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