Actual source code: mpiptap.c
petsc-master 2020-02-15
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,const PetscInt *map,PetscInt dof,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,offset;
707: PetscErrorCode ierr;
710: pcstart = P->cmap->rstart;
711: pcstart *= dof;
712: pcend = P->cmap->rend;
713: pcend *= dof;
714: /* diagonal portion: Ad[i,:]*P */
715: ai = ad->i;
716: nzi = ai[i+1] - ai[i];
717: aj = ad->j + ai[i];
718: for (j=0; j<nzi; j++) {
719: row = aj[j];
720: offset = row%dof;
721: row /= dof;
722: nzpi = pd->i[row+1] - pd->i[row];
723: pj = pd->j + pd->i[row];
724: for (k=0; k<nzpi; k++) {
725: PetscHSetIAdd(dht,pj[k]*dof+offset+pcstart);
726: }
727: }
728: /* off diag P */
729: for (j=0; j<nzi; j++) {
730: row = aj[j];
731: offset = row%dof;
732: row /= dof;
733: nzpi = po->i[row+1] - po->i[row];
734: pj = po->j + po->i[row];
735: for (k=0; k<nzpi; k++) {
736: PetscHSetIAdd(oht,p->garray[pj[k]]*dof+offset);
737: }
738: }
740: /* off diagonal part: Ao[i, :]*P_oth */
741: if (ao) {
742: ai = ao->i;
743: pi = p_oth->i;
744: nzi = ai[i+1] - ai[i];
745: aj = ao->j + ai[i];
746: for (j=0; j<nzi; j++) {
747: row = aj[j];
748: offset = a->garray[row]%dof;
749: row = map[row];
750: pnz = pi[row+1] - pi[row];
751: p_othcols = p_oth->j + pi[row];
752: for (col=0; col<pnz; col++) {
753: column = p_othcols[col] * dof + offset;
754: if (column>=pcstart && column<pcend) {
755: PetscHSetIAdd(dht,column);
756: } else {
757: PetscHSetIAdd(oht,column);
758: }
759: }
760: }
761: } /* end if (ao) */
762: return(0);
763: }
765: PETSC_STATIC_INLINE PetscErrorCode MatPtAPNumericComputeOneRowOfAP_private(Mat A,Mat P,Mat P_oth,const PetscInt *map,PetscInt dof,PetscInt i,PetscHMapIV hmap)
766: {
767: Mat_MPIAIJ *a=(Mat_MPIAIJ*)A->data,*p=(Mat_MPIAIJ*)P->data;
768: 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;
769: PetscInt *ai,nzi,j,*aj,row,col,*pi,pnz,*p_othcols,pcstart,*pj,k,nzpi,offset;
770: PetscScalar ra,*aa,*pa;
771: PetscErrorCode ierr;
774: pcstart = P->cmap->rstart;
775: pcstart *= dof;
777: /* diagonal portion: Ad[i,:]*P */
778: ai = ad->i;
779: nzi = ai[i+1] - ai[i];
780: aj = ad->j + ai[i];
781: aa = ad->a + ai[i];
782: for (j=0; j<nzi; j++) {
783: ra = aa[j];
784: row = aj[j];
785: offset = row%dof;
786: row /= dof;
787: nzpi = pd->i[row+1] - pd->i[row];
788: pj = pd->j + pd->i[row];
789: pa = pd->a + pd->i[row];
790: for (k=0; k<nzpi; k++) {
791: PetscHMapIVAddValue(hmap,pj[k]*dof+offset+pcstart,ra*pa[k]);
792: }
793: PetscLogFlops(2.0*nzpi);
794: }
795: for (j=0; j<nzi; j++) {
796: ra = aa[j];
797: row = aj[j];
798: offset = row%dof;
799: row /= dof;
800: nzpi = po->i[row+1] - po->i[row];
801: pj = po->j + po->i[row];
802: pa = po->a + po->i[row];
803: for (k=0; k<nzpi; k++) {
804: PetscHMapIVAddValue(hmap,p->garray[pj[k]]*dof+offset,ra*pa[k]);
805: }
806: PetscLogFlops(2.0*nzpi);
807: }
809: /* off diagonal part: Ao[i, :]*P_oth */
810: if (ao) {
811: ai = ao->i;
812: pi = p_oth->i;
813: nzi = ai[i+1] - ai[i];
814: aj = ao->j + ai[i];
815: aa = ao->a + ai[i];
816: for (j=0; j<nzi; j++) {
817: row = aj[j];
818: offset = a->garray[row]%dof;
819: row = map[row];
820: ra = aa[j];
821: pnz = pi[row+1] - pi[row];
822: p_othcols = p_oth->j + pi[row];
823: pa = p_oth->a + pi[row];
824: for (col=0; col<pnz; col++) {
825: PetscHMapIVAddValue(hmap,p_othcols[col]*dof+offset,ra*pa[col]);
826: }
827: PetscLogFlops(2.0*pnz);
828: }
829: } /* end if (ao) */
831: return(0);
832: }
834: PetscErrorCode MatGetBrowsOfAcols_MPIXAIJ(Mat,Mat,PetscInt dof,MatReuse,Mat*);
836: PetscErrorCode MatPtAPNumeric_MPIAIJ_MPIXAIJ_allatonce(Mat A,Mat P,PetscInt dof,Mat C)
837: {
838: PetscErrorCode ierr;
839: Mat_MPIAIJ *p=(Mat_MPIAIJ*)P->data,*c=(Mat_MPIAIJ*)C->data;
840: Mat_SeqAIJ *cd,*co,*po=(Mat_SeqAIJ*)p->B->data,*pd=(Mat_SeqAIJ*)p->A->data;
841: Mat_APMPI *ptap = c->ap;
842: PetscHMapIV hmap;
843: 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;
844: PetscScalar *c_rmta,*c_otha,*poa,*pda,*apvalues,*apvaluestmp,*c_rmtaa;
845: PetscInt offset,ii,pocol;
846: const PetscInt *mappingindices;
847: IS map;
848: MPI_Comm comm;
851: PetscObjectGetComm((PetscObject)A,&comm);
852: if (!ptap->P_oth) SETERRQ(comm,PETSC_ERR_ARG_WRONGSTATE,"PtAP cannot be reused. Do not call MatFreeIntermediateDataStructures() or use '-mat_freeintermediatedatastructures'");
854: MatZeroEntries(C);
856: /* Get P_oth = ptap->P_oth and P_loc = ptap->P_loc */
857: /*-----------------------------------------------------*/
858: if (ptap->reuse == MAT_REUSE_MATRIX) {
859: /* P_oth and P_loc are obtained in MatPtASymbolic() when reuse == MAT_INITIAL_MATRIX */
860: MatGetBrowsOfAcols_MPIXAIJ(A,P,dof,MAT_REUSE_MATRIX,&ptap->P_oth);
861: }
862: PetscObjectQuery((PetscObject)ptap->P_oth,"aoffdiagtopothmapping",(PetscObject*)&map);
864: MatGetLocalSize(p->B,NULL,&pon);
865: pon *= dof;
866: PetscCalloc2(ptap->c_rmti[pon],&c_rmtj,ptap->c_rmti[pon],&c_rmta);
867: MatGetLocalSize(A,&am,NULL);
868: cmaxr = 0;
869: for (i=0; i<pon; i++) {
870: cmaxr = PetscMax(cmaxr,ptap->c_rmti[i+1]-ptap->c_rmti[i]);
871: }
872: PetscCalloc4(cmaxr,&apindices,cmaxr,&apvalues,cmaxr,&apvaluestmp,pon,&c_rmtc);
873: PetscHMapIVCreate(&hmap);
874: PetscHMapIVResize(hmap,cmaxr);
875: ISGetIndices(map,&mappingindices);
876: for (i=0; i<am && pon; i++) {
877: PetscHMapIVClear(hmap);
878: offset = i%dof;
879: ii = i/dof;
880: nzi = po->i[ii+1] - po->i[ii];
881: if (!nzi) continue;
882: MatPtAPNumericComputeOneRowOfAP_private(A,P,ptap->P_oth,mappingindices,dof,i,hmap);
883: voff = 0;
884: PetscHMapIVGetPairs(hmap,&voff,apindices,apvalues);
885: if (!voff) continue;
887: /* Form C(ii, :) */
888: poj = po->j + po->i[ii];
889: poa = po->a + po->i[ii];
890: for (j=0; j<nzi; j++) {
891: pocol = poj[j]*dof+offset;
892: c_rmtjj = c_rmtj + ptap->c_rmti[pocol];
893: c_rmtaa = c_rmta + ptap->c_rmti[pocol];
894: for (jj=0; jj<voff; jj++) {
895: apvaluestmp[jj] = apvalues[jj]*poa[j];
896: /*If the row is empty */
897: if (!c_rmtc[pocol]) {
898: c_rmtjj[jj] = apindices[jj];
899: c_rmtaa[jj] = apvaluestmp[jj];
900: c_rmtc[pocol]++;
901: } else {
902: PetscFindInt(apindices[jj],c_rmtc[pocol],c_rmtjj,&loc);
903: if (loc>=0){ /* hit */
904: c_rmtaa[loc] += apvaluestmp[jj];
905: PetscLogFlops(1.0);
906: } else { /* new element */
907: loc = -(loc+1);
908: /* Move data backward */
909: for (kk=c_rmtc[pocol]; kk>loc; kk--) {
910: c_rmtjj[kk] = c_rmtjj[kk-1];
911: c_rmtaa[kk] = c_rmtaa[kk-1];
912: }/* End kk */
913: c_rmtjj[loc] = apindices[jj];
914: c_rmtaa[loc] = apvaluestmp[jj];
915: c_rmtc[pocol]++;
916: }
917: }
918: PetscLogFlops(voff);
919: } /* End jj */
920: } /* End j */
921: } /* End i */
923: PetscFree4(apindices,apvalues,apvaluestmp,c_rmtc);
924: PetscHMapIVDestroy(&hmap);
926: MatGetLocalSize(P,NULL,&pn);
927: pn *= dof;
928: PetscCalloc2(ptap->c_othi[pn],&c_othj,ptap->c_othi[pn],&c_otha);
930: PetscSFReduceBegin(ptap->sf,MPIU_INT,c_rmtj,c_othj,MPIU_REPLACE);
931: PetscSFReduceBegin(ptap->sf,MPIU_SCALAR,c_rmta,c_otha,MPIU_REPLACE);
932: MatGetOwnershipRangeColumn(P,&pcstart,&pcend);
933: pcstart = pcstart*dof;
934: pcend = pcend*dof;
935: cd = (Mat_SeqAIJ*)(c->A)->data;
936: co = (Mat_SeqAIJ*)(c->B)->data;
938: cmaxr = 0;
939: for (i=0; i<pn; i++) {
940: cmaxr = PetscMax(cmaxr,(cd->i[i+1]-cd->i[i])+(co->i[i+1]-co->i[i]));
941: }
942: PetscCalloc5(cmaxr,&apindices,cmaxr,&apvalues,cmaxr,&apvaluestmp,pn,&dcc,pn,&occ);
943: PetscHMapIVCreate(&hmap);
944: PetscHMapIVResize(hmap,cmaxr);
945: for (i=0; i<am && pn; i++) {
946: PetscHMapIVClear(hmap);
947: offset = i%dof;
948: ii = i/dof;
949: nzi = pd->i[ii+1] - pd->i[ii];
950: if (!nzi) continue;
951: MatPtAPNumericComputeOneRowOfAP_private(A,P,ptap->P_oth,mappingindices,dof,i,hmap);
952: voff = 0;
953: PetscHMapIVGetPairs(hmap,&voff,apindices,apvalues);
954: if (!voff) continue;
955: /* Form C(ii, :) */
956: pdj = pd->j + pd->i[ii];
957: pda = pd->a + pd->i[ii];
958: for (j=0; j<nzi; j++) {
959: row = pcstart + pdj[j] * dof + offset;
960: for (jj=0; jj<voff; jj++) {
961: apvaluestmp[jj] = apvalues[jj]*pda[j];
962: }
963: PetscLogFlops(voff);
964: MatSetValues(C,1,&row,voff,apindices,apvaluestmp,ADD_VALUES);
965: }
966: }
967: ISRestoreIndices(map,&mappingindices);
968: MatGetOwnershipRangeColumn(C,&ccstart,&ccend);
969: PetscFree5(apindices,apvalues,apvaluestmp,dcc,occ);
970: PetscHMapIVDestroy(&hmap);
971: PetscSFReduceEnd(ptap->sf,MPIU_INT,c_rmtj,c_othj,MPIU_REPLACE);
972: PetscSFReduceEnd(ptap->sf,MPIU_SCALAR,c_rmta,c_otha,MPIU_REPLACE);
973: PetscFree2(c_rmtj,c_rmta);
975: /* Add contributions from remote */
976: for (i = 0; i < pn; i++) {
977: row = i + pcstart;
978: 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);
979: }
980: PetscFree2(c_othj,c_otha);
982: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
983: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
985: ptap->reuse = MAT_REUSE_MATRIX;
987: /* supporting struct ptap consumes almost same amount of memory as C=PtAP, release it if C will not be updated by A and P */
988: if (ptap->freestruct) {
989: MatFreeIntermediateDataStructures(C);
990: }
991: return(0);
992: }
994: PetscErrorCode MatPtAPNumeric_MPIAIJ_MPIAIJ_allatonce(Mat A,Mat P,Mat C)
995: {
996: PetscErrorCode ierr;
1000: MatPtAPNumeric_MPIAIJ_MPIXAIJ_allatonce(A,P,1,C);
1001: return(0);
1002: }
1004: PetscErrorCode MatPtAPNumeric_MPIAIJ_MPIXAIJ_allatonce_merged(Mat A,Mat P,PetscInt dof,Mat C)
1005: {
1006: PetscErrorCode ierr;
1007: Mat_MPIAIJ *p=(Mat_MPIAIJ*)P->data,*c=(Mat_MPIAIJ*)C->data;
1008: Mat_SeqAIJ *cd,*co,*po=(Mat_SeqAIJ*)p->B->data,*pd=(Mat_SeqAIJ*)p->A->data;
1009: Mat_APMPI *ptap = c->ap;
1010: PetscHMapIV hmap;
1011: 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;
1012: PetscScalar *c_rmta,*c_otha,*poa,*pda,*apvalues,*apvaluestmp,*c_rmtaa;
1013: PetscInt offset,ii,pocol;
1014: const PetscInt *mappingindices;
1015: IS map;
1016: MPI_Comm comm;
1019: PetscObjectGetComm((PetscObject)A,&comm);
1020: if (!ptap->P_oth) SETERRQ(comm,PETSC_ERR_ARG_WRONGSTATE,"PtAP cannot be reused. Do not call MatFreeIntermediateDataStructures() or use '-mat_freeintermediatedatastructures'");
1022: MatZeroEntries(C);
1024: /* Get P_oth = ptap->P_oth and P_loc = ptap->P_loc */
1025: /*-----------------------------------------------------*/
1026: if (ptap->reuse == MAT_REUSE_MATRIX) {
1027: /* P_oth and P_loc are obtained in MatPtASymbolic() when reuse == MAT_INITIAL_MATRIX */
1028: MatGetBrowsOfAcols_MPIXAIJ(A,P,dof,MAT_REUSE_MATRIX,&ptap->P_oth);
1029: }
1030: PetscObjectQuery((PetscObject)ptap->P_oth,"aoffdiagtopothmapping",(PetscObject*)&map);
1031: MatGetLocalSize(p->B,NULL,&pon);
1032: pon *= dof;
1033: MatGetLocalSize(P,NULL,&pn);
1034: pn *= dof;
1036: PetscCalloc2(ptap->c_rmti[pon],&c_rmtj,ptap->c_rmti[pon],&c_rmta);
1037: MatGetLocalSize(A,&am,NULL);
1038: MatGetOwnershipRangeColumn(P,&pcstart,&pcend);
1039: pcstart *= dof;
1040: pcend *= dof;
1041: cmaxr = 0;
1042: for (i=0; i<pon; i++) {
1043: cmaxr = PetscMax(cmaxr,ptap->c_rmti[i+1]-ptap->c_rmti[i]);
1044: }
1045: cd = (Mat_SeqAIJ*)(c->A)->data;
1046: co = (Mat_SeqAIJ*)(c->B)->data;
1047: for (i=0; i<pn; i++) {
1048: cmaxr = PetscMax(cmaxr,(cd->i[i+1]-cd->i[i])+(co->i[i+1]-co->i[i]));
1049: }
1050: PetscCalloc4(cmaxr,&apindices,cmaxr,&apvalues,cmaxr,&apvaluestmp,pon,&c_rmtc);
1051: PetscHMapIVCreate(&hmap);
1052: PetscHMapIVResize(hmap,cmaxr);
1053: ISGetIndices(map,&mappingindices);
1054: for (i=0; i<am && (pon || pn); i++) {
1055: PetscHMapIVClear(hmap);
1056: offset = i%dof;
1057: ii = i/dof;
1058: nzi = po->i[ii+1] - po->i[ii];
1059: dnzi = pd->i[ii+1] - pd->i[ii];
1060: if (!nzi && !dnzi) continue;
1061: MatPtAPNumericComputeOneRowOfAP_private(A,P,ptap->P_oth,mappingindices,dof,i,hmap);
1062: voff = 0;
1063: PetscHMapIVGetPairs(hmap,&voff,apindices,apvalues);
1064: if (!voff) continue;
1066: /* Form remote C(ii, :) */
1067: poj = po->j + po->i[ii];
1068: poa = po->a + po->i[ii];
1069: for (j=0; j<nzi; j++) {
1070: pocol = poj[j]*dof+offset;
1071: c_rmtjj = c_rmtj + ptap->c_rmti[pocol];
1072: c_rmtaa = c_rmta + ptap->c_rmti[pocol];
1073: for (jj=0; jj<voff; jj++) {
1074: apvaluestmp[jj] = apvalues[jj]*poa[j];
1075: /*If the row is empty */
1076: if (!c_rmtc[pocol]) {
1077: c_rmtjj[jj] = apindices[jj];
1078: c_rmtaa[jj] = apvaluestmp[jj];
1079: c_rmtc[pocol]++;
1080: } else {
1081: PetscFindInt(apindices[jj],c_rmtc[pocol],c_rmtjj,&loc);
1082: if (loc>=0){ /* hit */
1083: c_rmtaa[loc] += apvaluestmp[jj];
1084: PetscLogFlops(1.0);
1085: } else { /* new element */
1086: loc = -(loc+1);
1087: /* Move data backward */
1088: for (kk=c_rmtc[pocol]; kk>loc; kk--) {
1089: c_rmtjj[kk] = c_rmtjj[kk-1];
1090: c_rmtaa[kk] = c_rmtaa[kk-1];
1091: }/* End kk */
1092: c_rmtjj[loc] = apindices[jj];
1093: c_rmtaa[loc] = apvaluestmp[jj];
1094: c_rmtc[pocol]++;
1095: }
1096: }
1097: } /* End jj */
1098: PetscLogFlops(voff);
1099: } /* End j */
1101: /* Form local C(ii, :) */
1102: pdj = pd->j + pd->i[ii];
1103: pda = pd->a + pd->i[ii];
1104: for (j=0; j<dnzi; j++) {
1105: row = pcstart + pdj[j] * dof + offset;
1106: for (jj=0; jj<voff; jj++) {
1107: apvaluestmp[jj] = apvalues[jj]*pda[j];
1108: }/* End kk */
1109: PetscLogFlops(voff);
1110: MatSetValues(C,1,&row,voff,apindices,apvaluestmp,ADD_VALUES);
1111: }/* End j */
1112: } /* End i */
1114: ISRestoreIndices(map,&mappingindices);
1115: PetscFree4(apindices,apvalues,apvaluestmp,c_rmtc);
1116: PetscHMapIVDestroy(&hmap);
1117: PetscCalloc2(ptap->c_othi[pn],&c_othj,ptap->c_othi[pn],&c_otha);
1119: PetscSFReduceBegin(ptap->sf,MPIU_INT,c_rmtj,c_othj,MPIU_REPLACE);
1120: PetscSFReduceBegin(ptap->sf,MPIU_SCALAR,c_rmta,c_otha,MPIU_REPLACE);
1121: PetscSFReduceEnd(ptap->sf,MPIU_INT,c_rmtj,c_othj,MPIU_REPLACE);
1122: PetscSFReduceEnd(ptap->sf,MPIU_SCALAR,c_rmta,c_otha,MPIU_REPLACE);
1123: PetscFree2(c_rmtj,c_rmta);
1125: /* Add contributions from remote */
1126: for (i = 0; i < pn; i++) {
1127: row = i + pcstart;
1128: 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);
1129: }
1130: PetscFree2(c_othj,c_otha);
1132: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
1133: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
1135: ptap->reuse = MAT_REUSE_MATRIX;
1137: /* supporting struct ptap consumes almost same amount of memory as C=PtAP, release it if C will not be updated by A and P */
1138: if (ptap->freestruct) {
1139: MatFreeIntermediateDataStructures(C);
1140: }
1141: return(0);
1142: }
1144: PetscErrorCode MatPtAPNumeric_MPIAIJ_MPIAIJ_allatonce_merged(Mat A,Mat P,Mat C)
1145: {
1146: PetscErrorCode ierr;
1150: MatPtAPNumeric_MPIAIJ_MPIXAIJ_allatonce_merged(A,P,1,C);
1151: return(0);
1152: }
1154: PetscErrorCode MatPtAPSymbolic_MPIAIJ_MPIXAIJ_allatonce(Mat A,Mat P,PetscInt dof,PetscReal fill,Mat *C)
1155: {
1156: Mat_APMPI *ptap;
1157: Mat_MPIAIJ *p=(Mat_MPIAIJ*)P->data,*c;
1158: MPI_Comm comm;
1159: Mat Cmpi;
1160: Mat_SeqAIJ *pd=(Mat_SeqAIJ*)p->A->data,*po=(Mat_SeqAIJ*)p->B->data;
1161: MatType mtype;
1162: PetscSF sf;
1163: PetscSFNode *iremote;
1164: PetscInt rootspacesize,*rootspace,*rootspaceoffsets,nleaves;
1165: const PetscInt *rootdegrees;
1166: PetscHSetI ht,oht,*hta,*hto;
1167: PetscInt pn,pon,*c_rmtc,i,j,nzi,htsize,htosize,*c_rmtj,off,*c_othj,rcvncols,sendncols,*c_rmtoffsets;
1168: PetscInt lidx,*rdj,col,pcstart,pcend,*dnz,*onz,am,arstart,arend,*poj,*pdj;
1169: PetscInt nalg=2,alg=0,offset,ii;
1170: PetscMPIInt owner;
1171: const PetscInt *mappingindices;
1172: PetscBool flg;
1173: const char *algTypes[2] = {"overlapping","merged"};
1174: IS map;
1175: PetscErrorCode ierr;
1178: PetscObjectGetComm((PetscObject)A,&comm);
1180: /* Create symbolic parallel matrix Cmpi */
1181: MatGetLocalSize(P,NULL,&pn);
1182: pn *= dof;
1183: MatCreate(comm,&Cmpi);
1184: MatGetType(A,&mtype);
1185: MatSetType(Cmpi,mtype);
1186: MatSetSizes(Cmpi,pn,pn,PETSC_DETERMINE,PETSC_DETERMINE);
1188: PetscNew(&ptap);
1189: ptap->reuse = MAT_INITIAL_MATRIX;
1190: ptap->algType = 2;
1192: /* Get P_oth by taking rows of P (= non-zero cols of local A) from other processors */
1193: MatGetBrowsOfAcols_MPIXAIJ(A,P,dof,MAT_INITIAL_MATRIX,&ptap->P_oth);
1194: PetscObjectQuery((PetscObject)ptap->P_oth,"aoffdiagtopothmapping",(PetscObject*)&map);
1195: /* This equals to the number of offdiag columns in P */
1196: MatGetLocalSize(p->B,NULL,&pon);
1197: pon *= dof;
1198: /* offsets */
1199: PetscMalloc1(pon+1,&ptap->c_rmti);
1200: /* The number of columns we will send to remote ranks */
1201: PetscMalloc1(pon,&c_rmtc);
1202: PetscMalloc1(pon,&hta);
1203: for (i=0; i<pon; i++) {
1204: PetscHSetICreate(&hta[i]);
1205: }
1206: MatGetLocalSize(A,&am,NULL);
1207: MatGetOwnershipRange(A,&arstart,&arend);
1208: /* Create hash table to merge all columns for C(i, :) */
1209: PetscHSetICreate(&ht);
1211: ISGetIndices(map,&mappingindices);
1212: ptap->c_rmti[0] = 0;
1213: /* 2) Pass 1: calculate the size for C_rmt (a matrix need to be sent to other processors) */
1214: for (i=0; i<am && pon; i++) {
1215: /* Form one row of AP */
1216: PetscHSetIClear(ht);
1217: offset = i%dof;
1218: ii = i/dof;
1219: /* If the off diag is empty, we should not do any calculation */
1220: nzi = po->i[ii+1] - po->i[ii];
1221: if (!nzi) continue;
1223: MatPtAPSymbolicComputeOneRowOfAP_private(A,P,ptap->P_oth,mappingindices,dof,i,ht,ht);
1224: PetscHSetIGetSize(ht,&htsize);
1225: /* If AP is empty, just continue */
1226: if (!htsize) continue;
1227: /* Form C(ii, :) */
1228: poj = po->j + po->i[ii];
1229: for (j=0; j<nzi; j++) {
1230: PetscHSetIUpdate(hta[poj[j]*dof+offset],ht);
1231: }
1232: }
1234: for (i=0; i<pon; i++) {
1235: PetscHSetIGetSize(hta[i],&htsize);
1236: ptap->c_rmti[i+1] = ptap->c_rmti[i] + htsize;
1237: c_rmtc[i] = htsize;
1238: }
1240: PetscMalloc1(ptap->c_rmti[pon],&c_rmtj);
1242: for (i=0; i<pon; i++) {
1243: off = 0;
1244: PetscHSetIGetElems(hta[i],&off,c_rmtj+ptap->c_rmti[i]);
1245: PetscHSetIDestroy(&hta[i]);
1246: }
1247: PetscFree(hta);
1249: PetscMalloc1(pon,&iremote);
1250: for (i=0; i<pon; i++) {
1251: owner = 0; lidx = 0;
1252: offset = i%dof;
1253: ii = i/dof;
1254: PetscLayoutFindOwnerIndex(P->cmap,p->garray[ii],&owner,&lidx);
1255: iremote[i].index = lidx*dof + offset;
1256: iremote[i].rank = owner;
1257: }
1259: PetscSFCreate(comm,&sf);
1260: PetscSFSetGraph(sf,pn,pon,NULL,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER);
1261: /* Reorder ranks properly so that the data handled by gather and scatter have the same order */
1262: PetscSFSetRankOrder(sf,PETSC_TRUE);
1263: PetscSFSetFromOptions(sf);
1264: PetscSFSetUp(sf);
1265: /* How many neighbors have contributions to my rows? */
1266: PetscSFComputeDegreeBegin(sf,&rootdegrees);
1267: PetscSFComputeDegreeEnd(sf,&rootdegrees);
1268: rootspacesize = 0;
1269: for (i = 0; i < pn; i++) {
1270: rootspacesize += rootdegrees[i];
1271: }
1272: PetscMalloc1(rootspacesize,&rootspace);
1273: PetscMalloc1(rootspacesize+1,&rootspaceoffsets);
1274: /* Get information from leaves
1275: * Number of columns other people contribute to my rows
1276: * */
1277: PetscSFGatherBegin(sf,MPIU_INT,c_rmtc,rootspace);
1278: PetscSFGatherEnd(sf,MPIU_INT,c_rmtc,rootspace);
1279: PetscFree(c_rmtc);
1280: PetscCalloc1(pn+1,&ptap->c_othi);
1281: /* The number of columns is received for each row */
1282: ptap->c_othi[0] = 0;
1283: rootspacesize = 0;
1284: rootspaceoffsets[0] = 0;
1285: for (i = 0; i < pn; i++) {
1286: rcvncols = 0;
1287: for (j = 0; j<rootdegrees[i]; j++) {
1288: rcvncols += rootspace[rootspacesize];
1289: rootspaceoffsets[rootspacesize+1] = rootspaceoffsets[rootspacesize] + rootspace[rootspacesize];
1290: rootspacesize++;
1291: }
1292: ptap->c_othi[i+1] = ptap->c_othi[i] + rcvncols;
1293: }
1294: PetscFree(rootspace);
1296: PetscMalloc1(pon,&c_rmtoffsets);
1297: PetscSFScatterBegin(sf,MPIU_INT,rootspaceoffsets,c_rmtoffsets);
1298: PetscSFScatterEnd(sf,MPIU_INT,rootspaceoffsets,c_rmtoffsets);
1299: PetscSFDestroy(&sf);
1300: PetscFree(rootspaceoffsets);
1302: PetscCalloc1(ptap->c_rmti[pon],&iremote);
1303: nleaves = 0;
1304: for (i = 0; i<pon; i++) {
1305: owner = 0;
1306: ii = i/dof;
1307: PetscLayoutFindOwnerIndex(P->cmap,p->garray[ii],&owner,NULL);
1308: sendncols = ptap->c_rmti[i+1] - ptap->c_rmti[i];
1309: for (j=0; j<sendncols; j++) {
1310: iremote[nleaves].rank = owner;
1311: iremote[nleaves++].index = c_rmtoffsets[i] + j;
1312: }
1313: }
1314: PetscFree(c_rmtoffsets);
1315: PetscCalloc1(ptap->c_othi[pn],&c_othj);
1317: PetscSFCreate(comm,&ptap->sf);
1318: PetscSFSetGraph(ptap->sf,ptap->c_othi[pn],nleaves,NULL,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER);
1319: PetscSFSetFromOptions(ptap->sf);
1320: /* One to one map */
1321: PetscSFReduceBegin(ptap->sf,MPIU_INT,c_rmtj,c_othj,MPIU_REPLACE);
1323: PetscMalloc2(pn,&dnz,pn,&onz);
1324: PetscHSetICreate(&oht);
1325: MatGetOwnershipRangeColumn(P,&pcstart,&pcend);
1326: pcstart *= dof;
1327: pcend *= dof;
1328: PetscMalloc2(pn,&hta,pn,&hto);
1329: for (i=0; i<pn; i++) {
1330: PetscHSetICreate(&hta[i]);
1331: PetscHSetICreate(&hto[i]);
1332: }
1333: /* Work on local part */
1334: /* 4) Pass 1: Estimate memory for C_loc */
1335: for (i=0; i<am && pn; i++) {
1336: PetscHSetIClear(ht);
1337: PetscHSetIClear(oht);
1338: offset = i%dof;
1339: ii = i/dof;
1340: nzi = pd->i[ii+1] - pd->i[ii];
1341: if (!nzi) continue;
1343: MatPtAPSymbolicComputeOneRowOfAP_private(A,P,ptap->P_oth,mappingindices,dof,i,ht,oht);
1344: PetscHSetIGetSize(ht,&htsize);
1345: PetscHSetIGetSize(oht,&htosize);
1346: if (!(htsize+htosize)) continue;
1347: /* Form C(ii, :) */
1348: pdj = pd->j + pd->i[ii];
1349: for (j=0; j<nzi; j++) {
1350: PetscHSetIUpdate(hta[pdj[j]*dof+offset],ht);
1351: PetscHSetIUpdate(hto[pdj[j]*dof+offset],oht);
1352: }
1353: }
1355: ISRestoreIndices(map,&mappingindices);
1357: PetscHSetIDestroy(&ht);
1358: PetscHSetIDestroy(&oht);
1360: /* Get remote data */
1361: PetscSFReduceEnd(ptap->sf,MPIU_INT,c_rmtj,c_othj,MPIU_REPLACE);
1362: PetscFree(c_rmtj);
1364: for (i = 0; i < pn; i++) {
1365: nzi = ptap->c_othi[i+1] - ptap->c_othi[i];
1366: rdj = c_othj + ptap->c_othi[i];
1367: for (j = 0; j < nzi; j++) {
1368: col = rdj[j];
1369: /* diag part */
1370: if (col>=pcstart && col<pcend) {
1371: PetscHSetIAdd(hta[i],col);
1372: } else { /* off diag */
1373: PetscHSetIAdd(hto[i],col);
1374: }
1375: }
1376: PetscHSetIGetSize(hta[i],&htsize);
1377: dnz[i] = htsize;
1378: PetscHSetIDestroy(&hta[i]);
1379: PetscHSetIGetSize(hto[i],&htsize);
1380: onz[i] = htsize;
1381: PetscHSetIDestroy(&hto[i]);
1382: }
1384: PetscFree2(hta,hto);
1385: PetscFree(c_othj);
1387: /* local sizes and preallocation */
1388: MatSetSizes(Cmpi,pn,pn,PETSC_DETERMINE,PETSC_DETERMINE);
1389: MatSetBlockSizes(Cmpi,dof>1? dof: P->cmap->bs,dof>1? dof: P->cmap->bs);
1390: MatMPIAIJSetPreallocation(Cmpi,0,dnz,0,onz);
1391: MatSetUp(Cmpi);
1392: PetscFree2(dnz,onz);
1394: /* attach the supporting struct to Cmpi for reuse */
1395: c = (Mat_MPIAIJ*)Cmpi->data;
1396: c->ap = ptap;
1397: ptap->duplicate = Cmpi->ops->duplicate;
1398: ptap->destroy = Cmpi->ops->destroy;
1399: ptap->view = Cmpi->ops->view;
1401: /* Cmpi is not ready for use - assembly will be done by MatPtAPNumeric() */
1402: Cmpi->assembled = PETSC_FALSE;
1403: /* pick an algorithm */
1404: PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MatPtAP","Mat");
1405: alg = 0;
1406: PetscOptionsEList("-matptap_allatonce_via","PtAP allatonce numeric approach","MatPtAP",algTypes,nalg,algTypes[alg],&alg,&flg);
1407: PetscOptionsEnd();
1408: switch (alg) {
1409: case 0:
1410: Cmpi->ops->ptapnumeric = MatPtAPNumeric_MPIAIJ_MPIAIJ_allatonce;
1411: break;
1412: case 1:
1413: Cmpi->ops->ptapnumeric = MatPtAPNumeric_MPIAIJ_MPIAIJ_allatonce_merged;
1414: break;
1415: default:
1416: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG," Unsupported allatonce numerical algorithm \n");
1417: }
1418: Cmpi->ops->destroy = MatDestroy_MPIAIJ_PtAP;
1419: Cmpi->ops->view = MatView_MPIAIJ_PtAP;
1420: Cmpi->ops->freeintermediatedatastructures = MatFreeIntermediateDataStructures_MPIAIJ_AP;
1421: *C = Cmpi;
1422: return(0);
1423: }
1425: PetscErrorCode MatPtAPSymbolic_MPIAIJ_MPIAIJ_allatonce(Mat A,Mat P,PetscReal fill,Mat *C)
1426: {
1427: PetscErrorCode ierr;
1431: MatPtAPSymbolic_MPIAIJ_MPIXAIJ_allatonce(A,P,1,fill,C);
1432: return(0);
1433: }
1435: PetscErrorCode MatPtAPSymbolic_MPIAIJ_MPIXAIJ_allatonce_merged(Mat A,Mat P,PetscInt dof,PetscReal fill,Mat *C)
1436: {
1437: Mat_APMPI *ptap;
1438: Mat_MPIAIJ *p=(Mat_MPIAIJ*)P->data,*c;
1439: MPI_Comm comm;
1440: Mat Cmpi;
1441: Mat_SeqAIJ *pd=(Mat_SeqAIJ*)p->A->data,*po=(Mat_SeqAIJ*)p->B->data;
1442: MatType mtype;
1443: PetscSF sf;
1444: PetscSFNode *iremote;
1445: PetscInt rootspacesize,*rootspace,*rootspaceoffsets,nleaves;
1446: const PetscInt *rootdegrees;
1447: PetscHSetI ht,oht,*hta,*hto,*htd;
1448: PetscInt pn,pon,*c_rmtc,i,j,nzi,dnzi,htsize,htosize,*c_rmtj,off,*c_othj,rcvncols,sendncols,*c_rmtoffsets;
1449: PetscInt lidx,*rdj,col,pcstart,pcend,*dnz,*onz,am,arstart,arend,*poj,*pdj;
1450: PetscInt nalg=2,alg=0,offset,ii;
1451: PetscMPIInt owner;
1452: PetscBool flg;
1453: const char *algTypes[2] = {"merged","overlapping"};
1454: const PetscInt *mappingindices;
1455: IS map;
1456: PetscErrorCode ierr;
1459: PetscObjectGetComm((PetscObject)A,&comm);
1461: /* Create symbolic parallel matrix Cmpi */
1462: MatGetLocalSize(P,NULL,&pn);
1463: pn *= dof;
1464: MatCreate(comm,&Cmpi);
1465: MatGetType(A,&mtype);
1466: MatSetType(Cmpi,mtype);
1467: MatSetSizes(Cmpi,pn,pn,PETSC_DETERMINE,PETSC_DETERMINE);
1469: PetscNew(&ptap);
1470: ptap->reuse = MAT_INITIAL_MATRIX;
1471: ptap->algType = 3;
1473: /* 0) Get P_oth by taking rows of P (= non-zero cols of local A) from other processors */
1474: MatGetBrowsOfAcols_MPIXAIJ(A,P,dof,MAT_INITIAL_MATRIX,&ptap->P_oth);
1475: PetscObjectQuery((PetscObject)ptap->P_oth,"aoffdiagtopothmapping",(PetscObject*)&map);
1477: /* This equals to the number of offdiag columns in P */
1478: MatGetLocalSize(p->B,NULL,&pon);
1479: pon *= dof;
1480: /* offsets */
1481: PetscMalloc1(pon+1,&ptap->c_rmti);
1482: /* The number of columns we will send to remote ranks */
1483: PetscMalloc1(pon,&c_rmtc);
1484: PetscMalloc1(pon,&hta);
1485: for (i=0; i<pon; i++) {
1486: PetscHSetICreate(&hta[i]);
1487: }
1488: MatGetLocalSize(A,&am,NULL);
1489: MatGetOwnershipRange(A,&arstart,&arend);
1490: /* Create hash table to merge all columns for C(i, :) */
1491: PetscHSetICreate(&ht);
1492: PetscHSetICreate(&oht);
1493: PetscMalloc2(pn,&htd,pn,&hto);
1494: for (i=0; i<pn; i++) {
1495: PetscHSetICreate(&htd[i]);
1496: PetscHSetICreate(&hto[i]);
1497: }
1499: ISGetIndices(map,&mappingindices);
1500: ptap->c_rmti[0] = 0;
1501: /* 2) Pass 1: calculate the size for C_rmt (a matrix need to be sent to other processors) */
1502: for (i=0; i<am && (pon || pn); i++) {
1503: /* Form one row of AP */
1504: PetscHSetIClear(ht);
1505: PetscHSetIClear(oht);
1506: offset = i%dof;
1507: ii = i/dof;
1508: /* If the off diag is empty, we should not do any calculation */
1509: nzi = po->i[ii+1] - po->i[ii];
1510: dnzi = pd->i[ii+1] - pd->i[ii];
1511: if (!nzi && !dnzi) continue;
1513: MatPtAPSymbolicComputeOneRowOfAP_private(A,P,ptap->P_oth,mappingindices,dof,i,ht,oht);
1514: PetscHSetIGetSize(ht,&htsize);
1515: PetscHSetIGetSize(oht,&htosize);
1516: /* If AP is empty, just continue */
1517: if (!(htsize+htosize)) continue;
1519: /* Form remote C(ii, :) */
1520: poj = po->j + po->i[ii];
1521: for (j=0; j<nzi; j++) {
1522: PetscHSetIUpdate(hta[poj[j]*dof+offset],ht);
1523: PetscHSetIUpdate(hta[poj[j]*dof+offset],oht);
1524: }
1526: /* Form local C(ii, :) */
1527: pdj = pd->j + pd->i[ii];
1528: for (j=0; j<dnzi; j++) {
1529: PetscHSetIUpdate(htd[pdj[j]*dof+offset],ht);
1530: PetscHSetIUpdate(hto[pdj[j]*dof+offset],oht);
1531: }
1532: }
1534: ISRestoreIndices(map,&mappingindices);
1536: PetscHSetIDestroy(&ht);
1537: PetscHSetIDestroy(&oht);
1539: for (i=0; i<pon; i++) {
1540: PetscHSetIGetSize(hta[i],&htsize);
1541: ptap->c_rmti[i+1] = ptap->c_rmti[i] + htsize;
1542: c_rmtc[i] = htsize;
1543: }
1545: PetscMalloc1(ptap->c_rmti[pon],&c_rmtj);
1547: for (i=0; i<pon; i++) {
1548: off = 0;
1549: PetscHSetIGetElems(hta[i],&off,c_rmtj+ptap->c_rmti[i]);
1550: PetscHSetIDestroy(&hta[i]);
1551: }
1552: PetscFree(hta);
1554: PetscMalloc1(pon,&iremote);
1555: for (i=0; i<pon; i++) {
1556: owner = 0; lidx = 0;
1557: offset = i%dof;
1558: ii = i/dof;
1559: PetscLayoutFindOwnerIndex(P->cmap,p->garray[ii],&owner,&lidx);
1560: iremote[i].index = lidx*dof+offset;
1561: iremote[i].rank = owner;
1562: }
1564: PetscSFCreate(comm,&sf);
1565: PetscSFSetGraph(sf,pn,pon,NULL,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER);
1566: /* Reorder ranks properly so that the data handled by gather and scatter have the same order */
1567: PetscSFSetRankOrder(sf,PETSC_TRUE);
1568: PetscSFSetFromOptions(sf);
1569: PetscSFSetUp(sf);
1570: /* How many neighbors have contributions to my rows? */
1571: PetscSFComputeDegreeBegin(sf,&rootdegrees);
1572: PetscSFComputeDegreeEnd(sf,&rootdegrees);
1573: rootspacesize = 0;
1574: for (i = 0; i < pn; i++) {
1575: rootspacesize += rootdegrees[i];
1576: }
1577: PetscMalloc1(rootspacesize,&rootspace);
1578: PetscMalloc1(rootspacesize+1,&rootspaceoffsets);
1579: /* Get information from leaves
1580: * Number of columns other people contribute to my rows
1581: * */
1582: PetscSFGatherBegin(sf,MPIU_INT,c_rmtc,rootspace);
1583: PetscSFGatherEnd(sf,MPIU_INT,c_rmtc,rootspace);
1584: PetscFree(c_rmtc);
1585: PetscMalloc1(pn+1,&ptap->c_othi);
1586: /* The number of columns is received for each row */
1587: ptap->c_othi[0] = 0;
1588: rootspacesize = 0;
1589: rootspaceoffsets[0] = 0;
1590: for (i = 0; i < pn; i++) {
1591: rcvncols = 0;
1592: for (j = 0; j<rootdegrees[i]; j++) {
1593: rcvncols += rootspace[rootspacesize];
1594: rootspaceoffsets[rootspacesize+1] = rootspaceoffsets[rootspacesize] + rootspace[rootspacesize];
1595: rootspacesize++;
1596: }
1597: ptap->c_othi[i+1] = ptap->c_othi[i] + rcvncols;
1598: }
1599: PetscFree(rootspace);
1601: PetscMalloc1(pon,&c_rmtoffsets);
1602: PetscSFScatterBegin(sf,MPIU_INT,rootspaceoffsets,c_rmtoffsets);
1603: PetscSFScatterEnd(sf,MPIU_INT,rootspaceoffsets,c_rmtoffsets);
1604: PetscSFDestroy(&sf);
1605: PetscFree(rootspaceoffsets);
1607: PetscCalloc1(ptap->c_rmti[pon],&iremote);
1608: nleaves = 0;
1609: for (i = 0; i<pon; i++) {
1610: owner = 0;
1611: ii = i/dof;
1612: PetscLayoutFindOwnerIndex(P->cmap,p->garray[ii],&owner,NULL);
1613: sendncols = ptap->c_rmti[i+1] - ptap->c_rmti[i];
1614: for (j=0; j<sendncols; j++) {
1615: iremote[nleaves].rank = owner;
1616: iremote[nleaves++].index = c_rmtoffsets[i] + j;
1617: }
1618: }
1619: PetscFree(c_rmtoffsets);
1620: PetscCalloc1(ptap->c_othi[pn],&c_othj);
1622: PetscSFCreate(comm,&ptap->sf);
1623: PetscSFSetGraph(ptap->sf,ptap->c_othi[pn],nleaves,NULL,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER);
1624: PetscSFSetFromOptions(ptap->sf);
1625: /* One to one map */
1626: PetscSFReduceBegin(ptap->sf,MPIU_INT,c_rmtj,c_othj,MPIU_REPLACE);
1627: /* Get remote data */
1628: PetscSFReduceEnd(ptap->sf,MPIU_INT,c_rmtj,c_othj,MPIU_REPLACE);
1629: PetscFree(c_rmtj);
1630: PetscMalloc2(pn,&dnz,pn,&onz);
1631: MatGetOwnershipRangeColumn(P,&pcstart,&pcend);
1632: pcstart *= dof;
1633: pcend *= dof;
1634: for (i = 0; i < pn; i++) {
1635: nzi = ptap->c_othi[i+1] - ptap->c_othi[i];
1636: rdj = c_othj + ptap->c_othi[i];
1637: for (j = 0; j < nzi; j++) {
1638: col = rdj[j];
1639: /* diag part */
1640: if (col>=pcstart && col<pcend) {
1641: PetscHSetIAdd(htd[i],col);
1642: } else { /* off diag */
1643: PetscHSetIAdd(hto[i],col);
1644: }
1645: }
1646: PetscHSetIGetSize(htd[i],&htsize);
1647: dnz[i] = htsize;
1648: PetscHSetIDestroy(&htd[i]);
1649: PetscHSetIGetSize(hto[i],&htsize);
1650: onz[i] = htsize;
1651: PetscHSetIDestroy(&hto[i]);
1652: }
1654: PetscFree2(htd,hto);
1655: PetscFree(c_othj);
1657: /* local sizes and preallocation */
1658: MatSetSizes(Cmpi,pn,pn,PETSC_DETERMINE,PETSC_DETERMINE);
1659: MatSetBlockSizes(Cmpi, dof>1? dof: P->cmap->bs,dof>1? dof: P->cmap->bs);
1660: MatMPIAIJSetPreallocation(Cmpi,0,dnz,0,onz);
1661: PetscFree2(dnz,onz);
1663: /* attach the supporting struct to Cmpi for reuse */
1664: c = (Mat_MPIAIJ*)Cmpi->data;
1665: c->ap = ptap;
1666: ptap->duplicate = Cmpi->ops->duplicate;
1667: ptap->destroy = Cmpi->ops->destroy;
1668: ptap->view = Cmpi->ops->view;
1670: /* Cmpi is not ready for use - assembly will be done by MatPtAPNumeric() */
1671: Cmpi->assembled = PETSC_FALSE;
1672: /* pick an algorithm */
1673: PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MatPtAP","Mat");
1674: alg = 0;
1675: PetscOptionsEList("-matptap_allatonce_via","PtAP allatonce numeric approach","MatPtAP",algTypes,nalg,algTypes[alg],&alg,&flg);
1676: PetscOptionsEnd();
1677: switch (alg) {
1678: case 0:
1679: Cmpi->ops->ptapnumeric = MatPtAPNumeric_MPIAIJ_MPIAIJ_allatonce_merged;
1680: break;
1681: case 1:
1682: Cmpi->ops->ptapnumeric = MatPtAPNumeric_MPIAIJ_MPIAIJ_allatonce;
1683: break;
1684: default:
1685: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG," Unsupported allatonce numerical algorithm \n");
1686: }
1687: Cmpi->ops->destroy = MatDestroy_MPIAIJ_PtAP;
1688: Cmpi->ops->view = MatView_MPIAIJ_PtAP;
1689: Cmpi->ops->freeintermediatedatastructures = MatFreeIntermediateDataStructures_MPIAIJ_AP;
1690: *C = Cmpi;
1691: return(0);
1692: }
1694: PetscErrorCode MatPtAPSymbolic_MPIAIJ_MPIAIJ_allatonce_merged(Mat A,Mat P,PetscReal fill,Mat *C)
1695: {
1696: PetscErrorCode ierr;
1700: MatPtAPSymbolic_MPIAIJ_MPIXAIJ_allatonce_merged(A,P,1,fill,C);
1701: return(0);
1702: }
1704: PetscErrorCode MatPtAPSymbolic_MPIAIJ_MPIAIJ(Mat A,Mat P,PetscReal fill,Mat *C)
1705: {
1706: PetscErrorCode ierr;
1707: Mat_APMPI *ptap;
1708: Mat_MPIAIJ *a=(Mat_MPIAIJ*)A->data,*p=(Mat_MPIAIJ*)P->data,*c;
1709: MPI_Comm comm;
1710: PetscMPIInt size,rank;
1711: Mat Cmpi;
1712: PetscFreeSpaceList free_space=NULL,current_space=NULL;
1713: PetscInt am=A->rmap->n,pm=P->rmap->n,pN=P->cmap->N,pn=P->cmap->n;
1714: PetscInt *lnk,i,k,pnz,row,nsend;
1715: PetscBT lnkbt;
1716: PetscMPIInt tagi,tagj,*len_si,*len_s,*len_ri,icompleted=0,nrecv;
1717: PetscInt **buf_rj,**buf_ri,**buf_ri_k;
1718: PetscInt len,proc,*dnz,*onz,*owners,nzi,nspacedouble;
1719: PetscInt nrows,*buf_s,*buf_si,*buf_si_i,**nextrow,**nextci;
1720: MPI_Request *swaits,*rwaits;
1721: MPI_Status *sstatus,rstatus;
1722: PetscLayout rowmap;
1723: PetscInt *owners_co,*coi,*coj; /* i and j array of (p->B)^T*A*P - used in the communication */
1724: PetscMPIInt *len_r,*id_r; /* array of length of comm->size, store send/recv matrix values */
1725: PetscInt *api,*apj,*Jptr,apnz,*prmap=p->garray,con,j,ap_rmax=0,Crmax,*aj,*ai,*pi;
1726: Mat_SeqAIJ *p_loc,*p_oth=NULL,*ad=(Mat_SeqAIJ*)(a->A)->data,*ao=NULL,*c_loc,*c_oth;
1727: PetscScalar *apv;
1728: PetscTable ta;
1729: MatType mtype;
1730: const char *prefix;
1731: #if defined(PETSC_USE_INFO)
1732: PetscReal apfill;
1733: #endif
1736: PetscObjectGetComm((PetscObject)A,&comm);
1737: MPI_Comm_size(comm,&size);
1738: MPI_Comm_rank(comm,&rank);
1740: if (size > 1) ao = (Mat_SeqAIJ*)(a->B)->data;
1742: /* create symbolic parallel matrix Cmpi */
1743: MatCreate(comm,&Cmpi);
1744: MatGetType(A,&mtype);
1745: MatSetType(Cmpi,mtype);
1747: /* Do dense axpy in MatPtAPNumeric_MPIAIJ_MPIAIJ() */
1748: Cmpi->ops->ptapnumeric = MatPtAPNumeric_MPIAIJ_MPIAIJ;
1750: /* create struct Mat_APMPI and attached it to C later */
1751: PetscNew(&ptap);
1752: ptap->reuse = MAT_INITIAL_MATRIX;
1753: ptap->algType = 1;
1755: /* get P_oth by taking rows of P (= non-zero cols of local A) from other processors */
1756: MatGetBrowsOfAoCols_MPIAIJ(A,P,MAT_INITIAL_MATRIX,&ptap->startsj_s,&ptap->startsj_r,&ptap->bufa,&ptap->P_oth);
1757: /* get P_loc by taking all local rows of P */
1758: MatMPIAIJGetLocalMat(P,MAT_INITIAL_MATRIX,&ptap->P_loc);
1760: /* (0) compute Rd = Pd^T, Ro = Po^T */
1761: /* --------------------------------- */
1762: MatTranspose(p->A,MAT_INITIAL_MATRIX,&ptap->Rd);
1763: MatTranspose(p->B,MAT_INITIAL_MATRIX,&ptap->Ro);
1765: /* (1) compute symbolic AP = A_loc*P = Ad*P_loc + Ao*P_oth (api,apj) */
1766: /* ----------------------------------------------------------------- */
1767: p_loc = (Mat_SeqAIJ*)(ptap->P_loc)->data;
1768: if (ptap->P_oth) p_oth = (Mat_SeqAIJ*)(ptap->P_oth)->data;
1770: /* create and initialize a linked list */
1771: PetscTableCreate(pn,pN,&ta); /* for compute AP_loc and Cmpi */
1772: MatRowMergeMax_SeqAIJ(p_loc,ptap->P_loc->rmap->N,ta);
1773: MatRowMergeMax_SeqAIJ(p_oth,ptap->P_oth->rmap->N,ta);
1774: PetscTableGetCount(ta,&Crmax); /* Crmax = nnz(sum of Prows) */
1775: /* printf("[%d] est %d, Crmax %d; pN %d\n",rank,5*(p_loc->rmax+p_oth->rmax + (PetscInt)(1.e-2*pN)),Crmax,pN); */
1777: PetscLLCondensedCreate(Crmax,pN,&lnk,&lnkbt);
1779: /* Initial FreeSpace size is fill*(nnz(A) + nnz(P)) */
1780: if (ao) {
1781: PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,PetscIntSumTruncate(ad->i[am],PetscIntSumTruncate(ao->i[am],p_loc->i[pm]))),&free_space);
1782: } else {
1783: PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,PetscIntSumTruncate(ad->i[am],p_loc->i[pm])),&free_space);
1784: }
1785: current_space = free_space;
1786: nspacedouble = 0;
1788: PetscMalloc1(am+1,&api);
1789: api[0] = 0;
1790: for (i=0; i<am; i++) {
1791: /* diagonal portion: Ad[i,:]*P */
1792: ai = ad->i; pi = p_loc->i;
1793: nzi = ai[i+1] - ai[i];
1794: aj = ad->j + ai[i];
1795: for (j=0; j<nzi; j++) {
1796: row = aj[j];
1797: pnz = pi[row+1] - pi[row];
1798: Jptr = p_loc->j + pi[row];
1799: /* add non-zero cols of P into the sorted linked list lnk */
1800: PetscLLCondensedAddSorted(pnz,Jptr,lnk,lnkbt);
1801: }
1802: /* off-diagonal portion: Ao[i,:]*P */
1803: if (ao) {
1804: ai = ao->i; pi = p_oth->i;
1805: nzi = ai[i+1] - ai[i];
1806: aj = ao->j + ai[i];
1807: for (j=0; j<nzi; j++) {
1808: row = aj[j];
1809: pnz = pi[row+1] - pi[row];
1810: Jptr = p_oth->j + pi[row];
1811: PetscLLCondensedAddSorted(pnz,Jptr,lnk,lnkbt);
1812: }
1813: }
1814: apnz = lnk[0];
1815: api[i+1] = api[i] + apnz;
1816: if (ap_rmax < apnz) ap_rmax = apnz;
1818: /* if free space is not available, double the total space in the list */
1819: if (current_space->local_remaining<apnz) {
1820: PetscFreeSpaceGet(PetscIntSumTruncate(apnz,current_space->total_array_size),¤t_space);
1821: nspacedouble++;
1822: }
1824: /* Copy data into free space, then initialize lnk */
1825: PetscLLCondensedClean(pN,apnz,current_space->array,lnk,lnkbt);
1827: current_space->array += apnz;
1828: current_space->local_used += apnz;
1829: current_space->local_remaining -= apnz;
1830: }
1831: /* Allocate space for apj and apv, initialize apj, and */
1832: /* destroy list of free space and other temporary array(s) */
1833: PetscMalloc2(api[am],&apj,api[am],&apv);
1834: PetscFreeSpaceContiguous(&free_space,apj);
1835: PetscLLDestroy(lnk,lnkbt);
1837: /* Create AP_loc for reuse */
1838: MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,am,pN,api,apj,apv,&ptap->AP_loc);
1840: #if defined(PETSC_USE_INFO)
1841: if (ao) {
1842: apfill = (PetscReal)api[am]/(ad->i[am]+ao->i[am]+p_loc->i[pm]+1);
1843: } else {
1844: apfill = (PetscReal)api[am]/(ad->i[am]+p_loc->i[pm]+1);
1845: }
1846: ptap->AP_loc->info.mallocs = nspacedouble;
1847: ptap->AP_loc->info.fill_ratio_given = fill;
1848: ptap->AP_loc->info.fill_ratio_needed = apfill;
1850: if (api[am]) {
1851: PetscInfo3(ptap->AP_loc,"Nonscalable algorithm, AP_loc reallocs %D; Fill ratio: given %g needed %g.\n",nspacedouble,(double)fill,(double)apfill);
1852: PetscInfo1(ptap->AP_loc,"Use MatPtAP(A,B,MatReuse,%g,&C) for best AP_loc performance.;\n",(double)apfill);
1853: } else {
1854: PetscInfo(ptap->AP_loc,"Nonscalable algorithm, AP_loc is empty \n");
1855: }
1856: #endif
1858: /* (2-1) compute symbolic Co = Ro*AP_loc */
1859: /* ------------------------------------ */
1860: MatGetOptionsPrefix(A,&prefix);
1861: MatSetOptionsPrefix(ptap->Ro,prefix);
1862: MatAppendOptionsPrefix(ptap->Ro,"inner_offdiag_");
1863: MatMatMultSymbolic_SeqAIJ_SeqAIJ(ptap->Ro,ptap->AP_loc,fill,&ptap->C_oth);
1865: /* (3) send coj of C_oth to other processors */
1866: /* ------------------------------------------ */
1867: /* determine row ownership */
1868: PetscLayoutCreate(comm,&rowmap);
1869: rowmap->n = pn;
1870: rowmap->bs = 1;
1871: PetscLayoutSetUp(rowmap);
1872: owners = rowmap->range;
1874: /* determine the number of messages to send, their lengths */
1875: PetscMalloc4(size,&len_s,size,&len_si,size,&sstatus,size+2,&owners_co);
1876: PetscArrayzero(len_s,size);
1877: PetscArrayzero(len_si,size);
1879: c_oth = (Mat_SeqAIJ*)ptap->C_oth->data;
1880: coi = c_oth->i; coj = c_oth->j;
1881: con = ptap->C_oth->rmap->n;
1882: proc = 0;
1883: for (i=0; i<con; i++) {
1884: while (prmap[i] >= owners[proc+1]) proc++;
1885: len_si[proc]++; /* num of rows in Co(=Pt*AP) to be sent to [proc] */
1886: len_s[proc] += coi[i+1] - coi[i]; /* num of nonzeros in Co to be sent to [proc] */
1887: }
1889: len = 0; /* max length of buf_si[], see (4) */
1890: owners_co[0] = 0;
1891: nsend = 0;
1892: for (proc=0; proc<size; proc++) {
1893: owners_co[proc+1] = owners_co[proc] + len_si[proc];
1894: if (len_s[proc]) {
1895: nsend++;
1896: len_si[proc] = 2*(len_si[proc] + 1); /* length of buf_si to be sent to [proc] */
1897: len += len_si[proc];
1898: }
1899: }
1901: /* determine the number and length of messages to receive for coi and coj */
1902: PetscGatherNumberOfMessages(comm,NULL,len_s,&nrecv);
1903: PetscGatherMessageLengths2(comm,nsend,nrecv,len_s,len_si,&id_r,&len_r,&len_ri);
1905: /* post the Irecv and Isend of coj */
1906: PetscCommGetNewTag(comm,&tagj);
1907: PetscPostIrecvInt(comm,tagj,nrecv,id_r,len_r,&buf_rj,&rwaits);
1908: PetscMalloc1(nsend+1,&swaits);
1909: for (proc=0, k=0; proc<size; proc++) {
1910: if (!len_s[proc]) continue;
1911: i = owners_co[proc];
1912: MPI_Isend(coj+coi[i],len_s[proc],MPIU_INT,proc,tagj,comm,swaits+k);
1913: k++;
1914: }
1916: /* (2-2) compute symbolic C_loc = Rd*AP_loc */
1917: /* ---------------------------------------- */
1918: MatSetOptionsPrefix(ptap->Rd,prefix);
1919: MatAppendOptionsPrefix(ptap->Rd,"inner_diag_");
1920: MatMatMultSymbolic_SeqAIJ_SeqAIJ(ptap->Rd,ptap->AP_loc,fill,&ptap->C_loc);
1921: c_loc = (Mat_SeqAIJ*)ptap->C_loc->data;
1923: /* receives coj are complete */
1924: for (i=0; i<nrecv; i++) {
1925: MPI_Waitany(nrecv,rwaits,&icompleted,&rstatus);
1926: }
1927: PetscFree(rwaits);
1928: if (nsend) {MPI_Waitall(nsend,swaits,sstatus);}
1930: /* add received column indices into ta to update Crmax */
1931: for (k=0; k<nrecv; k++) {/* k-th received message */
1932: Jptr = buf_rj[k];
1933: for (j=0; j<len_r[k]; j++) {
1934: PetscTableAdd(ta,*(Jptr+j)+1,1,INSERT_VALUES);
1935: }
1936: }
1937: PetscTableGetCount(ta,&Crmax);
1938: PetscTableDestroy(&ta);
1940: /* (4) send and recv coi */
1941: /*-----------------------*/
1942: PetscCommGetNewTag(comm,&tagi);
1943: PetscPostIrecvInt(comm,tagi,nrecv,id_r,len_ri,&buf_ri,&rwaits);
1944: PetscMalloc1(len+1,&buf_s);
1945: buf_si = buf_s; /* points to the beginning of k-th msg to be sent */
1946: for (proc=0,k=0; proc<size; proc++) {
1947: if (!len_s[proc]) continue;
1948: /* form outgoing message for i-structure:
1949: buf_si[0]: nrows to be sent
1950: [1:nrows]: row index (global)
1951: [nrows+1:2*nrows+1]: i-structure index
1952: */
1953: /*-------------------------------------------*/
1954: nrows = len_si[proc]/2 - 1; /* num of rows in Co to be sent to [proc] */
1955: buf_si_i = buf_si + nrows+1;
1956: buf_si[0] = nrows;
1957: buf_si_i[0] = 0;
1958: nrows = 0;
1959: for (i=owners_co[proc]; i<owners_co[proc+1]; i++) {
1960: nzi = coi[i+1] - coi[i];
1961: buf_si_i[nrows+1] = buf_si_i[nrows] + nzi; /* i-structure */
1962: buf_si[nrows+1] = prmap[i] -owners[proc]; /* local row index */
1963: nrows++;
1964: }
1965: MPI_Isend(buf_si,len_si[proc],MPIU_INT,proc,tagi,comm,swaits+k);
1966: k++;
1967: buf_si += len_si[proc];
1968: }
1969: for (i=0; i<nrecv; i++) {
1970: MPI_Waitany(nrecv,rwaits,&icompleted,&rstatus);
1971: }
1972: PetscFree(rwaits);
1973: if (nsend) {MPI_Waitall(nsend,swaits,sstatus);}
1975: PetscFree4(len_s,len_si,sstatus,owners_co);
1976: PetscFree(len_ri);
1977: PetscFree(swaits);
1978: PetscFree(buf_s);
1980: /* (5) compute the local portion of Cmpi */
1981: /* ------------------------------------------ */
1982: /* set initial free space to be Crmax, sufficient for holding nozeros in each row of Cmpi */
1983: PetscFreeSpaceGet(Crmax,&free_space);
1984: current_space = free_space;
1986: PetscMalloc3(nrecv,&buf_ri_k,nrecv,&nextrow,nrecv,&nextci);
1987: for (k=0; k<nrecv; k++) {
1988: buf_ri_k[k] = buf_ri[k]; /* beginning of k-th recved i-structure */
1989: nrows = *buf_ri_k[k];
1990: nextrow[k] = buf_ri_k[k] + 1; /* next row number of k-th recved i-structure */
1991: nextci[k] = buf_ri_k[k] + (nrows + 1); /* poins to the next i-structure of k-th recved i-structure */
1992: }
1994: MatPreallocateInitialize(comm,pn,pn,dnz,onz);
1995: PetscLLCondensedCreate(Crmax,pN,&lnk,&lnkbt);
1996: for (i=0; i<pn; i++) {
1997: /* add C_loc into Cmpi */
1998: nzi = c_loc->i[i+1] - c_loc->i[i];
1999: Jptr = c_loc->j + c_loc->i[i];
2000: PetscLLCondensedAddSorted(nzi,Jptr,lnk,lnkbt);
2002: /* add received col data into lnk */
2003: for (k=0; k<nrecv; k++) { /* k-th received message */
2004: if (i == *nextrow[k]) { /* i-th row */
2005: nzi = *(nextci[k]+1) - *nextci[k];
2006: Jptr = buf_rj[k] + *nextci[k];
2007: PetscLLCondensedAddSorted(nzi,Jptr,lnk,lnkbt);
2008: nextrow[k]++; nextci[k]++;
2009: }
2010: }
2011: nzi = lnk[0];
2013: /* copy data into free space, then initialize lnk */
2014: PetscLLCondensedClean(pN,nzi,current_space->array,lnk,lnkbt);
2015: MatPreallocateSet(i+owners[rank],nzi,current_space->array,dnz,onz);
2016: }
2017: PetscFree3(buf_ri_k,nextrow,nextci);
2018: PetscLLDestroy(lnk,lnkbt);
2019: PetscFreeSpaceDestroy(free_space);
2021: /* local sizes and preallocation */
2022: MatSetSizes(Cmpi,pn,pn,PETSC_DETERMINE,PETSC_DETERMINE);
2023: if (P->cmap->bs > 0) {
2024: PetscLayoutSetBlockSize(Cmpi->rmap,P->cmap->bs);
2025: PetscLayoutSetBlockSize(Cmpi->cmap,P->cmap->bs);
2026: }
2027: MatMPIAIJSetPreallocation(Cmpi,0,dnz,0,onz);
2028: MatPreallocateFinalize(dnz,onz);
2030: /* members in merge */
2031: PetscFree(id_r);
2032: PetscFree(len_r);
2033: PetscFree(buf_ri[0]);
2034: PetscFree(buf_ri);
2035: PetscFree(buf_rj[0]);
2036: PetscFree(buf_rj);
2037: PetscLayoutDestroy(&rowmap);
2039: /* attach the supporting struct to Cmpi for reuse */
2040: c = (Mat_MPIAIJ*)Cmpi->data;
2041: c->ap = ptap;
2042: ptap->duplicate = Cmpi->ops->duplicate;
2043: ptap->destroy = Cmpi->ops->destroy;
2044: ptap->view = Cmpi->ops->view;
2045: PetscCalloc1(pN,&ptap->apa);
2047: /* Cmpi is not ready for use - assembly will be done by MatPtAPNumeric() */
2048: Cmpi->assembled = PETSC_FALSE;
2049: Cmpi->ops->destroy = MatDestroy_MPIAIJ_PtAP;
2050: Cmpi->ops->view = MatView_MPIAIJ_PtAP;
2051: Cmpi->ops->freeintermediatedatastructures = MatFreeIntermediateDataStructures_MPIAIJ_AP;
2052: *C = Cmpi;
2053: return(0);
2054: }
2056: PetscErrorCode MatPtAPNumeric_MPIAIJ_MPIAIJ(Mat A,Mat P,Mat C)
2057: {
2058: PetscErrorCode ierr;
2059: Mat_MPIAIJ *a=(Mat_MPIAIJ*)A->data,*p=(Mat_MPIAIJ*)P->data,*c=(Mat_MPIAIJ*)C->data;
2060: Mat_SeqAIJ *ad=(Mat_SeqAIJ*)(a->A)->data,*ao=(Mat_SeqAIJ*)(a->B)->data;
2061: Mat_SeqAIJ *ap,*p_loc,*p_oth=NULL,*c_seq;
2062: Mat_APMPI *ptap = c->ap;
2063: Mat AP_loc,C_loc,C_oth;
2064: PetscInt i,rstart,rend,cm,ncols,row;
2065: PetscInt *api,*apj,am = A->rmap->n,j,col,apnz;
2066: PetscScalar *apa;
2067: const PetscInt *cols;
2068: const PetscScalar *vals;
2071: if (!ptap->AP_loc) {
2072: MPI_Comm comm;
2073: PetscObjectGetComm((PetscObject)C,&comm);
2074: SETERRQ(comm,PETSC_ERR_ARG_WRONGSTATE,"PtAP cannot be reused. Do not call MatFreeIntermediateDataStructures() or use '-mat_freeintermediatedatastructures'");
2075: }
2077: MatZeroEntries(C);
2078: /* 1) get R = Pd^T,Ro = Po^T */
2079: if (ptap->reuse == MAT_REUSE_MATRIX) {
2080: MatTranspose(p->A,MAT_REUSE_MATRIX,&ptap->Rd);
2081: MatTranspose(p->B,MAT_REUSE_MATRIX,&ptap->Ro);
2082: }
2084: /* 2) get AP_loc */
2085: AP_loc = ptap->AP_loc;
2086: ap = (Mat_SeqAIJ*)AP_loc->data;
2088: /* 2-1) get P_oth = ptap->P_oth and P_loc = ptap->P_loc */
2089: /*-----------------------------------------------------*/
2090: if (ptap->reuse == MAT_REUSE_MATRIX) {
2091: /* P_oth and P_loc are obtained in MatPtASymbolic() when reuse == MAT_INITIAL_MATRIX */
2092: MatGetBrowsOfAoCols_MPIAIJ(A,P,MAT_REUSE_MATRIX,&ptap->startsj_s,&ptap->startsj_r,&ptap->bufa,&ptap->P_oth);
2093: MatMPIAIJGetLocalMat(P,MAT_REUSE_MATRIX,&ptap->P_loc);
2094: }
2096: /* 2-2) compute numeric A_loc*P - dominating part */
2097: /* ---------------------------------------------- */
2098: /* get data from symbolic products */
2099: p_loc = (Mat_SeqAIJ*)(ptap->P_loc)->data;
2100: if (ptap->P_oth) {
2101: p_oth = (Mat_SeqAIJ*)(ptap->P_oth)->data;
2102: }
2103: apa = ptap->apa;
2104: api = ap->i;
2105: apj = ap->j;
2106: for (i=0; i<am; i++) {
2107: /* AP[i,:] = A[i,:]*P = Ad*P_loc Ao*P_oth */
2108: AProw_nonscalable(i,ad,ao,p_loc,p_oth,apa);
2109: apnz = api[i+1] - api[i];
2110: for (j=0; j<apnz; j++) {
2111: col = apj[j+api[i]];
2112: ap->a[j+ap->i[i]] = apa[col];
2113: apa[col] = 0.0;
2114: }
2115: }
2117: /* 3) C_loc = Rd*AP_loc, C_oth = Ro*AP_loc */
2118: ((ptap->C_loc)->ops->matmultnumeric)(ptap->Rd,AP_loc,ptap->C_loc);
2119: ((ptap->C_oth)->ops->matmultnumeric)(ptap->Ro,AP_loc,ptap->C_oth);
2120: C_loc = ptap->C_loc;
2121: C_oth = ptap->C_oth;
2123: /* add C_loc and Co to to C */
2124: MatGetOwnershipRange(C,&rstart,&rend);
2126: /* C_loc -> C */
2127: cm = C_loc->rmap->N;
2128: c_seq = (Mat_SeqAIJ*)C_loc->data;
2129: cols = c_seq->j;
2130: vals = c_seq->a;
2133: /* The (fast) MatSetValues_MPIAIJ_CopyFromCSRFormat function can only be used when C->was_assembled is PETSC_FALSE and */
2134: /* when there are no off-processor parts. */
2135: /* If was_assembled is true, then the statement aj[rowstart_diag+dnz_row] = mat_j[col] - cstart; in MatSetValues_MPIAIJ_CopyFromCSRFormat */
2136: /* is no longer true. Then the more complex function MatSetValues_MPIAIJ() has to be used, where the column index is looked up from */
2137: /* a table, and other, more complex stuff has to be done. */
2138: if (C->assembled) {
2139: C->was_assembled = PETSC_TRUE;
2140: C->assembled = PETSC_FALSE;
2141: }
2142: if (C->was_assembled) {
2143: for (i=0; i<cm; i++) {
2144: ncols = c_seq->i[i+1] - c_seq->i[i];
2145: row = rstart + i;
2146: MatSetValues_MPIAIJ(C,1,&row,ncols,cols,vals,ADD_VALUES);
2147: cols += ncols; vals += ncols;
2148: }
2149: } else {
2150: MatSetValues_MPIAIJ_CopyFromCSRFormat(C,c_seq->j,c_seq->i,c_seq->a);
2151: }
2153: /* Co -> C, off-processor part */
2154: cm = C_oth->rmap->N;
2155: c_seq = (Mat_SeqAIJ*)C_oth->data;
2156: cols = c_seq->j;
2157: vals = c_seq->a;
2158: for (i=0; i<cm; i++) {
2159: ncols = c_seq->i[i+1] - c_seq->i[i];
2160: row = p->garray[i];
2161: MatSetValues(C,1,&row,ncols,cols,vals,ADD_VALUES);
2162: cols += ncols; vals += ncols;
2163: }
2165: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
2166: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
2168: ptap->reuse = MAT_REUSE_MATRIX;
2170: /* supporting struct ptap consumes almost same amount of memory as C=PtAP, release it if C will not be updated by A and P */
2171: if (ptap->freestruct) {
2172: MatFreeIntermediateDataStructures(C);
2173: }
2174: return(0);
2175: }