xref: /petsc/src/mat/impls/aij/mpi/mpimatmatmult.c (revision 2205254efee3a00a594e5e2a3a70f74dcb40bc03) 
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 #undef __FUNCT__
13 #define __FUNCT__ "MatMatMult_MPIAIJ_MPIAIJ"
14 PetscErrorCode MatMatMult_MPIAIJ_MPIAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill, Mat *C)
15 {
16   PetscErrorCode ierr;
17 
18   PetscFunctionBegin;
19   if (scall == MAT_INITIAL_MATRIX) {
20     ierr = MatMatMultSymbolic_MPIAIJ_MPIAIJ(A,B,fill,C);CHKERRQ(ierr);
21   }
22 
23   ierr = (*(*C)->ops->matmultnumeric)(A,B,*C);CHKERRQ(ierr);
24   PetscFunctionReturn(0);
25 }
26 
27 #undef __FUNCT__
28 #define __FUNCT__ "MatDestroy_MPIAIJ_MatMatMult"
29 PetscErrorCode MatDestroy_MPIAIJ_MatMatMult(Mat A)
30 {
31   PetscErrorCode ierr;
32   Mat_MPIAIJ     *a    = (Mat_MPIAIJ*)A->data;
33   Mat_PtAPMPI    *ptap = a->ptap;
34 
35   PetscFunctionBegin;
36   ierr = PetscFree2(ptap->startsj_s,ptap->startsj_r);CHKERRQ(ierr);
37   ierr = PetscFree(ptap->bufa);CHKERRQ(ierr);
38   ierr = MatDestroy(&ptap->P_loc);CHKERRQ(ierr);
39   ierr = MatDestroy(&ptap->P_oth);CHKERRQ(ierr);
40   ierr = MatDestroy(&ptap->Pt);CHKERRQ(ierr);
41   ierr = PetscFree(ptap->api);CHKERRQ(ierr);
42   ierr = PetscFree(ptap->apj);CHKERRQ(ierr);
43   ierr = PetscFree(ptap->apa);CHKERRQ(ierr);
44   ierr = ptap->destroy(A);CHKERRQ(ierr);
45   ierr = PetscFree(ptap);CHKERRQ(ierr);
46   PetscFunctionReturn(0);
47 }
48 
49 #undef __FUNCT__
50 #define __FUNCT__ "MatDuplicate_MPIAIJ_MatMatMult"
51 PetscErrorCode MatDuplicate_MPIAIJ_MatMatMult(Mat A, MatDuplicateOption op, Mat *M)
52 {
53   PetscErrorCode ierr;
54   Mat_MPIAIJ     *a    = (Mat_MPIAIJ*)A->data;
55   Mat_PtAPMPI    *ptap = a->ptap;
56 
57   PetscFunctionBegin;
58   ierr = (*ptap->duplicate)(A,op,M);CHKERRQ(ierr);
59   (*M)->ops->destroy   = ptap->destroy;   /* = MatDestroy_MPIAIJ, *M doesn't duplicate A's special structure! */
60   (*M)->ops->duplicate = ptap->duplicate; /* = MatDuplicate_MPIAIJ */
61   PetscFunctionReturn(0);
62 }
63 
64 #undef __FUNCT__
65 #define __FUNCT__ "MatMatMultNumeric_MPIAIJ_MPIAIJ"
66 PetscErrorCode MatMatMultNumeric_MPIAIJ_MPIAIJ(Mat A,Mat P,Mat C)
67 {
68   PetscErrorCode ierr;
69   Mat_MPIAIJ     *a  =(Mat_MPIAIJ*)A->data,*c=(Mat_MPIAIJ*)C->data;
70   Mat_SeqAIJ     *ad =(Mat_SeqAIJ*)(a->A)->data,*ao=(Mat_SeqAIJ*)(a->B)->data;
71   Mat_SeqAIJ     *cd =(Mat_SeqAIJ*)(c->A)->data,*co=(Mat_SeqAIJ*)(c->B)->data;
72   PetscInt       *adi=ad->i,*adj,*aoi=ao->i,*aoj;
73   PetscScalar    *ada,*aoa,*cda=cd->a,*coa=co->a;
74   Mat_SeqAIJ     *p_loc,*p_oth;
75   PetscInt       *pi_loc,*pj_loc,*pi_oth,*pj_oth,*pj;
76   PetscScalar    *pa_loc,*pa_oth,*pa,*apa,valtmp,*ca;
77   PetscInt       cm   =C->rmap->n,anz,pnz;
78   Mat_PtAPMPI    *ptap=c->ptap;
79   PetscInt       *api,*apj,*apJ,i,j,k,row;
80   PetscInt       cstart=C->cmap->rstart;
81   PetscInt       cdnz,conz,k0,k1;
82 
83   PetscFunctionBegin;
84   /* 1) get P_oth = ptap->P_oth  and P_loc = ptap->P_loc */
85   /*-----------------------------------------------------*/
86   /* update numerical values of P_oth and P_loc */
87   ierr = MatGetBrowsOfAoCols_MPIAIJ(A,P,MAT_REUSE_MATRIX,&ptap->startsj_s,&ptap->startsj_r,&ptap->bufa,&ptap->P_oth);CHKERRQ(ierr);
88   ierr = MatMPIAIJGetLocalMat(P,MAT_REUSE_MATRIX,&ptap->P_loc);CHKERRQ(ierr);
89 
90   /* 2) compute numeric C_loc = A_loc*P = Ad*P_loc + Ao*P_oth */
91   /*----------------------------------------------------------*/
92   /* get data from symbolic products */
93   p_loc = (Mat_SeqAIJ*)(ptap->P_loc)->data;
94   p_oth = (Mat_SeqAIJ*)(ptap->P_oth)->data;
95   pi_loc=p_loc->i; pj_loc=p_loc->j; pa_loc=p_loc->a;
96   pi_oth=p_oth->i; pj_oth=p_oth->j; pa_oth=p_oth->a;
97 
98   /* get apa for storing dense row A[i,:]*P */
99   apa = ptap->apa;
100 
101   api = ptap->api;
102   apj = ptap->apj;
103   for (i=0; i<cm; i++) {
104     /* diagonal portion of A */
105     anz = adi[i+1] - adi[i];
106     adj = ad->j + adi[i];
107     ada = ad->a + adi[i];
108     for (j=0; j<anz; j++) {
109       row = adj[j];
110       pnz = pi_loc[row+1] - pi_loc[row];
111       pj  = pj_loc + pi_loc[row];
112       pa  = pa_loc + pi_loc[row];
113 
114       /* perform dense axpy */
115       valtmp = ada[j];
116       for (k=0; k<pnz; k++) {
117         apa[pj[k]] += valtmp*pa[k];
118       }
119       ierr = PetscLogFlops(2.0*pnz);CHKERRQ(ierr);
120     }
121 
122     /* off-diagonal portion of A */
123     anz = aoi[i+1] - aoi[i];
124     aoj = ao->j + aoi[i];
125     aoa = ao->a + aoi[i];
126     for (j=0; j<anz; j++) {
127       row = aoj[j];
128       pnz = pi_oth[row+1] - pi_oth[row];
129       pj  = pj_oth + pi_oth[row];
130       pa  = pa_oth + pi_oth[row];
131 
132       /* perform dense axpy */
133       valtmp = aoa[j];
134       for (k=0; k<pnz; k++) {
135         apa[pj[k]] += valtmp*pa[k];
136       }
137       ierr = PetscLogFlops(2.0*pnz);CHKERRQ(ierr);
138     }
139 
140     /* set values in C */
141     apJ  = apj + api[i];
142     cdnz = cd->i[i+1] - cd->i[i];
143     conz = co->i[i+1] - co->i[i];
144 
145     /* 1st off-diagoanl part of C */
146     ca = coa + co->i[i];
147     k  = 0;
148     for (k0=0; k0<conz; k0++) {
149       if (apJ[k] >= cstart) break;
150       ca[k0]      = apa[apJ[k]];
151       apa[apJ[k]] = 0.0;
152       k++;
153     }
154 
155     /* diagonal part of C */
156     ca = cda + cd->i[i];
157     for (k1=0; k1<cdnz; k1++) {
158       ca[k1]      = apa[apJ[k]];
159       apa[apJ[k]] = 0.0;
160       k++;
161     }
162 
163     /* 2nd off-diagoanl part of C */
164     ca = coa + co->i[i];
165     for (; k0<conz; k0++) {
166       ca[k0]      = apa[apJ[k]];
167       apa[apJ[k]] = 0.0;
168       k++;
169     }
170   }
171   ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
172   ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
173   PetscFunctionReturn(0);
174 }
175 
176 #undef __FUNCT__
177 #define __FUNCT__ "MatMatMultSymbolic_MPIAIJ_MPIAIJ"
178 PetscErrorCode MatMatMultSymbolic_MPIAIJ_MPIAIJ(Mat A,Mat P,PetscReal fill,Mat *C)
179 {
180   PetscErrorCode     ierr;
181   MPI_Comm           comm=((PetscObject)A)->comm;
182   Mat                Cmpi;
183   Mat_PtAPMPI        *ptap;
184   PetscFreeSpaceList free_space=PETSC_NULL,current_space=PETSC_NULL;
185   Mat_MPIAIJ         *a        =(Mat_MPIAIJ*)A->data,*c;
186   Mat_SeqAIJ         *ad       =(Mat_SeqAIJ*)(a->A)->data,*ao=(Mat_SeqAIJ*)(a->B)->data,*p_loc,*p_oth;
187   PetscInt           *pi_loc,*pj_loc,*pi_oth,*pj_oth,*dnz,*onz;
188   PetscInt           *adi=ad->i,*adj=ad->j,*aoi=ao->i,*aoj=ao->j,rstart=A->rmap->rstart;
189   PetscInt           *lnk,i,pnz,row,*api,*apj,*Jptr,apnz,nspacedouble=0,j,nzi;
190   PetscInt           am=A->rmap->n,pN=P->cmap->N,pn=P->cmap->n,pm=P->rmap->n;
191   PetscBT            lnkbt;
192   PetscScalar        *apa;
193   PetscReal          afill;
194   PetscBool          scalable=PETSC_TRUE;
195   PetscInt           nlnk_max,armax,prmax;
196 
197   PetscFunctionBegin;
198   if (A->cmap->rstart != P->rmap->rstart || A->cmap->rend != P->rmap->rend) {
199     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);
200   }
201 
202   ierr = PetscObjectOptionsBegin((PetscObject)A);CHKERRQ(ierr);
203   ierr = PetscOptionsBool("-matmatmult_scalable","Use a scalable but slower C=A*B","",scalable,&scalable,PETSC_NULL);CHKERRQ(ierr);
204   if (scalable) {
205     ierr = MatMatMultSymbolic_MPIAIJ_MPIAIJ_Scalable(A,P,fill,C);CHKERRQ(ierr);
206     PetscFunctionReturn(0);
207   }
208   ierr = PetscOptionsEnd();CHKERRQ(ierr);
209 
210   /* create struct Mat_PtAPMPI and attached it to C later */
211   ierr = PetscNew(Mat_PtAPMPI,&ptap);CHKERRQ(ierr);
212 
213   /* get P_oth by taking rows of P (= non-zero cols of local A) from other processors */
214   ierr = MatGetBrowsOfAoCols_MPIAIJ(A,P,MAT_INITIAL_MATRIX,&ptap->startsj_s,&ptap->startsj_r,&ptap->bufa,&ptap->P_oth);CHKERRQ(ierr);
215 
216   /* get P_loc by taking all local rows of P */
217   ierr = MatMPIAIJGetLocalMat(P,MAT_INITIAL_MATRIX,&ptap->P_loc);CHKERRQ(ierr);
218 
219   p_loc  = (Mat_SeqAIJ*)(ptap->P_loc)->data;
220   p_oth  = (Mat_SeqAIJ*)(ptap->P_oth)->data;
221   pi_loc = p_loc->i; pj_loc = p_loc->j;
222   pi_oth = p_oth->i; pj_oth = p_oth->j;
223 
224   /* first, compute symbolic AP = A_loc*P = A_diag*P_loc + A_off*P_oth */
225   /*-------------------------------------------------------------------*/
226   ierr      = PetscMalloc((am+2)*sizeof(PetscInt),&api);CHKERRQ(ierr);
227   ptap->api = api;
228   api[0]    = 0;
229 
230   /* create and initialize a linked list */
231   armax    = ad->rmax+ao->rmax;
232   prmax    = PetscMax(p_loc->rmax,p_oth->rmax);
233   nlnk_max = armax*prmax;
234   if (!nlnk_max || nlnk_max > pN) nlnk_max = pN;
235   ierr = PetscLLCondensedCreate(nlnk_max,pN,&lnk,&lnkbt);CHKERRQ(ierr);
236 
237   /* Initial FreeSpace size is fill*(nnz(A)+nnz(P)) */
238   ierr = PetscFreeSpaceGet((PetscInt)(fill*(adi[am]+aoi[am]+pi_loc[pm])),&free_space);CHKERRQ(ierr);
239 
240   current_space = free_space;
241 
242   ierr = MatPreallocateInitialize(comm,am,pn,dnz,onz);CHKERRQ(ierr);
243   for (i=0; i<am; i++) {
244     /* diagonal portion of A */
245     nzi = adi[i+1] - adi[i];
246     for (j=0; j<nzi; j++) {
247       row  = *adj++;
248       pnz  = pi_loc[row+1] - pi_loc[row];
249       Jptr = pj_loc + pi_loc[row];
250       /* add non-zero cols of P into the sorted linked list lnk */
251       ierr = PetscLLCondensedAddSorted(pnz,Jptr,lnk,lnkbt);CHKERRQ(ierr);
252     }
253     /* off-diagonal portion of A */
254     nzi = aoi[i+1] - aoi[i];
255     for (j=0; j<nzi; j++) {
256       row  = *aoj++;
257       pnz  = pi_oth[row+1] - pi_oth[row];
258       Jptr = pj_oth + pi_oth[row];
259       ierr = PetscLLCondensedAddSorted(pnz,Jptr,lnk,lnkbt);CHKERRQ(ierr);
260     }
261 
262     apnz     = lnk[0];
263     api[i+1] = api[i] + apnz;
264 
265     /* if free space is not available, double the total space in the list */
266     if (current_space->local_remaining<apnz) {
267       ierr = PetscFreeSpaceGet(apnz+current_space->total_array_size,&current_space);CHKERRQ(ierr);
268       nspacedouble++;
269     }
270 
271     /* Copy data into free space, then initialize lnk */
272     ierr = PetscLLCondensedClean(pN,apnz,current_space->array,lnk,lnkbt);CHKERRQ(ierr);
273     ierr = MatPreallocateSet(i+rstart,apnz,current_space->array,dnz,onz);CHKERRQ(ierr);
274 
275     current_space->array           += apnz;
276     current_space->local_used      += apnz;
277     current_space->local_remaining -= apnz;
278   }
279 
280   /* Allocate space for apj, initialize apj, and */
281   /* destroy list of free space and other temporary array(s) */
282   ierr = PetscMalloc((api[am]+1)*sizeof(PetscInt),&ptap->apj);CHKERRQ(ierr);
283   apj  = ptap->apj;
284   ierr = PetscFreeSpaceContiguous(&free_space,ptap->apj);CHKERRQ(ierr);
285   ierr = PetscLLDestroy(lnk,lnkbt);CHKERRQ(ierr);
286 
287   /* malloc apa to store dense row A[i,:]*P */
288   ierr = PetscMalloc(pN*sizeof(PetscScalar),&apa);CHKERRQ(ierr);
289   ierr = PetscMemzero(apa,pN*sizeof(PetscScalar));CHKERRQ(ierr);
290 
291   ptap->apa = apa;
292 
293   /* create and assemble symbolic parallel matrix Cmpi */
294   /*----------------------------------------------------*/
295   ierr = MatCreate(comm,&Cmpi);CHKERRQ(ierr);
296   ierr = MatSetSizes(Cmpi,am,pn,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr);
297   ierr = MatSetBlockSizes(Cmpi,A->rmap->bs,P->cmap->bs);CHKERRQ(ierr);
298 
299   ierr = MatSetType(Cmpi,MATMPIAIJ);CHKERRQ(ierr);
300   ierr = MatMPIAIJSetPreallocation(Cmpi,0,dnz,0,onz);CHKERRQ(ierr);
301   ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr);
302   for (i=0; i<am; i++) {
303     row  = i + rstart;
304     apnz = api[i+1] - api[i];
305     ierr = MatSetValues(Cmpi,1,&row,apnz,apj,apa,INSERT_VALUES);CHKERRQ(ierr);
306     apj += apnz;
307   }
308   ierr = MatAssemblyBegin(Cmpi,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
309   ierr = MatAssemblyEnd(Cmpi,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
310 
311   ptap->destroy        = Cmpi->ops->destroy;
312   ptap->duplicate      = Cmpi->ops->duplicate;
313   Cmpi->ops->destroy   = MatDestroy_MPIAIJ_MatMatMult;
314   Cmpi->ops->duplicate = MatDuplicate_MPIAIJ_MatMatMult;
315 
316   /* attach the supporting struct to Cmpi for reuse */
317   c       = (Mat_MPIAIJ*)Cmpi->data;
318   c->ptap = ptap;
319 
320   *C = Cmpi;
321 
322   /* set MatInfo */
323   afill = (PetscReal)api[am]/(adi[am]+aoi[am]+pi_loc[pm]+1) + 1.e-5;
324   if (afill < 1.0) afill = 1.0;
325   Cmpi->info.mallocs           = nspacedouble;
326   Cmpi->info.fill_ratio_given  = fill;
327   Cmpi->info.fill_ratio_needed = afill;
328 
329 #if defined(PETSC_USE_INFO)
330   if (api[am]) {
331     ierr = PetscInfo3(Cmpi,"Reallocs %D; Fill ratio: given %G needed %G.\n",nspacedouble,fill,afill);CHKERRQ(ierr);
332     ierr = PetscInfo1(Cmpi,"Use MatMatMult(A,B,MatReuse,%G,&C) for best performance.;\n",afill);CHKERRQ(ierr);
333   } else {
334     ierr = PetscInfo(Cmpi,"Empty matrix product\n");CHKERRQ(ierr);
335   }
336 #endif
337   PetscFunctionReturn(0);
338 }
339 
340 #undef __FUNCT__
341 #define __FUNCT__ "MatMatMult_MPIAIJ_MPIDense"
342 PetscErrorCode MatMatMult_MPIAIJ_MPIDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
343 {
344   PetscErrorCode ierr;
345 
346   PetscFunctionBegin;
347   if (scall == MAT_INITIAL_MATRIX) {
348     ierr = MatMatMultSymbolic_MPIAIJ_MPIDense(A,B,fill,C);CHKERRQ(ierr);
349   }
350   ierr = MatMatMultNumeric_MPIAIJ_MPIDense(A,B,*C);CHKERRQ(ierr);
351   PetscFunctionReturn(0);
352 }
353 
354 typedef struct {
355   Mat         workB;
356   PetscScalar *rvalues,*svalues;
357   MPI_Request *rwaits,*swaits;
358 } MPIAIJ_MPIDense;
359 
360 #undef __FUNCT__
361 #define __FUNCT__ "MatMPIAIJ_MPIDenseDestroy"
362 PetscErrorCode MatMPIAIJ_MPIDenseDestroy(void *ctx)
363 {
364   MPIAIJ_MPIDense *contents = (MPIAIJ_MPIDense*) ctx;
365   PetscErrorCode  ierr;
366 
367   PetscFunctionBegin;
368   ierr = MatDestroy(&contents->workB);CHKERRQ(ierr);
369   ierr = PetscFree4(contents->rvalues,contents->svalues,contents->rwaits,contents->swaits);CHKERRQ(ierr);
370   ierr = PetscFree(contents);CHKERRQ(ierr);
371   PetscFunctionReturn(0);
372 }
373 
374 #undef __FUNCT__
375 #define __FUNCT__ "MatMatMultSymbolic_MPIAIJ_MPIDense"
376 PetscErrorCode MatMatMultSymbolic_MPIAIJ_MPIDense(Mat A,Mat B,PetscReal fill,Mat *C)
377 {
378   PetscErrorCode         ierr;
379   Mat_MPIAIJ             *aij = (Mat_MPIAIJ*) A->data;
380   PetscInt               nz   = aij->B->cmap->n;
381   PetscContainer         container;
382   MPIAIJ_MPIDense        *contents;
383   VecScatter             ctx   = aij->Mvctx;
384   VecScatter_MPI_General *from = (VecScatter_MPI_General*) ctx->fromdata;
385   VecScatter_MPI_General *to   = (VecScatter_MPI_General*) ctx->todata;
386   PetscInt               m     = A->rmap->n,n=B->cmap->n;
387 
388   PetscFunctionBegin;
389   ierr = MatCreate(((PetscObject)B)->comm,C);CHKERRQ(ierr);
390   ierr = MatSetSizes(*C,m,n,A->rmap->N,B->cmap->N);CHKERRQ(ierr);
391   ierr = MatSetBlockSizes(*C,A->rmap->bs,B->cmap->bs);CHKERRQ(ierr);
392   ierr = MatSetType(*C,MATMPIDENSE);CHKERRQ(ierr);
393   ierr = MatMPIDenseSetPreallocation(*C,PETSC_NULL);CHKERRQ(ierr);
394   ierr = MatAssemblyBegin(*C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
395   ierr = MatAssemblyEnd(*C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
396 
397   (*C)->ops->matmultnumeric = MatMatMultNumeric_MPIAIJ_MPIDense;
398 
399   ierr = PetscNew(MPIAIJ_MPIDense,&contents);CHKERRQ(ierr);
400   /* Create work matrix used to store off processor rows of B needed for local product */
401   ierr = MatCreateSeqDense(PETSC_COMM_SELF,nz,B->cmap->N,PETSC_NULL,&contents->workB);CHKERRQ(ierr);
402   /* Create work arrays needed */
403   ierr = PetscMalloc4(B->cmap->N*from->starts[from->n],PetscScalar,&contents->rvalues,
404                       B->cmap->N*to->starts[to->n],PetscScalar,&contents->svalues,
405                       from->n,MPI_Request,&contents->rwaits,
406                       to->n,MPI_Request,&contents->swaits);CHKERRQ(ierr);
407 
408   ierr = PetscContainerCreate(((PetscObject)A)->comm,&container);CHKERRQ(ierr);
409   ierr = PetscContainerSetPointer(container,contents);CHKERRQ(ierr);
410   ierr = PetscContainerSetUserDestroy(container,MatMPIAIJ_MPIDenseDestroy);CHKERRQ(ierr);
411   ierr = PetscObjectCompose((PetscObject)(*C),"workB",(PetscObject)container);CHKERRQ(ierr);
412   ierr = PetscContainerDestroy(&container);CHKERRQ(ierr);
413   PetscFunctionReturn(0);
414 }
415 
416 #undef __FUNCT__
417 #define __FUNCT__ "MatMPIDenseScatter"
418 /*
419     Performs an efficient scatter on the rows of B needed by this process; this is
420     a modification of the VecScatterBegin_() routines.
421 */
422 PetscErrorCode MatMPIDenseScatter(Mat A,Mat B,Mat C,Mat *outworkB)
423 {
424   Mat_MPIAIJ             *aij = (Mat_MPIAIJ*)A->data;
425   PetscErrorCode         ierr;
426   PetscScalar            *b,*w,*svalues,*rvalues;
427   VecScatter             ctx   = aij->Mvctx;
428   VecScatter_MPI_General *from = (VecScatter_MPI_General*) ctx->fromdata;
429   VecScatter_MPI_General *to   = (VecScatter_MPI_General*) ctx->todata;
430   PetscInt               i,j,k;
431   PetscInt               *sindices,*sstarts,*rindices,*rstarts;
432   PetscMPIInt            *sprocs,*rprocs,nrecvs;
433   MPI_Request            *swaits,*rwaits;
434   MPI_Comm               comm = ((PetscObject)A)->comm;
435   PetscMPIInt            tag  = ((PetscObject)ctx)->tag,ncols = B->cmap->N, nrows = aij->B->cmap->n,imdex,nrowsB = B->rmap->n;
436   MPI_Status             status;
437   MPIAIJ_MPIDense        *contents;
438   PetscContainer         container;
439   Mat                    workB;
440 
441   PetscFunctionBegin;
442   ierr = PetscObjectQuery((PetscObject)C,"workB",(PetscObject*)&container);CHKERRQ(ierr);
443   if (!container) SETERRQ(comm,PETSC_ERR_PLIB,"Container does not exist");
444   ierr = PetscContainerGetPointer(container,(void**)&contents);CHKERRQ(ierr);
445 
446   workB = *outworkB = contents->workB;
447   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);
448   sindices = to->indices;
449   sstarts  = to->starts;
450   sprocs   = to->procs;
451   swaits   = contents->swaits;
452   svalues  = contents->svalues;
453 
454   rindices = from->indices;
455   rstarts  = from->starts;
456   rprocs   = from->procs;
457   rwaits   = contents->rwaits;
458   rvalues  = contents->rvalues;
459 
460   ierr = MatDenseGetArray(B,&b);CHKERRQ(ierr);
461   ierr = MatDenseGetArray(workB,&w);CHKERRQ(ierr);
462 
463   for (i=0; i<from->n; i++) {
464     ierr = MPI_Irecv(rvalues+ncols*rstarts[i],ncols*(rstarts[i+1]-rstarts[i]),MPIU_SCALAR,rprocs[i],tag,comm,rwaits+i);CHKERRQ(ierr);
465   }
466 
467   for (i=0; i<to->n; i++) {
468     /* pack a message at a time */
469     CHKMEMQ;
470     for (j=0; j<sstarts[i+1]-sstarts[i]; j++) {
471       for (k=0; k<ncols; k++) {
472         svalues[ncols*(sstarts[i] + j) + k] = b[sindices[sstarts[i]+j] + nrowsB*k];
473       }
474     }
475     CHKMEMQ;
476     ierr = MPI_Isend(svalues+ncols*sstarts[i],ncols*(sstarts[i+1]-sstarts[i]),MPIU_SCALAR,sprocs[i],tag,comm,swaits+i);CHKERRQ(ierr);
477   }
478 
479   nrecvs = from->n;
480   while (nrecvs) {
481     ierr = MPI_Waitany(from->n,rwaits,&imdex,&status);CHKERRQ(ierr);
482     nrecvs--;
483     /* unpack a message at a time */
484     CHKMEMQ;
485     for (j=0; j<rstarts[imdex+1]-rstarts[imdex]; j++) {
486       for (k=0; k<ncols; k++) {
487         w[rindices[rstarts[imdex]+j] + nrows*k] = rvalues[ncols*(rstarts[imdex] + j) + k];
488       }
489     }
490     CHKMEMQ;
491   }
492   if (to->n) {ierr = MPI_Waitall(to->n,swaits,to->sstatus);CHKERRQ(ierr);}
493 
494   ierr = MatDenseRestoreArray(B,&b);CHKERRQ(ierr);
495   ierr = MatDenseRestoreArray(workB,&w);CHKERRQ(ierr);
496   ierr = MatAssemblyBegin(workB,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
497   ierr = MatAssemblyEnd(workB,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
498   PetscFunctionReturn(0);
499 }
500 extern PetscErrorCode MatMatMultNumericAdd_SeqAIJ_SeqDense(Mat,Mat,Mat);
501 
502 #undef __FUNCT__
503 #define __FUNCT__ "MatMatMultNumeric_MPIAIJ_MPIDense"
504 PetscErrorCode MatMatMultNumeric_MPIAIJ_MPIDense(Mat A,Mat B,Mat C)
505 {
506   PetscErrorCode ierr;
507   Mat_MPIAIJ     *aij    = (Mat_MPIAIJ*)A->data;
508   Mat_MPIDense   *bdense = (Mat_MPIDense*)B->data;
509   Mat_MPIDense   *cdense = (Mat_MPIDense*)C->data;
510   Mat            workB;
511 
512   PetscFunctionBegin;
513   /* diagonal block of A times all local rows of B*/
514   ierr = MatMatMultNumeric_SeqAIJ_SeqDense(aij->A,bdense->A,cdense->A);CHKERRQ(ierr);
515 
516   /* get off processor parts of B needed to complete the product */
517   ierr = MatMPIDenseScatter(A,B,C,&workB);CHKERRQ(ierr);
518 
519   /* off-diagonal block of A times nonlocal rows of B */
520   ierr = MatMatMultNumericAdd_SeqAIJ_SeqDense(aij->B,workB,cdense->A);CHKERRQ(ierr);
521   ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
522   ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
523   PetscFunctionReturn(0);
524 }
525 
526 #undef __FUNCT__
527 #define __FUNCT__ "MatMatMultNumeric_MPIAIJ_MPIAIJ_Scalable"
528 PetscErrorCode MatMatMultNumeric_MPIAIJ_MPIAIJ_Scalable(Mat A,Mat P,Mat C)
529 {
530   PetscErrorCode ierr;
531   Mat_MPIAIJ     *a   = (Mat_MPIAIJ*)A->data,*c=(Mat_MPIAIJ*)C->data;
532   Mat_SeqAIJ     *ad  = (Mat_SeqAIJ*)(a->A)->data,*ao=(Mat_SeqAIJ*)(a->B)->data;
533   Mat_SeqAIJ     *cd  = (Mat_SeqAIJ*)(c->A)->data,*co=(Mat_SeqAIJ*)(c->B)->data;
534   PetscInt       *adi = ad->i,*adj,*aoi=ao->i,*aoj;
535   PetscScalar    *ada,*aoa,*cda=cd->a,*coa=co->a;
536   Mat_SeqAIJ     *p_loc,*p_oth;
537   PetscInt       *pi_loc,*pj_loc,*pi_oth,*pj_oth,*pj;
538   PetscScalar    *pa_loc,*pa_oth,*pa,valtmp,*ca;
539   PetscInt       cm          = C->rmap->n,anz,pnz;
540   Mat_PtAPMPI    *ptap       = c->ptap;
541   PetscScalar    *apa_sparse = ptap->apa;
542   PetscInt       *api,*apj,*apJ,i,j,k,row;
543   PetscInt       cstart = C->cmap->rstart;
544   PetscInt       cdnz,conz,k0,k1,nextp;
545 
546   PetscFunctionBegin;
547   /* 1) get P_oth = ptap->P_oth  and P_loc = ptap->P_loc */
548   /*-----------------------------------------------------*/
549   /* update numerical values of P_oth and P_loc */
550   ierr = MatGetBrowsOfAoCols_MPIAIJ(A,P,MAT_REUSE_MATRIX,&ptap->startsj_s,&ptap->startsj_r,&ptap->bufa,&ptap->P_oth);CHKERRQ(ierr);
551   ierr = MatMPIAIJGetLocalMat(P,MAT_REUSE_MATRIX,&ptap->P_loc);CHKERRQ(ierr);
552 
553   /* 2) compute numeric C_loc = A_loc*P = Ad*P_loc + Ao*P_oth */
554   /*----------------------------------------------------------*/
555   /* get data from symbolic products */
556   p_loc = (Mat_SeqAIJ*)(ptap->P_loc)->data;
557   p_oth = (Mat_SeqAIJ*)(ptap->P_oth)->data;
558   pi_loc=p_loc->i; pj_loc=p_loc->j; pa_loc=p_loc->a;
559   pi_oth=p_oth->i; pj_oth=p_oth->j; pa_oth=p_oth->a;
560 
561   api = ptap->api;
562   apj = ptap->apj;
563   for (i=0; i<cm; i++) {
564     apJ = apj + api[i];
565 
566     /* diagonal portion of A */
567     anz = adi[i+1] - adi[i];
568     adj = ad->j + adi[i];
569     ada = ad->a + adi[i];
570     for (j=0; j<anz; j++) {
571       row = adj[j];
572       pnz = pi_loc[row+1] - pi_loc[row];
573       pj  = pj_loc + pi_loc[row];
574       pa  = pa_loc + pi_loc[row];
575       /* perform sparse axpy */
576       valtmp = ada[j];
577       nextp  = 0;
578       for (k=0; nextp<pnz; k++) {
579         if (apJ[k] == pj[nextp]) { /* column of AP == column of P */
580           apa_sparse[k] += valtmp*pa[nextp++];
581         }
582       }
583       ierr = PetscLogFlops(2.0*pnz);CHKERRQ(ierr);
584     }
585 
586     /* off-diagonal portion of A */
587     anz = aoi[i+1] - aoi[i];
588     aoj = ao->j + aoi[i];
589     aoa = ao->a + aoi[i];
590     for (j=0; j<anz; j++) {
591       row = aoj[j];
592       pnz = pi_oth[row+1] - pi_oth[row];
593       pj  = pj_oth + pi_oth[row];
594       pa  = pa_oth + pi_oth[row];
595       /* perform sparse axpy */
596       valtmp = aoa[j];
597       nextp  = 0;
598       for (k=0; nextp<pnz; k++) {
599         if (apJ[k] == pj[nextp]) { /* column of AP == column of P */
600           apa_sparse[k] += valtmp*pa[nextp++];
601         }
602       }
603       ierr = PetscLogFlops(2.0*pnz);CHKERRQ(ierr);
604     }
605 
606     /* set values in C */
607     cdnz = cd->i[i+1] - cd->i[i];
608     conz = co->i[i+1] - co->i[i];
609 
610     /* 1st off-diagoanl part of C */
611     ca = coa + co->i[i];
612     k  = 0;
613     for (k0=0; k0<conz; k0++) {
614       if (apJ[k] >= cstart) break;
615       ca[k0]        = apa_sparse[k];
616       apa_sparse[k] = 0.0;
617       k++;
618     }
619 
620     /* diagonal part of C */
621     ca = cda + cd->i[i];
622     for (k1=0; k1<cdnz; k1++) {
623       ca[k1]        = apa_sparse[k];
624       apa_sparse[k] = 0.0;
625       k++;
626     }
627 
628     /* 2nd off-diagoanl part of C */
629     ca = coa + co->i[i];
630     for (; k0<conz; k0++) {
631       ca[k0]        = apa_sparse[k];
632       apa_sparse[k] = 0.0;
633       k++;
634     }
635   }
636   ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
637   ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
638   PetscFunctionReturn(0);
639 }
640 
641 /* same as MatMatMultSymbolic_MPIAIJ_MPIAIJ(), except using LLCondensed to avoid O(BN) memory requirement */
642 #undef __FUNCT__
643 #define __FUNCT__ "MatMatMultSymbolic_MPIAIJ_MPIAIJ_Scalable"
644 PetscErrorCode MatMatMultSymbolic_MPIAIJ_MPIAIJ_Scalable(Mat A,Mat P,PetscReal fill,Mat *C)
645 {
646   PetscErrorCode     ierr;
647   MPI_Comm           comm=((PetscObject)A)->comm;
648   Mat                Cmpi;
649   Mat_PtAPMPI        *ptap;
650   PetscFreeSpaceList free_space = PETSC_NULL,current_space=PETSC_NULL;
651   Mat_MPIAIJ         *a         = (Mat_MPIAIJ*)A->data,*c;
652   Mat_SeqAIJ         *ad        = (Mat_SeqAIJ*)(a->A)->data,*ao=(Mat_SeqAIJ*)(a->B)->data,*p_loc,*p_oth;
653   PetscInt           *pi_loc,*pj_loc,*pi_oth,*pj_oth,*dnz,*onz;
654   PetscInt           *adi=ad->i,*adj=ad->j,*aoi=ao->i,*aoj=ao->j,rstart=A->rmap->rstart;
655   PetscInt           i,pnz,row,*api,*apj,*Jptr,apnz,nspacedouble=0,j,nzi,*lnk,apnz_max=0;
656   PetscInt           am=A->rmap->n,pN=P->cmap->N,pn=P->cmap->n,pm=P->rmap->n;
657   PetscInt           nlnk_max,armax,prmax;
658   PetscReal          afill;
659   PetscScalar        *apa;
660 
661   PetscFunctionBegin;
662   /* create struct Mat_PtAPMPI and attached it to C later */
663   ierr = PetscNew(Mat_PtAPMPI,&ptap);CHKERRQ(ierr);
664 
665   /* get P_oth by taking rows of P (= non-zero cols of local A) from other processors */
666   ierr = MatGetBrowsOfAoCols_MPIAIJ(A,P,MAT_INITIAL_MATRIX,&ptap->startsj_s,&ptap->startsj_r,&ptap->bufa,&ptap->P_oth);CHKERRQ(ierr);
667 
668   /* get P_loc by taking all local rows of P */
669   ierr = MatMPIAIJGetLocalMat(P,MAT_INITIAL_MATRIX,&ptap->P_loc);CHKERRQ(ierr);
670 
671   p_loc  = (Mat_SeqAIJ*)(ptap->P_loc)->data;
672   p_oth  = (Mat_SeqAIJ*)(ptap->P_oth)->data;
673   pi_loc = p_loc->i; pj_loc = p_loc->j;
674   pi_oth = p_oth->i; pj_oth = p_oth->j;
675 
676   /* first, compute symbolic AP = A_loc*P = A_diag*P_loc + A_off*P_oth */
677   /*-------------------------------------------------------------------*/
678   ierr      = PetscMalloc((am+2)*sizeof(PetscInt),&api);CHKERRQ(ierr);
679   ptap->api = api;
680   api[0]    = 0;
681 
682   /* create and initialize a linked list */
683   armax    = ad->rmax+ao->rmax;
684   prmax    = PetscMax(p_loc->rmax,p_oth->rmax);
685   nlnk_max = armax*prmax;
686   if (!nlnk_max || nlnk_max > pN) nlnk_max = pN;
687   ierr = PetscLLCondensedCreate_Scalable(nlnk_max,&lnk);CHKERRQ(ierr);
688 
689   /* Initial FreeSpace size is fill*(nnz(A)+nnz(P)) */
690   ierr = PetscFreeSpaceGet((PetscInt)(fill*(adi[am]+aoi[am]+pi_loc[pm])),&free_space);CHKERRQ(ierr);
691 
692   current_space = free_space;
693 
694   ierr = MatPreallocateInitialize(comm,am,pn,dnz,onz);CHKERRQ(ierr);
695   for (i=0; i<am; i++) {
696     /* diagonal portion of A */
697     nzi = adi[i+1] - adi[i];
698     for (j=0; j<nzi; j++) {
699       row  = *adj++;
700       pnz  = pi_loc[row+1] - pi_loc[row];
701       Jptr = pj_loc + pi_loc[row];
702       /* add non-zero cols of P into the sorted linked list lnk */
703       ierr = PetscLLCondensedAddSorted_Scalable(pnz,Jptr,lnk);CHKERRQ(ierr);
704     }
705     /* off-diagonal portion of A */
706     nzi = aoi[i+1] - aoi[i];
707     for (j=0; j<nzi; j++) {
708       row  = *aoj++;
709       pnz  = pi_oth[row+1] - pi_oth[row];
710       Jptr = pj_oth + pi_oth[row];
711       ierr = PetscLLCondensedAddSorted_Scalable(pnz,Jptr,lnk);CHKERRQ(ierr);
712     }
713 
714     apnz     = *lnk;
715     api[i+1] = api[i] + apnz;
716     if (apnz > apnz_max) apnz_max = apnz;
717 
718     /* if free space is not available, double the total space in the list */
719     if (current_space->local_remaining<apnz) {
720       ierr = PetscFreeSpaceGet(apnz+current_space->total_array_size,&current_space);CHKERRQ(ierr);
721       nspacedouble++;
722     }
723 
724     /* Copy data into free space, then initialize lnk */
725     ierr = PetscLLCondensedClean_Scalable(apnz,current_space->array,lnk);CHKERRQ(ierr);
726     ierr = MatPreallocateSet(i+rstart,apnz,current_space->array,dnz,onz);CHKERRQ(ierr);
727 
728     current_space->array           += apnz;
729     current_space->local_used      += apnz;
730     current_space->local_remaining -= apnz;
731   }
732 
733   /* Allocate space for apj, initialize apj, and */
734   /* destroy list of free space and other temporary array(s) */
735   ierr = PetscMalloc((api[am]+1)*sizeof(PetscInt),&ptap->apj);CHKERRQ(ierr);
736   apj  = ptap->apj;
737   ierr = PetscFreeSpaceContiguous(&free_space,ptap->apj);CHKERRQ(ierr);
738   ierr = PetscLLCondensedDestroy_Scalable(lnk);CHKERRQ(ierr);
739 
740   /* create and assemble symbolic parallel matrix Cmpi */
741   /*----------------------------------------------------*/
742   ierr = MatCreate(comm,&Cmpi);CHKERRQ(ierr);
743   ierr = MatSetSizes(Cmpi,am,pn,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr);
744   ierr = MatSetBlockSizes(Cmpi,A->rmap->bs,P->cmap->bs);CHKERRQ(ierr);
745   ierr = MatSetType(Cmpi,MATMPIAIJ);CHKERRQ(ierr);
746   ierr = MatMPIAIJSetPreallocation(Cmpi,0,dnz,0,onz);CHKERRQ(ierr);
747   ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr);
748 
749   /* malloc apa for assembly Cmpi */
750   ierr = PetscMalloc(apnz_max*sizeof(PetscScalar),&apa);CHKERRQ(ierr);
751   ierr = PetscMemzero(apa,apnz_max*sizeof(PetscScalar));CHKERRQ(ierr);
752 
753   ptap->apa = apa;
754   for (i=0; i<am; i++) {
755     row  = i + rstart;
756     apnz = api[i+1] - api[i];
757     ierr = MatSetValues(Cmpi,1,&row,apnz,apj,apa,INSERT_VALUES);CHKERRQ(ierr);
758     apj += apnz;
759   }
760   ierr = MatAssemblyBegin(Cmpi,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
761   ierr = MatAssemblyEnd(Cmpi,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
762 
763   ptap->destroy             = Cmpi->ops->destroy;
764   ptap->duplicate           = Cmpi->ops->duplicate;
765   Cmpi->ops->matmultnumeric = MatMatMultNumeric_MPIAIJ_MPIAIJ_Scalable;
766   Cmpi->ops->destroy        = MatDestroy_MPIAIJ_MatMatMult;
767   Cmpi->ops->duplicate      = MatDuplicate_MPIAIJ_MatMatMult;
768 
769   /* attach the supporting struct to Cmpi for reuse */
770   c       = (Mat_MPIAIJ*)Cmpi->data;
771   c->ptap = ptap;
772 
773   *C = Cmpi;
774 
775   /* set MatInfo */
776   afill = (PetscReal)api[am]/(adi[am]+aoi[am]+pi_loc[pm]+1) + 1.e-5;
777   if (afill < 1.0) afill = 1.0;
778   Cmpi->info.mallocs           = nspacedouble;
779   Cmpi->info.fill_ratio_given  = fill;
780   Cmpi->info.fill_ratio_needed = afill;
781 
782 #if defined(PETSC_USE_INFO)
783   if (api[am]) {
784     ierr = PetscInfo3(Cmpi,"Reallocs %D; Fill ratio: given %G needed %G.\n",nspacedouble,fill,afill);CHKERRQ(ierr);
785     ierr = PetscInfo1(Cmpi,"Use MatMatMult(A,B,MatReuse,%G,&C) for best performance.;\n",afill);CHKERRQ(ierr);
786   } else {
787     ierr = PetscInfo(Cmpi,"Empty matrix product\n");CHKERRQ(ierr);
788   }
789 #endif
790   PetscFunctionReturn(0);
791 }
792 
793 /*-------------------------------------------------------------------------*/
794 #undef __FUNCT__
795 #define __FUNCT__ "MatTransposeMatMult_MPIAIJ_MPIAIJ"
796 PetscErrorCode MatTransposeMatMult_MPIAIJ_MPIAIJ(Mat P,Mat A,MatReuse scall,PetscReal fill,Mat *C)
797 {
798   PetscErrorCode ierr;
799   PetscBool      scalable=PETSC_TRUE,viamatmatmult=PETSC_FALSE;
800 
801   PetscFunctionBegin;
802   ierr = PetscOptionsBool("-mattransposematmult_viamatmatmult","Use R=Pt and C=R*A","",viamatmatmult,&viamatmatmult,PETSC_NULL);CHKERRQ(ierr);
803   if (viamatmatmult) {
804     Mat         Pt;
805     Mat_PtAPMPI *ptap;
806     Mat_MPIAIJ  *c;
807     if (scall == MAT_INITIAL_MATRIX) {
808       ierr = MatTranspose(P,MAT_INITIAL_MATRIX,&Pt);CHKERRQ(ierr);
809       ierr = MatMatMult(Pt,A,MAT_INITIAL_MATRIX,fill,C);CHKERRQ(ierr);
810 
811       c        = (Mat_MPIAIJ*)(*C)->data;
812       ptap     = c->ptap;
813       ptap->Pt = Pt;
814     } else if (scall == MAT_REUSE_MATRIX) {
815       c    = (Mat_MPIAIJ*)(*C)->data;
816       ptap = c->ptap;
817       Pt   = ptap->Pt;
818       ierr = MatTranspose(P,scall,&Pt);CHKERRQ(ierr);
819       ierr = MatMatMult(Pt,A,scall,fill,C);CHKERRQ(ierr);
820     } else SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not supported");
821     PetscFunctionReturn(0);
822   }
823 
824   if (scall == MAT_INITIAL_MATRIX) {
825     ierr = PetscObjectOptionsBegin((PetscObject)A);CHKERRQ(ierr);
826     ierr = PetscOptionsBool("-mattransposematmult_scalable","Use a scalable but slower C=Pt*A","",scalable,&scalable,PETSC_NULL);CHKERRQ(ierr);
827     if  (scalable) {
828       ierr = MatTransposeMatMultSymbolic_MPIAIJ_MPIAIJ_Scalable(P,A,fill,C);CHKERRQ(ierr);
829     } else {
830       ierr = MatTransposeMatMultSymbolic_MPIAIJ_MPIAIJ(P,A,fill,C);CHKERRQ(ierr);
831     }
832     ierr = PetscOptionsEnd();CHKERRQ(ierr);
833   }
834   ierr = (*(*C)->ops->mattransposemultnumeric)(P,A,*C);CHKERRQ(ierr);
835   PetscFunctionReturn(0);
836 }
837 
838 #undef __FUNCT__
839 #define __FUNCT__ "MatTransposeMatMultNumeric_MPIAIJ_MPIAIJ"
840 PetscErrorCode MatTransposeMatMultNumeric_MPIAIJ_MPIAIJ(Mat P,Mat A,Mat C)
841 {
842   PetscErrorCode      ierr;
843   Mat_Merge_SeqsToMPI *merge;
844   Mat_MPIAIJ          *p =(Mat_MPIAIJ*)P->data,*c=(Mat_MPIAIJ*)C->data;
845   Mat_SeqAIJ          *pd=(Mat_SeqAIJ*)(p->A)->data,*po=(Mat_SeqAIJ*)(p->B)->data;
846   Mat_PtAPMPI         *ptap;
847   PetscInt            *adj,*aJ;
848   PetscInt            i,j,k,anz,pnz,row,*cj;
849   MatScalar           *ada,*aval,*ca,valtmp;
850   PetscInt            am  =A->rmap->n,cm=C->rmap->n,pon=(p->B)->cmap->n;
851   MPI_Comm            comm=((PetscObject)C)->comm;
852   PetscMPIInt         size,rank,taga,*len_s;
853   PetscInt            *owners,proc,nrows,**buf_ri_k,**nextrow,**nextci;
854   PetscInt            **buf_ri,**buf_rj;
855   PetscInt            cnz=0,*bj_i,*bi,*bj,bnz,nextcj;  /* bi,bj,ba: local array of C(mpi mat) */
856   MPI_Request         *s_waits,*r_waits;
857   MPI_Status          *status;
858   MatScalar           **abuf_r,*ba_i,*pA,*coa,*ba;
859   PetscInt            *ai,*aj,*coi,*coj;
860   PetscInt            *poJ,*pdJ;
861   Mat                 A_loc;
862   Mat_SeqAIJ          *a_loc;
863 
864   PetscFunctionBegin;
865   ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
866   ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
867 
868   ptap  = c->ptap;
869   merge = ptap->merge;
870 
871   /* 2) compute numeric C_seq = P_loc^T*A_loc*P - dominating part */
872   /*--------------------------------------------------------------*/
873   /* get data from symbolic products */
874   coi  = merge->coi; coj = merge->coj;
875   ierr = PetscMalloc((coi[pon]+1)*sizeof(MatScalar),&coa);CHKERRQ(ierr);
876   ierr = PetscMemzero(coa,coi[pon]*sizeof(MatScalar));CHKERRQ(ierr);
877 
878   bi     = merge->bi; bj = merge->bj;
879   owners = merge->rowmap->range;
880   ierr   = PetscMalloc((bi[cm]+1)*sizeof(MatScalar),&ba);CHKERRQ(ierr);
881   ierr   = PetscMemzero(ba,bi[cm]*sizeof(MatScalar));CHKERRQ(ierr);
882 
883   /* get A_loc by taking all local rows of A */
884   A_loc = ptap->A_loc;
885   ierr  = MatMPIAIJGetLocalMat(A,MAT_REUSE_MATRIX,&A_loc);CHKERRQ(ierr);
886   a_loc = (Mat_SeqAIJ*)(A_loc)->data;
887   ai    = a_loc->i;
888   aj    = a_loc->j;
889 
890   ierr = PetscMalloc((A->cmap->N)*sizeof(PetscScalar),&aval);CHKERRQ(ierr); /* non-scalable!!! */
891   ierr = PetscMemzero(aval,A->cmap->N*sizeof(PetscScalar));CHKERRQ(ierr);
892 
893   for (i=0; i<am; i++) {
894     /* 2-a) put A[i,:] to dense array aval */
895     anz = ai[i+1] - ai[i];
896     adj = aj + ai[i];
897     ada = a_loc->a + ai[i];
898     for (j=0; j<anz; j++) {
899       aval[adj[j]] = ada[j];
900     }
901 
902     /* 2-b) Compute Cseq = P_loc[i,:]^T*A[i,:] using outer product */
903     /*--------------------------------------------------------------*/
904     /* put the value into Co=(p->B)^T*A (off-diagonal part, send to others) */
905     pnz = po->i[i+1] - po->i[i];
906     poJ = po->j + po->i[i];
907     pA  = po->a + po->i[i];
908     for (j=0; j<pnz; j++) {
909       row = poJ[j];
910       cnz = coi[row+1] - coi[row];
911       cj  = coj + coi[row];
912       ca  = coa + coi[row];
913       /* perform dense axpy */
914       valtmp = pA[j];
915       for (k=0; k<cnz; k++) {
916         ca[k] += valtmp*aval[cj[k]];
917       }
918       ierr = PetscLogFlops(2.0*cnz);CHKERRQ(ierr);
919     }
920 
921     /* put the value into Cd (diagonal part) */
922     pnz = pd->i[i+1] - pd->i[i];
923     pdJ = pd->j + pd->i[i];
924     pA  = pd->a + pd->i[i];
925     for (j=0; j<pnz; j++) {
926       row = pdJ[j];
927       cnz = bi[row+1] - bi[row];
928       cj  = bj + bi[row];
929       ca  = ba + bi[row];
930       /* perform dense axpy */
931       valtmp = pA[j];
932       for (k=0; k<cnz; k++) {
933         ca[k] += valtmp*aval[cj[k]];
934       }
935       ierr = PetscLogFlops(2.0*cnz);CHKERRQ(ierr);
936     }
937 
938     /* zero the current row of Pt*A */
939     aJ = aj + ai[i];
940     for (k=0; k<anz; k++) aval[aJ[k]] = 0.0;
941   }
942 
943   /* 3) send and recv matrix values coa */
944   /*------------------------------------*/
945   buf_ri = merge->buf_ri;
946   buf_rj = merge->buf_rj;
947   len_s  = merge->len_s;
948   ierr   = PetscCommGetNewTag(comm,&taga);CHKERRQ(ierr);
949   ierr   = PetscPostIrecvScalar(comm,taga,merge->nrecv,merge->id_r,merge->len_r,&abuf_r,&r_waits);CHKERRQ(ierr);
950 
951   ierr = PetscMalloc2(merge->nsend+1,MPI_Request,&s_waits,size,MPI_Status,&status);CHKERRQ(ierr);
952   for (proc=0,k=0; proc<size; proc++) {
953     if (!len_s[proc]) continue;
954     i    = merge->owners_co[proc];
955     ierr = MPI_Isend(coa+coi[i],len_s[proc],MPIU_MATSCALAR,proc,taga,comm,s_waits+k);CHKERRQ(ierr);
956     k++;
957   }
958   if (merge->nrecv) {ierr = MPI_Waitall(merge->nrecv,r_waits,status);CHKERRQ(ierr);}
959   if (merge->nsend) {ierr = MPI_Waitall(merge->nsend,s_waits,status);CHKERRQ(ierr);}
960 
961   ierr = PetscFree2(s_waits,status);CHKERRQ(ierr);
962   ierr = PetscFree(r_waits);CHKERRQ(ierr);
963   ierr = PetscFree(coa);CHKERRQ(ierr);
964 
965   /* 4) insert local Cseq and received values into Cmpi */
966   /*----------------------------------------------------*/
967   ierr = PetscMalloc3(merge->nrecv,PetscInt**,&buf_ri_k,merge->nrecv,PetscInt*,&nextrow,merge->nrecv,PetscInt*,&nextci);CHKERRQ(ierr);
968   for (k=0; k<merge->nrecv; k++) {
969     buf_ri_k[k] = buf_ri[k]; /* beginning of k-th recved i-structure */
970     nrows       = *(buf_ri_k[k]);
971     nextrow[k]  = buf_ri_k[k]+1;  /* next row number of k-th recved i-structure */
972     nextci[k]   = buf_ri_k[k] + (nrows + 1); /* poins to the next i-structure of k-th recved i-structure  */
973   }
974 
975   for (i=0; i<cm; i++) {
976     row  = owners[rank] + i; /* global row index of C_seq */
977     bj_i = bj + bi[i];  /* col indices of the i-th row of C */
978     ba_i = ba + bi[i];
979     bnz  = bi[i+1] - bi[i];
980     /* add received vals into ba */
981     for (k=0; k<merge->nrecv; k++) { /* k-th received message */
982       /* i-th row */
983       if (i == *nextrow[k]) {
984         cnz    = *(nextci[k]+1) - *nextci[k];
985         cj     = buf_rj[k] + *(nextci[k]);
986         ca     = abuf_r[k] + *(nextci[k]);
987         nextcj = 0;
988         for (j=0; nextcj<cnz; j++) {
989           if (bj_i[j] == cj[nextcj]) { /* bcol == ccol */
990             ba_i[j] += ca[nextcj++];
991           }
992         }
993         nextrow[k]++; nextci[k]++;
994         ierr = PetscLogFlops(2.0*cnz);CHKERRQ(ierr);
995       }
996     }
997     ierr = MatSetValues(C,1,&row,bnz,bj_i,ba_i,INSERT_VALUES);CHKERRQ(ierr);
998   }
999   ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1000   ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1001 
1002   ierr = PetscFree(ba);CHKERRQ(ierr);
1003   ierr = PetscFree(abuf_r[0]);CHKERRQ(ierr);
1004   ierr = PetscFree(abuf_r);CHKERRQ(ierr);
1005   ierr = PetscFree3(buf_ri_k,nextrow,nextci);CHKERRQ(ierr);
1006   ierr = PetscFree(aval);CHKERRQ(ierr);
1007   PetscFunctionReturn(0);
1008 }
1009 
1010 /* This routine is modified from MatPtAPSymbolic_MPIAIJ_MPIAIJ() */
1011 #undef __FUNCT__
1012 #define __FUNCT__ "MatTransposeMatMultSymbolic_MPIAIJ_MPIAIJ"
1013 PetscErrorCode MatTransposeMatMultSymbolic_MPIAIJ_MPIAIJ(Mat P,Mat A,PetscReal fill,Mat *C)
1014 {
1015   PetscErrorCode      ierr;
1016   Mat                 Cmpi,A_loc,POt,PDt;
1017   Mat_PtAPMPI         *ptap;
1018   PetscFreeSpaceList  free_space=PETSC_NULL,current_space=PETSC_NULL;
1019   Mat_MPIAIJ          *p        =(Mat_MPIAIJ*)P->data,*c;
1020   PetscInt            *pdti,*pdtj,*poti,*potj,*ptJ;
1021   PetscInt            nnz;
1022   PetscInt            *lnk,*owners_co,*coi,*coj,i,k,pnz,row;
1023   PetscInt            am=A->rmap->n,pn=P->cmap->n;
1024   PetscBT             lnkbt;
1025   MPI_Comm            comm=((PetscObject)A)->comm;
1026   PetscMPIInt         size,rank,tagi,tagj,*len_si,*len_s,*len_ri;
1027   PetscInt            **buf_rj,**buf_ri,**buf_ri_k;
1028   PetscInt            len,proc,*dnz,*onz,*owners;
1029   PetscInt            nzi,*bi,*bj;
1030   PetscInt            nrows,*buf_s,*buf_si,*buf_si_i,**nextrow,**nextci;
1031   MPI_Request         *swaits,*rwaits;
1032   MPI_Status          *sstatus,rstatus;
1033   Mat_Merge_SeqsToMPI *merge;
1034   PetscInt            *ai,*aj,*Jptr,anz,*prmap=p->garray,pon,nspacedouble=0,j;
1035   PetscReal           afill  =1.0,afill_tmp;
1036   PetscInt            rstart = P->cmap->rstart,rmax,aN=A->cmap->N,Crmax;
1037   PetscScalar         *vals;
1038   Mat_SeqAIJ          *a_loc, *pdt,*pot;
1039 
1040   PetscFunctionBegin;
1041   /* check if matrix local sizes are compatible */
1042   if (A->rmap->rstart != P->rmap->rstart || A->rmap->rend != P->rmap->rend) {
1043     SETERRQ4(comm,PETSC_ERR_ARG_SIZ,"Matrix local dimensions are incompatible, A (%D, %D) != P (%D,%D)",A->rmap->rstart,A->rmap->rend,P->rmap->rstart,P->rmap->rend);
1044   }
1045 
1046   ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
1047   ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
1048 
1049   /* create struct Mat_PtAPMPI and attached it to C later */
1050   ierr = PetscNew(Mat_PtAPMPI,&ptap);CHKERRQ(ierr);
1051 
1052   /* get A_loc by taking all local rows of A */
1053   ierr = MatMPIAIJGetLocalMat(A,MAT_INITIAL_MATRIX,&A_loc);CHKERRQ(ierr);
1054 
1055   ptap->A_loc = A_loc;
1056 
1057   a_loc = (Mat_SeqAIJ*)(A_loc)->data;
1058   ai    = a_loc->i;
1059   aj    = a_loc->j;
1060 
1061   /* determine symbolic Co=(p->B)^T*A - send to others */
1062   /*----------------------------------------------------*/
1063   ierr = MatTransposeSymbolic_SeqAIJ(p->A,&PDt);CHKERRQ(ierr);
1064   pdt  = (Mat_SeqAIJ*)PDt->data;
1065   pdti = pdt->i; pdtj = pdt->j;
1066 
1067   ierr = MatTransposeSymbolic_SeqAIJ(p->B,&POt);CHKERRQ(ierr);
1068   pot  = (Mat_SeqAIJ*)POt->data;
1069   poti = pot->i; potj = pot->j;
1070 
1071   /* then, compute symbolic Co = (p->B)^T*A */
1072   pon    = (p->B)->cmap->n; /* total num of rows to be sent to other processors >= (num of nonzero rows of C_seq) - pn */
1073   ierr   = PetscMalloc((pon+1)*sizeof(PetscInt),&coi);CHKERRQ(ierr);
1074   coi[0] = 0;
1075 
1076   /* set initial free space to be fill*(nnz(p->B) + nnz(A)) */
1077   nnz           = fill*(poti[pon] + ai[am]);
1078   ierr          = PetscFreeSpaceGet(nnz,&free_space);CHKERRQ(ierr);
1079   current_space = free_space;
1080 
1081   /* create and initialize a linked list */
1082   i     = PetscMax(pdt->rmax,pot->rmax);
1083   Crmax = i*a_loc->rmax*size;
1084   if (!Crmax || Crmax > aN) Crmax = aN;
1085   ierr = PetscLLCondensedCreate(Crmax,aN,&lnk,&lnkbt);CHKERRQ(ierr);
1086 
1087   for (i=0; i<pon; i++) {
1088     pnz = poti[i+1] - poti[i];
1089     ptJ = potj + poti[i];
1090     for (j=0; j<pnz; j++) {
1091       row  = ptJ[j]; /* row of A_loc == col of Pot */
1092       anz  = ai[row+1] - ai[row];
1093       Jptr = aj + ai[row];
1094       /* add non-zero cols of AP into the sorted linked list lnk */
1095       ierr = PetscLLCondensedAddSorted(anz,Jptr,lnk,lnkbt);CHKERRQ(ierr);
1096     }
1097     nnz = lnk[0];
1098 
1099     /* If free space is not available, double the total space in the list */
1100     if (current_space->local_remaining<nnz) {
1101       ierr = PetscFreeSpaceGet(nnz+current_space->total_array_size,&current_space);CHKERRQ(ierr);
1102       nspacedouble++;
1103     }
1104 
1105     /* Copy data into free space, and zero out denserows */
1106     ierr = PetscLLCondensedClean(aN,nnz,current_space->array,lnk,lnkbt);CHKERRQ(ierr);
1107 
1108     current_space->array           += nnz;
1109     current_space->local_used      += nnz;
1110     current_space->local_remaining -= nnz;
1111 
1112     coi[i+1] = coi[i] + nnz;
1113   }
1114 
1115   ierr = PetscMalloc((coi[pon]+1)*sizeof(PetscInt),&coj);CHKERRQ(ierr);
1116   ierr = PetscFreeSpaceContiguous(&free_space,coj);CHKERRQ(ierr);
1117 
1118   afill_tmp = (PetscReal)coi[pon]/(poti[pon] + ai[am]+1);
1119   if (afill_tmp > afill) afill = afill_tmp;
1120 
1121   /* send j-array (coj) of Co to other processors */
1122   /*----------------------------------------------*/
1123   /* determine row ownership */
1124   ierr = PetscNew(Mat_Merge_SeqsToMPI,&merge);CHKERRQ(ierr);
1125   ierr = PetscLayoutCreate(comm,&merge->rowmap);CHKERRQ(ierr);
1126 
1127   merge->rowmap->n  = pn;
1128   merge->rowmap->bs = 1;
1129 
1130   ierr   = PetscLayoutSetUp(merge->rowmap);CHKERRQ(ierr);
1131   owners = merge->rowmap->range;
1132 
1133   /* determine the number of messages to send, their lengths */
1134   ierr = PetscMalloc(size*sizeof(PetscMPIInt),&len_si);CHKERRQ(ierr);
1135   ierr = PetscMemzero(len_si,size*sizeof(PetscMPIInt));CHKERRQ(ierr);
1136   ierr = PetscMalloc(size*sizeof(PetscMPIInt),&merge->len_s);CHKERRQ(ierr);
1137 
1138   len_s        = merge->len_s;
1139   merge->nsend = 0;
1140 
1141   ierr = PetscMalloc((size+2)*sizeof(PetscInt),&owners_co);CHKERRQ(ierr);
1142   ierr = PetscMemzero(len_s,size*sizeof(PetscMPIInt));CHKERRQ(ierr);
1143 
1144   proc = 0;
1145   for (i=0; i<pon; i++) {
1146     while (prmap[i] >= owners[proc+1]) proc++;
1147     len_si[proc]++;  /* num of rows in Co to be sent to [proc] */
1148     len_s[proc] += coi[i+1] - coi[i];
1149   }
1150 
1151   len          = 0; /* max length of buf_si[] */
1152   owners_co[0] = 0;
1153   for (proc=0; proc<size; proc++) {
1154     owners_co[proc+1] = owners_co[proc] + len_si[proc];
1155     if (len_si[proc]) {
1156       merge->nsend++;
1157       len_si[proc] = 2*(len_si[proc] + 1);
1158       len         += len_si[proc];
1159     }
1160   }
1161 
1162   /* determine the number and length of messages to receive for coi and coj  */
1163   ierr = PetscGatherNumberOfMessages(comm,PETSC_NULL,len_s,&merge->nrecv);CHKERRQ(ierr);
1164   ierr = PetscGatherMessageLengths2(comm,merge->nsend,merge->nrecv,len_s,len_si,&merge->id_r,&merge->len_r,&len_ri);CHKERRQ(ierr);
1165 
1166   /* post the Irecv and Isend of coj */
1167   ierr = PetscCommGetNewTag(comm,&tagj);CHKERRQ(ierr);
1168   ierr = PetscPostIrecvInt(comm,tagj,merge->nrecv,merge->id_r,merge->len_r,&buf_rj,&rwaits);CHKERRQ(ierr);
1169   ierr = PetscMalloc((merge->nsend+1)*sizeof(MPI_Request),&swaits);CHKERRQ(ierr);
1170   for (proc=0, k=0; proc<size; proc++) {
1171     if (!len_s[proc]) continue;
1172     i    = owners_co[proc];
1173     ierr = MPI_Isend(coj+coi[i],len_s[proc],MPIU_INT,proc,tagj,comm,swaits+k);CHKERRQ(ierr);
1174     k++;
1175   }
1176 
1177   /* receives and sends of coj are complete */
1178   ierr = PetscMalloc(size*sizeof(MPI_Status),&sstatus);CHKERRQ(ierr);
1179   for (i=0; i<merge->nrecv; i++) {
1180     PetscMPIInt icompleted;
1181     ierr = MPI_Waitany(merge->nrecv,rwaits,&icompleted,&rstatus);CHKERRQ(ierr);
1182   }
1183   ierr = PetscFree(rwaits);CHKERRQ(ierr);
1184   if (merge->nsend) {ierr = MPI_Waitall(merge->nsend,swaits,sstatus);CHKERRQ(ierr);}
1185 
1186   /* send and recv coi */
1187   /*-------------------*/
1188   ierr   = PetscCommGetNewTag(comm,&tagi);CHKERRQ(ierr);
1189   ierr   = PetscPostIrecvInt(comm,tagi,merge->nrecv,merge->id_r,len_ri,&buf_ri,&rwaits);CHKERRQ(ierr);
1190   ierr   = PetscMalloc((len+1)*sizeof(PetscInt),&buf_s);CHKERRQ(ierr);
1191   buf_si = buf_s;  /* points to the beginning of k-th msg to be sent */
1192   for (proc=0,k=0; proc<size; proc++) {
1193     if (!len_s[proc]) continue;
1194     /* form outgoing message for i-structure:
1195          buf_si[0]:                 nrows to be sent
1196                [1:nrows]:           row index (global)
1197                [nrows+1:2*nrows+1]: i-structure index
1198     */
1199     /*-------------------------------------------*/
1200     nrows       = len_si[proc]/2 - 1;
1201     buf_si_i    = buf_si + nrows+1;
1202     buf_si[0]   = nrows;
1203     buf_si_i[0] = 0;
1204     nrows       = 0;
1205     for (i=owners_co[proc]; i<owners_co[proc+1]; i++) {
1206       nzi               = coi[i+1] - coi[i];
1207       buf_si_i[nrows+1] = buf_si_i[nrows] + nzi; /* i-structure */
1208       buf_si[nrows+1]   = prmap[i] -owners[proc]; /* local row index */
1209       nrows++;
1210     }
1211     ierr = MPI_Isend(buf_si,len_si[proc],MPIU_INT,proc,tagi,comm,swaits+k);CHKERRQ(ierr);
1212     k++;
1213     buf_si += len_si[proc];
1214   }
1215   i = merge->nrecv;
1216   while (i--) {
1217     PetscMPIInt icompleted;
1218     ierr = MPI_Waitany(merge->nrecv,rwaits,&icompleted,&rstatus);CHKERRQ(ierr);
1219   }
1220   ierr = PetscFree(rwaits);CHKERRQ(ierr);
1221   if (merge->nsend) {ierr = MPI_Waitall(merge->nsend,swaits,sstatus);CHKERRQ(ierr);}
1222   ierr = PetscFree(len_si);CHKERRQ(ierr);
1223   ierr = PetscFree(len_ri);CHKERRQ(ierr);
1224   ierr = PetscFree(swaits);CHKERRQ(ierr);
1225   ierr = PetscFree(sstatus);CHKERRQ(ierr);
1226   ierr = PetscFree(buf_s);CHKERRQ(ierr);
1227 
1228   /* compute the local portion of C (mpi mat) */
1229   /*------------------------------------------*/
1230   /* allocate bi array and free space for accumulating nonzero column info */
1231   ierr  = PetscMalloc((pn+1)*sizeof(PetscInt),&bi);CHKERRQ(ierr);
1232   bi[0] = 0;
1233 
1234   /* set initial free space to be fill*(nnz(P) + nnz(A)) */
1235   nnz           = fill*(pdti[pn] + poti[pon] + ai[am]);
1236   ierr          = PetscFreeSpaceGet(nnz,&free_space);CHKERRQ(ierr);
1237   current_space = free_space;
1238 
1239   ierr = PetscMalloc3(merge->nrecv,PetscInt**,&buf_ri_k,merge->nrecv,PetscInt*,&nextrow,merge->nrecv,PetscInt*,&nextci);CHKERRQ(ierr);
1240   for (k=0; k<merge->nrecv; k++) {
1241     buf_ri_k[k] = buf_ri[k]; /* beginning of k-th recved i-structure */
1242     nrows       = *buf_ri_k[k];
1243     nextrow[k]  = buf_ri_k[k] + 1;  /* next row number of k-th recved i-structure */
1244     nextci[k]   = buf_ri_k[k] + (nrows + 1); /* poins to the next i-structure of k-th recved i-structure  */
1245   }
1246 
1247   ierr = MatPreallocateInitialize(comm,pn,A->cmap->n,dnz,onz);CHKERRQ(ierr);
1248   rmax = 0;
1249   for (i=0; i<pn; i++) {
1250     /* add pdt[i,:]*AP into lnk */
1251     pnz = pdti[i+1] - pdti[i];
1252     ptJ = pdtj + pdti[i];
1253     for (j=0; j<pnz; j++) {
1254       row  = ptJ[j];  /* row of AP == col of Pt */
1255       anz  = ai[row+1] - ai[row];
1256       Jptr = aj + ai[row];
1257       /* add non-zero cols of AP into the sorted linked list lnk */
1258       ierr = PetscLLCondensedAddSorted(anz,Jptr,lnk,lnkbt);CHKERRQ(ierr);
1259     }
1260 
1261     /* add received col data into lnk */
1262     for (k=0; k<merge->nrecv; k++) { /* k-th received message */
1263       if (i == *nextrow[k]) { /* i-th row */
1264         nzi  = *(nextci[k]+1) - *nextci[k];
1265         Jptr = buf_rj[k] + *nextci[k];
1266         ierr = PetscLLCondensedAddSorted(nzi,Jptr,lnk,lnkbt);CHKERRQ(ierr);
1267         nextrow[k]++; nextci[k]++;
1268       }
1269     }
1270     nnz = lnk[0];
1271 
1272     /* if free space is not available, make more free space */
1273     if (current_space->local_remaining<nnz) {
1274       ierr = PetscFreeSpaceGet(nnz+current_space->total_array_size,&current_space);CHKERRQ(ierr);
1275       nspacedouble++;
1276     }
1277     /* copy data into free space, then initialize lnk */
1278     ierr = PetscLLCondensedClean(aN,nnz,current_space->array,lnk,lnkbt);CHKERRQ(ierr);
1279     ierr = MatPreallocateSet(i+owners[rank],nnz,current_space->array,dnz,onz);CHKERRQ(ierr);
1280 
1281     current_space->array           += nnz;
1282     current_space->local_used      += nnz;
1283     current_space->local_remaining -= nnz;
1284 
1285     bi[i+1] = bi[i] + nnz;
1286     if (nnz > rmax) rmax = nnz;
1287   }
1288   ierr = PetscFree3(buf_ri_k,nextrow,nextci);CHKERRQ(ierr);
1289 
1290   ierr = PetscMalloc((bi[pn]+1)*sizeof(PetscInt),&bj);CHKERRQ(ierr);
1291   ierr = PetscFreeSpaceContiguous(&free_space,bj);CHKERRQ(ierr);
1292 
1293   afill_tmp = (PetscReal)bi[pn]/(pdti[pn] + poti[pon] + ai[am]+1);
1294   if (afill_tmp > afill) afill = afill_tmp;
1295   ierr = PetscLLCondensedDestroy(lnk,lnkbt);CHKERRQ(ierr);
1296   ierr = MatDestroy(&POt);CHKERRQ(ierr);
1297   ierr = MatDestroy(&PDt);CHKERRQ(ierr);
1298 
1299   /* create symbolic parallel matrix Cmpi - why cannot be assembled in Numeric part   */
1300   /*----------------------------------------------------------------------------------*/
1301   ierr = PetscMalloc((rmax+1)*sizeof(PetscScalar),&vals);CHKERRQ(ierr);
1302   ierr = PetscMemzero(vals,rmax*sizeof(PetscScalar));CHKERRQ(ierr);
1303 
1304   ierr = MatCreate(comm,&Cmpi);CHKERRQ(ierr);
1305   ierr = MatSetSizes(Cmpi,pn,A->cmap->n,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr);
1306   ierr = MatSetBlockSizes(Cmpi,P->cmap->bs,A->cmap->bs);CHKERRQ(ierr);
1307   ierr = MatSetType(Cmpi,MATMPIAIJ);CHKERRQ(ierr);
1308   ierr = MatMPIAIJSetPreallocation(Cmpi,0,dnz,0,onz);CHKERRQ(ierr);
1309   ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr);
1310   ierr = MatSetBlockSize(Cmpi,1);CHKERRQ(ierr);
1311   for (i=0; i<pn; i++) {
1312     row  = i + rstart;
1313     nnz  = bi[i+1] - bi[i];
1314     Jptr = bj + bi[i];
1315     ierr = MatSetValues(Cmpi,1,&row,nnz,Jptr,vals,INSERT_VALUES);CHKERRQ(ierr);
1316   }
1317   ierr = MatAssemblyBegin(Cmpi,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1318   ierr = MatAssemblyEnd(Cmpi,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1319   ierr = PetscFree(vals);CHKERRQ(ierr);
1320 
1321   merge->bi        = bi;
1322   merge->bj        = bj;
1323   merge->coi       = coi;
1324   merge->coj       = coj;
1325   merge->buf_ri    = buf_ri;
1326   merge->buf_rj    = buf_rj;
1327   merge->owners_co = owners_co;
1328   merge->destroy   = Cmpi->ops->destroy;
1329   merge->duplicate = Cmpi->ops->duplicate;
1330 
1331   Cmpi->ops->mattransposemultnumeric = MatTransposeMatMultNumeric_MPIAIJ_MPIAIJ;
1332   Cmpi->ops->destroy                 = MatDestroy_MPIAIJ_PtAP;
1333 
1334   /* attach the supporting struct to Cmpi for reuse */
1335   c           = (Mat_MPIAIJ*)Cmpi->data;
1336   c->ptap     = ptap;
1337   ptap->api   = PETSC_NULL;
1338   ptap->apj   = PETSC_NULL;
1339   ptap->merge = merge;
1340   ptap->rmax  = rmax;
1341 
1342   *C = Cmpi;
1343 #if defined(PETSC_USE_INFO)
1344   if (bi[pn] != 0) {
1345     ierr = PetscInfo3(Cmpi,"Reallocs %D; Fill ratio: given %G needed %G.\n",nspacedouble,fill,afill);CHKERRQ(ierr);
1346     ierr = PetscInfo1(Cmpi,"Use MatTransposeMatMult(A,B,MatReuse,%G,&C) for best performance.\n",afill);CHKERRQ(ierr);
1347   } else {
1348     ierr = PetscInfo(Cmpi,"Empty matrix product\n");CHKERRQ(ierr);
1349   }
1350 #endif
1351   PetscFunctionReturn(0);
1352 }
1353 
1354 #undef __FUNCT__
1355 #define __FUNCT__ "MatTransposeMatMultNumeric_MPIAIJ_MPIAIJ_Scalable"
1356 PetscErrorCode MatTransposeMatMultNumeric_MPIAIJ_MPIAIJ_Scalable(Mat P,Mat A,Mat C)
1357 {
1358   PetscErrorCode      ierr;
1359   Mat_Merge_SeqsToMPI *merge;
1360   Mat_MPIAIJ          *p =(Mat_MPIAIJ*)P->data,*c=(Mat_MPIAIJ*)C->data;
1361   Mat_SeqAIJ          *pd=(Mat_SeqAIJ*)(p->A)->data,*po=(Mat_SeqAIJ*)(p->B)->data;
1362   Mat_PtAPMPI         *ptap;
1363   PetscInt            *adj;
1364   PetscInt            i,j,k,anz,pnz,row,*cj,nexta;
1365   MatScalar           *ada,*ca,valtmp;
1366   PetscInt            am  =A->rmap->n,cm=C->rmap->n,pon=(p->B)->cmap->n;
1367   MPI_Comm            comm=((PetscObject)C)->comm;
1368   PetscMPIInt         size,rank,taga,*len_s;
1369   PetscInt            *owners,proc,nrows,**buf_ri_k,**nextrow,**nextci;
1370   PetscInt            **buf_ri,**buf_rj;
1371   PetscInt            cnz=0,*bj_i,*bi,*bj,bnz,nextcj;  /* bi,bj,ba: local array of C(mpi mat) */
1372   MPI_Request         *s_waits,*r_waits;
1373   MPI_Status          *status;
1374   MatScalar           **abuf_r,*ba_i,*pA,*coa,*ba;
1375   PetscInt            *ai,*aj,*coi,*coj;
1376   PetscInt            *poJ,*pdJ;
1377   Mat                 A_loc;
1378   Mat_SeqAIJ          *a_loc;
1379 
1380   PetscFunctionBegin;
1381   ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
1382   ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
1383 
1384   ptap  = c->ptap;
1385   merge = ptap->merge;
1386 
1387   /* 2) compute numeric C_seq = P_loc^T*A_loc */
1388   /*------------------------------------------*/
1389   /* get data from symbolic products */
1390   coi    = merge->coi; coj = merge->coj;
1391   ierr   = PetscMalloc((coi[pon]+1)*sizeof(MatScalar),&coa);CHKERRQ(ierr);
1392   ierr   = PetscMemzero(coa,coi[pon]*sizeof(MatScalar));CHKERRQ(ierr);
1393   bi     = merge->bi; bj = merge->bj;
1394   owners = merge->rowmap->range;
1395   ierr   = PetscMalloc((bi[cm]+1)*sizeof(MatScalar),&ba);CHKERRQ(ierr);
1396   ierr   = PetscMemzero(ba,bi[cm]*sizeof(MatScalar));CHKERRQ(ierr);
1397 
1398   /* get A_loc by taking all local rows of A */
1399   A_loc = ptap->A_loc;
1400   ierr  = MatMPIAIJGetLocalMat(A,MAT_REUSE_MATRIX,&A_loc);CHKERRQ(ierr);
1401   a_loc = (Mat_SeqAIJ*)(A_loc)->data;
1402   ai    = a_loc->i;
1403   aj    = a_loc->j;
1404 
1405   for (i=0; i<am; i++) {
1406     anz = ai[i+1] - ai[i];
1407     adj = aj + ai[i];
1408     ada = a_loc->a + ai[i];
1409 
1410     /* 2-b) Compute Cseq = P_loc[i,:]^T*A[i,:] using outer product */
1411     /*-------------------------------------------------------------*/
1412     /* put the value into Co=(p->B)^T*A (off-diagonal part, send to others) */
1413     pnz = po->i[i+1] - po->i[i];
1414     poJ = po->j + po->i[i];
1415     pA  = po->a + po->i[i];
1416     for (j=0; j<pnz; j++) {
1417       row = poJ[j];
1418       cj  = coj + coi[row];
1419       ca  = coa + coi[row];
1420       /* perform sparse axpy */
1421       nexta  = 0;
1422       valtmp = pA[j];
1423       for (k=0; nexta<anz; k++) {
1424         if (cj[k] == adj[nexta]) {
1425           ca[k] += valtmp*ada[nexta];
1426           nexta++;
1427         }
1428       }
1429       ierr = PetscLogFlops(2.0*anz);CHKERRQ(ierr);
1430     }
1431 
1432     /* put the value into Cd (diagonal part) */
1433     pnz = pd->i[i+1] - pd->i[i];
1434     pdJ = pd->j + pd->i[i];
1435     pA  = pd->a + pd->i[i];
1436     for (j=0; j<pnz; j++) {
1437       row = pdJ[j];
1438       cj  = bj + bi[row];
1439       ca  = ba + bi[row];
1440       /* perform sparse axpy */
1441       nexta  = 0;
1442       valtmp = pA[j];
1443       for (k=0; nexta<anz; k++) {
1444         if (cj[k] == adj[nexta]) {
1445           ca[k] += valtmp*ada[nexta];
1446           nexta++;
1447         }
1448       }
1449       ierr = PetscLogFlops(2.0*anz);CHKERRQ(ierr);
1450     }
1451   }
1452 
1453   /* 3) send and recv matrix values coa */
1454   /*------------------------------------*/
1455   buf_ri = merge->buf_ri;
1456   buf_rj = merge->buf_rj;
1457   len_s  = merge->len_s;
1458   ierr   = PetscCommGetNewTag(comm,&taga);CHKERRQ(ierr);
1459   ierr   = PetscPostIrecvScalar(comm,taga,merge->nrecv,merge->id_r,merge->len_r,&abuf_r,&r_waits);CHKERRQ(ierr);
1460 
1461   ierr = PetscMalloc2(merge->nsend+1,MPI_Request,&s_waits,size,MPI_Status,&status);CHKERRQ(ierr);
1462   for (proc=0,k=0; proc<size; proc++) {
1463     if (!len_s[proc]) continue;
1464     i    = merge->owners_co[proc];
1465     ierr = MPI_Isend(coa+coi[i],len_s[proc],MPIU_MATSCALAR,proc,taga,comm,s_waits+k);CHKERRQ(ierr);
1466     k++;
1467   }
1468   if (merge->nrecv) {ierr = MPI_Waitall(merge->nrecv,r_waits,status);CHKERRQ(ierr);}
1469   if (merge->nsend) {ierr = MPI_Waitall(merge->nsend,s_waits,status);CHKERRQ(ierr);}
1470 
1471   ierr = PetscFree2(s_waits,status);CHKERRQ(ierr);
1472   ierr = PetscFree(r_waits);CHKERRQ(ierr);
1473   ierr = PetscFree(coa);CHKERRQ(ierr);
1474 
1475   /* 4) insert local Cseq and received values into Cmpi */
1476   /*----------------------------------------------------*/
1477   ierr = PetscMalloc3(merge->nrecv,PetscInt**,&buf_ri_k,merge->nrecv,PetscInt*,&nextrow,merge->nrecv,PetscInt*,&nextci);CHKERRQ(ierr);
1478   for (k=0; k<merge->nrecv; k++) {
1479     buf_ri_k[k] = buf_ri[k]; /* beginning of k-th recved i-structure */
1480     nrows       = *(buf_ri_k[k]);
1481     nextrow[k]  = buf_ri_k[k]+1;  /* next row number of k-th recved i-structure */
1482     nextci[k]   = buf_ri_k[k] + (nrows + 1); /* poins to the next i-structure of k-th recved i-structure  */
1483   }
1484 
1485   for (i=0; i<cm; i++) {
1486     row  = owners[rank] + i; /* global row index of C_seq */
1487     bj_i = bj + bi[i];  /* col indices of the i-th row of C */
1488     ba_i = ba + bi[i];
1489     bnz  = bi[i+1] - bi[i];
1490     /* add received vals into ba */
1491     for (k=0; k<merge->nrecv; k++) { /* k-th received message */
1492       /* i-th row */
1493       if (i == *nextrow[k]) {
1494         cnz    = *(nextci[k]+1) - *nextci[k];
1495         cj     = buf_rj[k] + *(nextci[k]);
1496         ca     = abuf_r[k] + *(nextci[k]);
1497         nextcj = 0;
1498         for (j=0; nextcj<cnz; j++) {
1499           if (bj_i[j] == cj[nextcj]) { /* bcol == ccol */
1500             ba_i[j] += ca[nextcj++];
1501           }
1502         }
1503         nextrow[k]++; nextci[k]++;
1504         ierr = PetscLogFlops(2.0*cnz);CHKERRQ(ierr);
1505       }
1506     }
1507     ierr = MatSetValues(C,1,&row,bnz,bj_i,ba_i,INSERT_VALUES);CHKERRQ(ierr);
1508   }
1509   ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1510   ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1511 
1512   ierr = PetscFree(ba);CHKERRQ(ierr);
1513   ierr = PetscFree(abuf_r[0]);CHKERRQ(ierr);
1514   ierr = PetscFree(abuf_r);CHKERRQ(ierr);
1515   ierr = PetscFree3(buf_ri_k,nextrow,nextci);CHKERRQ(ierr);
1516   PetscFunctionReturn(0);
1517 }
1518 
1519 /* This routine is modified from MatPtAPSymbolic_MPIAIJ_MPIAIJ() */
1520 #undef __FUNCT__
1521 #define __FUNCT__ "MatTransposeMatMultSymbolic_MPIAIJ_MPIAIJ_Scalable"
1522 PetscErrorCode MatTransposeMatMultSymbolic_MPIAIJ_MPIAIJ_Scalable(Mat P,Mat A,PetscReal fill,Mat *C)
1523 {
1524   PetscErrorCode      ierr;
1525   Mat                 Cmpi,A_loc,POt,PDt;
1526   Mat_PtAPMPI         *ptap;
1527   PetscFreeSpaceList  free_space=PETSC_NULL,current_space=PETSC_NULL;
1528   Mat_MPIAIJ          *p        =(Mat_MPIAIJ*)P->data,*c;
1529   PetscInt            *pdti,*pdtj,*poti,*potj,*ptJ;
1530   PetscInt            nnz;
1531   PetscInt            *lnk,*owners_co,*coi,*coj,i,k,pnz,row;
1532   PetscInt            am  =A->rmap->n,pn=P->cmap->n;
1533   MPI_Comm            comm=((PetscObject)A)->comm;
1534   PetscMPIInt         size,rank,tagi,tagj,*len_si,*len_s,*len_ri;
1535   PetscInt            **buf_rj,**buf_ri,**buf_ri_k;
1536   PetscInt            len,proc,*dnz,*onz,*owners;
1537   PetscInt            nzi,*bi,*bj;
1538   PetscInt            nrows,*buf_s,*buf_si,*buf_si_i,**nextrow,**nextci;
1539   MPI_Request         *swaits,*rwaits;
1540   MPI_Status          *sstatus,rstatus;
1541   Mat_Merge_SeqsToMPI *merge;
1542   PetscInt            *ai,*aj,*Jptr,anz,*prmap=p->garray,pon,nspacedouble=0,j;
1543   PetscReal           afill  =1.0,afill_tmp;
1544   PetscInt            rstart = P->cmap->rstart,rmax,aN=A->cmap->N,Crmax;
1545   PetscScalar         *vals;
1546   Mat_SeqAIJ          *a_loc, *pdt,*pot;
1547 
1548   PetscFunctionBegin;
1549   /* check if matrix local sizes are compatible */
1550   if (A->rmap->rstart != P->rmap->rstart || A->rmap->rend != P->rmap->rend) {
1551     SETERRQ4(comm,PETSC_ERR_ARG_SIZ,"Matrix local dimensions are incompatible, A (%D, %D) != P (%D,%D)",A->rmap->rstart,A->rmap->rend,P->rmap->rstart,P->rmap->rend);
1552   }
1553 
1554   ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
1555   ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
1556 
1557   /* create struct Mat_PtAPMPI and attached it to C later */
1558   ierr = PetscNew(Mat_PtAPMPI,&ptap);CHKERRQ(ierr);
1559 
1560   /* get A_loc by taking all local rows of A */
1561   ierr = MatMPIAIJGetLocalMat(A,MAT_INITIAL_MATRIX,&A_loc);CHKERRQ(ierr);
1562 
1563   ptap->A_loc = A_loc;
1564   a_loc       = (Mat_SeqAIJ*)(A_loc)->data;
1565   ai          = a_loc->i;
1566   aj          = a_loc->j;
1567 
1568   /* determine symbolic Co=(p->B)^T*A - send to others */
1569   /*----------------------------------------------------*/
1570   ierr = MatTransposeSymbolic_SeqAIJ(p->A,&PDt);CHKERRQ(ierr);
1571   pdt  = (Mat_SeqAIJ*)PDt->data;
1572   pdti = pdt->i; pdtj = pdt->j;
1573 
1574   ierr = MatTransposeSymbolic_SeqAIJ(p->B,&POt);CHKERRQ(ierr);
1575   pot  = (Mat_SeqAIJ*)POt->data;
1576   poti = pot->i; potj = pot->j;
1577 
1578   /* then, compute symbolic Co = (p->B)^T*A */
1579   pon    = (p->B)->cmap->n; /* total num of rows to be sent to other processors
1580                          >= (num of nonzero rows of C_seq) - pn */
1581   ierr   = PetscMalloc((pon+1)*sizeof(PetscInt),&coi);CHKERRQ(ierr);
1582   coi[0] = 0;
1583 
1584   /* set initial free space to be fill*(nnz(p->B) + nnz(A)) */
1585   nnz           = fill*(poti[pon] + ai[am]);
1586   ierr          = PetscFreeSpaceGet(nnz,&free_space);CHKERRQ(ierr);
1587   current_space = free_space;
1588 
1589   /* create and initialize a linked list */
1590   i     = PetscMax(pdt->rmax,pot->rmax);
1591   Crmax = i*a_loc->rmax*size; /* non-scalable! */
1592   if (!Crmax || Crmax > aN) Crmax = aN;
1593   ierr = PetscLLCondensedCreate_Scalable(Crmax,&lnk);CHKERRQ(ierr);
1594 
1595   for (i=0; i<pon; i++) {
1596     pnz = poti[i+1] - poti[i];
1597     ptJ = potj + poti[i];
1598     for (j=0; j<pnz; j++) {
1599       row  = ptJ[j]; /* row of A_loc == col of Pot */
1600       anz  = ai[row+1] - ai[row];
1601       Jptr = aj + ai[row];
1602       /* add non-zero cols of AP into the sorted linked list lnk */
1603       ierr = PetscLLCondensedAddSorted_Scalable(anz,Jptr,lnk);CHKERRQ(ierr);
1604     }
1605     nnz = lnk[0];
1606 
1607     /* If free space is not available, double the total space in the list */
1608     if (current_space->local_remaining<nnz) {
1609       ierr = PetscFreeSpaceGet(nnz+current_space->total_array_size,&current_space);CHKERRQ(ierr);
1610       nspacedouble++;
1611     }
1612 
1613     /* Copy data into free space, and zero out denserows */
1614     ierr = PetscLLCondensedClean_Scalable(nnz,current_space->array,lnk);CHKERRQ(ierr);
1615 
1616     current_space->array           += nnz;
1617     current_space->local_used      += nnz;
1618     current_space->local_remaining -= nnz;
1619 
1620     coi[i+1] = coi[i] + nnz;
1621   }
1622 
1623   ierr = PetscMalloc((coi[pon]+1)*sizeof(PetscInt),&coj);CHKERRQ(ierr);
1624   ierr = PetscFreeSpaceContiguous(&free_space,coj);CHKERRQ(ierr);
1625 
1626   afill_tmp = (PetscReal)coi[pon]/(poti[pon] + ai[am]+1);
1627   if (afill_tmp > afill) afill = afill_tmp;
1628 
1629   /* send j-array (coj) of Co to other processors */
1630   /*----------------------------------------------*/
1631   /* determine row ownership */
1632   ierr = PetscNew(Mat_Merge_SeqsToMPI,&merge);CHKERRQ(ierr);
1633   ierr = PetscLayoutCreate(comm,&merge->rowmap);CHKERRQ(ierr);
1634 
1635   merge->rowmap->n  = pn;
1636   merge->rowmap->bs = 1;
1637 
1638   ierr   = PetscLayoutSetUp(merge->rowmap);CHKERRQ(ierr);
1639   owners = merge->rowmap->range;
1640 
1641   /* determine the number of messages to send, their lengths */
1642   ierr = PetscMalloc(size*sizeof(PetscMPIInt),&len_si);CHKERRQ(ierr);
1643   ierr = PetscMemzero(len_si,size*sizeof(PetscMPIInt));CHKERRQ(ierr);
1644   ierr = PetscMalloc(size*sizeof(PetscMPIInt),&merge->len_s);CHKERRQ(ierr);
1645 
1646   len_s        = merge->len_s;
1647   merge->nsend = 0;
1648 
1649   ierr = PetscMalloc((size+2)*sizeof(PetscInt),&owners_co);CHKERRQ(ierr);
1650   ierr = PetscMemzero(len_s,size*sizeof(PetscMPIInt));CHKERRQ(ierr);
1651 
1652   proc = 0;
1653   for (i=0; i<pon; i++) {
1654     while (prmap[i] >= owners[proc+1]) proc++;
1655     len_si[proc]++;  /* num of rows in Co to be sent to [proc] */
1656     len_s[proc] += coi[i+1] - coi[i];
1657   }
1658 
1659   len          = 0; /* max length of buf_si[] */
1660   owners_co[0] = 0;
1661   for (proc=0; proc<size; proc++) {
1662     owners_co[proc+1] = owners_co[proc] + len_si[proc];
1663     if (len_si[proc]) {
1664       merge->nsend++;
1665       len_si[proc] = 2*(len_si[proc] + 1);
1666       len         += len_si[proc];
1667     }
1668   }
1669 
1670   /* determine the number and length of messages to receive for coi and coj  */
1671   ierr = PetscGatherNumberOfMessages(comm,PETSC_NULL,len_s,&merge->nrecv);CHKERRQ(ierr);
1672   ierr = PetscGatherMessageLengths2(comm,merge->nsend,merge->nrecv,len_s,len_si,&merge->id_r,&merge->len_r,&len_ri);CHKERRQ(ierr);
1673 
1674   /* post the Irecv and Isend of coj */
1675   ierr = PetscCommGetNewTag(comm,&tagj);CHKERRQ(ierr);
1676   ierr = PetscPostIrecvInt(comm,tagj,merge->nrecv,merge->id_r,merge->len_r,&buf_rj,&rwaits);CHKERRQ(ierr);
1677   ierr = PetscMalloc((merge->nsend+1)*sizeof(MPI_Request),&swaits);CHKERRQ(ierr);
1678   for (proc=0, k=0; proc<size; proc++) {
1679     if (!len_s[proc]) continue;
1680     i    = owners_co[proc];
1681     ierr = MPI_Isend(coj+coi[i],len_s[proc],MPIU_INT,proc,tagj,comm,swaits+k);CHKERRQ(ierr);
1682     k++;
1683   }
1684 
1685   /* receives and sends of coj are complete */
1686   ierr = PetscMalloc(size*sizeof(MPI_Status),&sstatus);CHKERRQ(ierr);
1687   for (i=0; i<merge->nrecv; i++) {
1688     PetscMPIInt icompleted;
1689     ierr = MPI_Waitany(merge->nrecv,rwaits,&icompleted,&rstatus);CHKERRQ(ierr);
1690   }
1691   ierr = PetscFree(rwaits);CHKERRQ(ierr);
1692   if (merge->nsend) {ierr = MPI_Waitall(merge->nsend,swaits,sstatus);CHKERRQ(ierr);}
1693 
1694   /* send and recv coi */
1695   /*-------------------*/
1696   ierr   = PetscCommGetNewTag(comm,&tagi);CHKERRQ(ierr);
1697   ierr   = PetscPostIrecvInt(comm,tagi,merge->nrecv,merge->id_r,len_ri,&buf_ri,&rwaits);CHKERRQ(ierr);
1698   ierr   = PetscMalloc((len+1)*sizeof(PetscInt),&buf_s);CHKERRQ(ierr);
1699   buf_si = buf_s;  /* points to the beginning of k-th msg to be sent */
1700   for (proc=0,k=0; proc<size; proc++) {
1701     if (!len_s[proc]) continue;
1702     /* form outgoing message for i-structure:
1703          buf_si[0]:                 nrows to be sent
1704                [1:nrows]:           row index (global)
1705                [nrows+1:2*nrows+1]: i-structure index
1706     */
1707     /*-------------------------------------------*/
1708     nrows       = len_si[proc]/2 - 1;
1709     buf_si_i    = buf_si + nrows+1;
1710     buf_si[0]   = nrows;
1711     buf_si_i[0] = 0;
1712     nrows       = 0;
1713     for (i=owners_co[proc]; i<owners_co[proc+1]; i++) {
1714       nzi               = coi[i+1] - coi[i];
1715       buf_si_i[nrows+1] = buf_si_i[nrows] + nzi;  /* i-structure */
1716       buf_si[nrows+1]   = prmap[i] -owners[proc]; /* local row index */
1717       nrows++;
1718     }
1719     ierr = MPI_Isend(buf_si,len_si[proc],MPIU_INT,proc,tagi,comm,swaits+k);CHKERRQ(ierr);
1720     k++;
1721     buf_si += len_si[proc];
1722   }
1723   i = merge->nrecv;
1724   while (i--) {
1725     PetscMPIInt icompleted;
1726     ierr = MPI_Waitany(merge->nrecv,rwaits,&icompleted,&rstatus);CHKERRQ(ierr);
1727   }
1728   ierr = PetscFree(rwaits);CHKERRQ(ierr);
1729   if (merge->nsend) {ierr = MPI_Waitall(merge->nsend,swaits,sstatus);CHKERRQ(ierr);}
1730   ierr = PetscFree(len_si);CHKERRQ(ierr);
1731   ierr = PetscFree(len_ri);CHKERRQ(ierr);
1732   ierr = PetscFree(swaits);CHKERRQ(ierr);
1733   ierr = PetscFree(sstatus);CHKERRQ(ierr);
1734   ierr = PetscFree(buf_s);CHKERRQ(ierr);
1735 
1736   /* compute the local portion of C (mpi mat) */
1737   /*------------------------------------------*/
1738   /* allocate bi array and free space for accumulating nonzero column info */
1739   ierr  = PetscMalloc((pn+1)*sizeof(PetscInt),&bi);CHKERRQ(ierr);
1740   bi[0] = 0;
1741 
1742   /* set initial free space to be fill*(nnz(P) + nnz(AP)) */
1743   nnz           = fill*(pdti[pn] + poti[pon] + ai[am]);
1744   ierr          = PetscFreeSpaceGet(nnz,&free_space);CHKERRQ(ierr);
1745   current_space = free_space;
1746 
1747   ierr = PetscMalloc3(merge->nrecv,PetscInt**,&buf_ri_k,merge->nrecv,PetscInt*,&nextrow,merge->nrecv,PetscInt*,&nextci);CHKERRQ(ierr);
1748   for (k=0; k<merge->nrecv; k++) {
1749     buf_ri_k[k] = buf_ri[k]; /* beginning of k-th recved i-structure */
1750     nrows       = *buf_ri_k[k];
1751     nextrow[k]  = buf_ri_k[k] + 1;  /* next row number of k-th recved i-structure */
1752     nextci[k]   = buf_ri_k[k] + (nrows + 1); /* points to the next i-structure of k-th recieved i-structure  */
1753   }
1754 
1755   ierr = MatPreallocateInitialize(comm,pn,A->cmap->n,dnz,onz);CHKERRQ(ierr);
1756   rmax = 0;
1757   for (i=0; i<pn; i++) {
1758     /* add pdt[i,:]*AP into lnk */
1759     pnz = pdti[i+1] - pdti[i];
1760     ptJ = pdtj + pdti[i];
1761     for (j=0; j<pnz; j++) {
1762       row  = ptJ[j];  /* row of AP == col of Pt */
1763       anz  = ai[row+1] - ai[row];
1764       Jptr = aj + ai[row];
1765       /* add non-zero cols of AP into the sorted linked list lnk */
1766       ierr = PetscLLCondensedAddSorted_Scalable(anz,Jptr,lnk);CHKERRQ(ierr);
1767     }
1768 
1769     /* add received col data into lnk */
1770     for (k=0; k<merge->nrecv; k++) { /* k-th received message */
1771       if (i == *nextrow[k]) { /* i-th row */
1772         nzi  = *(nextci[k]+1) - *nextci[k];
1773         Jptr = buf_rj[k] + *nextci[k];
1774         ierr = PetscLLCondensedAddSorted_Scalable(nzi,Jptr,lnk);CHKERRQ(ierr);
1775         nextrow[k]++; nextci[k]++;
1776       }
1777     }
1778     nnz = lnk[0];
1779 
1780     /* if free space is not available, make more free space */
1781     if (current_space->local_remaining<nnz) {
1782       ierr = PetscFreeSpaceGet(nnz+current_space->total_array_size,&current_space);CHKERRQ(ierr);
1783       nspacedouble++;
1784     }
1785     /* copy data into free space, then initialize lnk */
1786     ierr = PetscLLCondensedClean_Scalable(nnz,current_space->array,lnk);CHKERRQ(ierr);
1787     ierr = MatPreallocateSet(i+owners[rank],nnz,current_space->array,dnz,onz);CHKERRQ(ierr);
1788 
1789     current_space->array           += nnz;
1790     current_space->local_used      += nnz;
1791     current_space->local_remaining -= nnz;
1792 
1793     bi[i+1] = bi[i] + nnz;
1794     if (nnz > rmax) rmax = nnz;
1795   }
1796   ierr = PetscFree3(buf_ri_k,nextrow,nextci);CHKERRQ(ierr);
1797 
1798   ierr      = PetscMalloc((bi[pn]+1)*sizeof(PetscInt),&bj);CHKERRQ(ierr);
1799   ierr      = PetscFreeSpaceContiguous(&free_space,bj);CHKERRQ(ierr);
1800   afill_tmp = (PetscReal)bi[pn]/(pdti[pn] + poti[pon] + ai[am]+1);
1801   if (afill_tmp > afill) afill = afill_tmp;
1802   ierr = PetscLLCondensedDestroy_Scalable(lnk);CHKERRQ(ierr);
1803   ierr = MatDestroy(&POt);CHKERRQ(ierr);
1804   ierr = MatDestroy(&PDt);CHKERRQ(ierr);
1805 
1806   /* create symbolic parallel matrix Cmpi - why cannot be assembled in Numeric part   */
1807   /*----------------------------------------------------------------------------------*/
1808   ierr = PetscMalloc((rmax+1)*sizeof(PetscScalar),&vals);CHKERRQ(ierr);
1809   ierr = PetscMemzero(vals,rmax*sizeof(PetscScalar));CHKERRQ(ierr);
1810 
1811   ierr = MatCreate(comm,&Cmpi);CHKERRQ(ierr);
1812   ierr = MatSetSizes(Cmpi,pn,A->cmap->n,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr);
1813   ierr = MatSetBlockSizes(Cmpi,P->cmap->bs,A->cmap->bs);CHKERRQ(ierr);
1814   ierr = MatSetType(Cmpi,MATMPIAIJ);CHKERRQ(ierr);
1815   ierr = MatMPIAIJSetPreallocation(Cmpi,0,dnz,0,onz);CHKERRQ(ierr);
1816   ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr);
1817   ierr = MatSetBlockSize(Cmpi,1);CHKERRQ(ierr);
1818   for (i=0; i<pn; i++) {
1819     row  = i + rstart;
1820     nnz  = bi[i+1] - bi[i];
1821     Jptr = bj + bi[i];
1822     ierr = MatSetValues(Cmpi,1,&row,nnz,Jptr,vals,INSERT_VALUES);CHKERRQ(ierr);
1823   }
1824   ierr = MatAssemblyBegin(Cmpi,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1825   ierr = MatAssemblyEnd(Cmpi,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1826   ierr = PetscFree(vals);CHKERRQ(ierr);
1827 
1828   merge->bi        = bi;
1829   merge->bj        = bj;
1830   merge->coi       = coi;
1831   merge->coj       = coj;
1832   merge->buf_ri    = buf_ri;
1833   merge->buf_rj    = buf_rj;
1834   merge->owners_co = owners_co;
1835   merge->destroy   = Cmpi->ops->destroy;
1836   merge->duplicate = Cmpi->ops->duplicate;
1837 
1838   Cmpi->ops->mattransposemultnumeric = MatTransposeMatMultNumeric_MPIAIJ_MPIAIJ_Scalable;
1839   Cmpi->ops->destroy                 = MatDestroy_MPIAIJ_PtAP;
1840 
1841   /* attach the supporting struct to Cmpi for reuse */
1842   c = (Mat_MPIAIJ*)Cmpi->data;
1843 
1844   c->ptap     = ptap;
1845   ptap->api   = PETSC_NULL;
1846   ptap->apj   = PETSC_NULL;
1847   ptap->merge = merge;
1848   ptap->rmax  = rmax;
1849   ptap->apa   = PETSC_NULL;
1850 
1851   *C = Cmpi;
1852 #if defined(PETSC_USE_INFO)
1853   if (bi[pn] != 0) {
1854     ierr = PetscInfo3(Cmpi,"Reallocs %D; Fill ratio: given %G needed %G.\n",nspacedouble,fill,afill);CHKERRQ(ierr);
1855     ierr = PetscInfo1(Cmpi,"Use MatTransposeMatMult(A,B,MatReuse,%G,&C) for best performance.\n",afill);CHKERRQ(ierr);
1856   } else {
1857     ierr = PetscInfo(Cmpi,"Empty matrix product\n");CHKERRQ(ierr);
1858   }
1859 #endif
1860   PetscFunctionReturn(0);
1861 }
1862