xref: /petsc/src/mat/impls/aij/mpi/mpimatmatmult.c (revision d54338ec4fdf48b56b5efdd9bf04ca9d50d6ab04) 
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 {
814       SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not supported");
815     }
816     PetscFunctionReturn(0);
817   }
818 
819   if (scall == MAT_INITIAL_MATRIX) {
820     ierr = PetscObjectOptionsBegin((PetscObject)A);CHKERRQ(ierr);
821       ierr = PetscOptionsBool("-mattransposematmult_scalable","Use a scalable but slower C=Pt*A","",scalable,&scalable,PETSC_NULL);CHKERRQ(ierr);
822       if  (scalable) {
823         ierr = MatTransposeMatMultSymbolic_MPIAIJ_MPIAIJ_Scalable(P,A,fill,C);CHKERRQ(ierr);
824       } else {
825         ierr = MatTransposeMatMultSymbolic_MPIAIJ_MPIAIJ(P,A,fill,C);CHKERRQ(ierr);
826       }
827     ierr = PetscOptionsEnd();CHKERRQ(ierr);
828   }
829   ierr = (*(*C)->ops->mattransposemultnumeric)(P,A,*C);CHKERRQ(ierr);
830   PetscFunctionReturn(0);
831 }
832 
833 #undef __FUNCT__
834 #define __FUNCT__ "MatTransposeMatMultNumeric_MPIAIJ_MPIAIJ"
835 PetscErrorCode MatTransposeMatMultNumeric_MPIAIJ_MPIAIJ(Mat P,Mat A,Mat C)
836 {
837   PetscErrorCode       ierr;
838   Mat_Merge_SeqsToMPI  *merge;
839   Mat_MPIAIJ           *p=(Mat_MPIAIJ*)P->data,*c=(Mat_MPIAIJ*)C->data;
840   Mat_SeqAIJ           *pd=(Mat_SeqAIJ*)(p->A)->data,*po=(Mat_SeqAIJ*)(p->B)->data;
841   Mat_PtAPMPI          *ptap;
842   PetscInt             *adj,*aJ;
843   PetscInt             i,j,k,anz,pnz,row,*cj;
844   MatScalar            *ada,*aval,*ca,valtmp;
845   PetscInt             am=A->rmap->n,cm=C->rmap->n,pon=(p->B)->cmap->n;
846   MPI_Comm             comm=((PetscObject)C)->comm;
847   PetscMPIInt          size,rank,taga,*len_s;
848   PetscInt             *owners,proc,nrows,**buf_ri_k,**nextrow,**nextci;
849   PetscInt             **buf_ri,**buf_rj;
850   PetscInt             cnz=0,*bj_i,*bi,*bj,bnz,nextcj; /* bi,bj,ba: local array of C(mpi mat) */
851   MPI_Request          *s_waits,*r_waits;
852   MPI_Status           *status;
853   MatScalar            **abuf_r,*ba_i,*pA,*coa,*ba;
854   PetscInt             *ai,*aj,*coi,*coj;
855   PetscInt             *poJ,*pdJ;
856   Mat                  A_loc;
857   Mat_SeqAIJ           *a_loc;
858 
859   PetscFunctionBegin;
860   ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
861   ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
862 
863   ptap  = c->ptap;
864   merge = ptap->merge;
865 
866   /* 2) compute numeric C_seq = P_loc^T*A_loc*P - dominating part */
867   /*--------------------------------------------------------------*/
868   /* get data from symbolic products */
869   coi = merge->coi; coj = merge->coj;
870   ierr = PetscMalloc((coi[pon]+1)*sizeof(MatScalar),&coa);CHKERRQ(ierr);
871   ierr = PetscMemzero(coa,coi[pon]*sizeof(MatScalar));CHKERRQ(ierr);
872 
873   bi     = merge->bi; bj = merge->bj;
874   owners = merge->rowmap->range;
875   ierr   = PetscMalloc((bi[cm]+1)*sizeof(MatScalar),&ba);CHKERRQ(ierr);
876   ierr   = PetscMemzero(ba,bi[cm]*sizeof(MatScalar));CHKERRQ(ierr);
877 
878   /* get A_loc by taking all local rows of A */
879   A_loc = ptap->A_loc;
880   ierr = MatMPIAIJGetLocalMat(A,MAT_REUSE_MATRIX,&A_loc);CHKERRQ(ierr);
881   a_loc = (Mat_SeqAIJ*)(A_loc)->data;
882   ai   = a_loc->i;
883   aj   = a_loc->j;
884 
885   ierr = PetscMalloc((A->cmap->N)*sizeof(PetscScalar),&aval);CHKERRQ(ierr); /* non-scalable!!! */
886   ierr = PetscMemzero(aval,A->cmap->N*sizeof(PetscScalar));CHKERRQ(ierr);
887 
888     for (i=0; i<am; i++) {
889       /* 2-a) put A[i,:] to dense array aval */
890       anz = ai[i+1] - ai[i];
891       adj = aj + ai[i];
892       ada = a_loc->a + ai[i];
893       for (j=0; j<anz; j++) {
894         aval[adj[j]] = ada[j];
895       }
896 
897       /* 2-b) Compute Cseq = P_loc[i,:]^T*A[i,:] using outer product */
898       /*--------------------------------------------------------------*/
899       /* put the value into Co=(p->B)^T*A (off-diagonal part, send to others) */
900       pnz = po->i[i+1] - po->i[i];
901       poJ = po->j + po->i[i];
902       pA  = po->a + po->i[i];
903       for (j=0; j<pnz; j++) {
904         row = poJ[j];
905         cnz = coi[row+1] - coi[row];
906         cj  = coj + coi[row];
907         ca  = coa + coi[row];
908         /* perform dense axpy */
909         valtmp = pA[j];
910         for (k=0; k<cnz; k++) {
911           ca[k] += valtmp*aval[cj[k]];
912         }
913         ierr = PetscLogFlops(2.0*cnz);CHKERRQ(ierr);
914       }
915 
916       /* put the value into Cd (diagonal part) */
917       pnz = pd->i[i+1] - pd->i[i];
918       pdJ = pd->j + pd->i[i];
919       pA  = pd->a + pd->i[i];
920       for (j=0; j<pnz; j++) {
921         row = pdJ[j];
922         cnz = bi[row+1] - bi[row];
923         cj  = bj + bi[row];
924         ca  = ba + bi[row];
925         /* perform dense axpy */
926         valtmp = pA[j];
927         for (k=0; k<cnz; k++) {
928           ca[k] += valtmp*aval[cj[k]];
929         }
930         ierr = PetscLogFlops(2.0*cnz);CHKERRQ(ierr);
931       }
932 
933       /* zero the current row of Pt*A */
934       aJ = aj + ai[i];
935       for (k=0; k<anz; k++) aval[aJ[k]] = 0.0;
936     }
937 
938   /* 3) send and recv matrix values coa */
939   /*------------------------------------*/
940   buf_ri = merge->buf_ri;
941   buf_rj = merge->buf_rj;
942   len_s  = merge->len_s;
943   ierr = PetscCommGetNewTag(comm,&taga);CHKERRQ(ierr);
944   ierr = PetscPostIrecvScalar(comm,taga,merge->nrecv,merge->id_r,merge->len_r,&abuf_r,&r_waits);CHKERRQ(ierr);
945 
946   ierr = PetscMalloc2(merge->nsend+1,MPI_Request,&s_waits,size,MPI_Status,&status);CHKERRQ(ierr);
947   for (proc=0,k=0; proc<size; proc++) {
948     if (!len_s[proc]) continue;
949     i = merge->owners_co[proc];
950     ierr = MPI_Isend(coa+coi[i],len_s[proc],MPIU_MATSCALAR,proc,taga,comm,s_waits+k);CHKERRQ(ierr);
951     k++;
952   }
953   if (merge->nrecv) {ierr = MPI_Waitall(merge->nrecv,r_waits,status);CHKERRQ(ierr);}
954   if (merge->nsend) {ierr = MPI_Waitall(merge->nsend,s_waits,status);CHKERRQ(ierr);}
955 
956   ierr = PetscFree2(s_waits,status);CHKERRQ(ierr);
957   ierr = PetscFree(r_waits);CHKERRQ(ierr);
958   ierr = PetscFree(coa);CHKERRQ(ierr);
959 
960   /* 4) insert local Cseq and received values into Cmpi */
961   /*----------------------------------------------------*/
962   ierr = PetscMalloc3(merge->nrecv,PetscInt**,&buf_ri_k,merge->nrecv,PetscInt*,&nextrow,merge->nrecv,PetscInt*,&nextci);CHKERRQ(ierr);
963   for (k=0; k<merge->nrecv; k++) {
964     buf_ri_k[k] = buf_ri[k]; /* beginning of k-th recved i-structure */
965     nrows       = *(buf_ri_k[k]);
966     nextrow[k]  = buf_ri_k[k]+1;  /* next row number of k-th recved i-structure */
967     nextci[k]   = buf_ri_k[k] + (nrows + 1);/* poins to the next i-structure of k-th recved i-structure  */
968   }
969 
970   for (i=0; i<cm; i++) {
971     row = owners[rank] + i; /* global row index of C_seq */
972     bj_i = bj + bi[i];  /* col indices of the i-th row of C */
973     ba_i = ba + bi[i];
974     bnz  = bi[i+1] - bi[i];
975     /* add received vals into ba */
976     for (k=0; k<merge->nrecv; k++) { /* k-th received message */
977       /* i-th row */
978       if (i == *nextrow[k]) {
979         cnz = *(nextci[k]+1) - *nextci[k];
980         cj  = buf_rj[k] + *(nextci[k]);
981         ca  = abuf_r[k] + *(nextci[k]);
982         nextcj = 0;
983         for (j=0; nextcj<cnz; j++) {
984           if (bj_i[j] == cj[nextcj]) { /* bcol == ccol */
985             ba_i[j] += ca[nextcj++];
986           }
987         }
988         nextrow[k]++; nextci[k]++;
989         ierr = PetscLogFlops(2.0*cnz);CHKERRQ(ierr);
990       }
991     }
992     ierr = MatSetValues(C,1,&row,bnz,bj_i,ba_i,INSERT_VALUES);CHKERRQ(ierr);
993   }
994   ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
995   ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
996 
997   ierr = PetscFree(ba);CHKERRQ(ierr);
998   ierr = PetscFree(abuf_r[0]);CHKERRQ(ierr);
999   ierr = PetscFree(abuf_r);CHKERRQ(ierr);
1000   ierr = PetscFree3(buf_ri_k,nextrow,nextci);CHKERRQ(ierr);
1001   ierr = PetscFree(aval);CHKERRQ(ierr);
1002   PetscFunctionReturn(0);
1003 }
1004 
1005 /* This routine is modified from MatPtAPSymbolic_MPIAIJ_MPIAIJ() */
1006 #undef __FUNCT__
1007 #define __FUNCT__ "MatTransposeMatMultSymbolic_MPIAIJ_MPIAIJ"
1008 PetscErrorCode MatTransposeMatMultSymbolic_MPIAIJ_MPIAIJ(Mat P,Mat A,PetscReal fill,Mat *C)
1009 {
1010   PetscErrorCode       ierr;
1011   Mat                  Cmpi,A_loc,POt,PDt;
1012   Mat_PtAPMPI          *ptap;
1013   PetscFreeSpaceList   free_space=PETSC_NULL,current_space=PETSC_NULL;
1014   Mat_MPIAIJ           *p=(Mat_MPIAIJ*)P->data,*c;
1015   PetscInt             *pdti,*pdtj,*poti,*potj,*ptJ;
1016   PetscInt             nnz;
1017   PetscInt             *lnk,*owners_co,*coi,*coj,i,k,pnz,row;
1018   PetscInt             am=A->rmap->n,pn=P->cmap->n;
1019   PetscBT              lnkbt;
1020   MPI_Comm             comm=((PetscObject)A)->comm;
1021   PetscMPIInt          size,rank,tagi,tagj,*len_si,*len_s,*len_ri;
1022   PetscInt             **buf_rj,**buf_ri,**buf_ri_k;
1023   PetscInt             len,proc,*dnz,*onz,*owners;
1024   PetscInt             nzi,*bi,*bj;
1025   PetscInt             nrows,*buf_s,*buf_si,*buf_si_i,**nextrow,**nextci;
1026   MPI_Request          *swaits,*rwaits;
1027   MPI_Status           *sstatus,rstatus;
1028   Mat_Merge_SeqsToMPI  *merge;
1029   PetscInt             *ai,*aj,*Jptr,anz,*prmap=p->garray,pon,nspacedouble=0,j;
1030   PetscReal            afill=1.0,afill_tmp;
1031   PetscInt             rstart = P->cmap->rstart,rmax,aN=A->cmap->N,Crmax;
1032   PetscScalar          *vals;
1033   Mat_SeqAIJ           *a_loc, *pdt,*pot;
1034 
1035   PetscFunctionBegin;
1036   /* check if matrix local sizes are compatible */
1037   if (A->rmap->rstart != P->rmap->rstart || A->rmap->rend != P->rmap->rend) {
1038     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);
1039   }
1040 
1041   ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
1042   ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
1043 
1044   /* create struct Mat_PtAPMPI and attached it to C later */
1045   ierr = PetscNew(Mat_PtAPMPI,&ptap);CHKERRQ(ierr);
1046 
1047   /* get A_loc by taking all local rows of A */
1048   ierr = MatMPIAIJGetLocalMat(A,MAT_INITIAL_MATRIX,&A_loc);CHKERRQ(ierr);
1049   ptap->A_loc = A_loc;
1050   a_loc = (Mat_SeqAIJ*)(A_loc)->data;
1051   ai   = a_loc->i;
1052   aj   = a_loc->j;
1053 
1054   /* determine symbolic Co=(p->B)^T*A - send to others */
1055   /*----------------------------------------------------*/
1056   ierr = MatTransposeSymbolic_SeqAIJ(p->A,&PDt);CHKERRQ(ierr);
1057   pdt = (Mat_SeqAIJ*)PDt->data;
1058   pdti = pdt->i; pdtj = pdt->j;
1059 
1060   ierr = MatTransposeSymbolic_SeqAIJ(p->B,&POt);CHKERRQ(ierr);
1061   pot = (Mat_SeqAIJ*)POt->data;
1062   poti = pot->i; potj = pot->j;
1063 
1064   /* then, compute symbolic Co = (p->B)^T*A */
1065   pon = (p->B)->cmap->n; /* total num of rows to be sent to other processors
1066                          >= (num of nonzero rows of C_seq) - pn */
1067   ierr = PetscMalloc((pon+1)*sizeof(PetscInt),&coi);CHKERRQ(ierr);
1068   coi[0] = 0;
1069 
1070   /* set initial free space to be fill*(nnz(p->B) + nnz(A)) */
1071   nnz           = fill*(poti[pon] + ai[am]);
1072   ierr          = PetscFreeSpaceGet(nnz,&free_space);CHKERRQ(ierr);
1073   current_space = free_space;
1074 
1075   /* create and initialize a linked list */
1076   i = PetscMax(pdt->rmax,pot->rmax);
1077   Crmax = i*a_loc->rmax*size;
1078   if (!Crmax || Crmax > aN) Crmax = aN;
1079   ierr = PetscLLCondensedCreate(Crmax,aN,&lnk,&lnkbt);CHKERRQ(ierr);
1080 
1081   for (i=0; i<pon; i++) {
1082     pnz = poti[i+1] - poti[i];
1083     ptJ = potj + poti[i];
1084     for (j=0; j<pnz; j++) {
1085       row  = ptJ[j]; /* row of A_loc == col of Pot */
1086       anz  = ai[row+1] - ai[row];
1087       Jptr = aj + ai[row];
1088       /* add non-zero cols of AP into the sorted linked list lnk */
1089       ierr = PetscLLCondensedAddSorted(anz,Jptr,lnk,lnkbt);CHKERRQ(ierr);
1090     }
1091     nnz = lnk[0];
1092 
1093     /* If free space is not available, double the total space in the list */
1094     if (current_space->local_remaining<nnz) {
1095       ierr = PetscFreeSpaceGet(nnz+current_space->total_array_size,&current_space);CHKERRQ(ierr);
1096       nspacedouble++;
1097     }
1098 
1099     /* Copy data into free space, and zero out denserows */
1100     ierr = PetscLLCondensedClean(aN,nnz,current_space->array,lnk,lnkbt);CHKERRQ(ierr);
1101     current_space->array           += nnz;
1102     current_space->local_used      += nnz;
1103     current_space->local_remaining -= nnz;
1104     coi[i+1] = coi[i] + nnz;
1105   }
1106 
1107   ierr = PetscMalloc((coi[pon]+1)*sizeof(PetscInt),&coj);CHKERRQ(ierr);
1108   ierr = PetscFreeSpaceContiguous(&free_space,coj);CHKERRQ(ierr);
1109   afill_tmp = (PetscReal)coi[pon]/(poti[pon] + ai[am]+1);
1110   if (afill_tmp > afill) afill = afill_tmp;
1111 
1112   /* send j-array (coj) of Co to other processors */
1113   /*----------------------------------------------*/
1114   /* determine row ownership */
1115   ierr = PetscNew(Mat_Merge_SeqsToMPI,&merge);CHKERRQ(ierr);
1116   ierr = PetscLayoutCreate(comm,&merge->rowmap);CHKERRQ(ierr);
1117   merge->rowmap->n = pn;
1118   merge->rowmap->bs = 1;
1119   ierr = PetscLayoutSetUp(merge->rowmap);CHKERRQ(ierr);
1120   owners = merge->rowmap->range;
1121 
1122   /* determine the number of messages to send, their lengths */
1123   ierr = PetscMalloc(size*sizeof(PetscMPIInt),&len_si);CHKERRQ(ierr);
1124   ierr = PetscMemzero(len_si,size*sizeof(PetscMPIInt));CHKERRQ(ierr);
1125   ierr = PetscMalloc(size*sizeof(PetscMPIInt),&merge->len_s);CHKERRQ(ierr);
1126   len_s = merge->len_s;
1127   merge->nsend = 0;
1128 
1129   ierr = PetscMalloc((size+2)*sizeof(PetscInt),&owners_co);CHKERRQ(ierr);
1130   ierr = PetscMemzero(len_s,size*sizeof(PetscMPIInt));CHKERRQ(ierr);
1131 
1132   proc = 0;
1133   for (i=0; i<pon; i++) {
1134     while (prmap[i] >= owners[proc+1]) proc++;
1135     len_si[proc]++;  /* num of rows in Co to be sent to [proc] */
1136     len_s[proc] += coi[i+1] - coi[i];
1137   }
1138 
1139   len   = 0;  /* max length of buf_si[] */
1140   owners_co[0] = 0;
1141   for (proc=0; proc<size; proc++) {
1142     owners_co[proc+1] = owners_co[proc] + len_si[proc];
1143     if (len_si[proc]) {
1144       merge->nsend++;
1145       len_si[proc] = 2*(len_si[proc] + 1);
1146       len += len_si[proc];
1147     }
1148   }
1149 
1150   /* determine the number and length of messages to receive for coi and coj  */
1151   ierr = PetscGatherNumberOfMessages(comm,PETSC_NULL,len_s,&merge->nrecv);CHKERRQ(ierr);
1152   ierr = PetscGatherMessageLengths2(comm,merge->nsend,merge->nrecv,len_s,len_si,&merge->id_r,&merge->len_r,&len_ri);CHKERRQ(ierr);
1153 
1154   /* post the Irecv and Isend of coj */
1155   ierr = PetscCommGetNewTag(comm,&tagj);CHKERRQ(ierr);
1156   ierr = PetscPostIrecvInt(comm,tagj,merge->nrecv,merge->id_r,merge->len_r,&buf_rj,&rwaits);CHKERRQ(ierr);
1157   ierr = PetscMalloc((merge->nsend+1)*sizeof(MPI_Request),&swaits);CHKERRQ(ierr);
1158   for (proc=0, k=0; proc<size; proc++) {
1159     if (!len_s[proc]) continue;
1160     i = owners_co[proc];
1161     ierr = MPI_Isend(coj+coi[i],len_s[proc],MPIU_INT,proc,tagj,comm,swaits+k);CHKERRQ(ierr);
1162     k++;
1163   }
1164 
1165   /* receives and sends of coj are complete */
1166   ierr = PetscMalloc(size*sizeof(MPI_Status),&sstatus);CHKERRQ(ierr);
1167   for (i=0; i<merge->nrecv; i++) {
1168     PetscMPIInt icompleted;
1169     ierr = MPI_Waitany(merge->nrecv,rwaits,&icompleted,&rstatus);CHKERRQ(ierr);
1170   }
1171   ierr = PetscFree(rwaits);CHKERRQ(ierr);
1172   if (merge->nsend) {ierr = MPI_Waitall(merge->nsend,swaits,sstatus);CHKERRQ(ierr);}
1173 
1174   /* send and recv coi */
1175   /*-------------------*/
1176   ierr = PetscCommGetNewTag(comm,&tagi);CHKERRQ(ierr);
1177   ierr = PetscPostIrecvInt(comm,tagi,merge->nrecv,merge->id_r,len_ri,&buf_ri,&rwaits);CHKERRQ(ierr);
1178   ierr = PetscMalloc((len+1)*sizeof(PetscInt),&buf_s);CHKERRQ(ierr);
1179   buf_si = buf_s;  /* points to the beginning of k-th msg to be sent */
1180   for (proc=0,k=0; proc<size; proc++) {
1181     if (!len_s[proc]) continue;
1182     /* form outgoing message for i-structure:
1183          buf_si[0]:                 nrows to be sent
1184                [1:nrows]:           row index (global)
1185                [nrows+1:2*nrows+1]: i-structure index
1186     */
1187     /*-------------------------------------------*/
1188     nrows = len_si[proc]/2 - 1;
1189     buf_si_i    = buf_si + nrows+1;
1190     buf_si[0]   = nrows;
1191     buf_si_i[0] = 0;
1192     nrows = 0;
1193     for (i=owners_co[proc]; i<owners_co[proc+1]; i++) {
1194       nzi = coi[i+1] - coi[i];
1195       buf_si_i[nrows+1] = buf_si_i[nrows] + nzi; /* i-structure */
1196       buf_si[nrows+1] =prmap[i] -owners[proc]; /* local row index */
1197       nrows++;
1198     }
1199     ierr = MPI_Isend(buf_si,len_si[proc],MPIU_INT,proc,tagi,comm,swaits+k);CHKERRQ(ierr);
1200     k++;
1201     buf_si += len_si[proc];
1202   }
1203   i = merge->nrecv;
1204   while (i--) {
1205     PetscMPIInt icompleted;
1206     ierr = MPI_Waitany(merge->nrecv,rwaits,&icompleted,&rstatus);CHKERRQ(ierr);
1207   }
1208   ierr = PetscFree(rwaits);CHKERRQ(ierr);
1209   if (merge->nsend) {ierr = MPI_Waitall(merge->nsend,swaits,sstatus);CHKERRQ(ierr);}
1210   ierr = PetscFree(len_si);CHKERRQ(ierr);
1211   ierr = PetscFree(len_ri);CHKERRQ(ierr);
1212   ierr = PetscFree(swaits);CHKERRQ(ierr);
1213   ierr = PetscFree(sstatus);CHKERRQ(ierr);
1214   ierr = PetscFree(buf_s);CHKERRQ(ierr);
1215 
1216   /* compute the local portion of C (mpi mat) */
1217   /*------------------------------------------*/
1218   /* allocate bi array and free space for accumulating nonzero column info */
1219   ierr = PetscMalloc((pn+1)*sizeof(PetscInt),&bi);CHKERRQ(ierr);
1220   bi[0] = 0;
1221 
1222   /* set initial free space to be fill*(nnz(P) + nnz(A)) */
1223   nnz           = fill*(pdti[pn] + poti[pon] + ai[am]);
1224   ierr          = PetscFreeSpaceGet(nnz,&free_space);CHKERRQ(ierr);
1225   current_space = free_space;
1226 
1227   ierr = PetscMalloc3(merge->nrecv,PetscInt**,&buf_ri_k,merge->nrecv,PetscInt*,&nextrow,merge->nrecv,PetscInt*,&nextci);CHKERRQ(ierr);
1228   for (k=0; k<merge->nrecv; k++) {
1229     buf_ri_k[k] = buf_ri[k]; /* beginning of k-th recved i-structure */
1230     nrows       = *buf_ri_k[k];
1231     nextrow[k]  = buf_ri_k[k] + 1;  /* next row number of k-th recved i-structure */
1232     nextci[k]   = buf_ri_k[k] + (nrows + 1);/* poins to the next i-structure of k-th recved i-structure  */
1233   }
1234 
1235   ierr = MatPreallocateInitialize(comm,pn,A->cmap->n,dnz,onz);CHKERRQ(ierr);
1236   rmax = 0;
1237   for (i=0; i<pn; i++) {
1238     /* add pdt[i,:]*AP into lnk */
1239     pnz = pdti[i+1] - pdti[i];
1240     ptJ = pdtj + pdti[i];
1241     for (j=0; j<pnz; j++) {
1242       row  = ptJ[j];  /* row of AP == col of Pt */
1243       anz  = ai[row+1] - ai[row];
1244       Jptr = aj + ai[row];
1245       /* add non-zero cols of AP into the sorted linked list lnk */
1246       ierr = PetscLLCondensedAddSorted(anz,Jptr,lnk,lnkbt);CHKERRQ(ierr);
1247     }
1248 
1249     /* add received col data into lnk */
1250     for (k=0; k<merge->nrecv; k++) { /* k-th received message */
1251       if (i == *nextrow[k]) { /* i-th row */
1252         nzi = *(nextci[k]+1) - *nextci[k];
1253         Jptr  = buf_rj[k] + *nextci[k];
1254         ierr = PetscLLCondensedAddSorted(nzi,Jptr,lnk,lnkbt);CHKERRQ(ierr);
1255         nextrow[k]++; nextci[k]++;
1256       }
1257     }
1258     nnz = lnk[0];
1259 
1260     /* if free space is not available, make more free space */
1261     if (current_space->local_remaining<nnz) {
1262       ierr = PetscFreeSpaceGet(nnz+current_space->total_array_size,&current_space);CHKERRQ(ierr);
1263       nspacedouble++;
1264     }
1265     /* copy data into free space, then initialize lnk */
1266     ierr = PetscLLCondensedClean(aN,nnz,current_space->array,lnk,lnkbt);CHKERRQ(ierr);
1267     ierr = MatPreallocateSet(i+owners[rank],nnz,current_space->array,dnz,onz);CHKERRQ(ierr);
1268     current_space->array           += nnz;
1269     current_space->local_used      += nnz;
1270     current_space->local_remaining -= nnz;
1271     bi[i+1] = bi[i] + nnz;
1272     if (nnz > rmax) rmax = nnz;
1273   }
1274   ierr = PetscFree3(buf_ri_k,nextrow,nextci);CHKERRQ(ierr);
1275 
1276   ierr = PetscMalloc((bi[pn]+1)*sizeof(PetscInt),&bj);CHKERRQ(ierr);
1277   ierr = PetscFreeSpaceContiguous(&free_space,bj);CHKERRQ(ierr);
1278   afill_tmp = (PetscReal)bi[pn]/(pdti[pn] + poti[pon] + ai[am]+1);
1279   if (afill_tmp > afill) afill = afill_tmp;
1280   ierr = PetscLLCondensedDestroy(lnk,lnkbt);CHKERRQ(ierr);
1281   ierr = MatDestroy(&POt);CHKERRQ(ierr);
1282   ierr = MatDestroy(&PDt);CHKERRQ(ierr);
1283 
1284   /* create symbolic parallel matrix Cmpi - why cannot be assembled in Numeric part   */
1285   /*----------------------------------------------------------------------------------*/
1286   ierr = PetscMalloc((rmax+1)*sizeof(PetscScalar),&vals);CHKERRQ(ierr);
1287   ierr = PetscMemzero(vals,rmax*sizeof(PetscScalar));CHKERRQ(ierr);
1288 
1289   ierr = MatCreate(comm,&Cmpi);CHKERRQ(ierr);
1290   ierr = MatSetSizes(Cmpi,pn,A->cmap->n,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr);
1291   ierr = MatSetBlockSizes(Cmpi,P->cmap->bs,A->cmap->bs);CHKERRQ(ierr);
1292   ierr = MatSetType(Cmpi,MATMPIAIJ);CHKERRQ(ierr);
1293   ierr = MatMPIAIJSetPreallocation(Cmpi,0,dnz,0,onz);CHKERRQ(ierr);
1294   ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr);
1295   ierr = MatSetBlockSize(Cmpi,1);CHKERRQ(ierr);
1296   for (i=0; i<pn; i++) {
1297     row = i + rstart;
1298     nnz = bi[i+1] - bi[i];
1299     Jptr = bj + bi[i];
1300     ierr = MatSetValues(Cmpi,1,&row,nnz,Jptr,vals,INSERT_VALUES);CHKERRQ(ierr);
1301   }
1302   ierr = MatAssemblyBegin(Cmpi,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1303   ierr = MatAssemblyEnd(Cmpi,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1304   ierr = PetscFree(vals);CHKERRQ(ierr);
1305 
1306   merge->bi            = bi;
1307   merge->bj            = bj;
1308   merge->coi           = coi;
1309   merge->coj           = coj;
1310   merge->buf_ri        = buf_ri;
1311   merge->buf_rj        = buf_rj;
1312   merge->owners_co     = owners_co;
1313   merge->destroy       = Cmpi->ops->destroy;
1314   merge->duplicate     = Cmpi->ops->duplicate;
1315 
1316   Cmpi->ops->mattransposemultnumeric = MatTransposeMatMultNumeric_MPIAIJ_MPIAIJ;
1317   Cmpi->ops->destroy                 = MatDestroy_MPIAIJ_PtAP;
1318 
1319   /* attach the supporting struct to Cmpi for reuse */
1320   c = (Mat_MPIAIJ*)Cmpi->data;
1321   c->ptap        = ptap;
1322   ptap->api      = PETSC_NULL;
1323   ptap->apj      = PETSC_NULL;
1324   ptap->merge    = merge;
1325   ptap->rmax     = rmax;
1326 
1327   *C = Cmpi;
1328 #if defined(PETSC_USE_INFO)
1329   if (bi[pn] != 0) {
1330     ierr = PetscInfo3(Cmpi,"Reallocs %D; Fill ratio: given %G needed %G.\n",nspacedouble,fill,afill);CHKERRQ(ierr);
1331     ierr = PetscInfo1(Cmpi,"Use MatTransposeMatMult(A,B,MatReuse,%G,&C) for best performance.\n",afill);CHKERRQ(ierr);
1332   } else {
1333     ierr = PetscInfo(Cmpi,"Empty matrix product\n");CHKERRQ(ierr);
1334   }
1335 #endif
1336   PetscFunctionReturn(0);
1337 }
1338 
1339 #undef __FUNCT__
1340 #define __FUNCT__ "MatTransposeMatMultNumeric_MPIAIJ_MPIAIJ_Scalable"
1341 PetscErrorCode MatTransposeMatMultNumeric_MPIAIJ_MPIAIJ_Scalable(Mat P,Mat A,Mat C)
1342 {
1343   PetscErrorCode       ierr;
1344   Mat_Merge_SeqsToMPI  *merge;
1345   Mat_MPIAIJ           *p=(Mat_MPIAIJ*)P->data,*c=(Mat_MPIAIJ*)C->data;
1346   Mat_SeqAIJ           *pd=(Mat_SeqAIJ*)(p->A)->data,*po=(Mat_SeqAIJ*)(p->B)->data;
1347   Mat_PtAPMPI          *ptap;
1348   PetscInt             *adj;
1349   PetscInt             i,j,k,anz,pnz,row,*cj,nexta;
1350   MatScalar            *ada,*ca,valtmp;
1351   PetscInt             am=A->rmap->n,cm=C->rmap->n,pon=(p->B)->cmap->n;
1352   MPI_Comm             comm=((PetscObject)C)->comm;
1353   PetscMPIInt          size,rank,taga,*len_s;
1354   PetscInt             *owners,proc,nrows,**buf_ri_k,**nextrow,**nextci;
1355   PetscInt             **buf_ri,**buf_rj;
1356   PetscInt             cnz=0,*bj_i,*bi,*bj,bnz,nextcj; /* bi,bj,ba: local array of C(mpi mat) */
1357   MPI_Request          *s_waits,*r_waits;
1358   MPI_Status           *status;
1359   MatScalar            **abuf_r,*ba_i,*pA,*coa,*ba;
1360   PetscInt             *ai,*aj,*coi,*coj;
1361   PetscInt             *poJ,*pdJ;
1362   Mat                  A_loc;
1363   Mat_SeqAIJ           *a_loc;
1364 
1365   PetscFunctionBegin;
1366   ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
1367   ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
1368 
1369   ptap  = c->ptap;
1370   merge = ptap->merge;
1371 
1372   /* 2) compute numeric C_seq = P_loc^T*A_loc */
1373   /*------------------------------------------*/
1374   /* get data from symbolic products */
1375   coi = merge->coi; coj = merge->coj;
1376   ierr = PetscMalloc((coi[pon]+1)*sizeof(MatScalar),&coa);CHKERRQ(ierr);
1377   ierr = PetscMemzero(coa,coi[pon]*sizeof(MatScalar));CHKERRQ(ierr);
1378   bi     = merge->bi; bj = merge->bj;
1379   owners = merge->rowmap->range;
1380   ierr   = PetscMalloc((bi[cm]+1)*sizeof(MatScalar),&ba);CHKERRQ(ierr);
1381   ierr   = PetscMemzero(ba,bi[cm]*sizeof(MatScalar));CHKERRQ(ierr);
1382 
1383   /* get A_loc by taking all local rows of A */
1384   A_loc = ptap->A_loc;
1385   ierr = MatMPIAIJGetLocalMat(A,MAT_REUSE_MATRIX,&A_loc);CHKERRQ(ierr);
1386   a_loc = (Mat_SeqAIJ*)(A_loc)->data;
1387   ai   = a_loc->i;
1388   aj   = a_loc->j;
1389 
1390   for (i=0; i<am; i++) {
1391     anz = ai[i+1] - ai[i];
1392     adj = aj + ai[i];
1393     ada = a_loc->a + ai[i];
1394 
1395     /* 2-b) Compute Cseq = P_loc[i,:]^T*A[i,:] using outer product */
1396     /*-------------------------------------------------------------*/
1397     /* put the value into Co=(p->B)^T*A (off-diagonal part, send to others) */
1398     pnz = po->i[i+1] - po->i[i];
1399     poJ = po->j + po->i[i];
1400     pA  = po->a + po->i[i];
1401     for (j=0; j<pnz; j++) {
1402       row = poJ[j];
1403       cj  = coj + coi[row];
1404       ca  = coa + coi[row];
1405       /* perform sparse axpy */
1406       nexta  = 0;
1407       valtmp = pA[j];
1408       for (k=0; nexta<anz; k++) {
1409         if (cj[k] == adj[nexta]) {
1410           ca[k] += valtmp*ada[nexta];
1411           nexta++;
1412         }
1413       }
1414       ierr = PetscLogFlops(2.0*anz);CHKERRQ(ierr);
1415     }
1416 
1417     /* put the value into Cd (diagonal part) */
1418     pnz = pd->i[i+1] - pd->i[i];
1419     pdJ = pd->j + pd->i[i];
1420     pA  = pd->a + pd->i[i];
1421     for (j=0; j<pnz; j++) {
1422       row = pdJ[j];
1423       cj  = bj + bi[row];
1424       ca  = ba + bi[row];
1425       /* perform sparse axpy */
1426       nexta  = 0;
1427       valtmp = pA[j];
1428       for (k=0; nexta<anz; k++) {
1429         if (cj[k] == adj[nexta]) {
1430           ca[k] += valtmp*ada[nexta];
1431           nexta++;
1432         }
1433       }
1434       ierr = PetscLogFlops(2.0*anz);CHKERRQ(ierr);
1435     }
1436   }
1437 
1438   /* 3) send and recv matrix values coa */
1439   /*------------------------------------*/
1440   buf_ri = merge->buf_ri;
1441   buf_rj = merge->buf_rj;
1442   len_s  = merge->len_s;
1443   ierr = PetscCommGetNewTag(comm,&taga);CHKERRQ(ierr);
1444   ierr = PetscPostIrecvScalar(comm,taga,merge->nrecv,merge->id_r,merge->len_r,&abuf_r,&r_waits);CHKERRQ(ierr);
1445 
1446   ierr = PetscMalloc2(merge->nsend+1,MPI_Request,&s_waits,size,MPI_Status,&status);CHKERRQ(ierr);
1447   for (proc=0,k=0; proc<size; proc++) {
1448     if (!len_s[proc]) continue;
1449     i = merge->owners_co[proc];
1450     ierr = MPI_Isend(coa+coi[i],len_s[proc],MPIU_MATSCALAR,proc,taga,comm,s_waits+k);CHKERRQ(ierr);
1451     k++;
1452   }
1453   if (merge->nrecv) {ierr = MPI_Waitall(merge->nrecv,r_waits,status);CHKERRQ(ierr);}
1454   if (merge->nsend) {ierr = MPI_Waitall(merge->nsend,s_waits,status);CHKERRQ(ierr);}
1455 
1456   ierr = PetscFree2(s_waits,status);CHKERRQ(ierr);
1457   ierr = PetscFree(r_waits);CHKERRQ(ierr);
1458   ierr = PetscFree(coa);CHKERRQ(ierr);
1459 
1460   /* 4) insert local Cseq and received values into Cmpi */
1461   /*----------------------------------------------------*/
1462   ierr = PetscMalloc3(merge->nrecv,PetscInt**,&buf_ri_k,merge->nrecv,PetscInt*,&nextrow,merge->nrecv,PetscInt*,&nextci);CHKERRQ(ierr);
1463   for (k=0; k<merge->nrecv; k++) {
1464     buf_ri_k[k] = buf_ri[k]; /* beginning of k-th recved i-structure */
1465     nrows       = *(buf_ri_k[k]);
1466     nextrow[k]  = buf_ri_k[k]+1;  /* next row number of k-th recved i-structure */
1467     nextci[k]   = buf_ri_k[k] + (nrows + 1);/* poins to the next i-structure of k-th recved i-structure  */
1468   }
1469 
1470   for (i=0; i<cm; i++) {
1471     row = owners[rank] + i; /* global row index of C_seq */
1472     bj_i = bj + bi[i];  /* col indices of the i-th row of C */
1473     ba_i = ba + bi[i];
1474     bnz  = bi[i+1] - bi[i];
1475     /* add received vals into ba */
1476     for (k=0; k<merge->nrecv; k++) { /* k-th received message */
1477       /* i-th row */
1478       if (i == *nextrow[k]) {
1479         cnz = *(nextci[k]+1) - *nextci[k];
1480         cj  = buf_rj[k] + *(nextci[k]);
1481         ca  = abuf_r[k] + *(nextci[k]);
1482         nextcj = 0;
1483         for (j=0; nextcj<cnz; j++) {
1484           if (bj_i[j] == cj[nextcj]) { /* bcol == ccol */
1485             ba_i[j] += ca[nextcj++];
1486           }
1487         }
1488         nextrow[k]++; nextci[k]++;
1489         ierr = PetscLogFlops(2.0*cnz);CHKERRQ(ierr);
1490       }
1491     }
1492     ierr = MatSetValues(C,1,&row,bnz,bj_i,ba_i,INSERT_VALUES);CHKERRQ(ierr);
1493   }
1494   ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1495   ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1496 
1497   ierr = PetscFree(ba);CHKERRQ(ierr);
1498   ierr = PetscFree(abuf_r[0]);CHKERRQ(ierr);
1499   ierr = PetscFree(abuf_r);CHKERRQ(ierr);
1500   ierr = PetscFree3(buf_ri_k,nextrow,nextci);CHKERRQ(ierr);
1501   PetscFunctionReturn(0);
1502 }
1503 
1504 /* This routine is modified from MatPtAPSymbolic_MPIAIJ_MPIAIJ() */
1505 #undef __FUNCT__
1506 #define __FUNCT__ "MatTransposeMatMultSymbolic_MPIAIJ_MPIAIJ_Scalable"
1507 PetscErrorCode MatTransposeMatMultSymbolic_MPIAIJ_MPIAIJ_Scalable(Mat P,Mat A,PetscReal fill,Mat *C)
1508 {
1509   PetscErrorCode       ierr;
1510   Mat                  Cmpi,A_loc,POt,PDt;
1511   Mat_PtAPMPI          *ptap;
1512   PetscFreeSpaceList   free_space=PETSC_NULL,current_space=PETSC_NULL;
1513   Mat_MPIAIJ           *p=(Mat_MPIAIJ*)P->data,*c;
1514   PetscInt             *pdti,*pdtj,*poti,*potj,*ptJ;
1515   PetscInt             nnz;
1516   PetscInt             *lnk,*owners_co,*coi,*coj,i,k,pnz,row;
1517   PetscInt             am=A->rmap->n,pn=P->cmap->n;
1518   MPI_Comm             comm=((PetscObject)A)->comm;
1519   PetscMPIInt          size,rank,tagi,tagj,*len_si,*len_s,*len_ri;
1520   PetscInt             **buf_rj,**buf_ri,**buf_ri_k;
1521   PetscInt             len,proc,*dnz,*onz,*owners;
1522   PetscInt             nzi,*bi,*bj;
1523   PetscInt             nrows,*buf_s,*buf_si,*buf_si_i,**nextrow,**nextci;
1524   MPI_Request          *swaits,*rwaits;
1525   MPI_Status           *sstatus,rstatus;
1526   Mat_Merge_SeqsToMPI  *merge;
1527   PetscInt             *ai,*aj,*Jptr,anz,*prmap=p->garray,pon,nspacedouble=0,j;
1528   PetscReal            afill=1.0,afill_tmp;
1529   PetscInt             rstart = P->cmap->rstart,rmax,aN=A->cmap->N,Crmax;
1530   PetscScalar          *vals;
1531   Mat_SeqAIJ           *a_loc, *pdt,*pot;
1532 
1533   PetscFunctionBegin;
1534   /* check if matrix local sizes are compatible */
1535   if (A->rmap->rstart != P->rmap->rstart || A->rmap->rend != P->rmap->rend) {
1536     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);
1537   }
1538 
1539   ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
1540   ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
1541 
1542   /* create struct Mat_PtAPMPI and attached it to C later */
1543   ierr = PetscNew(Mat_PtAPMPI,&ptap);CHKERRQ(ierr);
1544 
1545   /* get A_loc by taking all local rows of A */
1546   ierr = MatMPIAIJGetLocalMat(A,MAT_INITIAL_MATRIX,&A_loc);CHKERRQ(ierr);
1547   ptap->A_loc = A_loc;
1548   a_loc = (Mat_SeqAIJ*)(A_loc)->data;
1549   ai   = a_loc->i;
1550   aj   = a_loc->j;
1551 
1552   /* determine symbolic Co=(p->B)^T*A - send to others */
1553   /*----------------------------------------------------*/
1554   ierr = MatTransposeSymbolic_SeqAIJ(p->A,&PDt);CHKERRQ(ierr);
1555   pdt = (Mat_SeqAIJ*)PDt->data;
1556   pdti = pdt->i; pdtj = pdt->j;
1557 
1558   ierr = MatTransposeSymbolic_SeqAIJ(p->B,&POt);CHKERRQ(ierr);
1559   pot = (Mat_SeqAIJ*)POt->data;
1560   poti = pot->i; potj = pot->j;
1561 
1562   /* then, compute symbolic Co = (p->B)^T*A */
1563   pon = (p->B)->cmap->n; /* total num of rows to be sent to other processors
1564                          >= (num of nonzero rows of C_seq) - pn */
1565   ierr = PetscMalloc((pon+1)*sizeof(PetscInt),&coi);CHKERRQ(ierr);
1566   coi[0] = 0;
1567 
1568   /* set initial free space to be fill*(nnz(p->B) + nnz(A)) */
1569   nnz           = fill*(poti[pon] + ai[am]);
1570   ierr          = PetscFreeSpaceGet(nnz,&free_space);CHKERRQ(ierr);
1571   current_space = free_space;
1572 
1573   /* create and initialize a linked list */
1574   i = PetscMax(pdt->rmax,pot->rmax);
1575   Crmax = i*a_loc->rmax*size; /* non-scalable! */
1576   if (!Crmax || Crmax > aN) Crmax = aN;
1577   ierr = PetscLLCondensedCreate_Scalable(Crmax,&lnk);CHKERRQ(ierr);
1578 
1579   for (i=0; i<pon; i++) {
1580     pnz = poti[i+1] - poti[i];
1581     ptJ = potj + poti[i];
1582     for (j=0; j<pnz; j++) {
1583       row  = ptJ[j]; /* row of A_loc == col of Pot */
1584       anz  = ai[row+1] - ai[row];
1585       Jptr = aj + ai[row];
1586       /* add non-zero cols of AP into the sorted linked list lnk */
1587       ierr = PetscLLCondensedAddSorted_Scalable(anz,Jptr,lnk);CHKERRQ(ierr);
1588     }
1589     nnz = lnk[0];
1590 
1591     /* If free space is not available, double the total space in the list */
1592     if (current_space->local_remaining<nnz) {
1593       ierr = PetscFreeSpaceGet(nnz+current_space->total_array_size,&current_space);CHKERRQ(ierr);
1594       nspacedouble++;
1595     }
1596 
1597     /* Copy data into free space, and zero out denserows */
1598     ierr = PetscLLCondensedClean_Scalable(nnz,current_space->array,lnk);CHKERRQ(ierr);
1599     current_space->array           += nnz;
1600     current_space->local_used      += nnz;
1601     current_space->local_remaining -= nnz;
1602     coi[i+1] = coi[i] + nnz;
1603   }
1604 
1605   ierr = PetscMalloc((coi[pon]+1)*sizeof(PetscInt),&coj);CHKERRQ(ierr);
1606   ierr = PetscFreeSpaceContiguous(&free_space,coj);CHKERRQ(ierr);
1607   afill_tmp = (PetscReal)coi[pon]/(poti[pon] + ai[am]+1);
1608   if (afill_tmp > afill) afill = afill_tmp;
1609 
1610   /* send j-array (coj) of Co to other processors */
1611   /*----------------------------------------------*/
1612   /* determine row ownership */
1613   ierr = PetscNew(Mat_Merge_SeqsToMPI,&merge);CHKERRQ(ierr);
1614   ierr = PetscLayoutCreate(comm,&merge->rowmap);CHKERRQ(ierr);
1615   merge->rowmap->n = pn;
1616   merge->rowmap->bs = 1;
1617   ierr = PetscLayoutSetUp(merge->rowmap);CHKERRQ(ierr);
1618   owners = merge->rowmap->range;
1619 
1620   /* determine the number of messages to send, their lengths */
1621   ierr = PetscMalloc(size*sizeof(PetscMPIInt),&len_si);CHKERRQ(ierr);
1622   ierr = PetscMemzero(len_si,size*sizeof(PetscMPIInt));CHKERRQ(ierr);
1623   ierr = PetscMalloc(size*sizeof(PetscMPIInt),&merge->len_s);CHKERRQ(ierr);
1624   len_s = merge->len_s;
1625   merge->nsend = 0;
1626 
1627   ierr = PetscMalloc((size+2)*sizeof(PetscInt),&owners_co);CHKERRQ(ierr);
1628   ierr = PetscMemzero(len_s,size*sizeof(PetscMPIInt));CHKERRQ(ierr);
1629 
1630   proc = 0;
1631   for (i=0; i<pon; i++) {
1632     while (prmap[i] >= owners[proc+1]) proc++;
1633     len_si[proc]++;  /* num of rows in Co to be sent to [proc] */
1634     len_s[proc] += coi[i+1] - coi[i];
1635   }
1636 
1637   len   = 0;  /* max length of buf_si[] */
1638   owners_co[0] = 0;
1639   for (proc=0; proc<size; proc++) {
1640     owners_co[proc+1] = owners_co[proc] + len_si[proc];
1641     if (len_si[proc]) {
1642       merge->nsend++;
1643       len_si[proc] = 2*(len_si[proc] + 1);
1644       len += len_si[proc];
1645     }
1646   }
1647 
1648   /* determine the number and length of messages to receive for coi and coj  */
1649   ierr = PetscGatherNumberOfMessages(comm,PETSC_NULL,len_s,&merge->nrecv);CHKERRQ(ierr);
1650   ierr = PetscGatherMessageLengths2(comm,merge->nsend,merge->nrecv,len_s,len_si,&merge->id_r,&merge->len_r,&len_ri);CHKERRQ(ierr);
1651 
1652   /* post the Irecv and Isend of coj */
1653   ierr = PetscCommGetNewTag(comm,&tagj);CHKERRQ(ierr);
1654   ierr = PetscPostIrecvInt(comm,tagj,merge->nrecv,merge->id_r,merge->len_r,&buf_rj,&rwaits);CHKERRQ(ierr);
1655   ierr = PetscMalloc((merge->nsend+1)*sizeof(MPI_Request),&swaits);CHKERRQ(ierr);
1656   for (proc=0, k=0; proc<size; proc++) {
1657     if (!len_s[proc]) continue;
1658     i = owners_co[proc];
1659     ierr = MPI_Isend(coj+coi[i],len_s[proc],MPIU_INT,proc,tagj,comm,swaits+k);CHKERRQ(ierr);
1660     k++;
1661   }
1662 
1663   /* receives and sends of coj are complete */
1664   ierr = PetscMalloc(size*sizeof(MPI_Status),&sstatus);CHKERRQ(ierr);
1665   for (i=0; i<merge->nrecv; i++) {
1666     PetscMPIInt icompleted;
1667     ierr = MPI_Waitany(merge->nrecv,rwaits,&icompleted,&rstatus);CHKERRQ(ierr);
1668   }
1669   ierr = PetscFree(rwaits);CHKERRQ(ierr);
1670   if (merge->nsend) {ierr = MPI_Waitall(merge->nsend,swaits,sstatus);CHKERRQ(ierr);}
1671 
1672   /* send and recv coi */
1673   /*-------------------*/
1674   ierr = PetscCommGetNewTag(comm,&tagi);CHKERRQ(ierr);
1675   ierr = PetscPostIrecvInt(comm,tagi,merge->nrecv,merge->id_r,len_ri,&buf_ri,&rwaits);CHKERRQ(ierr);
1676   ierr = PetscMalloc((len+1)*sizeof(PetscInt),&buf_s);CHKERRQ(ierr);
1677   buf_si = buf_s;  /* points to the beginning of k-th msg to be sent */
1678   for (proc=0,k=0; proc<size; proc++) {
1679     if (!len_s[proc]) continue;
1680     /* form outgoing message for i-structure:
1681          buf_si[0]:                 nrows to be sent
1682                [1:nrows]:           row index (global)
1683                [nrows+1:2*nrows+1]: i-structure index
1684     */
1685     /*-------------------------------------------*/
1686     nrows = len_si[proc]/2 - 1;
1687     buf_si_i    = buf_si + nrows+1;
1688     buf_si[0]   = nrows;
1689     buf_si_i[0] = 0;
1690     nrows = 0;
1691     for (i=owners_co[proc]; i<owners_co[proc+1]; i++) {
1692       nzi = coi[i+1] - coi[i];
1693       buf_si_i[nrows+1] = buf_si_i[nrows] + nzi; /* i-structure */
1694       buf_si[nrows+1] =prmap[i] -owners[proc]; /* local row index */
1695       nrows++;
1696     }
1697     ierr = MPI_Isend(buf_si,len_si[proc],MPIU_INT,proc,tagi,comm,swaits+k);CHKERRQ(ierr);
1698     k++;
1699     buf_si += len_si[proc];
1700   }
1701   i = merge->nrecv;
1702   while (i--) {
1703     PetscMPIInt icompleted;
1704     ierr = MPI_Waitany(merge->nrecv,rwaits,&icompleted,&rstatus);CHKERRQ(ierr);
1705   }
1706   ierr = PetscFree(rwaits);CHKERRQ(ierr);
1707   if (merge->nsend) {ierr = MPI_Waitall(merge->nsend,swaits,sstatus);CHKERRQ(ierr);}
1708   ierr = PetscFree(len_si);CHKERRQ(ierr);
1709   ierr = PetscFree(len_ri);CHKERRQ(ierr);
1710   ierr = PetscFree(swaits);CHKERRQ(ierr);
1711   ierr = PetscFree(sstatus);CHKERRQ(ierr);
1712   ierr = PetscFree(buf_s);CHKERRQ(ierr);
1713 
1714   /* compute the local portion of C (mpi mat) */
1715   /*------------------------------------------*/
1716   /* allocate bi array and free space for accumulating nonzero column info */
1717   ierr = PetscMalloc((pn+1)*sizeof(PetscInt),&bi);CHKERRQ(ierr);
1718   bi[0] = 0;
1719 
1720   /* set initial free space to be fill*(nnz(P) + nnz(AP)) */
1721   nnz           = fill*(pdti[pn] + poti[pon] + ai[am]);
1722   ierr          = PetscFreeSpaceGet(nnz,&free_space);CHKERRQ(ierr);
1723   current_space = free_space;
1724 
1725   ierr = PetscMalloc3(merge->nrecv,PetscInt**,&buf_ri_k,merge->nrecv,PetscInt*,&nextrow,merge->nrecv,PetscInt*,&nextci);CHKERRQ(ierr);
1726   for (k=0; k<merge->nrecv; k++) {
1727     buf_ri_k[k] = buf_ri[k]; /* beginning of k-th recved i-structure */
1728     nrows       = *buf_ri_k[k];
1729     nextrow[k]  = buf_ri_k[k] + 1;  /* next row number of k-th recved i-structure */
1730     nextci[k]   = buf_ri_k[k] + (nrows + 1);/* poins to the next i-structure of k-th recved i-structure  */
1731   }
1732 
1733   ierr = MatPreallocateInitialize(comm,pn,A->cmap->n,dnz,onz);CHKERRQ(ierr);
1734   rmax = 0;
1735   for (i=0; i<pn; i++) {
1736     /* add pdt[i,:]*AP into lnk */
1737     pnz = pdti[i+1] - pdti[i];
1738     ptJ = pdtj + pdti[i];
1739     for (j=0; j<pnz; j++) {
1740       row  = ptJ[j];  /* row of AP == col of Pt */
1741       anz  = ai[row+1] - ai[row];
1742       Jptr = aj + ai[row];
1743       /* add non-zero cols of AP into the sorted linked list lnk */
1744       ierr = PetscLLCondensedAddSorted_Scalable(anz,Jptr,lnk);CHKERRQ(ierr);
1745     }
1746 
1747     /* add received col data into lnk */
1748     for (k=0; k<merge->nrecv; k++) { /* k-th received message */
1749       if (i == *nextrow[k]) { /* i-th row */
1750         nzi = *(nextci[k]+1) - *nextci[k];
1751         Jptr  = buf_rj[k] + *nextci[k];
1752         ierr = PetscLLCondensedAddSorted_Scalable(nzi,Jptr,lnk);CHKERRQ(ierr);
1753         nextrow[k]++; nextci[k]++;
1754       }
1755     }
1756     nnz = lnk[0];
1757 
1758     /* if free space is not available, make more free space */
1759     if (current_space->local_remaining<nnz) {
1760       ierr = PetscFreeSpaceGet(nnz+current_space->total_array_size,&current_space);CHKERRQ(ierr);
1761       nspacedouble++;
1762     }
1763     /* copy data into free space, then initialize lnk */
1764     ierr = PetscLLCondensedClean_Scalable(nnz,current_space->array,lnk);CHKERRQ(ierr);
1765     ierr = MatPreallocateSet(i+owners[rank],nnz,current_space->array,dnz,onz);CHKERRQ(ierr);
1766     current_space->array           += nnz;
1767     current_space->local_used      += nnz;
1768     current_space->local_remaining -= nnz;
1769     bi[i+1] = bi[i] + nnz;
1770     if (nnz > rmax) rmax = nnz;
1771   }
1772   ierr = PetscFree3(buf_ri_k,nextrow,nextci);CHKERRQ(ierr);
1773 
1774   ierr = PetscMalloc((bi[pn]+1)*sizeof(PetscInt),&bj);CHKERRQ(ierr);
1775   ierr = PetscFreeSpaceContiguous(&free_space,bj);CHKERRQ(ierr);
1776   afill_tmp = (PetscReal)bi[pn]/(pdti[pn] + poti[pon] + ai[am]+1);
1777   if (afill_tmp > afill) afill = afill_tmp;
1778   ierr = PetscLLCondensedDestroy_Scalable(lnk);CHKERRQ(ierr);
1779   ierr = MatDestroy(&POt);CHKERRQ(ierr);
1780   ierr = MatDestroy(&PDt);CHKERRQ(ierr);
1781 
1782   /* create symbolic parallel matrix Cmpi - why cannot be assembled in Numeric part   */
1783   /*----------------------------------------------------------------------------------*/
1784   ierr = PetscMalloc((rmax+1)*sizeof(PetscScalar),&vals);CHKERRQ(ierr);
1785   ierr = PetscMemzero(vals,rmax*sizeof(PetscScalar));CHKERRQ(ierr);
1786 
1787   ierr = MatCreate(comm,&Cmpi);CHKERRQ(ierr);
1788   ierr = MatSetSizes(Cmpi,pn,A->cmap->n,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr);
1789   ierr = MatSetBlockSizes(Cmpi,P->cmap->bs,A->cmap->bs);CHKERRQ(ierr);
1790   ierr = MatSetType(Cmpi,MATMPIAIJ);CHKERRQ(ierr);
1791   ierr = MatMPIAIJSetPreallocation(Cmpi,0,dnz,0,onz);CHKERRQ(ierr);
1792   ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr);
1793   ierr = MatSetBlockSize(Cmpi,1);CHKERRQ(ierr);
1794   for (i=0; i<pn; i++) {
1795     row = i + rstart;
1796     nnz = bi[i+1] - bi[i];
1797     Jptr = bj + bi[i];
1798     ierr = MatSetValues(Cmpi,1,&row,nnz,Jptr,vals,INSERT_VALUES);CHKERRQ(ierr);
1799   }
1800   ierr = MatAssemblyBegin(Cmpi,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1801   ierr = MatAssemblyEnd(Cmpi,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1802   ierr = PetscFree(vals);CHKERRQ(ierr);
1803 
1804   merge->bi            = bi;
1805   merge->bj            = bj;
1806   merge->coi           = coi;
1807   merge->coj           = coj;
1808   merge->buf_ri        = buf_ri;
1809   merge->buf_rj        = buf_rj;
1810   merge->owners_co     = owners_co;
1811   merge->destroy       = Cmpi->ops->destroy;
1812   merge->duplicate     = Cmpi->ops->duplicate;
1813 
1814   Cmpi->ops->mattransposemultnumeric = MatTransposeMatMultNumeric_MPIAIJ_MPIAIJ_Scalable;
1815   Cmpi->ops->destroy                 = MatDestroy_MPIAIJ_PtAP;
1816 
1817   /* attach the supporting struct to Cmpi for reuse */
1818   c = (Mat_MPIAIJ*)Cmpi->data;
1819   c->ptap        = ptap;
1820   ptap->api      = PETSC_NULL;
1821   ptap->apj      = PETSC_NULL;
1822   ptap->merge    = merge;
1823   ptap->rmax     = rmax;
1824   ptap->apa      = PETSC_NULL;
1825 
1826   *C = Cmpi;
1827 #if defined(PETSC_USE_INFO)
1828   if (bi[pn] != 0) {
1829     ierr = PetscInfo3(Cmpi,"Reallocs %D; Fill ratio: given %G needed %G.\n",nspacedouble,fill,afill);CHKERRQ(ierr);
1830     ierr = PetscInfo1(Cmpi,"Use MatTransposeMatMult(A,B,MatReuse,%G,&C) for best performance.\n",afill);CHKERRQ(ierr);
1831   } else {
1832     ierr = PetscInfo(Cmpi,"Empty matrix product\n");CHKERRQ(ierr);
1833   }
1834 #endif
1835   PetscFunctionReturn(0);
1836 }
1837