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