Actual source code: mpiptap.c
petsc-master 2019-06-22
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),¤t_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: PetscCalloc1(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: PetscCalloc1(rootspacesize,&rootspace);
1210: PetscCalloc1(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: PetscCalloc1(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: PetscCalloc2(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: PetscCalloc1(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: PetscCalloc1(rootspacesize,&rootspace);
1486: PetscCalloc1(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: PetscCalloc1(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: PetscCalloc1(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; lidx = 0;
1519: PetscLayoutFindOwnerIndex(P->cmap,p->garray[i],&owner,&lidx);
1520: sendncols = ptap->c_rmti[i+1] - ptap->c_rmti[i];
1521: for (j=0; j<sendncols; j++) {
1522: iremote[nleaves].rank = owner;
1523: iremote[nleaves++].index = c_rmtoffsets[i] + j;
1524: }
1525: }
1526: PetscFree(c_rmtoffsets);
1527: PetscCalloc1(ptap->c_othi[pn],&c_othj);
1529: PetscSFCreate(comm,&ptap->sf);
1530: PetscSFSetGraph(ptap->sf,ptap->c_othi[pn],nleaves,NULL,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER);
1531: PetscSFSetFromOptions(ptap->sf);
1532: /* One to one map */
1533: PetscSFReduceBegin(ptap->sf,MPIU_INT,c_rmtj,c_othj,MPIU_REPLACE);
1534: /* Get remote data */
1535: PetscSFReduceEnd(ptap->sf,MPIU_INT,c_rmtj,c_othj,MPIU_REPLACE);
1536: PetscFree(c_rmtj);
1537: PetscCalloc2(pn,&dnz,pn,&onz);
1538: MatGetOwnershipRangeColumn(P,&pcstart,&pcend);
1540: for (i = 0; i < pn; i++) {
1541: nzi = ptap->c_othi[i+1] - ptap->c_othi[i];
1542: rdj = c_othj + ptap->c_othi[i];
1543: for (j = 0; j < nzi; j++) {
1544: col = rdj[j];
1545: /* diag part */
1546: if (col>=pcstart && col<pcend) {
1547: PetscHSetIAdd(htd[i],col);
1548: } else { /* off diag */
1549: PetscHSetIAdd(hto[i],col);
1550: }
1551: }
1552: PetscHSetIGetSize(htd[i],&htsize);
1553: dnz[i] = htsize;
1554: PetscHSetIDestroy(&htd[i]);
1555: PetscHSetIGetSize(hto[i],&htsize);
1556: onz[i] = htsize;
1557: PetscHSetIDestroy(&hto[i]);
1558: }
1560: PetscFree2(htd,hto);
1561: PetscFree(c_othj);
1563: /* local sizes and preallocation */
1564: MatSetSizes(Cmpi,pn,pn,PETSC_DETERMINE,PETSC_DETERMINE);
1565: MatSetBlockSizes(Cmpi,PetscAbs(P->cmap->bs),PetscAbs(P->cmap->bs));
1566: MatMPIAIJSetPreallocation(Cmpi,0,dnz,0,onz);
1567: PetscFree2(dnz,onz);
1569: /* attach the supporting struct to Cmpi for reuse */
1570: c = (Mat_MPIAIJ*)Cmpi->data;
1571: c->ap = ptap;
1572: ptap->duplicate = Cmpi->ops->duplicate;
1573: ptap->destroy = Cmpi->ops->destroy;
1574: ptap->view = Cmpi->ops->view;
1576: /* Cmpi is not ready for use - assembly will be done by MatPtAPNumeric() */
1577: Cmpi->assembled = PETSC_FALSE;
1578: /* pick an algorithm */
1579: PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MatPtAP","Mat");
1580: alg = 0;
1581: PetscOptionsEList("-matptap_allatonce_via","PtAP allatonce numeric approach","MatPtAP",algTypes,nalg,algTypes[alg],&alg,&flg);
1582: PetscOptionsEnd();
1583: switch (alg) {
1584: case 0:
1585: Cmpi->ops->ptapnumeric = MatPtAPNumeric_MPIAIJ_MPIAIJ_allatonce_merged;
1586: break;
1587: case 1:
1588: Cmpi->ops->ptapnumeric = MatPtAPNumeric_MPIAIJ_MPIAIJ_allatonce;
1589: break;
1590: default:
1591: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG," Unsupported allatonce numerical algorithm \n");
1592: }
1593: Cmpi->ops->destroy = MatDestroy_MPIAIJ_PtAP;
1594: Cmpi->ops->view = MatView_MPIAIJ_PtAP;
1595: Cmpi->ops->freeintermediatedatastructures = MatFreeIntermediateDataStructures_MPIAIJ_AP;
1596: *C = Cmpi;
1597: return(0);
1598: }
1600: PetscErrorCode MatPtAPSymbolic_MPIAIJ_MPIAIJ(Mat A,Mat P,PetscReal fill,Mat *C)
1601: {
1602: PetscErrorCode ierr;
1603: Mat_APMPI *ptap;
1604: Mat_MPIAIJ *a=(Mat_MPIAIJ*)A->data,*p=(Mat_MPIAIJ*)P->data,*c;
1605: MPI_Comm comm;
1606: PetscMPIInt size,rank;
1607: Mat Cmpi;
1608: PetscFreeSpaceList free_space=NULL,current_space=NULL;
1609: PetscInt am=A->rmap->n,pm=P->rmap->n,pN=P->cmap->N,pn=P->cmap->n;
1610: PetscInt *lnk,i,k,pnz,row,nsend;
1611: PetscBT lnkbt;
1612: PetscMPIInt tagi,tagj,*len_si,*len_s,*len_ri,icompleted=0,nrecv;
1613: PetscInt **buf_rj,**buf_ri,**buf_ri_k;
1614: PetscInt len,proc,*dnz,*onz,*owners,nzi,nspacedouble;
1615: PetscInt nrows,*buf_s,*buf_si,*buf_si_i,**nextrow,**nextci;
1616: MPI_Request *swaits,*rwaits;
1617: MPI_Status *sstatus,rstatus;
1618: PetscLayout rowmap;
1619: PetscInt *owners_co,*coi,*coj; /* i and j array of (p->B)^T*A*P - used in the communication */
1620: PetscMPIInt *len_r,*id_r; /* array of length of comm->size, store send/recv matrix values */
1621: PetscInt *api,*apj,*Jptr,apnz,*prmap=p->garray,con,j,ap_rmax=0,Crmax,*aj,*ai,*pi;
1622: Mat_SeqAIJ *p_loc,*p_oth=NULL,*ad=(Mat_SeqAIJ*)(a->A)->data,*ao=NULL,*c_loc,*c_oth;
1623: PetscScalar *apv;
1624: PetscTable ta;
1625: MatType mtype;
1626: const char *prefix;
1627: #if defined(PETSC_USE_INFO)
1628: PetscReal apfill;
1629: #endif
1632: PetscObjectGetComm((PetscObject)A,&comm);
1633: MPI_Comm_size(comm,&size);
1634: MPI_Comm_rank(comm,&rank);
1636: if (size > 1) ao = (Mat_SeqAIJ*)(a->B)->data;
1638: /* create symbolic parallel matrix Cmpi */
1639: MatCreate(comm,&Cmpi);
1640: MatGetType(A,&mtype);
1641: MatSetType(Cmpi,mtype);
1643: /* Do dense axpy in MatPtAPNumeric_MPIAIJ_MPIAIJ() */
1644: Cmpi->ops->ptapnumeric = MatPtAPNumeric_MPIAIJ_MPIAIJ;
1646: /* create struct Mat_APMPI and attached it to C later */
1647: PetscNew(&ptap);
1648: ptap->reuse = MAT_INITIAL_MATRIX;
1649: ptap->algType = 1;
1651: /* get P_oth by taking rows of P (= non-zero cols of local A) from other processors */
1652: MatGetBrowsOfAoCols_MPIAIJ(A,P,MAT_INITIAL_MATRIX,&ptap->startsj_s,&ptap->startsj_r,&ptap->bufa,&ptap->P_oth);
1653: /* get P_loc by taking all local rows of P */
1654: MatMPIAIJGetLocalMat(P,MAT_INITIAL_MATRIX,&ptap->P_loc);
1656: /* (0) compute Rd = Pd^T, Ro = Po^T */
1657: /* --------------------------------- */
1658: MatTranspose(p->A,MAT_INITIAL_MATRIX,&ptap->Rd);
1659: MatTranspose(p->B,MAT_INITIAL_MATRIX,&ptap->Ro);
1661: /* (1) compute symbolic AP = A_loc*P = Ad*P_loc + Ao*P_oth (api,apj) */
1662: /* ----------------------------------------------------------------- */
1663: p_loc = (Mat_SeqAIJ*)(ptap->P_loc)->data;
1664: if (ptap->P_oth) p_oth = (Mat_SeqAIJ*)(ptap->P_oth)->data;
1666: /* create and initialize a linked list */
1667: PetscTableCreate(pn,pN,&ta); /* for compute AP_loc and Cmpi */
1668: MatRowMergeMax_SeqAIJ(p_loc,ptap->P_loc->rmap->N,ta);
1669: MatRowMergeMax_SeqAIJ(p_oth,ptap->P_oth->rmap->N,ta);
1670: PetscTableGetCount(ta,&Crmax); /* Crmax = nnz(sum of Prows) */
1671: /* printf("[%d] est %d, Crmax %d; pN %d\n",rank,5*(p_loc->rmax+p_oth->rmax + (PetscInt)(1.e-2*pN)),Crmax,pN); */
1673: PetscLLCondensedCreate(Crmax,pN,&lnk,&lnkbt);
1675: /* Initial FreeSpace size is fill*(nnz(A) + nnz(P)) */
1676: if (ao) {
1677: PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,PetscIntSumTruncate(ad->i[am],PetscIntSumTruncate(ao->i[am],p_loc->i[pm]))),&free_space);
1678: } else {
1679: PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,PetscIntSumTruncate(ad->i[am],p_loc->i[pm])),&free_space);
1680: }
1681: current_space = free_space;
1682: nspacedouble = 0;
1684: PetscMalloc1(am+1,&api);
1685: api[0] = 0;
1686: for (i=0; i<am; i++) {
1687: /* diagonal portion: Ad[i,:]*P */
1688: ai = ad->i; pi = p_loc->i;
1689: nzi = ai[i+1] - ai[i];
1690: aj = ad->j + ai[i];
1691: for (j=0; j<nzi; j++) {
1692: row = aj[j];
1693: pnz = pi[row+1] - pi[row];
1694: Jptr = p_loc->j + pi[row];
1695: /* add non-zero cols of P into the sorted linked list lnk */
1696: PetscLLCondensedAddSorted(pnz,Jptr,lnk,lnkbt);
1697: }
1698: /* off-diagonal portion: Ao[i,:]*P */
1699: if (ao) {
1700: ai = ao->i; pi = p_oth->i;
1701: nzi = ai[i+1] - ai[i];
1702: aj = ao->j + ai[i];
1703: for (j=0; j<nzi; j++) {
1704: row = aj[j];
1705: pnz = pi[row+1] - pi[row];
1706: Jptr = p_oth->j + pi[row];
1707: PetscLLCondensedAddSorted(pnz,Jptr,lnk,lnkbt);
1708: }
1709: }
1710: apnz = lnk[0];
1711: api[i+1] = api[i] + apnz;
1712: if (ap_rmax < apnz) ap_rmax = apnz;
1714: /* if free space is not available, double the total space in the list */
1715: if (current_space->local_remaining<apnz) {
1716: PetscFreeSpaceGet(PetscIntSumTruncate(apnz,current_space->total_array_size),¤t_space);
1717: nspacedouble++;
1718: }
1720: /* Copy data into free space, then initialize lnk */
1721: PetscLLCondensedClean(pN,apnz,current_space->array,lnk,lnkbt);
1723: current_space->array += apnz;
1724: current_space->local_used += apnz;
1725: current_space->local_remaining -= apnz;
1726: }
1727: /* Allocate space for apj and apv, initialize apj, and */
1728: /* destroy list of free space and other temporary array(s) */
1729: PetscMalloc2(api[am],&apj,api[am],&apv);
1730: PetscFreeSpaceContiguous(&free_space,apj);
1731: PetscLLDestroy(lnk,lnkbt);
1733: /* Create AP_loc for reuse */
1734: MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,am,pN,api,apj,apv,&ptap->AP_loc);
1736: #if defined(PETSC_USE_INFO)
1737: if (ao) {
1738: apfill = (PetscReal)api[am]/(ad->i[am]+ao->i[am]+p_loc->i[pm]+1);
1739: } else {
1740: apfill = (PetscReal)api[am]/(ad->i[am]+p_loc->i[pm]+1);
1741: }
1742: ptap->AP_loc->info.mallocs = nspacedouble;
1743: ptap->AP_loc->info.fill_ratio_given = fill;
1744: ptap->AP_loc->info.fill_ratio_needed = apfill;
1746: if (api[am]) {
1747: PetscInfo3(ptap->AP_loc,"Nonscalable algorithm, AP_loc reallocs %D; Fill ratio: given %g needed %g.\n",nspacedouble,(double)fill,(double)apfill);
1748: PetscInfo1(ptap->AP_loc,"Use MatPtAP(A,B,MatReuse,%g,&C) for best AP_loc performance.;\n",(double)apfill);
1749: } else {
1750: PetscInfo(ptap->AP_loc,"Nonscalable algorithm, AP_loc is empty \n");
1751: }
1752: #endif
1754: /* (2-1) compute symbolic Co = Ro*AP_loc */
1755: /* ------------------------------------ */
1756: MatGetOptionsPrefix(A,&prefix);
1757: MatSetOptionsPrefix(ptap->Ro,prefix);
1758: MatAppendOptionsPrefix(ptap->Ro,"inner_offdiag_");
1759: MatMatMultSymbolic_SeqAIJ_SeqAIJ(ptap->Ro,ptap->AP_loc,fill,&ptap->C_oth);
1761: /* (3) send coj of C_oth to other processors */
1762: /* ------------------------------------------ */
1763: /* determine row ownership */
1764: PetscLayoutCreate(comm,&rowmap);
1765: rowmap->n = pn;
1766: rowmap->bs = 1;
1767: PetscLayoutSetUp(rowmap);
1768: owners = rowmap->range;
1770: /* determine the number of messages to send, their lengths */
1771: PetscMalloc4(size,&len_s,size,&len_si,size,&sstatus,size+2,&owners_co);
1772: PetscArrayzero(len_s,size);
1773: PetscArrayzero(len_si,size);
1775: c_oth = (Mat_SeqAIJ*)ptap->C_oth->data;
1776: coi = c_oth->i; coj = c_oth->j;
1777: con = ptap->C_oth->rmap->n;
1778: proc = 0;
1779: for (i=0; i<con; i++) {
1780: while (prmap[i] >= owners[proc+1]) proc++;
1781: len_si[proc]++; /* num of rows in Co(=Pt*AP) to be sent to [proc] */
1782: len_s[proc] += coi[i+1] - coi[i]; /* num of nonzeros in Co to be sent to [proc] */
1783: }
1785: len = 0; /* max length of buf_si[], see (4) */
1786: owners_co[0] = 0;
1787: nsend = 0;
1788: for (proc=0; proc<size; proc++) {
1789: owners_co[proc+1] = owners_co[proc] + len_si[proc];
1790: if (len_s[proc]) {
1791: nsend++;
1792: len_si[proc] = 2*(len_si[proc] + 1); /* length of buf_si to be sent to [proc] */
1793: len += len_si[proc];
1794: }
1795: }
1797: /* determine the number and length of messages to receive for coi and coj */
1798: PetscGatherNumberOfMessages(comm,NULL,len_s,&nrecv);
1799: PetscGatherMessageLengths2(comm,nsend,nrecv,len_s,len_si,&id_r,&len_r,&len_ri);
1801: /* post the Irecv and Isend of coj */
1802: PetscCommGetNewTag(comm,&tagj);
1803: PetscPostIrecvInt(comm,tagj,nrecv,id_r,len_r,&buf_rj,&rwaits);
1804: PetscMalloc1(nsend+1,&swaits);
1805: for (proc=0, k=0; proc<size; proc++) {
1806: if (!len_s[proc]) continue;
1807: i = owners_co[proc];
1808: MPI_Isend(coj+coi[i],len_s[proc],MPIU_INT,proc,tagj,comm,swaits+k);
1809: k++;
1810: }
1812: /* (2-2) compute symbolic C_loc = Rd*AP_loc */
1813: /* ---------------------------------------- */
1814: MatSetOptionsPrefix(ptap->Rd,prefix);
1815: MatAppendOptionsPrefix(ptap->Rd,"inner_diag_");
1816: MatMatMultSymbolic_SeqAIJ_SeqAIJ(ptap->Rd,ptap->AP_loc,fill,&ptap->C_loc);
1817: c_loc = (Mat_SeqAIJ*)ptap->C_loc->data;
1819: /* receives coj are complete */
1820: for (i=0; i<nrecv; i++) {
1821: MPI_Waitany(nrecv,rwaits,&icompleted,&rstatus);
1822: }
1823: PetscFree(rwaits);
1824: if (nsend) {MPI_Waitall(nsend,swaits,sstatus);}
1826: /* add received column indices into ta to update Crmax */
1827: for (k=0; k<nrecv; k++) {/* k-th received message */
1828: Jptr = buf_rj[k];
1829: for (j=0; j<len_r[k]; j++) {
1830: PetscTableAdd(ta,*(Jptr+j)+1,1,INSERT_VALUES);
1831: }
1832: }
1833: PetscTableGetCount(ta,&Crmax);
1834: PetscTableDestroy(&ta);
1836: /* (4) send and recv coi */
1837: /*-----------------------*/
1838: PetscCommGetNewTag(comm,&tagi);
1839: PetscPostIrecvInt(comm,tagi,nrecv,id_r,len_ri,&buf_ri,&rwaits);
1840: PetscMalloc1(len+1,&buf_s);
1841: buf_si = buf_s; /* points to the beginning of k-th msg to be sent */
1842: for (proc=0,k=0; proc<size; proc++) {
1843: if (!len_s[proc]) continue;
1844: /* form outgoing message for i-structure:
1845: buf_si[0]: nrows to be sent
1846: [1:nrows]: row index (global)
1847: [nrows+1:2*nrows+1]: i-structure index
1848: */
1849: /*-------------------------------------------*/
1850: nrows = len_si[proc]/2 - 1; /* num of rows in Co to be sent to [proc] */
1851: buf_si_i = buf_si + nrows+1;
1852: buf_si[0] = nrows;
1853: buf_si_i[0] = 0;
1854: nrows = 0;
1855: for (i=owners_co[proc]; i<owners_co[proc+1]; i++) {
1856: nzi = coi[i+1] - coi[i];
1857: buf_si_i[nrows+1] = buf_si_i[nrows] + nzi; /* i-structure */
1858: buf_si[nrows+1] = prmap[i] -owners[proc]; /* local row index */
1859: nrows++;
1860: }
1861: MPI_Isend(buf_si,len_si[proc],MPIU_INT,proc,tagi,comm,swaits+k);
1862: k++;
1863: buf_si += len_si[proc];
1864: }
1865: for (i=0; i<nrecv; i++) {
1866: MPI_Waitany(nrecv,rwaits,&icompleted,&rstatus);
1867: }
1868: PetscFree(rwaits);
1869: if (nsend) {MPI_Waitall(nsend,swaits,sstatus);}
1871: PetscFree4(len_s,len_si,sstatus,owners_co);
1872: PetscFree(len_ri);
1873: PetscFree(swaits);
1874: PetscFree(buf_s);
1876: /* (5) compute the local portion of Cmpi */
1877: /* ------------------------------------------ */
1878: /* set initial free space to be Crmax, sufficient for holding nozeros in each row of Cmpi */
1879: PetscFreeSpaceGet(Crmax,&free_space);
1880: current_space = free_space;
1882: PetscMalloc3(nrecv,&buf_ri_k,nrecv,&nextrow,nrecv,&nextci);
1883: for (k=0; k<nrecv; k++) {
1884: buf_ri_k[k] = buf_ri[k]; /* beginning of k-th recved i-structure */
1885: nrows = *buf_ri_k[k];
1886: nextrow[k] = buf_ri_k[k] + 1; /* next row number of k-th recved i-structure */
1887: nextci[k] = buf_ri_k[k] + (nrows + 1); /* poins to the next i-structure of k-th recved i-structure */
1888: }
1890: MatPreallocateInitialize(comm,pn,pn,dnz,onz);
1891: PetscLLCondensedCreate(Crmax,pN,&lnk,&lnkbt);
1892: for (i=0; i<pn; i++) {
1893: /* add C_loc into Cmpi */
1894: nzi = c_loc->i[i+1] - c_loc->i[i];
1895: Jptr = c_loc->j + c_loc->i[i];
1896: PetscLLCondensedAddSorted(nzi,Jptr,lnk,lnkbt);
1898: /* add received col data into lnk */
1899: for (k=0; k<nrecv; k++) { /* k-th received message */
1900: if (i == *nextrow[k]) { /* i-th row */
1901: nzi = *(nextci[k]+1) - *nextci[k];
1902: Jptr = buf_rj[k] + *nextci[k];
1903: PetscLLCondensedAddSorted(nzi,Jptr,lnk,lnkbt);
1904: nextrow[k]++; nextci[k]++;
1905: }
1906: }
1907: nzi = lnk[0];
1909: /* copy data into free space, then initialize lnk */
1910: PetscLLCondensedClean(pN,nzi,current_space->array,lnk,lnkbt);
1911: MatPreallocateSet(i+owners[rank],nzi,current_space->array,dnz,onz);
1912: }
1913: PetscFree3(buf_ri_k,nextrow,nextci);
1914: PetscLLDestroy(lnk,lnkbt);
1915: PetscFreeSpaceDestroy(free_space);
1917: /* local sizes and preallocation */
1918: MatSetSizes(Cmpi,pn,pn,PETSC_DETERMINE,PETSC_DETERMINE);
1919: if (P->cmap->bs > 0) {
1920: PetscLayoutSetBlockSize(Cmpi->rmap,P->cmap->bs);
1921: PetscLayoutSetBlockSize(Cmpi->cmap,P->cmap->bs);
1922: }
1923: MatMPIAIJSetPreallocation(Cmpi,0,dnz,0,onz);
1924: MatPreallocateFinalize(dnz,onz);
1926: /* members in merge */
1927: PetscFree(id_r);
1928: PetscFree(len_r);
1929: PetscFree(buf_ri[0]);
1930: PetscFree(buf_ri);
1931: PetscFree(buf_rj[0]);
1932: PetscFree(buf_rj);
1933: PetscLayoutDestroy(&rowmap);
1935: /* attach the supporting struct to Cmpi for reuse */
1936: c = (Mat_MPIAIJ*)Cmpi->data;
1937: c->ap = ptap;
1938: ptap->duplicate = Cmpi->ops->duplicate;
1939: ptap->destroy = Cmpi->ops->destroy;
1940: ptap->view = Cmpi->ops->view;
1941: PetscCalloc1(pN,&ptap->apa);
1943: /* Cmpi is not ready for use - assembly will be done by MatPtAPNumeric() */
1944: Cmpi->assembled = PETSC_FALSE;
1945: Cmpi->ops->destroy = MatDestroy_MPIAIJ_PtAP;
1946: Cmpi->ops->view = MatView_MPIAIJ_PtAP;
1947: Cmpi->ops->freeintermediatedatastructures = MatFreeIntermediateDataStructures_MPIAIJ_AP;
1948: *C = Cmpi;
1949: return(0);
1950: }
1952: PetscErrorCode MatPtAPNumeric_MPIAIJ_MPIAIJ(Mat A,Mat P,Mat C)
1953: {
1954: PetscErrorCode ierr;
1955: Mat_MPIAIJ *a=(Mat_MPIAIJ*)A->data,*p=(Mat_MPIAIJ*)P->data,*c=(Mat_MPIAIJ*)C->data;
1956: Mat_SeqAIJ *ad=(Mat_SeqAIJ*)(a->A)->data,*ao=(Mat_SeqAIJ*)(a->B)->data;
1957: Mat_SeqAIJ *ap,*p_loc,*p_oth=NULL,*c_seq;
1958: Mat_APMPI *ptap = c->ap;
1959: Mat AP_loc,C_loc,C_oth;
1960: PetscInt i,rstart,rend,cm,ncols,row;
1961: PetscInt *api,*apj,am = A->rmap->n,j,col,apnz;
1962: PetscScalar *apa;
1963: const PetscInt *cols;
1964: const PetscScalar *vals;
1967: if (!ptap->AP_loc) {
1968: MPI_Comm comm;
1969: PetscObjectGetComm((PetscObject)C,&comm);
1970: SETERRQ(comm,PETSC_ERR_ARG_WRONGSTATE,"PtAP cannot be reused. Do not call MatFreeIntermediateDataStructures() or use '-mat_freeintermediatedatastructures'");
1971: }
1973: MatZeroEntries(C);
1974: /* 1) get R = Pd^T,Ro = Po^T */
1975: if (ptap->reuse == MAT_REUSE_MATRIX) {
1976: MatTranspose(p->A,MAT_REUSE_MATRIX,&ptap->Rd);
1977: MatTranspose(p->B,MAT_REUSE_MATRIX,&ptap->Ro);
1978: }
1980: /* 2) get AP_loc */
1981: AP_loc = ptap->AP_loc;
1982: ap = (Mat_SeqAIJ*)AP_loc->data;
1984: /* 2-1) get P_oth = ptap->P_oth and P_loc = ptap->P_loc */
1985: /*-----------------------------------------------------*/
1986: if (ptap->reuse == MAT_REUSE_MATRIX) {
1987: /* P_oth and P_loc are obtained in MatPtASymbolic() when reuse == MAT_INITIAL_MATRIX */
1988: MatGetBrowsOfAoCols_MPIAIJ(A,P,MAT_REUSE_MATRIX,&ptap->startsj_s,&ptap->startsj_r,&ptap->bufa,&ptap->P_oth);
1989: MatMPIAIJGetLocalMat(P,MAT_REUSE_MATRIX,&ptap->P_loc);
1990: }
1992: /* 2-2) compute numeric A_loc*P - dominating part */
1993: /* ---------------------------------------------- */
1994: /* get data from symbolic products */
1995: p_loc = (Mat_SeqAIJ*)(ptap->P_loc)->data;
1996: if (ptap->P_oth) {
1997: p_oth = (Mat_SeqAIJ*)(ptap->P_oth)->data;
1998: }
1999: apa = ptap->apa;
2000: api = ap->i;
2001: apj = ap->j;
2002: for (i=0; i<am; i++) {
2003: /* AP[i,:] = A[i,:]*P = Ad*P_loc Ao*P_oth */
2004: AProw_nonscalable(i,ad,ao,p_loc,p_oth,apa);
2005: apnz = api[i+1] - api[i];
2006: for (j=0; j<apnz; j++) {
2007: col = apj[j+api[i]];
2008: ap->a[j+ap->i[i]] = apa[col];
2009: apa[col] = 0.0;
2010: }
2011: }
2013: /* 3) C_loc = Rd*AP_loc, C_oth = Ro*AP_loc */
2014: ((ptap->C_loc)->ops->matmultnumeric)(ptap->Rd,AP_loc,ptap->C_loc);
2015: ((ptap->C_oth)->ops->matmultnumeric)(ptap->Ro,AP_loc,ptap->C_oth);
2016: C_loc = ptap->C_loc;
2017: C_oth = ptap->C_oth;
2019: /* add C_loc and Co to to C */
2020: MatGetOwnershipRange(C,&rstart,&rend);
2022: /* C_loc -> C */
2023: cm = C_loc->rmap->N;
2024: c_seq = (Mat_SeqAIJ*)C_loc->data;
2025: cols = c_seq->j;
2026: vals = c_seq->a;
2029: /* The (fast) MatSetValues_MPIAIJ_CopyFromCSRFormat function can only be used when C->was_assembled is PETSC_FALSE and */
2030: /* when there are no off-processor parts. */
2031: /* If was_assembled is true, then the statement aj[rowstart_diag+dnz_row] = mat_j[col] - cstart; in MatSetValues_MPIAIJ_CopyFromCSRFormat */
2032: /* is no longer true. Then the more complex function MatSetValues_MPIAIJ() has to be used, where the column index is looked up from */
2033: /* a table, and other, more complex stuff has to be done. */
2034: if (C->assembled) {
2035: C->was_assembled = PETSC_TRUE;
2036: C->assembled = PETSC_FALSE;
2037: }
2038: if (C->was_assembled) {
2039: for (i=0; i<cm; i++) {
2040: ncols = c_seq->i[i+1] - c_seq->i[i];
2041: row = rstart + i;
2042: MatSetValues_MPIAIJ(C,1,&row,ncols,cols,vals,ADD_VALUES);
2043: cols += ncols; vals += ncols;
2044: }
2045: } else {
2046: MatSetValues_MPIAIJ_CopyFromCSRFormat(C,c_seq->j,c_seq->i,c_seq->a);
2047: }
2049: /* Co -> C, off-processor part */
2050: cm = C_oth->rmap->N;
2051: c_seq = (Mat_SeqAIJ*)C_oth->data;
2052: cols = c_seq->j;
2053: vals = c_seq->a;
2054: for (i=0; i<cm; i++) {
2055: ncols = c_seq->i[i+1] - c_seq->i[i];
2056: row = p->garray[i];
2057: MatSetValues(C,1,&row,ncols,cols,vals,ADD_VALUES);
2058: cols += ncols; vals += ncols;
2059: }
2061: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
2062: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
2064: ptap->reuse = MAT_REUSE_MATRIX;
2066: /* supporting struct ptap consumes almost same amount of memory as C=PtAP, release it if C will not be updated by A and P */
2067: if (ptap->freestruct) {
2068: MatFreeIntermediateDataStructures(C);
2069: }
2070: return(0);
2071: }