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;
843: PetscFunctionBegin;
844: PetscCall(MatSetRandom(sell->A, rctx));
845: PetscCall(MatSetRandom(sell->B, rctx));
846: PetscCall(MatAssemblyBegin(x, MAT_FINAL_ASSEMBLY));
847: PetscCall(MatAssemblyEnd(x, MAT_FINAL_ASSEMBLY));
848: PetscFunctionReturn(PETSC_SUCCESS);
849: }
851: static PetscErrorCode MatSetFromOptions_MPISELL(Mat A, PetscOptionItems PetscOptionsObject)
852: {
853: PetscFunctionBegin;
854: PetscOptionsHeadBegin(PetscOptionsObject, "MPISELL options");
855: PetscOptionsHeadEnd();
856: PetscFunctionReturn(PETSC_SUCCESS);
857: }
859: static PetscErrorCode MatShift_MPISELL(Mat Y, PetscScalar a)
860: {
861: Mat_MPISELL *msell = (Mat_MPISELL *)Y->data;
862: Mat_SeqSELL *sell = (Mat_SeqSELL *)msell->A->data;
864: PetscFunctionBegin;
865: if (!Y->preallocated) {
866: PetscCall(MatMPISELLSetPreallocation(Y, 1, NULL, 0, NULL));
867: } else if (!sell->nz) {
868: PetscInt nonew = sell->nonew;
869: PetscCall(MatSeqSELLSetPreallocation(msell->A, 1, NULL));
870: sell->nonew = nonew;
871: }
872: PetscCall(MatShift_Basic(Y, a));
873: PetscFunctionReturn(PETSC_SUCCESS);
874: }
876: static PetscErrorCode MatGetDiagonalBlock_MPISELL(Mat A, Mat *a)
877: {
878: PetscFunctionBegin;
879: *a = ((Mat_MPISELL *)A->data)->A;
880: PetscFunctionReturn(PETSC_SUCCESS);
881: }
883: static PetscErrorCode MatStoreValues_MPISELL(Mat mat)
884: {
885: Mat_MPISELL *sell = (Mat_MPISELL *)mat->data;
887: PetscFunctionBegin;
888: PetscCall(MatStoreValues(sell->A));
889: PetscCall(MatStoreValues(sell->B));
890: PetscFunctionReturn(PETSC_SUCCESS);
891: }
893: static PetscErrorCode MatRetrieveValues_MPISELL(Mat mat)
894: {
895: Mat_MPISELL *sell = (Mat_MPISELL *)mat->data;
897: PetscFunctionBegin;
898: PetscCall(MatRetrieveValues(sell->A));
899: PetscCall(MatRetrieveValues(sell->B));
900: PetscFunctionReturn(PETSC_SUCCESS);
901: }
903: static PetscErrorCode MatMPISELLSetPreallocation_MPISELL(Mat B, PetscInt d_rlenmax, const PetscInt d_rlen[], PetscInt o_rlenmax, const PetscInt o_rlen[])
904: {
905: Mat_MPISELL *b;
907: PetscFunctionBegin;
908: PetscCall(PetscLayoutSetUp(B->rmap));
909: PetscCall(PetscLayoutSetUp(B->cmap));
910: b = (Mat_MPISELL *)B->data;
912: if (!B->preallocated) {
913: /* Explicitly create 2 MATSEQSELL matrices. */
914: PetscCall(MatCreate(PETSC_COMM_SELF, &b->A));
915: PetscCall(MatSetSizes(b->A, B->rmap->n, B->cmap->n, B->rmap->n, B->cmap->n));
916: PetscCall(MatSetBlockSizesFromMats(b->A, B, B));
917: PetscCall(MatSetType(b->A, MATSEQSELL));
918: PetscCall(MatCreate(PETSC_COMM_SELF, &b->B));
919: PetscCall(MatSetSizes(b->B, B->rmap->n, B->cmap->N, B->rmap->n, B->cmap->N));
920: PetscCall(MatSetBlockSizesFromMats(b->B, B, B));
921: PetscCall(MatSetType(b->B, MATSEQSELL));
922: }
924: PetscCall(MatSeqSELLSetPreallocation(b->A, d_rlenmax, d_rlen));
925: PetscCall(MatSeqSELLSetPreallocation(b->B, o_rlenmax, o_rlen));
926: B->preallocated = PETSC_TRUE;
927: B->was_assembled = PETSC_FALSE;
928: /*
929: critical for MatAssemblyEnd to work.
930: MatAssemblyBegin checks it to set up was_assembled
931: and MatAssemblyEnd checks was_assembled to determine whether to build garray
932: */
933: B->assembled = PETSC_FALSE;
934: PetscFunctionReturn(PETSC_SUCCESS);
935: }
937: static PetscErrorCode MatDuplicate_MPISELL(Mat matin, MatDuplicateOption cpvalues, Mat *newmat)
938: {
939: Mat mat;
940: Mat_MPISELL *a, *oldmat = (Mat_MPISELL *)matin->data;
942: PetscFunctionBegin;
943: *newmat = NULL;
944: PetscCall(MatCreate(PetscObjectComm((PetscObject)matin), &mat));
945: PetscCall(MatSetSizes(mat, matin->rmap->n, matin->cmap->n, matin->rmap->N, matin->cmap->N));
946: PetscCall(MatSetBlockSizesFromMats(mat, matin, matin));
947: PetscCall(MatSetType(mat, ((PetscObject)matin)->type_name));
948: a = (Mat_MPISELL *)mat->data;
950: mat->factortype = matin->factortype;
951: mat->assembled = PETSC_TRUE;
952: mat->insertmode = NOT_SET_VALUES;
953: mat->preallocated = PETSC_TRUE;
955: a->size = oldmat->size;
956: a->rank = oldmat->rank;
957: a->donotstash = oldmat->donotstash;
958: a->roworiented = oldmat->roworiented;
959: a->rowindices = NULL;
960: a->rowvalues = NULL;
961: a->getrowactive = PETSC_FALSE;
963: PetscCall(PetscLayoutReference(matin->rmap, &mat->rmap));
964: PetscCall(PetscLayoutReference(matin->cmap, &mat->cmap));
966: if (oldmat->colmap) {
967: #if PetscDefined(USE_CTABLE)
968: PetscCall(PetscHMapIDuplicate(oldmat->colmap, &a->colmap));
969: #else
970: PetscCall(PetscMalloc1(mat->cmap->N, &a->colmap));
971: PetscCall(PetscArraycpy(a->colmap, oldmat->colmap, mat->cmap->N));
972: #endif
973: } else a->colmap = NULL;
974: if (oldmat->garray) {
975: PetscInt len;
976: len = oldmat->B->cmap->n;
977: PetscCall(PetscMalloc1(len + 1, &a->garray));
978: if (len) PetscCall(PetscArraycpy(a->garray, oldmat->garray, len));
979: } else a->garray = NULL;
981: PetscCall(VecDuplicate(oldmat->lvec, &a->lvec));
982: PetscCall(VecScatterCopy(oldmat->Mvctx, &a->Mvctx));
983: PetscCall(MatDuplicate(oldmat->A, cpvalues, &a->A));
984: PetscCall(MatDuplicate(oldmat->B, cpvalues, &a->B));
985: PetscCall(PetscFunctionListDuplicate(((PetscObject)matin)->qlist, &((PetscObject)mat)->qlist));
986: *newmat = mat;
987: PetscFunctionReturn(PETSC_SUCCESS);
988: }
990: static const struct _MatOps MatOps_Values = {MatSetValues_MPISELL,
991: NULL,
992: NULL,
993: MatMult_MPISELL,
994: /* 4*/ MatMultAdd_MPISELL,
995: MatMultTranspose_MPISELL,
996: MatMultTransposeAdd_MPISELL,
997: NULL,
998: NULL,
999: NULL,
1000: /*10*/ NULL,
1001: NULL,
1002: NULL,
1003: MatSOR_MPISELL,
1004: NULL,
1005: /*15*/ MatGetInfo_MPISELL,
1006: MatEqual_MPISELL,
1007: MatGetDiagonal_MPISELL,
1008: MatDiagonalScale_MPISELL,
1009: NULL,
1010: /*20*/ MatAssemblyBegin_MPISELL,
1011: MatAssemblyEnd_MPISELL,
1012: MatSetOption_MPISELL,
1013: MatZeroEntries_MPISELL,
1014: /*24*/ NULL,
1015: NULL,
1016: NULL,
1017: NULL,
1018: NULL,
1019: /*29*/ MatSetUp_MPISELL,
1020: NULL,
1021: NULL,
1022: MatGetDiagonalBlock_MPISELL,
1023: NULL,
1024: /*34*/ MatDuplicate_MPISELL,
1025: NULL,
1026: NULL,
1027: NULL,
1028: NULL,
1029: /*39*/ NULL,
1030: NULL,
1031: NULL,
1032: MatGetValues_MPISELL,
1033: MatCopy_MPISELL,
1034: /*44*/ NULL,
1035: MatScale_MPISELL,
1036: MatShift_MPISELL,
1037: MatDiagonalSet_MPISELL,
1038: NULL,
1039: /*49*/ MatSetRandom_MPISELL,
1040: NULL,
1041: NULL,
1042: NULL,
1043: NULL,
1044: /*54*/ MatFDColoringCreate_MPIXAIJ,
1045: NULL,
1046: MatSetUnfactored_MPISELL,
1047: NULL,
1048: NULL,
1049: /*59*/ NULL,
1050: MatDestroy_MPISELL,
1051: MatView_MPISELL,
1052: NULL,
1053: NULL,
1054: /*64*/ NULL,
1055: NULL,
1056: NULL,
1057: NULL,
1058: NULL,
1059: /*69*/ NULL,
1060: NULL,
1061: NULL,
1062: MatFDColoringApply_AIJ, /* reuse AIJ function */
1063: MatSetFromOptions_MPISELL,
1064: NULL,
1065: /*75*/ NULL,
1066: NULL,
1067: NULL,
1068: NULL,
1069: NULL,
1070: /*80*/ NULL,
1071: NULL,
1072: NULL,
1073: /*83*/ NULL,
1074: NULL,
1075: NULL,
1076: NULL,
1077: NULL,
1078: NULL,
1079: /*89*/ NULL,
1080: NULL,
1081: NULL,
1082: NULL,
1083: MatConjugate_MPISELL,
1084: /*94*/ NULL,
1085: NULL,
1086: NULL,
1087: NULL,
1088: NULL,
1089: /*99*/ NULL,
1090: NULL,
1091: NULL,
1092: NULL,
1093: NULL,
1094: /*104*/ NULL,
1095: NULL,
1096: MatGetGhosts_MPISELL,
1097: NULL,
1098: NULL,
1099: /*109*/ MatMultDiagonalBlock_MPISELL,
1100: NULL,
1101: NULL,
1102: NULL,
1103: NULL,
1104: /*114*/ NULL,
1105: NULL,
1106: MatInvertBlockDiagonal_MPISELL,
1107: NULL,
1108: /*119*/ NULL,
1109: NULL,
1110: NULL,
1111: NULL,
1112: NULL,
1113: /*124*/ NULL,
1114: NULL,
1115: NULL,
1116: NULL,
1117: MatFDColoringSetUp_MPIXAIJ,
1118: /*129*/ NULL,
1119: NULL,
1120: NULL,
1121: NULL,
1122: NULL,
1123: /*134*/ NULL,
1124: NULL,
1125: NULL,
1126: NULL,
1127: NULL,
1128: /*139*/ NULL,
1129: NULL,
1130: NULL,
1131: NULL,
1132: NULL,
1133: NULL,
1134: /*144*/ NULL,
1135: NULL,
1136: NULL,
1137: NULL};
1139: /*@
1140: MatMPISELLSetPreallocation - Preallocates memory for a `MATMPISELL` sparse parallel matrix in sell format.
1141: For good matrix assembly performance the user should preallocate the matrix storage by
1142: setting the parameters `d_nz` (or `d_nnz`) and `o_nz` (or `o_nnz`).
1144: Collective
1146: Input Parameters:
1147: + B - the matrix
1148: . d_nz - number of nonzeros per row in DIAGONAL portion of local submatrix
1149: (same value is used for all local rows)
1150: . d_nnz - array containing the number of nonzeros in the various rows of the
1151: DIAGONAL portion of the local submatrix (possibly different for each row)
1152: or NULL (`PETSC_NULL_INTEGER` in Fortran), if `d_nz` is used to specify the nonzero structure.
1153: The size of this array is equal to the number of local rows, i.e 'm'.
1154: For matrices that will be factored, you must leave room for (and set)
1155: the diagonal entry even if it is zero.
1156: . o_nz - number of nonzeros per row in the OFF-DIAGONAL portion of local
1157: submatrix (same value is used for all local rows).
1158: - o_nnz - array containing the number of nonzeros in the various rows of the
1159: OFF-DIAGONAL portion of the local submatrix (possibly different for
1160: each row) or NULL (`PETSC_NULL_INTEGER` in Fortran), if `o_nz` is used to specify the nonzero
1161: structure. The size of this array is equal to the number
1162: of local rows, i.e 'm'.
1164: Example usage:
1165: Consider the following 8x8 matrix with 34 non-zero values, that is
1166: assembled across 3 processors. Lets assume that proc0 owns 3 rows,
1167: proc1 owns 3 rows, proc2 owns 2 rows. This division can be shown
1168: as follows
1170: .vb
1171: 1 2 0 | 0 3 0 | 0 4
1172: Proc0 0 5 6 | 7 0 0 | 8 0
1173: 9 0 10 | 11 0 0 | 12 0
1174: -------------------------------------
1175: 13 0 14 | 15 16 17 | 0 0
1176: Proc1 0 18 0 | 19 20 21 | 0 0
1177: 0 0 0 | 22 23 0 | 24 0
1178: -------------------------------------
1179: Proc2 25 26 27 | 0 0 28 | 29 0
1180: 30 0 0 | 31 32 33 | 0 34
1181: .ve
1183: This can be represented as a collection of submatrices as
1185: .vb
1186: A B C
1187: D E F
1188: G H I
1189: .ve
1191: Where the submatrices A,B,C are owned by proc0, D,E,F are
1192: owned by proc1, G,H,I are owned by proc2.
1194: The 'm' parameters for proc0,proc1,proc2 are 3,3,2 respectively.
1195: The 'n' parameters for proc0,proc1,proc2 are 3,3,2 respectively.
1196: The 'M','N' parameters are 8,8, and have the same values on all procs.
1198: The DIAGONAL submatrices corresponding to proc0,proc1,proc2 are
1199: submatrices [A], [E], [I] respectively. The OFF-DIAGONAL submatrices
1200: corresponding to proc0,proc1,proc2 are [BC], [DF], [GH] respectively.
1201: Internally, each processor stores the DIAGONAL part, and the OFF-DIAGONAL
1202: part as `MATSEQSELL` matrices. For example, proc1 will store [E] as a `MATSEQSELL`
1203: matrix, and [DF] as another SeqSELL matrix.
1205: When `d_nz`, `o_nz` parameters are specified, `d_nz` storage elements are
1206: allocated for every row of the local DIAGONAL submatrix, and o_nz
1207: storage locations are allocated for every row of the OFF-DIAGONAL submatrix.
1208: One way to choose `d_nz` and `o_nz` is to use the maximum number of nonzeros over
1209: the local rows for each of the local DIAGONAL, and the OFF-DIAGONAL submatrices.
1210: In this case, the values of d_nz,o_nz are
1211: .vb
1212: proc0 dnz = 2, o_nz = 2
1213: proc1 dnz = 3, o_nz = 2
1214: proc2 dnz = 1, o_nz = 4
1215: .ve
1216: We are allocating m*(d_nz+o_nz) storage locations for every proc. This
1217: translates to 3*(2+2)=12 for proc0, 3*(3+2)=15 for proc1, 2*(1+4)=10
1218: for proc3. i.e we are using 12+15+10=37 storage locations to store
1219: 34 values.
1221: When `d_nnz`, `o_nnz` parameters are specified, the storage is specified
1222: for every row, corresponding to both DIAGONAL and OFF-DIAGONAL submatrices.
1223: In the above case the values for d_nnz,o_nnz are
1224: .vb
1225: proc0 d_nnz = [2,2,2] and o_nnz = [2,2,2]
1226: proc1 d_nnz = [3,3,2] and o_nnz = [2,1,1]
1227: proc2 d_nnz = [1,1] and o_nnz = [4,4]
1228: .ve
1229: Here the space allocated is according to nz (or maximum values in the nnz
1230: if nnz is provided) for DIAGONAL and OFF-DIAGONAL submatrices, i.e (2+2+3+2)*3+(1+4)*2=37
1232: Level: intermediate
1234: Notes:
1235: If the *_nnz parameter is given then the *_nz parameter is ignored
1237: The stored row and column indices begin with zero.
1239: The parallel matrix is partitioned such that the first m0 rows belong to
1240: process 0, the next m1 rows belong to process 1, the next m2 rows belong
1241: to process 2 etc.. where m0,m1,m2... are the input parameter 'm'.
1243: The DIAGONAL portion of the local submatrix of a processor can be defined
1244: as the submatrix which is obtained by extraction the part corresponding to
1245: the rows r1-r2 and columns c1-c2 of the global matrix, where r1 is the
1246: first row that belongs to the processor, r2 is the last row belonging to
1247: the this processor, and c1-c2 is range of indices of the local part of a
1248: vector suitable for applying the matrix to. This is an mxn matrix. In the
1249: common case of a square matrix, the row and column ranges are the same and
1250: the DIAGONAL part is also square. The remaining portion of the local
1251: submatrix (mxN) constitute the OFF-DIAGONAL portion.
1253: If `o_nnz`, `d_nnz` are specified, then `o_nz`, and `d_nz` are ignored.
1255: You can call `MatGetInfo()` to get information on how effective the preallocation was;
1256: for example the fields mallocs,nz_allocated,nz_used,nz_unneeded;
1257: You can also run with the option -info and look for messages with the string
1258: malloc in them to see if additional memory allocation was needed.
1260: .seealso: `Mat`, `MatCreate()`, `MatCreateSeqSELL()`, `MatSetValues()`, `MatCreateSELL()`,
1261: `MATMPISELL`, `MatGetInfo()`, `PetscSplitOwnership()`, `MATSELL`
1262: @*/
1263: PetscErrorCode MatMPISELLSetPreallocation(Mat B, PetscInt d_nz, const PetscInt d_nnz[], PetscInt o_nz, const PetscInt o_nnz[])
1264: {
1265: PetscFunctionBegin;
1268: PetscTryMethod(B, "MatMPISELLSetPreallocation_C", (Mat, PetscInt, const PetscInt[], PetscInt, const PetscInt[]), (B, d_nz, d_nnz, o_nz, o_nnz));
1269: PetscFunctionReturn(PETSC_SUCCESS);
1270: }
1272: /*MC
1273: MATMPISELL - MATMPISELL = "mpisell" - A matrix type to be used for MPI sparse matrices,
1274: based on the sliced Ellpack format
1276: Options Database Key:
1277: . -mat_type sell - sets the matrix type to `MATSELL` during a call to `MatSetFromOptions()`
1279: Level: beginner
1281: .seealso: `Mat`, `MatCreateSELL()`, `MATSEQSELL`, `MATSELL`, `MATSEQAIJ`, `MATAIJ`, `MATMPIAIJ`
1282: M*/
1284: /*@
1285: MatCreateSELL - Creates a sparse parallel matrix in `MATSELL` format.
1287: Collective
1289: Input Parameters:
1290: + comm - MPI communicator
1291: . m - number of local rows (or `PETSC_DECIDE` to have calculated if M is given)
1292: This value should be the same as the local size used in creating the
1293: y vector for the matrix-vector product y = Ax.
1294: . n - This value should be the same as the local size used in creating the
1295: x vector for the matrix-vector product y = Ax. (or `PETSC_DECIDE` to have
1296: calculated if `N` is given) For square matrices n is almost always `m`.
1297: . M - number of global rows (or `PETSC_DETERMINE` to have calculated if `m` is given)
1298: . N - number of global columns (or `PETSC_DETERMINE` to have calculated if `n` is given)
1299: . d_rlenmax - max number of nonzeros per row in DIAGONAL portion of local submatrix
1300: (same value is used for all local rows)
1301: . d_rlen - array containing the number of nonzeros in the various rows of the
1302: DIAGONAL portion of the local submatrix (possibly different for each row)
1303: or `NULL`, if d_rlenmax is used to specify the nonzero structure.
1304: The size of this array is equal to the number of local rows, i.e `m`.
1305: . o_rlenmax - max number of nonzeros per row in the OFF-DIAGONAL portion of local
1306: submatrix (same value is used for all local rows).
1307: - o_rlen - array containing the number of nonzeros in the various rows of the
1308: OFF-DIAGONAL portion of the local submatrix (possibly different for
1309: each row) or `NULL`, if `o_rlenmax` is used to specify the nonzero
1310: structure. The size of this array is equal to the number
1311: of local rows, i.e `m`.
1313: Output Parameter:
1314: . A - the matrix
1316: Options Database Key:
1317: . -mat_sell_oneindex - Internally use indexing starting at 1
1318: rather than 0. When calling `MatSetValues()`,
1319: the user still MUST index entries starting at 0!
1321: Example:
1322: Consider the following 8x8 matrix with 34 non-zero values, that is
1323: assembled across 3 processors. Lets assume that proc0 owns 3 rows,
1324: proc1 owns 3 rows, proc2 owns 2 rows. This division can be shown
1325: as follows
1327: .vb
1328: 1 2 0 | 0 3 0 | 0 4
1329: Proc0 0 5 6 | 7 0 0 | 8 0
1330: 9 0 10 | 11 0 0 | 12 0
1331: -------------------------------------
1332: 13 0 14 | 15 16 17 | 0 0
1333: Proc1 0 18 0 | 19 20 21 | 0 0
1334: 0 0 0 | 22 23 0 | 24 0
1335: -------------------------------------
1336: Proc2 25 26 27 | 0 0 28 | 29 0
1337: 30 0 0 | 31 32 33 | 0 34
1338: .ve
1340: This can be represented as a collection of submatrices as
1341: .vb
1342: A B C
1343: D E F
1344: G H I
1345: .ve
1347: Where the submatrices A,B,C are owned by proc0, D,E,F are
1348: owned by proc1, G,H,I are owned by proc2.
1350: The 'm' parameters for proc0,proc1,proc2 are 3,3,2 respectively.
1351: The 'n' parameters for proc0,proc1,proc2 are 3,3,2 respectively.
1352: The 'M','N' parameters are 8,8, and have the same values on all procs.
1354: The DIAGONAL submatrices corresponding to proc0,proc1,proc2 are
1355: submatrices [A], [E], [I] respectively. The OFF-DIAGONAL submatrices
1356: corresponding to proc0,proc1,proc2 are [BC], [DF], [GH] respectively.
1357: Internally, each processor stores the DIAGONAL part, and the OFF-DIAGONAL
1358: part as `MATSEQSELL` matrices. For example, proc1 will store [E] as a `MATSEQSELL`
1359: matrix, and [DF] as another `MATSEQSELL` matrix.
1361: When d_rlenmax, o_rlenmax parameters are specified, d_rlenmax storage elements are
1362: allocated for every row of the local DIAGONAL submatrix, and o_rlenmax
1363: storage locations are allocated for every row of the OFF-DIAGONAL submatrix.
1364: One way to choose `d_rlenmax` and `o_rlenmax` is to use the maximum number of nonzeros over
1365: the local rows for each of the local DIAGONAL, and the OFF-DIAGONAL submatrices.
1366: In this case, the values of d_rlenmax,o_rlenmax are
1367: .vb
1368: proc0 - d_rlenmax = 2, o_rlenmax = 2
1369: proc1 - d_rlenmax = 3, o_rlenmax = 2
1370: proc2 - d_rlenmax = 1, o_rlenmax = 4
1371: .ve
1372: We are allocating m*(d_rlenmax+o_rlenmax) storage locations for every proc. This
1373: translates to 3*(2+2)=12 for proc0, 3*(3+2)=15 for proc1, 2*(1+4)=10
1374: for proc3. i.e we are using 12+15+10=37 storage locations to store
1375: 34 values.
1377: When `d_rlen`, `o_rlen` parameters are specified, the storage is specified
1378: for every row, corresponding to both DIAGONAL and OFF-DIAGONAL submatrices.
1379: In the above case the values for `d_nnz`, `o_nnz` are
1380: .vb
1381: proc0 - d_nnz = [2,2,2] and o_nnz = [2,2,2]
1382: proc1 - d_nnz = [3,3,2] and o_nnz = [2,1,1]
1383: proc2 - d_nnz = [1,1] and o_nnz = [4,4]
1384: .ve
1385: Here the space allocated is still 37 though there are 34 nonzeros because
1386: the allocation is always done according to rlenmax.
1388: Level: intermediate
1390: Notes:
1391: It is recommended that one use the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()`,
1392: MatXXXXSetPreallocation() paradigm instead of this routine directly.
1393: [MatXXXXSetPreallocation() is, for example, `MatSeqSELLSetPreallocation()`]
1395: If the *_rlen parameter is given then the *_rlenmax parameter is ignored
1397: `m`, `n`, `M`, `N` parameters specify the size of the matrix, and its partitioning across
1398: processors, while `d_rlenmax`, `d_rlen`, `o_rlenmax` , `o_rlen` parameters specify the approximate
1399: storage requirements for this matrix.
1401: If `PETSC_DECIDE` or `PETSC_DETERMINE` is used for a particular argument on one
1402: processor than it must be used on all processors that share the object for
1403: that argument.
1405: The user MUST specify either the local or global matrix dimensions
1406: (possibly both).
1408: The parallel matrix is partitioned across processors such that the
1409: first m0 rows belong to process 0, the next m1 rows belong to
1410: process 1, the next m2 rows belong to process 2 etc.. where
1411: m0,m1,m2,.. are the input parameter 'm'. i.e each processor stores
1412: values corresponding to [`m` x `N`] submatrix.
1414: The columns are logically partitioned with the n0 columns belonging
1415: to 0th partition, the next n1 columns belonging to the next
1416: partition etc.. where n0,n1,n2... are the input parameter `n`.
1418: The DIAGONAL portion of the local submatrix on any given processor
1419: is the submatrix corresponding to the rows and columns `m`, `n`
1420: corresponding to the given processor. i.e diagonal matrix on
1421: process 0 is [m0 x n0], diagonal matrix on process 1 is [m1 x n1]
1422: etc. The remaining portion of the local submatrix [m x (N-n)]
1423: constitute the OFF-DIAGONAL portion. The example below better
1424: illustrates this concept.
1426: For a square global matrix we define each processor's diagonal portion
1427: to be its local rows and the corresponding columns (a square submatrix);
1428: each processor's off-diagonal portion encompasses the remainder of the
1429: local matrix (a rectangular submatrix).
1431: If `o_rlen`, `d_rlen` are specified, then `o_rlenmax`, and `d_rlenmax` are ignored.
1433: When calling this routine with a single process communicator, a matrix of
1434: type `MATSEQSELL` is returned. If a matrix of type `MATMPISELL` is desired for this
1435: type of communicator, use the construction mechanism
1436: .vb
1437: MatCreate(...,&A);
1438: MatSetType(A,MATMPISELL);
1439: MatSetSizes(A, m,n,M,N);
1440: MatMPISELLSetPreallocation(A,...);
1441: .ve
1443: .seealso: `Mat`, `MATSELL`, `MatCreate()`, `MatCreateSeqSELL()`, `MatSetValues()`, `MatMPISELLSetPreallocation()`, `MATMPISELL`
1444: @*/
1445: 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)
1446: {
1447: PetscMPIInt size;
1449: PetscFunctionBegin;
1450: PetscCall(MatCreate(comm, A));
1451: PetscCall(MatSetSizes(*A, m, n, M, N));
1452: PetscCallMPI(MPI_Comm_size(comm, &size));
1453: if (size > 1) {
1454: PetscCall(MatSetType(*A, MATMPISELL));
1455: PetscCall(MatMPISELLSetPreallocation(*A, d_rlenmax, d_rlen, o_rlenmax, o_rlen));
1456: } else {
1457: PetscCall(MatSetType(*A, MATSEQSELL));
1458: PetscCall(MatSeqSELLSetPreallocation(*A, d_rlenmax, d_rlen));
1459: }
1460: PetscFunctionReturn(PETSC_SUCCESS);
1461: }
1463: /*@
1464: MatMPISELLGetSeqSELL - Returns the local pieces of this distributed matrix
1466: Not Collective
1468: Input Parameter:
1469: . A - the `MATMPISELL` matrix
1471: Output Parameters:
1472: + Ad - The diagonal portion of `A`
1473: . Ao - The off-diagonal portion of `A`
1474: - colmap - An array mapping local column numbers of `Ao` to global column numbers of the parallel matrix
1476: Level: advanced
1478: .seealso: `Mat`, `MATSEQSELL`, `MATMPISELL`
1479: @*/
1480: PetscErrorCode MatMPISELLGetSeqSELL(Mat A, Mat *Ad, Mat *Ao, const PetscInt *colmap[])
1481: {
1482: Mat_MPISELL *a = (Mat_MPISELL *)A->data;
1483: PetscBool flg;
1485: PetscFunctionBegin;
1486: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATMPISELL, &flg));
1487: PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "This function requires a MATMPISELL matrix as input");
1488: if (Ad) *Ad = a->A;
1489: if (Ao) *Ao = a->B;
1490: if (colmap) *colmap = a->garray;
1491: PetscFunctionReturn(PETSC_SUCCESS);
1492: }
1494: /*@
1495: MatMPISELLGetLocalMatCondensed - Creates a `MATSEQSELL` matrix from an `MATMPISELL` matrix by
1496: taking all its local rows and NON-ZERO columns
1498: Not Collective
1500: Input Parameters:
1501: + A - the matrix
1502: . scall - either `MAT_INITIAL_MATRIX` or `MAT_REUSE_MATRIX`
1503: . row - index sets of rows to extract (or `NULL`)
1504: - col - index sets of columns to extract (or `NULL`)
1506: Output Parameter:
1507: . A_loc - the local sequential matrix generated
1509: Level: advanced
1511: .seealso: `Mat`, `MATSEQSELL`, `MATMPISELL`, `MatGetOwnershipRange()`, `MatMPISELLGetLocalMat()`
1512: @*/
1513: PetscErrorCode MatMPISELLGetLocalMatCondensed(Mat A, MatReuse scall, IS *row, IS *col, Mat *A_loc)
1514: {
1515: Mat_MPISELL *a = (Mat_MPISELL *)A->data;
1516: PetscInt i, start, end, ncols, nzA, nzB, *cmap, imark, *idx;
1517: IS isrowa, iscola;
1518: Mat *aloc;
1519: PetscBool match;
1521: PetscFunctionBegin;
1522: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATMPISELL, &match));
1523: PetscCheck(match, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Requires MATMPISELL matrix as input");
1524: PetscCall(PetscLogEventBegin(MAT_Getlocalmatcondensed, A, 0, 0, 0));
1525: if (!row) {
1526: start = A->rmap->rstart;
1527: end = A->rmap->rend;
1528: PetscCall(ISCreateStride(PETSC_COMM_SELF, end - start, start, 1, &isrowa));
1529: } else {
1530: isrowa = *row;
1531: }
1532: if (!col) {
1533: start = A->cmap->rstart;
1534: cmap = a->garray;
1535: nzA = a->A->cmap->n;
1536: nzB = a->B->cmap->n;
1537: PetscCall(PetscMalloc1(nzA + nzB, &idx));
1538: ncols = 0;
1539: for (i = 0; i < nzB; i++) {
1540: if (cmap[i] < start) idx[ncols++] = cmap[i];
1541: else break;
1542: }
1543: imark = i;
1544: for (i = 0; i < nzA; i++) idx[ncols++] = start + i;
1545: for (i = imark; i < nzB; i++) idx[ncols++] = cmap[i];
1546: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, ncols, idx, PETSC_OWN_POINTER, &iscola));
1547: } else {
1548: iscola = *col;
1549: }
1550: if (scall != MAT_INITIAL_MATRIX) {
1551: PetscCall(PetscMalloc1(1, &aloc));
1552: aloc[0] = *A_loc;
1553: }
1554: PetscCall(MatCreateSubMatrices(A, 1, &isrowa, &iscola, scall, &aloc));
1555: *A_loc = aloc[0];
1556: PetscCall(PetscFree(aloc));
1557: if (!row) PetscCall(ISDestroy(&isrowa));
1558: if (!col) PetscCall(ISDestroy(&iscola));
1559: PetscCall(PetscLogEventEnd(MAT_Getlocalmatcondensed, A, 0, 0, 0));
1560: PetscFunctionReturn(PETSC_SUCCESS);
1561: }
1563: #include <../src/mat/impls/aij/mpi/mpiaij.h>
1565: PetscErrorCode MatConvert_MPISELL_MPIAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
1566: {
1567: Mat_MPISELL *a = (Mat_MPISELL *)A->data;
1568: Mat B;
1569: Mat_MPIAIJ *b;
1571: PetscFunctionBegin;
1572: PetscCheck(A->assembled, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Matrix must be assembled");
1574: if (reuse == MAT_REUSE_MATRIX) {
1575: B = *newmat;
1576: } else {
1577: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
1578: PetscCall(MatSetType(B, MATMPIAIJ));
1579: PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
1580: PetscCall(MatSetBlockSizes(B, A->rmap->bs, A->cmap->bs));
1581: PetscCall(MatSeqAIJSetPreallocation(B, 0, NULL));
1582: PetscCall(MatMPIAIJSetPreallocation(B, 0, NULL, 0, NULL));
1583: }
1584: b = (Mat_MPIAIJ *)B->data;
1586: if (reuse == MAT_REUSE_MATRIX) {
1587: PetscCall(MatConvert_SeqSELL_SeqAIJ(a->A, MATSEQAIJ, MAT_REUSE_MATRIX, &b->A));
1588: PetscCall(MatConvert_SeqSELL_SeqAIJ(a->B, MATSEQAIJ, MAT_REUSE_MATRIX, &b->B));
1589: } else {
1590: PetscCall(MatDestroy(&b->A));
1591: PetscCall(MatDestroy(&b->B));
1592: PetscCall(MatDisAssemble_MPISELL(A));
1593: PetscCall(MatConvert_SeqSELL_SeqAIJ(a->A, MATSEQAIJ, MAT_INITIAL_MATRIX, &b->A));
1594: PetscCall(MatConvert_SeqSELL_SeqAIJ(a->B, MATSEQAIJ, MAT_INITIAL_MATRIX, &b->B));
1595: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
1596: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
1597: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
1598: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
1599: }
1601: if (reuse == MAT_INPLACE_MATRIX) {
1602: PetscCall(MatHeaderReplace(A, &B));
1603: } else {
1604: *newmat = B;
1605: }
1606: PetscFunctionReturn(PETSC_SUCCESS);
1607: }
1609: PetscErrorCode MatConvert_MPIAIJ_MPISELL(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
1610: {
1611: Mat_MPIAIJ *a = (Mat_MPIAIJ *)A->data;
1612: Mat B;
1613: Mat_MPISELL *b;
1615: PetscFunctionBegin;
1616: PetscCheck(A->assembled, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Matrix must be assembled");
1618: if (reuse == MAT_REUSE_MATRIX) {
1619: B = *newmat;
1620: } else {
1621: Mat_SeqAIJ *Aa = (Mat_SeqAIJ *)a->A->data, *Ba = (Mat_SeqAIJ *)a->B->data;
1622: PetscInt i, d_nz = 0, o_nz = 0, m = A->rmap->N, n = A->cmap->N, lm = A->rmap->n, ln = A->cmap->n;
1623: PetscInt *d_nnz, *o_nnz;
1624: PetscCall(PetscMalloc2(lm, &d_nnz, lm, &o_nnz));
1625: for (i = 0; i < lm; i++) {
1626: d_nnz[i] = Aa->i[i + 1] - Aa->i[i];
1627: o_nnz[i] = Ba->i[i + 1] - Ba->i[i];
1628: if (d_nnz[i] > d_nz) d_nz = d_nnz[i];
1629: if (o_nnz[i] > o_nz) o_nz = o_nnz[i];
1630: }
1631: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
1632: PetscCall(MatSetType(B, MATMPISELL));
1633: PetscCall(MatSetSizes(B, lm, ln, m, n));
1634: PetscCall(MatSetBlockSizes(B, A->rmap->bs, A->cmap->bs));
1635: PetscCall(MatSeqSELLSetPreallocation(B, d_nz, d_nnz));
1636: PetscCall(MatMPISELLSetPreallocation(B, d_nz, d_nnz, o_nz, o_nnz));
1637: PetscCall(PetscFree2(d_nnz, o_nnz));
1638: }
1639: b = (Mat_MPISELL *)B->data;
1641: if (reuse == MAT_REUSE_MATRIX) {
1642: PetscCall(MatConvert_SeqAIJ_SeqSELL(a->A, MATSEQSELL, MAT_REUSE_MATRIX, &b->A));
1643: PetscCall(MatConvert_SeqAIJ_SeqSELL(a->B, MATSEQSELL, MAT_REUSE_MATRIX, &b->B));
1644: } else {
1645: PetscBool nooffprocentries_A = A->nooffprocentries, nooffprocentries_B = B->nooffprocentries;
1647: PetscCall(MatDestroy(&b->A));
1648: PetscCall(MatDestroy(&b->B));
1649: /* Expand a->B from compacted local off-diag columns back to global columns so the new MPISELL's
1650: MatAssemblyEnd() builds the correct garray/Mvctx for its off-diagonal block. */
1651: PetscCall(MatDisAssemble_MPIAIJ(A, PETSC_FALSE));
1652: PetscCall(MatConvert_SeqAIJ_SeqSELL(a->A, MATSEQSELL, MAT_INITIAL_MATRIX, &b->A));
1653: PetscCall(MatConvert_SeqAIJ_SeqSELL(a->B, MATSEQSELL, MAT_INITIAL_MATRIX, &b->B));
1654: /* The locally-populated A and B have no stashed off-processor entries, so skip the stash scatter. */
1655: A->nooffprocentries = PETSC_TRUE;
1656: B->nooffprocentries = PETSC_TRUE;
1657: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
1658: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
1659: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
1660: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
1661: A->nooffprocentries = nooffprocentries_A;
1662: B->nooffprocentries = nooffprocentries_B;
1663: }
1665: if (reuse == MAT_INPLACE_MATRIX) {
1666: PetscCall(MatHeaderReplace(A, &B));
1667: } else {
1668: *newmat = B;
1669: }
1670: PetscFunctionReturn(PETSC_SUCCESS);
1671: }
1673: PetscErrorCode MatSOR_MPISELL(Mat matin, Vec bb, PetscReal omega, MatSORType flag, PetscReal fshift, PetscInt its, PetscInt lits, Vec xx)
1674: {
1675: Mat_MPISELL *mat = (Mat_MPISELL *)matin->data;
1676: Vec bb1 = NULL;
1678: PetscFunctionBegin;
1679: if (flag == SOR_APPLY_UPPER) {
1680: PetscUseTypeMethod(mat->A, sor, bb, omega, flag, fshift, lits, 1, xx);
1681: PetscFunctionReturn(PETSC_SUCCESS);
1682: }
1684: if (its > 1 || ~flag & SOR_ZERO_INITIAL_GUESS || flag & SOR_EISENSTAT) PetscCall(VecDuplicate(bb, &bb1));
1686: if ((flag & SOR_LOCAL_SYMMETRIC_SWEEP) == SOR_LOCAL_SYMMETRIC_SWEEP) {
1687: if (flag & SOR_ZERO_INITIAL_GUESS) {
1688: PetscUseTypeMethod(mat->A, sor, bb, omega, flag, fshift, lits, 1, xx);
1689: its--;
1690: }
1692: while (its--) {
1693: PetscCall(VecScatterBegin(mat->Mvctx, xx, mat->lvec, INSERT_VALUES, SCATTER_FORWARD));
1694: PetscCall(VecScatterEnd(mat->Mvctx, xx, mat->lvec, INSERT_VALUES, SCATTER_FORWARD));
1696: /* update rhs: bb1 = bb - B*x */
1697: PetscCall(VecScale(mat->lvec, -1.0));
1698: PetscUseTypeMethod(mat->B, multadd, mat->lvec, bb, bb1);
1700: /* local sweep */
1701: PetscUseTypeMethod(mat->A, sor, bb1, omega, SOR_SYMMETRIC_SWEEP, fshift, lits, 1, xx);
1702: }
1703: } else if (flag & SOR_LOCAL_FORWARD_SWEEP) {
1704: if (flag & SOR_ZERO_INITIAL_GUESS) {
1705: PetscUseTypeMethod(mat->A, sor, bb, omega, flag, fshift, lits, 1, xx);
1706: its--;
1707: }
1708: while (its--) {
1709: PetscCall(VecScatterBegin(mat->Mvctx, xx, mat->lvec, INSERT_VALUES, SCATTER_FORWARD));
1710: PetscCall(VecScatterEnd(mat->Mvctx, xx, mat->lvec, INSERT_VALUES, SCATTER_FORWARD));
1712: /* update rhs: bb1 = bb - B*x */
1713: PetscCall(VecScale(mat->lvec, -1.0));
1714: PetscUseTypeMethod(mat->B, multadd, mat->lvec, bb, bb1);
1716: /* local sweep */
1717: PetscUseTypeMethod(mat->A, sor, bb1, omega, SOR_FORWARD_SWEEP, fshift, lits, 1, xx);
1718: }
1719: } else if (flag & SOR_LOCAL_BACKWARD_SWEEP) {
1720: if (flag & SOR_ZERO_INITIAL_GUESS) {
1721: PetscUseTypeMethod(mat->A, sor, bb, omega, flag, fshift, lits, 1, xx);
1722: its--;
1723: }
1724: while (its--) {
1725: PetscCall(VecScatterBegin(mat->Mvctx, xx, mat->lvec, INSERT_VALUES, SCATTER_FORWARD));
1726: PetscCall(VecScatterEnd(mat->Mvctx, xx, mat->lvec, INSERT_VALUES, SCATTER_FORWARD));
1728: /* update rhs: bb1 = bb - B*x */
1729: PetscCall(VecScale(mat->lvec, -1.0));
1730: PetscUseTypeMethod(mat->B, multadd, mat->lvec, bb, bb1);
1732: /* local sweep */
1733: PetscUseTypeMethod(mat->A, sor, bb1, omega, SOR_BACKWARD_SWEEP, fshift, lits, 1, xx);
1734: }
1735: } else SETERRQ(PetscObjectComm((PetscObject)matin), PETSC_ERR_SUP, "Parallel SOR not supported");
1737: PetscCall(VecDestroy(&bb1));
1739: matin->factorerrortype = mat->A->factorerrortype;
1740: PetscFunctionReturn(PETSC_SUCCESS);
1741: }
1743: #if PetscDefined(HAVE_CUDA)
1744: PETSC_INTERN PetscErrorCode MatConvert_MPISELL_MPISELLCUDA(Mat, MatType, MatReuse, Mat *);
1745: #endif
1747: /*MC
1748: MATMPISELL - MATMPISELL = "MPISELL" - A matrix type to be used for parallel sparse matrices.
1750: Options Database Keys:
1751: . -mat_type mpisell - sets the matrix type to `MATMPISELL` during a call to `MatSetFromOptions()`
1753: Level: beginner
1755: .seealso: `Mat`, `MATSELL`, `MATSEQSELL`, `MatCreateSELL()`
1756: M*/
1757: PETSC_EXTERN PetscErrorCode MatCreate_MPISELL(Mat B)
1758: {
1759: Mat_MPISELL *b;
1760: PetscMPIInt size;
1762: PetscFunctionBegin;
1763: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)B), &size));
1764: PetscCall(PetscNew(&b));
1765: B->data = (void *)b;
1766: B->ops[0] = MatOps_Values;
1767: B->assembled = PETSC_FALSE;
1768: B->insertmode = NOT_SET_VALUES;
1769: b->size = size;
1770: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)B), &b->rank));
1771: /* build cache for off array entries formed */
1772: PetscCall(MatStashCreate_Private(PetscObjectComm((PetscObject)B), 1, &B->stash));
1774: b->donotstash = PETSC_FALSE;
1775: b->colmap = NULL;
1776: b->garray = NULL;
1777: b->roworiented = PETSC_TRUE;
1779: /* stuff used for matrix vector multiply */
1780: b->lvec = NULL;
1781: b->Mvctx = NULL;
1783: /* stuff for MatGetRow() */
1784: b->rowindices = NULL;
1785: b->rowvalues = NULL;
1786: b->getrowactive = PETSC_FALSE;
1788: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatStoreValues_C", MatStoreValues_MPISELL));
1789: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatRetrieveValues_C", MatRetrieveValues_MPISELL));
1790: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatIsTranspose_C", MatIsTranspose_MPISELL));
1791: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatMPISELLSetPreallocation_C", MatMPISELLSetPreallocation_MPISELL));
1792: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_mpisell_mpiaij_C", MatConvert_MPISELL_MPIAIJ));
1793: #if PetscDefined(HAVE_CUDA)
1794: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_mpisell_mpisellcuda_C", MatConvert_MPISELL_MPISELLCUDA));
1795: #endif
1796: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDiagonalScaleLocal_C", MatDiagonalScaleLocal_MPISELL));
1797: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatGetMultPetscSF_C", MatGetMultPetscSF_MPISELL));
1798: PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATMPISELL));
1799: PetscFunctionReturn(PETSC_SUCCESS);
1800: }