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