xref: /petsc/src/mat/impls/aij/mpi/mpimatmatmult.c (revision ee3a9e216a68c71a318bf6447fcbaabcae63b2e7) 
1 
2 /*
3   Defines matrix-matrix product routines for pairs of MPIAIJ matrices
4           C = A * B
5 */
6 #include <../src/mat/impls/aij/seq/aij.h> /*I "petscmat.h" I*/
7 #include <../src/mat/utils/freespace.h>
8 #include <../src/mat/impls/aij/mpi/mpiaij.h>
9 #include <petscbt.h>
10 #include <../src/mat/impls/dense/mpi/mpidense.h>
11 /*
12 #define DEBUG_MATMATMULT
13  */
14 
15 #undef __FUNCT__
16 #define __FUNCT__ "MatMatMult_MPIAIJ_MPIAIJ"
17 PetscErrorCode MatMatMult_MPIAIJ_MPIAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill, Mat *C)
18 {
19   PetscErrorCode ierr;
20 
21   PetscFunctionBegin;
22   if (scall == MAT_INITIAL_MATRIX){
23     ierr = PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr);
24     ierr = MatMatMultSymbolic_MPIAIJ_MPIAIJ(A,B,fill,C);CHKERRQ(ierr);
25     ierr = PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr);
26   }
27   ierr = PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr);
28   ierr = (*(*C)->ops->matmultnumeric)(A,B,*C);CHKERRQ(ierr);
29   ierr = PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr);
30   PetscFunctionReturn(0);
31 }
32 
33 #undef __FUNCT__
34 #define __FUNCT__ "PetscContainerDestroy_Mat_MatMatMultMPI"
35 PetscErrorCode PetscContainerDestroy_Mat_MatMatMultMPI(void *ptr)
36 {
37   PetscErrorCode       ierr;
38   Mat_MatMatMultMPI    *mult=(Mat_MatMatMultMPI*)ptr;
39 
40   PetscFunctionBegin;
41   ierr = ISDestroy(&mult->isrowa);CHKERRQ(ierr);
42   ierr = ISDestroy(&mult->isrowb);CHKERRQ(ierr);
43   ierr = ISDestroy(&mult->iscolb);CHKERRQ(ierr);
44   ierr = MatDestroy(&mult->C_seq);CHKERRQ(ierr);
45   ierr = MatDestroy(&mult->A_loc);CHKERRQ(ierr);
46   ierr = MatDestroy(&mult->B_seq);CHKERRQ(ierr);
47   ierr = PetscFree(mult);CHKERRQ(ierr);
48   PetscFunctionReturn(0);
49 }
50 
51 #undef __FUNCT__
52 #define __FUNCT__ "MatDestroy_MPIAIJ_MatMatMult"
53 PetscErrorCode MatDestroy_MPIAIJ_MatMatMult(Mat A)
54 {
55   PetscErrorCode ierr;
56   Mat_MPIAIJ     *a=(Mat_MPIAIJ*)A->data;
57   Mat_PtAPMPI    *ptap=a->ptap;
58 
59   PetscFunctionBegin;
60   ierr = PetscFree2(ptap->startsj,ptap->startsj_r);CHKERRQ(ierr);
61   ierr = PetscFree(ptap->bufa);CHKERRQ(ierr);
62   ierr = MatDestroy(&ptap->P_loc);CHKERRQ(ierr);
63   ierr = MatDestroy(&ptap->P_oth);CHKERRQ(ierr);
64   ierr = PetscFree(ptap->api);CHKERRQ(ierr);
65   ierr = PetscFree(ptap->apj);CHKERRQ(ierr);
66   ierr = PetscFree(ptap->apa);CHKERRQ(ierr);
67   ierr = ptap->destroy(A);CHKERRQ(ierr);
68   ierr = PetscFree(ptap);CHKERRQ(ierr);
69   PetscFunctionReturn(0);
70 }
71 
72 #undef __FUNCT__
73 #define __FUNCT__ "MatDestroy_MPIAIJ_MatMatMult_32"
74 PetscErrorCode MatDestroy_MPIAIJ_MatMatMult_32(Mat A)
75 {
76   PetscErrorCode     ierr;
77   PetscContainer     container;
78   Mat_MatMatMultMPI  *mult=PETSC_NULL;
79 
80   PetscFunctionBegin;
81   ierr = PetscObjectQuery((PetscObject)A,"Mat_MatMatMultMPI",(PetscObject *)&container);CHKERRQ(ierr);
82   if (!container) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_PLIB,"Container does not exist");
83   ierr = PetscContainerGetPointer(container,(void **)&mult);CHKERRQ(ierr);
84   A->ops->destroy   = mult->destroy;
85   A->ops->duplicate = mult->duplicate;
86   if (A->ops->destroy) {
87     ierr = (*A->ops->destroy)(A);CHKERRQ(ierr);
88   }
89   ierr = PetscObjectCompose((PetscObject)A,"Mat_MatMatMultMPI",0);CHKERRQ(ierr);
90   PetscFunctionReturn(0);
91 }
92 
93 #undef __FUNCT__
94 #define __FUNCT__ "MatDuplicate_MPIAIJ_MatMatMult_32"
95 PetscErrorCode MatDuplicate_MPIAIJ_MatMatMult_32(Mat A, MatDuplicateOption op, Mat *M)
96 {
97   PetscErrorCode     ierr;
98   Mat_MatMatMultMPI  *mult;
99   PetscContainer     container;
100 
101   PetscFunctionBegin;
102   ierr = PetscObjectQuery((PetscObject)A,"Mat_MatMatMultMPI",(PetscObject *)&container);CHKERRQ(ierr);
103   if (!container) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_PLIB,"Container does not exist");
104   ierr  = PetscContainerGetPointer(container,(void **)&mult);CHKERRQ(ierr);
105   /* Note: the container is not duplicated, because it requires deep copying of
106      several large data sets (see PetscContainerDestroy_Mat_MatMatMultMPI()).
107      These data sets are only used for repeated calling of MatMatMultNumeric().
108      *M is unlikely being used in this way. Thus we create *M with pure mpiaij format */
109   ierr = (*mult->duplicate)(A,op,M);CHKERRQ(ierr);
110   (*M)->ops->destroy   = mult->destroy;   /* = MatDestroy_MPIAIJ, *M doesn't duplicate A's container! */
111   (*M)->ops->duplicate = mult->duplicate; /* = MatDuplicate_MPIAIJ */
112   PetscFunctionReturn(0);
113 }
114 
115 #undef __FUNCT__
116 #define __FUNCT__ "MatDuplicate_MPIAIJ_MatMatMult"
117 PetscErrorCode MatDuplicate_MPIAIJ_MatMatMult(Mat A, MatDuplicateOption op, Mat *M)
118 {
119   PetscErrorCode     ierr;
120   Mat_MPIAIJ         *a=(Mat_MPIAIJ*)A->data;
121   Mat_PtAPMPI        *ptap=a->ptap;
122 
123   PetscFunctionBegin;
124   ierr = (*ptap->duplicate)(A,op,M);CHKERRQ(ierr);
125   (*M)->ops->destroy   = ptap->destroy;   /* = MatDestroy_MPIAIJ, *M doesn't duplicate A's special structure! */
126   (*M)->ops->duplicate = ptap->duplicate; /* = MatDuplicate_MPIAIJ */
127   PetscFunctionReturn(0);
128 }
129 
130 #undef __FUNCT__
131 #define __FUNCT__ "MatMatMultNumeric_MPIAIJ_MPIAIJ"
132 PetscErrorCode MatMatMultNumeric_MPIAIJ_MPIAIJ(Mat A,Mat P,Mat C)
133 {
134   PetscErrorCode     ierr;
135   Mat_MPIAIJ         *a=(Mat_MPIAIJ*)A->data,*c=(Mat_MPIAIJ*)C->data;
136   Mat_SeqAIJ         *ad=(Mat_SeqAIJ*)(a->A)->data,*ao=(Mat_SeqAIJ*)(a->B)->data;
137   Mat_SeqAIJ         *cd=(Mat_SeqAIJ*)(c->A)->data,*co=(Mat_SeqAIJ*)(c->B)->data;
138   PetscInt           *adi=ad->i,*adj,*aoi=ao->i,*aoj;
139   PetscScalar        *ada,*aoa,*cda=cd->a,*coa=co->a;
140   Mat_SeqAIJ         *p_loc,*p_oth;
141   PetscInt           *pi_loc,*pj_loc,*pi_oth,*pj_oth,*pj;
142   PetscScalar        *pa_loc,*pa_oth,*pa,*apa,valtmp,*ca;
143   PetscInt           cm=C->rmap->n,anz,pnz;
144   Mat_PtAPMPI        *ptap=c->ptap;
145   PetscInt           *api,*apj,*apJ,i,j,k,row;
146   PetscInt           cstart=C->cmap->rstart;
147   PetscInt           cdnz,conz,k0,k1;
148 #if defined(DEBUG_MATMATMULT)
149   PetscMPIInt        rank=a->rank;
150 #endif
151 
152   PetscFunctionBegin;
153 #if defined(DEBUG_MATMATMULT)
154   if (!rank) ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] call MatMatMultNumeric_MPIAIJ_MPIAIJ()...\n",rank);
155 #endif
156 
157   /* 1) get P_oth = ptap->P_oth  and P_loc = ptap->P_loc */
158   /*-----------------------------------------------------*/
159   /* update numerical values of P_oth and P_loc */
160   ierr = MatGetBrowsOfAoCols_MPIAIJ(A,P,MAT_REUSE_MATRIX,&ptap->startsj,&ptap->startsj_r,&ptap->bufa,&ptap->P_oth);CHKERRQ(ierr);
161   ierr = MatMPIAIJGetLocalMat(P,MAT_REUSE_MATRIX,&ptap->P_loc);CHKERRQ(ierr);
162 #if defined(DEBUG_MATMATMULT)
163   if (!rank) ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] got P_oth and P_loc...\n",rank);
164 #endif
165 
166   /* 2) compute numeric C_loc = A_loc*P = Ad*P_loc + Ao*P_oth */
167   /*----------------------------------------------------------*/
168   /* get data from symbolic products */
169   p_loc = (Mat_SeqAIJ*)(ptap->P_loc)->data;
170   p_oth = (Mat_SeqAIJ*)(ptap->P_oth)->data;
171   pi_loc=p_loc->i; pj_loc=p_loc->j; pa_loc=p_loc->a;
172   pi_oth=p_oth->i; pj_oth=p_oth->j; pa_oth=p_oth->a;
173 
174   /* get apa for storing dense row A[i,:]*P */
175   apa = ptap->apa;
176 
177   api = ptap->api;
178   apj = ptap->apj;
179   for (i=0; i<cm; i++) {
180     /* diagonal portion of A */
181     anz = adi[i+1] - adi[i];
182     adj = ad->j + adi[i];
183     ada = ad->a + adi[i];
184     for (j=0; j<anz; j++) {
185       row = adj[j];
186       pnz = pi_loc[row+1] - pi_loc[row];
187       pj  = pj_loc + pi_loc[row];
188       pa  = pa_loc + pi_loc[row];
189 
190       /* perform dense axpy */
191       valtmp = ada[j];
192       for (k=0; k<pnz; k++){
193         apa[pj[k]] += valtmp*pa[k];
194       }
195       ierr = PetscLogFlops(2.0*pnz);CHKERRQ(ierr);
196     }
197 
198     /* off-diagonal portion of A */
199     anz = aoi[i+1] - aoi[i];
200     aoj = ao->j + aoi[i];
201     aoa = ao->a + aoi[i];
202     for (j=0; j<anz; j++) {
203       row = aoj[j];
204       pnz = pi_oth[row+1] - pi_oth[row];
205       pj  = pj_oth + pi_oth[row];
206       pa  = pa_oth + pi_oth[row];
207 
208       /* perform dense axpy */
209       valtmp = aoa[j];
210       for (k=0; k<pnz; k++){
211         apa[pj[k]] += valtmp*pa[k];
212       }
213       ierr = PetscLogFlops(2.0*pnz);CHKERRQ(ierr);
214     }
215 
216     /* set values in C */
217     apJ = apj + api[i];
218     cdnz = cd->i[i+1] - cd->i[i];
219     conz = co->i[i+1] - co->i[i];
220 
221     /* 1st off-diagoanl part of C */
222     ca = coa + co->i[i];
223     k  = 0;
224     for (k0=0; k0<conz; k0++){
225       if (apJ[k] >= cstart) break;
226       ca[k0]      = apa[apJ[k]];
227       apa[apJ[k]] = 0.0;
228       k++;
229     }
230 
231     /* diagonal part of C */
232     ca = cda + cd->i[i];
233     for (k1=0; k1<cdnz; k1++){
234       ca[k1]      = apa[apJ[k]];
235       apa[apJ[k]] = 0.0;
236       k++;
237     }
238 
239     /* 2nd off-diagoanl part of C */
240     ca = coa + co->i[i];
241     for (; k0<conz; k0++){
242       ca[k0]      = apa[apJ[k]];
243       apa[apJ[k]] = 0.0;
244       k++;
245     }
246   }
247   ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
248   ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
249   PetscFunctionReturn(0);
250 }
251 
252 #undef __FUNCT__
253 #define __FUNCT__ "MatMatMultSymbolic_MPIAIJ_MPIAIJ"
254 PetscErrorCode MatMatMultSymbolic_MPIAIJ_MPIAIJ(Mat A,Mat P,PetscReal fill,Mat *C)
255 {
256   PetscErrorCode       ierr;
257   MPI_Comm             comm=((PetscObject)A)->comm;
258   Mat                  Cmpi;
259   Mat_PtAPMPI          *ptap;
260   PetscFreeSpaceList   free_space=PETSC_NULL,current_space=PETSC_NULL;
261   Mat_MPIAIJ           *a=(Mat_MPIAIJ*)A->data,*c;
262   Mat_SeqAIJ           *ad=(Mat_SeqAIJ*)(a->A)->data,*ao=(Mat_SeqAIJ*)(a->B)->data,*p_loc,*p_oth;
263   PetscInt             *pi_loc,*pj_loc,*pi_oth,*pj_oth,*dnz,*onz;
264   PetscInt             *adi=ad->i,*adj=ad->j,*aoi=ao->i,*aoj=ao->j,rstart=A->rmap->rstart;
265   PetscInt             nlnk,*lnk,i,pnz,row,*api,*apj,*Jptr,apnz,nspacedouble=0,j,nzi;
266   PetscInt             am=A->rmap->n,pN=P->cmap->N,pn=P->cmap->n,pm=P->rmap->n;
267   PetscBT              lnkbt;
268   PetscScalar          *apa;
269   PetscReal            afill;
270   PetscBool            scalable=PETSC_FALSE;
271 #if defined(DEBUG_MATMATMULT)
272   PetscMPIInt          rank=a->rank;
273 #endif
274 
275   PetscFunctionBegin;
276   if (A->cmap->rstart != P->rmap->rstart || A->cmap->rend != P->rmap->rend){
277     SETERRQ4(comm,PETSC_ERR_ARG_SIZ,"Matrix local dimensions are incompatible, (%D, %D) != (%D,%D)",A->cmap->rstart,A->cmap->rend,P->rmap->rstart,P->rmap->rend);
278   }
279 
280   ierr = PetscObjectOptionsBegin((PetscObject)A);CHKERRQ(ierr);
281 
282     ierr = PetscOptionsBool("-matmatmult_32","Use a scalable but slower C=A*B","",scalable,&scalable,PETSC_NULL);CHKERRQ(ierr);
283     if (scalable){
284       ierr = MatMatMultSymbolic_MPIAIJ_MPIAIJ_Scalable_32(A,P,fill,C);;CHKERRQ(ierr);
285       PetscFunctionReturn(0);
286     }
287     ierr = PetscOptionsBool("-matmatmult_scalable","Use a scalable but slower C=A*B","",scalable,&scalable,PETSC_NULL);CHKERRQ(ierr);
288     if (scalable){
289       ierr = MatMatMultSymbolic_MPIAIJ_MPIAIJ_Scalable(A,P,fill,C);;CHKERRQ(ierr);
290       PetscFunctionReturn(0);
291     }
292   ierr = PetscOptionsEnd();CHKERRQ(ierr);
293 
294 #if defined(DEBUG_MATMATMULT)
295   if (!rank) ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] call MatMatMultSymbolic_MPIAIJ_MPIAIJ()...\n",rank);
296 #endif
297 
298   /* create struct Mat_PtAPMPI and attached it to C later */
299   ierr = PetscNew(Mat_PtAPMPI,&ptap);CHKERRQ(ierr);
300 
301   /* malloc apa to store dense row A[i,:]*P */
302   ierr = PetscMalloc(pN*sizeof(PetscScalar),&apa);CHKERRQ(ierr);
303   ierr = PetscMemzero(apa,pN*sizeof(PetscScalar));CHKERRQ(ierr);
304   ptap->apa = apa;
305 #if defined(DEBUG_MATMATMULT)
306   if (!rank) ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] Malloc apa pN %D is done...\n",rank,pN);
307 #endif
308 
309   /* get P_oth by taking rows of P (= non-zero cols of local A) from other processors */
310   ierr = MatGetBrowsOfAoCols_MPIAIJ(A,P,MAT_INITIAL_MATRIX,&ptap->startsj,&ptap->startsj_r,&ptap->bufa,&ptap->P_oth);CHKERRQ(ierr);
311 #if defined(DEBUG_MATMATMULT)
312   if (!rank) ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] P_oth is done...\n",rank);
313 #endif
314   /* get P_loc by taking all local rows of P */
315   ierr = MatMPIAIJGetLocalMat(P,MAT_INITIAL_MATRIX,&ptap->P_loc);CHKERRQ(ierr);
316 
317   p_loc = (Mat_SeqAIJ*)(ptap->P_loc)->data;
318   p_oth = (Mat_SeqAIJ*)(ptap->P_oth)->data;
319   pi_loc = p_loc->i; pj_loc = p_loc->j;
320   pi_oth = p_oth->i; pj_oth = p_oth->j;
321 
322 #if defined(DEBUG_MATMATMULT)
323   if (!rank) ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] P_loc is done, Annz %D * P_locnnz %D = %D...\n",rank,ad->i[am]+ao->i[am],p_loc->rmax,(ad->i[am]+ao->i[am])*p_loc->rmax);
324 #endif
325 
326   /* first, compute symbolic AP = A_loc*P = A_diag*P_loc + A_off*P_oth */
327   /*-------------------------------------------------------------------*/
328   ierr  = PetscMalloc((am+2)*sizeof(PetscInt),&api);CHKERRQ(ierr);
329   ptap->api = api;
330   api[0]    = 0;
331 
332   /* create and initialize a linked list */
333   nlnk = pN+1;
334   ierr = PetscLLCreate(pN,pN,nlnk,lnk,lnkbt);CHKERRQ(ierr);
335 
336   /* Initial FreeSpace size is fill*(nnz(A)+nnz(P)) */
337   ierr = PetscFreeSpaceGet((PetscInt)(fill*(adi[am]+aoi[am]+pi_loc[pm])),&free_space);CHKERRQ(ierr);
338   current_space = free_space;
339 
340   ierr = MatPreallocateInitialize(comm,am,pn,dnz,onz);CHKERRQ(ierr);
341   for (i=0; i<am; i++) {
342     apnz = 0;
343     /* diagonal portion of A */
344     nzi = adi[i+1] - adi[i];
345     for (j=0; j<nzi; j++){
346       row = *adj++;
347       pnz = pi_loc[row+1] - pi_loc[row];
348       Jptr  = pj_loc + pi_loc[row];
349       /* add non-zero cols of P into the sorted linked list lnk */
350       ierr = PetscLLAddSorted(pnz,Jptr,pN,nlnk,lnk,lnkbt);CHKERRQ(ierr);
351       apnz += nlnk;
352     }
353     /* off-diagonal portion of A */
354     nzi = aoi[i+1] - aoi[i];
355     for (j=0; j<nzi; j++){
356       row = *aoj++;
357       pnz = pi_oth[row+1] - pi_oth[row];
358       Jptr  = pj_oth + pi_oth[row];
359       ierr = PetscLLAddSorted(pnz,Jptr,pN,nlnk,lnk,lnkbt);CHKERRQ(ierr);
360       apnz += nlnk;
361     }
362 
363     api[i+1] = api[i] + apnz;
364 
365     /* if free space is not available, double the total space in the list */
366     if (current_space->local_remaining<apnz) {
367       ierr = PetscFreeSpaceGet(apnz+current_space->total_array_size,&current_space);CHKERRQ(ierr);
368       nspacedouble++;
369     }
370 
371     /* Copy data into free space, then initialize lnk */
372     ierr = PetscLLClean(pN,pN,apnz,lnk,current_space->array,lnkbt);CHKERRQ(ierr);
373     ierr = MatPreallocateSet(i+rstart,apnz,current_space->array,dnz,onz);CHKERRQ(ierr);
374     current_space->array           += apnz;
375     current_space->local_used      += apnz;
376     current_space->local_remaining -= apnz;
377   }
378 
379   /* Allocate space for apj, initialize apj, and */
380   /* destroy list of free space and other temporary array(s) */
381   ierr = PetscMalloc((api[am]+1)*sizeof(PetscInt),&ptap->apj);CHKERRQ(ierr);
382   apj  = ptap->apj;
383   ierr = PetscFreeSpaceContiguous(&free_space,ptap->apj);CHKERRQ(ierr);
384   ierr = PetscLLDestroy(lnk,lnkbt);CHKERRQ(ierr);
385 #if defined(DEBUG_MATMATMULT)
386   if (!rank) ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] AP is done...\n",rank);
387 #endif
388 
389   /* create and assemble symbolic parallel matrix Cmpi */
390   /*----------------------------------------------------*/
391   ierr = MatCreate(comm,&Cmpi);CHKERRQ(ierr);
392   ierr = MatSetSizes(Cmpi,am,pn,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr);
393   ierr = MatSetType(Cmpi,MATMPIAIJ);CHKERRQ(ierr);
394   ierr = MatMPIAIJSetPreallocation(Cmpi,0,dnz,0,onz);CHKERRQ(ierr);
395   ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr);
396   ierr = MatSetBlockSize(Cmpi,1);CHKERRQ(ierr);
397   for (i=0; i<am; i++){
398     row  = i + rstart;
399     apnz = api[i+1] - api[i];
400     ierr = MatSetValues(Cmpi,1,&row,apnz,apj,apa,INSERT_VALUES);CHKERRQ(ierr);
401     apj += apnz;
402   }
403   ierr = MatAssemblyBegin(Cmpi,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
404   ierr = MatAssemblyEnd(Cmpi,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
405 
406   ptap->destroy             = Cmpi->ops->destroy;
407   ptap->duplicate           = Cmpi->ops->duplicate;
408   Cmpi->ops->matmultnumeric = MatMatMultNumeric_MPIAIJ_MPIAIJ;
409   Cmpi->ops->destroy        = MatDestroy_MPIAIJ_MatMatMult;
410   Cmpi->ops->duplicate      = MatDuplicate_MPIAIJ_MatMatMult;
411 
412   /* attach the supporting struct to Cmpi for reuse */
413   c = (Mat_MPIAIJ*)Cmpi->data;
414   c->ptap  = ptap;
415 
416   *C = Cmpi;
417 
418   /* set MatInfo */
419   afill = (PetscReal)api[am]/(adi[am]+aoi[am]+pi_loc[pm]) + 1.e-5;
420   if (afill < 1.0) afill = 1.0;
421   Cmpi->info.mallocs           = nspacedouble;
422   Cmpi->info.fill_ratio_given  = fill;
423   Cmpi->info.fill_ratio_needed = afill;
424 
425 #if defined(PETSC_USE_INFO)
426   if (api[am]) {
427     ierr = PetscInfo3(Cmpi,"Reallocs %D; Fill ratio: given %G needed %G.\n",nspacedouble,fill,afill);CHKERRQ(ierr);
428     ierr = PetscInfo1(Cmpi,"Use MatMatMult(A,B,MatReuse,%G,&C) for best performance.;\n",afill);CHKERRQ(ierr);
429   } else {
430     ierr = PetscInfo(Cmpi,"Empty matrix product\n");CHKERRQ(ierr);
431   }
432 #endif
433   PetscFunctionReturn(0);
434 }
435 
436 /* implementation used in PETSc-3.2 */
437 /* This routine is called ONLY in the case of reusing previously computed symbolic C */
438 #undef __FUNCT__
439 #define __FUNCT__ "MatMatMultNumeric_MPIAIJ_MPIAIJ_Scalable_32"
440 PetscErrorCode MatMatMultNumeric_MPIAIJ_MPIAIJ_Scalable_32(Mat A,Mat B,Mat C)
441 {
442   PetscErrorCode     ierr;
443   Mat                *seq;
444   Mat_MatMatMultMPI  *mult;
445   PetscContainer     container;
446 
447   PetscFunctionBegin;
448   ierr = PetscObjectQuery((PetscObject)C,"Mat_MatMatMultMPI",(PetscObject *)&container);CHKERRQ(ierr);
449   if (!container) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_PLIB,"Container does not exist");
450   ierr  = PetscContainerGetPointer(container,(void **)&mult);CHKERRQ(ierr);
451 
452   if (mult->skipNumeric){ /* first numeric product is done during symbolic product */
453     mult->skipNumeric = PETSC_FALSE;
454     PetscFunctionReturn(0);
455   }
456 #if defined(DEBUG_MATMATMULT)
457   PetscMPIInt rank;
458   MPI_Comm    comm = ((PetscObject)C)->comm;
459   ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
460   ierr = MPI_Barrier(comm);CHKERRQ(ierr);
461   if (!rank) ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] call MatMatMultNumeric_MPIAIJ_MPIAIJ_Scalable_32()...\n",rank);
462 #endif
463 
464   seq = &mult->B_seq;
465   ierr = MatGetSubMatrices(B,1,&mult->isrowb,&mult->iscolb,MAT_REUSE_MATRIX,&seq);CHKERRQ(ierr);
466   mult->B_seq = *seq;
467 
468   seq = &mult->A_loc;
469   ierr = MatGetSubMatrices(A,1,&mult->isrowa,&mult->isrowb,MAT_REUSE_MATRIX,&seq);CHKERRQ(ierr);
470   mult->A_loc = *seq;
471 
472   ierr = MatMatMultNumeric_SeqAIJ_SeqAIJ_Scalable(mult->A_loc,mult->B_seq,mult->C_seq);CHKERRQ(ierr);
473 
474   ierr = PetscObjectReference((PetscObject)mult->C_seq);CHKERRQ(ierr);
475   ierr = MatMerge(((PetscObject)A)->comm,mult->C_seq,B->cmap->n,MAT_REUSE_MATRIX,&C);CHKERRQ(ierr);
476   PetscFunctionReturn(0);
477 }
478 
479 /* numeric product is computed as well */
480 #undef __FUNCT__
481 #define __FUNCT__ "MatMatMultSymbolic_MPIAIJ_MPIAIJ_Scalable_32"
482 PetscErrorCode MatMatMultSymbolic_MPIAIJ_MPIAIJ_Scalable_32(Mat A,Mat B,PetscReal fill,Mat *C)
483 {
484   PetscErrorCode     ierr;
485   Mat_MatMatMultMPI  *mult;
486   PetscContainer     container;
487   Mat                AB,*seq;
488   Mat_MPIAIJ         *a=(Mat_MPIAIJ*)A->data;
489   PetscInt           *idx,i,start,ncols,nzA,nzB,*cmap,imark;
490 #if defined(DEBUG_MATMATMULT)
491   MPI_Comm             comm = ((PetscObject)A)->comm;
492   PetscMPIInt          rank=a->rank;
493 #endif
494 
495   PetscFunctionBegin;
496 #if defined(DEBUG_MATMATMULT)
497   if (!rank) ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] call MatMatMultSymbolic_MPIAIJ_MPIAIJ_Scalable_32()...\n",rank);
498 #endif
499   if (A->cmap->rstart != B->rmap->rstart || A->cmap->rend != B->rmap->rend){
500     SETERRQ4(((PetscObject)A)->comm,PETSC_ERR_ARG_SIZ,"Matrix local dimensions are incompatible, (%D, %D) != (%D,%D)",A->cmap->rstart,A->cmap->rend,B->rmap->rstart,B->rmap->rend);
501   }
502 
503   ierr = PetscNew(Mat_MatMatMultMPI,&mult);CHKERRQ(ierr);
504 
505   /* get isrowb: nonzero col of A */
506   start = A->cmap->rstart;
507   cmap  = a->garray;
508   nzA   = a->A->cmap->n;
509   nzB   = a->B->cmap->n;
510   ierr  = PetscMalloc((nzA+nzB)*sizeof(PetscInt), &idx);CHKERRQ(ierr);
511   ncols = 0;
512   for (i=0; i<nzB; i++) {  /* row < local row index */
513     if (cmap[i] < start) idx[ncols++] = cmap[i];
514     else break;
515   }
516   imark = i;
517   for (i=0; i<nzA; i++) idx[ncols++] = start + i;  /* local rows */
518   for (i=imark; i<nzB; i++) idx[ncols++] = cmap[i]; /* row > local row index */
519   ierr = ISCreateGeneral(PETSC_COMM_SELF,ncols,idx,PETSC_OWN_POINTER,&mult->isrowb);CHKERRQ(ierr);
520   ierr = ISCreateStride(PETSC_COMM_SELF,B->cmap->N,0,1,&mult->iscolb);CHKERRQ(ierr);
521 
522   /*  get isrowa: all local rows of A */
523   ierr = ISCreateStride(PETSC_COMM_SELF,A->rmap->n,A->rmap->rstart,1,&mult->isrowa);CHKERRQ(ierr);
524 
525   /* Below should go to MatMatMultNumeric_MPIAIJ_MPIAIJ() - How to generate C there? */
526   /* create a seq matrix B_seq = submatrix of B by taking rows of B that equal to nonzero col of A */
527   ierr = MatGetSubMatrices(B,1,&mult->isrowb,&mult->iscolb,MAT_INITIAL_MATRIX,&seq);CHKERRQ(ierr);
528   mult->B_seq = *seq;
529   ierr = PetscFree(seq);CHKERRQ(ierr);
530 #if defined(DEBUG_MATMATMULT)
531   if (!rank) ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] B_seq is done...\n",rank);
532 #endif
533 
534   /*  create a seq matrix A_seq = submatrix of A by taking all local rows of A */
535   ierr = MatGetSubMatrices(A,1,&mult->isrowa,&mult->isrowb,MAT_INITIAL_MATRIX,&seq);CHKERRQ(ierr);
536   mult->A_loc = *seq;
537   ierr = PetscFree(seq);CHKERRQ(ierr);
538 #if defined(DEBUG_MATMATMULT)
539   if (!rank) ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] A_loc is done...\n",rank);
540 #endif
541 
542   /* compute C_seq = A_seq * B_seq */
543   ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable(mult->A_loc,mult->B_seq,fill,&mult->C_seq);CHKERRQ(ierr);
544 #if defined(DEBUG_MATMATMULT)
545   ierr = MPI_Barrier(comm);CHKERRQ(ierr);
546   if (!rank) ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] C_seq Symbolic is done...\n",rank);
547 #endif
548   ierr = MatMatMultNumeric_SeqAIJ_SeqAIJ_Scalable(mult->A_loc,mult->B_seq,mult->C_seq);CHKERRQ(ierr);
549 #if defined(DEBUG_MATMATMULT)
550   if (!rank) ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] C_seq Numeric is done...\n",rank);
551 #endif
552 
553   /* create mpi matrix C by concatinating C_seq */
554   ierr = PetscObjectReference((PetscObject)mult->C_seq);CHKERRQ(ierr); /* prevent C_seq being destroyed by MatMerge() */
555   ierr = MatMergeSymbolic(((PetscObject)A)->comm,mult->C_seq,B->cmap->n,&AB);CHKERRQ(ierr);
556   ierr = MatMergeNumeric(((PetscObject)A)->comm,mult->C_seq,B->cmap->n,AB);CHKERRQ(ierr);
557 #if defined(DEBUG_MATMATMULT)
558   if (!rank) ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] Merge is done...\n",rank);
559 #endif
560 
561   /* attach the supporting struct to C for reuse of symbolic C */
562   ierr = PetscContainerCreate(PETSC_COMM_SELF,&container);CHKERRQ(ierr);
563   ierr = PetscContainerSetPointer(container,mult);CHKERRQ(ierr);
564   ierr = PetscContainerSetUserDestroy(container,PetscContainerDestroy_Mat_MatMatMultMPI);CHKERRQ(ierr);
565   ierr = PetscObjectCompose((PetscObject)AB,"Mat_MatMatMultMPI",(PetscObject)container);CHKERRQ(ierr);
566   ierr = PetscContainerDestroy(&container);CHKERRQ(ierr);
567   mult->skipNumeric  = PETSC_TRUE; /* a numeric product is done here */
568   mult->destroy      = AB->ops->destroy;
569   mult->duplicate    = AB->ops->duplicate;
570   AB->ops->matmultnumeric = MatMatMultNumeric_MPIAIJ_MPIAIJ_Scalable_32;
571   AB->ops->destroy        = MatDestroy_MPIAIJ_MatMatMult_32;
572   AB->ops->duplicate      = MatDuplicate_MPIAIJ_MatMatMult_32;
573   AB->ops->matmult        = MatMatMult_MPIAIJ_MPIAIJ;
574   *C                      = AB;
575   PetscFunctionReturn(0);
576 }
577 
578 #undef __FUNCT__
579 #define __FUNCT__ "MatMatMult_MPIAIJ_MPIDense"
580 PetscErrorCode MatMatMult_MPIAIJ_MPIDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
581 {
582   PetscErrorCode ierr;
583 
584   PetscFunctionBegin;
585   if (scall == MAT_INITIAL_MATRIX){
586     ierr = MatMatMultSymbolic_MPIAIJ_MPIDense(A,B,fill,C);CHKERRQ(ierr);
587   }
588   ierr = MatMatMultNumeric_MPIAIJ_MPIDense(A,B,*C);CHKERRQ(ierr);
589   PetscFunctionReturn(0);
590 }
591 
592 typedef struct {
593   Mat         workB;
594   PetscScalar *rvalues,*svalues;
595   MPI_Request *rwaits,*swaits;
596 } MPIAIJ_MPIDense;
597 
598 #undef __FUNCT__
599 #define __FUNCT__ "MPIAIJ_MPIDenseDestroy"
600 PetscErrorCode MPIAIJ_MPIDenseDestroy(void *ctx)
601 {
602   MPIAIJ_MPIDense *contents = (MPIAIJ_MPIDense*) ctx;
603   PetscErrorCode  ierr;
604 
605   PetscFunctionBegin;
606   ierr = MatDestroy(&contents->workB);CHKERRQ(ierr);
607   ierr = PetscFree4(contents->rvalues,contents->svalues,contents->rwaits,contents->swaits);CHKERRQ(ierr);
608   ierr = PetscFree(contents);CHKERRQ(ierr);
609   PetscFunctionReturn(0);
610 }
611 
612 #undef __FUNCT__
613 #define __FUNCT__ "MatMatMultSymbolic_MPIAIJ_MPIDense"
614 PetscErrorCode MatMatMultSymbolic_MPIAIJ_MPIDense(Mat A,Mat B,PetscReal fill,Mat *C)
615 {
616   PetscErrorCode         ierr;
617   Mat_MPIAIJ             *aij = (Mat_MPIAIJ*) A->data;
618   PetscInt               nz = aij->B->cmap->n;
619   PetscContainer         container;
620   MPIAIJ_MPIDense        *contents;
621   VecScatter             ctx = aij->Mvctx;
622   VecScatter_MPI_General *from = (VecScatter_MPI_General*) ctx->fromdata;
623   VecScatter_MPI_General *to   = ( VecScatter_MPI_General*) ctx->todata;
624   PetscInt               m=A->rmap->n,n=B->cmap->n;
625 
626   PetscFunctionBegin;
627   ierr = MatCreate(((PetscObject)B)->comm,C);CHKERRQ(ierr);
628   ierr = MatSetSizes(*C,m,n,A->rmap->N,B->cmap->N);CHKERRQ(ierr);
629   ierr = MatSetType(*C,MATMPIDENSE);CHKERRQ(ierr);
630   ierr = MatAssemblyBegin(*C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
631   ierr = MatAssemblyEnd(*C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
632   (*C)->ops->matmult = MatMatMult_MPIAIJ_MPIDense;
633 
634   ierr = PetscNew(MPIAIJ_MPIDense,&contents);CHKERRQ(ierr);
635   /* Create work matrix used to store off processor rows of B needed for local product */
636   ierr = MatCreateSeqDense(PETSC_COMM_SELF,nz,B->cmap->N,PETSC_NULL,&contents->workB);CHKERRQ(ierr);
637   /* Create work arrays needed */
638   ierr = PetscMalloc4(B->cmap->N*from->starts[from->n],PetscScalar,&contents->rvalues,
639                       B->cmap->N*to->starts[to->n],PetscScalar,&contents->svalues,
640                       from->n,MPI_Request,&contents->rwaits,
641                       to->n,MPI_Request,&contents->swaits);CHKERRQ(ierr);
642 
643   ierr = PetscContainerCreate(((PetscObject)A)->comm,&container);CHKERRQ(ierr);
644   ierr = PetscContainerSetPointer(container,contents);CHKERRQ(ierr);
645   ierr = PetscContainerSetUserDestroy(container,MPIAIJ_MPIDenseDestroy);CHKERRQ(ierr);
646   ierr = PetscObjectCompose((PetscObject)(*C),"workB",(PetscObject)container);CHKERRQ(ierr);
647   ierr = PetscContainerDestroy(&container);CHKERRQ(ierr);
648   PetscFunctionReturn(0);
649 }
650 
651 #undef __FUNCT__
652 #define __FUNCT__ "MatMPIDenseScatter"
653 /*
654     Performs an efficient scatter on the rows of B needed by this process; this is
655     a modification of the VecScatterBegin_() routines.
656 */
657 PetscErrorCode MatMPIDenseScatter(Mat A,Mat B,Mat C,Mat *outworkB)
658 {
659   Mat_MPIAIJ             *aij = (Mat_MPIAIJ*)A->data;
660   PetscErrorCode         ierr;
661   PetscScalar            *b,*w,*svalues,*rvalues;
662   VecScatter             ctx = aij->Mvctx;
663   VecScatter_MPI_General *from = (VecScatter_MPI_General*) ctx->fromdata;
664   VecScatter_MPI_General *to   = ( VecScatter_MPI_General*) ctx->todata;
665   PetscInt               i,j,k;
666   PetscInt               *sindices,*sstarts,*rindices,*rstarts;
667   PetscMPIInt            *sprocs,*rprocs,nrecvs;
668   MPI_Request            *swaits,*rwaits;
669   MPI_Comm               comm = ((PetscObject)A)->comm;
670   PetscMPIInt            tag = ((PetscObject)ctx)->tag,ncols = B->cmap->N, nrows = aij->B->cmap->n,imdex,nrowsB = B->rmap->n;
671   MPI_Status             status;
672   MPIAIJ_MPIDense        *contents;
673   PetscContainer         container;
674   Mat                    workB;
675 
676   PetscFunctionBegin;
677   ierr = PetscObjectQuery((PetscObject)C,"workB",(PetscObject*)&container);CHKERRQ(ierr);
678   if (!container) SETERRQ(comm,PETSC_ERR_PLIB,"Container does not exist");
679   ierr = PetscContainerGetPointer(container,(void**)&contents);CHKERRQ(ierr);
680 
681   workB = *outworkB = contents->workB;
682   if (nrows != workB->rmap->n) SETERRQ2(comm,PETSC_ERR_PLIB,"Number of rows of workB %D not equal to columns of aij->B %D",nrows,workB->cmap->n);
683   sindices  = to->indices;
684   sstarts   = to->starts;
685   sprocs    = to->procs;
686   swaits    = contents->swaits;
687   svalues   = contents->svalues;
688 
689   rindices  = from->indices;
690   rstarts   = from->starts;
691   rprocs    = from->procs;
692   rwaits    = contents->rwaits;
693   rvalues   = contents->rvalues;
694 
695   ierr = MatGetArray(B,&b);CHKERRQ(ierr);
696   ierr = MatGetArray(workB,&w);CHKERRQ(ierr);
697 
698   for (i=0; i<from->n; i++) {
699     ierr = MPI_Irecv(rvalues+ncols*rstarts[i],ncols*(rstarts[i+1]-rstarts[i]),MPIU_SCALAR,rprocs[i],tag,comm,rwaits+i);CHKERRQ(ierr);
700   }
701 
702   for (i=0; i<to->n; i++) {
703     /* pack a message at a time */
704     CHKMEMQ;
705     for (j=0; j<sstarts[i+1]-sstarts[i]; j++){
706       for (k=0; k<ncols; k++) {
707         svalues[ncols*(sstarts[i] + j) + k] = b[sindices[sstarts[i]+j] + nrowsB*k];
708       }
709     }
710     CHKMEMQ;
711     ierr = MPI_Isend(svalues+ncols*sstarts[i],ncols*(sstarts[i+1]-sstarts[i]),MPIU_SCALAR,sprocs[i],tag,comm,swaits+i);CHKERRQ(ierr);
712   }
713 
714   nrecvs = from->n;
715   while (nrecvs) {
716     ierr = MPI_Waitany(from->n,rwaits,&imdex,&status);CHKERRQ(ierr);
717     nrecvs--;
718     /* unpack a message at a time */
719     CHKMEMQ;
720     for (j=0; j<rstarts[imdex+1]-rstarts[imdex]; j++){
721       for (k=0; k<ncols; k++) {
722         w[rindices[rstarts[imdex]+j] + nrows*k] = rvalues[ncols*(rstarts[imdex] + j) + k];
723       }
724     }
725     CHKMEMQ;
726   }
727   if (to->n) {ierr = MPI_Waitall(to->n,swaits,to->sstatus);CHKERRQ(ierr);}
728 
729   ierr = MatRestoreArray(B,&b);CHKERRQ(ierr);
730   ierr = MatRestoreArray(workB,&w);CHKERRQ(ierr);
731   ierr = MatAssemblyBegin(workB,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
732   ierr = MatAssemblyEnd(workB,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
733   PetscFunctionReturn(0);
734 }
735 extern PetscErrorCode MatMatMultNumericAdd_SeqAIJ_SeqDense(Mat,Mat,Mat);
736 
737 #undef __FUNCT__
738 #define __FUNCT__ "MatMatMultNumeric_MPIAIJ_MPIDense"
739 PetscErrorCode MatMatMultNumeric_MPIAIJ_MPIDense(Mat A,Mat B,Mat C)
740 {
741   PetscErrorCode       ierr;
742   Mat_MPIAIJ           *aij = (Mat_MPIAIJ*)A->data;
743   Mat_MPIDense         *bdense = (Mat_MPIDense*)B->data;
744   Mat_MPIDense         *cdense = (Mat_MPIDense*)C->data;
745   Mat                  workB;
746 
747   PetscFunctionBegin;
748 
749   /* diagonal block of A times all local rows of B*/
750   ierr = MatMatMultNumeric_SeqAIJ_SeqDense(aij->A,bdense->A,cdense->A);CHKERRQ(ierr);
751 
752   /* get off processor parts of B needed to complete the product */
753   ierr = MatMPIDenseScatter(A,B,C,&workB);CHKERRQ(ierr);
754 
755   /* off-diagonal block of A times nonlocal rows of B */
756   ierr = MatMatMultNumericAdd_SeqAIJ_SeqDense(aij->B,workB,cdense->A);CHKERRQ(ierr);
757   ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
758   ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
759   PetscFunctionReturn(0);
760 }
761 
762 #undef __FUNCT__
763 #define __FUNCT__ "MatMatMultNumeric_MPIAIJ_MPIAIJ_Scalable"
764 PetscErrorCode MatMatMultNumeric_MPIAIJ_MPIAIJ_Scalable(Mat A,Mat P,Mat C)
765 {
766   PetscErrorCode     ierr;
767   Mat_MPIAIJ         *a=(Mat_MPIAIJ*)A->data,*c=(Mat_MPIAIJ*)C->data;
768   Mat_SeqAIJ         *ad=(Mat_SeqAIJ*)(a->A)->data,*ao=(Mat_SeqAIJ*)(a->B)->data;
769   Mat_SeqAIJ         *cd=(Mat_SeqAIJ*)(c->A)->data,*co=(Mat_SeqAIJ*)(c->B)->data;
770   PetscInt           *adi=ad->i,*adj,*aoi=ao->i,*aoj;
771   PetscScalar        *ada,*aoa,*cda=cd->a,*coa=co->a;
772   Mat_SeqAIJ         *p_loc,*p_oth;
773   PetscInt           *pi_loc,*pj_loc,*pi_oth,*pj_oth,*pj;
774   PetscScalar        *pa_loc,*pa_oth,*pa,valtmp,*ca;
775   PetscInt           cm=C->rmap->n,anz,pnz;
776   Mat_PtAPMPI        *ptap=c->ptap;
777   PetscScalar        *apa_sparse=ptap->apa;
778   PetscInt           *api,*apj,*apJ,i,j,k,row;
779   PetscInt           cstart=C->cmap->rstart;
780   PetscInt           cdnz,conz,k0,k1,nextp;
781 #if defined(DEBUG_MATMATMULT)
782   PetscMPIInt        rank=a->rank;
783 #endif
784 
785   PetscFunctionBegin;
786 #if defined(DEBUG_MATMATMULT)
787   if (!rank) ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] call MatMatMultNumeric_MPIAIJ_MPIAIJ_Scalable()...\n",rank);
788 #endif
789 
790   /* 1) get P_oth = ptap->P_oth  and P_loc = ptap->P_loc */
791   /*-----------------------------------------------------*/
792   /* update numerical values of P_oth and P_loc */
793   ierr = MatGetBrowsOfAoCols_MPIAIJ(A,P,MAT_REUSE_MATRIX,&ptap->startsj,&ptap->startsj_r,&ptap->bufa,&ptap->P_oth);CHKERRQ(ierr);
794   ierr = MatMPIAIJGetLocalMat(P,MAT_REUSE_MATRIX,&ptap->P_loc);CHKERRQ(ierr);
795 #if defined(DEBUG_MATMATMULT)
796   if (!rank) ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] got P_oth and P_loc...\n",rank);
797 #endif
798 
799   /* 2) compute numeric C_loc = A_loc*P = Ad*P_loc + Ao*P_oth */
800   /*----------------------------------------------------------*/
801   /* get data from symbolic products */
802   p_loc = (Mat_SeqAIJ*)(ptap->P_loc)->data;
803   p_oth = (Mat_SeqAIJ*)(ptap->P_oth)->data;
804   pi_loc=p_loc->i; pj_loc=p_loc->j; pa_loc=p_loc->a;
805   pi_oth=p_oth->i; pj_oth=p_oth->j; pa_oth=p_oth->a;
806 
807   api = ptap->api;
808   apj = ptap->apj;
809   for (i=0; i<cm; i++) {
810     apJ = apj + api[i];
811 
812     /* diagonal portion of A */
813     anz = adi[i+1] - adi[i];
814     adj = ad->j + adi[i];
815     ada = ad->a + adi[i];
816     for (j=0; j<anz; j++) {
817       row = adj[j];
818       pnz = pi_loc[row+1] - pi_loc[row];
819       pj  = pj_loc + pi_loc[row];
820       pa  = pa_loc + pi_loc[row];
821       /* perform sparse axpy */
822       valtmp = ada[j];
823       nextp  = 0;
824       for (k=0; nextp<pnz; k++) {
825         if (apJ[k] == pj[nextp]) { /* column of AP == column of P */
826           apa_sparse[k] += valtmp*pa[nextp++];
827         }
828       }
829       ierr = PetscLogFlops(2.0*pnz);CHKERRQ(ierr);
830     }
831 
832     /* off-diagonal portion of A */
833     anz = aoi[i+1] - aoi[i];
834     aoj = ao->j + aoi[i];
835     aoa = ao->a + aoi[i];
836     for (j=0; j<anz; j++) {
837       row = aoj[j];
838       pnz = pi_oth[row+1] - pi_oth[row];
839       pj  = pj_oth + pi_oth[row];
840       pa  = pa_oth + pi_oth[row];
841       /* perform sparse axpy */
842       valtmp = aoa[j];
843       nextp  = 0;
844       for (k=0; nextp<pnz; k++) {
845         if (apJ[k] == pj[nextp]) { /* column of AP == column of P */
846           apa_sparse[k] += valtmp*pa[nextp++];
847         }
848       }
849       ierr = PetscLogFlops(2.0*pnz);CHKERRQ(ierr);
850     }
851 
852     /* set values in C */
853     cdnz = cd->i[i+1] - cd->i[i];
854     conz = co->i[i+1] - co->i[i];
855 
856     /* 1st off-diagoanl part of C */
857     ca = coa + co->i[i];
858     k  = 0;
859     for (k0=0; k0<conz; k0++){
860       if (apJ[k] >= cstart) break;
861       ca[k0]      = apa_sparse[k];
862       apa_sparse[k] = 0.0;
863       k++;
864     }
865 
866     /* diagonal part of C */
867     ca = cda + cd->i[i];
868     for (k1=0; k1<cdnz; k1++){
869       ca[k1]      = apa_sparse[k];
870       apa_sparse[k] = 0.0;
871       k++;
872     }
873 
874     /* 2nd off-diagoanl part of C */
875     ca = coa + co->i[i];
876     for (; k0<conz; k0++){
877       ca[k0]      = apa_sparse[k];
878       apa_sparse[k] = 0.0;
879       k++;
880     }
881   }
882   ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
883   ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
884   PetscFunctionReturn(0);
885 }
886 
887 /* same as MatMatMultSymbolic_MPIAIJ_MPIAIJ(), except using LLCondensed to avoid O(BN) memory requirement */
888 #undef __FUNCT__
889 #define __FUNCT__ "MatMatMultSymbolic_MPIAIJ_MPIAIJ_Scalable"
890 PetscErrorCode MatMatMultSymbolic_MPIAIJ_MPIAIJ_Scalable(Mat A,Mat P,PetscReal fill,Mat *C)
891 {
892   PetscErrorCode       ierr;
893   MPI_Comm             comm=((PetscObject)A)->comm;
894   Mat                  Cmpi;
895   Mat_PtAPMPI          *ptap;
896   PetscFreeSpaceList   free_space=PETSC_NULL,current_space=PETSC_NULL;
897   Mat_MPIAIJ           *a=(Mat_MPIAIJ*)A->data,*c;
898   Mat_SeqAIJ           *ad=(Mat_SeqAIJ*)(a->A)->data,*ao=(Mat_SeqAIJ*)(a->B)->data,*p_loc,*p_oth;
899   PetscInt             *pi_loc,*pj_loc,*pi_oth,*pj_oth,*dnz,*onz;
900   PetscInt             *adi=ad->i,*adj=ad->j,*aoi=ao->i,*aoj=ao->j,rstart=A->rmap->rstart;
901   PetscInt             i,pnz,row,*api,*apj,*Jptr,apnz,nspacedouble=0,j,nzi,*nlnk,*lnk,apnz_max=0;
902   PetscInt             am=A->rmap->n,pN=P->cmap->N,pn=P->cmap->n,pm=P->rmap->n;
903   PetscInt             nlnk_max,armax,prmax;
904   PetscBT              lnkbt;
905   PetscReal            afill;
906   PetscScalar          *apa;
907 #if defined(DEBUG_MATMATMULT)
908   PetscMPIInt          rank=a->rank;
909 #endif
910 
911   PetscFunctionBegin;
912 #if defined(DEBUG_MATMATMULT)
913   ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
914   if (!rank) ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] call MatMatMultSymbolic_MPIAIJ_MPIAIJ_Scalable()...\n",rank);
915 #endif
916 
917   /* create struct Mat_PtAPMPI and attached it to C later */
918   ierr = PetscNew(Mat_PtAPMPI,&ptap);CHKERRQ(ierr);
919 
920   /* get P_oth by taking rows of P (= non-zero cols of local A) from other processors */
921   ierr = MatGetBrowsOfAoCols_MPIAIJ(A,P,MAT_INITIAL_MATRIX,&ptap->startsj,&ptap->startsj_r,&ptap->bufa,&ptap->P_oth);CHKERRQ(ierr);
922 #if defined(DEBUG_MATMATMULT)
923   if (!rank) ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] P_oth is done...\n",rank);
924 #endif
925   /* get P_loc by taking all local rows of P */
926   ierr = MatMPIAIJGetLocalMat(P,MAT_INITIAL_MATRIX,&ptap->P_loc);CHKERRQ(ierr);
927 
928   p_loc = (Mat_SeqAIJ*)(ptap->P_loc)->data;
929   p_oth = (Mat_SeqAIJ*)(ptap->P_oth)->data;
930   pi_loc = p_loc->i; pj_loc = p_loc->j;
931   pi_oth = p_oth->i; pj_oth = p_oth->j;
932 
933 #if defined(DEBUG_MATMATMULT)
934   if (!rank) ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] P_loc is done, Annz %D * P_locnnz %D = %D...\n",rank,ad->i[am]+ao->i[am],p_loc->rmax,(ad->i[am]+ao->i[am])*p_loc->rmax);
935 #endif
936 
937   /* first, compute symbolic AP = A_loc*P = A_diag*P_loc + A_off*P_oth */
938   /*-------------------------------------------------------------------*/
939   ierr  = PetscMalloc((am+2)*sizeof(PetscInt),&api);CHKERRQ(ierr);
940   ptap->api = api;
941   api[0]    = 0;
942 
943   /* create and initialize a linked list */
944   armax = ad->rmax+ao->rmax;
945   prmax = PetscMax(p_loc->rmax,p_oth->rmax);
946   nlnk_max = armax*prmax;
947   if (!nlnk_max || nlnk_max > pN) nlnk_max = pN;
948 #if defined(DEBUG_MATMATMULT)
949   if (!rank) ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] pN %d; nlnk_max %d; Armax %d+%d=%d; Prmax Max(%d,%d)=%d\n",rank,pN,nlnk_max,ad->rmax,ao->rmax,armax,p_loc->rmax,p_oth->rmax,prmax);
950 #endif
951   ierr = PetscLLCondensedCreate(nlnk_max,pN,lnk,nlnk,lnkbt);
952 
953   /* Initial FreeSpace size is fill*(nnz(A)+nnz(P)) */
954   ierr = PetscFreeSpaceGet((PetscInt)(fill*(adi[am]+aoi[am]+pi_loc[pm])),&free_space);CHKERRQ(ierr);
955   current_space = free_space;
956 
957   ierr = MatPreallocateInitialize(comm,am,pn,dnz,onz);CHKERRQ(ierr);
958   for (i=0; i<am; i++) {
959     apnz = 0;
960     /* diagonal portion of A */
961     nzi = adi[i+1] - adi[i];
962     for (j=0; j<nzi; j++){
963       row = *adj++;
964       pnz = pi_loc[row+1] - pi_loc[row];
965       Jptr  = pj_loc + pi_loc[row];
966       /* add non-zero cols of P into the sorted linked list lnk */
967       ierr = PetscLLCondensedAddSorted(nlnk_max,pN,pnz,Jptr,nlnk,lnk,lnkbt);CHKERRQ(ierr);
968     }
969     /* off-diagonal portion of A */
970     nzi = aoi[i+1] - aoi[i];
971     for (j=0; j<nzi; j++){
972       row = *aoj++;
973       pnz = pi_oth[row+1] - pi_oth[row];
974       Jptr  = pj_oth + pi_oth[row];
975       ierr = PetscLLCondensedAddSorted(nlnk_max,pN,pnz,Jptr,nlnk,lnk,lnkbt);CHKERRQ(ierr);
976     }
977 
978     apnz     = *nlnk;
979     api[i+1] = api[i] + apnz;
980     if (apnz > apnz_max) apnz_max = apnz;
981 
982     /* if free space is not available, double the total space in the list */
983     if (current_space->local_remaining<apnz) {
984       ierr = PetscFreeSpaceGet(apnz+current_space->total_array_size,&current_space);CHKERRQ(ierr);
985       nspacedouble++;
986     }
987 
988     /* Copy data into free space, then initialize lnk */
989     ierr = PetscLLCondensedClean(nlnk_max,pN,apnz,current_space->array,nlnk,lnk,lnkbt);CHKERRQ(ierr);
990     ierr = MatPreallocateSet(i+rstart,apnz,current_space->array,dnz,onz);CHKERRQ(ierr);
991     current_space->array           += apnz;
992     current_space->local_used      += apnz;
993     current_space->local_remaining -= apnz;
994   }
995 
996   /* Allocate space for apj, initialize apj, and */
997   /* destroy list of free space and other temporary array(s) */
998   ierr = PetscMalloc((api[am]+1)*sizeof(PetscInt),&ptap->apj);CHKERRQ(ierr);
999   apj  = ptap->apj;
1000   ierr = PetscFreeSpaceContiguous(&free_space,ptap->apj);CHKERRQ(ierr);
1001   ierr = PetscLLCondensedDestroy(lnk,lnkbt);CHKERRQ(ierr);
1002 #if defined(DEBUG_MATMATMULT)
1003   if (!rank) ierr = PetscPrintf(PETSC_COMM_SELF,"[%d] AP is done..., apnz_max %d\n",rank,apnz_max);
1004 #endif
1005 
1006   /* create and assemble symbolic parallel matrix Cmpi */
1007   /*----------------------------------------------------*/
1008   ierr = MatCreate(comm,&Cmpi);CHKERRQ(ierr);
1009   ierr = MatSetSizes(Cmpi,am,pn,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr);
1010   ierr = MatSetType(Cmpi,MATMPIAIJ);CHKERRQ(ierr);
1011   ierr = MatMPIAIJSetPreallocation(Cmpi,0,dnz,0,onz);CHKERRQ(ierr);
1012   ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr);
1013   ierr = MatSetBlockSize(Cmpi,1);CHKERRQ(ierr);
1014 
1015   /* malloc apa for assembly Cmpi */
1016   ierr = PetscMalloc(apnz_max*sizeof(PetscScalar),&apa);CHKERRQ(ierr);
1017   ierr = PetscMemzero(apa,apnz_max*sizeof(PetscScalar));CHKERRQ(ierr);
1018   ptap->apa = apa;
1019   for (i=0; i<am; i++){
1020     row  = i + rstart;
1021     apnz = api[i+1] - api[i];
1022     ierr = MatSetValues(Cmpi,1,&row,apnz,apj,apa,INSERT_VALUES);CHKERRQ(ierr);
1023     apj += apnz;
1024   }
1025   ierr = MatAssemblyBegin(Cmpi,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1026   ierr = MatAssemblyEnd(Cmpi,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1027 
1028   ptap->destroy             = Cmpi->ops->destroy;
1029   ptap->duplicate           = Cmpi->ops->duplicate;
1030   Cmpi->ops->matmultnumeric = MatMatMultNumeric_MPIAIJ_MPIAIJ_Scalable;
1031   Cmpi->ops->destroy        = MatDestroy_MPIAIJ_MatMatMult;
1032   Cmpi->ops->duplicate      = MatDuplicate_MPIAIJ_MatMatMult;
1033 
1034   /* attach the supporting struct to Cmpi for reuse */
1035   c = (Mat_MPIAIJ*)Cmpi->data;
1036   c->ptap  = ptap;
1037 
1038   *C = Cmpi;
1039 
1040   /* set MatInfo */
1041   afill = (PetscReal)api[am]/(adi[am]+aoi[am]+pi_loc[pm]) + 1.e-5;
1042   if (afill < 1.0) afill = 1.0;
1043   Cmpi->info.mallocs           = nspacedouble;
1044   Cmpi->info.fill_ratio_given  = fill;
1045   Cmpi->info.fill_ratio_needed = afill;
1046 
1047 #if defined(PETSC_USE_INFO)
1048   if (api[am]) {
1049     ierr = PetscInfo3(Cmpi,"Reallocs %D; Fill ratio: given %G needed %G.\n",nspacedouble,fill,afill);CHKERRQ(ierr);
1050     ierr = PetscInfo1(Cmpi,"Use MatMatMult(A,B,MatReuse,%G,&C) for best performance.;\n",afill);CHKERRQ(ierr);
1051   } else {
1052     ierr = PetscInfo(Cmpi,"Empty matrix product\n");CHKERRQ(ierr);
1053   }
1054 #endif
1055   PetscFunctionReturn(0);
1056 }
1057