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, &notme);
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: }