Actual source code: mpiptap.c
1: /*
2: Defines projective product routines where A is a MPIAIJ matrix
3: C = P^T * A * P
4: */
6: #include <../src/mat/impls/aij/seq/aij.h>
7: #include <../src/mat/utils/freespace.h>
8: #include <../src/mat/impls/aij/mpi/mpiaij.h>
9: #include <petscbt.h>
10: #include <petsctime.h>
11: #include <petsc/private/hashmapiv.h>
12: #include <petsc/private/hashseti.h>
13: #include <petscsf.h>
15: static PetscErrorCode MatView_MPIAIJ_PtAP(Mat A, PetscViewer viewer)
16: {
17: PetscBool isascii;
18: PetscViewerFormat format;
19: MatProductCtx_APMPI *ptap;
21: PetscFunctionBegin;
22: MatCheckProduct(A, 1);
23: ptap = (MatProductCtx_APMPI *)A->product->data;
24: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
25: if (isascii) {
26: PetscCall(PetscViewerGetFormat(viewer, &format));
27: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
28: if (ptap->algType == 0) {
29: PetscCall(PetscViewerASCIIPrintf(viewer, "using scalable MatPtAP() implementation\n"));
30: } else if (ptap->algType == 1) {
31: PetscCall(PetscViewerASCIIPrintf(viewer, "using nonscalable MatPtAP() implementation\n"));
32: } else if (ptap->algType == 2) {
33: PetscCall(PetscViewerASCIIPrintf(viewer, "using allatonce MatPtAP() implementation\n"));
34: } else if (ptap->algType == 3) {
35: PetscCall(PetscViewerASCIIPrintf(viewer, "using merged allatonce MatPtAP() implementation\n"));
36: }
37: }
38: }
39: PetscFunctionReturn(PETSC_SUCCESS);
40: }
42: PetscErrorCode MatProductCtxDestroy_MPIAIJ_PtAP(PetscCtxRt data)
43: {
44: MatProductCtx_APMPI *ptap = *(MatProductCtx_APMPI **)data;
45: MatMergeSeqsToMPI *merge;
47: PetscFunctionBegin;
48: PetscCall(PetscFree2(ptap->startsj_s, ptap->startsj_r));
49: PetscCall(PetscFree(ptap->bufa));
50: PetscCall(MatDestroy(&ptap->P_loc));
51: PetscCall(MatDestroy(&ptap->P_oth));
52: PetscCall(MatDestroy(&ptap->A_loc)); /* used by MatTransposeMatMult() */
53: PetscCall(MatDestroy(&ptap->Rd));
54: PetscCall(MatDestroy(&ptap->Ro));
55: if (ptap->AP_loc) { /* used by alg_rap */
56: Mat_SeqAIJ *ap = (Mat_SeqAIJ *)ptap->AP_loc->data;
57: PetscCall(PetscFree(ap->i));
58: PetscCall(PetscFree2(ap->j, ap->a));
59: PetscCall(MatDestroy(&ptap->AP_loc));
60: } else { /* used by alg_ptap */
61: PetscCall(PetscFree(ptap->api));
62: PetscCall(PetscFree(ptap->apj));
63: }
64: PetscCall(MatDestroy(&ptap->C_loc));
65: PetscCall(MatDestroy(&ptap->C_oth));
66: PetscCall(PetscFree(ptap->apa));
68: PetscCall(MatDestroy(&ptap->Pt));
70: merge = ptap->merge;
71: if (merge) { /* used by alg_ptap */
72: PetscCall(PetscFree(merge->id_r));
73: PetscCall(PetscFree(merge->len_s));
74: PetscCall(PetscFree(merge->len_r));
75: PetscCall(PetscFree(merge->bi));
76: PetscCall(PetscFree(merge->bj));
77: PetscCall(PetscFree(merge->buf_ri[0]));
78: PetscCall(PetscFree(merge->buf_ri));
79: PetscCall(PetscFree(merge->buf_rj[0]));
80: PetscCall(PetscFree(merge->buf_rj));
81: PetscCall(PetscFree(merge->coi));
82: PetscCall(PetscFree(merge->coj));
83: PetscCall(PetscFree(merge->owners_co));
84: PetscCall(PetscLayoutDestroy(&merge->rowmap));
85: PetscCall(PetscFree(ptap->merge));
86: }
87: PetscCall(ISLocalToGlobalMappingDestroy(&ptap->ltog));
89: PetscCall(PetscSFDestroy(&ptap->sf));
90: PetscCall(PetscFree(ptap->c_othi));
91: PetscCall(PetscFree(ptap->c_rmti));
92: PetscCall(PetscFree(ptap));
93: PetscFunctionReturn(PETSC_SUCCESS);
94: }
96: PetscErrorCode MatPtAPNumeric_MPIAIJ_MPIAIJ_scalable(Mat A, Mat P, Mat C)
97: {
98: Mat_MPIAIJ *a = (Mat_MPIAIJ *)A->data, *p = (Mat_MPIAIJ *)P->data;
99: Mat_SeqAIJ *ad = (Mat_SeqAIJ *)a->A->data, *ao = (Mat_SeqAIJ *)a->B->data;
100: Mat_SeqAIJ *ap, *p_loc, *p_oth = NULL, *c_seq;
101: MatProductCtx_APMPI *ptap;
102: Mat AP_loc, C_loc, C_oth;
103: PetscInt i, rstart, rend, cm, ncols, row, *api, *apj, am = A->rmap->n, apnz, nout;
104: PetscScalar *apa;
105: const PetscInt *cols;
106: const PetscScalar *vals;
108: PetscFunctionBegin;
109: MatCheckProduct(C, 3);
110: ptap = (MatProductCtx_APMPI *)C->product->data;
111: PetscCheck(ptap, PetscObjectComm((PetscObject)C), PETSC_ERR_ARG_WRONGSTATE, "PtAP cannot be computed. Missing data");
112: PetscCheck(ptap->AP_loc, PetscObjectComm((PetscObject)C), PETSC_ERR_ARG_WRONGSTATE, "PtAP cannot be reused. Do not call MatProductClear()");
114: PetscCall(MatZeroEntries(C));
116: /* 1) get R = Pd^T,Ro = Po^T */
117: if (ptap->reuse == MAT_REUSE_MATRIX) {
118: PetscCall(MatTranspose(p->A, MAT_REUSE_MATRIX, &ptap->Rd));
119: PetscCall(MatTranspose(p->B, MAT_REUSE_MATRIX, &ptap->Ro));
120: }
122: /* 2) get AP_loc */
123: AP_loc = ptap->AP_loc;
124: ap = (Mat_SeqAIJ *)AP_loc->data;
126: /* 2-1) get P_oth = ptap->P_oth and P_loc = ptap->P_loc */
127: if (ptap->reuse == MAT_REUSE_MATRIX) {
128: /* P_oth and P_loc are obtained in MatPtASymbolic() when reuse == MAT_INITIAL_MATRIX */
129: PetscCall(MatGetBrowsOfAoCols_MPIAIJ(A, P, MAT_REUSE_MATRIX, &ptap->startsj_s, &ptap->startsj_r, &ptap->bufa, &ptap->P_oth));
130: PetscCall(MatMPIAIJGetLocalMat(P, MAT_REUSE_MATRIX, &ptap->P_loc));
131: }
133: /* 2-2) compute numeric A_loc*P - dominating part */
134: /* get data from symbolic products */
135: p_loc = (Mat_SeqAIJ *)ptap->P_loc->data;
136: if (ptap->P_oth) p_oth = (Mat_SeqAIJ *)ptap->P_oth->data;
138: api = ap->i;
139: apj = ap->j;
140: PetscCall(ISLocalToGlobalMappingApply(ptap->ltog, api[AP_loc->rmap->n], apj, apj));
141: for (i = 0; i < am; i++) {
142: /* AP[i,:] = A[i,:]*P = Ad*P_loc Ao*P_oth */
143: apnz = api[i + 1] - api[i];
144: apa = ap->a + api[i];
145: PetscCall(PetscArrayzero(apa, apnz));
146: AProw_scalable(i, ad, ao, p_loc, p_oth, api, apj, apa);
147: }
148: PetscCall(ISGlobalToLocalMappingApply(ptap->ltog, IS_GTOLM_DROP, api[AP_loc->rmap->n], apj, &nout, apj));
149: PetscCheck(api[AP_loc->rmap->n] == nout, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Incorrect mapping %" PetscInt_FMT " != %" PetscInt_FMT, api[AP_loc->rmap->n], nout);
151: /* 3) C_loc = Rd*AP_loc, C_oth = Ro*AP_loc */
152: /* Always use scalable version since we are in the MPI scalable version */
153: PetscCall(MatMatMultNumeric_SeqAIJ_SeqAIJ_Scalable(ptap->Rd, AP_loc, ptap->C_loc));
154: PetscCall(MatMatMultNumeric_SeqAIJ_SeqAIJ_Scalable(ptap->Ro, AP_loc, ptap->C_oth));
156: C_loc = ptap->C_loc;
157: C_oth = ptap->C_oth;
159: /* add C_loc and Co to C */
160: PetscCall(MatGetOwnershipRange(C, &rstart, &rend));
162: /* C_loc -> C */
163: cm = C_loc->rmap->N;
164: c_seq = (Mat_SeqAIJ *)C_loc->data;
165: cols = c_seq->j;
166: vals = c_seq->a;
167: PetscCall(ISLocalToGlobalMappingApply(ptap->ltog, c_seq->i[C_loc->rmap->n], c_seq->j, c_seq->j));
169: /* The (fast) MatSetValues_MPIAIJ_CopyFromCSRFormat function can only be used when C->was_assembled is PETSC_FALSE and */
170: /* when there are no off-processor parts. */
171: /* If was_assembled is true, then the statement aj[rowstart_diag+dnz_row] = mat_j[col] - cstart; in MatSetValues_MPIAIJ_CopyFromCSRFormat */
172: /* is no longer true. Then the more complex function MatSetValues_MPIAIJ() has to be used, where the column index is looked up from */
173: /* a table, and other, more complex stuff has to be done. */
174: if (C->assembled) {
175: C->was_assembled = PETSC_TRUE;
176: C->assembled = PETSC_FALSE;
177: }
178: if (C->was_assembled) {
179: for (i = 0; i < cm; i++) {
180: ncols = c_seq->i[i + 1] - c_seq->i[i];
181: row = rstart + i;
182: PetscCall(MatSetValues_MPIAIJ(C, 1, &row, ncols, cols, vals, ADD_VALUES));
183: cols += ncols;
184: vals += ncols;
185: }
186: } else {
187: PetscCall(MatSetValues_MPIAIJ_CopyFromCSRFormat(C, c_seq->j, c_seq->i, c_seq->a));
188: }
189: PetscCall(ISGlobalToLocalMappingApply(ptap->ltog, IS_GTOLM_DROP, c_seq->i[C_loc->rmap->n], c_seq->j, &nout, c_seq->j));
190: PetscCheck(c_seq->i[C_loc->rmap->n] == nout, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Incorrect mapping %" PetscInt_FMT " != %" PetscInt_FMT, c_seq->i[C_loc->rmap->n], nout);
192: /* Co -> C, off-processor part */
193: cm = C_oth->rmap->N;
194: c_seq = (Mat_SeqAIJ *)C_oth->data;
195: cols = c_seq->j;
196: vals = c_seq->a;
197: PetscCall(ISLocalToGlobalMappingApply(ptap->ltog, c_seq->i[C_oth->rmap->n], c_seq->j, c_seq->j));
198: for (i = 0; i < cm; i++) {
199: ncols = c_seq->i[i + 1] - c_seq->i[i];
200: row = p->garray[i];
201: PetscCall(MatSetValues(C, 1, &row, ncols, cols, vals, ADD_VALUES));
202: cols += ncols;
203: vals += ncols;
204: }
205: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
206: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
208: ptap->reuse = MAT_REUSE_MATRIX;
210: PetscCall(ISGlobalToLocalMappingApply(ptap->ltog, IS_GTOLM_DROP, c_seq->i[C_oth->rmap->n], c_seq->j, &nout, c_seq->j));
211: PetscCheck(c_seq->i[C_oth->rmap->n] == nout, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Incorrect mapping %" PetscInt_FMT " != %" PetscInt_FMT, c_seq->i[C_loc->rmap->n], nout);
212: PetscFunctionReturn(PETSC_SUCCESS);
213: }
215: PetscErrorCode MatPtAPSymbolic_MPIAIJ_MPIAIJ_scalable(Mat A, Mat P, PetscReal fill, Mat Cmpi)
216: {
217: MatProductCtx_APMPI *ptap;
218: Mat_MPIAIJ *a = (Mat_MPIAIJ *)A->data, *p = (Mat_MPIAIJ *)P->data;
219: MPI_Comm comm;
220: PetscMPIInt size, rank;
221: Mat P_loc, P_oth;
222: PetscFreeSpaceList free_space = NULL, current_space = NULL;
223: PetscInt am = A->rmap->n, pm = P->rmap->n, pN = P->cmap->N, pn = P->cmap->n;
224: PetscInt *lnk, i, k, pnz, row;
225: PetscMPIInt tagi, tagj, *len_si, *len_s, *len_ri, nrecv, nsend, proc;
226: PETSC_UNUSED PetscMPIInt icompleted = 0;
227: PetscInt **buf_rj, **buf_ri, **buf_ri_k;
228: const PetscInt *owners;
229: PetscInt len, *dnz, *onz, nzi, nspacedouble;
230: PetscInt nrows, *buf_s, *buf_si, *buf_si_i, **nextrow, **nextci;
231: MPI_Request *swaits, *rwaits;
232: MPI_Status *sstatus, rstatus;
233: PetscLayout rowmap;
234: PetscInt *owners_co, *coi, *coj; /* i and j array of (p->B)^T*A*P - used in the communication */
235: PetscMPIInt *len_r, *id_r; /* array of length of comm->size, store send/recv matrix values */
236: PetscInt *api, *apj, *Jptr, apnz, *prmap = p->garray, con, j, Crmax, *aj, *ai, *pi, nout;
237: Mat_SeqAIJ *p_loc, *p_oth = NULL, *ad = (Mat_SeqAIJ *)a->A->data, *ao = NULL, *c_loc, *c_oth;
238: PetscScalar *apv;
239: PetscHMapI ta;
240: MatType mtype;
241: const char *prefix;
242: PetscReal apfill;
244: PetscFunctionBegin;
245: MatCheckProduct(Cmpi, 4);
246: PetscCheck(!Cmpi->product->data, PetscObjectComm((PetscObject)Cmpi), PETSC_ERR_PLIB, "Product data not empty");
247: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
248: PetscCallMPI(MPI_Comm_size(comm, &size));
249: PetscCallMPI(MPI_Comm_rank(comm, &rank));
251: if (size > 1) ao = (Mat_SeqAIJ *)a->B->data;
253: /* create symbolic parallel matrix Cmpi */
254: PetscCall(MatGetType(A, &mtype));
255: PetscCall(MatSetType(Cmpi, mtype));
257: /* create struct MatProductCtx_APMPI and attached it to C later */
258: PetscCall(PetscNew(&ptap));
259: ptap->reuse = MAT_INITIAL_MATRIX;
260: ptap->algType = 0;
262: /* get P_oth by taking rows of P (= non-zero cols of local A) from other processors */
263: PetscCall(MatGetBrowsOfAoCols_MPIAIJ(A, P, MAT_INITIAL_MATRIX, &ptap->startsj_s, &ptap->startsj_r, &ptap->bufa, &P_oth));
264: /* get P_loc by taking all local rows of P */
265: PetscCall(MatMPIAIJGetLocalMat(P, MAT_INITIAL_MATRIX, &P_loc));
267: ptap->P_loc = P_loc;
268: ptap->P_oth = P_oth;
270: /* (0) compute Rd = Pd^T, Ro = Po^T */
271: PetscCall(MatTranspose(p->A, MAT_INITIAL_MATRIX, &ptap->Rd));
272: PetscCall(MatTranspose(p->B, MAT_INITIAL_MATRIX, &ptap->Ro));
274: /* (1) compute symbolic AP = A_loc*P = Ad*P_loc + Ao*P_oth (api,apj) */
275: p_loc = (Mat_SeqAIJ *)P_loc->data;
276: if (P_oth) p_oth = (Mat_SeqAIJ *)P_oth->data;
278: /* create and initialize a linked list */
279: PetscCall(PetscHMapICreateWithSize(pn, &ta)); /* for compute AP_loc and Cmpi */
280: MatRowMergeMax_SeqAIJ(p_loc, P_loc->rmap->N, ta);
281: MatRowMergeMax_SeqAIJ(p_oth, P_oth->rmap->N, ta);
282: PetscCall(PetscHMapIGetSize(ta, &Crmax)); /* Crmax = nnz(sum of Prows) */
284: PetscCall(PetscLLCondensedCreate_Scalable(Crmax, &lnk));
286: /* Initial FreeSpace size is fill*(nnz(A) + nnz(P)) */
287: if (ao) {
288: PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(fill, PetscIntSumTruncate(ad->i[am], PetscIntSumTruncate(ao->i[am], p_loc->i[pm]))), &free_space));
289: } else {
290: PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(fill, PetscIntSumTruncate(ad->i[am], p_loc->i[pm])), &free_space));
291: }
292: current_space = free_space;
293: nspacedouble = 0;
295: PetscCall(PetscMalloc1(am + 1, &api));
296: api[0] = 0;
297: for (i = 0; i < am; i++) {
298: /* diagonal portion: Ad[i,:]*P */
299: ai = ad->i;
300: pi = p_loc->i;
301: nzi = ai[i + 1] - ai[i];
302: aj = PetscSafePointerPlusOffset(ad->j, ai[i]);
303: for (j = 0; j < nzi; j++) {
304: row = aj[j];
305: pnz = pi[row + 1] - pi[row];
306: Jptr = p_loc->j + pi[row];
307: /* add non-zero cols of P into the sorted linked list lnk */
308: PetscCall(PetscLLCondensedAddSorted_Scalable(pnz, Jptr, lnk));
309: }
310: /* off-diagonal portion: Ao[i,:]*P */
311: if (ao) {
312: ai = ao->i;
313: pi = p_oth->i;
314: nzi = ai[i + 1] - ai[i];
315: aj = PetscSafePointerPlusOffset(ao->j, ai[i]);
316: for (j = 0; j < nzi; j++) {
317: row = aj[j];
318: pnz = pi[row + 1] - pi[row];
319: Jptr = p_oth->j + pi[row];
320: PetscCall(PetscLLCondensedAddSorted_Scalable(pnz, Jptr, lnk));
321: }
322: }
323: apnz = lnk[0];
324: api[i + 1] = api[i] + apnz;
326: /* if free space is not available, double the total space in the list */
327: if (current_space->local_remaining < apnz) {
328: PetscCall(PetscFreeSpaceGet(PetscIntSumTruncate(apnz, current_space->total_array_size), ¤t_space));
329: nspacedouble++;
330: }
332: /* Copy data into free space, then initialize lnk */
333: PetscCall(PetscLLCondensedClean_Scalable(apnz, current_space->array, lnk));
335: current_space->array += apnz;
336: current_space->local_used += apnz;
337: current_space->local_remaining -= apnz;
338: }
339: /* Allocate space for apj and apv, initialize apj, and */
340: /* destroy list of free space and other temporary array(s) */
341: PetscCall(PetscCalloc2(api[am], &apj, api[am], &apv));
342: PetscCall(PetscFreeSpaceContiguous(&free_space, apj));
343: PetscCall(PetscLLCondensedDestroy_Scalable(lnk));
345: /* Create AP_loc for reuse */
346: PetscCall(MatCreateSeqAIJWithArrays(PETSC_COMM_SELF, am, pN, api, apj, apv, &ptap->AP_loc));
347: PetscCall(MatSeqAIJCompactOutExtraColumns_SeqAIJ(ptap->AP_loc, &ptap->ltog));
349: if (PetscDefined(USE_INFO)) {
350: if (ao) apfill = (PetscReal)api[am] / (ad->i[am] + ao->i[am] + p_loc->i[pm] + 1);
351: else apfill = (PetscReal)api[am] / (ad->i[am] + p_loc->i[pm] + 1);
352: ptap->AP_loc->info.mallocs = nspacedouble;
353: ptap->AP_loc->info.fill_ratio_given = fill;
354: ptap->AP_loc->info.fill_ratio_needed = apfill;
356: if (api[am]) {
357: PetscCall(PetscInfo(ptap->AP_loc, "Scalable algorithm, AP_loc reallocs %" PetscInt_FMT "; Fill ratio: given %g needed %g.\n", nspacedouble, (double)fill, (double)apfill));
358: PetscCall(PetscInfo(ptap->AP_loc, "Use MatPtAP(A,B,MatReuse,%g,&C) for best AP_loc performance.;\n", (double)apfill));
359: } else PetscCall(PetscInfo(ptap->AP_loc, "Scalable algorithm, AP_loc is empty \n"));
360: }
362: /* (2-1) compute symbolic Co = Ro*AP_loc */
363: PetscCall(MatProductCreate(ptap->Ro, ptap->AP_loc, NULL, &ptap->C_oth));
364: PetscCall(MatGetOptionsPrefix(A, &prefix));
365: PetscCall(MatSetOptionsPrefix(ptap->C_oth, prefix));
366: PetscCall(MatAppendOptionsPrefix(ptap->C_oth, "inner_offdiag_"));
368: PetscCall(MatProductSetType(ptap->C_oth, MATPRODUCT_AB));
369: PetscCall(MatProductSetAlgorithm(ptap->C_oth, "sorted"));
370: PetscCall(MatProductSetFill(ptap->C_oth, fill));
371: PetscCall(MatProductSetFromOptions(ptap->C_oth));
372: PetscCall(MatProductSymbolic(ptap->C_oth));
374: /* (3) send coj of C_oth to other processors */
375: /* determine row ownership */
376: PetscCall(PetscLayoutCreate(comm, &rowmap));
377: PetscCall(PetscLayoutSetLocalSize(rowmap, pn));
378: PetscCall(PetscLayoutSetBlockSize(rowmap, 1));
379: PetscCall(PetscLayoutSetUp(rowmap));
380: PetscCall(PetscLayoutGetRanges(rowmap, &owners));
382: /* determine the number of messages to send, their lengths */
383: PetscCall(PetscMalloc4(size, &len_s, size, &len_si, size, &sstatus, size + 2, &owners_co));
384: PetscCall(PetscArrayzero(len_s, size));
385: PetscCall(PetscArrayzero(len_si, size));
387: c_oth = (Mat_SeqAIJ *)ptap->C_oth->data;
388: coi = c_oth->i;
389: coj = c_oth->j;
390: con = ptap->C_oth->rmap->n;
391: proc = 0;
392: PetscCall(ISLocalToGlobalMappingApply(ptap->ltog, coi[con], coj, coj));
393: for (i = 0; i < con; i++) {
394: while (prmap[i] >= owners[proc + 1]) proc++;
395: len_si[proc]++; /* num of rows in Co(=Pt*AP) to be sent to [proc] */
396: len_s[proc] += coi[i + 1] - coi[i]; /* num of nonzeros in Co to be sent to [proc] */
397: }
399: len = 0; /* max length of buf_si[], see (4) */
400: owners_co[0] = 0;
401: nsend = 0;
402: for (proc = 0; proc < size; proc++) {
403: owners_co[proc + 1] = owners_co[proc] + len_si[proc];
404: if (len_s[proc]) {
405: nsend++;
406: len_si[proc] = 2 * (len_si[proc] + 1); /* length of buf_si to be sent to [proc] */
407: len += len_si[proc];
408: }
409: }
411: /* determine the number and length of messages to receive for coi and coj */
412: PetscCall(PetscGatherNumberOfMessages(comm, NULL, len_s, &nrecv));
413: PetscCall(PetscGatherMessageLengths2(comm, nsend, nrecv, len_s, len_si, &id_r, &len_r, &len_ri));
415: /* post the Irecv and Isend of coj */
416: PetscCall(PetscCommGetNewTag(comm, &tagj));
417: PetscCall(PetscPostIrecvInt(comm, tagj, nrecv, id_r, len_r, &buf_rj, &rwaits));
418: PetscCall(PetscMalloc1(nsend + 1, &swaits));
419: for (proc = 0, k = 0; proc < size; proc++) {
420: if (!len_s[proc]) continue;
421: i = owners_co[proc];
422: PetscCallMPI(MPIU_Isend(coj + coi[i], len_s[proc], MPIU_INT, proc, tagj, comm, swaits + k));
423: k++;
424: }
426: /* (2-2) compute symbolic C_loc = Rd*AP_loc */
427: PetscCall(MatProductCreate(ptap->Rd, ptap->AP_loc, NULL, &ptap->C_loc));
428: PetscCall(MatProductSetType(ptap->C_loc, MATPRODUCT_AB));
429: PetscCall(MatProductSetAlgorithm(ptap->C_loc, "default"));
430: PetscCall(MatProductSetFill(ptap->C_loc, fill));
432: PetscCall(MatSetOptionsPrefix(ptap->C_loc, prefix));
433: PetscCall(MatAppendOptionsPrefix(ptap->C_loc, "inner_diag_"));
435: PetscCall(MatProductSetFromOptions(ptap->C_loc));
436: PetscCall(MatProductSymbolic(ptap->C_loc));
438: c_loc = (Mat_SeqAIJ *)ptap->C_loc->data;
439: PetscCall(ISLocalToGlobalMappingApply(ptap->ltog, c_loc->i[ptap->C_loc->rmap->n], c_loc->j, c_loc->j));
441: /* receives coj are complete */
442: for (i = 0; i < nrecv; i++) PetscCallMPI(MPI_Waitany(nrecv, rwaits, &icompleted, &rstatus));
443: PetscCall(PetscFree(rwaits));
444: if (nsend) PetscCallMPI(MPI_Waitall(nsend, swaits, sstatus));
446: /* add received column indices into ta to update Crmax */
447: for (k = 0; k < nrecv; k++) { /* k-th received message */
448: Jptr = buf_rj[k];
449: for (j = 0; j < len_r[k]; j++) PetscCall(PetscHMapISet(ta, *(Jptr + j) + 1, 1));
450: }
451: PetscCall(PetscHMapIGetSize(ta, &Crmax));
452: PetscCall(PetscHMapIDestroy(&ta));
454: /* (4) send and recv coi */
455: PetscCall(PetscCommGetNewTag(comm, &tagi));
456: PetscCall(PetscPostIrecvInt(comm, tagi, nrecv, id_r, len_ri, &buf_ri, &rwaits));
457: PetscCall(PetscMalloc1(len + 1, &buf_s));
458: buf_si = buf_s; /* points to the beginning of k-th msg to be sent */
459: for (proc = 0, k = 0; proc < size; proc++) {
460: if (!len_s[proc]) continue;
461: /* form outgoing message for i-structure:
462: buf_si[0]: nrows to be sent
463: [1:nrows]: row index (global)
464: [nrows+1:2*nrows+1]: i-structure index
465: */
466: nrows = len_si[proc] / 2 - 1; /* num of rows in Co to be sent to [proc] */
467: buf_si_i = buf_si + nrows + 1;
468: buf_si[0] = nrows;
469: buf_si_i[0] = 0;
470: nrows = 0;
471: for (i = owners_co[proc]; i < owners_co[proc + 1]; i++) {
472: nzi = coi[i + 1] - coi[i];
473: buf_si_i[nrows + 1] = buf_si_i[nrows] + nzi; /* i-structure */
474: buf_si[nrows + 1] = prmap[i] - owners[proc]; /* local row index */
475: nrows++;
476: }
477: PetscCallMPI(MPIU_Isend(buf_si, len_si[proc], MPIU_INT, proc, tagi, comm, swaits + k));
478: k++;
479: buf_si += len_si[proc];
480: }
481: for (i = 0; i < nrecv; i++) PetscCallMPI(MPI_Waitany(nrecv, rwaits, &icompleted, &rstatus));
482: PetscCall(PetscFree(rwaits));
483: if (nsend) PetscCallMPI(MPI_Waitall(nsend, swaits, sstatus));
485: PetscCall(PetscFree4(len_s, len_si, sstatus, owners_co));
486: PetscCall(PetscFree(len_ri));
487: PetscCall(PetscFree(swaits));
488: PetscCall(PetscFree(buf_s));
490: /* (5) compute the local portion of Cmpi */
491: /* set initial free space to be Crmax, sufficient for holding nonzeros in each row of Cmpi */
492: PetscCall(PetscFreeSpaceGet(Crmax, &free_space));
493: current_space = free_space;
495: PetscCall(PetscMalloc3(nrecv, &buf_ri_k, nrecv, &nextrow, nrecv, &nextci));
496: for (k = 0; k < nrecv; k++) {
497: buf_ri_k[k] = buf_ri[k]; /* beginning of k-th recved i-structure */
498: nrows = *buf_ri_k[k];
499: nextrow[k] = buf_ri_k[k] + 1; /* next row number of k-th recved i-structure */
500: nextci[k] = buf_ri_k[k] + (nrows + 1); /* points to the next i-structure of k-th recved i-structure */
501: }
503: MatPreallocateBegin(comm, pn, pn, dnz, onz);
504: PetscCall(PetscLLCondensedCreate_Scalable(Crmax, &lnk));
505: for (i = 0; i < pn; i++) {
506: /* add C_loc into Cmpi */
507: nzi = c_loc->i[i + 1] - c_loc->i[i];
508: Jptr = c_loc->j + c_loc->i[i];
509: PetscCall(PetscLLCondensedAddSorted_Scalable(nzi, Jptr, lnk));
511: /* add received col data into lnk */
512: for (k = 0; k < nrecv; k++) { /* k-th received message */
513: if (i == *nextrow[k]) { /* i-th row */
514: nzi = *(nextci[k] + 1) - *nextci[k];
515: Jptr = buf_rj[k] + *nextci[k];
516: PetscCall(PetscLLCondensedAddSorted_Scalable(nzi, Jptr, lnk));
517: nextrow[k]++;
518: nextci[k]++;
519: }
520: }
521: nzi = lnk[0];
523: /* copy data into free space, then initialize lnk */
524: PetscCall(PetscLLCondensedClean_Scalable(nzi, current_space->array, lnk));
525: PetscCall(MatPreallocateSet(i + owners[rank], nzi, current_space->array, dnz, onz));
526: }
527: PetscCall(PetscFree3(buf_ri_k, nextrow, nextci));
528: PetscCall(PetscLLCondensedDestroy_Scalable(lnk));
529: PetscCall(PetscFreeSpaceDestroy(free_space));
531: /* local sizes and preallocation */
532: PetscCall(MatSetSizes(Cmpi, pn, pn, PETSC_DETERMINE, PETSC_DETERMINE));
533: if (P->cmap->bs > 0) {
534: PetscCall(PetscLayoutSetBlockSize(Cmpi->rmap, P->cmap->bs));
535: PetscCall(PetscLayoutSetBlockSize(Cmpi->cmap, P->cmap->bs));
536: }
537: PetscCall(MatMPIAIJSetPreallocation(Cmpi, 0, dnz, 0, onz));
538: MatPreallocateEnd(dnz, onz);
540: /* members in merge */
541: PetscCall(PetscFree(id_r));
542: PetscCall(PetscFree(len_r));
543: PetscCall(PetscFree(buf_ri[0]));
544: PetscCall(PetscFree(buf_ri));
545: PetscCall(PetscFree(buf_rj[0]));
546: PetscCall(PetscFree(buf_rj));
547: PetscCall(PetscLayoutDestroy(&rowmap));
549: nout = 0;
550: PetscCall(ISGlobalToLocalMappingApply(ptap->ltog, IS_GTOLM_DROP, c_oth->i[ptap->C_oth->rmap->n], c_oth->j, &nout, c_oth->j));
551: PetscCheck(c_oth->i[ptap->C_oth->rmap->n] == nout, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Incorrect mapping %" PetscInt_FMT " != %" PetscInt_FMT, c_oth->i[ptap->C_oth->rmap->n], nout);
552: PetscCall(ISGlobalToLocalMappingApply(ptap->ltog, IS_GTOLM_DROP, c_loc->i[ptap->C_loc->rmap->n], c_loc->j, &nout, c_loc->j));
553: PetscCheck(c_loc->i[ptap->C_loc->rmap->n] == nout, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Incorrect mapping %" PetscInt_FMT " != %" PetscInt_FMT, c_loc->i[ptap->C_loc->rmap->n], nout);
555: /* attach the supporting struct to Cmpi for reuse */
556: Cmpi->product->data = ptap;
557: Cmpi->product->view = MatView_MPIAIJ_PtAP;
558: Cmpi->product->destroy = MatProductCtxDestroy_MPIAIJ_PtAP;
560: /* Cmpi is not ready for use - assembly will be done by MatPtAPNumeric() */
561: Cmpi->assembled = PETSC_FALSE;
562: Cmpi->ops->ptapnumeric = MatPtAPNumeric_MPIAIJ_MPIAIJ_scalable;
563: PetscFunctionReturn(PETSC_SUCCESS);
564: }
566: static inline PetscErrorCode MatPtAPSymbolicComputeOneRowOfAP_private(Mat A, Mat P, Mat P_oth, const PetscInt *map, PetscInt dof, PetscInt i, PetscHSetI dht, PetscHSetI oht)
567: {
568: Mat_MPIAIJ *a = (Mat_MPIAIJ *)A->data, *p = (Mat_MPIAIJ *)P->data;
569: 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;
570: PetscInt *ai, nzi, j, *aj, row, col, *pi, *pj, pnz, nzpi, *p_othcols, k;
571: PetscInt pcstart, pcend, column, offset;
573: PetscFunctionBegin;
574: pcstart = P->cmap->rstart;
575: pcstart *= dof;
576: pcend = P->cmap->rend;
577: pcend *= dof;
578: /* diagonal portion: Ad[i,:]*P */
579: ai = ad->i;
580: nzi = ai[i + 1] - ai[i];
581: aj = ad->j + ai[i];
582: for (j = 0; j < nzi; j++) {
583: row = aj[j];
584: offset = row % dof;
585: row /= dof;
586: nzpi = pd->i[row + 1] - pd->i[row];
587: pj = pd->j + pd->i[row];
588: for (k = 0; k < nzpi; k++) PetscCall(PetscHSetIAdd(dht, pj[k] * dof + offset + pcstart));
589: }
590: /* off-diagonal P */
591: for (j = 0; j < nzi; j++) {
592: row = aj[j];
593: offset = row % dof;
594: row /= dof;
595: nzpi = po->i[row + 1] - po->i[row];
596: pj = PetscSafePointerPlusOffset(po->j, po->i[row]);
597: for (k = 0; k < nzpi; k++) PetscCall(PetscHSetIAdd(oht, p->garray[pj[k]] * dof + offset));
598: }
600: /* off-diagonal part: Ao[i, :]*P_oth */
601: if (ao) {
602: ai = ao->i;
603: pi = p_oth->i;
604: nzi = ai[i + 1] - ai[i];
605: aj = ao->j + ai[i];
606: for (j = 0; j < nzi; j++) {
607: row = aj[j];
608: offset = a->garray[row] % dof;
609: row = map[row];
610: pnz = pi[row + 1] - pi[row];
611: p_othcols = p_oth->j + pi[row];
612: for (col = 0; col < pnz; col++) {
613: column = p_othcols[col] * dof + offset;
614: if (column >= pcstart && column < pcend) {
615: PetscCall(PetscHSetIAdd(dht, column));
616: } else {
617: PetscCall(PetscHSetIAdd(oht, column));
618: }
619: }
620: }
621: } /* end if (ao) */
622: PetscFunctionReturn(PETSC_SUCCESS);
623: }
625: static inline PetscErrorCode MatPtAPNumericComputeOneRowOfAP_private(Mat A, Mat P, Mat P_oth, const PetscInt *map, PetscInt dof, PetscInt i, PetscHMapIV hmap)
626: {
627: Mat_MPIAIJ *a = (Mat_MPIAIJ *)A->data, *p = (Mat_MPIAIJ *)P->data;
628: 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;
629: PetscInt *ai, nzi, j, *aj, row, col, *pi, pnz, *p_othcols, pcstart, *pj, k, nzpi, offset;
630: PetscScalar ra, *aa, *pa;
632: PetscFunctionBegin;
633: pcstart = P->cmap->rstart;
634: pcstart *= dof;
636: /* diagonal portion: Ad[i,:]*P */
637: ai = ad->i;
638: nzi = ai[i + 1] - ai[i];
639: aj = ad->j + ai[i];
640: aa = ad->a + ai[i];
641: for (j = 0; j < nzi; j++) {
642: ra = aa[j];
643: row = aj[j];
644: offset = row % dof;
645: row /= dof;
646: nzpi = pd->i[row + 1] - pd->i[row];
647: pj = pd->j + pd->i[row];
648: pa = pd->a + pd->i[row];
649: for (k = 0; k < nzpi; k++) PetscCall(PetscHMapIVAddValue(hmap, pj[k] * dof + offset + pcstart, ra * pa[k]));
650: PetscCall(PetscLogFlops(2.0 * nzpi));
651: }
652: for (j = 0; j < nzi; j++) {
653: ra = aa[j];
654: row = aj[j];
655: offset = row % dof;
656: row /= dof;
657: nzpi = po->i[row + 1] - po->i[row];
658: pj = PetscSafePointerPlusOffset(po->j, po->i[row]);
659: pa = PetscSafePointerPlusOffset(po->a, po->i[row]);
660: for (k = 0; k < nzpi; k++) PetscCall(PetscHMapIVAddValue(hmap, p->garray[pj[k]] * dof + offset, ra * pa[k]));
661: PetscCall(PetscLogFlops(2.0 * nzpi));
662: }
664: /* off-diagonal part: Ao[i, :]*P_oth */
665: if (ao) {
666: ai = ao->i;
667: pi = p_oth->i;
668: nzi = ai[i + 1] - ai[i];
669: aj = ao->j + ai[i];
670: aa = ao->a + ai[i];
671: for (j = 0; j < nzi; j++) {
672: row = aj[j];
673: offset = a->garray[row] % dof;
674: row = map[row];
675: ra = aa[j];
676: pnz = pi[row + 1] - pi[row];
677: p_othcols = p_oth->j + pi[row];
678: pa = p_oth->a + pi[row];
679: for (col = 0; col < pnz; col++) PetscCall(PetscHMapIVAddValue(hmap, p_othcols[col] * dof + offset, ra * pa[col]));
680: PetscCall(PetscLogFlops(2.0 * pnz));
681: }
682: } /* end if (ao) */
683: PetscFunctionReturn(PETSC_SUCCESS);
684: }
686: PetscErrorCode MatGetBrowsOfAcols_MPIXAIJ(Mat, Mat, PetscInt dof, MatReuse, Mat *);
688: PetscErrorCode MatPtAPNumeric_MPIAIJ_MPIXAIJ_allatonce(Mat A, Mat P, PetscInt dof, Mat C)
689: {
690: Mat_MPIAIJ *p = (Mat_MPIAIJ *)P->data, *c = (Mat_MPIAIJ *)C->data;
691: Mat_SeqAIJ *cd, *co, *po = (Mat_SeqAIJ *)p->B->data, *pd = (Mat_SeqAIJ *)p->A->data;
692: MatProductCtx_APMPI *ptap;
693: PetscHMapIV hmap;
694: 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;
695: PetscScalar *c_rmta, *c_otha, *poa, *pda, *apvalues, *apvaluestmp, *c_rmtaa;
696: PetscInt offset, ii, pocol;
697: const PetscInt *mappingindices;
698: IS map;
700: PetscFunctionBegin;
701: MatCheckProduct(C, 4);
702: ptap = (MatProductCtx_APMPI *)C->product->data;
703: PetscCheck(ptap, PetscObjectComm((PetscObject)C), PETSC_ERR_ARG_WRONGSTATE, "PtAP cannot be computed. Missing data");
704: PetscCheck(ptap->P_oth, PetscObjectComm((PetscObject)C), PETSC_ERR_ARG_WRONGSTATE, "PtAP cannot be reused. Do not call MatProductClear()");
706: PetscCall(MatZeroEntries(C));
708: /* Get P_oth = ptap->P_oth and P_loc = ptap->P_loc */
709: if (ptap->reuse == MAT_REUSE_MATRIX) {
710: /* P_oth and P_loc are obtained in MatPtASymbolic() when reuse == MAT_INITIAL_MATRIX */
711: PetscCall(MatGetBrowsOfAcols_MPIXAIJ(A, P, dof, MAT_REUSE_MATRIX, &ptap->P_oth));
712: }
713: PetscCall(PetscObjectQuery((PetscObject)ptap->P_oth, "aoffdiagtopothmapping", (PetscObject *)&map));
715: PetscCall(MatGetLocalSize(p->B, NULL, &pon));
716: pon *= dof;
717: PetscCall(PetscCalloc2(ptap->c_rmti[pon], &c_rmtj, ptap->c_rmti[pon], &c_rmta));
718: PetscCall(MatGetLocalSize(A, &am, NULL));
719: cmaxr = 0;
720: for (i = 0; i < pon; i++) cmaxr = PetscMax(cmaxr, ptap->c_rmti[i + 1] - ptap->c_rmti[i]);
721: PetscCall(PetscCalloc4(cmaxr, &apindices, cmaxr, &apvalues, cmaxr, &apvaluestmp, pon, &c_rmtc));
722: PetscCall(PetscHMapIVCreateWithSize(cmaxr, &hmap));
723: PetscCall(ISGetIndices(map, &mappingindices));
724: for (i = 0; i < am && pon; i++) {
725: PetscCall(PetscHMapIVClear(hmap));
726: offset = i % dof;
727: ii = i / dof;
728: nzi = po->i[ii + 1] - po->i[ii];
729: if (!nzi) continue;
730: PetscCall(MatPtAPNumericComputeOneRowOfAP_private(A, P, ptap->P_oth, mappingindices, dof, i, hmap));
731: voff = 0;
732: PetscCall(PetscHMapIVGetPairs(hmap, &voff, apindices, apvalues));
733: if (!voff) continue;
735: /* Form C(ii, :) */
736: poj = po->j + po->i[ii];
737: poa = po->a + po->i[ii];
738: for (j = 0; j < nzi; j++) {
739: pocol = poj[j] * dof + offset;
740: c_rmtjj = c_rmtj + ptap->c_rmti[pocol];
741: c_rmtaa = c_rmta + ptap->c_rmti[pocol];
742: for (jj = 0; jj < voff; jj++) {
743: apvaluestmp[jj] = apvalues[jj] * poa[j];
744: /* If the row is empty */
745: if (!c_rmtc[pocol]) {
746: c_rmtjj[jj] = apindices[jj];
747: c_rmtaa[jj] = apvaluestmp[jj];
748: c_rmtc[pocol]++;
749: } else {
750: PetscCall(PetscFindInt(apindices[jj], c_rmtc[pocol], c_rmtjj, &loc));
751: if (loc >= 0) { /* hit */
752: c_rmtaa[loc] += apvaluestmp[jj];
753: PetscCall(PetscLogFlops(1.0));
754: } else { /* new element */
755: loc = -(loc + 1);
756: /* Move data backward */
757: for (kk = c_rmtc[pocol]; kk > loc; kk--) {
758: c_rmtjj[kk] = c_rmtjj[kk - 1];
759: c_rmtaa[kk] = c_rmtaa[kk - 1];
760: } /* End kk */
761: c_rmtjj[loc] = apindices[jj];
762: c_rmtaa[loc] = apvaluestmp[jj];
763: c_rmtc[pocol]++;
764: }
765: }
766: PetscCall(PetscLogFlops(voff));
767: } /* End jj */
768: } /* End j */
769: } /* End i */
771: PetscCall(PetscFree4(apindices, apvalues, apvaluestmp, c_rmtc));
772: PetscCall(PetscHMapIVDestroy(&hmap));
774: PetscCall(MatGetLocalSize(P, NULL, &pn));
775: pn *= dof;
776: PetscCall(PetscCalloc2(ptap->c_othi[pn], &c_othj, ptap->c_othi[pn], &c_otha));
778: PetscCall(PetscSFReduceBegin(ptap->sf, MPIU_INT, c_rmtj, c_othj, MPI_REPLACE));
779: PetscCall(PetscSFReduceBegin(ptap->sf, MPIU_SCALAR, c_rmta, c_otha, MPI_REPLACE));
780: PetscCall(MatGetOwnershipRangeColumn(P, &pcstart, &pcend));
781: pcstart = pcstart * dof;
782: pcend = pcend * dof;
783: cd = (Mat_SeqAIJ *)c->A->data;
784: co = (Mat_SeqAIJ *)c->B->data;
786: cmaxr = 0;
787: for (i = 0; i < pn; i++) cmaxr = PetscMax(cmaxr, (cd->i[i + 1] - cd->i[i]) + (co->i[i + 1] - co->i[i]));
788: PetscCall(PetscCalloc5(cmaxr, &apindices, cmaxr, &apvalues, cmaxr, &apvaluestmp, pn, &dcc, pn, &occ));
789: PetscCall(PetscHMapIVCreateWithSize(cmaxr, &hmap));
790: for (i = 0; i < am && pn; i++) {
791: PetscCall(PetscHMapIVClear(hmap));
792: offset = i % dof;
793: ii = i / dof;
794: nzi = pd->i[ii + 1] - pd->i[ii];
795: if (!nzi) continue;
796: PetscCall(MatPtAPNumericComputeOneRowOfAP_private(A, P, ptap->P_oth, mappingindices, dof, i, hmap));
797: voff = 0;
798: PetscCall(PetscHMapIVGetPairs(hmap, &voff, apindices, apvalues));
799: if (!voff) continue;
800: /* Form C(ii, :) */
801: pdj = pd->j + pd->i[ii];
802: pda = pd->a + pd->i[ii];
803: for (j = 0; j < nzi; j++) {
804: row = pcstart + pdj[j] * dof + offset;
805: for (jj = 0; jj < voff; jj++) apvaluestmp[jj] = apvalues[jj] * pda[j];
806: PetscCall(PetscLogFlops(voff));
807: PetscCall(MatSetValues(C, 1, &row, voff, apindices, apvaluestmp, ADD_VALUES));
808: }
809: }
810: PetscCall(ISRestoreIndices(map, &mappingindices));
811: PetscCall(MatGetOwnershipRangeColumn(C, &ccstart, &ccend));
812: PetscCall(PetscFree5(apindices, apvalues, apvaluestmp, dcc, occ));
813: PetscCall(PetscHMapIVDestroy(&hmap));
814: PetscCall(PetscSFReduceEnd(ptap->sf, MPIU_INT, c_rmtj, c_othj, MPI_REPLACE));
815: PetscCall(PetscSFReduceEnd(ptap->sf, MPIU_SCALAR, c_rmta, c_otha, MPI_REPLACE));
816: PetscCall(PetscFree2(c_rmtj, c_rmta));
818: /* Add contributions from remote */
819: for (i = 0; i < pn; i++) {
820: row = i + pcstart;
821: PetscCall(MatSetValues(C, 1, &row, ptap->c_othi[i + 1] - ptap->c_othi[i], PetscSafePointerPlusOffset(c_othj, ptap->c_othi[i]), PetscSafePointerPlusOffset(c_otha, ptap->c_othi[i]), ADD_VALUES));
822: }
823: PetscCall(PetscFree2(c_othj, c_otha));
825: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
826: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
828: ptap->reuse = MAT_REUSE_MATRIX;
829: PetscFunctionReturn(PETSC_SUCCESS);
830: }
832: PetscErrorCode MatPtAPNumeric_MPIAIJ_MPIAIJ_allatonce(Mat A, Mat P, Mat C)
833: {
834: PetscFunctionBegin;
835: PetscCall(MatPtAPNumeric_MPIAIJ_MPIXAIJ_allatonce(A, P, 1, C));
836: PetscFunctionReturn(PETSC_SUCCESS);
837: }
839: PetscErrorCode MatPtAPNumeric_MPIAIJ_MPIXAIJ_allatonce_merged(Mat A, Mat P, PetscInt dof, Mat C)
840: {
841: Mat_MPIAIJ *p = (Mat_MPIAIJ *)P->data, *c = (Mat_MPIAIJ *)C->data;
842: Mat_SeqAIJ *cd, *co, *po = (Mat_SeqAIJ *)p->B->data, *pd = (Mat_SeqAIJ *)p->A->data;
843: MatProductCtx_APMPI *ptap;
844: PetscHMapIV hmap;
845: 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;
846: PetscScalar *c_rmta, *c_otha, *poa, *pda, *apvalues, *apvaluestmp, *c_rmtaa;
847: PetscInt offset, ii, pocol;
848: const PetscInt *mappingindices;
849: IS map;
851: PetscFunctionBegin;
852: MatCheckProduct(C, 4);
853: ptap = (MatProductCtx_APMPI *)C->product->data;
854: PetscCheck(ptap, PetscObjectComm((PetscObject)C), PETSC_ERR_ARG_WRONGSTATE, "PtAP cannot be computed. Missing data");
855: PetscCheck(ptap->P_oth, PetscObjectComm((PetscObject)C), PETSC_ERR_ARG_WRONGSTATE, "PtAP cannot be reused. Do not call MatProductClear()");
857: PetscCall(MatZeroEntries(C));
859: /* Get P_oth = ptap->P_oth and P_loc = ptap->P_loc */
860: if (ptap->reuse == MAT_REUSE_MATRIX) {
861: /* P_oth and P_loc are obtained in MatPtASymbolic() when reuse == MAT_INITIAL_MATRIX */
862: PetscCall(MatGetBrowsOfAcols_MPIXAIJ(A, P, dof, MAT_REUSE_MATRIX, &ptap->P_oth));
863: }
864: PetscCall(PetscObjectQuery((PetscObject)ptap->P_oth, "aoffdiagtopothmapping", (PetscObject *)&map));
865: PetscCall(MatGetLocalSize(p->B, NULL, &pon));
866: pon *= dof;
867: PetscCall(MatGetLocalSize(P, NULL, &pn));
868: pn *= dof;
870: PetscCall(PetscCalloc2(ptap->c_rmti[pon], &c_rmtj, ptap->c_rmti[pon], &c_rmta));
871: PetscCall(MatGetLocalSize(A, &am, NULL));
872: PetscCall(MatGetOwnershipRangeColumn(P, &pcstart, &pcend));
873: pcstart *= dof;
874: pcend *= dof;
875: cmaxr = 0;
876: for (i = 0; i < pon; i++) cmaxr = PetscMax(cmaxr, ptap->c_rmti[i + 1] - ptap->c_rmti[i]);
877: cd = (Mat_SeqAIJ *)c->A->data;
878: co = (Mat_SeqAIJ *)c->B->data;
879: for (i = 0; i < pn; i++) cmaxr = PetscMax(cmaxr, (cd->i[i + 1] - cd->i[i]) + (co->i[i + 1] - co->i[i]));
880: PetscCall(PetscCalloc4(cmaxr, &apindices, cmaxr, &apvalues, cmaxr, &apvaluestmp, pon, &c_rmtc));
881: PetscCall(PetscHMapIVCreateWithSize(cmaxr, &hmap));
882: PetscCall(ISGetIndices(map, &mappingindices));
883: for (i = 0; i < am && (pon || pn); i++) {
884: PetscCall(PetscHMapIVClear(hmap));
885: offset = i % dof;
886: ii = i / dof;
887: nzi = po->i[ii + 1] - po->i[ii];
888: dnzi = pd->i[ii + 1] - pd->i[ii];
889: if (!nzi && !dnzi) continue;
890: PetscCall(MatPtAPNumericComputeOneRowOfAP_private(A, P, ptap->P_oth, mappingindices, dof, i, hmap));
891: voff = 0;
892: PetscCall(PetscHMapIVGetPairs(hmap, &voff, apindices, apvalues));
893: if (!voff) continue;
895: /* Form remote C(ii, :) */
896: poj = PetscSafePointerPlusOffset(po->j, po->i[ii]);
897: poa = PetscSafePointerPlusOffset(po->a, po->i[ii]);
898: for (j = 0; j < nzi; j++) {
899: pocol = poj[j] * dof + offset;
900: c_rmtjj = c_rmtj + ptap->c_rmti[pocol];
901: c_rmtaa = c_rmta + ptap->c_rmti[pocol];
902: for (jj = 0; jj < voff; jj++) {
903: apvaluestmp[jj] = apvalues[jj] * poa[j];
904: /* If the row is empty */
905: if (!c_rmtc[pocol]) {
906: c_rmtjj[jj] = apindices[jj];
907: c_rmtaa[jj] = apvaluestmp[jj];
908: c_rmtc[pocol]++;
909: } else {
910: PetscCall(PetscFindInt(apindices[jj], c_rmtc[pocol], c_rmtjj, &loc));
911: if (loc >= 0) { /* hit */
912: c_rmtaa[loc] += apvaluestmp[jj];
913: PetscCall(PetscLogFlops(1.0));
914: } else { /* new element */
915: loc = -(loc + 1);
916: /* Move data backward */
917: for (kk = c_rmtc[pocol]; kk > loc; kk--) {
918: c_rmtjj[kk] = c_rmtjj[kk - 1];
919: c_rmtaa[kk] = c_rmtaa[kk - 1];
920: } /* End kk */
921: c_rmtjj[loc] = apindices[jj];
922: c_rmtaa[loc] = apvaluestmp[jj];
923: c_rmtc[pocol]++;
924: }
925: }
926: } /* End jj */
927: PetscCall(PetscLogFlops(voff));
928: } /* End j */
930: /* Form local C(ii, :) */
931: pdj = pd->j + pd->i[ii];
932: pda = pd->a + pd->i[ii];
933: for (j = 0; j < dnzi; j++) {
934: row = pcstart + pdj[j] * dof + offset;
935: for (jj = 0; jj < voff; jj++) apvaluestmp[jj] = apvalues[jj] * pda[j]; /* End kk */
936: PetscCall(PetscLogFlops(voff));
937: PetscCall(MatSetValues(C, 1, &row, voff, apindices, apvaluestmp, ADD_VALUES));
938: } /* End j */
939: } /* End i */
941: PetscCall(ISRestoreIndices(map, &mappingindices));
942: PetscCall(PetscFree4(apindices, apvalues, apvaluestmp, c_rmtc));
943: PetscCall(PetscHMapIVDestroy(&hmap));
944: PetscCall(PetscCalloc2(ptap->c_othi[pn], &c_othj, ptap->c_othi[pn], &c_otha));
946: PetscCall(PetscSFReduceBegin(ptap->sf, MPIU_INT, c_rmtj, c_othj, MPI_REPLACE));
947: PetscCall(PetscSFReduceBegin(ptap->sf, MPIU_SCALAR, c_rmta, c_otha, MPI_REPLACE));
948: PetscCall(PetscSFReduceEnd(ptap->sf, MPIU_INT, c_rmtj, c_othj, MPI_REPLACE));
949: PetscCall(PetscSFReduceEnd(ptap->sf, MPIU_SCALAR, c_rmta, c_otha, MPI_REPLACE));
950: PetscCall(PetscFree2(c_rmtj, c_rmta));
952: /* Add contributions from remote */
953: for (i = 0; i < pn; i++) {
954: row = i + pcstart;
955: PetscCall(MatSetValues(C, 1, &row, ptap->c_othi[i + 1] - ptap->c_othi[i], PetscSafePointerPlusOffset(c_othj, ptap->c_othi[i]), PetscSafePointerPlusOffset(c_otha, ptap->c_othi[i]), ADD_VALUES));
956: }
957: PetscCall(PetscFree2(c_othj, c_otha));
959: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
960: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
962: ptap->reuse = MAT_REUSE_MATRIX;
963: PetscFunctionReturn(PETSC_SUCCESS);
964: }
966: PetscErrorCode MatPtAPNumeric_MPIAIJ_MPIAIJ_allatonce_merged(Mat A, Mat P, Mat C)
967: {
968: PetscFunctionBegin;
969: PetscCall(MatPtAPNumeric_MPIAIJ_MPIXAIJ_allatonce_merged(A, P, 1, C));
970: PetscFunctionReturn(PETSC_SUCCESS);
971: }
973: /* TODO: move algorithm selection to MatProductSetFromOptions */
974: PetscErrorCode MatPtAPSymbolic_MPIAIJ_MPIXAIJ_allatonce(Mat A, Mat P, PetscInt dof, PetscReal fill, Mat Cmpi)
975: {
976: MatProductCtx_APMPI *ptap;
977: Mat_MPIAIJ *p = (Mat_MPIAIJ *)P->data;
978: MPI_Comm comm;
979: Mat_SeqAIJ *pd = (Mat_SeqAIJ *)p->A->data, *po = (Mat_SeqAIJ *)p->B->data;
980: MatType mtype;
981: PetscSF sf;
982: PetscSFNode *iremote;
983: PetscInt rootspacesize, *rootspace, *rootspaceoffsets, nleaves;
984: const PetscInt *rootdegrees;
985: PetscHSetI ht, oht, *hta, *hto;
986: PetscInt pn, pon, *c_rmtc, i, j, nzi, htsize, htosize, *c_rmtj, off, *c_othj, rcvncols, sendncols, *c_rmtoffsets;
987: PetscInt lidx, *rdj, col, pcstart, pcend, *dnz, *onz, am, arstart, arend, *poj, *pdj;
988: PetscInt nalg = 2, alg = 0, offset, ii;
989: PetscMPIInt owner;
990: const PetscInt *mappingindices;
991: PetscBool flg;
992: const char *algTypes[2] = {"overlapping", "merged"};
993: IS map;
995: PetscFunctionBegin;
996: MatCheckProduct(Cmpi, 5);
997: PetscCheck(!Cmpi->product->data, PetscObjectComm((PetscObject)Cmpi), PETSC_ERR_PLIB, "Product data not empty");
998: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
1000: /* Create symbolic parallel matrix Cmpi */
1001: PetscCall(MatGetLocalSize(P, NULL, &pn));
1002: pn *= dof;
1003: PetscCall(MatGetType(A, &mtype));
1004: PetscCall(MatSetType(Cmpi, mtype));
1005: PetscCall(MatSetSizes(Cmpi, pn, pn, PETSC_DETERMINE, PETSC_DETERMINE));
1007: PetscCall(PetscNew(&ptap));
1008: ptap->reuse = MAT_INITIAL_MATRIX;
1009: ptap->algType = 2;
1011: /* Get P_oth by taking rows of P (= non-zero cols of local A) from other processors */
1012: PetscCall(MatGetBrowsOfAcols_MPIXAIJ(A, P, dof, MAT_INITIAL_MATRIX, &ptap->P_oth));
1013: PetscCall(PetscObjectQuery((PetscObject)ptap->P_oth, "aoffdiagtopothmapping", (PetscObject *)&map));
1014: /* This equals to the number of offdiag columns in P */
1015: PetscCall(MatGetLocalSize(p->B, NULL, &pon));
1016: pon *= dof;
1017: /* offsets */
1018: PetscCall(PetscMalloc1(pon + 1, &ptap->c_rmti));
1019: /* The number of columns we will send to remote ranks */
1020: PetscCall(PetscMalloc1(pon, &c_rmtc));
1021: PetscCall(PetscMalloc1(pon, &hta));
1022: for (i = 0; i < pon; i++) PetscCall(PetscHSetICreate(&hta[i]));
1023: PetscCall(MatGetLocalSize(A, &am, NULL));
1024: PetscCall(MatGetOwnershipRange(A, &arstart, &arend));
1025: /* Create hash table to merge all columns for C(i, :) */
1026: PetscCall(PetscHSetICreate(&ht));
1028: PetscCall(ISGetIndices(map, &mappingindices));
1029: ptap->c_rmti[0] = 0;
1030: /* 2) Pass 1: calculate the size for C_rmt (a matrix need to be sent to other processors) */
1031: for (i = 0; i < am && pon; i++) {
1032: /* Form one row of AP */
1033: PetscCall(PetscHSetIClear(ht));
1034: offset = i % dof;
1035: ii = i / dof;
1036: /* If the off-diagonal is empty, do not do any calculation */
1037: nzi = po->i[ii + 1] - po->i[ii];
1038: if (!nzi) continue;
1040: PetscCall(MatPtAPSymbolicComputeOneRowOfAP_private(A, P, ptap->P_oth, mappingindices, dof, i, ht, ht));
1041: PetscCall(PetscHSetIGetSize(ht, &htsize));
1042: /* If AP is empty, just continue */
1043: if (!htsize) continue;
1044: /* Form C(ii, :) */
1045: poj = po->j + po->i[ii];
1046: for (j = 0; j < nzi; j++) PetscCall(PetscHSetIUpdate(hta[poj[j] * dof + offset], ht));
1047: }
1049: for (i = 0; i < pon; i++) {
1050: PetscCall(PetscHSetIGetSize(hta[i], &htsize));
1051: ptap->c_rmti[i + 1] = ptap->c_rmti[i] + htsize;
1052: c_rmtc[i] = htsize;
1053: }
1055: PetscCall(PetscMalloc1(ptap->c_rmti[pon], &c_rmtj));
1057: for (i = 0; i < pon; i++) {
1058: off = 0;
1059: PetscCall(PetscHSetIGetElems(hta[i], &off, c_rmtj + ptap->c_rmti[i]));
1060: PetscCall(PetscHSetIDestroy(&hta[i]));
1061: }
1062: PetscCall(PetscFree(hta));
1064: PetscCall(PetscMalloc1(pon, &iremote));
1065: for (i = 0; i < pon; i++) {
1066: owner = 0;
1067: lidx = 0;
1068: offset = i % dof;
1069: ii = i / dof;
1070: PetscCall(PetscLayoutFindOwnerIndex(P->cmap, p->garray[ii], &owner, &lidx));
1071: iremote[i].index = lidx * dof + offset;
1072: iremote[i].rank = owner;
1073: }
1075: PetscCall(PetscSFCreate(comm, &sf));
1076: PetscCall(PetscSFSetGraph(sf, pn, pon, NULL, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER));
1077: /* Reorder ranks properly so that the data handled by gather and scatter have the same order */
1078: PetscCall(PetscSFSetRankOrder(sf, PETSC_TRUE));
1079: PetscCall(PetscSFSetFromOptions(sf));
1080: PetscCall(PetscSFSetUp(sf));
1081: /* How many neighbors have contributions to my rows? */
1082: PetscCall(PetscSFComputeDegreeBegin(sf, &rootdegrees));
1083: PetscCall(PetscSFComputeDegreeEnd(sf, &rootdegrees));
1084: rootspacesize = 0;
1085: for (i = 0; i < pn; i++) rootspacesize += rootdegrees[i];
1086: PetscCall(PetscMalloc1(rootspacesize, &rootspace));
1087: PetscCall(PetscMalloc1(rootspacesize + 1, &rootspaceoffsets));
1088: /* Get information from leaves
1089: * Number of columns other people contribute to my rows
1090: * */
1091: PetscCall(PetscSFGatherBegin(sf, MPIU_INT, c_rmtc, rootspace));
1092: PetscCall(PetscSFGatherEnd(sf, MPIU_INT, c_rmtc, rootspace));
1093: PetscCall(PetscFree(c_rmtc));
1094: PetscCall(PetscCalloc1(pn + 1, &ptap->c_othi));
1095: /* The number of columns is received for each row */
1096: ptap->c_othi[0] = 0;
1097: rootspacesize = 0;
1098: rootspaceoffsets[0] = 0;
1099: for (i = 0; i < pn; i++) {
1100: rcvncols = 0;
1101: for (j = 0; j < rootdegrees[i]; j++) {
1102: rcvncols += rootspace[rootspacesize];
1103: rootspaceoffsets[rootspacesize + 1] = rootspaceoffsets[rootspacesize] + rootspace[rootspacesize];
1104: rootspacesize++;
1105: }
1106: ptap->c_othi[i + 1] = ptap->c_othi[i] + rcvncols;
1107: }
1108: PetscCall(PetscFree(rootspace));
1110: PetscCall(PetscMalloc1(pon, &c_rmtoffsets));
1111: PetscCall(PetscSFScatterBegin(sf, MPIU_INT, rootspaceoffsets, c_rmtoffsets));
1112: PetscCall(PetscSFScatterEnd(sf, MPIU_INT, rootspaceoffsets, c_rmtoffsets));
1113: PetscCall(PetscSFDestroy(&sf));
1114: PetscCall(PetscFree(rootspaceoffsets));
1116: PetscCall(PetscCalloc1(ptap->c_rmti[pon], &iremote));
1117: nleaves = 0;
1118: for (i = 0; i < pon; i++) {
1119: owner = 0;
1120: ii = i / dof;
1121: PetscCall(PetscLayoutFindOwnerIndex(P->cmap, p->garray[ii], &owner, NULL));
1122: sendncols = ptap->c_rmti[i + 1] - ptap->c_rmti[i];
1123: for (j = 0; j < sendncols; j++) {
1124: iremote[nleaves].rank = owner;
1125: iremote[nleaves++].index = c_rmtoffsets[i] + j;
1126: }
1127: }
1128: PetscCall(PetscFree(c_rmtoffsets));
1129: PetscCall(PetscCalloc1(ptap->c_othi[pn], &c_othj));
1131: PetscCall(PetscSFCreate(comm, &ptap->sf));
1132: PetscCall(PetscSFSetGraph(ptap->sf, ptap->c_othi[pn], nleaves, NULL, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER));
1133: PetscCall(PetscSFSetFromOptions(ptap->sf));
1134: /* One to one map */
1135: PetscCall(PetscSFReduceBegin(ptap->sf, MPIU_INT, c_rmtj, c_othj, MPI_REPLACE));
1137: PetscCall(PetscMalloc2(pn, &dnz, pn, &onz));
1138: PetscCall(PetscHSetICreate(&oht));
1139: PetscCall(MatGetOwnershipRangeColumn(P, &pcstart, &pcend));
1140: pcstart *= dof;
1141: pcend *= dof;
1142: PetscCall(PetscMalloc2(pn, &hta, pn, &hto));
1143: for (i = 0; i < pn; i++) {
1144: PetscCall(PetscHSetICreate(&hta[i]));
1145: PetscCall(PetscHSetICreate(&hto[i]));
1146: }
1147: /* Work on local part */
1148: /* 4) Pass 1: Estimate memory for C_loc */
1149: for (i = 0; i < am && pn; i++) {
1150: PetscCall(PetscHSetIClear(ht));
1151: PetscCall(PetscHSetIClear(oht));
1152: offset = i % dof;
1153: ii = i / dof;
1154: nzi = pd->i[ii + 1] - pd->i[ii];
1155: if (!nzi) continue;
1157: PetscCall(MatPtAPSymbolicComputeOneRowOfAP_private(A, P, ptap->P_oth, mappingindices, dof, i, ht, oht));
1158: PetscCall(PetscHSetIGetSize(ht, &htsize));
1159: PetscCall(PetscHSetIGetSize(oht, &htosize));
1160: if (!(htsize + htosize)) continue;
1161: /* Form C(ii, :) */
1162: pdj = pd->j + pd->i[ii];
1163: for (j = 0; j < nzi; j++) {
1164: PetscCall(PetscHSetIUpdate(hta[pdj[j] * dof + offset], ht));
1165: PetscCall(PetscHSetIUpdate(hto[pdj[j] * dof + offset], oht));
1166: }
1167: }
1169: PetscCall(ISRestoreIndices(map, &mappingindices));
1171: PetscCall(PetscHSetIDestroy(&ht));
1172: PetscCall(PetscHSetIDestroy(&oht));
1174: /* Get remote data */
1175: PetscCall(PetscSFReduceEnd(ptap->sf, MPIU_INT, c_rmtj, c_othj, MPI_REPLACE));
1176: PetscCall(PetscFree(c_rmtj));
1178: for (i = 0; i < pn; i++) {
1179: nzi = ptap->c_othi[i + 1] - ptap->c_othi[i];
1180: rdj = PetscSafePointerPlusOffset(c_othj, ptap->c_othi[i]);
1181: for (j = 0; j < nzi; j++) {
1182: col = rdj[j];
1183: /* diag part */
1184: if (col >= pcstart && col < pcend) {
1185: PetscCall(PetscHSetIAdd(hta[i], col));
1186: } else { /* off-diagonal */
1187: PetscCall(PetscHSetIAdd(hto[i], col));
1188: }
1189: }
1190: PetscCall(PetscHSetIGetSize(hta[i], &htsize));
1191: dnz[i] = htsize;
1192: PetscCall(PetscHSetIDestroy(&hta[i]));
1193: PetscCall(PetscHSetIGetSize(hto[i], &htsize));
1194: onz[i] = htsize;
1195: PetscCall(PetscHSetIDestroy(&hto[i]));
1196: }
1198: PetscCall(PetscFree2(hta, hto));
1199: PetscCall(PetscFree(c_othj));
1201: /* local sizes and preallocation */
1202: PetscCall(MatSetSizes(Cmpi, pn, pn, PETSC_DETERMINE, PETSC_DETERMINE));
1203: PetscCall(MatSetBlockSizes(Cmpi, dof > 1 ? dof : P->cmap->bs, dof > 1 ? dof : P->cmap->bs));
1204: PetscCall(MatMPIAIJSetPreallocation(Cmpi, 0, dnz, 0, onz));
1205: PetscCall(MatSetUp(Cmpi));
1206: PetscCall(PetscFree2(dnz, onz));
1208: /* attach the supporting struct to Cmpi for reuse */
1209: Cmpi->product->data = ptap;
1210: Cmpi->product->destroy = MatProductCtxDestroy_MPIAIJ_PtAP;
1211: Cmpi->product->view = MatView_MPIAIJ_PtAP;
1213: /* Cmpi is not ready for use - assembly will be done by MatPtAPNumeric() */
1214: Cmpi->assembled = PETSC_FALSE;
1215: /* pick an algorithm */
1216: PetscOptionsBegin(PetscObjectComm((PetscObject)A), ((PetscObject)A)->prefix, "MatPtAP", "Mat");
1217: alg = 0;
1218: PetscCall(PetscOptionsEList("-matptap_allatonce_via", "PtAP allatonce numeric approach", "MatPtAP", algTypes, nalg, algTypes[alg], &alg, &flg));
1219: PetscOptionsEnd();
1220: switch (alg) {
1221: case 0:
1222: Cmpi->ops->ptapnumeric = MatPtAPNumeric_MPIAIJ_MPIAIJ_allatonce;
1223: break;
1224: case 1:
1225: Cmpi->ops->ptapnumeric = MatPtAPNumeric_MPIAIJ_MPIAIJ_allatonce_merged;
1226: break;
1227: default:
1228: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, " Unsupported allatonce numerical algorithm ");
1229: }
1230: PetscFunctionReturn(PETSC_SUCCESS);
1231: }
1233: PetscErrorCode MatPtAPSymbolic_MPIAIJ_MPIAIJ_allatonce(Mat A, Mat P, PetscReal fill, Mat C)
1234: {
1235: PetscFunctionBegin;
1236: PetscCall(MatPtAPSymbolic_MPIAIJ_MPIXAIJ_allatonce(A, P, 1, fill, C));
1237: PetscFunctionReturn(PETSC_SUCCESS);
1238: }
1240: PetscErrorCode MatPtAPSymbolic_MPIAIJ_MPIXAIJ_allatonce_merged(Mat A, Mat P, PetscInt dof, PetscReal fill, Mat Cmpi)
1241: {
1242: MatProductCtx_APMPI *ptap;
1243: Mat_MPIAIJ *p = (Mat_MPIAIJ *)P->data;
1244: MPI_Comm comm;
1245: Mat_SeqAIJ *pd = (Mat_SeqAIJ *)p->A->data, *po = (Mat_SeqAIJ *)p->B->data;
1246: MatType mtype;
1247: PetscSF sf;
1248: PetscSFNode *iremote;
1249: PetscInt rootspacesize, *rootspace, *rootspaceoffsets, nleaves;
1250: const PetscInt *rootdegrees;
1251: PetscHSetI ht, oht, *hta, *hto, *htd;
1252: PetscInt pn, pon, *c_rmtc, i, j, nzi, dnzi, htsize, htosize, *c_rmtj, off, *c_othj, rcvncols, sendncols, *c_rmtoffsets;
1253: PetscInt lidx, *rdj, col, pcstart, pcend, *dnz, *onz, am, arstart, arend, *poj, *pdj;
1254: PetscInt nalg = 2, alg = 0, offset, ii;
1255: PetscMPIInt owner;
1256: PetscBool flg;
1257: const char *algTypes[2] = {"merged", "overlapping"};
1258: const PetscInt *mappingindices;
1259: IS map;
1261: PetscFunctionBegin;
1262: MatCheckProduct(Cmpi, 5);
1263: PetscCheck(!Cmpi->product->data, PetscObjectComm((PetscObject)Cmpi), PETSC_ERR_PLIB, "Product data not empty");
1264: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
1266: /* Create symbolic parallel matrix Cmpi */
1267: PetscCall(MatGetLocalSize(P, NULL, &pn));
1268: pn *= dof;
1269: PetscCall(MatGetType(A, &mtype));
1270: PetscCall(MatSetType(Cmpi, mtype));
1271: PetscCall(MatSetSizes(Cmpi, pn, pn, PETSC_DETERMINE, PETSC_DETERMINE));
1273: PetscCall(PetscNew(&ptap));
1274: ptap->reuse = MAT_INITIAL_MATRIX;
1275: ptap->algType = 3;
1277: /* 0) Get P_oth by taking rows of P (= non-zero cols of local A) from other processors */
1278: PetscCall(MatGetBrowsOfAcols_MPIXAIJ(A, P, dof, MAT_INITIAL_MATRIX, &ptap->P_oth));
1279: PetscCall(PetscObjectQuery((PetscObject)ptap->P_oth, "aoffdiagtopothmapping", (PetscObject *)&map));
1281: /* This equals to the number of offdiag columns in P */
1282: PetscCall(MatGetLocalSize(p->B, NULL, &pon));
1283: pon *= dof;
1284: /* offsets */
1285: PetscCall(PetscMalloc1(pon + 1, &ptap->c_rmti));
1286: /* The number of columns we will send to remote ranks */
1287: PetscCall(PetscMalloc1(pon, &c_rmtc));
1288: PetscCall(PetscMalloc1(pon, &hta));
1289: for (i = 0; i < pon; i++) PetscCall(PetscHSetICreate(&hta[i]));
1290: PetscCall(MatGetLocalSize(A, &am, NULL));
1291: PetscCall(MatGetOwnershipRange(A, &arstart, &arend));
1292: /* Create hash table to merge all columns for C(i, :) */
1293: PetscCall(PetscHSetICreate(&ht));
1294: PetscCall(PetscHSetICreate(&oht));
1295: PetscCall(PetscMalloc2(pn, &htd, pn, &hto));
1296: for (i = 0; i < pn; i++) {
1297: PetscCall(PetscHSetICreate(&htd[i]));
1298: PetscCall(PetscHSetICreate(&hto[i]));
1299: }
1301: PetscCall(ISGetIndices(map, &mappingindices));
1302: ptap->c_rmti[0] = 0;
1303: /* 2) Pass 1: calculate the size for C_rmt (a matrix need to be sent to other processors) */
1304: for (i = 0; i < am && (pon || pn); i++) {
1305: /* Form one row of AP */
1306: PetscCall(PetscHSetIClear(ht));
1307: PetscCall(PetscHSetIClear(oht));
1308: offset = i % dof;
1309: ii = i / dof;
1310: /* If the off-diagonal is empty, do not do any calculation */
1311: nzi = po->i[ii + 1] - po->i[ii];
1312: dnzi = pd->i[ii + 1] - pd->i[ii];
1313: if (!nzi && !dnzi) continue;
1315: PetscCall(MatPtAPSymbolicComputeOneRowOfAP_private(A, P, ptap->P_oth, mappingindices, dof, i, ht, oht));
1316: PetscCall(PetscHSetIGetSize(ht, &htsize));
1317: PetscCall(PetscHSetIGetSize(oht, &htosize));
1318: /* If AP is empty, just continue */
1319: if (!(htsize + htosize)) continue;
1321: /* Form remote C(ii, :) */
1322: poj = PetscSafePointerPlusOffset(po->j, po->i[ii]);
1323: for (j = 0; j < nzi; j++) {
1324: PetscCall(PetscHSetIUpdate(hta[poj[j] * dof + offset], ht));
1325: PetscCall(PetscHSetIUpdate(hta[poj[j] * dof + offset], oht));
1326: }
1328: /* Form local C(ii, :) */
1329: pdj = pd->j + pd->i[ii];
1330: for (j = 0; j < dnzi; j++) {
1331: PetscCall(PetscHSetIUpdate(htd[pdj[j] * dof + offset], ht));
1332: PetscCall(PetscHSetIUpdate(hto[pdj[j] * dof + offset], oht));
1333: }
1334: }
1336: PetscCall(ISRestoreIndices(map, &mappingindices));
1338: PetscCall(PetscHSetIDestroy(&ht));
1339: PetscCall(PetscHSetIDestroy(&oht));
1341: for (i = 0; i < pon; i++) {
1342: PetscCall(PetscHSetIGetSize(hta[i], &htsize));
1343: ptap->c_rmti[i + 1] = ptap->c_rmti[i] + htsize;
1344: c_rmtc[i] = htsize;
1345: }
1347: PetscCall(PetscMalloc1(ptap->c_rmti[pon], &c_rmtj));
1349: for (i = 0; i < pon; i++) {
1350: off = 0;
1351: PetscCall(PetscHSetIGetElems(hta[i], &off, c_rmtj + ptap->c_rmti[i]));
1352: PetscCall(PetscHSetIDestroy(&hta[i]));
1353: }
1354: PetscCall(PetscFree(hta));
1356: PetscCall(PetscMalloc1(pon, &iremote));
1357: for (i = 0; i < pon; i++) {
1358: owner = 0;
1359: lidx = 0;
1360: offset = i % dof;
1361: ii = i / dof;
1362: PetscCall(PetscLayoutFindOwnerIndex(P->cmap, p->garray[ii], &owner, &lidx));
1363: iremote[i].index = lidx * dof + offset;
1364: iremote[i].rank = owner;
1365: }
1367: PetscCall(PetscSFCreate(comm, &sf));
1368: PetscCall(PetscSFSetGraph(sf, pn, pon, NULL, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER));
1369: /* Reorder ranks properly so that the data handled by gather and scatter have the same order */
1370: PetscCall(PetscSFSetRankOrder(sf, PETSC_TRUE));
1371: PetscCall(PetscSFSetFromOptions(sf));
1372: PetscCall(PetscSFSetUp(sf));
1373: /* How many neighbors have contributions to my rows? */
1374: PetscCall(PetscSFComputeDegreeBegin(sf, &rootdegrees));
1375: PetscCall(PetscSFComputeDegreeEnd(sf, &rootdegrees));
1376: rootspacesize = 0;
1377: for (i = 0; i < pn; i++) rootspacesize += rootdegrees[i];
1378: PetscCall(PetscMalloc1(rootspacesize, &rootspace));
1379: PetscCall(PetscMalloc1(rootspacesize + 1, &rootspaceoffsets));
1380: /* Get information from leaves
1381: * Number of columns other people contribute to my rows
1382: * */
1383: PetscCall(PetscSFGatherBegin(sf, MPIU_INT, c_rmtc, rootspace));
1384: PetscCall(PetscSFGatherEnd(sf, MPIU_INT, c_rmtc, rootspace));
1385: PetscCall(PetscFree(c_rmtc));
1386: PetscCall(PetscMalloc1(pn + 1, &ptap->c_othi));
1387: /* The number of columns is received for each row */
1388: ptap->c_othi[0] = 0;
1389: rootspacesize = 0;
1390: rootspaceoffsets[0] = 0;
1391: for (i = 0; i < pn; i++) {
1392: rcvncols = 0;
1393: for (j = 0; j < rootdegrees[i]; j++) {
1394: rcvncols += rootspace[rootspacesize];
1395: rootspaceoffsets[rootspacesize + 1] = rootspaceoffsets[rootspacesize] + rootspace[rootspacesize];
1396: rootspacesize++;
1397: }
1398: ptap->c_othi[i + 1] = ptap->c_othi[i] + rcvncols;
1399: }
1400: PetscCall(PetscFree(rootspace));
1402: PetscCall(PetscMalloc1(pon, &c_rmtoffsets));
1403: PetscCall(PetscSFScatterBegin(sf, MPIU_INT, rootspaceoffsets, c_rmtoffsets));
1404: PetscCall(PetscSFScatterEnd(sf, MPIU_INT, rootspaceoffsets, c_rmtoffsets));
1405: PetscCall(PetscSFDestroy(&sf));
1406: PetscCall(PetscFree(rootspaceoffsets));
1408: PetscCall(PetscCalloc1(ptap->c_rmti[pon], &iremote));
1409: nleaves = 0;
1410: for (i = 0; i < pon; i++) {
1411: owner = 0;
1412: ii = i / dof;
1413: PetscCall(PetscLayoutFindOwnerIndex(P->cmap, p->garray[ii], &owner, NULL));
1414: sendncols = ptap->c_rmti[i + 1] - ptap->c_rmti[i];
1415: for (j = 0; j < sendncols; j++) {
1416: iremote[nleaves].rank = owner;
1417: iremote[nleaves++].index = c_rmtoffsets[i] + j;
1418: }
1419: }
1420: PetscCall(PetscFree(c_rmtoffsets));
1421: PetscCall(PetscCalloc1(ptap->c_othi[pn], &c_othj));
1423: PetscCall(PetscSFCreate(comm, &ptap->sf));
1424: PetscCall(PetscSFSetGraph(ptap->sf, ptap->c_othi[pn], nleaves, NULL, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER));
1425: PetscCall(PetscSFSetFromOptions(ptap->sf));
1426: /* One to one map */
1427: PetscCall(PetscSFReduceBegin(ptap->sf, MPIU_INT, c_rmtj, c_othj, MPI_REPLACE));
1428: /* Get remote data */
1429: PetscCall(PetscSFReduceEnd(ptap->sf, MPIU_INT, c_rmtj, c_othj, MPI_REPLACE));
1430: PetscCall(PetscFree(c_rmtj));
1431: PetscCall(PetscMalloc2(pn, &dnz, pn, &onz));
1432: PetscCall(MatGetOwnershipRangeColumn(P, &pcstart, &pcend));
1433: pcstart *= dof;
1434: pcend *= dof;
1435: for (i = 0; i < pn; i++) {
1436: nzi = ptap->c_othi[i + 1] - ptap->c_othi[i];
1437: rdj = PetscSafePointerPlusOffset(c_othj, ptap->c_othi[i]);
1438: for (j = 0; j < nzi; j++) {
1439: col = rdj[j];
1440: /* diagonal part */
1441: if (col >= pcstart && col < pcend) {
1442: PetscCall(PetscHSetIAdd(htd[i], col));
1443: } else { /* off-diagonal */
1444: PetscCall(PetscHSetIAdd(hto[i], col));
1445: }
1446: }
1447: PetscCall(PetscHSetIGetSize(htd[i], &htsize));
1448: dnz[i] = htsize;
1449: PetscCall(PetscHSetIDestroy(&htd[i]));
1450: PetscCall(PetscHSetIGetSize(hto[i], &htsize));
1451: onz[i] = htsize;
1452: PetscCall(PetscHSetIDestroy(&hto[i]));
1453: }
1455: PetscCall(PetscFree2(htd, hto));
1456: PetscCall(PetscFree(c_othj));
1458: /* local sizes and preallocation */
1459: PetscCall(MatSetSizes(Cmpi, pn, pn, PETSC_DETERMINE, PETSC_DETERMINE));
1460: PetscCall(MatSetBlockSizes(Cmpi, dof > 1 ? dof : P->cmap->bs, dof > 1 ? dof : P->cmap->bs));
1461: PetscCall(MatMPIAIJSetPreallocation(Cmpi, 0, dnz, 0, onz));
1462: PetscCall(PetscFree2(dnz, onz));
1464: /* attach the supporting struct to Cmpi for reuse */
1465: Cmpi->product->data = ptap;
1466: Cmpi->product->destroy = MatProductCtxDestroy_MPIAIJ_PtAP;
1467: Cmpi->product->view = MatView_MPIAIJ_PtAP;
1469: /* Cmpi is not ready for use - assembly will be done by MatPtAPNumeric() */
1470: Cmpi->assembled = PETSC_FALSE;
1471: /* pick an algorithm */
1472: PetscOptionsBegin(PetscObjectComm((PetscObject)A), ((PetscObject)A)->prefix, "MatPtAP", "Mat");
1473: alg = 0;
1474: PetscCall(PetscOptionsEList("-matptap_allatonce_via", "PtAP allatonce numeric approach", "MatPtAP", algTypes, nalg, algTypes[alg], &alg, &flg));
1475: PetscOptionsEnd();
1476: switch (alg) {
1477: case 0:
1478: Cmpi->ops->ptapnumeric = MatPtAPNumeric_MPIAIJ_MPIAIJ_allatonce_merged;
1479: break;
1480: case 1:
1481: Cmpi->ops->ptapnumeric = MatPtAPNumeric_MPIAIJ_MPIAIJ_allatonce;
1482: break;
1483: default:
1484: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, " Unsupported allatonce numerical algorithm ");
1485: }
1486: PetscFunctionReturn(PETSC_SUCCESS);
1487: }
1489: PetscErrorCode MatPtAPSymbolic_MPIAIJ_MPIAIJ_allatonce_merged(Mat A, Mat P, PetscReal fill, Mat C)
1490: {
1491: PetscFunctionBegin;
1492: PetscCall(MatPtAPSymbolic_MPIAIJ_MPIXAIJ_allatonce_merged(A, P, 1, fill, C));
1493: PetscFunctionReturn(PETSC_SUCCESS);
1494: }
1496: PetscErrorCode MatPtAPSymbolic_MPIAIJ_MPIAIJ(Mat A, Mat P, PetscReal fill, Mat Cmpi)
1497: {
1498: MatProductCtx_APMPI *ptap;
1499: Mat_MPIAIJ *a = (Mat_MPIAIJ *)A->data, *p = (Mat_MPIAIJ *)P->data;
1500: MPI_Comm comm;
1501: PetscMPIInt size, rank, nsend, proc;
1502: PetscFreeSpaceList free_space = NULL, current_space = NULL;
1503: PetscInt am = A->rmap->n, pm = P->rmap->n, pN = P->cmap->N, pn = P->cmap->n;
1504: PetscInt *lnk, i, k, pnz, row;
1505: PetscBT lnkbt;
1506: PetscMPIInt tagi, tagj, *len_si, *len_s, *len_ri, nrecv;
1507: PETSC_UNUSED PetscMPIInt icompleted = 0;
1508: PetscInt **buf_rj, **buf_ri, **buf_ri_k;
1509: PetscInt len, *dnz, *onz, *owners, nzi, nspacedouble;
1510: PetscInt nrows, *buf_s, *buf_si, *buf_si_i, **nextrow, **nextci;
1511: MPI_Request *swaits, *rwaits;
1512: MPI_Status *sstatus, rstatus;
1513: PetscLayout rowmap;
1514: PetscInt *owners_co, *coi, *coj; /* i and j array of (p->B)^T*A*P - used in the communication */
1515: PetscMPIInt *len_r, *id_r; /* array of length of comm->size, store send/recv matrix values */
1516: PetscInt *api, *apj, *Jptr, apnz, *prmap = p->garray, con, j, ap_rmax = 0, Crmax, *aj, *ai, *pi;
1517: Mat_SeqAIJ *p_loc, *p_oth = NULL, *ad = (Mat_SeqAIJ *)a->A->data, *ao = NULL, *c_loc, *c_oth;
1518: PetscScalar *apv;
1519: PetscHMapI ta;
1520: MatType mtype;
1521: const char *prefix;
1522: PetscReal apfill;
1524: PetscFunctionBegin;
1525: MatCheckProduct(Cmpi, 4);
1526: PetscCheck(!Cmpi->product->data, PetscObjectComm((PetscObject)Cmpi), PETSC_ERR_PLIB, "Product data not empty");
1527: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
1528: PetscCallMPI(MPI_Comm_size(comm, &size));
1529: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1531: if (size > 1) ao = (Mat_SeqAIJ *)a->B->data;
1533: /* create symbolic parallel matrix Cmpi */
1534: PetscCall(MatGetType(A, &mtype));
1535: PetscCall(MatSetType(Cmpi, mtype));
1537: /* Do dense axpy in MatPtAPNumeric_MPIAIJ_MPIAIJ() */
1538: Cmpi->ops->ptapnumeric = MatPtAPNumeric_MPIAIJ_MPIAIJ;
1540: /* create struct MatProductCtx_APMPI and attached it to C later */
1541: PetscCall(PetscNew(&ptap));
1542: ptap->reuse = MAT_INITIAL_MATRIX;
1543: ptap->algType = 1;
1545: /* get P_oth by taking rows of P (= non-zero cols of local A) from other processors */
1546: PetscCall(MatGetBrowsOfAoCols_MPIAIJ(A, P, MAT_INITIAL_MATRIX, &ptap->startsj_s, &ptap->startsj_r, &ptap->bufa, &ptap->P_oth));
1547: /* get P_loc by taking all local rows of P */
1548: PetscCall(MatMPIAIJGetLocalMat(P, MAT_INITIAL_MATRIX, &ptap->P_loc));
1550: /* (0) compute Rd = Pd^T, Ro = Po^T */
1551: PetscCall(MatTranspose(p->A, MAT_INITIAL_MATRIX, &ptap->Rd));
1552: PetscCall(MatTranspose(p->B, MAT_INITIAL_MATRIX, &ptap->Ro));
1554: /* (1) compute symbolic AP = A_loc*P = Ad*P_loc + Ao*P_oth (api,apj) */
1555: p_loc = (Mat_SeqAIJ *)ptap->P_loc->data;
1556: if (ptap->P_oth) p_oth = (Mat_SeqAIJ *)ptap->P_oth->data;
1558: /* create and initialize a linked list */
1559: PetscCall(PetscHMapICreateWithSize(pn, &ta)); /* for compute AP_loc and Cmpi */
1560: MatRowMergeMax_SeqAIJ(p_loc, ptap->P_loc->rmap->N, ta);
1561: MatRowMergeMax_SeqAIJ(p_oth, ptap->P_oth->rmap->N, ta);
1562: PetscCall(PetscHMapIGetSize(ta, &Crmax)); /* Crmax = nnz(sum of Prows) */
1563: /* printf("[%d] est %d, Crmax %d; pN %d\n",rank,5*(p_loc->rmax+p_oth->rmax + (PetscInt)(1.e-2*pN)),Crmax,pN); */
1565: PetscCall(PetscLLCondensedCreate(Crmax, pN, &lnk, &lnkbt));
1567: /* Initial FreeSpace size is fill*(nnz(A) + nnz(P)) */
1568: if (ao) {
1569: PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(fill, PetscIntSumTruncate(ad->i[am], PetscIntSumTruncate(ao->i[am], p_loc->i[pm]))), &free_space));
1570: } else {
1571: PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(fill, PetscIntSumTruncate(ad->i[am], p_loc->i[pm])), &free_space));
1572: }
1573: current_space = free_space;
1574: nspacedouble = 0;
1576: PetscCall(PetscMalloc1(am + 1, &api));
1577: api[0] = 0;
1578: for (i = 0; i < am; i++) {
1579: /* diagonal portion: Ad[i,:]*P */
1580: ai = ad->i;
1581: pi = p_loc->i;
1582: nzi = ai[i + 1] - ai[i];
1583: aj = PetscSafePointerPlusOffset(ad->j, ai[i]);
1584: for (j = 0; j < nzi; j++) {
1585: row = aj[j];
1586: pnz = pi[row + 1] - pi[row];
1587: Jptr = p_loc->j + pi[row];
1588: /* add non-zero cols of P into the sorted linked list lnk */
1589: PetscCall(PetscLLCondensedAddSorted(pnz, Jptr, lnk, lnkbt));
1590: }
1591: /* off-diagonal portion: Ao[i,:]*P */
1592: if (ao) {
1593: ai = ao->i;
1594: pi = p_oth->i;
1595: nzi = ai[i + 1] - ai[i];
1596: aj = PetscSafePointerPlusOffset(ao->j, ai[i]);
1597: for (j = 0; j < nzi; j++) {
1598: row = aj[j];
1599: pnz = pi[row + 1] - pi[row];
1600: Jptr = PetscSafePointerPlusOffset(p_oth->j, pi[row]);
1601: PetscCall(PetscLLCondensedAddSorted(pnz, Jptr, lnk, lnkbt));
1602: }
1603: }
1604: apnz = lnk[0];
1605: api[i + 1] = api[i] + apnz;
1606: if (ap_rmax < apnz) ap_rmax = apnz;
1608: /* if free space is not available, double the total space in the list */
1609: if (current_space->local_remaining < apnz) {
1610: PetscCall(PetscFreeSpaceGet(PetscIntSumTruncate(apnz, current_space->total_array_size), ¤t_space));
1611: nspacedouble++;
1612: }
1614: /* Copy data into free space, then initialize lnk */
1615: PetscCall(PetscLLCondensedClean(pN, apnz, current_space->array, lnk, lnkbt));
1617: current_space->array = PetscSafePointerPlusOffset(current_space->array, apnz);
1618: current_space->local_used += apnz;
1619: current_space->local_remaining -= apnz;
1620: }
1621: /* Allocate space for apj and apv, initialize apj, and */
1622: /* destroy list of free space and other temporary array(s) */
1623: PetscCall(PetscMalloc2(api[am], &apj, api[am], &apv));
1624: PetscCall(PetscFreeSpaceContiguous(&free_space, apj));
1625: PetscCall(PetscLLDestroy(lnk, lnkbt));
1627: /* Create AP_loc for reuse */
1628: PetscCall(MatCreateSeqAIJWithArrays(PETSC_COMM_SELF, am, pN, api, apj, apv, &ptap->AP_loc));
1629: PetscCall(MatSetType(ptap->AP_loc, ((PetscObject)p->A)->type_name));
1630: if (PetscDefined(USE_INFO)) {
1631: if (ao) apfill = (PetscReal)api[am] / (ad->i[am] + ao->i[am] + p_loc->i[pm] + 1);
1632: else apfill = (PetscReal)api[am] / (ad->i[am] + p_loc->i[pm] + 1);
1633: ptap->AP_loc->info.mallocs = nspacedouble;
1634: ptap->AP_loc->info.fill_ratio_given = fill;
1635: ptap->AP_loc->info.fill_ratio_needed = apfill;
1637: if (api[am]) {
1638: PetscCall(PetscInfo(ptap->AP_loc, "Nonscalable algorithm, AP_loc reallocs %" PetscInt_FMT "; Fill ratio: given %g needed %g.\n", nspacedouble, (double)fill, (double)apfill));
1639: PetscCall(PetscInfo(ptap->AP_loc, "Use MatPtAP(A,B,MatReuse,%g,&C) for best AP_loc performance.;\n", (double)apfill));
1640: } else PetscCall(PetscInfo(ptap->AP_loc, "Nonscalable algorithm, AP_loc is empty \n"));
1641: }
1643: /* (2-1) compute symbolic Co = Ro*AP_loc */
1644: PetscCall(MatGetOptionsPrefix(A, &prefix));
1645: PetscCall(MatSetOptionsPrefix(ptap->Ro, prefix));
1646: PetscCall(MatAppendOptionsPrefix(ptap->Ro, "inner_offdiag_"));
1647: PetscCall(MatProductCreate(ptap->Ro, ptap->AP_loc, NULL, &ptap->C_oth));
1648: PetscCall(MatGetOptionsPrefix(Cmpi, &prefix));
1649: PetscCall(MatSetOptionsPrefix(ptap->C_oth, prefix));
1650: PetscCall(MatAppendOptionsPrefix(ptap->C_oth, "inner_C_oth_"));
1651: PetscCall(MatProductSetType(ptap->C_oth, MATPRODUCT_AB));
1652: PetscCall(MatProductSetAlgorithm(ptap->C_oth, "default"));
1653: PetscCall(MatProductSetFill(ptap->C_oth, fill));
1654: PetscCall(MatProductSetFromOptions(ptap->C_oth));
1655: PetscCall(MatProductSymbolic(ptap->C_oth));
1657: /* (3) send coj of C_oth to other processors */
1658: /* determine row ownership */
1659: PetscCall(PetscLayoutCreate(comm, &rowmap));
1660: rowmap->n = pn;
1661: rowmap->bs = 1;
1662: PetscCall(PetscLayoutSetUp(rowmap));
1663: owners = rowmap->range;
1665: /* determine the number of messages to send, their lengths */
1666: PetscCall(PetscMalloc4(size, &len_s, size, &len_si, size, &sstatus, size + 2, &owners_co));
1667: PetscCall(PetscArrayzero(len_s, size));
1668: PetscCall(PetscArrayzero(len_si, size));
1670: c_oth = (Mat_SeqAIJ *)ptap->C_oth->data;
1671: coi = c_oth->i;
1672: coj = c_oth->j;
1673: con = ptap->C_oth->rmap->n;
1674: proc = 0;
1675: for (i = 0; i < con; i++) {
1676: while (prmap[i] >= owners[proc + 1]) proc++;
1677: len_si[proc]++; /* num of rows in Co(=Pt*AP) to be sent to [proc] */
1678: len_s[proc] += coi[i + 1] - coi[i]; /* num of nonzeros in Co to be sent to [proc] */
1679: }
1681: len = 0; /* max length of buf_si[], see (4) */
1682: owners_co[0] = 0;
1683: nsend = 0;
1684: for (proc = 0; proc < size; proc++) {
1685: owners_co[proc + 1] = owners_co[proc] + len_si[proc];
1686: if (len_s[proc]) {
1687: nsend++;
1688: len_si[proc] = 2 * (len_si[proc] + 1); /* length of buf_si to be sent to [proc] */
1689: len += len_si[proc];
1690: }
1691: }
1693: /* determine the number and length of messages to receive for coi and coj */
1694: PetscCall(PetscGatherNumberOfMessages(comm, NULL, len_s, &nrecv));
1695: PetscCall(PetscGatherMessageLengths2(comm, nsend, nrecv, len_s, len_si, &id_r, &len_r, &len_ri));
1697: /* post the Irecv and Isend of coj */
1698: PetscCall(PetscCommGetNewTag(comm, &tagj));
1699: PetscCall(PetscPostIrecvInt(comm, tagj, nrecv, id_r, len_r, &buf_rj, &rwaits));
1700: PetscCall(PetscMalloc1(nsend + 1, &swaits));
1701: for (proc = 0, k = 0; proc < size; proc++) {
1702: if (!len_s[proc]) continue;
1703: i = owners_co[proc];
1704: PetscCallMPI(MPIU_Isend(coj + coi[i], len_s[proc], MPIU_INT, proc, tagj, comm, swaits + k));
1705: k++;
1706: }
1708: /* (2-2) compute symbolic C_loc = Rd*AP_loc */
1709: PetscCall(MatSetOptionsPrefix(ptap->Rd, prefix));
1710: PetscCall(MatAppendOptionsPrefix(ptap->Rd, "inner_diag_"));
1711: PetscCall(MatProductCreate(ptap->Rd, ptap->AP_loc, NULL, &ptap->C_loc));
1712: PetscCall(MatGetOptionsPrefix(Cmpi, &prefix));
1713: PetscCall(MatSetOptionsPrefix(ptap->C_loc, prefix));
1714: PetscCall(MatAppendOptionsPrefix(ptap->C_loc, "inner_C_loc_"));
1715: PetscCall(MatProductSetType(ptap->C_loc, MATPRODUCT_AB));
1716: PetscCall(MatProductSetAlgorithm(ptap->C_loc, "default"));
1717: PetscCall(MatProductSetFill(ptap->C_loc, fill));
1718: PetscCall(MatProductSetFromOptions(ptap->C_loc));
1719: PetscCall(MatProductSymbolic(ptap->C_loc));
1721: c_loc = (Mat_SeqAIJ *)ptap->C_loc->data;
1723: /* receives coj are complete */
1724: for (i = 0; i < nrecv; i++) PetscCallMPI(MPI_Waitany(nrecv, rwaits, &icompleted, &rstatus));
1725: PetscCall(PetscFree(rwaits));
1726: if (nsend) PetscCallMPI(MPI_Waitall(nsend, swaits, sstatus));
1728: /* add received column indices into ta to update Crmax */
1729: for (k = 0; k < nrecv; k++) { /* k-th received message */
1730: Jptr = buf_rj[k];
1731: for (j = 0; j < len_r[k]; j++) PetscCall(PetscHMapISet(ta, *(Jptr + j) + 1, 1));
1732: }
1733: PetscCall(PetscHMapIGetSize(ta, &Crmax));
1734: PetscCall(PetscHMapIDestroy(&ta));
1736: /* (4) send and recv coi */
1737: PetscCall(PetscCommGetNewTag(comm, &tagi));
1738: PetscCall(PetscPostIrecvInt(comm, tagi, nrecv, id_r, len_ri, &buf_ri, &rwaits));
1739: PetscCall(PetscMalloc1(len + 1, &buf_s));
1740: buf_si = buf_s; /* points to the beginning of k-th msg to be sent */
1741: for (proc = 0, k = 0; proc < size; proc++) {
1742: if (!len_s[proc]) continue;
1743: /* form outgoing message for i-structure:
1744: buf_si[0]: nrows to be sent
1745: [1:nrows]: row index (global)
1746: [nrows+1:2*nrows+1]: i-structure index
1747: */
1748: nrows = len_si[proc] / 2 - 1; /* num of rows in Co to be sent to [proc] */
1749: buf_si_i = buf_si + nrows + 1;
1750: buf_si[0] = nrows;
1751: buf_si_i[0] = 0;
1752: nrows = 0;
1753: for (i = owners_co[proc]; i < owners_co[proc + 1]; i++) {
1754: nzi = coi[i + 1] - coi[i];
1755: buf_si_i[nrows + 1] = buf_si_i[nrows] + nzi; /* i-structure */
1756: buf_si[nrows + 1] = prmap[i] - owners[proc]; /* local row index */
1757: nrows++;
1758: }
1759: PetscCallMPI(MPIU_Isend(buf_si, len_si[proc], MPIU_INT, proc, tagi, comm, swaits + k));
1760: k++;
1761: buf_si += len_si[proc];
1762: }
1763: for (i = 0; i < nrecv; i++) PetscCallMPI(MPI_Waitany(nrecv, rwaits, &icompleted, &rstatus));
1764: PetscCall(PetscFree(rwaits));
1765: if (nsend) PetscCallMPI(MPI_Waitall(nsend, swaits, sstatus));
1767: PetscCall(PetscFree4(len_s, len_si, sstatus, owners_co));
1768: PetscCall(PetscFree(len_ri));
1769: PetscCall(PetscFree(swaits));
1770: PetscCall(PetscFree(buf_s));
1772: /* (5) compute the local portion of Cmpi */
1773: /* set initial free space to be Crmax, sufficient for holding nonzeros in each row of Cmpi */
1774: PetscCall(PetscFreeSpaceGet(Crmax, &free_space));
1775: current_space = free_space;
1777: PetscCall(PetscMalloc3(nrecv, &buf_ri_k, nrecv, &nextrow, nrecv, &nextci));
1778: for (k = 0; k < nrecv; k++) {
1779: buf_ri_k[k] = buf_ri[k]; /* beginning of k-th recved i-structure */
1780: nrows = *buf_ri_k[k];
1781: nextrow[k] = buf_ri_k[k] + 1; /* next row number of k-th recved i-structure */
1782: nextci[k] = buf_ri_k[k] + (nrows + 1); /* points to the next i-structure of k-th recved i-structure */
1783: }
1785: MatPreallocateBegin(comm, pn, pn, dnz, onz);
1786: PetscCall(PetscLLCondensedCreate(Crmax, pN, &lnk, &lnkbt));
1787: for (i = 0; i < pn; i++) {
1788: /* add C_loc into Cmpi */
1789: nzi = c_loc->i[i + 1] - c_loc->i[i];
1790: Jptr = PetscSafePointerPlusOffset(c_loc->j, c_loc->i[i]);
1791: PetscCall(PetscLLCondensedAddSorted(nzi, Jptr, lnk, lnkbt));
1793: /* add received col data into lnk */
1794: for (k = 0; k < nrecv; k++) { /* k-th received message */
1795: if (i == *nextrow[k]) { /* i-th row */
1796: nzi = *(nextci[k] + 1) - *nextci[k];
1797: Jptr = buf_rj[k] + *nextci[k];
1798: PetscCall(PetscLLCondensedAddSorted(nzi, Jptr, lnk, lnkbt));
1799: nextrow[k]++;
1800: nextci[k]++;
1801: }
1802: }
1803: nzi = lnk[0];
1805: /* copy data into free space, then initialize lnk */
1806: PetscCall(PetscLLCondensedClean(pN, nzi, current_space->array, lnk, lnkbt));
1807: PetscCall(MatPreallocateSet(i + owners[rank], nzi, current_space->array, dnz, onz));
1808: }
1809: PetscCall(PetscFree3(buf_ri_k, nextrow, nextci));
1810: PetscCall(PetscLLDestroy(lnk, lnkbt));
1811: PetscCall(PetscFreeSpaceDestroy(free_space));
1813: /* local sizes and preallocation */
1814: PetscCall(MatSetSizes(Cmpi, pn, pn, PETSC_DETERMINE, PETSC_DETERMINE));
1815: if (P->cmap->bs > 0) {
1816: PetscCall(PetscLayoutSetBlockSize(Cmpi->rmap, P->cmap->bs));
1817: PetscCall(PetscLayoutSetBlockSize(Cmpi->cmap, P->cmap->bs));
1818: }
1819: PetscCall(MatMPIAIJSetPreallocation(Cmpi, 0, dnz, 0, onz));
1820: MatPreallocateEnd(dnz, onz);
1822: /* members in merge */
1823: PetscCall(PetscFree(id_r));
1824: PetscCall(PetscFree(len_r));
1825: PetscCall(PetscFree(buf_ri[0]));
1826: PetscCall(PetscFree(buf_ri));
1827: PetscCall(PetscFree(buf_rj[0]));
1828: PetscCall(PetscFree(buf_rj));
1829: PetscCall(PetscLayoutDestroy(&rowmap));
1831: PetscCall(PetscCalloc1(pN, &ptap->apa));
1833: /* attach the supporting struct to Cmpi for reuse */
1834: Cmpi->product->data = ptap;
1835: Cmpi->product->destroy = MatProductCtxDestroy_MPIAIJ_PtAP;
1836: Cmpi->product->view = MatView_MPIAIJ_PtAP;
1838: /* Cmpi is not ready for use - assembly will be done by MatPtAPNumeric() */
1839: Cmpi->assembled = PETSC_FALSE;
1840: PetscFunctionReturn(PETSC_SUCCESS);
1841: }
1843: PetscErrorCode MatPtAPNumeric_MPIAIJ_MPIAIJ(Mat A, Mat P, Mat C)
1844: {
1845: Mat_MPIAIJ *a = (Mat_MPIAIJ *)A->data, *p = (Mat_MPIAIJ *)P->data;
1846: Mat_SeqAIJ *ad = (Mat_SeqAIJ *)a->A->data, *ao = (Mat_SeqAIJ *)a->B->data;
1847: Mat_SeqAIJ *ap, *p_loc, *p_oth = NULL, *c_seq;
1848: MatProductCtx_APMPI *ptap;
1849: Mat AP_loc, C_loc, C_oth;
1850: PetscInt i, rstart, rend, cm, ncols, row;
1851: PetscInt *api, *apj, am = A->rmap->n, j, col, apnz;
1852: PetscScalar *apa;
1853: const PetscInt *cols;
1854: const PetscScalar *vals, *array, *dummy1, *dummy2, *dummy3, *dummy4;
1856: PetscFunctionBegin;
1857: MatCheckProduct(C, 3);
1858: ptap = (MatProductCtx_APMPI *)C->product->data;
1859: PetscCheck(ptap, PetscObjectComm((PetscObject)C), PETSC_ERR_ARG_WRONGSTATE, "PtAP cannot be computed. Missing data");
1860: PetscCheck(ptap->AP_loc, PetscObjectComm((PetscObject)C), PETSC_ERR_ARG_WRONGSTATE, "PtAP cannot be reused. Do not call MatProductClear()");
1862: PetscCall(MatZeroEntries(C));
1863: /* 1) get R = Pd^T,Ro = Po^T */
1864: if (ptap->reuse == MAT_REUSE_MATRIX) {
1865: PetscCall(MatTranspose(p->A, MAT_REUSE_MATRIX, &ptap->Rd));
1866: PetscCall(MatTranspose(p->B, MAT_REUSE_MATRIX, &ptap->Ro));
1867: }
1869: /* 2) get AP_loc */
1870: AP_loc = ptap->AP_loc;
1871: ap = (Mat_SeqAIJ *)AP_loc->data;
1873: /* 2-1) get P_oth = ptap->P_oth and P_loc = ptap->P_loc */
1874: if (ptap->reuse == MAT_REUSE_MATRIX) {
1875: /* P_oth and P_loc are obtained in MatPtASymbolic() when reuse == MAT_INITIAL_MATRIX */
1876: PetscCall(MatGetBrowsOfAoCols_MPIAIJ(A, P, MAT_REUSE_MATRIX, &ptap->startsj_s, &ptap->startsj_r, &ptap->bufa, &ptap->P_oth));
1877: PetscCall(MatMPIAIJGetLocalMat(P, MAT_REUSE_MATRIX, &ptap->P_loc));
1878: }
1880: /* 2-2) compute numeric A_loc*P - dominating part */
1881: /* get data from symbolic products */
1882: p_loc = (Mat_SeqAIJ *)ptap->P_loc->data;
1883: if (ptap->P_oth) p_oth = (Mat_SeqAIJ *)ptap->P_oth->data;
1884: api = ap->i;
1885: apj = ap->j;
1887: // Use MatSeqAIJGetArrayXXX to sync the matrix on host before accessing its values in the style of (Mat_SeqAIJ *)dat->a on host
1888: PetscCall(MatSeqAIJGetArrayWrite(AP_loc, &apa));
1889: PetscCall(MatSeqAIJGetArrayRead(a->A, &dummy1));
1890: PetscCall(MatSeqAIJGetArrayRead(a->B, &dummy2));
1891: PetscCall(MatSeqAIJGetArrayRead(ptap->P_loc, &dummy3));
1892: if (ptap->P_oth) PetscCall(MatSeqAIJGetArrayRead(ptap->P_oth, &dummy4));
1893: for (i = 0; i < am; i++) {
1894: /* AP[i,:] = A[i,:]*P = Ad*P_loc + Ao*P_oth */
1895: AProw_nonscalable(i, ad, ao, p_loc, p_oth, ptap->apa); // Directly access the value arrays from the Mat_SeqAIJ structs
1896: apnz = api[i + 1] - api[i];
1897: for (j = 0; j < apnz; j++) {
1898: col = apj[j + api[i]];
1899: apa[j + ap->i[i]] = ptap->apa[col];
1900: ptap->apa[col] = 0.0;
1901: }
1902: }
1903: PetscCall(MatSeqAIJRestoreArrayWrite(AP_loc, &apa));
1904: PetscCall(MatSeqAIJRestoreArrayRead(a->A, &dummy1));
1905: PetscCall(MatSeqAIJRestoreArrayRead(a->B, &dummy2));
1906: PetscCall(MatSeqAIJRestoreArrayRead(ptap->P_loc, &dummy3));
1907: if (ptap->P_oth) PetscCall(MatSeqAIJRestoreArrayRead(ptap->P_oth, &dummy4));
1909: /* We have modified the contents of local matrix AP_loc and must increase its ObjectState, since we are not doing AssemblyBegin/End on it. */
1910: PetscCall(PetscObjectStateIncrease((PetscObject)AP_loc));
1912: /* 3) C_loc = Rd*AP_loc, C_oth = Ro*AP_loc */
1913: PetscCall(MatProductNumeric(ptap->C_loc));
1914: PetscCall(MatProductNumeric(ptap->C_oth));
1915: C_loc = ptap->C_loc;
1916: C_oth = ptap->C_oth;
1918: /* add C_loc and Co to C */
1919: PetscCall(MatGetOwnershipRange(C, &rstart, &rend));
1921: /* C_loc -> C */
1922: cm = C_loc->rmap->N;
1923: c_seq = (Mat_SeqAIJ *)C_loc->data;
1924: cols = c_seq->j;
1926: PetscCall(MatSeqAIJGetArrayRead(C_loc, &array));
1927: vals = array;
1929: /* The (fast) MatSetValues_MPIAIJ_CopyFromCSRFormat function can only be used when C->was_assembled is PETSC_FALSE and */
1930: /* when there are no off-processor parts. */
1931: /* If was_assembled is true, then the statement aj[rowstart_diag+dnz_row] = mat_j[col] - cstart; in MatSetValues_MPIAIJ_CopyFromCSRFormat */
1932: /* is no longer true. Then the more complex function MatSetValues_MPIAIJ() has to be used, where the column index is looked up from */
1933: /* a table, and other, more complex stuff has to be done. */
1934: if (C->assembled) {
1935: C->was_assembled = PETSC_TRUE;
1936: C->assembled = PETSC_FALSE;
1937: }
1938: if (C->was_assembled) {
1939: for (i = 0; i < cm; i++) {
1940: ncols = c_seq->i[i + 1] - c_seq->i[i];
1941: row = rstart + i;
1942: PetscCall(MatSetValues_MPIAIJ(C, 1, &row, ncols, cols, vals, ADD_VALUES));
1943: cols = PetscSafePointerPlusOffset(cols, ncols);
1944: vals += ncols;
1945: }
1946: } else {
1947: PetscCall(MatSetValues_MPIAIJ_CopyFromCSRFormat(C, c_seq->j, c_seq->i, vals));
1948: }
1949: PetscCall(MatSeqAIJRestoreArrayRead(C_loc, &array));
1951: /* Co -> C, off-processor part */
1952: cm = C_oth->rmap->N;
1953: c_seq = (Mat_SeqAIJ *)C_oth->data;
1954: cols = c_seq->j;
1955: PetscCall(MatSeqAIJGetArrayRead(C_oth, &array));
1956: vals = array;
1957: for (i = 0; i < cm; i++) {
1958: ncols = c_seq->i[i + 1] - c_seq->i[i];
1959: row = p->garray[i];
1960: PetscCall(MatSetValues(C, 1, &row, ncols, cols, vals, ADD_VALUES));
1961: cols += ncols;
1962: vals += ncols;
1963: }
1964: PetscCall(MatSeqAIJRestoreArrayRead(C_oth, &array));
1966: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
1967: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
1969: ptap->reuse = MAT_REUSE_MATRIX;
1970: PetscFunctionReturn(PETSC_SUCCESS);
1971: }
1973: PETSC_INTERN PetscErrorCode MatProductSymbolic_PtAP_MPIAIJ_MPIAIJ(Mat C)
1974: {
1975: Mat_Product *product = C->product;
1976: Mat A = product->A, P = product->B;
1977: MatProductAlgorithm alg = product->alg;
1978: PetscReal fill = product->fill;
1979: PetscBool flg;
1981: PetscFunctionBegin;
1982: /* scalable: do R=P^T locally, then C=R*A*P */
1983: PetscCall(PetscStrcmp(alg, "scalable", &flg));
1984: if (flg) {
1985: PetscCall(MatPtAPSymbolic_MPIAIJ_MPIAIJ_scalable(A, P, product->fill, C));
1986: C->ops->productnumeric = MatProductNumeric_PtAP;
1987: goto next;
1988: }
1990: /* nonscalable: do R=P^T locally, then C=R*A*P */
1991: PetscCall(PetscStrcmp(alg, "nonscalable", &flg));
1992: if (flg) {
1993: PetscCall(MatPtAPSymbolic_MPIAIJ_MPIAIJ(A, P, fill, C));
1994: goto next;
1995: }
1997: /* allatonce */
1998: PetscCall(PetscStrcmp(alg, "allatonce", &flg));
1999: if (flg) {
2000: PetscCall(MatPtAPSymbolic_MPIAIJ_MPIAIJ_allatonce(A, P, fill, C));
2001: goto next;
2002: }
2004: /* allatonce_merged */
2005: PetscCall(PetscStrcmp(alg, "allatonce_merged", &flg));
2006: if (flg) {
2007: PetscCall(MatPtAPSymbolic_MPIAIJ_MPIAIJ_allatonce_merged(A, P, fill, C));
2008: goto next;
2009: }
2011: /* backend general code */
2012: PetscCall(PetscStrcmp(alg, "backend", &flg));
2013: if (flg) {
2014: PetscCall(MatProductSymbolic_MPIAIJBACKEND(C));
2015: PetscFunctionReturn(PETSC_SUCCESS);
2016: }
2018: /* hypre */
2019: #if PetscDefined(HAVE_HYPRE)
2020: PetscCall(PetscStrcmp(alg, "hypre", &flg));
2021: if (flg) {
2022: PetscCall(MatPtAPSymbolic_AIJ_AIJ_wHYPRE(A, P, fill, C));
2023: C->ops->productnumeric = MatProductNumeric_PtAP;
2024: PetscFunctionReturn(PETSC_SUCCESS);
2025: }
2026: #endif
2027: SETERRQ(PetscObjectComm((PetscObject)C), PETSC_ERR_SUP, "Mat Product Algorithm is not supported");
2029: next:
2030: C->ops->productnumeric = MatProductNumeric_PtAP;
2031: PetscFunctionReturn(PETSC_SUCCESS);
2032: }