Actual source code: mpisell.c
1: #include <../src/mat/impls/aij/mpi/mpiaij.h>
2: #include <../src/mat/impls/sell/mpi/mpisell.h>
3: #include <petsc/private/vecimpl.h>
4: #include <petsc/private/isimpl.h>
5: #include <petscblaslapack.h>
6: #include <petscsf.h>
8: /*MC
9: MATSELL - MATSELL = "sell" - A matrix type to be used for sparse matrices.
11: This matrix type is identical to `MATSEQSELL` when constructed with a single process communicator,
12: and `MATMPISELL` otherwise. As a result, for single process communicators,
13: `MatSeqSELLSetPreallocation()` is supported, and similarly `MatMPISELLSetPreallocation()` is supported
14: for communicators controlling multiple processes. It is recommended that you call both of
15: the above preallocation routines for simplicity.
17: Options Database Keys:
18: . -mat_type sell - sets the matrix type to `MATSELL` during a call to `MatSetFromOptions()`
20: Level: beginner
22: .seealso: `Mat`, `MATAIJ`, `MATBAIJ`, `MATSBAIJ`, `MatCreateSELL()`, `MatCreateSeqSELL()`, `MATSEQSELL`, `MATMPISELL`
23: M*/
25: static PetscErrorCode MatDiagonalSet_MPISELL(Mat Y, Vec D, InsertMode is)
26: {
27: Mat_MPISELL *sell = (Mat_MPISELL *)Y->data;
29: PetscFunctionBegin;
30: if (Y->assembled && Y->rmap->rstart == Y->cmap->rstart && Y->rmap->rend == Y->cmap->rend) {
31: PetscCall(MatDiagonalSet(sell->A, D, is));
32: } else {
33: PetscCall(MatDiagonalSet_Default(Y, D, is));
34: }
35: PetscFunctionReturn(PETSC_SUCCESS);
36: }
38: /*
39: Local utility routine that creates a mapping from the global column
40: number to the local number in the off-diagonal part of the local
41: storage of the matrix. When PETSC_USE_CTABLE is used this is scalable at
42: a slightly higher hash table cost; without it it is not scalable (each processor
43: has an order N integer array but is fast to access.
44: */
45: PetscErrorCode MatCreateColmap_MPISELL_Private(Mat mat)
46: {
47: Mat_MPISELL *sell = (Mat_MPISELL *)mat->data;
48: PetscInt n = sell->B->cmap->n, i;
50: PetscFunctionBegin;
51: PetscCheck(sell->garray, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MPISELL Matrix was assembled but is missing garray");
52: #if PetscDefined(USE_CTABLE)
53: PetscCall(PetscHMapICreateWithSize(n, &sell->colmap));
54: for (i = 0; i < n; i++) PetscCall(PetscHMapISet(sell->colmap, sell->garray[i] + 1, i + 1));
55: #else
56: PetscCall(PetscCalloc1(mat->cmap->N + 1, &sell->colmap));
57: for (i = 0; i < n; i++) sell->colmap[sell->garray[i]] = i + 1;
58: #endif
59: PetscFunctionReturn(PETSC_SUCCESS);
60: }
62: #define MatSetValues_SeqSELL_A_Private(row, col, value, addv, orow, ocol) \
63: { \
64: if (col <= lastcol1) low1 = 0; \
65: else high1 = nrow1; \
66: lastcol1 = col; \
67: while (high1 - low1 > 5) { \
68: t = (low1 + high1) / 2; \
69: if (cp1[sliceheight * t] > col) high1 = t; \
70: else low1 = t; \
71: } \
72: for (_i = low1; _i < high1; _i++) { \
73: if (cp1[sliceheight * _i] > col) break; \
74: if (cp1[sliceheight * _i] == col) { \
75: if (addv == ADD_VALUES) vp1[sliceheight * _i] += value; \
76: else vp1[sliceheight * _i] = value; \
77: inserted = PETSC_TRUE; \
78: goto a_noinsert; \
79: } \
80: } \
81: if (value == 0.0 && ignorezeroentries) { \
82: low1 = 0; \
83: high1 = nrow1; \
84: goto a_noinsert; \
85: } \
86: if (nonew == 1) { \
87: low1 = 0; \
88: high1 = nrow1; \
89: goto a_noinsert; \
90: } \
91: PetscCheck(nonew != -1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Inserting a new nonzero at global row/column (%" PetscInt_FMT ", %" PetscInt_FMT ") into matrix", orow, ocol); \
92: MatSeqXSELLReallocateSELL(A, am, 1, nrow1, a->sliidx, a->sliceheight, row / sliceheight, row, col, a->colidx, a->val, cp1, vp1, nonew, MatScalar); \
93: /* shift up all the later entries in this row */ \
94: for (ii = nrow1 - 1; ii >= _i; ii--) { \
95: cp1[sliceheight * (ii + 1)] = cp1[sliceheight * ii]; \
96: vp1[sliceheight * (ii + 1)] = vp1[sliceheight * ii]; \
97: } \
98: cp1[sliceheight * _i] = col; \
99: vp1[sliceheight * _i] = value; \
100: a->nz++; \
101: nrow1++; \
102: a_noinsert:; \
103: a->rlen[row] = nrow1; \
104: }
106: #define MatSetValues_SeqSELL_B_Private(row, col, value, addv, orow, ocol) \
107: { \
108: if (col <= lastcol2) low2 = 0; \
109: else high2 = nrow2; \
110: lastcol2 = col; \
111: while (high2 - low2 > 5) { \
112: t = (low2 + high2) / 2; \
113: if (cp2[sliceheight * t] > col) high2 = t; \
114: else low2 = t; \
115: } \
116: for (_i = low2; _i < high2; _i++) { \
117: if (cp2[sliceheight * _i] > col) break; \
118: if (cp2[sliceheight * _i] == col) { \
119: if (addv == ADD_VALUES) vp2[sliceheight * _i] += value; \
120: else vp2[sliceheight * _i] = value; \
121: inserted = PETSC_TRUE; \
122: goto b_noinsert; \
123: } \
124: } \
125: if (value == 0.0 && ignorezeroentries) { \
126: low2 = 0; \
127: high2 = nrow2; \
128: goto b_noinsert; \
129: } \
130: if (nonew == 1) { \
131: low2 = 0; \
132: high2 = nrow2; \
133: goto b_noinsert; \
134: } \
135: PetscCheck(nonew != -1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Inserting a new nonzero at global row/column (%" PetscInt_FMT ", %" PetscInt_FMT ") into matrix", orow, ocol); \
136: MatSeqXSELLReallocateSELL(B, bm, 1, nrow2, b->sliidx, b->sliceheight, row / sliceheight, row, col, b->colidx, b->val, cp2, vp2, nonew, MatScalar); \
137: /* shift up all the later entries in this row */ \
138: for (ii = nrow2 - 1; ii >= _i; ii--) { \
139: cp2[sliceheight * (ii + 1)] = cp2[sliceheight * ii]; \
140: vp2[sliceheight * (ii + 1)] = vp2[sliceheight * ii]; \
141: } \
142: cp2[sliceheight * _i] = col; \
143: vp2[sliceheight * _i] = value; \
144: b->nz++; \
145: nrow2++; \
146: b_noinsert:; \
147: b->rlen[row] = nrow2; \
148: }
150: static PetscErrorCode MatSetValues_MPISELL(Mat mat, PetscInt m, const PetscInt im[], PetscInt n, const PetscInt in[], const PetscScalar v[], InsertMode addv)
151: {
152: Mat_MPISELL *sell = (Mat_MPISELL *)mat->data;
153: PetscScalar value;
154: PetscInt i, j, rstart = mat->rmap->rstart, rend = mat->rmap->rend, shift1, shift2;
155: PetscInt cstart = mat->cmap->rstart, cend = mat->cmap->rend, row, col;
156: PetscBool roworiented = sell->roworiented;
158: /* Some Variables required in the macro */
159: Mat A = sell->A;
160: Mat_SeqSELL *a = (Mat_SeqSELL *)A->data;
161: PetscBool ignorezeroentries = a->ignorezeroentries, found;
162: Mat B = sell->B;
163: Mat_SeqSELL *b = (Mat_SeqSELL *)B->data;
164: PetscInt *cp1, *cp2, ii, _i, nrow1, nrow2, low1, high1, low2, high2, t, lastcol1, lastcol2, sliceheight = a->sliceheight;
165: MatScalar *vp1, *vp2;
167: PetscFunctionBegin;
168: for (i = 0; i < m; i++) {
169: if (im[i] < 0) continue;
170: PetscCheck(im[i] < mat->rmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, im[i], mat->rmap->N - 1);
171: if (im[i] >= rstart && im[i] < rend) {
172: row = im[i] - rstart;
173: lastcol1 = -1;
174: shift1 = a->sliidx[row / sliceheight] + (row % sliceheight); /* starting index of the row */
175: cp1 = PetscSafePointerPlusOffset(a->colidx, shift1);
176: vp1 = PetscSafePointerPlusOffset(a->val, shift1);
177: nrow1 = a->rlen[row];
178: low1 = 0;
179: high1 = nrow1;
180: lastcol2 = -1;
181: shift2 = b->sliidx[row / sliceheight] + (row % sliceheight); /* starting index of the row */
182: cp2 = PetscSafePointerPlusOffset(b->colidx, shift2);
183: vp2 = PetscSafePointerPlusOffset(b->val, shift2);
184: nrow2 = b->rlen[row];
185: low2 = 0;
186: high2 = nrow2;
188: for (j = 0; j < n; j++) {
189: if (roworiented) value = v[i * n + j];
190: else value = v[i + j * m];
191: if (ignorezeroentries && value == 0.0 && (addv == ADD_VALUES)) continue;
192: if (in[j] >= cstart && in[j] < cend) {
193: col = in[j] - cstart;
194: MatSetValue_SeqSELL_Private(A, row, col, value, addv, im[i], in[j], cp1, vp1, lastcol1, low1, high1); /* set one value */
195: #if PetscDefined(HAVE_CUDA)
196: if (A->offloadmask != PETSC_OFFLOAD_UNALLOCATED && found) A->offloadmask = PETSC_OFFLOAD_CPU;
197: #endif
198: } else if (in[j] < 0) {
199: continue;
200: } else {
201: PetscCheck(in[j] < mat->cmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Column too large: col %" PetscInt_FMT " max %" PetscInt_FMT, in[j], mat->cmap->N - 1);
202: if (mat->was_assembled) {
203: if (!sell->colmap) PetscCall(MatCreateColmap_MPISELL_Private(mat));
204: #if PetscDefined(USE_CTABLE)
205: PetscCall(PetscHMapIGetWithDefault(sell->colmap, in[j] + 1, 0, &col));
206: col--;
207: #else
208: col = sell->colmap[in[j]] - 1;
209: #endif
210: if (col < 0 && !((Mat_SeqSELL *)sell->B->data)->nonew) {
211: PetscCall(MatDisAssemble_MPISELL(mat));
212: col = in[j];
213: /* Reinitialize the variables required by MatSetValues_SeqSELL_B_Private() */
214: B = sell->B;
215: b = (Mat_SeqSELL *)B->data;
216: shift2 = b->sliidx[row / sliceheight] + (row % sliceheight); /* starting index of the row */
217: cp2 = b->colidx + shift2;
218: vp2 = b->val + shift2;
219: nrow2 = b->rlen[row];
220: low2 = 0;
221: high2 = nrow2;
222: found = PETSC_FALSE;
223: } else {
224: PetscCheck(col >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Inserting a new nonzero at global row/column (%" PetscInt_FMT ", %" PetscInt_FMT ") into matrix", im[i], in[j]);
225: }
226: } else col = in[j];
227: MatSetValue_SeqSELL_Private(B, row, col, value, addv, im[i], in[j], cp2, vp2, lastcol2, low2, high2); /* set one value */
228: #if PetscDefined(HAVE_CUDA)
229: if (B->offloadmask != PETSC_OFFLOAD_UNALLOCATED && found) B->offloadmask = PETSC_OFFLOAD_CPU;
230: #endif
231: }
232: }
233: } else {
234: PetscCheck(!mat->nooffprocentries, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Setting off process row %" PetscInt_FMT " even though MatSetOption(,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE) was set", im[i]);
235: if (!sell->donotstash) {
236: mat->assembled = PETSC_FALSE;
237: if (roworiented) {
238: PetscCall(MatStashValuesRow_Private(&mat->stash, im[i], n, in, v + i * n, (PetscBool)(ignorezeroentries && (addv == ADD_VALUES))));
239: } else {
240: PetscCall(MatStashValuesCol_Private(&mat->stash, im[i], n, in, v + i, m, (PetscBool)(ignorezeroentries && (addv == ADD_VALUES))));
241: }
242: }
243: }
244: }
245: PetscFunctionReturn(PETSC_SUCCESS);
246: }
248: static PetscErrorCode MatGetValues_MPISELL(Mat mat, PetscInt m, const PetscInt idxm[], PetscInt n, const PetscInt idxn[], PetscScalar v[])
249: {
250: Mat_MPISELL *sell = (Mat_MPISELL *)mat->data;
251: PetscInt i, j, rstart = mat->rmap->rstart, rend = mat->rmap->rend;
252: PetscInt cstart = mat->cmap->rstart, cend = mat->cmap->rend, row, col;
253: PetscBool roworiented = sell->roworiented;
254: PetscScalar *value;
256: PetscFunctionBegin;
257: for (i = 0; i < m; i++) {
258: if (idxm[i] < 0) continue; /* negative row */
259: PetscCheck(idxm[i] < mat->rmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, idxm[i], mat->rmap->N - 1);
260: PetscCheck(idxm[i] >= rstart && idxm[i] < rend, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only local values currently supported");
261: row = idxm[i] - rstart;
262: for (j = 0; j < n; j++) {
263: if (idxn[j] < 0) continue; /* negative column */
264: PetscCheck(idxn[j] < mat->cmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Column too large: col %" PetscInt_FMT " max %" PetscInt_FMT, idxn[j], mat->cmap->N - 1);
265: value = roworiented ? &v[j + i * n] : &v[i + j * m];
266: if (idxn[j] >= cstart && idxn[j] < cend) {
267: col = idxn[j] - cstart;
268: PetscCall(MatGetValues(sell->A, 1, &row, 1, &col, value));
269: } else {
270: if (!sell->colmap) PetscCall(MatCreateColmap_MPISELL_Private(mat));
271: #if PetscDefined(USE_CTABLE)
272: PetscCall(PetscHMapIGetWithDefault(sell->colmap, idxn[j] + 1, 0, &col));
273: col--;
274: #else
275: col = sell->colmap[idxn[j]] - 1;
276: #endif
277: if (col < 0 || sell->garray[col] != idxn[j]) *value = 0.0;
278: else PetscCall(MatGetValues(sell->B, 1, &row, 1, &col, value));
279: }
280: }
281: }
282: PetscFunctionReturn(PETSC_SUCCESS);
283: }
285: static PetscErrorCode MatAssemblyBegin_MPISELL(Mat mat, MatAssemblyType mode)
286: {
287: Mat_MPISELL *sell = (Mat_MPISELL *)mat->data;
288: PetscInt nstash, reallocs;
290: PetscFunctionBegin;
291: if (sell->donotstash || mat->nooffprocentries) PetscFunctionReturn(PETSC_SUCCESS);
293: PetscCall(MatStashScatterBegin_Private(mat, &mat->stash, mat->rmap->range));
294: PetscCall(MatStashGetInfo_Private(&mat->stash, &nstash, &reallocs));
295: PetscCall(PetscInfo(sell->A, "Stash has %" PetscInt_FMT " entries, uses %" PetscInt_FMT " mallocs.\n", nstash, reallocs));
296: PetscFunctionReturn(PETSC_SUCCESS);
297: }
299: PetscErrorCode MatAssemblyEnd_MPISELL(Mat mat, MatAssemblyType mode)
300: {
301: Mat_MPISELL *sell = (Mat_MPISELL *)mat->data;
302: PetscMPIInt n;
303: PetscInt i, flg;
304: PetscInt *row, *col;
305: PetscScalar *val;
306: PetscBool all_assembled;
307: /* do not use 'b = (Mat_SeqSELL*)sell->B->data' as B can be reset in disassembly */
308: PetscFunctionBegin;
309: if (!sell->donotstash && !mat->nooffprocentries) {
310: while (1) {
311: PetscCall(MatStashScatterGetMesg_Private(&mat->stash, &n, &row, &col, &val, &flg));
312: if (!flg) break;
314: for (i = 0; i < n; i++) { /* assemble one by one */
315: PetscCall(MatSetValues_MPISELL(mat, 1, row + i, 1, col + i, val + i, mat->insertmode));
316: }
317: }
318: PetscCall(MatStashScatterEnd_Private(&mat->stash));
319: }
320: #if PetscDefined(HAVE_CUDA)
321: if (mat->offloadmask == PETSC_OFFLOAD_CPU) sell->A->offloadmask = PETSC_OFFLOAD_CPU;
322: #endif
323: PetscCall(MatAssemblyBegin(sell->A, mode));
324: PetscCall(MatAssemblyEnd(sell->A, mode));
326: /*
327: determine if any process has disassembled, if so we must
328: also disassemble ourselves, in order that we may reassemble.
329: */
330: /*
331: if nonzero structure of submatrix B cannot change then we know that
332: no process disassembled thus we can skip this stuff
333: */
334: if (!((Mat_SeqSELL *)sell->B->data)->nonew) {
335: PetscCallMPI(MPIU_Allreduce(&mat->was_assembled, &all_assembled, 1, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)mat)));
336: if (mat->was_assembled && !all_assembled) PetscCall(MatDisAssemble_MPISELL(mat));
337: }
338: if (!mat->was_assembled && mode == MAT_FINAL_ASSEMBLY) PetscCall(MatSetUpMultiply_MPISELL(mat));
339: #if PetscDefined(HAVE_CUDA)
340: if (mat->offloadmask == PETSC_OFFLOAD_CPU && sell->B->offloadmask != PETSC_OFFLOAD_UNALLOCATED) sell->B->offloadmask = PETSC_OFFLOAD_CPU;
341: #endif
342: PetscCall(MatAssemblyBegin(sell->B, mode));
343: PetscCall(MatAssemblyEnd(sell->B, mode));
344: PetscCall(PetscFree2(sell->rowvalues, sell->rowindices));
345: sell->rowvalues = NULL;
346: PetscCall(VecDestroy(&sell->diag));
348: /* if no new nonzero locations are allowed in matrix then only set the matrix state the first time through */
349: if ((!mat->was_assembled && mode == MAT_FINAL_ASSEMBLY) || !((Mat_SeqSELL *)sell->A->data)->nonew) {
350: mat->nonzerostate = sell->A->nonzerostate + sell->B->nonzerostate;
351: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &mat->nonzerostate, 1, MPIU_INT64, MPI_SUM, PetscObjectComm((PetscObject)mat)));
352: }
353: #if PetscDefined(HAVE_CUDA)
354: mat->offloadmask = PETSC_OFFLOAD_BOTH;
355: #endif
356: PetscFunctionReturn(PETSC_SUCCESS);
357: }
359: static PetscErrorCode MatZeroEntries_MPISELL(Mat A)
360: {
361: Mat_MPISELL *l = (Mat_MPISELL *)A->data;
363: PetscFunctionBegin;
364: PetscCall(MatZeroEntries(l->A));
365: PetscCall(MatZeroEntries(l->B));
366: PetscFunctionReturn(PETSC_SUCCESS);
367: }
369: static PetscErrorCode MatMult_MPISELL(Mat A, Vec xx, Vec yy)
370: {
371: Mat_MPISELL *a = (Mat_MPISELL *)A->data;
372: PetscInt nt;
374: PetscFunctionBegin;
375: PetscCall(VecGetLocalSize(xx, &nt));
376: PetscCheck(nt == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Incompatible partition of A (%" PetscInt_FMT ") and xx (%" PetscInt_FMT ")", A->cmap->n, nt);
377: PetscCall(VecScatterBegin(a->Mvctx, xx, a->lvec, INSERT_VALUES, SCATTER_FORWARD));
378: PetscUseTypeMethod(a->A, mult, xx, yy);
379: PetscCall(VecScatterEnd(a->Mvctx, xx, a->lvec, INSERT_VALUES, SCATTER_FORWARD));
380: PetscUseTypeMethod(a->B, multadd, a->lvec, yy, yy);
381: PetscFunctionReturn(PETSC_SUCCESS);
382: }
384: static PetscErrorCode MatGetMultPetscSF_MPISELL(Mat A, PetscSF *sf)
385: {
386: Mat_MPISELL *a = (Mat_MPISELL *)A->data;
388: PetscFunctionBegin;
389: *sf = a->Mvctx;
390: PetscFunctionReturn(PETSC_SUCCESS);
391: }
393: static PetscErrorCode MatMultDiagonalBlock_MPISELL(Mat A, Vec bb, Vec xx)
394: {
395: Mat_MPISELL *a = (Mat_MPISELL *)A->data;
397: PetscFunctionBegin;
398: PetscCall(MatMultDiagonalBlock(a->A, bb, xx));
399: PetscFunctionReturn(PETSC_SUCCESS);
400: }
402: static PetscErrorCode MatMultAdd_MPISELL(Mat A, Vec xx, Vec yy, Vec zz)
403: {
404: Mat_MPISELL *a = (Mat_MPISELL *)A->data;
406: PetscFunctionBegin;
407: PetscCall(VecScatterBegin(a->Mvctx, xx, a->lvec, INSERT_VALUES, SCATTER_FORWARD));
408: PetscUseTypeMethod(a->A, multadd, xx, yy, zz);
409: PetscCall(VecScatterEnd(a->Mvctx, xx, a->lvec, INSERT_VALUES, SCATTER_FORWARD));
410: PetscUseTypeMethod(a->B, multadd, a->lvec, zz, zz);
411: PetscFunctionReturn(PETSC_SUCCESS);
412: }
414: static PetscErrorCode MatMultTranspose_MPISELL(Mat A, Vec xx, Vec yy)
415: {
416: Mat_MPISELL *a = (Mat_MPISELL *)A->data;
418: PetscFunctionBegin;
419: /* do nondiagonal part */
420: PetscUseTypeMethod(a->B, multtranspose, xx, a->lvec);
421: /* do local part */
422: PetscUseTypeMethod(a->A, multtranspose, xx, yy);
423: /* add partial results together */
424: PetscCall(VecScatterBegin(a->Mvctx, a->lvec, yy, ADD_VALUES, SCATTER_REVERSE));
425: PetscCall(VecScatterEnd(a->Mvctx, a->lvec, yy, ADD_VALUES, SCATTER_REVERSE));
426: PetscFunctionReturn(PETSC_SUCCESS);
427: }
429: static PetscErrorCode MatIsTranspose_MPISELL(Mat Amat, Mat Bmat, PetscReal tol, PetscBool *f)
430: {
431: MPI_Comm comm;
432: Mat_MPISELL *Asell = (Mat_MPISELL *)Amat->data, *Bsell;
433: Mat Adia = Asell->A, Bdia, Aoff, Boff, *Aoffs, *Boffs;
434: IS Me, Notme;
435: PetscInt M, N, first, last, *notme, i;
436: PetscMPIInt size;
438: PetscFunctionBegin;
439: /* Easy test: symmetric diagonal block */
440: Bsell = (Mat_MPISELL *)Bmat->data;
441: Bdia = Bsell->A;
442: PetscCall(MatIsTranspose(Adia, Bdia, tol, f));
443: if (!*f) PetscFunctionReturn(PETSC_SUCCESS);
444: PetscCall(PetscObjectGetComm((PetscObject)Amat, &comm));
445: PetscCallMPI(MPI_Comm_size(comm, &size));
446: if (size == 1) PetscFunctionReturn(PETSC_SUCCESS);
448: /* Hard test: off-diagonal block. This takes a MatCreateSubMatrix. */
449: PetscCall(MatGetSize(Amat, &M, &N));
450: PetscCall(MatGetOwnershipRange(Amat, &first, &last));
451: PetscCall(PetscMalloc1(N - last + first, ¬me));
452: for (i = 0; i < first; i++) notme[i] = i;
453: for (i = last; i < M; i++) notme[i - last + first] = i;
454: PetscCall(ISCreateGeneral(MPI_COMM_SELF, N - last + first, notme, PETSC_COPY_VALUES, &Notme));
455: PetscCall(ISCreateStride(MPI_COMM_SELF, last - first, first, 1, &Me));
456: PetscCall(MatCreateSubMatrices(Amat, 1, &Me, &Notme, MAT_INITIAL_MATRIX, &Aoffs));
457: Aoff = Aoffs[0];
458: PetscCall(MatCreateSubMatrices(Bmat, 1, &Notme, &Me, MAT_INITIAL_MATRIX, &Boffs));
459: Boff = Boffs[0];
460: PetscCall(MatIsTranspose(Aoff, Boff, tol, f));
461: PetscCall(MatDestroyMatrices(1, &Aoffs));
462: PetscCall(MatDestroyMatrices(1, &Boffs));
463: PetscCall(ISDestroy(&Me));
464: PetscCall(ISDestroy(&Notme));
465: PetscCall(PetscFree(notme));
466: PetscFunctionReturn(PETSC_SUCCESS);
467: }
469: static PetscErrorCode MatMultTransposeAdd_MPISELL(Mat A, Vec xx, Vec yy, Vec zz)
470: {
471: Mat_MPISELL *a = (Mat_MPISELL *)A->data;
473: PetscFunctionBegin;
474: /* do nondiagonal part */
475: PetscUseTypeMethod(a->B, multtranspose, xx, a->lvec);
476: /* do local part */
477: PetscUseTypeMethod(a->A, multtransposeadd, xx, yy, zz);
478: /* add partial results together */
479: PetscCall(VecScatterBegin(a->Mvctx, a->lvec, zz, ADD_VALUES, SCATTER_REVERSE));
480: PetscCall(VecScatterEnd(a->Mvctx, a->lvec, zz, ADD_VALUES, SCATTER_REVERSE));
481: PetscFunctionReturn(PETSC_SUCCESS);
482: }
484: /*
485: This only works correctly for square matrices where the subblock A->A is the
486: diagonal block
487: */
488: static PetscErrorCode MatGetDiagonal_MPISELL(Mat A, Vec v)
489: {
490: Mat_MPISELL *a = (Mat_MPISELL *)A->data;
492: PetscFunctionBegin;
493: PetscCheck(A->rmap->N == A->cmap->N, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Supports only square matrix where A->A is diag block");
494: PetscCheck(A->rmap->rstart == A->cmap->rstart && A->rmap->rend == A->cmap->rend, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "row partition must equal col partition");
495: PetscCall(MatGetDiagonal(a->A, v));
496: PetscFunctionReturn(PETSC_SUCCESS);
497: }
499: static PetscErrorCode MatScale_MPISELL(Mat A, PetscScalar aa)
500: {
501: Mat_MPISELL *a = (Mat_MPISELL *)A->data;
503: PetscFunctionBegin;
504: PetscCall(MatScale(a->A, aa));
505: PetscCall(MatScale(a->B, aa));
506: PetscFunctionReturn(PETSC_SUCCESS);
507: }
509: PetscErrorCode MatDestroy_MPISELL(Mat mat)
510: {
511: Mat_MPISELL *sell = (Mat_MPISELL *)mat->data;
513: PetscFunctionBegin;
514: PetscCall(PetscLogObjectState((PetscObject)mat, "Rows=%" PetscInt_FMT ", Cols=%" PetscInt_FMT, mat->rmap->N, mat->cmap->N));
515: PetscCall(MatStashDestroy_Private(&mat->stash));
516: PetscCall(VecDestroy(&sell->diag));
517: PetscCall(MatDestroy(&sell->A));
518: PetscCall(MatDestroy(&sell->B));
519: #if PetscDefined(USE_CTABLE)
520: PetscCall(PetscHMapIDestroy(&sell->colmap));
521: #else
522: PetscCall(PetscFree(sell->colmap));
523: #endif
524: PetscCall(PetscFree(sell->garray));
525: PetscCall(VecDestroy(&sell->lvec));
526: PetscCall(VecScatterDestroy(&sell->Mvctx));
527: PetscCall(PetscFree2(sell->rowvalues, sell->rowindices));
528: PetscCall(PetscFree(sell->ld));
529: PetscCall(PetscFree(mat->data));
531: PetscCall(PetscObjectChangeTypeName((PetscObject)mat, NULL));
532: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatStoreValues_C", NULL));
533: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatRetrieveValues_C", NULL));
534: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatIsTranspose_C", NULL));
535: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMPISELLSetPreallocation_C", NULL));
536: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_mpisell_mpiaij_C", NULL));
537: #if PetscDefined(HAVE_CUDA)
538: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_mpisell_mpisellcuda_C", NULL));
539: #endif
540: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDiagonalScaleLocal_C", NULL));
541: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatGetMultPetscSF_C", NULL));
542: PetscFunctionReturn(PETSC_SUCCESS);
543: }
545: #include <petscdraw.h>
546: static PetscErrorCode MatView_MPISELL_ASCIIorDraworSocket(Mat mat, PetscViewer viewer)
547: {
548: Mat_MPISELL *sell = (Mat_MPISELL *)mat->data;
549: PetscMPIInt rank = sell->rank, size = sell->size;
550: PetscBool isdraw, isascii, isbinary;
551: PetscViewer sviewer;
552: PetscViewerFormat format;
554: PetscFunctionBegin;
555: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
556: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
557: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
558: if (isascii) {
559: PetscCall(PetscViewerGetFormat(viewer, &format));
560: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
561: MatInfo info;
562: PetscInt *inodes;
564: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)mat), &rank));
565: PetscCall(MatGetInfo(mat, MAT_LOCAL, &info));
566: PetscCall(MatInodeGetInodeSizes(sell->A, NULL, &inodes, NULL));
567: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
568: if (!inodes) {
569: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Local rows %" PetscInt_FMT " nz %" PetscInt_FMT " nz alloced %" PetscInt_FMT " mem %" PetscInt_FMT ", not using I-node routines\n", rank, mat->rmap->n, (PetscInt)info.nz_used,
570: (PetscInt)info.nz_allocated, (PetscInt)info.memory));
571: } else {
572: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Local rows %" PetscInt_FMT " nz %" PetscInt_FMT " nz alloced %" PetscInt_FMT " mem %" PetscInt_FMT ", using I-node routines\n", rank, mat->rmap->n, (PetscInt)info.nz_used,
573: (PetscInt)info.nz_allocated, (PetscInt)info.memory));
574: }
575: PetscCall(MatGetInfo(sell->A, MAT_LOCAL, &info));
576: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] on-diagonal part: nz %" PetscInt_FMT " \n", rank, (PetscInt)info.nz_used));
577: PetscCall(MatGetInfo(sell->B, MAT_LOCAL, &info));
578: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] off-diagonal part: nz %" PetscInt_FMT " \n", rank, (PetscInt)info.nz_used));
579: PetscCall(PetscViewerFlush(viewer));
580: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
581: PetscCall(PetscViewerASCIIPrintf(viewer, "Information on VecScatter used in matrix-vector product: \n"));
582: PetscCall(VecScatterView(sell->Mvctx, viewer));
583: PetscFunctionReturn(PETSC_SUCCESS);
584: } else if (format == PETSC_VIEWER_ASCII_INFO) {
585: PetscInt inodecount, inodelimit, *inodes;
586: PetscCall(MatInodeGetInodeSizes(sell->A, &inodecount, &inodes, &inodelimit));
587: if (inodes) {
588: PetscCall(PetscViewerASCIIPrintf(viewer, "using I-node (on process 0) routines: found %" PetscInt_FMT " nodes, limit used is %" PetscInt_FMT "\n", inodecount, inodelimit));
589: } else {
590: PetscCall(PetscViewerASCIIPrintf(viewer, "not using I-node (on process 0) routines\n"));
591: }
592: PetscFunctionReturn(PETSC_SUCCESS);
593: } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) {
594: PetscFunctionReturn(PETSC_SUCCESS);
595: }
596: } else if (isbinary) {
597: if (size == 1) {
598: PetscCall(PetscObjectSetName((PetscObject)sell->A, ((PetscObject)mat)->name));
599: PetscCall(MatView(sell->A, viewer));
600: } else {
601: /* PetscCall(MatView_MPISELL_Binary(mat,viewer)); */
602: }
603: PetscFunctionReturn(PETSC_SUCCESS);
604: } else if (isdraw) {
605: PetscDraw draw;
606: PetscBool isnull;
607: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
608: PetscCall(PetscDrawIsNull(draw, &isnull));
609: if (isnull) PetscFunctionReturn(PETSC_SUCCESS);
610: }
612: {
613: /* assemble the entire matrix onto first processor. */
614: Mat A;
615: Mat_SeqSELL *Aloc;
616: PetscInt M = mat->rmap->N, N = mat->cmap->N, *acolidx, row, col, i, j;
617: MatScalar *aval;
618: PetscBool isnonzero;
620: PetscCall(MatCreate(PetscObjectComm((PetscObject)mat), &A));
621: if (rank == 0) {
622: PetscCall(MatSetSizes(A, M, N, M, N));
623: } else {
624: PetscCall(MatSetSizes(A, 0, 0, M, N));
625: }
626: /* This is just a temporary matrix, so explicitly using MATMPISELL is probably best */
627: PetscCall(MatSetType(A, MATMPISELL));
628: PetscCall(MatMPISELLSetPreallocation(A, 0, NULL, 0, NULL));
629: PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_FALSE));
631: /* copy over the A part */
632: Aloc = (Mat_SeqSELL *)sell->A->data;
633: acolidx = Aloc->colidx;
634: aval = Aloc->val;
635: for (i = 0; i < Aloc->totalslices; i++) { /* loop over slices */
636: for (j = Aloc->sliidx[i]; j < Aloc->sliidx[i + 1]; j++) {
637: isnonzero = (PetscBool)((j - Aloc->sliidx[i]) / Aloc->sliceheight < Aloc->rlen[i * Aloc->sliceheight + j % Aloc->sliceheight]);
638: if (isnonzero) { /* check the mask bit */
639: row = i * Aloc->sliceheight + j % Aloc->sliceheight + mat->rmap->rstart;
640: col = *acolidx + mat->rmap->rstart;
641: PetscCall(MatSetValues(A, 1, &row, 1, &col, aval, INSERT_VALUES));
642: }
643: aval++;
644: acolidx++;
645: }
646: }
648: /* copy over the B part */
649: Aloc = (Mat_SeqSELL *)sell->B->data;
650: acolidx = Aloc->colidx;
651: aval = Aloc->val;
652: for (i = 0; i < Aloc->totalslices; i++) {
653: for (j = Aloc->sliidx[i]; j < Aloc->sliidx[i + 1]; j++) {
654: isnonzero = (PetscBool)((j - Aloc->sliidx[i]) / Aloc->sliceheight < Aloc->rlen[i * Aloc->sliceheight + j % Aloc->sliceheight]);
655: if (isnonzero) {
656: row = i * Aloc->sliceheight + j % Aloc->sliceheight + mat->rmap->rstart;
657: col = sell->garray[*acolidx];
658: PetscCall(MatSetValues(A, 1, &row, 1, &col, aval, INSERT_VALUES));
659: }
660: aval++;
661: acolidx++;
662: }
663: }
665: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
666: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
667: /*
668: Everyone has to call to draw the matrix since the graphics waits are
669: synchronized across all processors that share the PetscDraw object
670: */
671: PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
672: if (rank == 0) {
673: PetscCall(PetscObjectSetName((PetscObject)((Mat_MPISELL *)A->data)->A, ((PetscObject)mat)->name));
674: PetscCall(MatView_SeqSELL(((Mat_MPISELL *)A->data)->A, sviewer));
675: }
676: PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
677: PetscCall(MatDestroy(&A));
678: }
679: PetscFunctionReturn(PETSC_SUCCESS);
680: }
682: static PetscErrorCode MatView_MPISELL(Mat mat, PetscViewer viewer)
683: {
684: PetscBool isascii, isdraw, issocket, isbinary;
686: PetscFunctionBegin;
687: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
688: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
689: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
690: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSOCKET, &issocket));
691: if (isascii || isdraw || isbinary || issocket) PetscCall(MatView_MPISELL_ASCIIorDraworSocket(mat, viewer));
692: PetscFunctionReturn(PETSC_SUCCESS);
693: }
695: static PetscErrorCode MatGetGhosts_MPISELL(Mat mat, PetscInt *nghosts, const PetscInt *ghosts[])
696: {
697: Mat_MPISELL *sell = (Mat_MPISELL *)mat->data;
699: PetscFunctionBegin;
700: PetscCall(MatGetSize(sell->B, NULL, nghosts));
701: if (ghosts) *ghosts = sell->garray;
702: PetscFunctionReturn(PETSC_SUCCESS);
703: }
705: static PetscErrorCode MatGetInfo_MPISELL(Mat matin, MatInfoType flag, MatInfo *info)
706: {
707: Mat_MPISELL *mat = (Mat_MPISELL *)matin->data;
708: Mat A = mat->A, B = mat->B;
709: PetscLogDouble irecv[5];
711: PetscFunctionBegin;
712: info->block_size = 1.0;
713: PetscCall(MatGetInfo(A, MAT_LOCAL, info));
715: irecv[0] = info->nz_used;
716: irecv[1] = info->nz_allocated;
717: irecv[2] = info->nz_unneeded;
718: irecv[3] = info->memory;
719: irecv[4] = info->mallocs;
721: PetscCall(MatGetInfo(B, MAT_LOCAL, info));
723: irecv[0] += info->nz_used;
724: irecv[1] += info->nz_allocated;
725: irecv[2] += info->nz_unneeded;
726: irecv[3] += info->memory;
727: irecv[4] += info->mallocs;
728: if (flag == MAT_LOCAL) {
729: info->nz_used = irecv[0];
730: info->nz_allocated = irecv[1];
731: info->nz_unneeded = irecv[2];
732: info->memory = irecv[3];
733: info->mallocs = irecv[4];
734: } else if (flag == MAT_GLOBAL_MAX) {
735: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, irecv, 5, MPIU_PETSCLOGDOUBLE, MPI_MAX, PetscObjectComm((PetscObject)matin)));
737: info->nz_used = irecv[0];
738: info->nz_allocated = irecv[1];
739: info->nz_unneeded = irecv[2];
740: info->memory = irecv[3];
741: info->mallocs = irecv[4];
742: } else if (flag == MAT_GLOBAL_SUM) {
743: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, irecv, 5, MPIU_PETSCLOGDOUBLE, MPI_SUM, PetscObjectComm((PetscObject)matin)));
745: info->nz_used = irecv[0];
746: info->nz_allocated = irecv[1];
747: info->nz_unneeded = irecv[2];
748: info->memory = irecv[3];
749: info->mallocs = irecv[4];
750: }
751: info->fill_ratio_given = 0; /* no parallel LU/ILU/Cholesky */
752: info->fill_ratio_needed = 0;
753: info->factor_mallocs = 0;
754: PetscFunctionReturn(PETSC_SUCCESS);
755: }
757: static PetscErrorCode MatSetOption_MPISELL(Mat A, MatOption op, PetscBool flg)
758: {
759: Mat_MPISELL *a = (Mat_MPISELL *)A->data;
761: PetscFunctionBegin;
762: switch (op) {
763: case MAT_NEW_NONZERO_LOCATIONS:
764: case MAT_NEW_NONZERO_ALLOCATION_ERR:
765: case MAT_UNUSED_NONZERO_LOCATION_ERR:
766: case MAT_KEEP_NONZERO_PATTERN:
767: case MAT_NEW_NONZERO_LOCATION_ERR:
768: case MAT_USE_INODES:
769: case MAT_IGNORE_ZERO_ENTRIES:
770: MatCheckPreallocated(A, 1);
771: PetscCall(MatSetOption(a->A, op, flg));
772: PetscCall(MatSetOption(a->B, op, flg));
773: break;
774: case MAT_ROW_ORIENTED:
775: MatCheckPreallocated(A, 1);
776: a->roworiented = flg;
778: PetscCall(MatSetOption(a->A, op, flg));
779: PetscCall(MatSetOption(a->B, op, flg));
780: break;
781: case MAT_IGNORE_OFF_PROC_ENTRIES:
782: a->donotstash = flg;
783: break;
784: case MAT_SYMMETRIC:
785: MatCheckPreallocated(A, 1);
786: PetscCall(MatSetOption(a->A, op, flg));
787: break;
788: case MAT_STRUCTURALLY_SYMMETRIC:
789: MatCheckPreallocated(A, 1);
790: PetscCall(MatSetOption(a->A, op, flg));
791: break;
792: case MAT_HERMITIAN:
793: MatCheckPreallocated(A, 1);
794: PetscCall(MatSetOption(a->A, op, flg));
795: break;
796: case MAT_SYMMETRY_ETERNAL:
797: MatCheckPreallocated(A, 1);
798: PetscCall(MatSetOption(a->A, op, flg));
799: break;
800: case MAT_STRUCTURAL_SYMMETRY_ETERNAL:
801: MatCheckPreallocated(A, 1);
802: PetscCall(MatSetOption(a->A, op, flg));
803: break;
804: default:
805: break;
806: }
807: PetscFunctionReturn(PETSC_SUCCESS);
808: }
810: static PetscErrorCode MatDiagonalScale_MPISELL(Mat mat, Vec ll, Vec rr)
811: {
812: Mat_MPISELL *sell = (Mat_MPISELL *)mat->data;
813: Mat a = sell->A, b = sell->B;
814: PetscInt s1, s2, s3;
816: PetscFunctionBegin;
817: PetscCall(MatGetLocalSize(mat, &s2, &s3));
818: if (rr) {
819: PetscCall(VecGetLocalSize(rr, &s1));
820: PetscCheck(s1 == s3, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "right vector non-conforming local size");
821: /* Overlap communication with computation. */
822: PetscCall(VecScatterBegin(sell->Mvctx, rr, sell->lvec, INSERT_VALUES, SCATTER_FORWARD));
823: }
824: if (ll) {
825: PetscCall(VecGetLocalSize(ll, &s1));
826: PetscCheck(s1 == s2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "left vector non-conforming local size");
827: PetscUseTypeMethod(b, diagonalscale, ll, NULL);
828: }
829: /* scale the diagonal block */
830: PetscUseTypeMethod(a, diagonalscale, ll, rr);
832: if (rr) {
833: /* Do a scatter end and then right scale the off-diagonal block */
834: PetscCall(VecScatterEnd(sell->Mvctx, rr, sell->lvec, INSERT_VALUES, SCATTER_FORWARD));
835: PetscUseTypeMethod(b, diagonalscale, NULL, sell->lvec);
836: }
837: PetscFunctionReturn(PETSC_SUCCESS);
838: }
840: static PetscErrorCode MatSetUnfactored_MPISELL(Mat A)
841: {
842: Mat_MPISELL *a = (Mat_MPISELL *)A->data;
844: PetscFunctionBegin;
845: PetscCall(MatSetUnfactored(a->A));
846: PetscFunctionReturn(PETSC_SUCCESS);
847: }
849: static PetscErrorCode MatEqual_MPISELL(Mat A, Mat B, PetscBool *flag)
850: {
851: Mat_MPISELL *matB = (Mat_MPISELL *)B->data, *matA = (Mat_MPISELL *)A->data;
852: Mat a, b, c, d;
854: PetscFunctionBegin;
855: a = matA->A;
856: b = matA->B;
857: c = matB->A;
858: d = matB->B;
860: PetscCall(MatEqual(a, c, flag));
861: if (*flag) PetscCall(MatEqual(b, d, flag));
862: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, flag, 1, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A)));
863: PetscFunctionReturn(PETSC_SUCCESS);
864: }
866: static PetscErrorCode MatCopy_MPISELL(Mat A, Mat B, MatStructure str)
867: {
868: Mat_MPISELL *a = (Mat_MPISELL *)A->data;
869: Mat_MPISELL *b = (Mat_MPISELL *)B->data;
871: PetscFunctionBegin;
872: /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */
873: if ((str != SAME_NONZERO_PATTERN) || (A->ops->copy != B->ops->copy)) {
874: /* because of the column compression in the off-processor part of the matrix a->B,
875: the number of columns in a->B and b->B may be different, hence we cannot call
876: the MatCopy() directly on the two parts. If need be, we can provide a more
877: efficient copy than the MatCopy_Basic() by first uncompressing the a->B matrices
878: then copying the submatrices */
879: PetscCall(MatCopy_Basic(A, B, str));
880: } else {
881: PetscCall(MatCopy(a->A, b->A, str));
882: PetscCall(MatCopy(a->B, b->B, str));
883: }
884: PetscFunctionReturn(PETSC_SUCCESS);
885: }
887: static PetscErrorCode MatSetUp_MPISELL(Mat A)
888: {
889: PetscFunctionBegin;
890: PetscCall(MatMPISELLSetPreallocation(A, PETSC_DEFAULT, NULL, PETSC_DEFAULT, NULL));
891: PetscFunctionReturn(PETSC_SUCCESS);
892: }
894: static PetscErrorCode MatConjugate_MPISELL(Mat mat)
895: {
896: Mat_MPISELL *sell = (Mat_MPISELL *)mat->data;
898: PetscFunctionBegin;
899: PetscCall(MatConjugate_SeqSELL(sell->A));
900: PetscCall(MatConjugate_SeqSELL(sell->B));
901: PetscFunctionReturn(PETSC_SUCCESS);
902: }
904: static PetscErrorCode MatInvertBlockDiagonal_MPISELL(Mat A, const PetscScalar **values)
905: {
906: Mat_MPISELL *a = (Mat_MPISELL *)A->data;
908: PetscFunctionBegin;
909: PetscCall(MatInvertBlockDiagonal(a->A, values));
910: A->factorerrortype = a->A->factorerrortype;
911: PetscFunctionReturn(PETSC_SUCCESS);
912: }
914: static PetscErrorCode MatSetRandom_MPISELL(Mat x, PetscRandom rctx)
915: {
916: Mat_MPISELL *sell = (Mat_MPISELL *)x->data;
918: PetscFunctionBegin;
919: PetscCall(MatSetRandom(sell->A, rctx));
920: PetscCall(MatSetRandom(sell->B, rctx));
921: PetscCall(MatAssemblyBegin(x, MAT_FINAL_ASSEMBLY));
922: PetscCall(MatAssemblyEnd(x, MAT_FINAL_ASSEMBLY));
923: PetscFunctionReturn(PETSC_SUCCESS);
924: }
926: static PetscErrorCode MatSetFromOptions_MPISELL(Mat A, PetscOptionItems PetscOptionsObject)
927: {
928: PetscFunctionBegin;
929: PetscOptionsHeadBegin(PetscOptionsObject, "MPISELL options");
930: PetscOptionsHeadEnd();
931: PetscFunctionReturn(PETSC_SUCCESS);
932: }
934: static PetscErrorCode MatShift_MPISELL(Mat Y, PetscScalar a)
935: {
936: Mat_MPISELL *msell = (Mat_MPISELL *)Y->data;
937: Mat_SeqSELL *sell = (Mat_SeqSELL *)msell->A->data;
939: PetscFunctionBegin;
940: if (!Y->preallocated) {
941: PetscCall(MatMPISELLSetPreallocation(Y, 1, NULL, 0, NULL));
942: } else if (!sell->nz) {
943: PetscInt nonew = sell->nonew;
944: PetscCall(MatSeqSELLSetPreallocation(msell->A, 1, NULL));
945: sell->nonew = nonew;
946: }
947: PetscCall(MatShift_Basic(Y, a));
948: PetscFunctionReturn(PETSC_SUCCESS);
949: }
951: static PetscErrorCode MatGetDiagonalBlock_MPISELL(Mat A, Mat *a)
952: {
953: PetscFunctionBegin;
954: *a = ((Mat_MPISELL *)A->data)->A;
955: PetscFunctionReturn(PETSC_SUCCESS);
956: }
958: static PetscErrorCode MatStoreValues_MPISELL(Mat mat)
959: {
960: Mat_MPISELL *sell = (Mat_MPISELL *)mat->data;
962: PetscFunctionBegin;
963: PetscCall(MatStoreValues(sell->A));
964: PetscCall(MatStoreValues(sell->B));
965: PetscFunctionReturn(PETSC_SUCCESS);
966: }
968: static PetscErrorCode MatRetrieveValues_MPISELL(Mat mat)
969: {
970: Mat_MPISELL *sell = (Mat_MPISELL *)mat->data;
972: PetscFunctionBegin;
973: PetscCall(MatRetrieveValues(sell->A));
974: PetscCall(MatRetrieveValues(sell->B));
975: PetscFunctionReturn(PETSC_SUCCESS);
976: }
978: static PetscErrorCode MatMPISELLSetPreallocation_MPISELL(Mat B, PetscInt d_rlenmax, const PetscInt d_rlen[], PetscInt o_rlenmax, const PetscInt o_rlen[])
979: {
980: Mat_MPISELL *b;
982: PetscFunctionBegin;
983: PetscCall(PetscLayoutSetUp(B->rmap));
984: PetscCall(PetscLayoutSetUp(B->cmap));
985: b = (Mat_MPISELL *)B->data;
987: if (!B->preallocated) {
988: /* Explicitly create 2 MATSEQSELL matrices. */
989: PetscCall(MatCreate(PETSC_COMM_SELF, &b->A));
990: PetscCall(MatSetSizes(b->A, B->rmap->n, B->cmap->n, B->rmap->n, B->cmap->n));
991: PetscCall(MatSetBlockSizesFromMats(b->A, B, B));
992: PetscCall(MatSetType(b->A, MATSEQSELL));
993: PetscCall(MatCreate(PETSC_COMM_SELF, &b->B));
994: PetscCall(MatSetSizes(b->B, B->rmap->n, B->cmap->N, B->rmap->n, B->cmap->N));
995: PetscCall(MatSetBlockSizesFromMats(b->B, B, B));
996: PetscCall(MatSetType(b->B, MATSEQSELL));
997: }
999: PetscCall(MatSeqSELLSetPreallocation(b->A, d_rlenmax, d_rlen));
1000: PetscCall(MatSeqSELLSetPreallocation(b->B, o_rlenmax, o_rlen));
1001: B->preallocated = PETSC_TRUE;
1002: B->was_assembled = PETSC_FALSE;
1003: /*
1004: critical for MatAssemblyEnd to work.
1005: MatAssemblyBegin checks it to set up was_assembled
1006: and MatAssemblyEnd checks was_assembled to determine whether to build garray
1007: */
1008: B->assembled = PETSC_FALSE;
1009: PetscFunctionReturn(PETSC_SUCCESS);
1010: }
1012: static PetscErrorCode MatDuplicate_MPISELL(Mat matin, MatDuplicateOption cpvalues, Mat *newmat)
1013: {
1014: Mat mat;
1015: Mat_MPISELL *a, *oldmat = (Mat_MPISELL *)matin->data;
1017: PetscFunctionBegin;
1018: *newmat = NULL;
1019: PetscCall(MatCreate(PetscObjectComm((PetscObject)matin), &mat));
1020: PetscCall(MatSetSizes(mat, matin->rmap->n, matin->cmap->n, matin->rmap->N, matin->cmap->N));
1021: PetscCall(MatSetBlockSizesFromMats(mat, matin, matin));
1022: PetscCall(MatSetType(mat, ((PetscObject)matin)->type_name));
1023: a = (Mat_MPISELL *)mat->data;
1025: mat->factortype = matin->factortype;
1026: mat->assembled = PETSC_TRUE;
1027: mat->insertmode = NOT_SET_VALUES;
1028: mat->preallocated = PETSC_TRUE;
1030: a->size = oldmat->size;
1031: a->rank = oldmat->rank;
1032: a->donotstash = oldmat->donotstash;
1033: a->roworiented = oldmat->roworiented;
1034: a->rowindices = NULL;
1035: a->rowvalues = NULL;
1036: a->getrowactive = PETSC_FALSE;
1038: PetscCall(PetscLayoutReference(matin->rmap, &mat->rmap));
1039: PetscCall(PetscLayoutReference(matin->cmap, &mat->cmap));
1041: if (oldmat->colmap) {
1042: #if PetscDefined(USE_CTABLE)
1043: PetscCall(PetscHMapIDuplicate(oldmat->colmap, &a->colmap));
1044: #else
1045: PetscCall(PetscMalloc1(mat->cmap->N, &a->colmap));
1046: PetscCall(PetscArraycpy(a->colmap, oldmat->colmap, mat->cmap->N));
1047: #endif
1048: } else a->colmap = NULL;
1049: if (oldmat->garray) {
1050: PetscInt len;
1051: len = oldmat->B->cmap->n;
1052: PetscCall(PetscMalloc1(len + 1, &a->garray));
1053: if (len) PetscCall(PetscArraycpy(a->garray, oldmat->garray, len));
1054: } else a->garray = NULL;
1056: PetscCall(VecDuplicate(oldmat->lvec, &a->lvec));
1057: PetscCall(VecScatterCopy(oldmat->Mvctx, &a->Mvctx));
1058: PetscCall(MatDuplicate(oldmat->A, cpvalues, &a->A));
1059: PetscCall(MatDuplicate(oldmat->B, cpvalues, &a->B));
1060: PetscCall(PetscFunctionListDuplicate(((PetscObject)matin)->qlist, &((PetscObject)mat)->qlist));
1061: *newmat = mat;
1062: PetscFunctionReturn(PETSC_SUCCESS);
1063: }
1065: static const struct _MatOps MatOps_Values = {MatSetValues_MPISELL,
1066: NULL,
1067: NULL,
1068: MatMult_MPISELL,
1069: /* 4*/ MatMultAdd_MPISELL,
1070: MatMultTranspose_MPISELL,
1071: MatMultTransposeAdd_MPISELL,
1072: NULL,
1073: NULL,
1074: NULL,
1075: /*10*/ NULL,
1076: NULL,
1077: NULL,
1078: MatSOR_MPISELL,
1079: NULL,
1080: /*15*/ MatGetInfo_MPISELL,
1081: MatEqual_MPISELL,
1082: MatGetDiagonal_MPISELL,
1083: MatDiagonalScale_MPISELL,
1084: NULL,
1085: /*20*/ MatAssemblyBegin_MPISELL,
1086: MatAssemblyEnd_MPISELL,
1087: MatSetOption_MPISELL,
1088: MatZeroEntries_MPISELL,
1089: /*24*/ NULL,
1090: NULL,
1091: NULL,
1092: NULL,
1093: NULL,
1094: /*29*/ MatSetUp_MPISELL,
1095: NULL,
1096: NULL,
1097: MatGetDiagonalBlock_MPISELL,
1098: NULL,
1099: /*34*/ MatDuplicate_MPISELL,
1100: NULL,
1101: NULL,
1102: NULL,
1103: NULL,
1104: /*39*/ NULL,
1105: NULL,
1106: NULL,
1107: MatGetValues_MPISELL,
1108: MatCopy_MPISELL,
1109: /*44*/ NULL,
1110: MatScale_MPISELL,
1111: MatShift_MPISELL,
1112: MatDiagonalSet_MPISELL,
1113: NULL,
1114: /*49*/ MatSetRandom_MPISELL,
1115: NULL,
1116: NULL,
1117: NULL,
1118: NULL,
1119: /*54*/ MatFDColoringCreate_MPIXAIJ,
1120: NULL,
1121: MatSetUnfactored_MPISELL,
1122: NULL,
1123: NULL,
1124: /*59*/ NULL,
1125: MatDestroy_MPISELL,
1126: MatView_MPISELL,
1127: NULL,
1128: NULL,
1129: /*64*/ NULL,
1130: NULL,
1131: NULL,
1132: NULL,
1133: NULL,
1134: /*69*/ NULL,
1135: NULL,
1136: NULL,
1137: MatFDColoringApply_AIJ, /* reuse AIJ function */
1138: MatSetFromOptions_MPISELL,
1139: NULL,
1140: /*75*/ NULL,
1141: NULL,
1142: NULL,
1143: NULL,
1144: NULL,
1145: /*80*/ NULL,
1146: NULL,
1147: NULL,
1148: /*83*/ NULL,
1149: NULL,
1150: NULL,
1151: NULL,
1152: NULL,
1153: NULL,
1154: /*89*/ NULL,
1155: NULL,
1156: NULL,
1157: NULL,
1158: MatConjugate_MPISELL,
1159: /*94*/ NULL,
1160: NULL,
1161: NULL,
1162: NULL,
1163: NULL,
1164: /*99*/ NULL,
1165: NULL,
1166: NULL,
1167: NULL,
1168: NULL,
1169: /*104*/ NULL,
1170: NULL,
1171: MatGetGhosts_MPISELL,
1172: NULL,
1173: NULL,
1174: /*109*/ MatMultDiagonalBlock_MPISELL,
1175: NULL,
1176: NULL,
1177: NULL,
1178: NULL,
1179: /*114*/ NULL,
1180: NULL,
1181: MatInvertBlockDiagonal_MPISELL,
1182: NULL,
1183: /*119*/ NULL,
1184: NULL,
1185: NULL,
1186: NULL,
1187: NULL,
1188: /*124*/ NULL,
1189: NULL,
1190: NULL,
1191: NULL,
1192: MatFDColoringSetUp_MPIXAIJ,
1193: /*129*/ NULL,
1194: NULL,
1195: NULL,
1196: NULL,
1197: NULL,
1198: /*134*/ NULL,
1199: NULL,
1200: NULL,
1201: NULL,
1202: NULL,
1203: /*139*/ NULL,
1204: NULL,
1205: NULL,
1206: NULL,
1207: NULL,
1208: MatADot_Default,
1209: /*144*/ MatANorm_Default,
1210: NULL,
1211: NULL,
1212: NULL};
1214: /*@C
1215: MatMPISELLSetPreallocation - Preallocates memory for a `MATMPISELL` sparse parallel matrix in sell format.
1216: For good matrix assembly performance the user should preallocate the matrix storage by
1217: setting the parameters `d_nz` (or `d_nnz`) and `o_nz` (or `o_nnz`).
1219: Collective
1221: Input Parameters:
1222: + B - the matrix
1223: . d_nz - number of nonzeros per row in DIAGONAL portion of local submatrix
1224: (same value is used for all local rows)
1225: . d_nnz - array containing the number of nonzeros in the various rows of the
1226: DIAGONAL portion of the local submatrix (possibly different for each row)
1227: or NULL (`PETSC_NULL_INTEGER` in Fortran), if `d_nz` is used to specify the nonzero structure.
1228: The size of this array is equal to the number of local rows, i.e 'm'.
1229: For matrices that will be factored, you must leave room for (and set)
1230: the diagonal entry even if it is zero.
1231: . o_nz - number of nonzeros per row in the OFF-DIAGONAL portion of local
1232: submatrix (same value is used for all local rows).
1233: - o_nnz - array containing the number of nonzeros in the various rows of the
1234: OFF-DIAGONAL portion of the local submatrix (possibly different for
1235: each row) or NULL (`PETSC_NULL_INTEGER` in Fortran), if `o_nz` is used to specify the nonzero
1236: structure. The size of this array is equal to the number
1237: of local rows, i.e 'm'.
1239: Example usage:
1240: Consider the following 8x8 matrix with 34 non-zero values, that is
1241: assembled across 3 processors. Lets assume that proc0 owns 3 rows,
1242: proc1 owns 3 rows, proc2 owns 2 rows. This division can be shown
1243: as follows
1245: .vb
1246: 1 2 0 | 0 3 0 | 0 4
1247: Proc0 0 5 6 | 7 0 0 | 8 0
1248: 9 0 10 | 11 0 0 | 12 0
1249: -------------------------------------
1250: 13 0 14 | 15 16 17 | 0 0
1251: Proc1 0 18 0 | 19 20 21 | 0 0
1252: 0 0 0 | 22 23 0 | 24 0
1253: -------------------------------------
1254: Proc2 25 26 27 | 0 0 28 | 29 0
1255: 30 0 0 | 31 32 33 | 0 34
1256: .ve
1258: This can be represented as a collection of submatrices as
1260: .vb
1261: A B C
1262: D E F
1263: G H I
1264: .ve
1266: Where the submatrices A,B,C are owned by proc0, D,E,F are
1267: owned by proc1, G,H,I are owned by proc2.
1269: The 'm' parameters for proc0,proc1,proc2 are 3,3,2 respectively.
1270: The 'n' parameters for proc0,proc1,proc2 are 3,3,2 respectively.
1271: The 'M','N' parameters are 8,8, and have the same values on all procs.
1273: The DIAGONAL submatrices corresponding to proc0,proc1,proc2 are
1274: submatrices [A], [E], [I] respectively. The OFF-DIAGONAL submatrices
1275: corresponding to proc0,proc1,proc2 are [BC], [DF], [GH] respectively.
1276: Internally, each processor stores the DIAGONAL part, and the OFF-DIAGONAL
1277: part as `MATSEQSELL` matrices. For example, proc1 will store [E] as a `MATSEQSELL`
1278: matrix, and [DF] as another SeqSELL matrix.
1280: When `d_nz`, `o_nz` parameters are specified, `d_nz` storage elements are
1281: allocated for every row of the local DIAGONAL submatrix, and o_nz
1282: storage locations are allocated for every row of the OFF-DIAGONAL submatrix.
1283: One way to choose `d_nz` and `o_nz` is to use the maximum number of nonzeros over
1284: the local rows for each of the local DIAGONAL, and the OFF-DIAGONAL submatrices.
1285: In this case, the values of d_nz,o_nz are
1286: .vb
1287: proc0 dnz = 2, o_nz = 2
1288: proc1 dnz = 3, o_nz = 2
1289: proc2 dnz = 1, o_nz = 4
1290: .ve
1291: We are allocating m*(d_nz+o_nz) storage locations for every proc. This
1292: translates to 3*(2+2)=12 for proc0, 3*(3+2)=15 for proc1, 2*(1+4)=10
1293: for proc3. i.e we are using 12+15+10=37 storage locations to store
1294: 34 values.
1296: When `d_nnz`, `o_nnz` parameters are specified, the storage is specified
1297: for every row, corresponding to both DIAGONAL and OFF-DIAGONAL submatrices.
1298: In the above case the values for d_nnz,o_nnz are
1299: .vb
1300: proc0 d_nnz = [2,2,2] and o_nnz = [2,2,2]
1301: proc1 d_nnz = [3,3,2] and o_nnz = [2,1,1]
1302: proc2 d_nnz = [1,1] and o_nnz = [4,4]
1303: .ve
1304: Here the space allocated is according to nz (or maximum values in the nnz
1305: if nnz is provided) for DIAGONAL and OFF-DIAGONAL submatrices, i.e (2+2+3+2)*3+(1+4)*2=37
1307: Level: intermediate
1309: Notes:
1310: If the *_nnz parameter is given then the *_nz parameter is ignored
1312: The stored row and column indices begin with zero.
1314: The parallel matrix is partitioned such that the first m0 rows belong to
1315: process 0, the next m1 rows belong to process 1, the next m2 rows belong
1316: to process 2 etc.. where m0,m1,m2... are the input parameter 'm'.
1318: The DIAGONAL portion of the local submatrix of a processor can be defined
1319: as the submatrix which is obtained by extraction the part corresponding to
1320: the rows r1-r2 and columns c1-c2 of the global matrix, where r1 is the
1321: first row that belongs to the processor, r2 is the last row belonging to
1322: the this processor, and c1-c2 is range of indices of the local part of a
1323: vector suitable for applying the matrix to. This is an mxn matrix. In the
1324: common case of a square matrix, the row and column ranges are the same and
1325: the DIAGONAL part is also square. The remaining portion of the local
1326: submatrix (mxN) constitute the OFF-DIAGONAL portion.
1328: If `o_nnz`, `d_nnz` are specified, then `o_nz`, and `d_nz` are ignored.
1330: You can call `MatGetInfo()` to get information on how effective the preallocation was;
1331: for example the fields mallocs,nz_allocated,nz_used,nz_unneeded;
1332: You can also run with the option -info and look for messages with the string
1333: malloc in them to see if additional memory allocation was needed.
1335: .seealso: `Mat`, `MatCreate()`, `MatCreateSeqSELL()`, `MatSetValues()`, `MatCreateSELL()`,
1336: `MATMPISELL`, `MatGetInfo()`, `PetscSplitOwnership()`, `MATSELL`
1337: @*/
1338: PetscErrorCode MatMPISELLSetPreallocation(Mat B, PetscInt d_nz, const PetscInt d_nnz[], PetscInt o_nz, const PetscInt o_nnz[])
1339: {
1340: PetscFunctionBegin;
1343: PetscTryMethod(B, "MatMPISELLSetPreallocation_C", (Mat, PetscInt, const PetscInt[], PetscInt, const PetscInt[]), (B, d_nz, d_nnz, o_nz, o_nnz));
1344: PetscFunctionReturn(PETSC_SUCCESS);
1345: }
1347: /*MC
1348: MATMPISELL - MATMPISELL = "mpisell" - A matrix type to be used for MPI sparse matrices,
1349: based on the sliced Ellpack format
1351: Options Database Key:
1352: . -mat_type sell - sets the matrix type to `MATSELL` during a call to `MatSetFromOptions()`
1354: Level: beginner
1356: .seealso: `Mat`, `MatCreateSELL()`, `MATSEQSELL`, `MATSELL`, `MATSEQAIJ`, `MATAIJ`, `MATMPIAIJ`
1357: M*/
1359: /*@C
1360: MatCreateSELL - Creates a sparse parallel matrix in `MATSELL` format.
1362: Collective
1364: Input Parameters:
1365: + comm - MPI communicator
1366: . m - number of local rows (or `PETSC_DECIDE` to have calculated if M is given)
1367: This value should be the same as the local size used in creating the
1368: y vector for the matrix-vector product y = Ax.
1369: . n - This value should be the same as the local size used in creating the
1370: x vector for the matrix-vector product y = Ax. (or `PETSC_DECIDE` to have
1371: calculated if `N` is given) For square matrices n is almost always `m`.
1372: . M - number of global rows (or `PETSC_DETERMINE` to have calculated if `m` is given)
1373: . N - number of global columns (or `PETSC_DETERMINE` to have calculated if `n` is given)
1374: . d_rlenmax - max number of nonzeros per row in DIAGONAL portion of local submatrix
1375: (same value is used for all local rows)
1376: . d_rlen - array containing the number of nonzeros in the various rows of the
1377: DIAGONAL portion of the local submatrix (possibly different for each row)
1378: or `NULL`, if d_rlenmax is used to specify the nonzero structure.
1379: The size of this array is equal to the number of local rows, i.e `m`.
1380: . o_rlenmax - max number of nonzeros per row in the OFF-DIAGONAL portion of local
1381: submatrix (same value is used for all local rows).
1382: - o_rlen - array containing the number of nonzeros in the various rows of the
1383: OFF-DIAGONAL portion of the local submatrix (possibly different for
1384: each row) or `NULL`, if `o_rlenmax` is used to specify the nonzero
1385: structure. The size of this array is equal to the number
1386: of local rows, i.e `m`.
1388: Output Parameter:
1389: . A - the matrix
1391: Options Database Key:
1392: . -mat_sell_oneindex - Internally use indexing starting at 1
1393: rather than 0. When calling `MatSetValues()`,
1394: the user still MUST index entries starting at 0!
1396: Example:
1397: Consider the following 8x8 matrix with 34 non-zero values, that is
1398: assembled across 3 processors. Lets assume that proc0 owns 3 rows,
1399: proc1 owns 3 rows, proc2 owns 2 rows. This division can be shown
1400: as follows
1402: .vb
1403: 1 2 0 | 0 3 0 | 0 4
1404: Proc0 0 5 6 | 7 0 0 | 8 0
1405: 9 0 10 | 11 0 0 | 12 0
1406: -------------------------------------
1407: 13 0 14 | 15 16 17 | 0 0
1408: Proc1 0 18 0 | 19 20 21 | 0 0
1409: 0 0 0 | 22 23 0 | 24 0
1410: -------------------------------------
1411: Proc2 25 26 27 | 0 0 28 | 29 0
1412: 30 0 0 | 31 32 33 | 0 34
1413: .ve
1415: This can be represented as a collection of submatrices as
1416: .vb
1417: A B C
1418: D E F
1419: G H I
1420: .ve
1422: Where the submatrices A,B,C are owned by proc0, D,E,F are
1423: owned by proc1, G,H,I are owned by proc2.
1425: The 'm' parameters for proc0,proc1,proc2 are 3,3,2 respectively.
1426: The 'n' parameters for proc0,proc1,proc2 are 3,3,2 respectively.
1427: The 'M','N' parameters are 8,8, and have the same values on all procs.
1429: The DIAGONAL submatrices corresponding to proc0,proc1,proc2 are
1430: submatrices [A], [E], [I] respectively. The OFF-DIAGONAL submatrices
1431: corresponding to proc0,proc1,proc2 are [BC], [DF], [GH] respectively.
1432: Internally, each processor stores the DIAGONAL part, and the OFF-DIAGONAL
1433: part as `MATSEQSELL` matrices. For example, proc1 will store [E] as a `MATSEQSELL`
1434: matrix, and [DF] as another `MATSEQSELL` matrix.
1436: When d_rlenmax, o_rlenmax parameters are specified, d_rlenmax storage elements are
1437: allocated for every row of the local DIAGONAL submatrix, and o_rlenmax
1438: storage locations are allocated for every row of the OFF-DIAGONAL submatrix.
1439: One way to choose `d_rlenmax` and `o_rlenmax` is to use the maximum number of nonzeros over
1440: the local rows for each of the local DIAGONAL, and the OFF-DIAGONAL submatrices.
1441: In this case, the values of d_rlenmax,o_rlenmax are
1442: .vb
1443: proc0 - d_rlenmax = 2, o_rlenmax = 2
1444: proc1 - d_rlenmax = 3, o_rlenmax = 2
1445: proc2 - d_rlenmax = 1, o_rlenmax = 4
1446: .ve
1447: We are allocating m*(d_rlenmax+o_rlenmax) storage locations for every proc. This
1448: translates to 3*(2+2)=12 for proc0, 3*(3+2)=15 for proc1, 2*(1+4)=10
1449: for proc3. i.e we are using 12+15+10=37 storage locations to store
1450: 34 values.
1452: When `d_rlen`, `o_rlen` parameters are specified, the storage is specified
1453: for every row, corresponding to both DIAGONAL and OFF-DIAGONAL submatrices.
1454: In the above case the values for `d_nnz`, `o_nnz` are
1455: .vb
1456: proc0 - d_nnz = [2,2,2] and o_nnz = [2,2,2]
1457: proc1 - d_nnz = [3,3,2] and o_nnz = [2,1,1]
1458: proc2 - d_nnz = [1,1] and o_nnz = [4,4]
1459: .ve
1460: Here the space allocated is still 37 though there are 34 nonzeros because
1461: the allocation is always done according to rlenmax.
1463: Level: intermediate
1465: Notes:
1466: It is recommended that one use the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()`,
1467: MatXXXXSetPreallocation() paradigm instead of this routine directly.
1468: [MatXXXXSetPreallocation() is, for example, `MatSeqSELLSetPreallocation()`]
1470: If the *_rlen parameter is given then the *_rlenmax parameter is ignored
1472: `m`, `n`, `M`, `N` parameters specify the size of the matrix, and its partitioning across
1473: processors, while `d_rlenmax`, `d_rlen`, `o_rlenmax` , `o_rlen` parameters specify the approximate
1474: storage requirements for this matrix.
1476: If `PETSC_DECIDE` or `PETSC_DETERMINE` is used for a particular argument on one
1477: processor than it must be used on all processors that share the object for
1478: that argument.
1480: The user MUST specify either the local or global matrix dimensions
1481: (possibly both).
1483: The parallel matrix is partitioned across processors such that the
1484: first m0 rows belong to process 0, the next m1 rows belong to
1485: process 1, the next m2 rows belong to process 2 etc.. where
1486: m0,m1,m2,.. are the input parameter 'm'. i.e each processor stores
1487: values corresponding to [`m` x `N`] submatrix.
1489: The columns are logically partitioned with the n0 columns belonging
1490: to 0th partition, the next n1 columns belonging to the next
1491: partition etc.. where n0,n1,n2... are the input parameter `n`.
1493: The DIAGONAL portion of the local submatrix on any given processor
1494: is the submatrix corresponding to the rows and columns `m`, `n`
1495: corresponding to the given processor. i.e diagonal matrix on
1496: process 0 is [m0 x n0], diagonal matrix on process 1 is [m1 x n1]
1497: etc. The remaining portion of the local submatrix [m x (N-n)]
1498: constitute the OFF-DIAGONAL portion. The example below better
1499: illustrates this concept.
1501: For a square global matrix we define each processor's diagonal portion
1502: to be its local rows and the corresponding columns (a square submatrix);
1503: each processor's off-diagonal portion encompasses the remainder of the
1504: local matrix (a rectangular submatrix).
1506: If `o_rlen`, `d_rlen` are specified, then `o_rlenmax`, and `d_rlenmax` are ignored.
1508: When calling this routine with a single process communicator, a matrix of
1509: type `MATSEQSELL` is returned. If a matrix of type `MATMPISELL` is desired for this
1510: type of communicator, use the construction mechanism
1511: .vb
1512: MatCreate(...,&A);
1513: MatSetType(A,MATMPISELL);
1514: MatSetSizes(A, m,n,M,N);
1515: MatMPISELLSetPreallocation(A,...);
1516: .ve
1518: .seealso: `Mat`, `MATSELL`, `MatCreate()`, `MatCreateSeqSELL()`, `MatSetValues()`, `MatMPISELLSetPreallocation()`, `MATMPISELL`
1519: @*/
1520: PetscErrorCode MatCreateSELL(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, PetscInt N, PetscInt d_rlenmax, const PetscInt d_rlen[], PetscInt o_rlenmax, const PetscInt o_rlen[], Mat *A)
1521: {
1522: PetscMPIInt size;
1524: PetscFunctionBegin;
1525: PetscCall(MatCreate(comm, A));
1526: PetscCall(MatSetSizes(*A, m, n, M, N));
1527: PetscCallMPI(MPI_Comm_size(comm, &size));
1528: if (size > 1) {
1529: PetscCall(MatSetType(*A, MATMPISELL));
1530: PetscCall(MatMPISELLSetPreallocation(*A, d_rlenmax, d_rlen, o_rlenmax, o_rlen));
1531: } else {
1532: PetscCall(MatSetType(*A, MATSEQSELL));
1533: PetscCall(MatSeqSELLSetPreallocation(*A, d_rlenmax, d_rlen));
1534: }
1535: PetscFunctionReturn(PETSC_SUCCESS);
1536: }
1538: /*@C
1539: MatMPISELLGetSeqSELL - Returns the local pieces of this distributed matrix
1541: Not Collective
1543: Input Parameter:
1544: . A - the `MATMPISELL` matrix
1546: Output Parameters:
1547: + Ad - The diagonal portion of `A`
1548: . Ao - The off-diagonal portion of `A`
1549: - colmap - An array mapping local column numbers of `Ao` to global column numbers of the parallel matrix
1551: Level: advanced
1553: .seealso: `Mat`, `MATSEQSELL`, `MATMPISELL`
1554: @*/
1555: PetscErrorCode MatMPISELLGetSeqSELL(Mat A, Mat *Ad, Mat *Ao, const PetscInt *colmap[])
1556: {
1557: Mat_MPISELL *a = (Mat_MPISELL *)A->data;
1558: PetscBool flg;
1560: PetscFunctionBegin;
1561: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATMPISELL, &flg));
1562: PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "This function requires a MATMPISELL matrix as input");
1563: if (Ad) *Ad = a->A;
1564: if (Ao) *Ao = a->B;
1565: if (colmap) *colmap = a->garray;
1566: PetscFunctionReturn(PETSC_SUCCESS);
1567: }
1569: /*@C
1570: MatMPISELLGetLocalMatCondensed - Creates a `MATSEQSELL` matrix from an `MATMPISELL` matrix by
1571: taking all its local rows and NON-ZERO columns
1573: Not Collective
1575: Input Parameters:
1576: + A - the matrix
1577: . scall - either `MAT_INITIAL_MATRIX` or `MAT_REUSE_MATRIX`
1578: . row - index sets of rows to extract (or `NULL`)
1579: - col - index sets of columns to extract (or `NULL`)
1581: Output Parameter:
1582: . A_loc - the local sequential matrix generated
1584: Level: advanced
1586: .seealso: `Mat`, `MATSEQSELL`, `MATMPISELL`, `MatGetOwnershipRange()`, `MatMPISELLGetLocalMat()`
1587: @*/
1588: PetscErrorCode MatMPISELLGetLocalMatCondensed(Mat A, MatReuse scall, IS *row, IS *col, Mat *A_loc)
1589: {
1590: Mat_MPISELL *a = (Mat_MPISELL *)A->data;
1591: PetscInt i, start, end, ncols, nzA, nzB, *cmap, imark, *idx;
1592: IS isrowa, iscola;
1593: Mat *aloc;
1594: PetscBool match;
1596: PetscFunctionBegin;
1597: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATMPISELL, &match));
1598: PetscCheck(match, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Requires MATMPISELL matrix as input");
1599: PetscCall(PetscLogEventBegin(MAT_Getlocalmatcondensed, A, 0, 0, 0));
1600: if (!row) {
1601: start = A->rmap->rstart;
1602: end = A->rmap->rend;
1603: PetscCall(ISCreateStride(PETSC_COMM_SELF, end - start, start, 1, &isrowa));
1604: } else {
1605: isrowa = *row;
1606: }
1607: if (!col) {
1608: start = A->cmap->rstart;
1609: cmap = a->garray;
1610: nzA = a->A->cmap->n;
1611: nzB = a->B->cmap->n;
1612: PetscCall(PetscMalloc1(nzA + nzB, &idx));
1613: ncols = 0;
1614: for (i = 0; i < nzB; i++) {
1615: if (cmap[i] < start) idx[ncols++] = cmap[i];
1616: else break;
1617: }
1618: imark = i;
1619: for (i = 0; i < nzA; i++) idx[ncols++] = start + i;
1620: for (i = imark; i < nzB; i++) idx[ncols++] = cmap[i];
1621: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, ncols, idx, PETSC_OWN_POINTER, &iscola));
1622: } else {
1623: iscola = *col;
1624: }
1625: if (scall != MAT_INITIAL_MATRIX) {
1626: PetscCall(PetscMalloc1(1, &aloc));
1627: aloc[0] = *A_loc;
1628: }
1629: PetscCall(MatCreateSubMatrices(A, 1, &isrowa, &iscola, scall, &aloc));
1630: *A_loc = aloc[0];
1631: PetscCall(PetscFree(aloc));
1632: if (!row) PetscCall(ISDestroy(&isrowa));
1633: if (!col) PetscCall(ISDestroy(&iscola));
1634: PetscCall(PetscLogEventEnd(MAT_Getlocalmatcondensed, A, 0, 0, 0));
1635: PetscFunctionReturn(PETSC_SUCCESS);
1636: }
1638: #include <../src/mat/impls/aij/mpi/mpiaij.h>
1640: PetscErrorCode MatConvert_MPISELL_MPIAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
1641: {
1642: Mat_MPISELL *a = (Mat_MPISELL *)A->data;
1643: Mat B;
1644: Mat_MPIAIJ *b;
1646: PetscFunctionBegin;
1647: PetscCheck(A->assembled, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Matrix must be assembled");
1649: if (reuse == MAT_REUSE_MATRIX) {
1650: B = *newmat;
1651: } else {
1652: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
1653: PetscCall(MatSetType(B, MATMPIAIJ));
1654: PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
1655: PetscCall(MatSetBlockSizes(B, A->rmap->bs, A->cmap->bs));
1656: PetscCall(MatSeqAIJSetPreallocation(B, 0, NULL));
1657: PetscCall(MatMPIAIJSetPreallocation(B, 0, NULL, 0, NULL));
1658: }
1659: b = (Mat_MPIAIJ *)B->data;
1661: if (reuse == MAT_REUSE_MATRIX) {
1662: PetscCall(MatConvert_SeqSELL_SeqAIJ(a->A, MATSEQAIJ, MAT_REUSE_MATRIX, &b->A));
1663: PetscCall(MatConvert_SeqSELL_SeqAIJ(a->B, MATSEQAIJ, MAT_REUSE_MATRIX, &b->B));
1664: } else {
1665: PetscCall(MatDestroy(&b->A));
1666: PetscCall(MatDestroy(&b->B));
1667: PetscCall(MatDisAssemble_MPISELL(A));
1668: PetscCall(MatConvert_SeqSELL_SeqAIJ(a->A, MATSEQAIJ, MAT_INITIAL_MATRIX, &b->A));
1669: PetscCall(MatConvert_SeqSELL_SeqAIJ(a->B, MATSEQAIJ, MAT_INITIAL_MATRIX, &b->B));
1670: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
1671: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
1672: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
1673: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
1674: }
1676: if (reuse == MAT_INPLACE_MATRIX) {
1677: PetscCall(MatHeaderReplace(A, &B));
1678: } else {
1679: *newmat = B;
1680: }
1681: PetscFunctionReturn(PETSC_SUCCESS);
1682: }
1684: PetscErrorCode MatConvert_MPIAIJ_MPISELL(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
1685: {
1686: Mat_MPIAIJ *a = (Mat_MPIAIJ *)A->data;
1687: Mat B;
1688: Mat_MPISELL *b;
1690: PetscFunctionBegin;
1691: PetscCheck(A->assembled, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Matrix must be assembled");
1693: if (reuse == MAT_REUSE_MATRIX) {
1694: B = *newmat;
1695: } else {
1696: Mat_SeqAIJ *Aa = (Mat_SeqAIJ *)a->A->data, *Ba = (Mat_SeqAIJ *)a->B->data;
1697: PetscInt i, d_nz = 0, o_nz = 0, m = A->rmap->N, n = A->cmap->N, lm = A->rmap->n, ln = A->cmap->n;
1698: PetscInt *d_nnz, *o_nnz;
1699: PetscCall(PetscMalloc2(lm, &d_nnz, lm, &o_nnz));
1700: for (i = 0; i < lm; i++) {
1701: d_nnz[i] = Aa->i[i + 1] - Aa->i[i];
1702: o_nnz[i] = Ba->i[i + 1] - Ba->i[i];
1703: if (d_nnz[i] > d_nz) d_nz = d_nnz[i];
1704: if (o_nnz[i] > o_nz) o_nz = o_nnz[i];
1705: }
1706: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
1707: PetscCall(MatSetType(B, MATMPISELL));
1708: PetscCall(MatSetSizes(B, lm, ln, m, n));
1709: PetscCall(MatSetBlockSizes(B, A->rmap->bs, A->cmap->bs));
1710: PetscCall(MatSeqSELLSetPreallocation(B, d_nz, d_nnz));
1711: PetscCall(MatMPISELLSetPreallocation(B, d_nz, d_nnz, o_nz, o_nnz));
1712: PetscCall(PetscFree2(d_nnz, o_nnz));
1713: }
1714: b = (Mat_MPISELL *)B->data;
1716: if (reuse == MAT_REUSE_MATRIX) {
1717: PetscCall(MatConvert_SeqAIJ_SeqSELL(a->A, MATSEQSELL, MAT_REUSE_MATRIX, &b->A));
1718: PetscCall(MatConvert_SeqAIJ_SeqSELL(a->B, MATSEQSELL, MAT_REUSE_MATRIX, &b->B));
1719: } else {
1720: PetscBool nooffprocentries_A = A->nooffprocentries, nooffprocentries_B = B->nooffprocentries;
1722: PetscCall(MatDestroy(&b->A));
1723: PetscCall(MatDestroy(&b->B));
1724: /* Expand a->B from compacted local off-diag columns back to global columns so the new MPISELL's
1725: MatAssemblyEnd() builds the correct garray/Mvctx for its off-diagonal block. */
1726: PetscCall(MatDisAssemble_MPIAIJ(A, PETSC_FALSE));
1727: PetscCall(MatConvert_SeqAIJ_SeqSELL(a->A, MATSEQSELL, MAT_INITIAL_MATRIX, &b->A));
1728: PetscCall(MatConvert_SeqAIJ_SeqSELL(a->B, MATSEQSELL, MAT_INITIAL_MATRIX, &b->B));
1729: /* The locally-populated A and B have no stashed off-processor entries, so skip the stash scatter. */
1730: A->nooffprocentries = PETSC_TRUE;
1731: B->nooffprocentries = PETSC_TRUE;
1732: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
1733: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
1734: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
1735: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
1736: A->nooffprocentries = nooffprocentries_A;
1737: B->nooffprocentries = nooffprocentries_B;
1738: }
1740: if (reuse == MAT_INPLACE_MATRIX) {
1741: PetscCall(MatHeaderReplace(A, &B));
1742: } else {
1743: *newmat = B;
1744: }
1745: PetscFunctionReturn(PETSC_SUCCESS);
1746: }
1748: PetscErrorCode MatSOR_MPISELL(Mat matin, Vec bb, PetscReal omega, MatSORType flag, PetscReal fshift, PetscInt its, PetscInt lits, Vec xx)
1749: {
1750: Mat_MPISELL *mat = (Mat_MPISELL *)matin->data;
1751: Vec bb1 = NULL;
1753: PetscFunctionBegin;
1754: if (flag == SOR_APPLY_UPPER) {
1755: PetscUseTypeMethod(mat->A, sor, bb, omega, flag, fshift, lits, 1, xx);
1756: PetscFunctionReturn(PETSC_SUCCESS);
1757: }
1759: if (its > 1 || ~flag & SOR_ZERO_INITIAL_GUESS || flag & SOR_EISENSTAT) PetscCall(VecDuplicate(bb, &bb1));
1761: if ((flag & SOR_LOCAL_SYMMETRIC_SWEEP) == SOR_LOCAL_SYMMETRIC_SWEEP) {
1762: if (flag & SOR_ZERO_INITIAL_GUESS) {
1763: PetscUseTypeMethod(mat->A, sor, bb, omega, flag, fshift, lits, 1, xx);
1764: its--;
1765: }
1767: while (its--) {
1768: PetscCall(VecScatterBegin(mat->Mvctx, xx, mat->lvec, INSERT_VALUES, SCATTER_FORWARD));
1769: PetscCall(VecScatterEnd(mat->Mvctx, xx, mat->lvec, INSERT_VALUES, SCATTER_FORWARD));
1771: /* update rhs: bb1 = bb - B*x */
1772: PetscCall(VecScale(mat->lvec, -1.0));
1773: PetscUseTypeMethod(mat->B, multadd, mat->lvec, bb, bb1);
1775: /* local sweep */
1776: PetscUseTypeMethod(mat->A, sor, bb1, omega, SOR_SYMMETRIC_SWEEP, fshift, lits, 1, xx);
1777: }
1778: } else if (flag & SOR_LOCAL_FORWARD_SWEEP) {
1779: if (flag & SOR_ZERO_INITIAL_GUESS) {
1780: PetscUseTypeMethod(mat->A, sor, bb, omega, flag, fshift, lits, 1, xx);
1781: its--;
1782: }
1783: while (its--) {
1784: PetscCall(VecScatterBegin(mat->Mvctx, xx, mat->lvec, INSERT_VALUES, SCATTER_FORWARD));
1785: PetscCall(VecScatterEnd(mat->Mvctx, xx, mat->lvec, INSERT_VALUES, SCATTER_FORWARD));
1787: /* update rhs: bb1 = bb - B*x */
1788: PetscCall(VecScale(mat->lvec, -1.0));
1789: PetscUseTypeMethod(mat->B, multadd, mat->lvec, bb, bb1);
1791: /* local sweep */
1792: PetscUseTypeMethod(mat->A, sor, bb1, omega, SOR_FORWARD_SWEEP, fshift, lits, 1, xx);
1793: }
1794: } else if (flag & SOR_LOCAL_BACKWARD_SWEEP) {
1795: if (flag & SOR_ZERO_INITIAL_GUESS) {
1796: PetscUseTypeMethod(mat->A, sor, bb, omega, flag, fshift, lits, 1, xx);
1797: its--;
1798: }
1799: while (its--) {
1800: PetscCall(VecScatterBegin(mat->Mvctx, xx, mat->lvec, INSERT_VALUES, SCATTER_FORWARD));
1801: PetscCall(VecScatterEnd(mat->Mvctx, xx, mat->lvec, INSERT_VALUES, SCATTER_FORWARD));
1803: /* update rhs: bb1 = bb - B*x */
1804: PetscCall(VecScale(mat->lvec, -1.0));
1805: PetscUseTypeMethod(mat->B, multadd, mat->lvec, bb, bb1);
1807: /* local sweep */
1808: PetscUseTypeMethod(mat->A, sor, bb1, omega, SOR_BACKWARD_SWEEP, fshift, lits, 1, xx);
1809: }
1810: } else SETERRQ(PetscObjectComm((PetscObject)matin), PETSC_ERR_SUP, "Parallel SOR not supported");
1812: PetscCall(VecDestroy(&bb1));
1814: matin->factorerrortype = mat->A->factorerrortype;
1815: PetscFunctionReturn(PETSC_SUCCESS);
1816: }
1818: #if PetscDefined(HAVE_CUDA)
1819: PETSC_INTERN PetscErrorCode MatConvert_MPISELL_MPISELLCUDA(Mat, MatType, MatReuse, Mat *);
1820: #endif
1822: /*MC
1823: MATMPISELL - MATMPISELL = "MPISELL" - A matrix type to be used for parallel sparse matrices.
1825: Options Database Keys:
1826: . -mat_type mpisell - sets the matrix type to `MATMPISELL` during a call to `MatSetFromOptions()`
1828: Level: beginner
1830: .seealso: `Mat`, `MATSELL`, `MATSEQSELL`, `MatCreateSELL()`
1831: M*/
1832: PETSC_EXTERN PetscErrorCode MatCreate_MPISELL(Mat B)
1833: {
1834: Mat_MPISELL *b;
1835: PetscMPIInt size;
1837: PetscFunctionBegin;
1838: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)B), &size));
1839: PetscCall(PetscNew(&b));
1840: B->data = (void *)b;
1841: B->ops[0] = MatOps_Values;
1842: B->assembled = PETSC_FALSE;
1843: B->insertmode = NOT_SET_VALUES;
1844: b->size = size;
1845: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)B), &b->rank));
1846: /* build cache for off array entries formed */
1847: PetscCall(MatStashCreate_Private(PetscObjectComm((PetscObject)B), 1, &B->stash));
1849: b->donotstash = PETSC_FALSE;
1850: b->colmap = NULL;
1851: b->garray = NULL;
1852: b->roworiented = PETSC_TRUE;
1854: /* stuff used for matrix vector multiply */
1855: b->lvec = NULL;
1856: b->Mvctx = NULL;
1858: /* stuff for MatGetRow() */
1859: b->rowindices = NULL;
1860: b->rowvalues = NULL;
1861: b->getrowactive = PETSC_FALSE;
1863: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatStoreValues_C", MatStoreValues_MPISELL));
1864: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatRetrieveValues_C", MatRetrieveValues_MPISELL));
1865: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatIsTranspose_C", MatIsTranspose_MPISELL));
1866: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatMPISELLSetPreallocation_C", MatMPISELLSetPreallocation_MPISELL));
1867: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_mpisell_mpiaij_C", MatConvert_MPISELL_MPIAIJ));
1868: #if PetscDefined(HAVE_CUDA)
1869: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_mpisell_mpisellcuda_C", MatConvert_MPISELL_MPISELLCUDA));
1870: #endif
1871: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDiagonalScaleLocal_C", MatDiagonalScaleLocal_MPISELL));
1872: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatGetMultPetscSF_C", MatGetMultPetscSF_MPISELL));
1873: PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATMPISELL));
1874: PetscFunctionReturn(PETSC_SUCCESS);
1875: }