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