Actual source code: mpiptap.c

petsc-master 2019-07-18
Report Typos and Errors

  2: /*
  3:   Defines projective product routines where A is a MPIAIJ matrix
  4:           C = P^T * A * P
  5: */

  7:  #include <../src/mat/impls/aij/seq/aij.h>
  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:  #include <petsc/private/hashmapiv.h>
 13:  #include <petsc/private/hashseti.h>
 14:  #include <petscsf.h>


 17: PetscErrorCode MatView_MPIAIJ_PtAP(Mat A,PetscViewer viewer)
 18: {
 19:   PetscErrorCode    ierr;
 20:   Mat_MPIAIJ        *a=(Mat_MPIAIJ*)A->data;
 21:   Mat_APMPI         *ptap=a->ap;
 22:   PetscBool         iascii;
 23:   PetscViewerFormat format;

 26:   if (!ptap) {
 27:     /* hack: MatDuplicate() sets oldmat->ops->view to newmat which is a base mat class with null ptpa! */
 28:     A->ops->view = MatView_MPIAIJ;
 29:     (A->ops->view)(A,viewer);
 30:     return(0);
 31:   }

 33:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
 34:   if (iascii) {
 35:     PetscViewerGetFormat(viewer,&format);
 36:     if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
 37:       if (ptap->algType == 0) {
 38:         PetscViewerASCIIPrintf(viewer,"using scalable MatPtAP() implementation\n");
 39:       } else if (ptap->algType == 1) {
 40:         PetscViewerASCIIPrintf(viewer,"using nonscalable MatPtAP() implementation\n");
 41:       } else if (ptap->algType == 2) {
 42:         PetscViewerASCIIPrintf(viewer,"using allatonce MatPtAP() implementation\n");
 43:       } else if (ptap->algType == 3) {
 44:         PetscViewerASCIIPrintf(viewer,"using merged allatonce MatPtAP() implementation\n");
 45:       }
 46:     }
 47:   }
 48:   (ptap->view)(A,viewer);
 49:   return(0);
 50: }

 52: PetscErrorCode MatFreeIntermediateDataStructures_MPIAIJ_AP(Mat A)
 53: {
 54:   PetscErrorCode      ierr;
 55:   Mat_MPIAIJ          *a=(Mat_MPIAIJ*)A->data;
 56:   Mat_APMPI           *ptap=a->ap;
 57:   Mat_Merge_SeqsToMPI *merge;

 60:   if (!ptap) return(0);

 62:   PetscFree2(ptap->startsj_s,ptap->startsj_r);
 63:   PetscFree(ptap->bufa);
 64:   MatDestroy(&ptap->P_loc);
 65:   MatDestroy(&ptap->P_oth);
 66:   MatDestroy(&ptap->A_loc); /* used by MatTransposeMatMult() */
 67:   MatDestroy(&ptap->Rd);
 68:   MatDestroy(&ptap->Ro);
 69:   if (ptap->AP_loc) { /* used by alg_rap */
 70:     Mat_SeqAIJ *ap = (Mat_SeqAIJ*)(ptap->AP_loc)->data;
 71:     PetscFree(ap->i);
 72:     PetscFree2(ap->j,ap->a);
 73:     MatDestroy(&ptap->AP_loc);
 74:   } else { /* used by alg_ptap */
 75:     PetscFree(ptap->api);
 76:     PetscFree(ptap->apj);
 77:   }
 78:   MatDestroy(&ptap->C_loc);
 79:   MatDestroy(&ptap->C_oth);
 80:   if (ptap->apa) {PetscFree(ptap->apa);}

 82:   MatDestroy(&ptap->Pt);

 84:   merge=ptap->merge;
 85:   if (merge) { /* used by alg_ptap */
 86:     PetscFree(merge->id_r);
 87:     PetscFree(merge->len_s);
 88:     PetscFree(merge->len_r);
 89:     PetscFree(merge->bi);
 90:     PetscFree(merge->bj);
 91:     PetscFree(merge->buf_ri[0]);
 92:     PetscFree(merge->buf_ri);
 93:     PetscFree(merge->buf_rj[0]);
 94:     PetscFree(merge->buf_rj);
 95:     PetscFree(merge->coi);
 96:     PetscFree(merge->coj);
 97:     PetscFree(merge->owners_co);
 98:     PetscLayoutDestroy(&merge->rowmap);
 99:     PetscFree(ptap->merge);
100:   }
101:   ISLocalToGlobalMappingDestroy(&ptap->ltog);

103:   PetscSFDestroy(&ptap->sf);
104:   PetscFree(ptap->c_othi);
105:   PetscFree(ptap->c_rmti);
106:   return(0);
107: }

109: PetscErrorCode MatDestroy_MPIAIJ_PtAP(Mat A)
110: {
112:   Mat_MPIAIJ     *a=(Mat_MPIAIJ*)A->data;
113:   Mat_APMPI      *ptap=a->ap;

116:   (*A->ops->freeintermediatedatastructures)(A);
117:   (*ptap->destroy)(A); /* MatDestroy_MPIAIJ(A) */
118:   PetscFree(ptap);
119:   return(0);
120: }

122: PETSC_INTERN PetscErrorCode MatPtAP_MPIAIJ_MPIAIJ(Mat A,Mat P,MatReuse scall,PetscReal fill,Mat *C)
123: {
125:   PetscBool      flg;
126:   MPI_Comm       comm;
127: #if !defined(PETSC_HAVE_HYPRE)
128:   const char          *algTypes[4] = {"scalable","nonscalable","allatonce","allatonce_merged"};
129:   PetscInt            nalg=4;
130: #else
131:   const char          *algTypes[5] = {"scalable","nonscalable","allatonce","allatonce_merged","hypre"};
132:   PetscInt            nalg=5;
133: #endif
134:   PetscInt            pN=P->cmap->N,alg=1; /* set default algorithm */

137:   /* check if matrix local sizes are compatible */
138:   PetscObjectGetComm((PetscObject)A,&comm);
139:   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);
140:   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);

142:   if (scall == MAT_INITIAL_MATRIX) {
143:     /* pick an algorithm */
144:     PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MatPtAP","Mat");
145:     PetscOptionsEList("-matptap_via","Algorithmic approach","MatPtAP",algTypes,nalg,algTypes[alg],&alg,&flg);
146:     PetscOptionsEnd();

148:     if (!flg && pN > 100000) { /* may switch to scalable algorithm as default */
149:       MatInfo     Ainfo,Pinfo;
150:       PetscInt    nz_local;
151:       PetscBool   alg_scalable_loc=PETSC_FALSE,alg_scalable;

153:       MatGetInfo(A,MAT_LOCAL,&Ainfo);
154:       MatGetInfo(P,MAT_LOCAL,&Pinfo);
155:       nz_local = (PetscInt)(Ainfo.nz_allocated + Pinfo.nz_allocated);

157:       if (pN > fill*nz_local) alg_scalable_loc = PETSC_TRUE;
158:       MPIU_Allreduce(&alg_scalable_loc,&alg_scalable,1,MPIU_BOOL,MPI_LOR,comm);

160:       if (alg_scalable) {
161:         alg = 0; /* scalable algorithm would 50% slower than nonscalable algorithm */
162:       }
163:     }

165:     switch (alg) {
166:     case 1:
167:       /* do R=P^T locally, then C=R*A*P -- nonscalable */
168:       PetscLogEventBegin(MAT_PtAPSymbolic,A,P,0,0);
169:       MatPtAPSymbolic_MPIAIJ_MPIAIJ(A,P,fill,C);
170:       PetscLogEventEnd(MAT_PtAPSymbolic,A,P,0,0);
171:       break;
172:     case 2:
173:       /* compute C=P^T*A*P allatonce */
174:       PetscLogEventBegin(MAT_PtAPSymbolic,A,P,0,0);
175:       MatPtAPSymbolic_MPIAIJ_MPIAIJ_allatonce(A,P,fill,C);
176:       PetscLogEventEnd(MAT_PtAPSymbolic,A,P,0,0);
177:       break;
178:     case 3:
179:       /* compute C=P^T*A*P allatonce */
180:       PetscLogEventBegin(MAT_PtAPSymbolic,A,P,0,0);
181:       MatPtAPSymbolic_MPIAIJ_MPIAIJ_allatonce_merged(A,P,fill,C);
182:       PetscLogEventEnd(MAT_PtAPSymbolic,A,P,0,0);
183:       break;
184: #if defined(PETSC_HAVE_HYPRE)
185:     case 4:
186:       /* Use boomerAMGBuildCoarseOperator */
187:       PetscLogEventBegin(MAT_PtAPSymbolic,A,P,0,0);
188:       MatPtAPSymbolic_AIJ_AIJ_wHYPRE(A,P,fill,C);
189:       PetscLogEventEnd(MAT_PtAPSymbolic,A,P,0,0);
190:       break;
191: #endif
192:     default:
193:       /* do R=P^T locally, then C=R*A*P */
194:       PetscLogEventBegin(MAT_PtAPSymbolic,A,P,0,0);
195:       MatPtAPSymbolic_MPIAIJ_MPIAIJ_scalable(A,P,fill,C);
196:       PetscLogEventEnd(MAT_PtAPSymbolic,A,P,0,0);
197:       break;
198:     }

200:     if (alg == 0 || alg == 1 || alg == 2 || alg == 3) {
201:       Mat_MPIAIJ *c  = (Mat_MPIAIJ*)(*C)->data;
202:       Mat_APMPI  *ap = c->ap;
203:       PetscOptionsBegin(PetscObjectComm((PetscObject)(*C)),((PetscObject)(*C))->prefix,"MatFreeIntermediateDataStructures","Mat");
204:       ap->freestruct = PETSC_FALSE;
205:       PetscOptionsBool("-mat_freeintermediatedatastructures","Free intermediate data structures", "MatFreeIntermediateDataStructures",ap->freestruct,&ap->freestruct, NULL);
206:       PetscOptionsEnd();
207:     }
208:   }

210:   PetscLogEventBegin(MAT_PtAPNumeric,A,P,0,0);
211:   (*(*C)->ops->ptapnumeric)(A,P,*C);
212:   PetscLogEventEnd(MAT_PtAPNumeric,A,P,0,0);
213:   return(0);
214: }

216: PetscErrorCode MatPtAPNumeric_MPIAIJ_MPIAIJ_scalable(Mat A,Mat P,Mat C)
217: {
218:   PetscErrorCode    ierr;
219:   Mat_MPIAIJ        *a=(Mat_MPIAIJ*)A->data,*p=(Mat_MPIAIJ*)P->data,*c=(Mat_MPIAIJ*)C->data;
220:   Mat_SeqAIJ        *ad=(Mat_SeqAIJ*)(a->A)->data,*ao=(Mat_SeqAIJ*)(a->B)->data;
221:   Mat_SeqAIJ        *ap,*p_loc,*p_oth=NULL,*c_seq;
222:   Mat_APMPI         *ptap = c->ap;
223:   Mat               AP_loc,C_loc,C_oth;
224:   PetscInt          i,rstart,rend,cm,ncols,row,*api,*apj,am = A->rmap->n,apnz,nout;
225:   PetscScalar       *apa;
226:   const PetscInt    *cols;
227:   const PetscScalar *vals;

230:   if (!ptap->AP_loc) {
231:     MPI_Comm comm;
232:     PetscObjectGetComm((PetscObject)C,&comm);
233:     SETERRQ(comm,PETSC_ERR_ARG_WRONGSTATE,"PtAP cannot be reused. Do not call MatFreeIntermediateDataStructures() or use '-mat_freeintermediatedatastructures'");
234:   }

236:   MatZeroEntries(C);

238:   /* 1) get R = Pd^T,Ro = Po^T */
239:   if (ptap->reuse == MAT_REUSE_MATRIX) {
240:     MatTranspose(p->A,MAT_REUSE_MATRIX,&ptap->Rd);
241:     MatTranspose(p->B,MAT_REUSE_MATRIX,&ptap->Ro);
242:   }

244:   /* 2) get AP_loc */
245:   AP_loc = ptap->AP_loc;
246:   ap = (Mat_SeqAIJ*)AP_loc->data;

248:   /* 2-1) get P_oth = ptap->P_oth  and P_loc = ptap->P_loc */
249:   /*-----------------------------------------------------*/
250:   if (ptap->reuse == MAT_REUSE_MATRIX) {
251:     /* P_oth and P_loc are obtained in MatPtASymbolic() when reuse == MAT_INITIAL_MATRIX */
252:     MatGetBrowsOfAoCols_MPIAIJ(A,P,MAT_REUSE_MATRIX,&ptap->startsj_s,&ptap->startsj_r,&ptap->bufa,&ptap->P_oth);
253:     MatMPIAIJGetLocalMat(P,MAT_REUSE_MATRIX,&ptap->P_loc);
254:   }

256:   /* 2-2) compute numeric A_loc*P - dominating part */
257:   /* ---------------------------------------------- */
258:   /* get data from symbolic products */
259:   p_loc = (Mat_SeqAIJ*)(ptap->P_loc)->data;
260:   if (ptap->P_oth) p_oth = (Mat_SeqAIJ*)(ptap->P_oth)->data;

262:   api   = ap->i;
263:   apj   = ap->j;
264:   ISLocalToGlobalMappingApply(ptap->ltog,api[AP_loc->rmap->n],apj,apj);
265:   for (i=0; i<am; i++) {
266:     /* AP[i,:] = A[i,:]*P = Ad*P_loc Ao*P_oth */
267:     apnz = api[i+1] - api[i];
268:     apa = ap->a + api[i];
269:     PetscArrayzero(apa,apnz);
270:     AProw_scalable(i,ad,ao,p_loc,p_oth,api,apj,apa);
271:   }
272:   ISGlobalToLocalMappingApply(ptap->ltog,IS_GTOLM_DROP,api[AP_loc->rmap->n],apj,&nout,apj);
273:   if (api[AP_loc->rmap->n] != nout) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incorrect mapping %D != %D\n",api[AP_loc->rmap->n],nout);

275:   /* 3) C_loc = Rd*AP_loc, C_oth = Ro*AP_loc */
276:   /* Always use scalable version since we are in the MPI scalable version */
277:   MatMatMultNumeric_SeqAIJ_SeqAIJ_Scalable(ptap->Rd,AP_loc,ptap->C_loc);
278:   MatMatMultNumeric_SeqAIJ_SeqAIJ_Scalable(ptap->Ro,AP_loc,ptap->C_oth);

280:   C_loc = ptap->C_loc;
281:   C_oth = ptap->C_oth;

283:   /* add C_loc and Co to to C */
284:   MatGetOwnershipRange(C,&rstart,&rend);

286:   /* C_loc -> C */
287:   cm    = C_loc->rmap->N;
288:   c_seq = (Mat_SeqAIJ*)C_loc->data;
289:   cols = c_seq->j;
290:   vals = c_seq->a;
291:   ISLocalToGlobalMappingApply(ptap->ltog,c_seq->i[C_loc->rmap->n],c_seq->j,c_seq->j);

293:   /* The (fast) MatSetValues_MPIAIJ_CopyFromCSRFormat function can only be used when C->was_assembled is PETSC_FALSE and */
294:   /* when there are no off-processor parts.  */
295:   /* If was_assembled is true, then the statement aj[rowstart_diag+dnz_row] = mat_j[col] - cstart; in MatSetValues_MPIAIJ_CopyFromCSRFormat */
296:   /* is no longer true. Then the more complex function MatSetValues_MPIAIJ() has to be used, where the column index is looked up from */
297:   /* a table, and other, more complex stuff has to be done. */
298:   if (C->assembled) {
299:     C->was_assembled = PETSC_TRUE;
300:     C->assembled     = PETSC_FALSE;
301:   }
302:   if (C->was_assembled) {
303:     for (i=0; i<cm; i++) {
304:       ncols = c_seq->i[i+1] - c_seq->i[i];
305:       row = rstart + i;
306:       MatSetValues_MPIAIJ(C,1,&row,ncols,cols,vals,ADD_VALUES);
307:       cols += ncols; vals += ncols;
308:     }
309:   } else {
310:     MatSetValues_MPIAIJ_CopyFromCSRFormat(C,c_seq->j,c_seq->i,c_seq->a);
311:   }
312:   ISGlobalToLocalMappingApply(ptap->ltog,IS_GTOLM_DROP,c_seq->i[C_loc->rmap->n],c_seq->j,&nout,c_seq->j);
313:   if (c_seq->i[C_loc->rmap->n] != nout) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incorrect mapping %D != %D\n",c_seq->i[C_loc->rmap->n],nout);

315:   /* Co -> C, off-processor part */
316:   cm = C_oth->rmap->N;
317:   c_seq = (Mat_SeqAIJ*)C_oth->data;
318:   cols = c_seq->j;
319:   vals = c_seq->a;
320:   ISLocalToGlobalMappingApply(ptap->ltog,c_seq->i[C_oth->rmap->n],c_seq->j,c_seq->j);
321:   for (i=0; i<cm; i++) {
322:     ncols = c_seq->i[i+1] - c_seq->i[i];
323:     row = p->garray[i];
324:     MatSetValues(C,1,&row,ncols,cols,vals,ADD_VALUES);
325:     cols += ncols; vals += ncols;
326:   }
327:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
328:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);

330:   ptap->reuse = MAT_REUSE_MATRIX;

332:   ISGlobalToLocalMappingApply(ptap->ltog,IS_GTOLM_DROP,c_seq->i[C_oth->rmap->n],c_seq->j,&nout,c_seq->j);
333:   if (c_seq->i[C_oth->rmap->n] != nout) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incorrect mapping %D != %D\n",c_seq->i[C_loc->rmap->n],nout);

335:   /* supporting struct ptap consumes almost same amount of memory as C=PtAP, release it if C will not be updated by A and P */
336:   if (ptap->freestruct) {
337:     MatFreeIntermediateDataStructures(C);
338:   }
339:   return(0);
340: }

342: PetscErrorCode MatPtAPSymbolic_MPIAIJ_MPIAIJ_scalable(Mat A,Mat P,PetscReal fill,Mat *C)
343: {
344:   PetscErrorCode      ierr;
345:   Mat_APMPI           *ptap;
346:   Mat_MPIAIJ          *a=(Mat_MPIAIJ*)A->data,*p=(Mat_MPIAIJ*)P->data,*c;
347:   MPI_Comm            comm;
348:   PetscMPIInt         size,rank;
349:   Mat                 Cmpi,P_loc,P_oth;
350:   PetscFreeSpaceList  free_space=NULL,current_space=NULL;
351:   PetscInt            am=A->rmap->n,pm=P->rmap->n,pN=P->cmap->N,pn=P->cmap->n;
352:   PetscInt            *lnk,i,k,pnz,row,nsend;
353:   PetscMPIInt         tagi,tagj,*len_si,*len_s,*len_ri,icompleted=0,nrecv;
354:   PetscInt            **buf_rj,**buf_ri,**buf_ri_k;
355:   PetscInt            len,proc,*dnz,*onz,*owners,nzi,nspacedouble;
356:   PetscInt            nrows,*buf_s,*buf_si,*buf_si_i,**nextrow,**nextci;
357:   MPI_Request         *swaits,*rwaits;
358:   MPI_Status          *sstatus,rstatus;
359:   PetscLayout         rowmap;
360:   PetscInt            *owners_co,*coi,*coj;    /* i and j array of (p->B)^T*A*P - used in the communication */
361:   PetscMPIInt         *len_r,*id_r;    /* array of length of comm->size, store send/recv matrix values */
362:   PetscInt            *api,*apj,*Jptr,apnz,*prmap=p->garray,con,j,Crmax,*aj,*ai,*pi,nout;
363:   Mat_SeqAIJ          *p_loc,*p_oth=NULL,*ad=(Mat_SeqAIJ*)(a->A)->data,*ao=NULL,*c_loc,*c_oth;
364:   PetscScalar         *apv;
365:   PetscTable          ta;
366:   MatType             mtype;
367:   const char          *prefix;
368: #if defined(PETSC_USE_INFO)
369:   PetscReal           apfill;
370: #endif

373:   PetscObjectGetComm((PetscObject)A,&comm);
374:   MPI_Comm_size(comm,&size);
375:   MPI_Comm_rank(comm,&rank);

377:   if (size > 1) ao = (Mat_SeqAIJ*)(a->B)->data;

379:   /* create symbolic parallel matrix Cmpi */
380:   MatCreate(comm,&Cmpi);
381:   MatGetType(A,&mtype);
382:   MatSetType(Cmpi,mtype);

384:   /* create struct Mat_APMPI and attached it to C later */
385:   PetscNew(&ptap);
386:   ptap->reuse = MAT_INITIAL_MATRIX;
387:   ptap->algType = 0;

389:   /* get P_oth by taking rows of P (= non-zero cols of local A) from other processors */
390:   MatGetBrowsOfAoCols_MPIAIJ(A,P,MAT_INITIAL_MATRIX,&ptap->startsj_s,&ptap->startsj_r,&ptap->bufa,&P_oth);
391:   /* get P_loc by taking all local rows of P */
392:   MatMPIAIJGetLocalMat(P,MAT_INITIAL_MATRIX,&P_loc);

394:   ptap->P_loc = P_loc;
395:   ptap->P_oth = P_oth;

397:   /* (0) compute Rd = Pd^T, Ro = Po^T  */
398:   /* --------------------------------- */
399:   MatTranspose(p->A,MAT_INITIAL_MATRIX,&ptap->Rd);
400:   MatTranspose(p->B,MAT_INITIAL_MATRIX,&ptap->Ro);

402:   /* (1) compute symbolic AP = A_loc*P = Ad*P_loc + Ao*P_oth (api,apj) */
403:   /* ----------------------------------------------------------------- */
404:   p_loc  = (Mat_SeqAIJ*)P_loc->data;
405:   if (P_oth) p_oth = (Mat_SeqAIJ*)P_oth->data;

407:   /* create and initialize a linked list */
408:   PetscTableCreate(pn,pN,&ta); /* for compute AP_loc and Cmpi */
409:   MatRowMergeMax_SeqAIJ(p_loc,P_loc->rmap->N,ta);
410:   MatRowMergeMax_SeqAIJ(p_oth,P_oth->rmap->N,ta);
411:   PetscTableGetCount(ta,&Crmax); /* Crmax = nnz(sum of Prows) */

413:   PetscLLCondensedCreate_Scalable(Crmax,&lnk);

415:   /* Initial FreeSpace size is fill*(nnz(A) + nnz(P)) */
416:   if (ao) {
417:     PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,PetscIntSumTruncate(ad->i[am],PetscIntSumTruncate(ao->i[am],p_loc->i[pm]))),&free_space);
418:   } else {
419:     PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,PetscIntSumTruncate(ad->i[am],p_loc->i[pm])),&free_space);
420:   }
421:   current_space = free_space;
422:   nspacedouble  = 0;

424:   PetscMalloc1(am+1,&api);
425:   api[0] = 0;
426:   for (i=0; i<am; i++) {
427:     /* diagonal portion: Ad[i,:]*P */
428:     ai = ad->i; pi = p_loc->i;
429:     nzi = ai[i+1] - ai[i];
430:     aj  = ad->j + ai[i];
431:     for (j=0; j<nzi; j++) {
432:       row  = aj[j];
433:       pnz  = pi[row+1] - pi[row];
434:       Jptr = p_loc->j + pi[row];
435:       /* add non-zero cols of P into the sorted linked list lnk */
436:       PetscLLCondensedAddSorted_Scalable(pnz,Jptr,lnk);
437:     }
438:     /* off-diagonal portion: Ao[i,:]*P */
439:     if (ao) {
440:       ai = ao->i; pi = p_oth->i;
441:       nzi = ai[i+1] - ai[i];
442:       aj  = ao->j + ai[i];
443:       for (j=0; j<nzi; j++) {
444:         row  = aj[j];
445:         pnz  = pi[row+1] - pi[row];
446:         Jptr = p_oth->j + pi[row];
447:         PetscLLCondensedAddSorted_Scalable(pnz,Jptr,lnk);
448:       }
449:     }
450:     apnz     = lnk[0];
451:     api[i+1] = api[i] + apnz;

453:     /* if free space is not available, double the total space in the list */
454:     if (current_space->local_remaining<apnz) {
455:       PetscFreeSpaceGet(PetscIntSumTruncate(apnz,current_space->total_array_size),&current_space);
456:       nspacedouble++;
457:     }

459:     /* Copy data into free space, then initialize lnk */
460:     PetscLLCondensedClean_Scalable(apnz,current_space->array,lnk);

462:     current_space->array           += apnz;
463:     current_space->local_used      += apnz;
464:     current_space->local_remaining -= apnz;
465:   }
466:   /* Allocate space for apj and apv, initialize apj, and */
467:   /* destroy list of free space and other temporary array(s) */
468:   PetscCalloc2(api[am],&apj,api[am],&apv);
469:   PetscFreeSpaceContiguous(&free_space,apj);
470:   PetscLLCondensedDestroy_Scalable(lnk);

472:   /* Create AP_loc for reuse */
473:   MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,am,pN,api,apj,apv,&ptap->AP_loc);
474:   MatSeqAIJCompactOutExtraColumns_SeqAIJ(ptap->AP_loc, &ptap->ltog);

476: #if defined(PETSC_USE_INFO)
477:   if (ao) {
478:     apfill = (PetscReal)api[am]/(ad->i[am]+ao->i[am]+p_loc->i[pm]+1);
479:   } else {
480:     apfill = (PetscReal)api[am]/(ad->i[am]+p_loc->i[pm]+1);
481:   }
482:   ptap->AP_loc->info.mallocs           = nspacedouble;
483:   ptap->AP_loc->info.fill_ratio_given  = fill;
484:   ptap->AP_loc->info.fill_ratio_needed = apfill;

486:   if (api[am]) {
487:     PetscInfo3(ptap->AP_loc,"Scalable algorithm, AP_loc reallocs %D; Fill ratio: given %g needed %g.\n",nspacedouble,(double)fill,(double)apfill);
488:     PetscInfo1(ptap->AP_loc,"Use MatPtAP(A,B,MatReuse,%g,&C) for best AP_loc performance.;\n",(double)apfill);
489:   } else {
490:     PetscInfo(ptap->AP_loc,"Scalable algorithm, AP_loc is empty \n");
491:   }
492: #endif

494:   /* (2-1) compute symbolic Co = Ro*AP_loc  */
495:   /* ------------------------------------ */
496:   MatGetOptionsPrefix(A,&prefix);
497:   MatSetOptionsPrefix(ptap->Ro,prefix);
498:   MatAppendOptionsPrefix(ptap->Ro,"inner_offdiag_");
499:   MatMatMultSymbolic_SeqAIJ_SeqAIJ(ptap->Ro,ptap->AP_loc,fill,&ptap->C_oth);

501:   /* (3) send coj of C_oth to other processors  */
502:   /* ------------------------------------------ */
503:   /* determine row ownership */
504:   PetscLayoutCreate(comm,&rowmap);
505:   rowmap->n  = pn;
506:   rowmap->bs = 1;
507:   PetscLayoutSetUp(rowmap);
508:   owners = rowmap->range;

510:   /* determine the number of messages to send, their lengths */
511:   PetscMalloc4(size,&len_s,size,&len_si,size,&sstatus,size+2,&owners_co);
512:   PetscArrayzero(len_s,size);
513:   PetscArrayzero(len_si,size);

515:   c_oth = (Mat_SeqAIJ*)ptap->C_oth->data;
516:   coi   = c_oth->i; coj = c_oth->j;
517:   con   = ptap->C_oth->rmap->n;
518:   proc  = 0;
519:   ISLocalToGlobalMappingApply(ptap->ltog,coi[con],coj,coj);
520:   for (i=0; i<con; i++) {
521:     while (prmap[i] >= owners[proc+1]) proc++;
522:     len_si[proc]++;               /* num of rows in Co(=Pt*AP) to be sent to [proc] */
523:     len_s[proc] += coi[i+1] - coi[i]; /* num of nonzeros in Co to be sent to [proc] */
524:   }

526:   len          = 0; /* max length of buf_si[], see (4) */
527:   owners_co[0] = 0;
528:   nsend        = 0;
529:   for (proc=0; proc<size; proc++) {
530:     owners_co[proc+1] = owners_co[proc] + len_si[proc];
531:     if (len_s[proc]) {
532:       nsend++;
533:       len_si[proc] = 2*(len_si[proc] + 1); /* length of buf_si to be sent to [proc] */
534:       len         += len_si[proc];
535:     }
536:   }

538:   /* determine the number and length of messages to receive for coi and coj  */
539:   PetscGatherNumberOfMessages(comm,NULL,len_s,&nrecv);
540:   PetscGatherMessageLengths2(comm,nsend,nrecv,len_s,len_si,&id_r,&len_r,&len_ri);

542:   /* post the Irecv and Isend of coj */
543:   PetscCommGetNewTag(comm,&tagj);
544:   PetscPostIrecvInt(comm,tagj,nrecv,id_r,len_r,&buf_rj,&rwaits);
545:   PetscMalloc1(nsend+1,&swaits);
546:   for (proc=0, k=0; proc<size; proc++) {
547:     if (!len_s[proc]) continue;
548:     i    = owners_co[proc];
549:     MPI_Isend(coj+coi[i],len_s[proc],MPIU_INT,proc,tagj,comm,swaits+k);
550:     k++;
551:   }

553:   /* (2-2) compute symbolic C_loc = Rd*AP_loc */
554:   /* ---------------------------------------- */
555:   MatSetOptionsPrefix(ptap->Rd,prefix);
556:   MatAppendOptionsPrefix(ptap->Rd,"inner_diag_");
557:   MatMatMultSymbolic_SeqAIJ_SeqAIJ(ptap->Rd,ptap->AP_loc,fill,&ptap->C_loc);
558:   c_loc = (Mat_SeqAIJ*)ptap->C_loc->data;
559:   ISLocalToGlobalMappingApply(ptap->ltog,c_loc->i[ptap->C_loc->rmap->n],c_loc->j,c_loc->j);

561:   /* receives coj are complete */
562:   for (i=0; i<nrecv; i++) {
563:     MPI_Waitany(nrecv,rwaits,&icompleted,&rstatus);
564:   }
565:   PetscFree(rwaits);
566:   if (nsend) {MPI_Waitall(nsend,swaits,sstatus);}

568:   /* add received column indices into ta to update Crmax */
569:   for (k=0; k<nrecv; k++) {/* k-th received message */
570:     Jptr = buf_rj[k];
571:     for (j=0; j<len_r[k]; j++) {
572:       PetscTableAdd(ta,*(Jptr+j)+1,1,INSERT_VALUES);
573:     }
574:   }
575:   PetscTableGetCount(ta,&Crmax);
576:   PetscTableDestroy(&ta);

578:   /* (4) send and recv coi */
579:   /*-----------------------*/
580:   PetscCommGetNewTag(comm,&tagi);
581:   PetscPostIrecvInt(comm,tagi,nrecv,id_r,len_ri,&buf_ri,&rwaits);
582:   PetscMalloc1(len+1,&buf_s);
583:   buf_si = buf_s;  /* points to the beginning of k-th msg to be sent */
584:   for (proc=0,k=0; proc<size; proc++) {
585:     if (!len_s[proc]) continue;
586:     /* form outgoing message for i-structure:
587:          buf_si[0]:                 nrows to be sent
588:                [1:nrows]:           row index (global)
589:                [nrows+1:2*nrows+1]: i-structure index
590:     */
591:     /*-------------------------------------------*/
592:     nrows       = len_si[proc]/2 - 1; /* num of rows in Co to be sent to [proc] */
593:     buf_si_i    = buf_si + nrows+1;
594:     buf_si[0]   = nrows;
595:     buf_si_i[0] = 0;
596:     nrows       = 0;
597:     for (i=owners_co[proc]; i<owners_co[proc+1]; i++) {
598:       nzi = coi[i+1] - coi[i];
599:       buf_si_i[nrows+1] = buf_si_i[nrows] + nzi;  /* i-structure */
600:       buf_si[nrows+1]   = prmap[i] -owners[proc]; /* local row index */
601:       nrows++;
602:     }
603:     MPI_Isend(buf_si,len_si[proc],MPIU_INT,proc,tagi,comm,swaits+k);
604:     k++;
605:     buf_si += len_si[proc];
606:   }
607:   for (i=0; i<nrecv; i++) {
608:     MPI_Waitany(nrecv,rwaits,&icompleted,&rstatus);
609:   }
610:   PetscFree(rwaits);
611:   if (nsend) {MPI_Waitall(nsend,swaits,sstatus);}

613:   PetscFree4(len_s,len_si,sstatus,owners_co);
614:   PetscFree(len_ri);
615:   PetscFree(swaits);
616:   PetscFree(buf_s);

618:   /* (5) compute the local portion of Cmpi      */
619:   /* ------------------------------------------ */
620:   /* set initial free space to be Crmax, sufficient for holding nozeros in each row of Cmpi */
621:   PetscFreeSpaceGet(Crmax,&free_space);
622:   current_space = free_space;

624:   PetscMalloc3(nrecv,&buf_ri_k,nrecv,&nextrow,nrecv,&nextci);
625:   for (k=0; k<nrecv; k++) {
626:     buf_ri_k[k] = buf_ri[k]; /* beginning of k-th recved i-structure */
627:     nrows       = *buf_ri_k[k];
628:     nextrow[k]  = buf_ri_k[k] + 1;  /* next row number of k-th recved i-structure */
629:     nextci[k]   = buf_ri_k[k] + (nrows + 1); /* poins to the next i-structure of k-th recved i-structure  */
630:   }

632:   MatPreallocateInitialize(comm,pn,pn,dnz,onz);
633:   PetscLLCondensedCreate_Scalable(Crmax,&lnk);
634:   for (i=0; i<pn; i++) {
635:     /* add C_loc into Cmpi */
636:     nzi  = c_loc->i[i+1] - c_loc->i[i];
637:     Jptr = c_loc->j + c_loc->i[i];
638:     PetscLLCondensedAddSorted_Scalable(nzi,Jptr,lnk);

640:     /* add received col data into lnk */
641:     for (k=0; k<nrecv; k++) { /* k-th received message */
642:       if (i == *nextrow[k]) { /* i-th row */
643:         nzi  = *(nextci[k]+1) - *nextci[k];
644:         Jptr = buf_rj[k] + *nextci[k];
645:         PetscLLCondensedAddSorted_Scalable(nzi,Jptr,lnk);
646:         nextrow[k]++; nextci[k]++;
647:       }
648:     }
649:     nzi = lnk[0];

651:     /* copy data into free space, then initialize lnk */
652:     PetscLLCondensedClean_Scalable(nzi,current_space->array,lnk);
653:     MatPreallocateSet(i+owners[rank],nzi,current_space->array,dnz,onz);
654:   }
655:   PetscFree3(buf_ri_k,nextrow,nextci);
656:   PetscLLCondensedDestroy_Scalable(lnk);
657:   PetscFreeSpaceDestroy(free_space);

659:   /* local sizes and preallocation */
660:   MatSetSizes(Cmpi,pn,pn,PETSC_DETERMINE,PETSC_DETERMINE);
661:   if (P->cmap->bs > 0) {
662:     PetscLayoutSetBlockSize(Cmpi->rmap,P->cmap->bs);
663:     PetscLayoutSetBlockSize(Cmpi->cmap,P->cmap->bs);
664:   }
665:   MatMPIAIJSetPreallocation(Cmpi,0,dnz,0,onz);
666:   MatPreallocateFinalize(dnz,onz);

668:   /* members in merge */
669:   PetscFree(id_r);
670:   PetscFree(len_r);
671:   PetscFree(buf_ri[0]);
672:   PetscFree(buf_ri);
673:   PetscFree(buf_rj[0]);
674:   PetscFree(buf_rj);
675:   PetscLayoutDestroy(&rowmap);

677:   /* attach the supporting struct to Cmpi for reuse */
678:   c = (Mat_MPIAIJ*)Cmpi->data;
679:   c->ap           = ptap;
680:   ptap->duplicate = Cmpi->ops->duplicate;
681:   ptap->destroy   = Cmpi->ops->destroy;
682:   ptap->view      = Cmpi->ops->view;

684:   /* Cmpi is not ready for use - assembly will be done by MatPtAPNumeric() */
685:   Cmpi->assembled        = PETSC_FALSE;
686:   Cmpi->ops->ptapnumeric = MatPtAPNumeric_MPIAIJ_MPIAIJ_scalable;
687:   Cmpi->ops->destroy     = MatDestroy_MPIAIJ_PtAP;
688:   Cmpi->ops->view        = MatView_MPIAIJ_PtAP;
689:   Cmpi->ops->freeintermediatedatastructures = MatFreeIntermediateDataStructures_MPIAIJ_AP;
690:   *C                     = Cmpi;

692:    nout = 0;
693:    ISGlobalToLocalMappingApply(ptap->ltog,IS_GTOLM_DROP,c_oth->i[ptap->C_oth->rmap->n],c_oth->j,&nout,c_oth->j);
694:    if (c_oth->i[ptap->C_oth->rmap->n] != nout) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incorrect mapping %D != %D\n",c_oth->i[ptap->C_oth->rmap->n],nout);
695:    ISGlobalToLocalMappingApply(ptap->ltog,IS_GTOLM_DROP,c_loc->i[ptap->C_loc->rmap->n],c_loc->j,&nout,c_loc->j);
696:    if (c_loc->i[ptap->C_loc->rmap->n] != nout) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incorrect mapping %D != %D\n",c_loc->i[ptap->C_loc->rmap->n],nout);

698:   return(0);
699: }

701: PETSC_STATIC_INLINE PetscErrorCode MatPtAPSymbolicComputeOneRowOfAP_private(Mat A,Mat P,Mat P_oth,PetscInt i,PetscHSetI dht,PetscHSetI oht)
702: {
703:   Mat_MPIAIJ           *a=(Mat_MPIAIJ*)A->data,*p=(Mat_MPIAIJ*)P->data;
704:   Mat_SeqAIJ           *ad=(Mat_SeqAIJ*)(a->A)->data,*ao=(Mat_SeqAIJ*)(a->B)->data,*p_oth=(Mat_SeqAIJ*)P_oth->data,*pd=(Mat_SeqAIJ*)p->A->data,*po=(Mat_SeqAIJ*)p->B->data;
705:   PetscInt             *ai,nzi,j,*aj,row,col,*pi,*pj,pnz,nzpi,*p_othcols,k;
706:   PetscInt             pcstart,pcend,column;
707:   PetscErrorCode       ierr;


711:   pcstart = P->cmap->rstart;
712:   pcend   = P->cmap->rend;
713:   /* diagonal portion: Ad[i,:]*P */
714:   ai = ad->i;
715:   nzi = ai[i+1] - ai[i];
716:   aj  = ad->j + ai[i];
717:   for (j=0; j<nzi; j++) {
718:     row  = aj[j];
719:     nzpi = pd->i[row+1] - pd->i[row];
720:     pj  = pd->j + pd->i[row];
721:     for (k=0; k<nzpi; k++) {
722:       PetscHSetIAdd(dht,pj[k]+pcstart);
723:     }
724:   }
725:   for (j=0; j<nzi; j++) {
726:     row  = aj[j];
727:     nzpi = po->i[row+1] - po->i[row];
728:     pj  = po->j + po->i[row];
729:     for (k=0; k<nzpi; k++) {
730:       PetscHSetIAdd(oht,p->garray[pj[k]]);
731:     }
732:   }

734:   /* off diagonal part: Ao[i, :]*P_oth */
735:   if (ao) {
736:     ai = ao->i;
737:     pi = p_oth->i;
738:     nzi = ai[i+1] - ai[i];
739:     aj  = ao->j + ai[i];
740:     for (j=0; j<nzi; j++) {
741:       row  = aj[j];
742:       pnz  = pi[row+1] - pi[row];
743:       p_othcols = p_oth->j + pi[row];
744:       for (col=0; col<pnz; col++) {
745:         column = p_othcols[col];
746:         if (column>=pcstart && column<pcend) {
747:           PetscHSetIAdd(dht,column);
748:         } else {
749:           PetscHSetIAdd(oht,column);
750:         }
751:       }
752:     }
753:   } /* end if (ao) */
754:   return(0);
755: }

757: PETSC_STATIC_INLINE PetscErrorCode MatPtAPNumericComputeOneRowOfAP_private(Mat A,Mat P,Mat P_oth,PetscInt i,PetscHMapIV hmap)
758: {
759:   Mat_MPIAIJ           *a=(Mat_MPIAIJ*)A->data,*p=(Mat_MPIAIJ*)P->data;
760:   Mat_SeqAIJ           *ad=(Mat_SeqAIJ*)(a->A)->data,*ao=(Mat_SeqAIJ*)(a->B)->data,*p_oth=(Mat_SeqAIJ*)P_oth->data,*pd=(Mat_SeqAIJ*)p->A->data,*po=(Mat_SeqAIJ*)p->B->data;
761:   PetscInt             *ai,nzi,j,*aj,row,col,*pi,pnz,*p_othcols,pcstart,*pj,k,nzpi;
762:   PetscScalar          ra,*aa,*pa;
763:   PetscErrorCode       ierr;

766:   pcstart = P->cmap->rstart;
767:   /* diagonal portion: Ad[i,:]*P */
768:   ai  = ad->i;
769:   nzi = ai[i+1] - ai[i];
770:   aj  = ad->j + ai[i];
771:   aa  = ad->a + ai[i];

773:   for (j=0; j<nzi; j++) {
774:     ra   = aa[j];
775:     row  = aj[j];
776:     nzpi = pd->i[row+1] - pd->i[row];
777:     pj = pd->j + pd->i[row];
778:     pa = pd->a + pd->i[row];
779:     for (k=0; k<nzpi; k++) {
780:       PetscHMapIVAddValue(hmap,pj[k]+pcstart,ra*pa[k]);
781:     }
782:     PetscLogFlops(2.0*nzpi);
783:   }
784:   for (j=0; j<nzi; j++) {
785:     ra   = aa[j];
786:     row  = aj[j];
787:     nzpi = po->i[row+1] - po->i[row];
788:     pj = po->j + po->i[row];
789:     pa = po->a + po->i[row];
790:     for (k=0; k<nzpi; k++) {
791:       PetscHMapIVAddValue(hmap,p->garray[pj[k]],ra*pa[k]);
792:     }
793:     PetscLogFlops(2.0*nzpi);
794:   }


797:   /* off diagonal part: Ao[i, :]*P_oth */
798:   if (ao) {
799:     ai = ao->i;
800:     pi = p_oth->i;
801:     nzi = ai[i+1] - ai[i];
802:     aj  = ao->j + ai[i];
803:     aa  = ao->a + ai[i];
804:     for (j=0; j<nzi; j++) {
805:       row  = aj[j];
806:       ra   = aa[j];
807:       pnz  = pi[row+1] - pi[row];
808:       p_othcols = p_oth->j + pi[row];
809:       pa   = p_oth->a + pi[row];
810:       for (col=0; col<pnz; col++) {
811:         PetscHMapIVAddValue(hmap,p_othcols[col],ra*pa[col]);
812:       }
813:       PetscLogFlops(2.0*pnz);
814:     }
815:   } /* end if (ao) */

817:   return(0);
818: }

820: PetscErrorCode MatPtAPNumeric_MPIAIJ_MPIAIJ_allatonce(Mat A,Mat P,Mat C)
821: {
822:   PetscErrorCode    ierr;
823:   Mat_MPIAIJ        *p=(Mat_MPIAIJ*)P->data,*c=(Mat_MPIAIJ*)C->data;
824:   Mat_SeqAIJ        *cd,*co,*po,*pd;
825:   Mat_APMPI         *ptap = c->ap;
826:   PetscHMapIV       hmap;
827:   PetscInt          i,j,jj,kk,nzi,*c_rmtj,voff,*c_othj,pn,pon,pcstart,pcend,ccstart,ccend,row,am,*poj,*pdj,*apindices,cmaxr,*c_rmtc,*c_rmtjj,*dcc,*occ,loc;
828:   PetscScalar       *c_rmta,*c_otha,*poa,*pda,*apvalues,*apvaluestmp,*c_rmtaa;
829:   MPI_Comm          comm;

832:   PetscObjectGetComm((PetscObject)A,&comm);
833:   if (!ptap->P_oth) SETERRQ(comm,PETSC_ERR_ARG_WRONGSTATE,"PtAP cannot be reused. Do not call MatFreeIntermediateDataStructures() or use '-mat_freeintermediatedatastructures'");

835:   MatZeroEntries(C);

837:   /* Get P_oth = ptap->P_oth  and P_loc = ptap->P_loc */
838:   /*-----------------------------------------------------*/
839:   if (ptap->reuse == MAT_REUSE_MATRIX) {
840:     /* P_oth and P_loc are obtained in MatPtASymbolic() when reuse == MAT_INITIAL_MATRIX */
841:     MatGetBrowsOfAoCols_MPIAIJ(A,P,MAT_REUSE_MATRIX,&ptap->startsj_s,&ptap->startsj_r,&ptap->bufa,&ptap->P_oth);
842:   }

844:   po = (Mat_SeqAIJ*) p->B->data;
845:   pd = (Mat_SeqAIJ*) p->A->data;

847:   MatGetLocalSize(p->B,NULL,&pon);

849:   PetscCalloc2(ptap->c_rmti[pon],&c_rmtj,ptap->c_rmti[pon],&c_rmta);
850:   MatGetLocalSize(A,&am,NULL);
851:   cmaxr = 0;
852:   for (i=0; i<pon; i++) {
853:     cmaxr = PetscMax(cmaxr,ptap->c_rmti[i+1]-ptap->c_rmti[i]);
854:   }
855:   PetscCalloc4(cmaxr,&apindices,cmaxr,&apvalues,cmaxr,&apvaluestmp,pon,&c_rmtc);
856:   PetscHMapIVCreate(&hmap);
857:   PetscHMapIVResize(hmap,cmaxr);
858:   for (i=0; i<am && pon; i++) {
859:     PetscHMapIVClear(hmap);
860:     nzi = po->i[i+1] - po->i[i];
861:     if (!nzi) continue;
862:     MatPtAPNumericComputeOneRowOfAP_private(A,P,ptap->P_oth,i,hmap);
863:     voff = 0;
864:     PetscHMapIVGetPairs(hmap,&voff,apindices,apvalues);
865:     if (!voff) continue;

867:     /* Form C(ii, :) */
868:     poj = po->j + po->i[i];
869:     poa = po->a + po->i[i];
870:     for (j=0; j<nzi; j++) {
871:       c_rmtjj = c_rmtj + ptap->c_rmti[poj[j]];
872:       c_rmtaa = c_rmta + ptap->c_rmti[poj[j]];
873:       for (jj=0; jj<voff; jj++) {
874:         apvaluestmp[jj] = apvalues[jj]*poa[j];
875:         /*If the row is empty */
876:         if (!c_rmtc[poj[j]]) {
877:           c_rmtjj[jj] = apindices[jj];
878:           c_rmtaa[jj] = apvaluestmp[jj];
879:           c_rmtc[poj[j]]++;
880:         } else {
881:           PetscFindInt(apindices[jj],c_rmtc[poj[j]],c_rmtjj,&loc);
882:           if (loc>=0){ /* hit */
883:             c_rmtaa[loc] += apvaluestmp[jj];
884:             PetscLogFlops(1.0);
885:           } else { /* new element */
886:             loc = -(loc+1);
887:             /* Move data backward */
888:             for (kk=c_rmtc[poj[j]]; kk>loc; kk--) {
889:               c_rmtjj[kk] = c_rmtjj[kk-1];
890:               c_rmtaa[kk] = c_rmtaa[kk-1];
891:             }/* End kk */
892:             c_rmtjj[loc] = apindices[jj];
893:             c_rmtaa[loc] = apvaluestmp[jj];
894:             c_rmtc[poj[j]]++;
895:           }
896:         }
897:         PetscLogFlops(voff);
898:       } /* End jj */
899:     } /* End j */
900:   } /* End i */

902:   PetscFree4(apindices,apvalues,apvaluestmp,c_rmtc);

904:   MatGetLocalSize(P,NULL,&pn);
905:   PetscCalloc2(ptap->c_othi[pn],&c_othj,ptap->c_othi[pn],&c_otha);

907:   PetscSFReduceBegin(ptap->sf,MPIU_INT,c_rmtj,c_othj,MPIU_REPLACE);
908:   PetscSFReduceBegin(ptap->sf,MPIU_SCALAR,c_rmta,c_otha,MPIU_REPLACE);
909:   MatGetOwnershipRangeColumn(P,&pcstart,&pcend);
910:   cd = (Mat_SeqAIJ*)(c->A)->data;
911:   co = (Mat_SeqAIJ*)(c->B)->data;

913:   cmaxr = 0;
914:   for (i=0; i<pn; i++) {
915:     cmaxr = PetscMax(cmaxr,(cd->i[i+1]-cd->i[i])+(co->i[i+1]-co->i[i]));
916:   }
917:   PetscCalloc5(cmaxr,&apindices,cmaxr,&apvalues,cmaxr,&apvaluestmp,pn,&dcc,pn,&occ);
918:   PetscHMapIVCreate(&hmap);
919:   PetscHMapIVResize(hmap,cmaxr);
920:   for (i=0; i<am && pn; i++) {
921:     PetscHMapIVClear(hmap);
922:     nzi = pd->i[i+1] - pd->i[i];
923:     if (!nzi) continue;
924:     MatPtAPNumericComputeOneRowOfAP_private(A,P,ptap->P_oth,i,hmap);
925:     voff = 0;
926:     PetscHMapIVGetPairs(hmap,&voff,apindices,apvalues);
927:     if (!voff) continue;
928:     /* Form C(ii, :) */
929:     pdj = pd->j + pd->i[i];
930:     pda = pd->a + pd->i[i];
931:     for (j=0; j<nzi; j++) {
932:       row = pcstart + pdj[j];
933:       for (jj=0; jj<voff; jj++) {
934:         apvaluestmp[jj] = apvalues[jj]*pda[j];
935:       }
936:       PetscLogFlops(voff);
937:       MatSetValues(C,1,&row,voff,apindices,apvaluestmp,ADD_VALUES);
938:     }
939:   }

941:   MatGetOwnershipRangeColumn(C,&ccstart,&ccend);
942:   PetscFree5(apindices,apvalues,apvaluestmp,dcc,occ);
943:   PetscHMapIVDestroy(&hmap);
944:   PetscSFReduceEnd(ptap->sf,MPIU_INT,c_rmtj,c_othj,MPI_SUM);
945:   PetscSFReduceEnd(ptap->sf,MPIU_SCALAR,c_rmta,c_otha,MPIU_REPLACE);
946:   PetscFree2(c_rmtj,c_rmta);

948:   /* Add contributions from remote */
949:   for (i = 0; i < pn; i++) {
950:     row = i + pcstart;
951:     MatSetValues(C,1,&row,ptap->c_othi[i+1]-ptap->c_othi[i],c_othj+ptap->c_othi[i],c_otha+ptap->c_othi[i],ADD_VALUES);
952:   }
953:   PetscFree2(c_othj,c_otha);

955:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
956:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);

958:   ptap->reuse = MAT_REUSE_MATRIX;

960:   /* supporting struct ptap consumes almost same amount of memory as C=PtAP, release it if C will not be updated by A and P */
961:   if (ptap->freestruct) {
962:     MatFreeIntermediateDataStructures(C);
963:   }
964:   return(0);
965: }

967: PetscErrorCode MatPtAPNumeric_MPIAIJ_MPIAIJ_allatonce_merged(Mat A,Mat P,Mat C)
968: {
969:   PetscErrorCode    ierr;
970:   Mat_MPIAIJ        *p=(Mat_MPIAIJ*)P->data,*c=(Mat_MPIAIJ*)C->data;
971:   Mat_SeqAIJ        *cd,*co,*po,*pd;
972:   Mat_APMPI         *ptap = c->ap;
973:   PetscHMapIV       hmap;
974:   PetscInt          i,j,jj,kk,nzi,dnzi,*c_rmtj,voff,*c_othj,pn,pon,pcstart,pcend,row,am,*poj,*pdj,*apindices,cmaxr,*c_rmtc,*c_rmtjj,loc;
975:   PetscScalar       *c_rmta,*c_otha,*poa,*pda,*apvalues,*apvaluestmp,*c_rmtaa;
976:   MPI_Comm          comm;

979:   PetscObjectGetComm((PetscObject)A,&comm);
980:   if (!ptap->P_oth) SETERRQ(comm,PETSC_ERR_ARG_WRONGSTATE,"PtAP cannot be reused. Do not call MatFreeIntermediateDataStructures() or use '-mat_freeintermediatedatastructures'");

982:   MatZeroEntries(C);

984:   /* Get P_oth = ptap->P_oth  and P_loc = ptap->P_loc */
985:   /*-----------------------------------------------------*/
986:   if (ptap->reuse == MAT_REUSE_MATRIX) {
987:     /* P_oth and P_loc are obtained in MatPtASymbolic() when reuse == MAT_INITIAL_MATRIX */
988:     MatGetBrowsOfAoCols_MPIAIJ(A,P,MAT_REUSE_MATRIX,&ptap->startsj_s,&ptap->startsj_r,&ptap->bufa,&ptap->P_oth);
989:   }

991:   po = (Mat_SeqAIJ*) p->B->data;
992:   pd = (Mat_SeqAIJ*) p->A->data;

994:   MatGetLocalSize(p->B,NULL,&pon);
995:   MatGetLocalSize(P,NULL,&pn);

997:   PetscCalloc2(ptap->c_rmti[pon],&c_rmtj,ptap->c_rmti[pon],&c_rmta);
998:   MatGetLocalSize(A,&am,NULL);
999:   MatGetOwnershipRangeColumn(P,&pcstart,&pcend);
1000:   cmaxr = 0;
1001:   for (i=0; i<pon; i++) {
1002:     cmaxr = PetscMax(cmaxr,ptap->c_rmti[i+1]-ptap->c_rmti[i]);
1003:   }
1004:   cd = (Mat_SeqAIJ*)(c->A)->data;
1005:   co = (Mat_SeqAIJ*)(c->B)->data;
1006:   for (i=0; i<pn; i++) {
1007:     cmaxr = PetscMax(cmaxr,(cd->i[i+1]-cd->i[i])+(co->i[i+1]-co->i[i]));
1008:   }
1009:   PetscCalloc4(cmaxr,&apindices,cmaxr,&apvalues,cmaxr,&apvaluestmp,pon,&c_rmtc);
1010:   PetscHMapIVCreate(&hmap);
1011:   PetscHMapIVResize(hmap,cmaxr);
1012:   for (i=0; i<am && (pon || pn); i++) {
1013:     PetscHMapIVClear(hmap);
1014:     nzi  = po->i[i+1] - po->i[i];
1015:     dnzi = pd->i[i+1] - pd->i[i];
1016:     if (!nzi && !dnzi) continue;
1017:     MatPtAPNumericComputeOneRowOfAP_private(A,P,ptap->P_oth,i,hmap);
1018:     voff = 0;
1019:     PetscHMapIVGetPairs(hmap,&voff,apindices,apvalues);
1020:     if (!voff) continue;

1022:     /* Form remote C(ii, :) */
1023:     poj = po->j + po->i[i];
1024:     poa = po->a + po->i[i];
1025:     for (j=0; j<nzi; j++) {
1026:       c_rmtjj = c_rmtj + ptap->c_rmti[poj[j]];
1027:       c_rmtaa = c_rmta + ptap->c_rmti[poj[j]];
1028:       for (jj=0; jj<voff; jj++) {
1029:         apvaluestmp[jj] = apvalues[jj]*poa[j];
1030:         /*If the row is empty */
1031:         if (!c_rmtc[poj[j]]) {
1032:           c_rmtjj[jj] = apindices[jj];
1033:           c_rmtaa[jj] = apvaluestmp[jj];
1034:           c_rmtc[poj[j]]++;
1035:         } else {
1036:           PetscFindInt(apindices[jj],c_rmtc[poj[j]],c_rmtjj,&loc);
1037:           if (loc>=0){ /* hit */
1038:             c_rmtaa[loc] += apvaluestmp[jj];
1039:             PetscLogFlops(1.0);
1040:           } else { /* new element */
1041:             loc = -(loc+1);
1042:             /* Move data backward */
1043:             for (kk=c_rmtc[poj[j]]; kk>loc; kk--) {
1044:               c_rmtjj[kk] = c_rmtjj[kk-1];
1045:               c_rmtaa[kk] = c_rmtaa[kk-1];
1046:             }/* End kk */
1047:             c_rmtjj[loc] = apindices[jj];
1048:             c_rmtaa[loc] = apvaluestmp[jj];
1049:             c_rmtc[poj[j]]++;
1050:           }
1051:         }
1052:       } /* End jj */
1053:       PetscLogFlops(voff);
1054:     } /* End j */

1056:     /* Form local C(ii, :) */
1057:     pdj = pd->j + pd->i[i];
1058:     pda = pd->a + pd->i[i];
1059:     for (j=0; j<dnzi; j++) {
1060:       row = pcstart + pdj[j];
1061:       for (jj=0; jj<voff; jj++) {
1062:         apvaluestmp[jj] = apvalues[jj]*pda[j];
1063:       }/* End kk */
1064:       PetscLogFlops(voff);
1065:       MatSetValues(C,1,&row,voff,apindices,apvaluestmp,ADD_VALUES);
1066:     }/* End j */
1067:   } /* End i */

1069:   PetscFree4(apindices,apvalues,apvaluestmp,c_rmtc);
1070:   PetscHMapIVDestroy(&hmap);
1071:   PetscCalloc2(ptap->c_othi[pn],&c_othj,ptap->c_othi[pn],&c_otha);

1073:   PetscSFReduceBegin(ptap->sf,MPIU_INT,c_rmtj,c_othj,MPIU_REPLACE);
1074:   PetscSFReduceBegin(ptap->sf,MPIU_SCALAR,c_rmta,c_otha,MPIU_REPLACE);
1075:   PetscSFReduceEnd(ptap->sf,MPIU_INT,c_rmtj,c_othj,MPI_SUM);
1076:   PetscSFReduceEnd(ptap->sf,MPIU_SCALAR,c_rmta,c_otha,MPIU_REPLACE);
1077:   PetscFree2(c_rmtj,c_rmta);

1079:   /* Add contributions from remote */
1080:   for (i = 0; i < pn; i++) {
1081:     row = i + pcstart;
1082:     MatSetValues(C,1,&row,ptap->c_othi[i+1]-ptap->c_othi[i],c_othj+ptap->c_othi[i],c_otha+ptap->c_othi[i],ADD_VALUES);
1083:   }
1084:   PetscFree2(c_othj,c_otha);

1086:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
1087:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);

1089:   ptap->reuse = MAT_REUSE_MATRIX;

1091:   /* supporting struct ptap consumes almost same amount of memory as C=PtAP, release it if C will not be updated by A and P */
1092:   if (ptap->freestruct) {
1093:     MatFreeIntermediateDataStructures(C);
1094:   }
1095:   return(0);
1096: }

1098: PetscErrorCode MatPtAPSymbolic_MPIAIJ_MPIAIJ_allatonce(Mat A,Mat P,PetscReal fill,Mat *C)
1099: {
1100:   Mat_APMPI           *ptap;
1101:   Mat_MPIAIJ          *p=(Mat_MPIAIJ*)P->data,*c;
1102:   MPI_Comm            comm;
1103:   Mat                 Cmpi;
1104:   Mat_SeqAIJ          *pd,*po;
1105:   MatType             mtype;
1106:   PetscSF             sf;
1107:   PetscSFNode         *iremote;
1108:   PetscInt            rootspacesize,*rootspace,*rootspaceoffsets,nleaves;
1109:   const PetscInt      *rootdegrees;
1110:   PetscHSetI          ht,oht,*hta,*hto;
1111:   PetscInt            pn,pon,*c_rmtc,i,j,nzi,htsize,htosize,*c_rmtj,off,*c_othj,rcvncols,sendncols,*c_rmtoffsets;
1112:   PetscInt            owner,lidx,*rdj,col,pcstart,pcend,*dnz,*onz,am,arstart,arend,*poj,*pdj;
1113:   PetscInt            nalg=2,alg=0;
1114:   PetscBool           flg;
1115:   const char          *algTypes[2] = {"overlapping","merged"};
1116:   PetscErrorCode      ierr;

1119:   PetscObjectGetComm((PetscObject)A,&comm);

1121:   /* Create symbolic parallel matrix Cmpi */
1122:   MatGetLocalSize(P,NULL,&pn);
1123:   MatCreate(comm,&Cmpi);
1124:   MatGetType(A,&mtype);
1125:   MatSetType(Cmpi,mtype);
1126:   MatSetSizes(Cmpi,pn,pn,PETSC_DETERMINE,PETSC_DETERMINE);

1128:   PetscNew(&ptap);
1129:   ptap->reuse = MAT_INITIAL_MATRIX;
1130:   ptap->algType = 2;

1132:   /* Get P_oth by taking rows of P (= non-zero cols of local A) from other processors */
1133:   MatGetBrowsOfAoCols_MPIAIJ(A,P,MAT_INITIAL_MATRIX,&ptap->startsj_s,&ptap->startsj_r,&ptap->bufa,&ptap->P_oth);

1135:   po = (Mat_SeqAIJ*)p->B->data;
1136:   pd = (Mat_SeqAIJ*)p->A->data;

1138:   /* This equals to the number of offdiag columns in P */
1139:   MatGetLocalSize(p->B,NULL,&pon);
1140:   /* offsets */
1141:   PetscMalloc1(pon+1,&ptap->c_rmti);
1142:   /* The number of columns we will send to remote ranks */
1143:   PetscMalloc1(pon,&c_rmtc);
1144:   PetscMalloc1(pon,&hta);
1145:   for (i=0; i<pon; i++) {
1146:     PetscHSetICreate(&hta[i]);
1147:   }
1148:   MatGetLocalSize(A,&am,NULL);
1149:   MatGetOwnershipRange(A,&arstart,&arend);
1150:   /* Create hash table to merge all columns for C(i, :) */
1151:   PetscHSetICreate(&ht);

1153:   ptap->c_rmti[0] = 0;
1154:   /* 2) Pass 1: calculate the size for C_rmt (a matrix need to be sent to other processors)  */
1155:   for (i=0; i<am && pon; i++) {
1156:     /* Form one row of AP */
1157:     PetscHSetIClear(ht);
1158:     /* If the off diag is empty, we should not do any calculation */
1159:     nzi = po->i[i+1] - po->i[i];
1160:     if (!nzi) continue;

1162:     MatPtAPSymbolicComputeOneRowOfAP_private(A,P,ptap->P_oth,i,ht,ht);
1163:     PetscHSetIGetSize(ht,&htsize);
1164:     /* If AP is empty, just continue */
1165:     if (!htsize) continue;
1166:     /* Form C(ii, :) */
1167:     poj = po->j + po->i[i];
1168:     for (j=0; j<nzi; j++) {
1169:       PetscHSetIUpdate(hta[poj[j]],ht);
1170:     }
1171:   }

1173:   for (i=0; i<pon; i++) {
1174:     PetscHSetIGetSize(hta[i],&htsize);
1175:     ptap->c_rmti[i+1] = ptap->c_rmti[i] + htsize;
1176:     c_rmtc[i] = htsize;
1177:   }

1179:   PetscMalloc1(ptap->c_rmti[pon],&c_rmtj);

1181:   for (i=0; i<pon; i++) {
1182:     off = 0;
1183:     PetscHSetIGetElems(hta[i],&off,c_rmtj+ptap->c_rmti[i]);
1184:     PetscHSetIDestroy(&hta[i]);
1185:   }
1186:   PetscFree(hta);

1188:   PetscMalloc1(pon,&iremote);
1189:   for (i=0; i<pon; i++) {
1190:     owner = 0; lidx = 0;
1191:     PetscLayoutFindOwnerIndex(P->cmap,p->garray[i],&owner,&lidx);
1192:     iremote[i].index = lidx;
1193:     iremote[i].rank  = owner;
1194:   }

1196:   PetscSFCreate(comm,&sf);
1197:   PetscSFSetGraph(sf,pn,pon,NULL,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER);
1198:   /* Reorder ranks properly so that the data handled by gather and scatter have the same order */
1199:   PetscSFSetRankOrder(sf,PETSC_TRUE);
1200:   PetscSFSetFromOptions(sf);
1201:   PetscSFSetUp(sf);
1202:   /* How many neighbors have contributions to my rows? */
1203:   PetscSFComputeDegreeBegin(sf,&rootdegrees);
1204:   PetscSFComputeDegreeEnd(sf,&rootdegrees);
1205:   rootspacesize = 0;
1206:   for (i = 0; i < pn; i++) {
1207:     rootspacesize += rootdegrees[i];
1208:   }
1209:   PetscMalloc1(rootspacesize,&rootspace);
1210:   PetscMalloc1(rootspacesize+1,&rootspaceoffsets);
1211:   /* Get information from leaves
1212:    * Number of columns other people contribute to my rows
1213:    * */
1214:   PetscSFGatherBegin(sf,MPIU_INT,c_rmtc,rootspace);
1215:   PetscSFGatherEnd(sf,MPIU_INT,c_rmtc,rootspace);
1216:   PetscFree(c_rmtc);
1217:   PetscCalloc1(pn+1,&ptap->c_othi);
1218:   /* The number of columns is received for each row */
1219:   ptap->c_othi[0] = 0;
1220:   rootspacesize = 0;
1221:   rootspaceoffsets[0] = 0;
1222:   for (i = 0; i < pn; i++) {
1223:     rcvncols = 0;
1224:     for (j = 0; j<rootdegrees[i]; j++) {
1225:       rcvncols += rootspace[rootspacesize];
1226:       rootspaceoffsets[rootspacesize+1] = rootspaceoffsets[rootspacesize] + rootspace[rootspacesize];
1227:       rootspacesize++;
1228:     }
1229:     ptap->c_othi[i+1] = ptap->c_othi[i] + rcvncols;
1230:   }
1231:   PetscFree(rootspace);

1233:   PetscMalloc1(pon,&c_rmtoffsets);
1234:   PetscSFScatterBegin(sf,MPIU_INT,rootspaceoffsets,c_rmtoffsets);
1235:   PetscSFScatterEnd(sf,MPIU_INT,rootspaceoffsets,c_rmtoffsets);
1236:   PetscSFDestroy(&sf);
1237:   PetscFree(rootspaceoffsets);

1239:   PetscCalloc1(ptap->c_rmti[pon],&iremote);
1240:   nleaves = 0;
1241:   for (i = 0; i<pon; i++) {
1242:     owner = 0; lidx = 0;
1243:     PetscLayoutFindOwnerIndex(P->cmap,p->garray[i],&owner,&lidx);
1244:     sendncols = ptap->c_rmti[i+1] - ptap->c_rmti[i];
1245:     for (j=0; j<sendncols; j++) {
1246:       iremote[nleaves].rank = owner;
1247:       iremote[nleaves++].index = c_rmtoffsets[i] + j;
1248:     }
1249:   }
1250:   PetscFree(c_rmtoffsets);
1251:   PetscCalloc1(ptap->c_othi[pn],&c_othj);

1253:   PetscSFCreate(comm,&ptap->sf);
1254:   PetscSFSetGraph(ptap->sf,ptap->c_othi[pn],nleaves,NULL,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER);
1255:   PetscSFSetFromOptions(ptap->sf);
1256:   /* One to one map */
1257:   PetscSFReduceBegin(ptap->sf,MPIU_INT,c_rmtj,c_othj,MPIU_REPLACE);

1259:   PetscMalloc2(pn,&dnz,pn,&onz);
1260:   PetscHSetICreate(&oht);
1261:   MatGetOwnershipRangeColumn(P,&pcstart,&pcend);
1262:   PetscMalloc2(pn,&hta,pn,&hto);
1263:   for (i=0; i<pn; i++) {
1264:     PetscHSetICreate(&hta[i]);
1265:     PetscHSetICreate(&hto[i]);
1266:   }
1267:   /* Work on local part */
1268:   /* 4) Pass 1: Estimate memory for C_loc */
1269:   for (i=0; i<am && pn; i++) {
1270:     PetscHSetIClear(ht);
1271:     PetscHSetIClear(oht);
1272:     nzi = pd->i[i+1] - pd->i[i];
1273:     if (!nzi) continue;

1275:     MatPtAPSymbolicComputeOneRowOfAP_private(A,P,ptap->P_oth,i,ht,oht);
1276:     PetscHSetIGetSize(ht,&htsize);
1277:     PetscHSetIGetSize(oht,&htosize);
1278:     if (!(htsize+htosize)) continue;
1279:     /* Form C(ii, :) */
1280:     pdj = pd->j + pd->i[i];
1281:     for (j=0; j<nzi; j++) {
1282:       PetscHSetIUpdate(hta[pdj[j]],ht);
1283:       PetscHSetIUpdate(hto[pdj[j]],oht);
1284:     }
1285:   }

1287:   PetscHSetIDestroy(&ht);
1288:   PetscHSetIDestroy(&oht);

1290:   /* Get remote data */
1291:   PetscSFReduceEnd(ptap->sf,MPIU_INT,c_rmtj,c_othj,MPIU_REPLACE);
1292:   PetscFree(c_rmtj);

1294:   for (i = 0; i < pn; i++) {
1295:     nzi = ptap->c_othi[i+1] - ptap->c_othi[i];
1296:     rdj = c_othj + ptap->c_othi[i];
1297:     for (j = 0; j < nzi; j++) {
1298:       col =  rdj[j];
1299:       /* diag part */
1300:       if (col>=pcstart && col<pcend) {
1301:         PetscHSetIAdd(hta[i],col);
1302:       } else { /* off diag */
1303:         PetscHSetIAdd(hto[i],col);
1304:       }
1305:     }
1306:     PetscHSetIGetSize(hta[i],&htsize);
1307:     dnz[i] = htsize;
1308:     PetscHSetIDestroy(&hta[i]);
1309:     PetscHSetIGetSize(hto[i],&htsize);
1310:     onz[i] = htsize;
1311:     PetscHSetIDestroy(&hto[i]);
1312:   }

1314:   PetscFree2(hta,hto);
1315:   PetscFree(c_othj);

1317:   /* local sizes and preallocation */
1318:   MatSetSizes(Cmpi,pn,pn,PETSC_DETERMINE,PETSC_DETERMINE);
1319:   MatSetBlockSizes(Cmpi,PetscAbs(P->cmap->bs),PetscAbs(P->cmap->bs));
1320:   MatMPIAIJSetPreallocation(Cmpi,0,dnz,0,onz);
1321:   PetscFree2(dnz,onz);

1323:   /* attach the supporting struct to Cmpi for reuse */
1324:   c = (Mat_MPIAIJ*)Cmpi->data;
1325:   c->ap           = ptap;
1326:   ptap->duplicate = Cmpi->ops->duplicate;
1327:   ptap->destroy   = Cmpi->ops->destroy;
1328:   ptap->view      = Cmpi->ops->view;

1330:   /* Cmpi is not ready for use - assembly will be done by MatPtAPNumeric() */
1331:   Cmpi->assembled        = PETSC_FALSE;
1332:   /* pick an algorithm */
1333:   PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MatPtAP","Mat");
1334:   alg = 0;
1335:   PetscOptionsEList("-matptap_allatonce_via","PtAP allatonce numeric approach","MatPtAP",algTypes,nalg,algTypes[alg],&alg,&flg);
1336:   PetscOptionsEnd();
1337:   switch (alg) {
1338:     case 0:
1339:       Cmpi->ops->ptapnumeric = MatPtAPNumeric_MPIAIJ_MPIAIJ_allatonce;
1340:       break;
1341:     case 1:
1342:       Cmpi->ops->ptapnumeric = MatPtAPNumeric_MPIAIJ_MPIAIJ_allatonce_merged;
1343:       break;
1344:     default:
1345:       SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG," Unsupported allatonce numerical algorithm \n");
1346:   }
1347:   Cmpi->ops->destroy     = MatDestroy_MPIAIJ_PtAP;
1348:   Cmpi->ops->view        = MatView_MPIAIJ_PtAP;
1349:   Cmpi->ops->freeintermediatedatastructures = MatFreeIntermediateDataStructures_MPIAIJ_AP;
1350:   *C                     = Cmpi;
1351:   return(0);
1352: }

1354: PetscErrorCode MatPtAPSymbolic_MPIAIJ_MPIAIJ_allatonce_merged(Mat A,Mat P,PetscReal fill,Mat *C)
1355: {
1356:   Mat_APMPI           *ptap;
1357:   Mat_MPIAIJ          *p=(Mat_MPIAIJ*)P->data,*c;
1358:   MPI_Comm            comm;
1359:   Mat                 Cmpi;
1360:   Mat_SeqAIJ          *pd,*po;
1361:   MatType             mtype;
1362:   PetscSF             sf;
1363:   PetscSFNode         *iremote;
1364:   PetscInt            rootspacesize,*rootspace,*rootspaceoffsets,nleaves;
1365:   const PetscInt      *rootdegrees;
1366:   PetscHSetI          ht,oht,*hta,*hto,*htd;
1367:   PetscInt            pn,pon,*c_rmtc,i,j,nzi,dnzi,htsize,htosize,*c_rmtj,off,*c_othj,rcvncols,sendncols,*c_rmtoffsets;
1368:   PetscInt            owner,lidx,*rdj,col,pcstart,pcend,*dnz,*onz,am,arstart,arend,*poj,*pdj;
1369:   PetscInt            nalg=2,alg=0;
1370:   PetscBool           flg;
1371:   const char          *algTypes[2] = {"merged","overlapping"};
1372:   PetscErrorCode      ierr;

1375:   PetscObjectGetComm((PetscObject)A,&comm);

1377:   /* Create symbolic parallel matrix Cmpi */
1378:   MatGetLocalSize(P,NULL,&pn);
1379:   MatCreate(comm,&Cmpi);
1380:   MatGetType(A,&mtype);
1381:   MatSetType(Cmpi,mtype);
1382:   MatSetSizes(Cmpi,pn,pn,PETSC_DETERMINE,PETSC_DETERMINE);

1384:   PetscNew(&ptap);
1385:   ptap->reuse = MAT_INITIAL_MATRIX;
1386:   ptap->algType = 3;

1388:   /* 0) Get P_oth by taking rows of P (= non-zero cols of local A) from other processors */
1389:   MatGetBrowsOfAoCols_MPIAIJ(A,P,MAT_INITIAL_MATRIX,&ptap->startsj_s,&ptap->startsj_r,&ptap->bufa,&ptap->P_oth);

1391:   po = (Mat_SeqAIJ*)p->B->data;
1392:   pd = (Mat_SeqAIJ*)p->A->data;

1394:   /* This equals to the number of offdiag columns in P */
1395:   MatGetLocalSize(p->B,NULL,&pon);
1396:   /* offsets */
1397:   PetscMalloc1(pon+1,&ptap->c_rmti);
1398:   /* The number of columns we will send to remote ranks */
1399:   PetscMalloc1(pon,&c_rmtc);
1400:   PetscMalloc1(pon,&hta);
1401:   for (i=0; i<pon; i++) {
1402:     PetscHSetICreate(&hta[i]);
1403:   }
1404:   MatGetLocalSize(A,&am,NULL);
1405:   MatGetOwnershipRange(A,&arstart,&arend);
1406:   /* Create hash table to merge all columns for C(i, :) */
1407:   PetscHSetICreate(&ht);
1408:   PetscHSetICreate(&oht);
1409:   PetscMalloc2(pn,&htd,pn,&hto);
1410:   for (i=0; i<pn; i++) {
1411:     PetscHSetICreate(&htd[i]);
1412:     PetscHSetICreate(&hto[i]);
1413:   }
1414:   ptap->c_rmti[0] = 0;
1415:   /* 2) Pass 1: calculate the size for C_rmt (a matrix need to be sent to other processors)  */
1416:   for (i=0; i<am && (pon || pn); i++) {
1417:     /* Form one row of AP */
1418:     PetscHSetIClear(ht);
1419:     PetscHSetIClear(oht);
1420:     /* If the off diag is empty, we should not do any calculation */
1421:     nzi = po->i[i+1] - po->i[i];
1422:     dnzi = pd->i[i+1] - pd->i[i];
1423:     if (!nzi && !dnzi) continue;

1425:     MatPtAPSymbolicComputeOneRowOfAP_private(A,P,ptap->P_oth,i,ht,oht);
1426:     PetscHSetIGetSize(ht,&htsize);
1427:     PetscHSetIGetSize(oht,&htosize);
1428:     /* If AP is empty, just continue */
1429:     if (!(htsize+htosize)) continue;

1431:     /* Form remote C(ii, :) */
1432:     poj = po->j + po->i[i];
1433:     for (j=0; j<nzi; j++) {
1434:       PetscHSetIUpdate(hta[poj[j]],ht);
1435:       PetscHSetIUpdate(hta[poj[j]],oht);
1436:     }

1438:     /* Form local C(ii, :) */
1439:     pdj = pd->j + pd->i[i];
1440:     for (j=0; j<dnzi; j++) {
1441:       PetscHSetIUpdate(htd[pdj[j]],ht);
1442:       PetscHSetIUpdate(hto[pdj[j]],oht);
1443:     }
1444:   }

1446:   PetscHSetIDestroy(&ht);
1447:   PetscHSetIDestroy(&oht);

1449:   for (i=0; i<pon; i++) {
1450:     PetscHSetIGetSize(hta[i],&htsize);
1451:     ptap->c_rmti[i+1] = ptap->c_rmti[i] + htsize;
1452:     c_rmtc[i] = htsize;
1453:   }

1455:   PetscMalloc1(ptap->c_rmti[pon],&c_rmtj);

1457:   for (i=0; i<pon; i++) {
1458:     off = 0;
1459:     PetscHSetIGetElems(hta[i],&off,c_rmtj+ptap->c_rmti[i]);
1460:     PetscHSetIDestroy(&hta[i]);
1461:   }
1462:   PetscFree(hta);

1464:   PetscMalloc1(pon,&iremote);
1465:   for (i=0; i<pon; i++) {
1466:     owner = 0; lidx = 0;
1467:     PetscLayoutFindOwnerIndex(P->cmap,p->garray[i],&owner,&lidx);
1468:     iremote[i].index = lidx;
1469:     iremote[i].rank  = owner;
1470:   }

1472:   PetscSFCreate(comm,&sf);
1473:   PetscSFSetGraph(sf,pn,pon,NULL,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER);
1474:   /* Reorder ranks properly so that the data handled by gather and scatter have the same order */
1475:   PetscSFSetRankOrder(sf,PETSC_TRUE);
1476:   PetscSFSetFromOptions(sf);
1477:   PetscSFSetUp(sf);
1478:   /* How many neighbors have contributions to my rows? */
1479:   PetscSFComputeDegreeBegin(sf,&rootdegrees);
1480:   PetscSFComputeDegreeEnd(sf,&rootdegrees);
1481:   rootspacesize = 0;
1482:   for (i = 0; i < pn; i++) {
1483:     rootspacesize += rootdegrees[i];
1484:   }
1485:   PetscMalloc1(rootspacesize,&rootspace);
1486:   PetscMalloc1(rootspacesize+1,&rootspaceoffsets);
1487:   /* Get information from leaves
1488:    * Number of columns other people contribute to my rows
1489:    * */
1490:   PetscSFGatherBegin(sf,MPIU_INT,c_rmtc,rootspace);
1491:   PetscSFGatherEnd(sf,MPIU_INT,c_rmtc,rootspace);
1492:   PetscFree(c_rmtc);
1493:   PetscMalloc1(pn+1,&ptap->c_othi);
1494:   /* The number of columns is received for each row */
1495:   ptap->c_othi[0]     = 0;
1496:   rootspacesize       = 0;
1497:   rootspaceoffsets[0] = 0;
1498:   for (i = 0; i < pn; i++) {
1499:     rcvncols = 0;
1500:     for (j = 0; j<rootdegrees[i]; j++) {
1501:       rcvncols += rootspace[rootspacesize];
1502:       rootspaceoffsets[rootspacesize+1] = rootspaceoffsets[rootspacesize] + rootspace[rootspacesize];
1503:       rootspacesize++;
1504:     }
1505:     ptap->c_othi[i+1] = ptap->c_othi[i] + rcvncols;
1506:   }
1507:   PetscFree(rootspace);

1509:   PetscMalloc1(pon,&c_rmtoffsets);
1510:   PetscSFScatterBegin(sf,MPIU_INT,rootspaceoffsets,c_rmtoffsets);
1511:   PetscSFScatterEnd(sf,MPIU_INT,rootspaceoffsets,c_rmtoffsets);
1512:   PetscSFDestroy(&sf);
1513:   PetscFree(rootspaceoffsets);

1515:   PetscCalloc1(ptap->c_rmti[pon],&iremote);
1516:   nleaves = 0;
1517:   for (i = 0; i<pon; i++) {
1518:     owner = 0;
1519:     lidx  = 0;
1520:     PetscLayoutFindOwnerIndex(P->cmap,p->garray[i],&owner,&lidx);
1521:     sendncols = ptap->c_rmti[i+1] - ptap->c_rmti[i];
1522:     for (j=0; j<sendncols; j++) {
1523:       iremote[nleaves].rank    = owner;
1524:       iremote[nleaves++].index = c_rmtoffsets[i] + j;
1525:     }
1526:   }
1527:   PetscFree(c_rmtoffsets);
1528:   PetscCalloc1(ptap->c_othi[pn],&c_othj);

1530:   PetscSFCreate(comm,&ptap->sf);
1531:   PetscSFSetGraph(ptap->sf,ptap->c_othi[pn],nleaves,NULL,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER);
1532:   PetscSFSetFromOptions(ptap->sf);
1533:   /* One to one map */
1534:   PetscSFReduceBegin(ptap->sf,MPIU_INT,c_rmtj,c_othj,MPIU_REPLACE);
1535:   /* Get remote data */
1536:   PetscSFReduceEnd(ptap->sf,MPIU_INT,c_rmtj,c_othj,MPIU_REPLACE);
1537:   PetscFree(c_rmtj);
1538:   PetscMalloc2(pn,&dnz,pn,&onz);
1539:   MatGetOwnershipRangeColumn(P,&pcstart,&pcend);

1541:   for (i = 0; i < pn; i++) {
1542:     nzi = ptap->c_othi[i+1] - ptap->c_othi[i];
1543:     rdj = c_othj + ptap->c_othi[i];
1544:     for (j = 0; j < nzi; j++) {
1545:       col =  rdj[j];
1546:       /* diag part */
1547:       if (col>=pcstart && col<pcend) {
1548:         PetscHSetIAdd(htd[i],col);
1549:       } else { /* off diag */
1550:         PetscHSetIAdd(hto[i],col);
1551:       }
1552:     }
1553:     PetscHSetIGetSize(htd[i],&htsize);
1554:     dnz[i] = htsize;
1555:     PetscHSetIDestroy(&htd[i]);
1556:     PetscHSetIGetSize(hto[i],&htsize);
1557:     onz[i] = htsize;
1558:     PetscHSetIDestroy(&hto[i]);
1559:   }

1561:   PetscFree2(htd,hto);
1562:   PetscFree(c_othj);

1564:   /* local sizes and preallocation */
1565:   MatSetSizes(Cmpi,pn,pn,PETSC_DETERMINE,PETSC_DETERMINE);
1566:   MatSetBlockSizes(Cmpi,PetscAbs(P->cmap->bs),PetscAbs(P->cmap->bs));
1567:   MatMPIAIJSetPreallocation(Cmpi,0,dnz,0,onz);
1568:   PetscFree2(dnz,onz);

1570:   /* attach the supporting struct to Cmpi for reuse */
1571:   c = (Mat_MPIAIJ*)Cmpi->data;
1572:   c->ap           = ptap;
1573:   ptap->duplicate = Cmpi->ops->duplicate;
1574:   ptap->destroy   = Cmpi->ops->destroy;
1575:   ptap->view      = Cmpi->ops->view;

1577:   /* Cmpi is not ready for use - assembly will be done by MatPtAPNumeric() */
1578:   Cmpi->assembled        = PETSC_FALSE;
1579:   /* pick an algorithm */
1580:   PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MatPtAP","Mat");
1581:   alg = 0;
1582:   PetscOptionsEList("-matptap_allatonce_via","PtAP allatonce numeric approach","MatPtAP",algTypes,nalg,algTypes[alg],&alg,&flg);
1583:   PetscOptionsEnd();
1584:   switch (alg) {
1585:     case 0:
1586:       Cmpi->ops->ptapnumeric = MatPtAPNumeric_MPIAIJ_MPIAIJ_allatonce_merged;
1587:       break;
1588:     case 1:
1589:       Cmpi->ops->ptapnumeric = MatPtAPNumeric_MPIAIJ_MPIAIJ_allatonce;
1590:       break;
1591:     default:
1592:       SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG," Unsupported allatonce numerical algorithm \n");
1593:   }
1594:   Cmpi->ops->destroy     = MatDestroy_MPIAIJ_PtAP;
1595:   Cmpi->ops->view        = MatView_MPIAIJ_PtAP;
1596:   Cmpi->ops->freeintermediatedatastructures = MatFreeIntermediateDataStructures_MPIAIJ_AP;
1597:   *C                     = Cmpi;
1598:   return(0);
1599: }

1601: PetscErrorCode MatPtAPSymbolic_MPIAIJ_MPIAIJ(Mat A,Mat P,PetscReal fill,Mat *C)
1602: {
1603:   PetscErrorCode      ierr;
1604:   Mat_APMPI           *ptap;
1605:   Mat_MPIAIJ          *a=(Mat_MPIAIJ*)A->data,*p=(Mat_MPIAIJ*)P->data,*c;
1606:   MPI_Comm            comm;
1607:   PetscMPIInt         size,rank;
1608:   Mat                 Cmpi;
1609:   PetscFreeSpaceList  free_space=NULL,current_space=NULL;
1610:   PetscInt            am=A->rmap->n,pm=P->rmap->n,pN=P->cmap->N,pn=P->cmap->n;
1611:   PetscInt            *lnk,i,k,pnz,row,nsend;
1612:   PetscBT             lnkbt;
1613:   PetscMPIInt         tagi,tagj,*len_si,*len_s,*len_ri,icompleted=0,nrecv;
1614:   PetscInt            **buf_rj,**buf_ri,**buf_ri_k;
1615:   PetscInt            len,proc,*dnz,*onz,*owners,nzi,nspacedouble;
1616:   PetscInt            nrows,*buf_s,*buf_si,*buf_si_i,**nextrow,**nextci;
1617:   MPI_Request         *swaits,*rwaits;
1618:   MPI_Status          *sstatus,rstatus;
1619:   PetscLayout         rowmap;
1620:   PetscInt            *owners_co,*coi,*coj;    /* i and j array of (p->B)^T*A*P - used in the communication */
1621:   PetscMPIInt         *len_r,*id_r;    /* array of length of comm->size, store send/recv matrix values */
1622:   PetscInt            *api,*apj,*Jptr,apnz,*prmap=p->garray,con,j,ap_rmax=0,Crmax,*aj,*ai,*pi;
1623:   Mat_SeqAIJ          *p_loc,*p_oth=NULL,*ad=(Mat_SeqAIJ*)(a->A)->data,*ao=NULL,*c_loc,*c_oth;
1624:   PetscScalar         *apv;
1625:   PetscTable          ta;
1626:   MatType             mtype;
1627:   const char          *prefix;
1628: #if defined(PETSC_USE_INFO)
1629:   PetscReal           apfill;
1630: #endif

1633:   PetscObjectGetComm((PetscObject)A,&comm);
1634:   MPI_Comm_size(comm,&size);
1635:   MPI_Comm_rank(comm,&rank);

1637:   if (size > 1) ao = (Mat_SeqAIJ*)(a->B)->data;

1639:   /* create symbolic parallel matrix Cmpi */
1640:   MatCreate(comm,&Cmpi);
1641:   MatGetType(A,&mtype);
1642:   MatSetType(Cmpi,mtype);

1644:   /* Do dense axpy in MatPtAPNumeric_MPIAIJ_MPIAIJ() */
1645:   Cmpi->ops->ptapnumeric = MatPtAPNumeric_MPIAIJ_MPIAIJ;

1647:   /* create struct Mat_APMPI and attached it to C later */
1648:   PetscNew(&ptap);
1649:   ptap->reuse = MAT_INITIAL_MATRIX;
1650:   ptap->algType = 1;

1652:   /* get P_oth by taking rows of P (= non-zero cols of local A) from other processors */
1653:   MatGetBrowsOfAoCols_MPIAIJ(A,P,MAT_INITIAL_MATRIX,&ptap->startsj_s,&ptap->startsj_r,&ptap->bufa,&ptap->P_oth);
1654:   /* get P_loc by taking all local rows of P */
1655:   MatMPIAIJGetLocalMat(P,MAT_INITIAL_MATRIX,&ptap->P_loc);

1657:   /* (0) compute Rd = Pd^T, Ro = Po^T  */
1658:   /* --------------------------------- */
1659:   MatTranspose(p->A,MAT_INITIAL_MATRIX,&ptap->Rd);
1660:   MatTranspose(p->B,MAT_INITIAL_MATRIX,&ptap->Ro);

1662:   /* (1) compute symbolic AP = A_loc*P = Ad*P_loc + Ao*P_oth (api,apj) */
1663:   /* ----------------------------------------------------------------- */
1664:   p_loc  = (Mat_SeqAIJ*)(ptap->P_loc)->data;
1665:   if (ptap->P_oth) p_oth  = (Mat_SeqAIJ*)(ptap->P_oth)->data;

1667:   /* create and initialize a linked list */
1668:   PetscTableCreate(pn,pN,&ta); /* for compute AP_loc and Cmpi */
1669:   MatRowMergeMax_SeqAIJ(p_loc,ptap->P_loc->rmap->N,ta);
1670:   MatRowMergeMax_SeqAIJ(p_oth,ptap->P_oth->rmap->N,ta);
1671:   PetscTableGetCount(ta,&Crmax); /* Crmax = nnz(sum of Prows) */
1672:   /* printf("[%d] est %d, Crmax %d; pN %d\n",rank,5*(p_loc->rmax+p_oth->rmax + (PetscInt)(1.e-2*pN)),Crmax,pN); */

1674:   PetscLLCondensedCreate(Crmax,pN,&lnk,&lnkbt);

1676:   /* Initial FreeSpace size is fill*(nnz(A) + nnz(P)) */
1677:   if (ao) {
1678:     PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,PetscIntSumTruncate(ad->i[am],PetscIntSumTruncate(ao->i[am],p_loc->i[pm]))),&free_space);
1679:   } else {
1680:     PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,PetscIntSumTruncate(ad->i[am],p_loc->i[pm])),&free_space);
1681:   }
1682:   current_space = free_space;
1683:   nspacedouble  = 0;

1685:   PetscMalloc1(am+1,&api);
1686:   api[0] = 0;
1687:   for (i=0; i<am; i++) {
1688:     /* diagonal portion: Ad[i,:]*P */
1689:     ai = ad->i; pi = p_loc->i;
1690:     nzi = ai[i+1] - ai[i];
1691:     aj  = ad->j + ai[i];
1692:     for (j=0; j<nzi; j++) {
1693:       row  = aj[j];
1694:       pnz  = pi[row+1] - pi[row];
1695:       Jptr = p_loc->j + pi[row];
1696:       /* add non-zero cols of P into the sorted linked list lnk */
1697:       PetscLLCondensedAddSorted(pnz,Jptr,lnk,lnkbt);
1698:     }
1699:     /* off-diagonal portion: Ao[i,:]*P */
1700:     if (ao) {
1701:       ai = ao->i; pi = p_oth->i;
1702:       nzi = ai[i+1] - ai[i];
1703:       aj  = ao->j + ai[i];
1704:       for (j=0; j<nzi; j++) {
1705:         row  = aj[j];
1706:         pnz  = pi[row+1] - pi[row];
1707:         Jptr = p_oth->j + pi[row];
1708:         PetscLLCondensedAddSorted(pnz,Jptr,lnk,lnkbt);
1709:       }
1710:     }
1711:     apnz     = lnk[0];
1712:     api[i+1] = api[i] + apnz;
1713:     if (ap_rmax < apnz) ap_rmax = apnz;

1715:     /* if free space is not available, double the total space in the list */
1716:     if (current_space->local_remaining<apnz) {
1717:       PetscFreeSpaceGet(PetscIntSumTruncate(apnz,current_space->total_array_size),&current_space);
1718:       nspacedouble++;
1719:     }

1721:     /* Copy data into free space, then initialize lnk */
1722:     PetscLLCondensedClean(pN,apnz,current_space->array,lnk,lnkbt);

1724:     current_space->array           += apnz;
1725:     current_space->local_used      += apnz;
1726:     current_space->local_remaining -= apnz;
1727:   }
1728:   /* Allocate space for apj and apv, initialize apj, and */
1729:   /* destroy list of free space and other temporary array(s) */
1730:   PetscMalloc2(api[am],&apj,api[am],&apv);
1731:   PetscFreeSpaceContiguous(&free_space,apj);
1732:   PetscLLDestroy(lnk,lnkbt);

1734:   /* Create AP_loc for reuse */
1735:   MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,am,pN,api,apj,apv,&ptap->AP_loc);

1737: #if defined(PETSC_USE_INFO)
1738:   if (ao) {
1739:     apfill = (PetscReal)api[am]/(ad->i[am]+ao->i[am]+p_loc->i[pm]+1);
1740:   } else {
1741:     apfill = (PetscReal)api[am]/(ad->i[am]+p_loc->i[pm]+1);
1742:   }
1743:   ptap->AP_loc->info.mallocs           = nspacedouble;
1744:   ptap->AP_loc->info.fill_ratio_given  = fill;
1745:   ptap->AP_loc->info.fill_ratio_needed = apfill;

1747:   if (api[am]) {
1748:     PetscInfo3(ptap->AP_loc,"Nonscalable algorithm, AP_loc reallocs %D; Fill ratio: given %g needed %g.\n",nspacedouble,(double)fill,(double)apfill);
1749:     PetscInfo1(ptap->AP_loc,"Use MatPtAP(A,B,MatReuse,%g,&C) for best AP_loc performance.;\n",(double)apfill);
1750:   } else {
1751:     PetscInfo(ptap->AP_loc,"Nonscalable algorithm, AP_loc is empty \n");
1752:   }
1753: #endif

1755:   /* (2-1) compute symbolic Co = Ro*AP_loc  */
1756:   /* ------------------------------------ */
1757:   MatGetOptionsPrefix(A,&prefix);
1758:   MatSetOptionsPrefix(ptap->Ro,prefix);
1759:   MatAppendOptionsPrefix(ptap->Ro,"inner_offdiag_");
1760:   MatMatMultSymbolic_SeqAIJ_SeqAIJ(ptap->Ro,ptap->AP_loc,fill,&ptap->C_oth);

1762:   /* (3) send coj of C_oth to other processors  */
1763:   /* ------------------------------------------ */
1764:   /* determine row ownership */
1765:   PetscLayoutCreate(comm,&rowmap);
1766:   rowmap->n  = pn;
1767:   rowmap->bs = 1;
1768:   PetscLayoutSetUp(rowmap);
1769:   owners = rowmap->range;

1771:   /* determine the number of messages to send, their lengths */
1772:   PetscMalloc4(size,&len_s,size,&len_si,size,&sstatus,size+2,&owners_co);
1773:   PetscArrayzero(len_s,size);
1774:   PetscArrayzero(len_si,size);

1776:   c_oth = (Mat_SeqAIJ*)ptap->C_oth->data;
1777:   coi   = c_oth->i; coj = c_oth->j;
1778:   con   = ptap->C_oth->rmap->n;
1779:   proc  = 0;
1780:   for (i=0; i<con; i++) {
1781:     while (prmap[i] >= owners[proc+1]) proc++;
1782:     len_si[proc]++;               /* num of rows in Co(=Pt*AP) to be sent to [proc] */
1783:     len_s[proc] += coi[i+1] - coi[i]; /* num of nonzeros in Co to be sent to [proc] */
1784:   }

1786:   len          = 0; /* max length of buf_si[], see (4) */
1787:   owners_co[0] = 0;
1788:   nsend        = 0;
1789:   for (proc=0; proc<size; proc++) {
1790:     owners_co[proc+1] = owners_co[proc] + len_si[proc];
1791:     if (len_s[proc]) {
1792:       nsend++;
1793:       len_si[proc] = 2*(len_si[proc] + 1); /* length of buf_si to be sent to [proc] */
1794:       len         += len_si[proc];
1795:     }
1796:   }

1798:   /* determine the number and length of messages to receive for coi and coj  */
1799:   PetscGatherNumberOfMessages(comm,NULL,len_s,&nrecv);
1800:   PetscGatherMessageLengths2(comm,nsend,nrecv,len_s,len_si,&id_r,&len_r,&len_ri);

1802:   /* post the Irecv and Isend of coj */
1803:   PetscCommGetNewTag(comm,&tagj);
1804:   PetscPostIrecvInt(comm,tagj,nrecv,id_r,len_r,&buf_rj,&rwaits);
1805:   PetscMalloc1(nsend+1,&swaits);
1806:   for (proc=0, k=0; proc<size; proc++) {
1807:     if (!len_s[proc]) continue;
1808:     i    = owners_co[proc];
1809:     MPI_Isend(coj+coi[i],len_s[proc],MPIU_INT,proc,tagj,comm,swaits+k);
1810:     k++;
1811:   }

1813:   /* (2-2) compute symbolic C_loc = Rd*AP_loc */
1814:   /* ---------------------------------------- */
1815:   MatSetOptionsPrefix(ptap->Rd,prefix);
1816:   MatAppendOptionsPrefix(ptap->Rd,"inner_diag_");
1817:   MatMatMultSymbolic_SeqAIJ_SeqAIJ(ptap->Rd,ptap->AP_loc,fill,&ptap->C_loc);
1818:   c_loc = (Mat_SeqAIJ*)ptap->C_loc->data;

1820:   /* receives coj are complete */
1821:   for (i=0; i<nrecv; i++) {
1822:     MPI_Waitany(nrecv,rwaits,&icompleted,&rstatus);
1823:   }
1824:   PetscFree(rwaits);
1825:   if (nsend) {MPI_Waitall(nsend,swaits,sstatus);}

1827:   /* add received column indices into ta to update Crmax */
1828:   for (k=0; k<nrecv; k++) {/* k-th received message */
1829:     Jptr = buf_rj[k];
1830:     for (j=0; j<len_r[k]; j++) {
1831:       PetscTableAdd(ta,*(Jptr+j)+1,1,INSERT_VALUES);
1832:     }
1833:   }
1834:   PetscTableGetCount(ta,&Crmax);
1835:   PetscTableDestroy(&ta);

1837:   /* (4) send and recv coi */
1838:   /*-----------------------*/
1839:   PetscCommGetNewTag(comm,&tagi);
1840:   PetscPostIrecvInt(comm,tagi,nrecv,id_r,len_ri,&buf_ri,&rwaits);
1841:   PetscMalloc1(len+1,&buf_s);
1842:   buf_si = buf_s;  /* points to the beginning of k-th msg to be sent */
1843:   for (proc=0,k=0; proc<size; proc++) {
1844:     if (!len_s[proc]) continue;
1845:     /* form outgoing message for i-structure:
1846:          buf_si[0]:                 nrows to be sent
1847:                [1:nrows]:           row index (global)
1848:                [nrows+1:2*nrows+1]: i-structure index
1849:     */
1850:     /*-------------------------------------------*/
1851:     nrows       = len_si[proc]/2 - 1; /* num of rows in Co to be sent to [proc] */
1852:     buf_si_i    = buf_si + nrows+1;
1853:     buf_si[0]   = nrows;
1854:     buf_si_i[0] = 0;
1855:     nrows       = 0;
1856:     for (i=owners_co[proc]; i<owners_co[proc+1]; i++) {
1857:       nzi = coi[i+1] - coi[i];
1858:       buf_si_i[nrows+1] = buf_si_i[nrows] + nzi;  /* i-structure */
1859:       buf_si[nrows+1]   = prmap[i] -owners[proc]; /* local row index */
1860:       nrows++;
1861:     }
1862:     MPI_Isend(buf_si,len_si[proc],MPIU_INT,proc,tagi,comm,swaits+k);
1863:     k++;
1864:     buf_si += len_si[proc];
1865:   }
1866:   for (i=0; i<nrecv; i++) {
1867:     MPI_Waitany(nrecv,rwaits,&icompleted,&rstatus);
1868:   }
1869:   PetscFree(rwaits);
1870:   if (nsend) {MPI_Waitall(nsend,swaits,sstatus);}

1872:   PetscFree4(len_s,len_si,sstatus,owners_co);
1873:   PetscFree(len_ri);
1874:   PetscFree(swaits);
1875:   PetscFree(buf_s);

1877:   /* (5) compute the local portion of Cmpi      */
1878:   /* ------------------------------------------ */
1879:   /* set initial free space to be Crmax, sufficient for holding nozeros in each row of Cmpi */
1880:   PetscFreeSpaceGet(Crmax,&free_space);
1881:   current_space = free_space;

1883:   PetscMalloc3(nrecv,&buf_ri_k,nrecv,&nextrow,nrecv,&nextci);
1884:   for (k=0; k<nrecv; k++) {
1885:     buf_ri_k[k] = buf_ri[k]; /* beginning of k-th recved i-structure */
1886:     nrows       = *buf_ri_k[k];
1887:     nextrow[k]  = buf_ri_k[k] + 1;  /* next row number of k-th recved i-structure */
1888:     nextci[k]   = buf_ri_k[k] + (nrows + 1); /* poins to the next i-structure of k-th recved i-structure  */
1889:   }

1891:   MatPreallocateInitialize(comm,pn,pn,dnz,onz);
1892:   PetscLLCondensedCreate(Crmax,pN,&lnk,&lnkbt);
1893:   for (i=0; i<pn; i++) {
1894:     /* add C_loc into Cmpi */
1895:     nzi  = c_loc->i[i+1] - c_loc->i[i];
1896:     Jptr = c_loc->j + c_loc->i[i];
1897:     PetscLLCondensedAddSorted(nzi,Jptr,lnk,lnkbt);

1899:     /* add received col data into lnk */
1900:     for (k=0; k<nrecv; k++) { /* k-th received message */
1901:       if (i == *nextrow[k]) { /* i-th row */
1902:         nzi  = *(nextci[k]+1) - *nextci[k];
1903:         Jptr = buf_rj[k] + *nextci[k];
1904:         PetscLLCondensedAddSorted(nzi,Jptr,lnk,lnkbt);
1905:         nextrow[k]++; nextci[k]++;
1906:       }
1907:     }
1908:     nzi = lnk[0];

1910:     /* copy data into free space, then initialize lnk */
1911:     PetscLLCondensedClean(pN,nzi,current_space->array,lnk,lnkbt);
1912:     MatPreallocateSet(i+owners[rank],nzi,current_space->array,dnz,onz);
1913:   }
1914:   PetscFree3(buf_ri_k,nextrow,nextci);
1915:   PetscLLDestroy(lnk,lnkbt);
1916:   PetscFreeSpaceDestroy(free_space);

1918:   /* local sizes and preallocation */
1919:   MatSetSizes(Cmpi,pn,pn,PETSC_DETERMINE,PETSC_DETERMINE);
1920:   if (P->cmap->bs > 0) {
1921:     PetscLayoutSetBlockSize(Cmpi->rmap,P->cmap->bs);
1922:     PetscLayoutSetBlockSize(Cmpi->cmap,P->cmap->bs);
1923:   }
1924:   MatMPIAIJSetPreallocation(Cmpi,0,dnz,0,onz);
1925:   MatPreallocateFinalize(dnz,onz);

1927:   /* members in merge */
1928:   PetscFree(id_r);
1929:   PetscFree(len_r);
1930:   PetscFree(buf_ri[0]);
1931:   PetscFree(buf_ri);
1932:   PetscFree(buf_rj[0]);
1933:   PetscFree(buf_rj);
1934:   PetscLayoutDestroy(&rowmap);

1936:   /* attach the supporting struct to Cmpi for reuse */
1937:   c = (Mat_MPIAIJ*)Cmpi->data;
1938:   c->ap           = ptap;
1939:   ptap->duplicate = Cmpi->ops->duplicate;
1940:   ptap->destroy   = Cmpi->ops->destroy;
1941:   ptap->view      = Cmpi->ops->view;
1942:   PetscCalloc1(pN,&ptap->apa);

1944:   /* Cmpi is not ready for use - assembly will be done by MatPtAPNumeric() */
1945:   Cmpi->assembled        = PETSC_FALSE;
1946:   Cmpi->ops->destroy     = MatDestroy_MPIAIJ_PtAP;
1947:   Cmpi->ops->view        = MatView_MPIAIJ_PtAP;
1948:   Cmpi->ops->freeintermediatedatastructures = MatFreeIntermediateDataStructures_MPIAIJ_AP;
1949:   *C                     = Cmpi;
1950:   return(0);
1951: }

1953: PetscErrorCode MatPtAPNumeric_MPIAIJ_MPIAIJ(Mat A,Mat P,Mat C)
1954: {
1955:   PetscErrorCode    ierr;
1956:   Mat_MPIAIJ        *a=(Mat_MPIAIJ*)A->data,*p=(Mat_MPIAIJ*)P->data,*c=(Mat_MPIAIJ*)C->data;
1957:   Mat_SeqAIJ        *ad=(Mat_SeqAIJ*)(a->A)->data,*ao=(Mat_SeqAIJ*)(a->B)->data;
1958:   Mat_SeqAIJ        *ap,*p_loc,*p_oth=NULL,*c_seq;
1959:   Mat_APMPI         *ptap = c->ap;
1960:   Mat               AP_loc,C_loc,C_oth;
1961:   PetscInt          i,rstart,rend,cm,ncols,row;
1962:   PetscInt          *api,*apj,am = A->rmap->n,j,col,apnz;
1963:   PetscScalar       *apa;
1964:   const PetscInt    *cols;
1965:   const PetscScalar *vals;

1968:   if (!ptap->AP_loc) {
1969:     MPI_Comm comm;
1970:     PetscObjectGetComm((PetscObject)C,&comm);
1971:     SETERRQ(comm,PETSC_ERR_ARG_WRONGSTATE,"PtAP cannot be reused. Do not call MatFreeIntermediateDataStructures() or use '-mat_freeintermediatedatastructures'");
1972:   }

1974:   MatZeroEntries(C);
1975:   /* 1) get R = Pd^T,Ro = Po^T */
1976:   if (ptap->reuse == MAT_REUSE_MATRIX) {
1977:     MatTranspose(p->A,MAT_REUSE_MATRIX,&ptap->Rd);
1978:     MatTranspose(p->B,MAT_REUSE_MATRIX,&ptap->Ro);
1979:   }

1981:   /* 2) get AP_loc */
1982:   AP_loc = ptap->AP_loc;
1983:   ap = (Mat_SeqAIJ*)AP_loc->data;

1985:   /* 2-1) get P_oth = ptap->P_oth  and P_loc = ptap->P_loc */
1986:   /*-----------------------------------------------------*/
1987:   if (ptap->reuse == MAT_REUSE_MATRIX) {
1988:     /* P_oth and P_loc are obtained in MatPtASymbolic() when reuse == MAT_INITIAL_MATRIX */
1989:     MatGetBrowsOfAoCols_MPIAIJ(A,P,MAT_REUSE_MATRIX,&ptap->startsj_s,&ptap->startsj_r,&ptap->bufa,&ptap->P_oth);
1990:     MatMPIAIJGetLocalMat(P,MAT_REUSE_MATRIX,&ptap->P_loc);
1991:   }

1993:   /* 2-2) compute numeric A_loc*P - dominating part */
1994:   /* ---------------------------------------------- */
1995:   /* get data from symbolic products */
1996:   p_loc = (Mat_SeqAIJ*)(ptap->P_loc)->data;
1997:   if (ptap->P_oth) {
1998:     p_oth = (Mat_SeqAIJ*)(ptap->P_oth)->data;
1999:   }
2000:   apa   = ptap->apa;
2001:   api   = ap->i;
2002:   apj   = ap->j;
2003:   for (i=0; i<am; i++) {
2004:     /* AP[i,:] = A[i,:]*P = Ad*P_loc Ao*P_oth */
2005:     AProw_nonscalable(i,ad,ao,p_loc,p_oth,apa);
2006:     apnz = api[i+1] - api[i];
2007:     for (j=0; j<apnz; j++) {
2008:       col = apj[j+api[i]];
2009:       ap->a[j+ap->i[i]] = apa[col];
2010:       apa[col] = 0.0;
2011:     }
2012:   }

2014:   /* 3) C_loc = Rd*AP_loc, C_oth = Ro*AP_loc */
2015:   ((ptap->C_loc)->ops->matmultnumeric)(ptap->Rd,AP_loc,ptap->C_loc);
2016:   ((ptap->C_oth)->ops->matmultnumeric)(ptap->Ro,AP_loc,ptap->C_oth);
2017:   C_loc = ptap->C_loc;
2018:   C_oth = ptap->C_oth;

2020:   /* add C_loc and Co to to C */
2021:   MatGetOwnershipRange(C,&rstart,&rend);

2023:   /* C_loc -> C */
2024:   cm    = C_loc->rmap->N;
2025:   c_seq = (Mat_SeqAIJ*)C_loc->data;
2026:   cols = c_seq->j;
2027:   vals = c_seq->a;


2030:   /* The (fast) MatSetValues_MPIAIJ_CopyFromCSRFormat function can only be used when C->was_assembled is PETSC_FALSE and */
2031:   /* when there are no off-processor parts.  */
2032:   /* If was_assembled is true, then the statement aj[rowstart_diag+dnz_row] = mat_j[col] - cstart; in MatSetValues_MPIAIJ_CopyFromCSRFormat */
2033:   /* is no longer true. Then the more complex function MatSetValues_MPIAIJ() has to be used, where the column index is looked up from */
2034:   /* a table, and other, more complex stuff has to be done. */
2035:   if (C->assembled) {
2036:     C->was_assembled = PETSC_TRUE;
2037:     C->assembled     = PETSC_FALSE;
2038:   }
2039:   if (C->was_assembled) {
2040:     for (i=0; i<cm; i++) {
2041:       ncols = c_seq->i[i+1] - c_seq->i[i];
2042:       row = rstart + i;
2043:       MatSetValues_MPIAIJ(C,1,&row,ncols,cols,vals,ADD_VALUES);
2044:       cols += ncols; vals += ncols;
2045:     }
2046:   } else {
2047:     MatSetValues_MPIAIJ_CopyFromCSRFormat(C,c_seq->j,c_seq->i,c_seq->a);
2048:   }

2050:   /* Co -> C, off-processor part */
2051:   cm = C_oth->rmap->N;
2052:   c_seq = (Mat_SeqAIJ*)C_oth->data;
2053:   cols = c_seq->j;
2054:   vals = c_seq->a;
2055:   for (i=0; i<cm; i++) {
2056:     ncols = c_seq->i[i+1] - c_seq->i[i];
2057:     row = p->garray[i];
2058:     MatSetValues(C,1,&row,ncols,cols,vals,ADD_VALUES);
2059:     cols += ncols; vals += ncols;
2060:   }

2062:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
2063:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);

2065:   ptap->reuse = MAT_REUSE_MATRIX;

2067:   /* supporting struct ptap consumes almost same amount of memory as C=PtAP, release it if C will not be updated by A and P */
2068:   if (ptap->freestruct) {
2069:     MatFreeIntermediateDataStructures(C);
2070:   }
2071:   return(0);
2072: }