Actual source code: mpibaij.c

  1: #include <../src/mat/impls/baij/mpi/mpibaij.h>

  3: #include <petsc/private/hashseti.h>
  4: #include <petscblaslapack.h>
  5: #include <petscsf.h>

  7: static PetscErrorCode MatDestroy_MPIBAIJ(Mat mat)
  8: {
  9:   Mat_MPIBAIJ *baij = (Mat_MPIBAIJ *)mat->data;

 11:   PetscFunctionBegin;
 12:   PetscCall(PetscLogObjectState((PetscObject)mat, "Rows=%" PetscInt_FMT ",Cols=%" PetscInt_FMT, mat->rmap->N, mat->cmap->N));
 13:   PetscCall(MatStashDestroy_Private(&mat->stash));
 14:   PetscCall(MatStashDestroy_Private(&mat->bstash));
 15:   PetscCall(MatDestroy(&baij->A));
 16:   PetscCall(MatDestroy(&baij->B));
 17: #if defined(PETSC_USE_CTABLE)
 18:   PetscCall(PetscHMapIDestroy(&baij->colmap));
 19: #else
 20:   PetscCall(PetscFree(baij->colmap));
 21: #endif
 22:   PetscCall(PetscFree(baij->garray));
 23:   PetscCall(VecDestroy(&baij->lvec));
 24:   PetscCall(VecScatterDestroy(&baij->Mvctx));
 25:   PetscCall(PetscFree2(baij->rowvalues, baij->rowindices));
 26:   PetscCall(PetscFree(baij->barray));
 27:   PetscCall(PetscFree2(baij->hd, baij->ht));
 28:   PetscCall(PetscFree(baij->rangebs));
 29:   PetscCall(PetscFree(mat->data));

 31:   PetscCall(PetscObjectChangeTypeName((PetscObject)mat, NULL));
 32:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatStoreValues_C", NULL));
 33:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatRetrieveValues_C", NULL));
 34:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMPIBAIJSetPreallocation_C", NULL));
 35:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMPIBAIJSetPreallocationCSR_C", NULL));
 36:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDiagonalScaleLocal_C", NULL));
 37:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatSetHashTableFactor_C", NULL));
 38:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_mpibaij_mpisbaij_C", NULL));
 39:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_mpibaij_mpiadj_C", NULL));
 40:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_mpibaij_mpiaij_C", NULL));
 41: #if defined(PETSC_HAVE_HYPRE)
 42:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_mpibaij_hypre_C", NULL));
 43: #endif
 44:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_mpibaij_is_C", NULL));
 45:   PetscFunctionReturn(PETSC_SUCCESS);
 46: }

 48: /* defines MatSetValues_MPI_Hash(), MatAssemblyBegin_MPI_Hash(), and  MatAssemblyEnd_MPI_Hash() */
 49: #define TYPE BAIJ
 50: #include "../src/mat/impls/aij/mpi/mpihashmat.h"
 51: #undef TYPE

 53: #if defined(PETSC_HAVE_HYPRE)
 54: PETSC_INTERN PetscErrorCode MatConvert_AIJ_HYPRE(Mat, MatType, MatReuse, Mat *);
 55: #endif

 57: static PetscErrorCode MatGetRowMaxAbs_MPIBAIJ(Mat A, Vec v, PetscInt idx[])
 58: {
 59:   Mat_MPIBAIJ       *a = (Mat_MPIBAIJ *)A->data;
 60:   PetscInt           i, *idxb = NULL, m = A->rmap->n, bs = A->cmap->bs;
 61:   PetscScalar       *va, *vv;
 62:   Vec                vB, vA;
 63:   const PetscScalar *vb;

 65:   PetscFunctionBegin;
 66:   PetscCall(MatCreateVecs(a->A, NULL, &vA));
 67:   PetscCall(MatGetRowMaxAbs(a->A, vA, idx));

 69:   PetscCall(VecGetArrayWrite(vA, &va));
 70:   if (idx) {
 71:     for (i = 0; i < m; i++) {
 72:       if (PetscAbsScalar(va[i])) idx[i] += A->cmap->rstart;
 73:     }
 74:   }

 76:   PetscCall(MatCreateVecs(a->B, NULL, &vB));
 77:   PetscCall(PetscMalloc1(m, &idxb));
 78:   PetscCall(MatGetRowMaxAbs(a->B, vB, idxb));

 80:   PetscCall(VecGetArrayWrite(v, &vv));
 81:   PetscCall(VecGetArrayRead(vB, &vb));
 82:   for (i = 0; i < m; i++) {
 83:     if (PetscAbsScalar(va[i]) < PetscAbsScalar(vb[i])) {
 84:       vv[i] = vb[i];
 85:       if (idx) idx[i] = bs * a->garray[idxb[i] / bs] + (idxb[i] % bs);
 86:     } else {
 87:       vv[i] = va[i];
 88:       if (idx && PetscAbsScalar(va[i]) == PetscAbsScalar(vb[i]) && idxb[i] != -1 && idx[i] > bs * a->garray[idxb[i] / bs] + (idxb[i] % bs)) idx[i] = bs * a->garray[idxb[i] / bs] + (idxb[i] % bs);
 89:     }
 90:   }
 91:   PetscCall(VecRestoreArrayWrite(vA, &vv));
 92:   PetscCall(VecRestoreArrayWrite(vA, &va));
 93:   PetscCall(VecRestoreArrayRead(vB, &vb));
 94:   PetscCall(PetscFree(idxb));
 95:   PetscCall(VecDestroy(&vA));
 96:   PetscCall(VecDestroy(&vB));
 97:   PetscFunctionReturn(PETSC_SUCCESS);
 98: }

100: static PetscErrorCode MatGetRowSumAbs_MPIBAIJ(Mat A, Vec v)
101: {
102:   Mat_MPIBAIJ *a = (Mat_MPIBAIJ *)A->data;
103:   Vec          vB, vA;

105:   PetscFunctionBegin;
106:   PetscCall(MatCreateVecs(a->A, NULL, &vA));
107:   PetscCall(MatGetRowSumAbs(a->A, vA));
108:   PetscCall(MatCreateVecs(a->B, NULL, &vB));
109:   PetscCall(MatGetRowSumAbs(a->B, vB));
110:   PetscCall(VecAXPY(vA, 1.0, vB));
111:   PetscCall(VecDestroy(&vB));
112:   PetscCall(VecCopy(vA, v));
113:   PetscCall(VecDestroy(&vA));
114:   PetscFunctionReturn(PETSC_SUCCESS);
115: }

117: static PetscErrorCode MatStoreValues_MPIBAIJ(Mat mat)
118: {
119:   Mat_MPIBAIJ *aij = (Mat_MPIBAIJ *)mat->data;

121:   PetscFunctionBegin;
122:   PetscCall(MatStoreValues(aij->A));
123:   PetscCall(MatStoreValues(aij->B));
124:   PetscFunctionReturn(PETSC_SUCCESS);
125: }

127: static PetscErrorCode MatRetrieveValues_MPIBAIJ(Mat mat)
128: {
129:   Mat_MPIBAIJ *aij = (Mat_MPIBAIJ *)mat->data;

131:   PetscFunctionBegin;
132:   PetscCall(MatRetrieveValues(aij->A));
133:   PetscCall(MatRetrieveValues(aij->B));
134:   PetscFunctionReturn(PETSC_SUCCESS);
135: }

137: /*
138:      Local utility routine that creates a mapping from the global column
139:    number to the local number in the off-diagonal part of the local
140:    storage of the matrix.  This is done in a non scalable way since the
141:    length of colmap equals the global matrix length.
142: */
143: PetscErrorCode MatCreateColmap_MPIBAIJ_Private(Mat mat)
144: {
145:   Mat_MPIBAIJ *baij = (Mat_MPIBAIJ *)mat->data;
146:   Mat_SeqBAIJ *B    = (Mat_SeqBAIJ *)baij->B->data;
147:   PetscInt     nbs = B->nbs, i, bs = mat->rmap->bs;

149:   PetscFunctionBegin;
150: #if defined(PETSC_USE_CTABLE)
151:   PetscCall(PetscHMapICreateWithSize(baij->nbs, &baij->colmap));
152:   for (i = 0; i < nbs; i++) PetscCall(PetscHMapISet(baij->colmap, baij->garray[i] + 1, i * bs + 1));
153: #else
154:   PetscCall(PetscCalloc1(baij->Nbs + 1, &baij->colmap));
155:   for (i = 0; i < nbs; i++) baij->colmap[baij->garray[i]] = i * bs + 1;
156: #endif
157:   PetscFunctionReturn(PETSC_SUCCESS);
158: }

160: #define MatSetValues_SeqBAIJ_A_Private(row, col, value, addv, orow, ocol) \
161:   do { \
162:     brow = row / bs; \
163:     rp   = PetscSafePointerPlusOffset(aj, ai[brow]); \
164:     ap   = PetscSafePointerPlusOffset(aa, bs2 * ai[brow]); \
165:     rmax = aimax[brow]; \
166:     nrow = ailen[brow]; \
167:     bcol = col / bs; \
168:     ridx = row % bs; \
169:     cidx = col % bs; \
170:     low  = 0; \
171:     high = nrow; \
172:     while (high - low > 3) { \
173:       t = (low + high) / 2; \
174:       if (rp[t] > bcol) high = t; \
175:       else low = t; \
176:     } \
177:     for (_i = low; _i < high; _i++) { \
178:       if (rp[_i] > bcol) break; \
179:       if (rp[_i] == bcol) { \
180:         bap = ap + bs2 * _i + bs * cidx + ridx; \
181:         if (addv == ADD_VALUES) *bap += value; \
182:         else *bap = value; \
183:         goto a_noinsert; \
184:       } \
185:     } \
186:     if (a->nonew == 1) goto a_noinsert; \
187:     PetscCheck(a->nonew != -1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Inserting a new nonzero at global row/column (%" PetscInt_FMT ", %" PetscInt_FMT ") into matrix", orow, ocol); \
188:     MatSeqXAIJReallocateAIJ(A, a->mbs, bs2, nrow, brow, bcol, rmax, aa, ai, aj, rp, ap, aimax, a->nonew, MatScalar); \
189:     N = nrow++ - 1; \
190:     /* shift up all the later entries in this row */ \
191:     PetscCall(PetscArraymove(rp + _i + 1, rp + _i, N - _i + 1)); \
192:     PetscCall(PetscArraymove(ap + bs2 * (_i + 1), ap + bs2 * _i, bs2 * (N - _i + 1))); \
193:     PetscCall(PetscArrayzero(ap + bs2 * _i, bs2)); \
194:     rp[_i]                          = bcol; \
195:     ap[bs2 * _i + bs * cidx + ridx] = value; \
196:   a_noinsert:; \
197:     ailen[brow] = nrow; \
198:   } while (0)

200: #define MatSetValues_SeqBAIJ_B_Private(row, col, value, addv, orow, ocol) \
201:   do { \
202:     brow = row / bs; \
203:     rp   = PetscSafePointerPlusOffset(bj, bi[brow]); \
204:     ap   = PetscSafePointerPlusOffset(ba, bs2 * bi[brow]); \
205:     rmax = bimax[brow]; \
206:     nrow = bilen[brow]; \
207:     bcol = col / bs; \
208:     ridx = row % bs; \
209:     cidx = col % bs; \
210:     low  = 0; \
211:     high = nrow; \
212:     while (high - low > 3) { \
213:       t = (low + high) / 2; \
214:       if (rp[t] > bcol) high = t; \
215:       else low = t; \
216:     } \
217:     for (_i = low; _i < high; _i++) { \
218:       if (rp[_i] > bcol) break; \
219:       if (rp[_i] == bcol) { \
220:         bap = ap + bs2 * _i + bs * cidx + ridx; \
221:         if (addv == ADD_VALUES) *bap += value; \
222:         else *bap = value; \
223:         goto b_noinsert; \
224:       } \
225:     } \
226:     if (b->nonew == 1) goto b_noinsert; \
227:     PetscCheck(b->nonew != -1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Inserting a new nonzero at global row/column  (%" PetscInt_FMT ", %" PetscInt_FMT ") into matrix", orow, ocol); \
228:     MatSeqXAIJReallocateAIJ(B, b->mbs, bs2, nrow, brow, bcol, rmax, ba, bi, bj, rp, ap, bimax, b->nonew, MatScalar); \
229:     N = nrow++ - 1; \
230:     /* shift up all the later entries in this row */ \
231:     PetscCall(PetscArraymove(rp + _i + 1, rp + _i, N - _i + 1)); \
232:     PetscCall(PetscArraymove(ap + bs2 * (_i + 1), ap + bs2 * _i, bs2 * (N - _i + 1))); \
233:     PetscCall(PetscArrayzero(ap + bs2 * _i, bs2)); \
234:     rp[_i]                          = bcol; \
235:     ap[bs2 * _i + bs * cidx + ridx] = value; \
236:   b_noinsert:; \
237:     bilen[brow] = nrow; \
238:   } while (0)

240: PetscErrorCode MatSetValues_MPIBAIJ(Mat mat, PetscInt m, const PetscInt im[], PetscInt n, const PetscInt in[], const PetscScalar v[], InsertMode addv)
241: {
242:   Mat_MPIBAIJ *baij = (Mat_MPIBAIJ *)mat->data;
243:   MatScalar    value;
244:   PetscBool    roworiented = baij->roworiented;
245:   PetscInt     i, j, row, col;
246:   PetscInt     rstart_orig = mat->rmap->rstart;
247:   PetscInt     rend_orig = mat->rmap->rend, cstart_orig = mat->cmap->rstart;
248:   PetscInt     cend_orig = mat->cmap->rend, bs = mat->rmap->bs;

250:   /* Some Variables required in the macro */
251:   Mat          A     = baij->A;
252:   Mat_SeqBAIJ *a     = (Mat_SeqBAIJ *)A->data;
253:   PetscInt    *aimax = a->imax, *ai = a->i, *ailen = a->ilen, *aj = a->j;
254:   MatScalar   *aa = a->a;

256:   Mat          B     = baij->B;
257:   Mat_SeqBAIJ *b     = (Mat_SeqBAIJ *)B->data;
258:   PetscInt    *bimax = b->imax, *bi = b->i, *bilen = b->ilen, *bj = b->j;
259:   MatScalar   *ba = b->a;

261:   PetscInt  *rp, ii, nrow, _i, rmax, N, brow, bcol;
262:   PetscInt   low, high, t, ridx, cidx, bs2 = a->bs2;
263:   MatScalar *ap, *bap;

265:   PetscFunctionBegin;
266:   for (i = 0; i < m; i++) {
267:     if (im[i] < 0) continue;
268:     PetscCheck(im[i] < mat->rmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, im[i], mat->rmap->N - 1);
269:     if (im[i] >= rstart_orig && im[i] < rend_orig) {
270:       row = im[i] - rstart_orig;
271:       for (j = 0; j < n; j++) {
272:         if (in[j] >= cstart_orig && in[j] < cend_orig) {
273:           col = in[j] - cstart_orig;
274:           if (roworiented) value = v[i * n + j];
275:           else value = v[i + j * m];
276:           MatSetValues_SeqBAIJ_A_Private(row, col, value, addv, im[i], in[j]);
277:         } else if (in[j] < 0) {
278:           continue;
279:         } else {
280:           PetscCheck(in[j] < mat->cmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Column too large: col %" PetscInt_FMT " max %" PetscInt_FMT, in[j], mat->cmap->N - 1);
281:           if (mat->was_assembled) {
282:             if (!baij->colmap) PetscCall(MatCreateColmap_MPIBAIJ_Private(mat));
283: #if defined(PETSC_USE_CTABLE)
284:             PetscCall(PetscHMapIGetWithDefault(baij->colmap, in[j] / bs + 1, 0, &col));
285:             col = col - 1;
286: #else
287:             col = baij->colmap[in[j] / bs] - 1;
288: #endif
289:             if (col < 0 && !((Mat_SeqBAIJ *)baij->B->data)->nonew) {
290:               PetscCall(MatDisAssemble_MPIBAIJ(mat));
291:               col = in[j];
292:               /* Reinitialize the variables required by MatSetValues_SeqBAIJ_B_Private() */
293:               B     = baij->B;
294:               b     = (Mat_SeqBAIJ *)B->data;
295:               bimax = b->imax;
296:               bi    = b->i;
297:               bilen = b->ilen;
298:               bj    = b->j;
299:               ba    = b->a;
300:             } else {
301:               PetscCheck(col >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Inserting a new nonzero (%" PetscInt_FMT ", %" PetscInt_FMT ") into matrix", im[i], in[j]);
302:               col += in[j] % bs;
303:             }
304:           } else col = in[j];
305:           if (roworiented) value = v[i * n + j];
306:           else value = v[i + j * m];
307:           MatSetValues_SeqBAIJ_B_Private(row, col, value, addv, im[i], in[j]);
308:           /* PetscCall(MatSetValues_SeqBAIJ(baij->B,1,&row,1,&col,&value,addv)); */
309:         }
310:       }
311:     } else {
312:       PetscCheck(!mat->nooffprocentries, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Setting off process row %" PetscInt_FMT " even though MatSetOption(,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE) was set", im[i]);
313:       if (!baij->donotstash) {
314:         mat->assembled = PETSC_FALSE;
315:         if (roworiented) {
316:           PetscCall(MatStashValuesRow_Private(&mat->stash, im[i], n, in, v + i * n, PETSC_FALSE));
317:         } else {
318:           PetscCall(MatStashValuesCol_Private(&mat->stash, im[i], n, in, v + i, m, PETSC_FALSE));
319:         }
320:       }
321:     }
322:   }
323:   PetscFunctionReturn(PETSC_SUCCESS);
324: }

326: static inline PetscErrorCode MatSetValuesBlocked_SeqBAIJ_Inlined(Mat A, PetscInt row, PetscInt col, const PetscScalar v[], InsertMode is, PetscInt orow, PetscInt ocol)
327: {
328:   Mat_SeqBAIJ       *a = (Mat_SeqBAIJ *)A->data;
329:   PetscInt          *rp, low, high, t, ii, jj, nrow, i, rmax, N;
330:   PetscInt          *imax = a->imax, *ai = a->i, *ailen = a->ilen;
331:   PetscInt          *aj = a->j, nonew = a->nonew, bs2 = a->bs2, bs = A->rmap->bs;
332:   PetscBool          roworiented = a->roworiented;
333:   const PetscScalar *value       = v;
334:   MatScalar         *ap, *aa = a->a, *bap;

336:   PetscFunctionBegin;
337:   rp    = aj + ai[row];
338:   ap    = aa + bs2 * ai[row];
339:   rmax  = imax[row];
340:   nrow  = ailen[row];
341:   value = v;
342:   low   = 0;
343:   high  = nrow;
344:   while (high - low > 7) {
345:     t = (low + high) / 2;
346:     if (rp[t] > col) high = t;
347:     else low = t;
348:   }
349:   for (i = low; i < high; i++) {
350:     if (rp[i] > col) break;
351:     if (rp[i] == col) {
352:       bap = ap + bs2 * i;
353:       if (roworiented) {
354:         if (is == ADD_VALUES) {
355:           for (ii = 0; ii < bs; ii++) {
356:             for (jj = ii; jj < bs2; jj += bs) bap[jj] += *value++;
357:           }
358:         } else {
359:           for (ii = 0; ii < bs; ii++) {
360:             for (jj = ii; jj < bs2; jj += bs) bap[jj] = *value++;
361:           }
362:         }
363:       } else {
364:         if (is == ADD_VALUES) {
365:           for (ii = 0; ii < bs; ii++, value += bs) {
366:             for (jj = 0; jj < bs; jj++) bap[jj] += value[jj];
367:             bap += bs;
368:           }
369:         } else {
370:           for (ii = 0; ii < bs; ii++, value += bs) {
371:             for (jj = 0; jj < bs; jj++) bap[jj] = value[jj];
372:             bap += bs;
373:           }
374:         }
375:       }
376:       goto noinsert2;
377:     }
378:   }
379:   if (nonew == 1) goto noinsert2;
380:   PetscCheck(nonew != -1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Inserting a new global block indexed nonzero block (%" PetscInt_FMT ", %" PetscInt_FMT ") in the matrix", orow, ocol);
381:   MatSeqXAIJReallocateAIJ(A, a->mbs, bs2, nrow, row, col, rmax, aa, ai, aj, rp, ap, imax, nonew, MatScalar);
382:   N = nrow++ - 1;
383:   high++;
384:   /* shift up all the later entries in this row */
385:   PetscCall(PetscArraymove(rp + i + 1, rp + i, N - i + 1));
386:   PetscCall(PetscArraymove(ap + bs2 * (i + 1), ap + bs2 * i, bs2 * (N - i + 1)));
387:   rp[i] = col;
388:   bap   = ap + bs2 * i;
389:   if (roworiented) {
390:     for (ii = 0; ii < bs; ii++) {
391:       for (jj = ii; jj < bs2; jj += bs) bap[jj] = *value++;
392:     }
393:   } else {
394:     for (ii = 0; ii < bs; ii++) {
395:       for (jj = 0; jj < bs; jj++) *bap++ = *value++;
396:     }
397:   }
398: noinsert2:;
399:   ailen[row] = nrow;
400:   PetscFunctionReturn(PETSC_SUCCESS);
401: }

403: /*
404:     This routine should be optimized so that the block copy at ** Here a copy is required ** below is not needed
405:     by passing additional stride information into the MatSetValuesBlocked_SeqBAIJ_Inlined() routine
406: */
407: static PetscErrorCode MatSetValuesBlocked_MPIBAIJ(Mat mat, PetscInt m, const PetscInt im[], PetscInt n, const PetscInt in[], const PetscScalar v[], InsertMode addv)
408: {
409:   Mat_MPIBAIJ       *baij = (Mat_MPIBAIJ *)mat->data;
410:   const PetscScalar *value;
411:   MatScalar         *barray      = baij->barray;
412:   PetscBool          roworiented = baij->roworiented;
413:   PetscInt           i, j, ii, jj, row, col, rstart = baij->rstartbs;
414:   PetscInt           rend = baij->rendbs, cstart = baij->cstartbs, stepval;
415:   PetscInt           cend = baij->cendbs, bs = mat->rmap->bs, bs2 = baij->bs2;

417:   PetscFunctionBegin;
418:   if (!barray) {
419:     PetscCall(PetscMalloc1(bs2, &barray));
420:     baij->barray = barray;
421:   }

423:   if (roworiented) stepval = (n - 1) * bs;
424:   else stepval = (m - 1) * bs;

426:   for (i = 0; i < m; i++) {
427:     if (im[i] < 0) continue;
428:     PetscCheck(im[i] < baij->Mbs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Block indexed row too large %" PetscInt_FMT " max %" PetscInt_FMT, im[i], baij->Mbs - 1);
429:     if (im[i] >= rstart && im[i] < rend) {
430:       row = im[i] - rstart;
431:       for (j = 0; j < n; j++) {
432:         /* If NumCol = 1 then a copy is not required */
433:         if ((roworiented) && (n == 1)) {
434:           barray = (MatScalar *)v + i * bs2;
435:         } else if ((!roworiented) && (m == 1)) {
436:           barray = (MatScalar *)v + j * bs2;
437:         } else { /* Here a copy is required */
438:           if (roworiented) {
439:             value = v + (i * (stepval + bs) + j) * bs;
440:           } else {
441:             value = v + (j * (stepval + bs) + i) * bs;
442:           }
443:           for (ii = 0; ii < bs; ii++, value += bs + stepval) {
444:             for (jj = 0; jj < bs; jj++) barray[jj] = value[jj];
445:             barray += bs;
446:           }
447:           barray -= bs2;
448:         }

450:         if (in[j] >= cstart && in[j] < cend) {
451:           col = in[j] - cstart;
452:           PetscCall(MatSetValuesBlocked_SeqBAIJ_Inlined(baij->A, row, col, barray, addv, im[i], in[j]));
453:         } else if (in[j] < 0) {
454:           continue;
455:         } else {
456:           PetscCheck(in[j] < baij->Nbs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Block indexed column too large %" PetscInt_FMT " max %" PetscInt_FMT, in[j], baij->Nbs - 1);
457:           if (mat->was_assembled) {
458:             if (!baij->colmap) PetscCall(MatCreateColmap_MPIBAIJ_Private(mat));

460: #if defined(PETSC_USE_CTABLE)
461:             PetscCall(PetscHMapIGetWithDefault(baij->colmap, in[j] + 1, 0, &col));
462:             col = col < 1 ? -1 : (col - 1) / bs;
463: #else
464:             col = baij->colmap[in[j]] < 1 ? -1 : (baij->colmap[in[j]] - 1) / bs;
465: #endif
466:             if (col < 0 && !((Mat_SeqBAIJ *)baij->B->data)->nonew) {
467:               PetscCall(MatDisAssemble_MPIBAIJ(mat));
468:               col = in[j];
469:             } else PetscCheck(col >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Inserting a new blocked indexed nonzero block (%" PetscInt_FMT ", %" PetscInt_FMT ") into matrix", im[i], in[j]);
470:           } else col = in[j];
471:           PetscCall(MatSetValuesBlocked_SeqBAIJ_Inlined(baij->B, row, col, barray, addv, im[i], in[j]));
472:         }
473:       }
474:     } else {
475:       PetscCheck(!mat->nooffprocentries, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Setting off process block indexed row %" PetscInt_FMT " even though MatSetOption(,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE) was set", im[i]);
476:       if (!baij->donotstash) {
477:         if (roworiented) {
478:           PetscCall(MatStashValuesRowBlocked_Private(&mat->bstash, im[i], n, in, v, m, n, i));
479:         } else {
480:           PetscCall(MatStashValuesColBlocked_Private(&mat->bstash, im[i], n, in, v, m, n, i));
481:         }
482:       }
483:     }
484:   }
485:   PetscFunctionReturn(PETSC_SUCCESS);
486: }

488: #define HASH_KEY             0.6180339887
489: #define HASH(size, key, tmp) (tmp = (key) * HASH_KEY, (PetscInt)((size) * (tmp - (PetscInt)tmp)))
490: /* #define HASH(size,key) ((PetscInt)((size)*fmod(((key)*HASH_KEY),1))) */
491: /* #define HASH(size,key,tmp) ((PetscInt)((size)*fmod(((key)*HASH_KEY),1))) */
492: static PetscErrorCode MatSetValues_MPIBAIJ_HT(Mat mat, PetscInt m, const PetscInt im[], PetscInt n, const PetscInt in[], const PetscScalar v[], InsertMode addv)
493: {
494:   Mat_MPIBAIJ *baij        = (Mat_MPIBAIJ *)mat->data;
495:   PetscBool    roworiented = baij->roworiented;
496:   PetscInt     i, j, row, col;
497:   PetscInt     rstart_orig = mat->rmap->rstart;
498:   PetscInt     rend_orig = mat->rmap->rend, Nbs = baij->Nbs;
499:   PetscInt     h1, key, size = baij->ht_size, bs = mat->rmap->bs, *HT = baij->ht, idx;
500:   PetscReal    tmp;
501:   MatScalar  **HD       = baij->hd, value;
502:   PetscInt     total_ct = baij->ht_total_ct, insert_ct = baij->ht_insert_ct;

504:   PetscFunctionBegin;
505:   for (i = 0; i < m; i++) {
506:     if (PetscDefined(USE_DEBUG)) {
507:       PetscCheck(im[i] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative row");
508:       PetscCheck(im[i] < mat->rmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, im[i], mat->rmap->N - 1);
509:     }
510:     row = im[i];
511:     if (row >= rstart_orig && row < rend_orig) {
512:       for (j = 0; j < n; j++) {
513:         col = in[j];
514:         if (roworiented) value = v[i * n + j];
515:         else value = v[i + j * m];
516:         /* Look up PetscInto the Hash Table */
517:         key = (row / bs) * Nbs + (col / bs) + 1;
518:         h1  = HASH(size, key, tmp);

520:         idx = h1;
521:         if (PetscDefined(USE_DEBUG)) {
522:           insert_ct++;
523:           total_ct++;
524:           if (HT[idx] != key) {
525:             for (idx = h1; (idx < size) && (HT[idx] != key); idx++, total_ct++);
526:             if (idx == size) {
527:               for (idx = 0; (idx < h1) && (HT[idx] != key); idx++, total_ct++);
528:               PetscCheck(idx != h1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "(%" PetscInt_FMT ",%" PetscInt_FMT ") has no entry in the hash table", row, col);
529:             }
530:           }
531:         } else if (HT[idx] != key) {
532:           for (idx = h1; (idx < size) && (HT[idx] != key); idx++);
533:           if (idx == size) {
534:             for (idx = 0; (idx < h1) && (HT[idx] != key); idx++);
535:             PetscCheck(idx != h1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "(%" PetscInt_FMT ",%" PetscInt_FMT ") has no entry in the hash table", row, col);
536:           }
537:         }
538:         /* A HASH table entry is found, so insert the values at the correct address */
539:         if (addv == ADD_VALUES) *(HD[idx] + (col % bs) * bs + (row % bs)) += value;
540:         else *(HD[idx] + (col % bs) * bs + (row % bs)) = value;
541:       }
542:     } else if (!baij->donotstash) {
543:       if (roworiented) {
544:         PetscCall(MatStashValuesRow_Private(&mat->stash, im[i], n, in, v + i * n, PETSC_FALSE));
545:       } else {
546:         PetscCall(MatStashValuesCol_Private(&mat->stash, im[i], n, in, v + i, m, PETSC_FALSE));
547:       }
548:     }
549:   }
550:   if (PetscDefined(USE_DEBUG)) {
551:     baij->ht_total_ct += total_ct;
552:     baij->ht_insert_ct += insert_ct;
553:   }
554:   PetscFunctionReturn(PETSC_SUCCESS);
555: }

557: static PetscErrorCode MatSetValuesBlocked_MPIBAIJ_HT(Mat mat, PetscInt m, const PetscInt im[], PetscInt n, const PetscInt in[], const PetscScalar v[], InsertMode addv)
558: {
559:   Mat_MPIBAIJ       *baij        = (Mat_MPIBAIJ *)mat->data;
560:   PetscBool          roworiented = baij->roworiented;
561:   PetscInt           i, j, ii, jj, row, col;
562:   PetscInt           rstart = baij->rstartbs;
563:   PetscInt           rend = mat->rmap->rend, stepval, bs = mat->rmap->bs, bs2 = baij->bs2, nbs2 = n * bs2;
564:   PetscInt           h1, key, size = baij->ht_size, idx, *HT = baij->ht, Nbs = baij->Nbs;
565:   PetscReal          tmp;
566:   MatScalar        **HD = baij->hd, *baij_a;
567:   const PetscScalar *v_t, *value;
568:   PetscInt           total_ct = baij->ht_total_ct, insert_ct = baij->ht_insert_ct;

570:   PetscFunctionBegin;
571:   if (roworiented) stepval = (n - 1) * bs;
572:   else stepval = (m - 1) * bs;

574:   for (i = 0; i < m; i++) {
575:     if (PetscDefined(USE_DEBUG)) {
576:       PetscCheck(im[i] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative row: %" PetscInt_FMT, im[i]);
577:       PetscCheck(im[i] < baij->Mbs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, im[i], baij->Mbs - 1);
578:     }
579:     row = im[i];
580:     v_t = v + i * nbs2;
581:     if (row >= rstart && row < rend) {
582:       for (j = 0; j < n; j++) {
583:         col = in[j];

585:         /* Look up into the Hash Table */
586:         key = row * Nbs + col + 1;
587:         h1  = HASH(size, key, tmp);

589:         idx = h1;
590:         if (PetscDefined(USE_DEBUG)) {
591:           total_ct++;
592:           insert_ct++;
593:           if (HT[idx] != key) {
594:             for (idx = h1; (idx < size) && (HT[idx] != key); idx++, total_ct++);
595:             if (idx == size) {
596:               for (idx = 0; (idx < h1) && (HT[idx] != key); idx++, total_ct++);
597:               PetscCheck(idx != h1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "(%" PetscInt_FMT ",%" PetscInt_FMT ") has no entry in the hash table", row, col);
598:             }
599:           }
600:         } else if (HT[idx] != key) {
601:           for (idx = h1; (idx < size) && (HT[idx] != key); idx++);
602:           if (idx == size) {
603:             for (idx = 0; (idx < h1) && (HT[idx] != key); idx++);
604:             PetscCheck(idx != h1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "(%" PetscInt_FMT ",%" PetscInt_FMT ") has no entry in the hash table", row, col);
605:           }
606:         }
607:         baij_a = HD[idx];
608:         if (roworiented) {
609:           /*value = v + i*(stepval+bs)*bs + j*bs;*/
610:           /* value = v + (i*(stepval+bs)+j)*bs; */
611:           value = v_t;
612:           v_t += bs;
613:           if (addv == ADD_VALUES) {
614:             for (ii = 0; ii < bs; ii++, value += stepval) {
615:               for (jj = ii; jj < bs2; jj += bs) baij_a[jj] += *value++;
616:             }
617:           } else {
618:             for (ii = 0; ii < bs; ii++, value += stepval) {
619:               for (jj = ii; jj < bs2; jj += bs) baij_a[jj] = *value++;
620:             }
621:           }
622:         } else {
623:           value = v + j * (stepval + bs) * bs + i * bs;
624:           if (addv == ADD_VALUES) {
625:             for (ii = 0; ii < bs; ii++, value += stepval, baij_a += bs) {
626:               for (jj = 0; jj < bs; jj++) baij_a[jj] += *value++;
627:             }
628:           } else {
629:             for (ii = 0; ii < bs; ii++, value += stepval, baij_a += bs) {
630:               for (jj = 0; jj < bs; jj++) baij_a[jj] = *value++;
631:             }
632:           }
633:         }
634:       }
635:     } else {
636:       if (!baij->donotstash) {
637:         if (roworiented) {
638:           PetscCall(MatStashValuesRowBlocked_Private(&mat->bstash, im[i], n, in, v, m, n, i));
639:         } else {
640:           PetscCall(MatStashValuesColBlocked_Private(&mat->bstash, im[i], n, in, v, m, n, i));
641:         }
642:       }
643:     }
644:   }
645:   if (PetscDefined(USE_DEBUG)) {
646:     baij->ht_total_ct += total_ct;
647:     baij->ht_insert_ct += insert_ct;
648:   }
649:   PetscFunctionReturn(PETSC_SUCCESS);
650: }

652: static PetscErrorCode MatGetValues_MPIBAIJ(Mat mat, PetscInt m, const PetscInt idxm[], PetscInt n, const PetscInt idxn[], PetscScalar v[])
653: {
654:   Mat_MPIBAIJ *baij = (Mat_MPIBAIJ *)mat->data;
655:   PetscInt     bs = mat->rmap->bs, i, j, bsrstart = mat->rmap->rstart, bsrend = mat->rmap->rend;
656:   PetscInt     bscstart = mat->cmap->rstart, bscend = mat->cmap->rend, row, col, data;

658:   PetscFunctionBegin;
659:   for (i = 0; i < m; i++) {
660:     if (idxm[i] < 0) continue; /* negative row */
661:     PetscCheck(idxm[i] < mat->rmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, idxm[i], mat->rmap->N - 1);
662:     if (idxm[i] >= bsrstart && idxm[i] < bsrend) {
663:       row = idxm[i] - bsrstart;
664:       for (j = 0; j < n; j++) {
665:         if (idxn[j] < 0) continue; /* negative column */
666:         PetscCheck(idxn[j] < mat->cmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Column too large: col %" PetscInt_FMT " max %" PetscInt_FMT, idxn[j], mat->cmap->N - 1);
667:         if (idxn[j] >= bscstart && idxn[j] < bscend) {
668:           col = idxn[j] - bscstart;
669:           PetscCall(MatGetValues_SeqBAIJ(baij->A, 1, &row, 1, &col, v + i * n + j));
670:         } else {
671:           if (!baij->colmap) PetscCall(MatCreateColmap_MPIBAIJ_Private(mat));
672: #if defined(PETSC_USE_CTABLE)
673:           PetscCall(PetscHMapIGetWithDefault(baij->colmap, idxn[j] / bs + 1, 0, &data));
674:           data--;
675: #else
676:           data = baij->colmap[idxn[j] / bs] - 1;
677: #endif
678:           if ((data < 0) || (baij->garray[data / bs] != idxn[j] / bs)) *(v + i * n + j) = 0.0;
679:           else {
680:             col = data + idxn[j] % bs;
681:             PetscCall(MatGetValues_SeqBAIJ(baij->B, 1, &row, 1, &col, v + i * n + j));
682:           }
683:         }
684:       }
685:     } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only local values currently supported");
686:   }
687:   PetscFunctionReturn(PETSC_SUCCESS);
688: }

690: static PetscErrorCode MatNorm_MPIBAIJ(Mat mat, NormType type, PetscReal *nrm)
691: {
692:   Mat_MPIBAIJ *baij = (Mat_MPIBAIJ *)mat->data;
693:   Mat_SeqBAIJ *amat = (Mat_SeqBAIJ *)baij->A->data, *bmat = (Mat_SeqBAIJ *)baij->B->data;
694:   PetscInt     i, j, bs2 = baij->bs2, bs = baij->A->rmap->bs, nz, row, col;
695:   PetscReal    sum = 0.0;
696:   MatScalar   *v;

698:   PetscFunctionBegin;
699:   if (baij->size == 1) {
700:     PetscCall(MatNorm(baij->A, type, nrm));
701:   } else {
702:     if (type == NORM_FROBENIUS) {
703:       v  = amat->a;
704:       nz = amat->nz * bs2;
705:       for (i = 0; i < nz; i++) {
706:         sum += PetscRealPart(PetscConj(*v) * (*v));
707:         v++;
708:       }
709:       v  = bmat->a;
710:       nz = bmat->nz * bs2;
711:       for (i = 0; i < nz; i++) {
712:         sum += PetscRealPart(PetscConj(*v) * (*v));
713:         v++;
714:       }
715:       PetscCallMPI(MPIU_Allreduce(&sum, nrm, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)mat)));
716:       *nrm = PetscSqrtReal(*nrm);
717:     } else if (type == NORM_1) { /* max column sum */
718:       PetscReal *tmp;
719:       PetscInt  *jj, *garray = baij->garray, cstart = baij->rstartbs;

721:       PetscCall(PetscCalloc1(mat->cmap->N, &tmp));
722:       v  = amat->a;
723:       jj = amat->j;
724:       for (i = 0; i < amat->nz; i++) {
725:         for (j = 0; j < bs; j++) {
726:           col = bs * (cstart + *jj) + j; /* column index */
727:           for (row = 0; row < bs; row++) {
728:             tmp[col] += PetscAbsScalar(*v);
729:             v++;
730:           }
731:         }
732:         jj++;
733:       }
734:       v  = bmat->a;
735:       jj = bmat->j;
736:       for (i = 0; i < bmat->nz; i++) {
737:         for (j = 0; j < bs; j++) {
738:           col = bs * garray[*jj] + j;
739:           for (row = 0; row < bs; row++) {
740:             tmp[col] += PetscAbsScalar(*v);
741:             v++;
742:           }
743:         }
744:         jj++;
745:       }
746:       PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, tmp, mat->cmap->N, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)mat)));
747:       *nrm = 0.0;
748:       for (j = 0; j < mat->cmap->N; j++) {
749:         if (tmp[j] > *nrm) *nrm = tmp[j];
750:       }
751:       PetscCall(PetscFree(tmp));
752:     } else if (type == NORM_INFINITY) { /* max row sum */
753:       PetscReal *sums;
754:       PetscCall(PetscMalloc1(bs, &sums));
755:       sum = 0.0;
756:       for (j = 0; j < amat->mbs; j++) {
757:         for (row = 0; row < bs; row++) sums[row] = 0.0;
758:         v  = amat->a + bs2 * amat->i[j];
759:         nz = amat->i[j + 1] - amat->i[j];
760:         for (i = 0; i < nz; i++) {
761:           for (col = 0; col < bs; col++) {
762:             for (row = 0; row < bs; row++) {
763:               sums[row] += PetscAbsScalar(*v);
764:               v++;
765:             }
766:           }
767:         }
768:         v  = bmat->a + bs2 * bmat->i[j];
769:         nz = bmat->i[j + 1] - bmat->i[j];
770:         for (i = 0; i < nz; i++) {
771:           for (col = 0; col < bs; col++) {
772:             for (row = 0; row < bs; row++) {
773:               sums[row] += PetscAbsScalar(*v);
774:               v++;
775:             }
776:           }
777:         }
778:         for (row = 0; row < bs; row++) {
779:           if (sums[row] > sum) sum = sums[row];
780:         }
781:       }
782:       PetscCallMPI(MPIU_Allreduce(&sum, nrm, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)mat)));
783:       PetscCall(PetscFree(sums));
784:     } else SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "No support for this norm yet");
785:   }
786:   PetscFunctionReturn(PETSC_SUCCESS);
787: }

789: /*
790:   Creates the hash table, and sets the table
791:   This table is created only once.
792:   If new entried need to be added to the matrix
793:   then the hash table has to be destroyed and
794:   recreated.
795: */
796: static PetscErrorCode MatCreateHashTable_MPIBAIJ_Private(Mat mat, PetscReal factor)
797: {
798:   Mat_MPIBAIJ *baij = (Mat_MPIBAIJ *)mat->data;
799:   Mat          A = baij->A, B = baij->B;
800:   Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data, *b = (Mat_SeqBAIJ *)B->data;
801:   PetscInt     i, j, k, nz = a->nz + b->nz, h1, *ai = a->i, *aj = a->j, *bi = b->i, *bj = b->j;
802:   PetscInt     ht_size, bs2 = baij->bs2, rstart = baij->rstartbs;
803:   PetscInt     cstart = baij->cstartbs, *garray = baij->garray, row, col, Nbs = baij->Nbs;
804:   PetscInt    *HT, key;
805:   MatScalar  **HD;
806:   PetscReal    tmp;
807: #if defined(PETSC_USE_INFO)
808:   PetscInt ct = 0, max = 0;
809: #endif

811:   PetscFunctionBegin;
812:   if (baij->ht) PetscFunctionReturn(PETSC_SUCCESS);

814:   baij->ht_size = (PetscInt)(factor * nz);
815:   ht_size       = baij->ht_size;

817:   /* Allocate Memory for Hash Table */
818:   PetscCall(PetscCalloc2(ht_size, &baij->hd, ht_size, &baij->ht));
819:   HD = baij->hd;
820:   HT = baij->ht;

822:   /* Loop Over A */
823:   for (i = 0; i < a->mbs; i++) {
824:     for (j = ai[i]; j < ai[i + 1]; j++) {
825:       row = i + rstart;
826:       col = aj[j] + cstart;

828:       key = row * Nbs + col + 1;
829:       h1  = HASH(ht_size, key, tmp);
830:       for (k = 0; k < ht_size; k++) {
831:         if (!HT[(h1 + k) % ht_size]) {
832:           HT[(h1 + k) % ht_size] = key;
833:           HD[(h1 + k) % ht_size] = a->a + j * bs2;
834:           break;
835: #if defined(PETSC_USE_INFO)
836:         } else {
837:           ct++;
838: #endif
839:         }
840:       }
841: #if defined(PETSC_USE_INFO)
842:       if (k > max) max = k;
843: #endif
844:     }
845:   }
846:   /* Loop Over B */
847:   for (i = 0; i < b->mbs; i++) {
848:     for (j = bi[i]; j < bi[i + 1]; j++) {
849:       row = i + rstart;
850:       col = garray[bj[j]];
851:       key = row * Nbs + col + 1;
852:       h1  = HASH(ht_size, key, tmp);
853:       for (k = 0; k < ht_size; k++) {
854:         if (!HT[(h1 + k) % ht_size]) {
855:           HT[(h1 + k) % ht_size] = key;
856:           HD[(h1 + k) % ht_size] = b->a + j * bs2;
857:           break;
858: #if defined(PETSC_USE_INFO)
859:         } else {
860:           ct++;
861: #endif
862:         }
863:       }
864: #if defined(PETSC_USE_INFO)
865:       if (k > max) max = k;
866: #endif
867:     }
868:   }

870:   /* Print Summary */
871: #if defined(PETSC_USE_INFO)
872:   for (i = 0, j = 0; i < ht_size; i++) {
873:     if (HT[i]) j++;
874:   }
875:   PetscCall(PetscInfo(mat, "Average Search = %5.2g,max search = %" PetscInt_FMT "\n", (!j) ? 0.0 : (double)(((PetscReal)(ct + j)) / j), max));
876: #endif
877:   PetscFunctionReturn(PETSC_SUCCESS);
878: }

880: static PetscErrorCode MatAssemblyBegin_MPIBAIJ(Mat mat, MatAssemblyType mode)
881: {
882:   Mat_MPIBAIJ *baij = (Mat_MPIBAIJ *)mat->data;
883:   PetscInt     nstash, reallocs;

885:   PetscFunctionBegin;
886:   if (baij->donotstash || mat->nooffprocentries) PetscFunctionReturn(PETSC_SUCCESS);

888:   PetscCall(MatStashScatterBegin_Private(mat, &mat->stash, mat->rmap->range));
889:   PetscCall(MatStashScatterBegin_Private(mat, &mat->bstash, baij->rangebs));
890:   PetscCall(MatStashGetInfo_Private(&mat->stash, &nstash, &reallocs));
891:   PetscCall(PetscInfo(mat, "Stash has %" PetscInt_FMT " entries,uses %" PetscInt_FMT " mallocs.\n", nstash, reallocs));
892:   PetscCall(MatStashGetInfo_Private(&mat->bstash, &nstash, &reallocs));
893:   PetscCall(PetscInfo(mat, "Block-Stash has %" PetscInt_FMT " entries, uses %" PetscInt_FMT " mallocs.\n", nstash, reallocs));
894:   PetscFunctionReturn(PETSC_SUCCESS);
895: }

897: static PetscErrorCode MatAssemblyEnd_MPIBAIJ(Mat mat, MatAssemblyType mode)
898: {
899:   Mat_MPIBAIJ *baij = (Mat_MPIBAIJ *)mat->data;
900:   Mat_SeqBAIJ *a    = (Mat_SeqBAIJ *)baij->A->data;
901:   PetscInt     i, j, rstart, ncols, flg, bs2 = baij->bs2;
902:   PetscInt    *row, *col;
903:   PetscBool    r1, r2, r3, other_disassembled;
904:   MatScalar   *val;
905:   PetscMPIInt  n;

907:   PetscFunctionBegin;
908:   /* do not use 'b=(Mat_SeqBAIJ*)baij->B->data' as B can be reset in disassembly */
909:   if (!baij->donotstash && !mat->nooffprocentries) {
910:     while (1) {
911:       PetscCall(MatStashScatterGetMesg_Private(&mat->stash, &n, &row, &col, &val, &flg));
912:       if (!flg) break;

914:       for (i = 0; i < n;) {
915:         /* Now identify the consecutive vals belonging to the same row */
916:         for (j = i, rstart = row[j]; j < n; j++) {
917:           if (row[j] != rstart) break;
918:         }
919:         if (j < n) ncols = j - i;
920:         else ncols = n - i;
921:         /* Now assemble all these values with a single function call */
922:         PetscCall(MatSetValues_MPIBAIJ(mat, 1, row + i, ncols, col + i, val + i, mat->insertmode));
923:         i = j;
924:       }
925:     }
926:     PetscCall(MatStashScatterEnd_Private(&mat->stash));
927:     /* Now process the block-stash. Since the values are stashed column-oriented,
928:        set the row-oriented flag to column-oriented, and after MatSetValues()
929:        restore the original flags */
930:     r1 = baij->roworiented;
931:     r2 = a->roworiented;
932:     r3 = ((Mat_SeqBAIJ *)baij->B->data)->roworiented;

934:     baij->roworiented                           = PETSC_FALSE;
935:     a->roworiented                              = PETSC_FALSE;
936:     ((Mat_SeqBAIJ *)baij->B->data)->roworiented = PETSC_FALSE;
937:     while (1) {
938:       PetscCall(MatStashScatterGetMesg_Private(&mat->bstash, &n, &row, &col, &val, &flg));
939:       if (!flg) break;

941:       for (i = 0; i < n;) {
942:         /* Now identify the consecutive vals belonging to the same row */
943:         for (j = i, rstart = row[j]; j < n; j++) {
944:           if (row[j] != rstart) break;
945:         }
946:         if (j < n) ncols = j - i;
947:         else ncols = n - i;
948:         PetscCall(MatSetValuesBlocked_MPIBAIJ(mat, 1, row + i, ncols, col + i, val + i * bs2, mat->insertmode));
949:         i = j;
950:       }
951:     }
952:     PetscCall(MatStashScatterEnd_Private(&mat->bstash));

954:     baij->roworiented                           = r1;
955:     a->roworiented                              = r2;
956:     ((Mat_SeqBAIJ *)baij->B->data)->roworiented = r3;
957:   }

959:   PetscCall(MatAssemblyBegin(baij->A, mode));
960:   PetscCall(MatAssemblyEnd(baij->A, mode));

962:   /* determine if any processor has disassembled, if so we must
963:      also disassemble ourselves, in order that we may reassemble. */
964:   /*
965:      if nonzero structure of submatrix B cannot change then we know that
966:      no processor disassembled thus we can skip this stuff
967:   */
968:   if (!((Mat_SeqBAIJ *)baij->B->data)->nonew) {
969:     PetscCallMPI(MPIU_Allreduce(&mat->was_assembled, &other_disassembled, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)mat)));
970:     if (mat->was_assembled && !other_disassembled) PetscCall(MatDisAssemble_MPIBAIJ(mat));
971:   }

973:   if (!mat->was_assembled && mode == MAT_FINAL_ASSEMBLY) PetscCall(MatSetUpMultiply_MPIBAIJ(mat));
974:   PetscCall(MatAssemblyBegin(baij->B, mode));
975:   PetscCall(MatAssemblyEnd(baij->B, mode));

977: #if defined(PETSC_USE_INFO)
978:   if (baij->ht && mode == MAT_FINAL_ASSEMBLY) {
979:     PetscCall(PetscInfo(mat, "Average Hash Table Search in MatSetValues = %5.2f\n", (double)((PetscReal)baij->ht_total_ct) / baij->ht_insert_ct));

981:     baij->ht_total_ct  = 0;
982:     baij->ht_insert_ct = 0;
983:   }
984: #endif
985:   if (baij->ht_flag && !baij->ht && mode == MAT_FINAL_ASSEMBLY) {
986:     PetscCall(MatCreateHashTable_MPIBAIJ_Private(mat, baij->ht_fact));

988:     mat->ops->setvalues        = MatSetValues_MPIBAIJ_HT;
989:     mat->ops->setvaluesblocked = MatSetValuesBlocked_MPIBAIJ_HT;
990:   }

992:   PetscCall(PetscFree2(baij->rowvalues, baij->rowindices));

994:   baij->rowvalues = NULL;

996:   /* if no new nonzero locations are allowed in matrix then only set the matrix state the first time through */
997:   if ((!mat->was_assembled && mode == MAT_FINAL_ASSEMBLY) || !((Mat_SeqBAIJ *)baij->A->data)->nonew) {
998:     PetscObjectState state = baij->A->nonzerostate + baij->B->nonzerostate;
999:     PetscCallMPI(MPIU_Allreduce(&state, &mat->nonzerostate, 1, MPIU_INT64, MPI_SUM, PetscObjectComm((PetscObject)mat)));
1000:   }
1001:   PetscFunctionReturn(PETSC_SUCCESS);
1002: }

1004: extern PetscErrorCode MatView_SeqBAIJ(Mat, PetscViewer);
1005: #include <petscdraw.h>
1006: static PetscErrorCode MatView_MPIBAIJ_ASCIIorDraworSocket(Mat mat, PetscViewer viewer)
1007: {
1008:   Mat_MPIBAIJ      *baij = (Mat_MPIBAIJ *)mat->data;
1009:   PetscMPIInt       rank = baij->rank;
1010:   PetscInt          bs   = mat->rmap->bs;
1011:   PetscBool         iascii, isdraw;
1012:   PetscViewer       sviewer;
1013:   PetscViewerFormat format;

1015:   PetscFunctionBegin;
1016:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
1017:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
1018:   if (iascii) {
1019:     PetscCall(PetscViewerGetFormat(viewer, &format));
1020:     if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1021:       MatInfo info;
1022:       PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)mat), &rank));
1023:       PetscCall(MatGetInfo(mat, MAT_LOCAL, &info));
1024:       PetscCall(PetscViewerASCIIPushSynchronized(viewer));
1025:       PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Local rows %" PetscInt_FMT " nz %" PetscInt_FMT " nz alloced %" PetscInt_FMT " bs %" PetscInt_FMT " mem %g\n", rank, mat->rmap->n, (PetscInt)info.nz_used, (PetscInt)info.nz_allocated,
1026:                                                    mat->rmap->bs, info.memory));
1027:       PetscCall(MatGetInfo(baij->A, MAT_LOCAL, &info));
1028:       PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] on-diagonal part: nz %" PetscInt_FMT " \n", rank, (PetscInt)info.nz_used));
1029:       PetscCall(MatGetInfo(baij->B, MAT_LOCAL, &info));
1030:       PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] off-diagonal part: nz %" PetscInt_FMT " \n", rank, (PetscInt)info.nz_used));
1031:       PetscCall(PetscViewerFlush(viewer));
1032:       PetscCall(PetscViewerASCIIPopSynchronized(viewer));
1033:       PetscCall(PetscViewerASCIIPrintf(viewer, "Information on VecScatter used in matrix-vector product: \n"));
1034:       PetscCall(VecScatterView(baij->Mvctx, viewer));
1035:       PetscFunctionReturn(PETSC_SUCCESS);
1036:     } else if (format == PETSC_VIEWER_ASCII_INFO) {
1037:       PetscCall(PetscViewerASCIIPrintf(viewer, "  block size is %" PetscInt_FMT "\n", bs));
1038:       PetscFunctionReturn(PETSC_SUCCESS);
1039:     } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) {
1040:       PetscFunctionReturn(PETSC_SUCCESS);
1041:     }
1042:   }

1044:   if (isdraw) {
1045:     PetscDraw draw;
1046:     PetscBool isnull;
1047:     PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
1048:     PetscCall(PetscDrawIsNull(draw, &isnull));
1049:     if (isnull) PetscFunctionReturn(PETSC_SUCCESS);
1050:   }

1052:   {
1053:     /* assemble the entire matrix onto first processor. */
1054:     Mat          A;
1055:     Mat_SeqBAIJ *Aloc;
1056:     PetscInt     M = mat->rmap->N, N = mat->cmap->N, *ai, *aj, col, i, j, k, *rvals, mbs = baij->mbs;
1057:     MatScalar   *a;
1058:     const char  *matname;

1060:     /* Here we are creating a temporary matrix, so will assume MPIBAIJ is acceptable */
1061:     /* Perhaps this should be the type of mat? */
1062:     PetscCall(MatCreate(PetscObjectComm((PetscObject)mat), &A));
1063:     if (rank == 0) {
1064:       PetscCall(MatSetSizes(A, M, N, M, N));
1065:     } else {
1066:       PetscCall(MatSetSizes(A, 0, 0, M, N));
1067:     }
1068:     PetscCall(MatSetType(A, MATMPIBAIJ));
1069:     PetscCall(MatMPIBAIJSetPreallocation(A, mat->rmap->bs, 0, NULL, 0, NULL));
1070:     PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_FALSE));

1072:     /* copy over the A part */
1073:     Aloc = (Mat_SeqBAIJ *)baij->A->data;
1074:     ai   = Aloc->i;
1075:     aj   = Aloc->j;
1076:     a    = Aloc->a;
1077:     PetscCall(PetscMalloc1(bs, &rvals));

1079:     for (i = 0; i < mbs; i++) {
1080:       rvals[0] = bs * (baij->rstartbs + i);
1081:       for (j = 1; j < bs; j++) rvals[j] = rvals[j - 1] + 1;
1082:       for (j = ai[i]; j < ai[i + 1]; j++) {
1083:         col = (baij->cstartbs + aj[j]) * bs;
1084:         for (k = 0; k < bs; k++) {
1085:           PetscCall(MatSetValues_MPIBAIJ(A, bs, rvals, 1, &col, a, INSERT_VALUES));
1086:           col++;
1087:           a += bs;
1088:         }
1089:       }
1090:     }
1091:     /* copy over the B part */
1092:     Aloc = (Mat_SeqBAIJ *)baij->B->data;
1093:     ai   = Aloc->i;
1094:     aj   = Aloc->j;
1095:     a    = Aloc->a;
1096:     for (i = 0; i < mbs; i++) {
1097:       rvals[0] = bs * (baij->rstartbs + i);
1098:       for (j = 1; j < bs; j++) rvals[j] = rvals[j - 1] + 1;
1099:       for (j = ai[i]; j < ai[i + 1]; j++) {
1100:         col = baij->garray[aj[j]] * bs;
1101:         for (k = 0; k < bs; k++) {
1102:           PetscCall(MatSetValues_MPIBAIJ(A, bs, rvals, 1, &col, a, INSERT_VALUES));
1103:           col++;
1104:           a += bs;
1105:         }
1106:       }
1107:     }
1108:     PetscCall(PetscFree(rvals));
1109:     PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
1110:     PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
1111:     /*
1112:        Everyone has to call to draw the matrix since the graphics waits are
1113:        synchronized across all processors that share the PetscDraw object
1114:     */
1115:     PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
1116:     if (((PetscObject)mat)->name) PetscCall(PetscObjectGetName((PetscObject)mat, &matname));
1117:     if (rank == 0) {
1118:       if (((PetscObject)mat)->name) PetscCall(PetscObjectSetName((PetscObject)((Mat_MPIBAIJ *)A->data)->A, matname));
1119:       PetscCall(MatView_SeqBAIJ(((Mat_MPIBAIJ *)A->data)->A, sviewer));
1120:     }
1121:     PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
1122:     PetscCall(MatDestroy(&A));
1123:   }
1124:   PetscFunctionReturn(PETSC_SUCCESS);
1125: }

1127: /* Used for both MPIBAIJ and MPISBAIJ matrices */
1128: PetscErrorCode MatView_MPIBAIJ_Binary(Mat mat, PetscViewer viewer)
1129: {
1130:   Mat_MPIBAIJ    *aij    = (Mat_MPIBAIJ *)mat->data;
1131:   Mat_SeqBAIJ    *A      = (Mat_SeqBAIJ *)aij->A->data;
1132:   Mat_SeqBAIJ    *B      = (Mat_SeqBAIJ *)aij->B->data;
1133:   const PetscInt *garray = aij->garray;
1134:   PetscInt        header[4], M, N, m, rs, cs, bs, cnt, i, j, ja, jb, k, l;
1135:   PetscCount      nz, hnz;
1136:   PetscInt       *rowlens, *colidxs;
1137:   PetscScalar    *matvals;
1138:   PetscMPIInt     rank;

1140:   PetscFunctionBegin;
1141:   PetscCall(PetscViewerSetUp(viewer));

1143:   M  = mat->rmap->N;
1144:   N  = mat->cmap->N;
1145:   m  = mat->rmap->n;
1146:   rs = mat->rmap->rstart;
1147:   cs = mat->cmap->rstart;
1148:   bs = mat->rmap->bs;
1149:   nz = bs * bs * (A->nz + B->nz);

1151:   /* write matrix header */
1152:   header[0] = MAT_FILE_CLASSID;
1153:   header[1] = M;
1154:   header[2] = N;
1155:   PetscCallMPI(MPI_Reduce(&nz, &hnz, 1, MPIU_COUNT, MPI_SUM, 0, PetscObjectComm((PetscObject)mat)));
1156:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)mat), &rank));
1157:   if (rank == 0) PetscCall(PetscIntCast(hnz, &header[3]));
1158:   PetscCall(PetscViewerBinaryWrite(viewer, header, 4, PETSC_INT));

1160:   /* fill in and store row lengths */
1161:   PetscCall(PetscMalloc1(m, &rowlens));
1162:   for (cnt = 0, i = 0; i < A->mbs; i++)
1163:     for (j = 0; j < bs; j++) rowlens[cnt++] = bs * (A->i[i + 1] - A->i[i] + B->i[i + 1] - B->i[i]);
1164:   PetscCall(PetscViewerBinaryWriteAll(viewer, rowlens, m, rs, M, PETSC_INT));
1165:   PetscCall(PetscFree(rowlens));

1167:   /* fill in and store column indices */
1168:   PetscCall(PetscMalloc1(nz, &colidxs));
1169:   for (cnt = 0, i = 0; i < A->mbs; i++) {
1170:     for (k = 0; k < bs; k++) {
1171:       for (jb = B->i[i]; jb < B->i[i + 1]; jb++) {
1172:         if (garray[B->j[jb]] > cs / bs) break;
1173:         for (l = 0; l < bs; l++) colidxs[cnt++] = bs * garray[B->j[jb]] + l;
1174:       }
1175:       for (ja = A->i[i]; ja < A->i[i + 1]; ja++)
1176:         for (l = 0; l < bs; l++) colidxs[cnt++] = bs * A->j[ja] + l + cs;
1177:       for (; jb < B->i[i + 1]; jb++)
1178:         for (l = 0; l < bs; l++) colidxs[cnt++] = bs * garray[B->j[jb]] + l;
1179:     }
1180:   }
1181:   PetscCheck(cnt == nz, PETSC_COMM_SELF, PETSC_ERR_LIB, "Internal PETSc error: cnt = %" PetscInt_FMT " nz = %" PetscCount_FMT, cnt, nz);
1182:   PetscCall(PetscViewerBinaryWriteAll(viewer, colidxs, nz, PETSC_DECIDE, PETSC_DECIDE, PETSC_INT));
1183:   PetscCall(PetscFree(colidxs));

1185:   /* fill in and store nonzero values */
1186:   PetscCall(PetscMalloc1(nz, &matvals));
1187:   for (cnt = 0, i = 0; i < A->mbs; i++) {
1188:     for (k = 0; k < bs; k++) {
1189:       for (jb = B->i[i]; jb < B->i[i + 1]; jb++) {
1190:         if (garray[B->j[jb]] > cs / bs) break;
1191:         for (l = 0; l < bs; l++) matvals[cnt++] = B->a[bs * (bs * jb + l) + k];
1192:       }
1193:       for (ja = A->i[i]; ja < A->i[i + 1]; ja++)
1194:         for (l = 0; l < bs; l++) matvals[cnt++] = A->a[bs * (bs * ja + l) + k];
1195:       for (; jb < B->i[i + 1]; jb++)
1196:         for (l = 0; l < bs; l++) matvals[cnt++] = B->a[bs * (bs * jb + l) + k];
1197:     }
1198:   }
1199:   PetscCall(PetscViewerBinaryWriteAll(viewer, matvals, nz, PETSC_DECIDE, PETSC_DECIDE, PETSC_SCALAR));
1200:   PetscCall(PetscFree(matvals));

1202:   /* write block size option to the viewer's .info file */
1203:   PetscCall(MatView_Binary_BlockSizes(mat, viewer));
1204:   PetscFunctionReturn(PETSC_SUCCESS);
1205: }

1207: PetscErrorCode MatView_MPIBAIJ(Mat mat, PetscViewer viewer)
1208: {
1209:   PetscBool iascii, isdraw, issocket, isbinary;

1211:   PetscFunctionBegin;
1212:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
1213:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
1214:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSOCKET, &issocket));
1215:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1216:   if (iascii || isdraw || issocket) {
1217:     PetscCall(MatView_MPIBAIJ_ASCIIorDraworSocket(mat, viewer));
1218:   } else if (isbinary) PetscCall(MatView_MPIBAIJ_Binary(mat, viewer));
1219:   PetscFunctionReturn(PETSC_SUCCESS);
1220: }

1222: static PetscErrorCode MatMult_MPIBAIJ(Mat A, Vec xx, Vec yy)
1223: {
1224:   Mat_MPIBAIJ *a = (Mat_MPIBAIJ *)A->data;
1225:   PetscInt     nt;

1227:   PetscFunctionBegin;
1228:   PetscCall(VecGetLocalSize(xx, &nt));
1229:   PetscCheck(nt == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Incompatible partition of A and xx");
1230:   PetscCall(VecGetLocalSize(yy, &nt));
1231:   PetscCheck(nt == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Incompatible partition of A and yy");
1232:   PetscCall(VecScatterBegin(a->Mvctx, xx, a->lvec, INSERT_VALUES, SCATTER_FORWARD));
1233:   PetscCall((*a->A->ops->mult)(a->A, xx, yy));
1234:   PetscCall(VecScatterEnd(a->Mvctx, xx, a->lvec, INSERT_VALUES, SCATTER_FORWARD));
1235:   PetscCall((*a->B->ops->multadd)(a->B, a->lvec, yy, yy));
1236:   PetscFunctionReturn(PETSC_SUCCESS);
1237: }

1239: static PetscErrorCode MatMultAdd_MPIBAIJ(Mat A, Vec xx, Vec yy, Vec zz)
1240: {
1241:   Mat_MPIBAIJ *a = (Mat_MPIBAIJ *)A->data;

1243:   PetscFunctionBegin;
1244:   PetscCall(VecScatterBegin(a->Mvctx, xx, a->lvec, INSERT_VALUES, SCATTER_FORWARD));
1245:   PetscCall((*a->A->ops->multadd)(a->A, xx, yy, zz));
1246:   PetscCall(VecScatterEnd(a->Mvctx, xx, a->lvec, INSERT_VALUES, SCATTER_FORWARD));
1247:   PetscCall((*a->B->ops->multadd)(a->B, a->lvec, zz, zz));
1248:   PetscFunctionReturn(PETSC_SUCCESS);
1249: }

1251: static PetscErrorCode MatMultTranspose_MPIBAIJ(Mat A, Vec xx, Vec yy)
1252: {
1253:   Mat_MPIBAIJ *a = (Mat_MPIBAIJ *)A->data;

1255:   PetscFunctionBegin;
1256:   /* do nondiagonal part */
1257:   PetscCall((*a->B->ops->multtranspose)(a->B, xx, a->lvec));
1258:   /* do local part */
1259:   PetscCall((*a->A->ops->multtranspose)(a->A, xx, yy));
1260:   /* add partial results together */
1261:   PetscCall(VecScatterBegin(a->Mvctx, a->lvec, yy, ADD_VALUES, SCATTER_REVERSE));
1262:   PetscCall(VecScatterEnd(a->Mvctx, a->lvec, yy, ADD_VALUES, SCATTER_REVERSE));
1263:   PetscFunctionReturn(PETSC_SUCCESS);
1264: }

1266: static PetscErrorCode MatMultTransposeAdd_MPIBAIJ(Mat A, Vec xx, Vec yy, Vec zz)
1267: {
1268:   Mat_MPIBAIJ *a = (Mat_MPIBAIJ *)A->data;

1270:   PetscFunctionBegin;
1271:   /* do nondiagonal part */
1272:   PetscCall((*a->B->ops->multtranspose)(a->B, xx, a->lvec));
1273:   /* do local part */
1274:   PetscCall((*a->A->ops->multtransposeadd)(a->A, xx, yy, zz));
1275:   /* add partial results together */
1276:   PetscCall(VecScatterBegin(a->Mvctx, a->lvec, zz, ADD_VALUES, SCATTER_REVERSE));
1277:   PetscCall(VecScatterEnd(a->Mvctx, a->lvec, zz, ADD_VALUES, SCATTER_REVERSE));
1278:   PetscFunctionReturn(PETSC_SUCCESS);
1279: }

1281: /*
1282:   This only works correctly for square matrices where the subblock A->A is the
1283:    diagonal block
1284: */
1285: static PetscErrorCode MatGetDiagonal_MPIBAIJ(Mat A, Vec v)
1286: {
1287:   PetscFunctionBegin;
1288:   PetscCheck(A->rmap->N == A->cmap->N, PETSC_COMM_SELF, PETSC_ERR_SUP, "Supports only square matrix where A->A is diag block");
1289:   PetscCall(MatGetDiagonal(((Mat_MPIBAIJ *)A->data)->A, v));
1290:   PetscFunctionReturn(PETSC_SUCCESS);
1291: }

1293: static PetscErrorCode MatScale_MPIBAIJ(Mat A, PetscScalar aa)
1294: {
1295:   Mat_MPIBAIJ *a = (Mat_MPIBAIJ *)A->data;

1297:   PetscFunctionBegin;
1298:   PetscCall(MatScale(a->A, aa));
1299:   PetscCall(MatScale(a->B, aa));
1300:   PetscFunctionReturn(PETSC_SUCCESS);
1301: }

1303: PetscErrorCode MatGetRow_MPIBAIJ(Mat matin, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v)
1304: {
1305:   Mat_MPIBAIJ *mat = (Mat_MPIBAIJ *)matin->data;
1306:   PetscScalar *vworkA, *vworkB, **pvA, **pvB, *v_p;
1307:   PetscInt     bs = matin->rmap->bs, bs2 = mat->bs2, i, *cworkA, *cworkB, **pcA, **pcB;
1308:   PetscInt     nztot, nzA, nzB, lrow, brstart = matin->rmap->rstart, brend = matin->rmap->rend;
1309:   PetscInt    *cmap, *idx_p, cstart = mat->cstartbs;

1311:   PetscFunctionBegin;
1312:   PetscCheck(row >= brstart && row < brend, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only local rows");
1313:   PetscCheck(!mat->getrowactive, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Already active");
1314:   mat->getrowactive = PETSC_TRUE;

1316:   if (!mat->rowvalues && (idx || v)) {
1317:     /*
1318:         allocate enough space to hold information from the longest row.
1319:     */
1320:     Mat_SeqBAIJ *Aa = (Mat_SeqBAIJ *)mat->A->data, *Ba = (Mat_SeqBAIJ *)mat->B->data;
1321:     PetscInt     max = 1, mbs = mat->mbs, tmp;
1322:     for (i = 0; i < mbs; i++) {
1323:       tmp = Aa->i[i + 1] - Aa->i[i] + Ba->i[i + 1] - Ba->i[i];
1324:       if (max < tmp) max = tmp;
1325:     }
1326:     PetscCall(PetscMalloc2(max * bs2, &mat->rowvalues, max * bs2, &mat->rowindices));
1327:   }
1328:   lrow = row - brstart;

1330:   pvA = &vworkA;
1331:   pcA = &cworkA;
1332:   pvB = &vworkB;
1333:   pcB = &cworkB;
1334:   if (!v) {
1335:     pvA = NULL;
1336:     pvB = NULL;
1337:   }
1338:   if (!idx) {
1339:     pcA = NULL;
1340:     if (!v) pcB = NULL;
1341:   }
1342:   PetscCall((*mat->A->ops->getrow)(mat->A, lrow, &nzA, pcA, pvA));
1343:   PetscCall((*mat->B->ops->getrow)(mat->B, lrow, &nzB, pcB, pvB));
1344:   nztot = nzA + nzB;

1346:   cmap = mat->garray;
1347:   if (v || idx) {
1348:     if (nztot) {
1349:       /* Sort by increasing column numbers, assuming A and B already sorted */
1350:       PetscInt imark = -1;
1351:       if (v) {
1352:         *v = v_p = mat->rowvalues;
1353:         for (i = 0; i < nzB; i++) {
1354:           if (cmap[cworkB[i] / bs] < cstart) v_p[i] = vworkB[i];
1355:           else break;
1356:         }
1357:         imark = i;
1358:         for (i = 0; i < nzA; i++) v_p[imark + i] = vworkA[i];
1359:         for (i = imark; i < nzB; i++) v_p[nzA + i] = vworkB[i];
1360:       }
1361:       if (idx) {
1362:         *idx = idx_p = mat->rowindices;
1363:         if (imark > -1) {
1364:           for (i = 0; i < imark; i++) idx_p[i] = cmap[cworkB[i] / bs] * bs + cworkB[i] % bs;
1365:         } else {
1366:           for (i = 0; i < nzB; i++) {
1367:             if (cmap[cworkB[i] / bs] < cstart) idx_p[i] = cmap[cworkB[i] / bs] * bs + cworkB[i] % bs;
1368:             else break;
1369:           }
1370:           imark = i;
1371:         }
1372:         for (i = 0; i < nzA; i++) idx_p[imark + i] = cstart * bs + cworkA[i];
1373:         for (i = imark; i < nzB; i++) idx_p[nzA + i] = cmap[cworkB[i] / bs] * bs + cworkB[i] % bs;
1374:       }
1375:     } else {
1376:       if (idx) *idx = NULL;
1377:       if (v) *v = NULL;
1378:     }
1379:   }
1380:   *nz = nztot;
1381:   PetscCall((*mat->A->ops->restorerow)(mat->A, lrow, &nzA, pcA, pvA));
1382:   PetscCall((*mat->B->ops->restorerow)(mat->B, lrow, &nzB, pcB, pvB));
1383:   PetscFunctionReturn(PETSC_SUCCESS);
1384: }

1386: PetscErrorCode MatRestoreRow_MPIBAIJ(Mat mat, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v)
1387: {
1388:   Mat_MPIBAIJ *baij = (Mat_MPIBAIJ *)mat->data;

1390:   PetscFunctionBegin;
1391:   PetscCheck(baij->getrowactive, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "MatGetRow not called");
1392:   baij->getrowactive = PETSC_FALSE;
1393:   PetscFunctionReturn(PETSC_SUCCESS);
1394: }

1396: static PetscErrorCode MatZeroEntries_MPIBAIJ(Mat A)
1397: {
1398:   Mat_MPIBAIJ *l = (Mat_MPIBAIJ *)A->data;

1400:   PetscFunctionBegin;
1401:   PetscCall(MatZeroEntries(l->A));
1402:   PetscCall(MatZeroEntries(l->B));
1403:   PetscFunctionReturn(PETSC_SUCCESS);
1404: }

1406: static PetscErrorCode MatGetInfo_MPIBAIJ(Mat matin, MatInfoType flag, MatInfo *info)
1407: {
1408:   Mat_MPIBAIJ   *a = (Mat_MPIBAIJ *)matin->data;
1409:   Mat            A = a->A, B = a->B;
1410:   PetscLogDouble isend[5], irecv[5];

1412:   PetscFunctionBegin;
1413:   info->block_size = (PetscReal)matin->rmap->bs;

1415:   PetscCall(MatGetInfo(A, MAT_LOCAL, info));

1417:   isend[0] = info->nz_used;
1418:   isend[1] = info->nz_allocated;
1419:   isend[2] = info->nz_unneeded;
1420:   isend[3] = info->memory;
1421:   isend[4] = info->mallocs;

1423:   PetscCall(MatGetInfo(B, MAT_LOCAL, info));

1425:   isend[0] += info->nz_used;
1426:   isend[1] += info->nz_allocated;
1427:   isend[2] += info->nz_unneeded;
1428:   isend[3] += info->memory;
1429:   isend[4] += info->mallocs;

1431:   if (flag == MAT_LOCAL) {
1432:     info->nz_used      = isend[0];
1433:     info->nz_allocated = isend[1];
1434:     info->nz_unneeded  = isend[2];
1435:     info->memory       = isend[3];
1436:     info->mallocs      = isend[4];
1437:   } else if (flag == MAT_GLOBAL_MAX) {
1438:     PetscCallMPI(MPIU_Allreduce(isend, irecv, 5, MPIU_PETSCLOGDOUBLE, MPI_MAX, PetscObjectComm((PetscObject)matin)));

1440:     info->nz_used      = irecv[0];
1441:     info->nz_allocated = irecv[1];
1442:     info->nz_unneeded  = irecv[2];
1443:     info->memory       = irecv[3];
1444:     info->mallocs      = irecv[4];
1445:   } else if (flag == MAT_GLOBAL_SUM) {
1446:     PetscCallMPI(MPIU_Allreduce(isend, irecv, 5, MPIU_PETSCLOGDOUBLE, MPI_SUM, PetscObjectComm((PetscObject)matin)));

1448:     info->nz_used      = irecv[0];
1449:     info->nz_allocated = irecv[1];
1450:     info->nz_unneeded  = irecv[2];
1451:     info->memory       = irecv[3];
1452:     info->mallocs      = irecv[4];
1453:   } else SETERRQ(PetscObjectComm((PetscObject)matin), PETSC_ERR_ARG_WRONG, "Unknown MatInfoType argument %d", (int)flag);
1454:   info->fill_ratio_given  = 0; /* no parallel LU/ILU/Cholesky */
1455:   info->fill_ratio_needed = 0;
1456:   info->factor_mallocs    = 0;
1457:   PetscFunctionReturn(PETSC_SUCCESS);
1458: }

1460: static PetscErrorCode MatSetOption_MPIBAIJ(Mat A, MatOption op, PetscBool flg)
1461: {
1462:   Mat_MPIBAIJ *a = (Mat_MPIBAIJ *)A->data;

1464:   PetscFunctionBegin;
1465:   switch (op) {
1466:   case MAT_NEW_NONZERO_LOCATIONS:
1467:   case MAT_NEW_NONZERO_ALLOCATION_ERR:
1468:   case MAT_UNUSED_NONZERO_LOCATION_ERR:
1469:   case MAT_KEEP_NONZERO_PATTERN:
1470:   case MAT_NEW_NONZERO_LOCATION_ERR:
1471:     MatCheckPreallocated(A, 1);
1472:     PetscCall(MatSetOption(a->A, op, flg));
1473:     PetscCall(MatSetOption(a->B, op, flg));
1474:     break;
1475:   case MAT_ROW_ORIENTED:
1476:     MatCheckPreallocated(A, 1);
1477:     a->roworiented = flg;

1479:     PetscCall(MatSetOption(a->A, op, flg));
1480:     PetscCall(MatSetOption(a->B, op, flg));
1481:     break;
1482:   case MAT_FORCE_DIAGONAL_ENTRIES:
1483:   case MAT_SORTED_FULL:
1484:     PetscCall(PetscInfo(A, "Option %s ignored\n", MatOptions[op]));
1485:     break;
1486:   case MAT_IGNORE_OFF_PROC_ENTRIES:
1487:     a->donotstash = flg;
1488:     break;
1489:   case MAT_USE_HASH_TABLE:
1490:     a->ht_flag = flg;
1491:     a->ht_fact = 1.39;
1492:     break;
1493:   case MAT_SYMMETRIC:
1494:   case MAT_STRUCTURALLY_SYMMETRIC:
1495:   case MAT_HERMITIAN:
1496:   case MAT_SUBMAT_SINGLEIS:
1497:   case MAT_SYMMETRY_ETERNAL:
1498:   case MAT_STRUCTURAL_SYMMETRY_ETERNAL:
1499:   case MAT_SPD_ETERNAL:
1500:     /* if the diagonal matrix is square it inherits some of the properties above */
1501:     break;
1502:   default:
1503:     SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "unknown option %d", op);
1504:   }
1505:   PetscFunctionReturn(PETSC_SUCCESS);
1506: }

1508: static PetscErrorCode MatTranspose_MPIBAIJ(Mat A, MatReuse reuse, Mat *matout)
1509: {
1510:   Mat_MPIBAIJ *baij = (Mat_MPIBAIJ *)A->data;
1511:   Mat_SeqBAIJ *Aloc;
1512:   Mat          B;
1513:   PetscInt     M = A->rmap->N, N = A->cmap->N, *ai, *aj, i, *rvals, j, k, col;
1514:   PetscInt     bs = A->rmap->bs, mbs = baij->mbs;
1515:   MatScalar   *a;

1517:   PetscFunctionBegin;
1518:   if (reuse == MAT_REUSE_MATRIX) PetscCall(MatTransposeCheckNonzeroState_Private(A, *matout));
1519:   if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) {
1520:     PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
1521:     PetscCall(MatSetSizes(B, A->cmap->n, A->rmap->n, N, M));
1522:     PetscCall(MatSetType(B, ((PetscObject)A)->type_name));
1523:     /* Do not know preallocation information, but must set block size */
1524:     PetscCall(MatMPIBAIJSetPreallocation(B, A->rmap->bs, PETSC_DECIDE, NULL, PETSC_DECIDE, NULL));
1525:   } else {
1526:     B = *matout;
1527:   }

1529:   /* copy over the A part */
1530:   Aloc = (Mat_SeqBAIJ *)baij->A->data;
1531:   ai   = Aloc->i;
1532:   aj   = Aloc->j;
1533:   a    = Aloc->a;
1534:   PetscCall(PetscMalloc1(bs, &rvals));

1536:   for (i = 0; i < mbs; i++) {
1537:     rvals[0] = bs * (baij->rstartbs + i);
1538:     for (j = 1; j < bs; j++) rvals[j] = rvals[j - 1] + 1;
1539:     for (j = ai[i]; j < ai[i + 1]; j++) {
1540:       col = (baij->cstartbs + aj[j]) * bs;
1541:       for (k = 0; k < bs; k++) {
1542:         PetscCall(MatSetValues_MPIBAIJ(B, 1, &col, bs, rvals, a, INSERT_VALUES));

1544:         col++;
1545:         a += bs;
1546:       }
1547:     }
1548:   }
1549:   /* copy over the B part */
1550:   Aloc = (Mat_SeqBAIJ *)baij->B->data;
1551:   ai   = Aloc->i;
1552:   aj   = Aloc->j;
1553:   a    = Aloc->a;
1554:   for (i = 0; i < mbs; i++) {
1555:     rvals[0] = bs * (baij->rstartbs + i);
1556:     for (j = 1; j < bs; j++) rvals[j] = rvals[j - 1] + 1;
1557:     for (j = ai[i]; j < ai[i + 1]; j++) {
1558:       col = baij->garray[aj[j]] * bs;
1559:       for (k = 0; k < bs; k++) {
1560:         PetscCall(MatSetValues_MPIBAIJ(B, 1, &col, bs, rvals, a, INSERT_VALUES));
1561:         col++;
1562:         a += bs;
1563:       }
1564:     }
1565:   }
1566:   PetscCall(PetscFree(rvals));
1567:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
1568:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));

1570:   if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) *matout = B;
1571:   else PetscCall(MatHeaderMerge(A, &B));
1572:   PetscFunctionReturn(PETSC_SUCCESS);
1573: }

1575: static PetscErrorCode MatDiagonalScale_MPIBAIJ(Mat mat, Vec ll, Vec rr)
1576: {
1577:   Mat_MPIBAIJ *baij = (Mat_MPIBAIJ *)mat->data;
1578:   Mat          a = baij->A, b = baij->B;
1579:   PetscInt     s1, s2, s3;

1581:   PetscFunctionBegin;
1582:   PetscCall(MatGetLocalSize(mat, &s2, &s3));
1583:   if (rr) {
1584:     PetscCall(VecGetLocalSize(rr, &s1));
1585:     PetscCheck(s1 == s3, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "right vector non-conforming local size");
1586:     /* Overlap communication with computation. */
1587:     PetscCall(VecScatterBegin(baij->Mvctx, rr, baij->lvec, INSERT_VALUES, SCATTER_FORWARD));
1588:   }
1589:   if (ll) {
1590:     PetscCall(VecGetLocalSize(ll, &s1));
1591:     PetscCheck(s1 == s2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "left vector non-conforming local size");
1592:     PetscUseTypeMethod(b, diagonalscale, ll, NULL);
1593:   }
1594:   /* scale  the diagonal block */
1595:   PetscUseTypeMethod(a, diagonalscale, ll, rr);

1597:   if (rr) {
1598:     /* Do a scatter end and then right scale the off-diagonal block */
1599:     PetscCall(VecScatterEnd(baij->Mvctx, rr, baij->lvec, INSERT_VALUES, SCATTER_FORWARD));
1600:     PetscUseTypeMethod(b, diagonalscale, NULL, baij->lvec);
1601:   }
1602:   PetscFunctionReturn(PETSC_SUCCESS);
1603: }

1605: static PetscErrorCode MatZeroRows_MPIBAIJ(Mat A, PetscInt N, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
1606: {
1607:   Mat_MPIBAIJ *l = (Mat_MPIBAIJ *)A->data;
1608:   PetscInt    *lrows;
1609:   PetscInt     r, len;
1610:   PetscBool    cong;

1612:   PetscFunctionBegin;
1613:   /* get locally owned rows */
1614:   PetscCall(MatZeroRowsMapLocal_Private(A, N, rows, &len, &lrows));
1615:   /* fix right-hand side if needed */
1616:   if (x && b) {
1617:     const PetscScalar *xx;
1618:     PetscScalar       *bb;

1620:     PetscCall(VecGetArrayRead(x, &xx));
1621:     PetscCall(VecGetArray(b, &bb));
1622:     for (r = 0; r < len; ++r) bb[lrows[r]] = diag * xx[lrows[r]];
1623:     PetscCall(VecRestoreArrayRead(x, &xx));
1624:     PetscCall(VecRestoreArray(b, &bb));
1625:   }

1627:   /* actually zap the local rows */
1628:   /*
1629:         Zero the required rows. If the "diagonal block" of the matrix
1630:      is square and the user wishes to set the diagonal we use separate
1631:      code so that MatSetValues() is not called for each diagonal allocating
1632:      new memory, thus calling lots of mallocs and slowing things down.

1634:   */
1635:   /* must zero l->B before l->A because the (diag) case below may put values into l->B*/
1636:   PetscCall(MatZeroRows_SeqBAIJ(l->B, len, lrows, 0.0, NULL, NULL));
1637:   PetscCall(MatHasCongruentLayouts(A, &cong));
1638:   if ((diag != 0.0) && cong) {
1639:     PetscCall(MatZeroRows_SeqBAIJ(l->A, len, lrows, diag, NULL, NULL));
1640:   } else if (diag != 0.0) {
1641:     PetscCall(MatZeroRows_SeqBAIJ(l->A, len, lrows, 0.0, NULL, NULL));
1642:     PetscCheck(!((Mat_SeqBAIJ *)l->A->data)->nonew, PETSC_COMM_SELF, PETSC_ERR_SUP, "MatZeroRows() on rectangular matrices cannot be used with the Mat options MAT_NEW_NONZERO_LOCATIONS, MAT_NEW_NONZERO_LOCATION_ERR, and MAT_NEW_NONZERO_ALLOCATION_ERR");
1643:     for (r = 0; r < len; ++r) {
1644:       const PetscInt row = lrows[r] + A->rmap->rstart;
1645:       PetscCall(MatSetValues(A, 1, &row, 1, &row, &diag, INSERT_VALUES));
1646:     }
1647:     PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
1648:     PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
1649:   } else {
1650:     PetscCall(MatZeroRows_SeqBAIJ(l->A, len, lrows, 0.0, NULL, NULL));
1651:   }
1652:   PetscCall(PetscFree(lrows));

1654:   /* only change matrix nonzero state if pattern was allowed to be changed */
1655:   if (!((Mat_SeqBAIJ *)l->A->data)->keepnonzeropattern || !((Mat_SeqBAIJ *)l->A->data)->nonew) {
1656:     PetscObjectState state = l->A->nonzerostate + l->B->nonzerostate;
1657:     PetscCallMPI(MPIU_Allreduce(&state, &A->nonzerostate, 1, MPIU_INT64, MPI_SUM, PetscObjectComm((PetscObject)A)));
1658:   }
1659:   PetscFunctionReturn(PETSC_SUCCESS);
1660: }

1662: static PetscErrorCode MatZeroRowsColumns_MPIBAIJ(Mat A, PetscInt N, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
1663: {
1664:   Mat_MPIBAIJ       *l = (Mat_MPIBAIJ *)A->data;
1665:   PetscMPIInt        n, p = 0;
1666:   PetscInt           i, j, k, r, len = 0, row, col, count;
1667:   PetscInt          *lrows, *owners = A->rmap->range;
1668:   PetscSFNode       *rrows;
1669:   PetscSF            sf;
1670:   const PetscScalar *xx;
1671:   PetscScalar       *bb, *mask;
1672:   Vec                xmask, lmask;
1673:   Mat_SeqBAIJ       *baij = (Mat_SeqBAIJ *)l->B->data;
1674:   PetscInt           bs = A->rmap->bs, bs2 = baij->bs2;
1675:   PetscScalar       *aa;

1677:   PetscFunctionBegin;
1678:   PetscCall(PetscMPIIntCast(A->rmap->n, &n));
1679:   /* Create SF where leaves are input rows and roots are owned rows */
1680:   PetscCall(PetscMalloc1(n, &lrows));
1681:   for (r = 0; r < n; ++r) lrows[r] = -1;
1682:   PetscCall(PetscMalloc1(N, &rrows));
1683:   for (r = 0; r < N; ++r) {
1684:     const PetscInt idx = rows[r];
1685:     PetscCheck(idx >= 0 && A->rmap->N > idx, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row %" PetscInt_FMT " out of range [0,%" PetscInt_FMT ")", idx, A->rmap->N);
1686:     if (idx < owners[p] || owners[p + 1] <= idx) { /* short-circuit the search if the last p owns this row too */
1687:       PetscCall(PetscLayoutFindOwner(A->rmap, idx, &p));
1688:     }
1689:     rrows[r].rank  = p;
1690:     rrows[r].index = rows[r] - owners[p];
1691:   }
1692:   PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)A), &sf));
1693:   PetscCall(PetscSFSetGraph(sf, n, N, NULL, PETSC_OWN_POINTER, rrows, PETSC_OWN_POINTER));
1694:   /* Collect flags for rows to be zeroed */
1695:   PetscCall(PetscSFReduceBegin(sf, MPIU_INT, (PetscInt *)rows, lrows, MPI_LOR));
1696:   PetscCall(PetscSFReduceEnd(sf, MPIU_INT, (PetscInt *)rows, lrows, MPI_LOR));
1697:   PetscCall(PetscSFDestroy(&sf));
1698:   /* Compress and put in row numbers */
1699:   for (r = 0; r < n; ++r)
1700:     if (lrows[r] >= 0) lrows[len++] = r;
1701:   /* zero diagonal part of matrix */
1702:   PetscCall(MatZeroRowsColumns(l->A, len, lrows, diag, x, b));
1703:   /* handle off-diagonal part of matrix */
1704:   PetscCall(MatCreateVecs(A, &xmask, NULL));
1705:   PetscCall(VecDuplicate(l->lvec, &lmask));
1706:   PetscCall(VecGetArray(xmask, &bb));
1707:   for (i = 0; i < len; i++) bb[lrows[i]] = 1;
1708:   PetscCall(VecRestoreArray(xmask, &bb));
1709:   PetscCall(VecScatterBegin(l->Mvctx, xmask, lmask, ADD_VALUES, SCATTER_FORWARD));
1710:   PetscCall(VecScatterEnd(l->Mvctx, xmask, lmask, ADD_VALUES, SCATTER_FORWARD));
1711:   PetscCall(VecDestroy(&xmask));
1712:   if (x) {
1713:     PetscCall(VecScatterBegin(l->Mvctx, x, l->lvec, INSERT_VALUES, SCATTER_FORWARD));
1714:     PetscCall(VecScatterEnd(l->Mvctx, x, l->lvec, INSERT_VALUES, SCATTER_FORWARD));
1715:     PetscCall(VecGetArrayRead(l->lvec, &xx));
1716:     PetscCall(VecGetArray(b, &bb));
1717:   }
1718:   PetscCall(VecGetArray(lmask, &mask));
1719:   /* remove zeroed rows of off-diagonal matrix */
1720:   for (i = 0; i < len; ++i) {
1721:     row   = lrows[i];
1722:     count = (baij->i[row / bs + 1] - baij->i[row / bs]) * bs;
1723:     aa    = baij->a + baij->i[row / bs] * bs2 + (row % bs);
1724:     for (k = 0; k < count; ++k) {
1725:       aa[0] = 0.0;
1726:       aa += bs;
1727:     }
1728:   }
1729:   /* loop over all elements of off process part of matrix zeroing removed columns*/
1730:   for (i = 0; i < l->B->rmap->N; ++i) {
1731:     row = i / bs;
1732:     for (j = baij->i[row]; j < baij->i[row + 1]; ++j) {
1733:       for (k = 0; k < bs; ++k) {
1734:         col = bs * baij->j[j] + k;
1735:         if (PetscAbsScalar(mask[col])) {
1736:           aa = baij->a + j * bs2 + (i % bs) + bs * k;
1737:           if (x) bb[i] -= aa[0] * xx[col];
1738:           aa[0] = 0.0;
1739:         }
1740:       }
1741:     }
1742:   }
1743:   if (x) {
1744:     PetscCall(VecRestoreArray(b, &bb));
1745:     PetscCall(VecRestoreArrayRead(l->lvec, &xx));
1746:   }
1747:   PetscCall(VecRestoreArray(lmask, &mask));
1748:   PetscCall(VecDestroy(&lmask));
1749:   PetscCall(PetscFree(lrows));

1751:   /* only change matrix nonzero state if pattern was allowed to be changed */
1752:   if (!((Mat_SeqBAIJ *)l->A->data)->nonew) {
1753:     PetscObjectState state = l->A->nonzerostate + l->B->nonzerostate;
1754:     PetscCallMPI(MPIU_Allreduce(&state, &A->nonzerostate, 1, MPIU_INT64, MPI_SUM, PetscObjectComm((PetscObject)A)));
1755:   }
1756:   PetscFunctionReturn(PETSC_SUCCESS);
1757: }

1759: static PetscErrorCode MatSetUnfactored_MPIBAIJ(Mat A)
1760: {
1761:   Mat_MPIBAIJ *a = (Mat_MPIBAIJ *)A->data;

1763:   PetscFunctionBegin;
1764:   PetscCall(MatSetUnfactored(a->A));
1765:   PetscFunctionReturn(PETSC_SUCCESS);
1766: }

1768: static PetscErrorCode MatDuplicate_MPIBAIJ(Mat, MatDuplicateOption, Mat *);

1770: static PetscErrorCode MatEqual_MPIBAIJ(Mat A, Mat B, PetscBool *flag)
1771: {
1772:   Mat_MPIBAIJ *matB = (Mat_MPIBAIJ *)B->data, *matA = (Mat_MPIBAIJ *)A->data;
1773:   Mat          a, b, c, d;
1774:   PetscBool    flg;

1776:   PetscFunctionBegin;
1777:   a = matA->A;
1778:   b = matA->B;
1779:   c = matB->A;
1780:   d = matB->B;

1782:   PetscCall(MatEqual(a, c, &flg));
1783:   if (flg) PetscCall(MatEqual(b, d, &flg));
1784:   PetscCallMPI(MPIU_Allreduce(&flg, flag, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A)));
1785:   PetscFunctionReturn(PETSC_SUCCESS);
1786: }

1788: static PetscErrorCode MatCopy_MPIBAIJ(Mat A, Mat B, MatStructure str)
1789: {
1790:   Mat_MPIBAIJ *a = (Mat_MPIBAIJ *)A->data;
1791:   Mat_MPIBAIJ *b = (Mat_MPIBAIJ *)B->data;

1793:   PetscFunctionBegin;
1794:   /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */
1795:   if ((str != SAME_NONZERO_PATTERN) || (A->ops->copy != B->ops->copy)) {
1796:     PetscCall(MatCopy_Basic(A, B, str));
1797:   } else {
1798:     PetscCall(MatCopy(a->A, b->A, str));
1799:     PetscCall(MatCopy(a->B, b->B, str));
1800:   }
1801:   PetscCall(PetscObjectStateIncrease((PetscObject)B));
1802:   PetscFunctionReturn(PETSC_SUCCESS);
1803: }

1805: PetscErrorCode MatAXPYGetPreallocation_MPIBAIJ(Mat Y, const PetscInt *yltog, Mat X, const PetscInt *xltog, PetscInt *nnz)
1806: {
1807:   PetscInt     bs = Y->rmap->bs, m = Y->rmap->N / bs;
1808:   Mat_SeqBAIJ *x = (Mat_SeqBAIJ *)X->data;
1809:   Mat_SeqBAIJ *y = (Mat_SeqBAIJ *)Y->data;

1811:   PetscFunctionBegin;
1812:   PetscCall(MatAXPYGetPreallocation_MPIX_private(m, x->i, x->j, xltog, y->i, y->j, yltog, nnz));
1813:   PetscFunctionReturn(PETSC_SUCCESS);
1814: }

1816: static PetscErrorCode MatAXPY_MPIBAIJ(Mat Y, PetscScalar a, Mat X, MatStructure str)
1817: {
1818:   Mat_MPIBAIJ *xx = (Mat_MPIBAIJ *)X->data, *yy = (Mat_MPIBAIJ *)Y->data;
1819:   PetscBLASInt bnz, one                         = 1;
1820:   Mat_SeqBAIJ *x, *y;
1821:   PetscInt     bs2 = Y->rmap->bs * Y->rmap->bs;

1823:   PetscFunctionBegin;
1824:   if (str == SAME_NONZERO_PATTERN) {
1825:     PetscScalar alpha = a;
1826:     x                 = (Mat_SeqBAIJ *)xx->A->data;
1827:     y                 = (Mat_SeqBAIJ *)yy->A->data;
1828:     PetscCall(PetscBLASIntCast(x->nz * bs2, &bnz));
1829:     PetscCallBLAS("BLASaxpy", BLASaxpy_(&bnz, &alpha, x->a, &one, y->a, &one));
1830:     x = (Mat_SeqBAIJ *)xx->B->data;
1831:     y = (Mat_SeqBAIJ *)yy->B->data;
1832:     PetscCall(PetscBLASIntCast(x->nz * bs2, &bnz));
1833:     PetscCallBLAS("BLASaxpy", BLASaxpy_(&bnz, &alpha, x->a, &one, y->a, &one));
1834:     PetscCall(PetscObjectStateIncrease((PetscObject)Y));
1835:   } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */
1836:     PetscCall(MatAXPY_Basic(Y, a, X, str));
1837:   } else {
1838:     Mat       B;
1839:     PetscInt *nnz_d, *nnz_o, bs = Y->rmap->bs;
1840:     PetscCall(PetscMalloc1(yy->A->rmap->N, &nnz_d));
1841:     PetscCall(PetscMalloc1(yy->B->rmap->N, &nnz_o));
1842:     PetscCall(MatCreate(PetscObjectComm((PetscObject)Y), &B));
1843:     PetscCall(PetscObjectSetName((PetscObject)B, ((PetscObject)Y)->name));
1844:     PetscCall(MatSetSizes(B, Y->rmap->n, Y->cmap->n, Y->rmap->N, Y->cmap->N));
1845:     PetscCall(MatSetBlockSizesFromMats(B, Y, Y));
1846:     PetscCall(MatSetType(B, MATMPIBAIJ));
1847:     PetscCall(MatAXPYGetPreallocation_SeqBAIJ(yy->A, xx->A, nnz_d));
1848:     PetscCall(MatAXPYGetPreallocation_MPIBAIJ(yy->B, yy->garray, xx->B, xx->garray, nnz_o));
1849:     PetscCall(MatMPIBAIJSetPreallocation(B, bs, 0, nnz_d, 0, nnz_o));
1850:     /* MatAXPY_BasicWithPreallocation() for BAIJ matrix is much slower than AIJ, even for bs=1 ! */
1851:     PetscCall(MatAXPY_BasicWithPreallocation(B, Y, a, X, str));
1852:     PetscCall(MatHeaderMerge(Y, &B));
1853:     PetscCall(PetscFree(nnz_d));
1854:     PetscCall(PetscFree(nnz_o));
1855:   }
1856:   PetscFunctionReturn(PETSC_SUCCESS);
1857: }

1859: static PetscErrorCode MatConjugate_MPIBAIJ(Mat mat)
1860: {
1861:   PetscFunctionBegin;
1862:   if (PetscDefined(USE_COMPLEX)) {
1863:     Mat_MPIBAIJ *a = (Mat_MPIBAIJ *)mat->data;

1865:     PetscCall(MatConjugate_SeqBAIJ(a->A));
1866:     PetscCall(MatConjugate_SeqBAIJ(a->B));
1867:   }
1868:   PetscFunctionReturn(PETSC_SUCCESS);
1869: }

1871: static PetscErrorCode MatRealPart_MPIBAIJ(Mat A)
1872: {
1873:   Mat_MPIBAIJ *a = (Mat_MPIBAIJ *)A->data;

1875:   PetscFunctionBegin;
1876:   PetscCall(MatRealPart(a->A));
1877:   PetscCall(MatRealPart(a->B));
1878:   PetscFunctionReturn(PETSC_SUCCESS);
1879: }

1881: static PetscErrorCode MatImaginaryPart_MPIBAIJ(Mat A)
1882: {
1883:   Mat_MPIBAIJ *a = (Mat_MPIBAIJ *)A->data;

1885:   PetscFunctionBegin;
1886:   PetscCall(MatImaginaryPart(a->A));
1887:   PetscCall(MatImaginaryPart(a->B));
1888:   PetscFunctionReturn(PETSC_SUCCESS);
1889: }

1891: static PetscErrorCode MatCreateSubMatrix_MPIBAIJ(Mat mat, IS isrow, IS iscol, MatReuse call, Mat *newmat)
1892: {
1893:   IS       iscol_local;
1894:   PetscInt csize;

1896:   PetscFunctionBegin;
1897:   PetscCall(ISGetLocalSize(iscol, &csize));
1898:   if (call == MAT_REUSE_MATRIX) {
1899:     PetscCall(PetscObjectQuery((PetscObject)*newmat, "ISAllGather", (PetscObject *)&iscol_local));
1900:     PetscCheck(iscol_local, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Submatrix passed in was not used before, cannot reuse");
1901:   } else {
1902:     PetscCall(ISAllGather(iscol, &iscol_local));
1903:   }
1904:   PetscCall(MatCreateSubMatrix_MPIBAIJ_Private(mat, isrow, iscol_local, csize, call, newmat, PETSC_FALSE));
1905:   if (call == MAT_INITIAL_MATRIX) {
1906:     PetscCall(PetscObjectCompose((PetscObject)*newmat, "ISAllGather", (PetscObject)iscol_local));
1907:     PetscCall(ISDestroy(&iscol_local));
1908:   }
1909:   PetscFunctionReturn(PETSC_SUCCESS);
1910: }

1912: /*
1913:   Not great since it makes two copies of the submatrix, first an SeqBAIJ
1914:   in local and then by concatenating the local matrices the end result.
1915:   Writing it directly would be much like MatCreateSubMatrices_MPIBAIJ().
1916:   This routine is used for BAIJ and SBAIJ matrices (unfortunate dependency).
1917: */
1918: PetscErrorCode MatCreateSubMatrix_MPIBAIJ_Private(Mat mat, IS isrow, IS iscol, PetscInt csize, MatReuse call, Mat *newmat, PetscBool sym)
1919: {
1920:   PetscMPIInt  rank, size;
1921:   PetscInt     i, m, n, rstart, row, rend, nz, *cwork, j, bs;
1922:   PetscInt    *ii, *jj, nlocal, *dlens, *olens, dlen, olen, jend, mglobal;
1923:   Mat          M, Mreuse;
1924:   MatScalar   *vwork, *aa;
1925:   MPI_Comm     comm;
1926:   IS           isrow_new, iscol_new;
1927:   Mat_SeqBAIJ *aij;

1929:   PetscFunctionBegin;
1930:   PetscCall(PetscObjectGetComm((PetscObject)mat, &comm));
1931:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
1932:   PetscCallMPI(MPI_Comm_size(comm, &size));
1933:   /* The compression and expansion should be avoided. Doesn't point
1934:      out errors, might change the indices, hence buggey */
1935:   PetscCall(ISCompressIndicesGeneral(mat->rmap->N, mat->rmap->n, mat->rmap->bs, 1, &isrow, &isrow_new));
1936:   if (isrow == iscol) {
1937:     iscol_new = isrow_new;
1938:     PetscCall(PetscObjectReference((PetscObject)iscol_new));
1939:   } else PetscCall(ISCompressIndicesGeneral(mat->cmap->N, mat->cmap->n, mat->cmap->bs, 1, &iscol, &iscol_new));

1941:   if (call == MAT_REUSE_MATRIX) {
1942:     PetscCall(PetscObjectQuery((PetscObject)*newmat, "SubMatrix", (PetscObject *)&Mreuse));
1943:     PetscCheck(Mreuse, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Submatrix passed in was not used before, cannot reuse");
1944:     PetscCall(MatCreateSubMatrices_MPIBAIJ_local(mat, 1, &isrow_new, &iscol_new, MAT_REUSE_MATRIX, &Mreuse, sym));
1945:   } else {
1946:     PetscCall(MatCreateSubMatrices_MPIBAIJ_local(mat, 1, &isrow_new, &iscol_new, MAT_INITIAL_MATRIX, &Mreuse, sym));
1947:   }
1948:   PetscCall(ISDestroy(&isrow_new));
1949:   PetscCall(ISDestroy(&iscol_new));
1950:   /*
1951:       m - number of local rows
1952:       n - number of columns (same on all processors)
1953:       rstart - first row in new global matrix generated
1954:   */
1955:   PetscCall(MatGetBlockSize(mat, &bs));
1956:   PetscCall(MatGetSize(Mreuse, &m, &n));
1957:   m = m / bs;
1958:   n = n / bs;

1960:   if (call == MAT_INITIAL_MATRIX) {
1961:     aij = (Mat_SeqBAIJ *)Mreuse->data;
1962:     ii  = aij->i;
1963:     jj  = aij->j;

1965:     /*
1966:         Determine the number of non-zeros in the diagonal and off-diagonal
1967:         portions of the matrix in order to do correct preallocation
1968:     */

1970:     /* first get start and end of "diagonal" columns */
1971:     if (csize == PETSC_DECIDE) {
1972:       PetscCall(ISGetSize(isrow, &mglobal));
1973:       if (mglobal == n * bs) { /* square matrix */
1974:         nlocal = m;
1975:       } else {
1976:         nlocal = n / size + ((n % size) > rank);
1977:       }
1978:     } else {
1979:       nlocal = csize / bs;
1980:     }
1981:     PetscCallMPI(MPI_Scan(&nlocal, &rend, 1, MPIU_INT, MPI_SUM, comm));
1982:     rstart = rend - nlocal;
1983:     PetscCheck(rank != size - 1 || rend == n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local column sizes %" PetscInt_FMT " do not add up to total number of columns %" PetscInt_FMT, rend, n);

1985:     /* next, compute all the lengths */
1986:     PetscCall(PetscMalloc2(m + 1, &dlens, m + 1, &olens));
1987:     for (i = 0; i < m; i++) {
1988:       jend = ii[i + 1] - ii[i];
1989:       olen = 0;
1990:       dlen = 0;
1991:       for (j = 0; j < jend; j++) {
1992:         if (*jj < rstart || *jj >= rend) olen++;
1993:         else dlen++;
1994:         jj++;
1995:       }
1996:       olens[i] = olen;
1997:       dlens[i] = dlen;
1998:     }
1999:     PetscCall(MatCreate(comm, &M));
2000:     PetscCall(MatSetSizes(M, bs * m, bs * nlocal, PETSC_DECIDE, bs * n));
2001:     PetscCall(MatSetType(M, sym ? ((PetscObject)mat)->type_name : MATMPIBAIJ));
2002:     PetscCall(MatMPIBAIJSetPreallocation(M, bs, 0, dlens, 0, olens));
2003:     PetscCall(MatMPISBAIJSetPreallocation(M, bs, 0, dlens, 0, olens));
2004:     PetscCall(PetscFree2(dlens, olens));
2005:   } else {
2006:     PetscInt ml, nl;

2008:     M = *newmat;
2009:     PetscCall(MatGetLocalSize(M, &ml, &nl));
2010:     PetscCheck(ml == m, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Previous matrix must be same size/layout as request");
2011:     PetscCall(MatZeroEntries(M));
2012:     /*
2013:          The next two lines are needed so we may call MatSetValues_MPIAIJ() below directly,
2014:        rather than the slower MatSetValues().
2015:     */
2016:     M->was_assembled = PETSC_TRUE;
2017:     M->assembled     = PETSC_FALSE;
2018:   }
2019:   PetscCall(MatSetOption(M, MAT_ROW_ORIENTED, PETSC_FALSE));
2020:   PetscCall(MatGetOwnershipRange(M, &rstart, &rend));
2021:   aij = (Mat_SeqBAIJ *)Mreuse->data;
2022:   ii  = aij->i;
2023:   jj  = aij->j;
2024:   aa  = aij->a;
2025:   for (i = 0; i < m; i++) {
2026:     row   = rstart / bs + i;
2027:     nz    = ii[i + 1] - ii[i];
2028:     cwork = jj;
2029:     jj    = PetscSafePointerPlusOffset(jj, nz);
2030:     vwork = aa;
2031:     aa    = PetscSafePointerPlusOffset(aa, nz * bs * bs);
2032:     PetscCall(MatSetValuesBlocked_MPIBAIJ(M, 1, &row, nz, cwork, vwork, INSERT_VALUES));
2033:   }

2035:   PetscCall(MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY));
2036:   PetscCall(MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY));
2037:   *newmat = M;

2039:   /* save submatrix used in processor for next request */
2040:   if (call == MAT_INITIAL_MATRIX) {
2041:     PetscCall(PetscObjectCompose((PetscObject)M, "SubMatrix", (PetscObject)Mreuse));
2042:     PetscCall(PetscObjectDereference((PetscObject)Mreuse));
2043:   }
2044:   PetscFunctionReturn(PETSC_SUCCESS);
2045: }

2047: static PetscErrorCode MatPermute_MPIBAIJ(Mat A, IS rowp, IS colp, Mat *B)
2048: {
2049:   MPI_Comm        comm, pcomm;
2050:   PetscInt        clocal_size, nrows;
2051:   const PetscInt *rows;
2052:   PetscMPIInt     size;
2053:   IS              crowp, lcolp;

2055:   PetscFunctionBegin;
2056:   PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
2057:   /* make a collective version of 'rowp' */
2058:   PetscCall(PetscObjectGetComm((PetscObject)rowp, &pcomm));
2059:   if (pcomm == comm) {
2060:     crowp = rowp;
2061:   } else {
2062:     PetscCall(ISGetSize(rowp, &nrows));
2063:     PetscCall(ISGetIndices(rowp, &rows));
2064:     PetscCall(ISCreateGeneral(comm, nrows, rows, PETSC_COPY_VALUES, &crowp));
2065:     PetscCall(ISRestoreIndices(rowp, &rows));
2066:   }
2067:   PetscCall(ISSetPermutation(crowp));
2068:   /* make a local version of 'colp' */
2069:   PetscCall(PetscObjectGetComm((PetscObject)colp, &pcomm));
2070:   PetscCallMPI(MPI_Comm_size(pcomm, &size));
2071:   if (size == 1) {
2072:     lcolp = colp;
2073:   } else {
2074:     PetscCall(ISAllGather(colp, &lcolp));
2075:   }
2076:   PetscCall(ISSetPermutation(lcolp));
2077:   /* now we just get the submatrix */
2078:   PetscCall(MatGetLocalSize(A, NULL, &clocal_size));
2079:   PetscCall(MatCreateSubMatrix_MPIBAIJ_Private(A, crowp, lcolp, clocal_size, MAT_INITIAL_MATRIX, B, PETSC_FALSE));
2080:   /* clean up */
2081:   if (pcomm != comm) PetscCall(ISDestroy(&crowp));
2082:   if (size > 1) PetscCall(ISDestroy(&lcolp));
2083:   PetscFunctionReturn(PETSC_SUCCESS);
2084: }

2086: static PetscErrorCode MatGetGhosts_MPIBAIJ(Mat mat, PetscInt *nghosts, const PetscInt *ghosts[])
2087: {
2088:   Mat_MPIBAIJ *baij = (Mat_MPIBAIJ *)mat->data;
2089:   Mat_SeqBAIJ *B    = (Mat_SeqBAIJ *)baij->B->data;

2091:   PetscFunctionBegin;
2092:   if (nghosts) *nghosts = B->nbs;
2093:   if (ghosts) *ghosts = baij->garray;
2094:   PetscFunctionReturn(PETSC_SUCCESS);
2095: }

2097: static PetscErrorCode MatGetSeqNonzeroStructure_MPIBAIJ(Mat A, Mat *newmat)
2098: {
2099:   Mat          B;
2100:   Mat_MPIBAIJ *a  = (Mat_MPIBAIJ *)A->data;
2101:   Mat_SeqBAIJ *ad = (Mat_SeqBAIJ *)a->A->data, *bd = (Mat_SeqBAIJ *)a->B->data;
2102:   Mat_SeqAIJ  *b;
2103:   PetscMPIInt  size, rank, *recvcounts = NULL, *displs = NULL;
2104:   PetscInt     sendcount, i, *rstarts = A->rmap->range, n, cnt, j, bs = A->rmap->bs;
2105:   PetscInt     m, *garray = a->garray, *lens, *jsendbuf, *a_jsendbuf, *b_jsendbuf;

2107:   PetscFunctionBegin;
2108:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
2109:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank));

2111:   /*   Tell every processor the number of nonzeros per row  */
2112:   PetscCall(PetscMalloc1(A->rmap->N / bs, &lens));
2113:   for (i = A->rmap->rstart / bs; i < A->rmap->rend / bs; i++) lens[i] = ad->i[i - A->rmap->rstart / bs + 1] - ad->i[i - A->rmap->rstart / bs] + bd->i[i - A->rmap->rstart / bs + 1] - bd->i[i - A->rmap->rstart / bs];
2114:   PetscCall(PetscMalloc1(2 * size, &recvcounts));
2115:   displs = recvcounts + size;
2116:   for (i = 0; i < size; i++) {
2117:     PetscCall(PetscMPIIntCast(A->rmap->range[i + 1] / bs - A->rmap->range[i] / bs, &recvcounts[i]));
2118:     PetscCall(PetscMPIIntCast(A->rmap->range[i] / bs, &displs[i]));
2119:   }
2120:   PetscCallMPI(MPI_Allgatherv(MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, lens, recvcounts, displs, MPIU_INT, PetscObjectComm((PetscObject)A)));
2121:   /* Create the sequential matrix of the same type as the local block diagonal  */
2122:   PetscCall(MatCreate(PETSC_COMM_SELF, &B));
2123:   PetscCall(MatSetSizes(B, A->rmap->N / bs, A->cmap->N / bs, PETSC_DETERMINE, PETSC_DETERMINE));
2124:   PetscCall(MatSetType(B, MATSEQAIJ));
2125:   PetscCall(MatSeqAIJSetPreallocation(B, 0, lens));
2126:   b = (Mat_SeqAIJ *)B->data;

2128:   /*     Copy my part of matrix column indices over  */
2129:   sendcount  = ad->nz + bd->nz;
2130:   jsendbuf   = b->j + b->i[rstarts[rank] / bs];
2131:   a_jsendbuf = ad->j;
2132:   b_jsendbuf = bd->j;
2133:   n          = A->rmap->rend / bs - A->rmap->rstart / bs;
2134:   cnt        = 0;
2135:   for (i = 0; i < n; i++) {
2136:     /* put in lower diagonal portion */
2137:     m = bd->i[i + 1] - bd->i[i];
2138:     while (m > 0) {
2139:       /* is it above diagonal (in bd (compressed) numbering) */
2140:       if (garray[*b_jsendbuf] > A->rmap->rstart / bs + i) break;
2141:       jsendbuf[cnt++] = garray[*b_jsendbuf++];
2142:       m--;
2143:     }

2145:     /* put in diagonal portion */
2146:     for (j = ad->i[i]; j < ad->i[i + 1]; j++) jsendbuf[cnt++] = A->rmap->rstart / bs + *a_jsendbuf++;

2148:     /* put in upper diagonal portion */
2149:     while (m-- > 0) jsendbuf[cnt++] = garray[*b_jsendbuf++];
2150:   }
2151:   PetscCheck(cnt == sendcount, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Corrupted PETSc matrix: nz given %" PetscInt_FMT " actual nz %" PetscInt_FMT, sendcount, cnt);

2153:   /*  Gather all column indices to all processors  */
2154:   for (i = 0; i < size; i++) {
2155:     recvcounts[i] = 0;
2156:     for (j = A->rmap->range[i] / bs; j < A->rmap->range[i + 1] / bs; j++) recvcounts[i] += lens[j];
2157:   }
2158:   displs[0] = 0;
2159:   for (i = 1; i < size; i++) displs[i] = displs[i - 1] + recvcounts[i - 1];
2160:   PetscCallMPI(MPI_Allgatherv(MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, b->j, recvcounts, displs, MPIU_INT, PetscObjectComm((PetscObject)A)));
2161:   /*  Assemble the matrix into usable form (note numerical values not yet set)  */
2162:   /* set the b->ilen (length of each row) values */
2163:   PetscCall(PetscArraycpy(b->ilen, lens, A->rmap->N / bs));
2164:   /* set the b->i indices */
2165:   b->i[0] = 0;
2166:   for (i = 1; i <= A->rmap->N / bs; i++) b->i[i] = b->i[i - 1] + lens[i - 1];
2167:   PetscCall(PetscFree(lens));
2168:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
2169:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
2170:   PetscCall(PetscFree(recvcounts));

2172:   PetscCall(MatPropagateSymmetryOptions(A, B));
2173:   *newmat = B;
2174:   PetscFunctionReturn(PETSC_SUCCESS);
2175: }

2177: static PetscErrorCode MatSOR_MPIBAIJ(Mat matin, Vec bb, PetscReal omega, MatSORType flag, PetscReal fshift, PetscInt its, PetscInt lits, Vec xx)
2178: {
2179:   Mat_MPIBAIJ *mat = (Mat_MPIBAIJ *)matin->data;
2180:   Vec          bb1 = NULL;

2182:   PetscFunctionBegin;
2183:   if (flag == SOR_APPLY_UPPER) {
2184:     PetscCall((*mat->A->ops->sor)(mat->A, bb, omega, flag, fshift, lits, 1, xx));
2185:     PetscFunctionReturn(PETSC_SUCCESS);
2186:   }

2188:   if (its > 1 || ~flag & SOR_ZERO_INITIAL_GUESS) PetscCall(VecDuplicate(bb, &bb1));

2190:   if ((flag & SOR_LOCAL_SYMMETRIC_SWEEP) == SOR_LOCAL_SYMMETRIC_SWEEP) {
2191:     if (flag & SOR_ZERO_INITIAL_GUESS) {
2192:       PetscCall((*mat->A->ops->sor)(mat->A, bb, omega, flag, fshift, lits, 1, xx));
2193:       its--;
2194:     }

2196:     while (its--) {
2197:       PetscCall(VecScatterBegin(mat->Mvctx, xx, mat->lvec, INSERT_VALUES, SCATTER_FORWARD));
2198:       PetscCall(VecScatterEnd(mat->Mvctx, xx, mat->lvec, INSERT_VALUES, SCATTER_FORWARD));

2200:       /* update rhs: bb1 = bb - B*x */
2201:       PetscCall(VecScale(mat->lvec, -1.0));
2202:       PetscCall((*mat->B->ops->multadd)(mat->B, mat->lvec, bb, bb1));

2204:       /* local sweep */
2205:       PetscCall((*mat->A->ops->sor)(mat->A, bb1, omega, SOR_SYMMETRIC_SWEEP, fshift, lits, 1, xx));
2206:     }
2207:   } else if (flag & SOR_LOCAL_FORWARD_SWEEP) {
2208:     if (flag & SOR_ZERO_INITIAL_GUESS) {
2209:       PetscCall((*mat->A->ops->sor)(mat->A, bb, omega, flag, fshift, lits, 1, xx));
2210:       its--;
2211:     }
2212:     while (its--) {
2213:       PetscCall(VecScatterBegin(mat->Mvctx, xx, mat->lvec, INSERT_VALUES, SCATTER_FORWARD));
2214:       PetscCall(VecScatterEnd(mat->Mvctx, xx, mat->lvec, INSERT_VALUES, SCATTER_FORWARD));

2216:       /* update rhs: bb1 = bb - B*x */
2217:       PetscCall(VecScale(mat->lvec, -1.0));
2218:       PetscCall((*mat->B->ops->multadd)(mat->B, mat->lvec, bb, bb1));

2220:       /* local sweep */
2221:       PetscCall((*mat->A->ops->sor)(mat->A, bb1, omega, SOR_FORWARD_SWEEP, fshift, lits, 1, xx));
2222:     }
2223:   } else if (flag & SOR_LOCAL_BACKWARD_SWEEP) {
2224:     if (flag & SOR_ZERO_INITIAL_GUESS) {
2225:       PetscCall((*mat->A->ops->sor)(mat->A, bb, omega, flag, fshift, lits, 1, xx));
2226:       its--;
2227:     }
2228:     while (its--) {
2229:       PetscCall(VecScatterBegin(mat->Mvctx, xx, mat->lvec, INSERT_VALUES, SCATTER_FORWARD));
2230:       PetscCall(VecScatterEnd(mat->Mvctx, xx, mat->lvec, INSERT_VALUES, SCATTER_FORWARD));

2232:       /* update rhs: bb1 = bb - B*x */
2233:       PetscCall(VecScale(mat->lvec, -1.0));
2234:       PetscCall((*mat->B->ops->multadd)(mat->B, mat->lvec, bb, bb1));

2236:       /* local sweep */
2237:       PetscCall((*mat->A->ops->sor)(mat->A, bb1, omega, SOR_BACKWARD_SWEEP, fshift, lits, 1, xx));
2238:     }
2239:   } else SETERRQ(PetscObjectComm((PetscObject)matin), PETSC_ERR_SUP, "Parallel version of SOR requested not supported");

2241:   PetscCall(VecDestroy(&bb1));
2242:   PetscFunctionReturn(PETSC_SUCCESS);
2243: }

2245: static PetscErrorCode MatGetColumnReductions_MPIBAIJ(Mat A, PetscInt type, PetscReal *reductions)
2246: {
2247:   Mat_MPIBAIJ *aij = (Mat_MPIBAIJ *)A->data;
2248:   PetscInt     m, N, i, *garray = aij->garray;
2249:   PetscInt     ib, jb, bs = A->rmap->bs;
2250:   Mat_SeqBAIJ *a_aij = (Mat_SeqBAIJ *)aij->A->data;
2251:   MatScalar   *a_val = a_aij->a;
2252:   Mat_SeqBAIJ *b_aij = (Mat_SeqBAIJ *)aij->B->data;
2253:   MatScalar   *b_val = b_aij->a;
2254:   PetscReal   *work;

2256:   PetscFunctionBegin;
2257:   PetscCall(MatGetSize(A, &m, &N));
2258:   PetscCall(PetscCalloc1(N, &work));
2259:   if (type == NORM_2) {
2260:     for (i = a_aij->i[0]; i < a_aij->i[aij->A->rmap->n / bs]; i++) {
2261:       for (jb = 0; jb < bs; jb++) {
2262:         for (ib = 0; ib < bs; ib++) {
2263:           work[A->cmap->rstart + a_aij->j[i] * bs + jb] += PetscAbsScalar(*a_val * *a_val);
2264:           a_val++;
2265:         }
2266:       }
2267:     }
2268:     for (i = b_aij->i[0]; i < b_aij->i[aij->B->rmap->n / bs]; i++) {
2269:       for (jb = 0; jb < bs; jb++) {
2270:         for (ib = 0; ib < bs; ib++) {
2271:           work[garray[b_aij->j[i]] * bs + jb] += PetscAbsScalar(*b_val * *b_val);
2272:           b_val++;
2273:         }
2274:       }
2275:     }
2276:   } else if (type == NORM_1) {
2277:     for (i = a_aij->i[0]; i < a_aij->i[aij->A->rmap->n / bs]; i++) {
2278:       for (jb = 0; jb < bs; jb++) {
2279:         for (ib = 0; ib < bs; ib++) {
2280:           work[A->cmap->rstart + a_aij->j[i] * bs + jb] += PetscAbsScalar(*a_val);
2281:           a_val++;
2282:         }
2283:       }
2284:     }
2285:     for (i = b_aij->i[0]; i < b_aij->i[aij->B->rmap->n / bs]; i++) {
2286:       for (jb = 0; jb < bs; jb++) {
2287:         for (ib = 0; ib < bs; ib++) {
2288:           work[garray[b_aij->j[i]] * bs + jb] += PetscAbsScalar(*b_val);
2289:           b_val++;
2290:         }
2291:       }
2292:     }
2293:   } else if (type == NORM_INFINITY) {
2294:     for (i = a_aij->i[0]; i < a_aij->i[aij->A->rmap->n / bs]; i++) {
2295:       for (jb = 0; jb < bs; jb++) {
2296:         for (ib = 0; ib < bs; ib++) {
2297:           PetscInt col = A->cmap->rstart + a_aij->j[i] * bs + jb;
2298:           work[col]    = PetscMax(PetscAbsScalar(*a_val), work[col]);
2299:           a_val++;
2300:         }
2301:       }
2302:     }
2303:     for (i = b_aij->i[0]; i < b_aij->i[aij->B->rmap->n / bs]; i++) {
2304:       for (jb = 0; jb < bs; jb++) {
2305:         for (ib = 0; ib < bs; ib++) {
2306:           PetscInt col = garray[b_aij->j[i]] * bs + jb;
2307:           work[col]    = PetscMax(PetscAbsScalar(*b_val), work[col]);
2308:           b_val++;
2309:         }
2310:       }
2311:     }
2312:   } else if (type == REDUCTION_SUM_REALPART || type == REDUCTION_MEAN_REALPART) {
2313:     for (i = a_aij->i[0]; i < a_aij->i[aij->A->rmap->n / bs]; i++) {
2314:       for (jb = 0; jb < bs; jb++) {
2315:         for (ib = 0; ib < bs; ib++) {
2316:           work[A->cmap->rstart + a_aij->j[i] * bs + jb] += PetscRealPart(*a_val);
2317:           a_val++;
2318:         }
2319:       }
2320:     }
2321:     for (i = b_aij->i[0]; i < b_aij->i[aij->B->rmap->n / bs]; i++) {
2322:       for (jb = 0; jb < bs; jb++) {
2323:         for (ib = 0; ib < bs; ib++) {
2324:           work[garray[b_aij->j[i]] * bs + jb] += PetscRealPart(*b_val);
2325:           b_val++;
2326:         }
2327:       }
2328:     }
2329:   } else if (type == REDUCTION_SUM_IMAGINARYPART || type == REDUCTION_MEAN_IMAGINARYPART) {
2330:     for (i = a_aij->i[0]; i < a_aij->i[aij->A->rmap->n / bs]; i++) {
2331:       for (jb = 0; jb < bs; jb++) {
2332:         for (ib = 0; ib < bs; ib++) {
2333:           work[A->cmap->rstart + a_aij->j[i] * bs + jb] += PetscImaginaryPart(*a_val);
2334:           a_val++;
2335:         }
2336:       }
2337:     }
2338:     for (i = b_aij->i[0]; i < b_aij->i[aij->B->rmap->n / bs]; i++) {
2339:       for (jb = 0; jb < bs; jb++) {
2340:         for (ib = 0; ib < bs; ib++) {
2341:           work[garray[b_aij->j[i]] * bs + jb] += PetscImaginaryPart(*b_val);
2342:           b_val++;
2343:         }
2344:       }
2345:     }
2346:   } else SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Unknown reduction type");
2347:   if (type == NORM_INFINITY) {
2348:     PetscCallMPI(MPIU_Allreduce(work, reductions, N, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)A)));
2349:   } else {
2350:     PetscCallMPI(MPIU_Allreduce(work, reductions, N, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)A)));
2351:   }
2352:   PetscCall(PetscFree(work));
2353:   if (type == NORM_2) {
2354:     for (i = 0; i < N; i++) reductions[i] = PetscSqrtReal(reductions[i]);
2355:   } else if (type == REDUCTION_MEAN_REALPART || type == REDUCTION_MEAN_IMAGINARYPART) {
2356:     for (i = 0; i < N; i++) reductions[i] /= m;
2357:   }
2358:   PetscFunctionReturn(PETSC_SUCCESS);
2359: }

2361: static PetscErrorCode MatInvertBlockDiagonal_MPIBAIJ(Mat A, const PetscScalar **values)
2362: {
2363:   Mat_MPIBAIJ *a = (Mat_MPIBAIJ *)A->data;

2365:   PetscFunctionBegin;
2366:   PetscCall(MatInvertBlockDiagonal(a->A, values));
2367:   A->factorerrortype             = a->A->factorerrortype;
2368:   A->factorerror_zeropivot_value = a->A->factorerror_zeropivot_value;
2369:   A->factorerror_zeropivot_row   = a->A->factorerror_zeropivot_row;
2370:   PetscFunctionReturn(PETSC_SUCCESS);
2371: }

2373: static PetscErrorCode MatShift_MPIBAIJ(Mat Y, PetscScalar a)
2374: {
2375:   Mat_MPIBAIJ *maij = (Mat_MPIBAIJ *)Y->data;
2376:   Mat_SeqBAIJ *aij  = (Mat_SeqBAIJ *)maij->A->data;

2378:   PetscFunctionBegin;
2379:   if (!Y->preallocated) {
2380:     PetscCall(MatMPIBAIJSetPreallocation(Y, Y->rmap->bs, 1, NULL, 0, NULL));
2381:   } else if (!aij->nz) {
2382:     PetscInt nonew = aij->nonew;
2383:     PetscCall(MatSeqBAIJSetPreallocation(maij->A, Y->rmap->bs, 1, NULL));
2384:     aij->nonew = nonew;
2385:   }
2386:   PetscCall(MatShift_Basic(Y, a));
2387:   PetscFunctionReturn(PETSC_SUCCESS);
2388: }

2390: static PetscErrorCode MatMissingDiagonal_MPIBAIJ(Mat A, PetscBool *missing, PetscInt *d)
2391: {
2392:   Mat_MPIBAIJ *a = (Mat_MPIBAIJ *)A->data;

2394:   PetscFunctionBegin;
2395:   PetscCheck(A->rmap->n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only works for square matrices");
2396:   PetscCall(MatMissingDiagonal(a->A, missing, d));
2397:   if (d) {
2398:     PetscInt rstart;
2399:     PetscCall(MatGetOwnershipRange(A, &rstart, NULL));
2400:     *d += rstart / A->rmap->bs;
2401:   }
2402:   PetscFunctionReturn(PETSC_SUCCESS);
2403: }

2405: static PetscErrorCode MatGetDiagonalBlock_MPIBAIJ(Mat A, Mat *a)
2406: {
2407:   PetscFunctionBegin;
2408:   *a = ((Mat_MPIBAIJ *)A->data)->A;
2409:   PetscFunctionReturn(PETSC_SUCCESS);
2410: }

2412: static PetscErrorCode MatEliminateZeros_MPIBAIJ(Mat A, PetscBool keep)
2413: {
2414:   Mat_MPIBAIJ *a = (Mat_MPIBAIJ *)A->data;

2416:   PetscFunctionBegin;
2417:   PetscCall(MatEliminateZeros_SeqBAIJ(a->A, keep));        // possibly keep zero diagonal coefficients
2418:   PetscCall(MatEliminateZeros_SeqBAIJ(a->B, PETSC_FALSE)); // never keep zero diagonal coefficients
2419:   PetscFunctionReturn(PETSC_SUCCESS);
2420: }

2422: static struct _MatOps MatOps_Values = {MatSetValues_MPIBAIJ,
2423:                                        MatGetRow_MPIBAIJ,
2424:                                        MatRestoreRow_MPIBAIJ,
2425:                                        MatMult_MPIBAIJ,
2426:                                        /* 4*/ MatMultAdd_MPIBAIJ,
2427:                                        MatMultTranspose_MPIBAIJ,
2428:                                        MatMultTransposeAdd_MPIBAIJ,
2429:                                        NULL,
2430:                                        NULL,
2431:                                        NULL,
2432:                                        /*10*/ NULL,
2433:                                        NULL,
2434:                                        NULL,
2435:                                        MatSOR_MPIBAIJ,
2436:                                        MatTranspose_MPIBAIJ,
2437:                                        /*15*/ MatGetInfo_MPIBAIJ,
2438:                                        MatEqual_MPIBAIJ,
2439:                                        MatGetDiagonal_MPIBAIJ,
2440:                                        MatDiagonalScale_MPIBAIJ,
2441:                                        MatNorm_MPIBAIJ,
2442:                                        /*20*/ MatAssemblyBegin_MPIBAIJ,
2443:                                        MatAssemblyEnd_MPIBAIJ,
2444:                                        MatSetOption_MPIBAIJ,
2445:                                        MatZeroEntries_MPIBAIJ,
2446:                                        /*24*/ MatZeroRows_MPIBAIJ,
2447:                                        NULL,
2448:                                        NULL,
2449:                                        NULL,
2450:                                        NULL,
2451:                                        /*29*/ MatSetUp_MPI_Hash,
2452:                                        NULL,
2453:                                        NULL,
2454:                                        MatGetDiagonalBlock_MPIBAIJ,
2455:                                        NULL,
2456:                                        /*34*/ MatDuplicate_MPIBAIJ,
2457:                                        NULL,
2458:                                        NULL,
2459:                                        NULL,
2460:                                        NULL,
2461:                                        /*39*/ MatAXPY_MPIBAIJ,
2462:                                        MatCreateSubMatrices_MPIBAIJ,
2463:                                        MatIncreaseOverlap_MPIBAIJ,
2464:                                        MatGetValues_MPIBAIJ,
2465:                                        MatCopy_MPIBAIJ,
2466:                                        /*44*/ NULL,
2467:                                        MatScale_MPIBAIJ,
2468:                                        MatShift_MPIBAIJ,
2469:                                        NULL,
2470:                                        MatZeroRowsColumns_MPIBAIJ,
2471:                                        /*49*/ NULL,
2472:                                        NULL,
2473:                                        NULL,
2474:                                        NULL,
2475:                                        NULL,
2476:                                        /*54*/ MatFDColoringCreate_MPIXAIJ,
2477:                                        NULL,
2478:                                        MatSetUnfactored_MPIBAIJ,
2479:                                        MatPermute_MPIBAIJ,
2480:                                        MatSetValuesBlocked_MPIBAIJ,
2481:                                        /*59*/ MatCreateSubMatrix_MPIBAIJ,
2482:                                        MatDestroy_MPIBAIJ,
2483:                                        MatView_MPIBAIJ,
2484:                                        NULL,
2485:                                        NULL,
2486:                                        /*64*/ NULL,
2487:                                        NULL,
2488:                                        NULL,
2489:                                        NULL,
2490:                                        NULL,
2491:                                        /*69*/ MatGetRowMaxAbs_MPIBAIJ,
2492:                                        NULL,
2493:                                        NULL,
2494:                                        NULL,
2495:                                        NULL,
2496:                                        /*74*/ NULL,
2497:                                        MatFDColoringApply_BAIJ,
2498:                                        NULL,
2499:                                        NULL,
2500:                                        NULL,
2501:                                        /*79*/ NULL,
2502:                                        NULL,
2503:                                        NULL,
2504:                                        NULL,
2505:                                        MatLoad_MPIBAIJ,
2506:                                        /*84*/ NULL,
2507:                                        NULL,
2508:                                        NULL,
2509:                                        NULL,
2510:                                        NULL,
2511:                                        /*89*/ NULL,
2512:                                        NULL,
2513:                                        NULL,
2514:                                        NULL,
2515:                                        NULL,
2516:                                        /*94*/ NULL,
2517:                                        NULL,
2518:                                        NULL,
2519:                                        NULL,
2520:                                        NULL,
2521:                                        /*99*/ NULL,
2522:                                        NULL,
2523:                                        NULL,
2524:                                        MatConjugate_MPIBAIJ,
2525:                                        NULL,
2526:                                        /*104*/ NULL,
2527:                                        MatRealPart_MPIBAIJ,
2528:                                        MatImaginaryPart_MPIBAIJ,
2529:                                        NULL,
2530:                                        NULL,
2531:                                        /*109*/ NULL,
2532:                                        NULL,
2533:                                        NULL,
2534:                                        NULL,
2535:                                        MatMissingDiagonal_MPIBAIJ,
2536:                                        /*114*/ MatGetSeqNonzeroStructure_MPIBAIJ,
2537:                                        NULL,
2538:                                        MatGetGhosts_MPIBAIJ,
2539:                                        NULL,
2540:                                        NULL,
2541:                                        /*119*/ NULL,
2542:                                        NULL,
2543:                                        NULL,
2544:                                        NULL,
2545:                                        MatGetMultiProcBlock_MPIBAIJ,
2546:                                        /*124*/ NULL,
2547:                                        MatGetColumnReductions_MPIBAIJ,
2548:                                        MatInvertBlockDiagonal_MPIBAIJ,
2549:                                        NULL,
2550:                                        NULL,
2551:                                        /*129*/ NULL,
2552:                                        NULL,
2553:                                        NULL,
2554:                                        NULL,
2555:                                        NULL,
2556:                                        /*134*/ NULL,
2557:                                        NULL,
2558:                                        NULL,
2559:                                        NULL,
2560:                                        NULL,
2561:                                        /*139*/ MatSetBlockSizes_Default,
2562:                                        NULL,
2563:                                        NULL,
2564:                                        MatFDColoringSetUp_MPIXAIJ,
2565:                                        NULL,
2566:                                        /*144*/ MatCreateMPIMatConcatenateSeqMat_MPIBAIJ,
2567:                                        NULL,
2568:                                        NULL,
2569:                                        NULL,
2570:                                        NULL,
2571:                                        NULL,
2572:                                        /*150*/ NULL,
2573:                                        MatEliminateZeros_MPIBAIJ,
2574:                                        MatGetRowSumAbs_MPIBAIJ,
2575:                                        NULL,
2576:                                        NULL,
2577:                                        /*155*/ NULL,
2578:                                        MatCopyHashToXAIJ_MPI_Hash};

2580: PETSC_INTERN PetscErrorCode MatConvert_MPIBAIJ_MPISBAIJ(Mat, MatType, MatReuse, Mat *);
2581: PETSC_INTERN PetscErrorCode MatConvert_XAIJ_IS(Mat, MatType, MatReuse, Mat *);

2583: static PetscErrorCode MatMPIBAIJSetPreallocationCSR_MPIBAIJ(Mat B, PetscInt bs, const PetscInt ii[], const PetscInt jj[], const PetscScalar V[])
2584: {
2585:   PetscInt        m, rstart, cstart, cend;
2586:   PetscInt        i, j, dlen, olen, nz, nz_max = 0, *d_nnz = NULL, *o_nnz = NULL;
2587:   const PetscInt *JJ          = NULL;
2588:   PetscScalar    *values      = NULL;
2589:   PetscBool       roworiented = ((Mat_MPIBAIJ *)B->data)->roworiented;
2590:   PetscBool       nooffprocentries;

2592:   PetscFunctionBegin;
2593:   PetscCall(PetscLayoutSetBlockSize(B->rmap, bs));
2594:   PetscCall(PetscLayoutSetBlockSize(B->cmap, bs));
2595:   PetscCall(PetscLayoutSetUp(B->rmap));
2596:   PetscCall(PetscLayoutSetUp(B->cmap));
2597:   PetscCall(PetscLayoutGetBlockSize(B->rmap, &bs));
2598:   m      = B->rmap->n / bs;
2599:   rstart = B->rmap->rstart / bs;
2600:   cstart = B->cmap->rstart / bs;
2601:   cend   = B->cmap->rend / bs;

2603:   PetscCheck(!ii[0], PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "ii[0] must be 0 but it is %" PetscInt_FMT, ii[0]);
2604:   PetscCall(PetscMalloc2(m, &d_nnz, m, &o_nnz));
2605:   for (i = 0; i < m; i++) {
2606:     nz = ii[i + 1] - ii[i];
2607:     PetscCheck(nz >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local row %" PetscInt_FMT " has a negative number of columns %" PetscInt_FMT, i, nz);
2608:     nz_max = PetscMax(nz_max, nz);
2609:     dlen   = 0;
2610:     olen   = 0;
2611:     JJ     = jj + ii[i];
2612:     for (j = 0; j < nz; j++) {
2613:       if (*JJ < cstart || *JJ >= cend) olen++;
2614:       else dlen++;
2615:       JJ++;
2616:     }
2617:     d_nnz[i] = dlen;
2618:     o_nnz[i] = olen;
2619:   }
2620:   PetscCall(MatMPIBAIJSetPreallocation(B, bs, 0, d_nnz, 0, o_nnz));
2621:   PetscCall(PetscFree2(d_nnz, o_nnz));

2623:   values = (PetscScalar *)V;
2624:   if (!values) PetscCall(PetscCalloc1(bs * bs * nz_max, &values));
2625:   for (i = 0; i < m; i++) {
2626:     PetscInt        row   = i + rstart;
2627:     PetscInt        ncols = ii[i + 1] - ii[i];
2628:     const PetscInt *icols = jj + ii[i];
2629:     if (bs == 1 || !roworiented) { /* block ordering matches the non-nested layout of MatSetValues so we can insert entire rows */
2630:       const PetscScalar *svals = values + (V ? (bs * bs * ii[i]) : 0);
2631:       PetscCall(MatSetValuesBlocked_MPIBAIJ(B, 1, &row, ncols, icols, svals, INSERT_VALUES));
2632:     } else { /* block ordering does not match so we can only insert one block at a time. */
2633:       PetscInt j;
2634:       for (j = 0; j < ncols; j++) {
2635:         const PetscScalar *svals = values + (V ? (bs * bs * (ii[i] + j)) : 0);
2636:         PetscCall(MatSetValuesBlocked_MPIBAIJ(B, 1, &row, 1, &icols[j], svals, INSERT_VALUES));
2637:       }
2638:     }
2639:   }

2641:   if (!V) PetscCall(PetscFree(values));
2642:   nooffprocentries    = B->nooffprocentries;
2643:   B->nooffprocentries = PETSC_TRUE;
2644:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
2645:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
2646:   B->nooffprocentries = nooffprocentries;

2648:   PetscCall(MatSetOption(B, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
2649:   PetscFunctionReturn(PETSC_SUCCESS);
2650: }

2652: /*@C
2653:   MatMPIBAIJSetPreallocationCSR - Creates a sparse parallel matrix in `MATBAIJ` format using the given nonzero structure and (optional) numerical values

2655:   Collective

2657:   Input Parameters:
2658: + B  - the matrix
2659: . bs - the block size
2660: . i  - the indices into `j` for the start of each local row (starts with zero)
2661: . j  - the column indices for each local row (starts with zero) these must be sorted for each row
2662: - v  - optional values in the matrix, use `NULL` if not provided

2664:   Level: advanced

2666:   Notes:
2667:   The `i`, `j`, and `v` arrays ARE copied by this routine into the internal format used by PETSc;
2668:   thus you CANNOT change the matrix entries by changing the values of `v` after you have
2669:   called this routine.

2671:   The order of the entries in values is specified by the `MatOption` `MAT_ROW_ORIENTED`.  For example, C programs
2672:   may want to use the default `MAT_ROW_ORIENTED` with value `PETSC_TRUE` and use an array v[nnz][bs][bs] where the second index is
2673:   over rows within a block and the last index is over columns within a block row.  Fortran programs will likely set
2674:   `MAT_ROW_ORIENTED` with value `PETSC_FALSE` and use a Fortran array v(bs,bs,nnz) in which the first index is over rows within a
2675:   block column and the second index is over columns within a block.

2677:   Though this routine has Preallocation() in the name it also sets the exact nonzero locations of the matrix entries and usually the numerical values as well

2679: .seealso: `Mat`, `MatCreate()`, `MatCreateSeqAIJ()`, `MatSetValues()`, `MatMPIBAIJSetPreallocation()`, `MatCreateAIJ()`, `MATMPIAIJ`, `MatCreateMPIBAIJWithArrays()`, `MATMPIBAIJ`
2680: @*/
2681: PetscErrorCode MatMPIBAIJSetPreallocationCSR(Mat B, PetscInt bs, const PetscInt i[], const PetscInt j[], const PetscScalar v[])
2682: {
2683:   PetscFunctionBegin;
2687:   PetscTryMethod(B, "MatMPIBAIJSetPreallocationCSR_C", (Mat, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[]), (B, bs, i, j, v));
2688:   PetscFunctionReturn(PETSC_SUCCESS);
2689: }

2691: PetscErrorCode MatMPIBAIJSetPreallocation_MPIBAIJ(Mat B, PetscInt bs, PetscInt d_nz, const PetscInt *d_nnz, PetscInt o_nz, const PetscInt *o_nnz)
2692: {
2693:   Mat_MPIBAIJ *b = (Mat_MPIBAIJ *)B->data;
2694:   PetscInt     i;
2695:   PetscMPIInt  size;

2697:   PetscFunctionBegin;
2698:   if (B->hash_active) {
2699:     B->ops[0]      = b->cops;
2700:     B->hash_active = PETSC_FALSE;
2701:   }
2702:   if (!B->preallocated) PetscCall(MatStashCreate_Private(PetscObjectComm((PetscObject)B), bs, &B->bstash));
2703:   PetscCall(MatSetBlockSize(B, PetscAbs(bs)));
2704:   PetscCall(PetscLayoutSetUp(B->rmap));
2705:   PetscCall(PetscLayoutSetUp(B->cmap));
2706:   PetscCall(PetscLayoutGetBlockSize(B->rmap, &bs));

2708:   if (d_nnz) {
2709:     for (i = 0; i < B->rmap->n / bs; i++) PetscCheck(d_nnz[i] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "d_nnz cannot be less than -1: local row %" PetscInt_FMT " value %" PetscInt_FMT, i, d_nnz[i]);
2710:   }
2711:   if (o_nnz) {
2712:     for (i = 0; i < B->rmap->n / bs; i++) PetscCheck(o_nnz[i] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "o_nnz cannot be less than -1: local row %" PetscInt_FMT " value %" PetscInt_FMT, i, o_nnz[i]);
2713:   }

2715:   b->bs2 = bs * bs;
2716:   b->mbs = B->rmap->n / bs;
2717:   b->nbs = B->cmap->n / bs;
2718:   b->Mbs = B->rmap->N / bs;
2719:   b->Nbs = B->cmap->N / bs;

2721:   for (i = 0; i <= b->size; i++) b->rangebs[i] = B->rmap->range[i] / bs;
2722:   b->rstartbs = B->rmap->rstart / bs;
2723:   b->rendbs   = B->rmap->rend / bs;
2724:   b->cstartbs = B->cmap->rstart / bs;
2725:   b->cendbs   = B->cmap->rend / bs;

2727: #if defined(PETSC_USE_CTABLE)
2728:   PetscCall(PetscHMapIDestroy(&b->colmap));
2729: #else
2730:   PetscCall(PetscFree(b->colmap));
2731: #endif
2732:   PetscCall(PetscFree(b->garray));
2733:   PetscCall(VecDestroy(&b->lvec));
2734:   PetscCall(VecScatterDestroy(&b->Mvctx));

2736:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)B), &size));

2738:   MatSeqXAIJGetOptions_Private(b->B);
2739:   PetscCall(MatDestroy(&b->B));
2740:   PetscCall(MatCreate(PETSC_COMM_SELF, &b->B));
2741:   PetscCall(MatSetSizes(b->B, B->rmap->n, size > 1 ? B->cmap->N : 0, B->rmap->n, size > 1 ? B->cmap->N : 0));
2742:   PetscCall(MatSetType(b->B, MATSEQBAIJ));
2743:   MatSeqXAIJRestoreOptions_Private(b->B);

2745:   MatSeqXAIJGetOptions_Private(b->A);
2746:   PetscCall(MatDestroy(&b->A));
2747:   PetscCall(MatCreate(PETSC_COMM_SELF, &b->A));
2748:   PetscCall(MatSetSizes(b->A, B->rmap->n, B->cmap->n, B->rmap->n, B->cmap->n));
2749:   PetscCall(MatSetType(b->A, MATSEQBAIJ));
2750:   MatSeqXAIJRestoreOptions_Private(b->A);

2752:   PetscCall(MatSeqBAIJSetPreallocation(b->A, bs, d_nz, d_nnz));
2753:   PetscCall(MatSeqBAIJSetPreallocation(b->B, bs, o_nz, o_nnz));
2754:   B->preallocated  = PETSC_TRUE;
2755:   B->was_assembled = PETSC_FALSE;
2756:   B->assembled     = PETSC_FALSE;
2757:   PetscFunctionReturn(PETSC_SUCCESS);
2758: }

2760: extern PetscErrorCode MatDiagonalScaleLocal_MPIBAIJ(Mat, Vec);
2761: extern PetscErrorCode MatSetHashTableFactor_MPIBAIJ(Mat, PetscReal);

2763: PETSC_INTERN PetscErrorCode MatConvert_MPIBAIJ_MPIAdj(Mat B, MatType newtype, MatReuse reuse, Mat *adj)
2764: {
2765:   Mat_MPIBAIJ    *b = (Mat_MPIBAIJ *)B->data;
2766:   Mat_SeqBAIJ    *d = (Mat_SeqBAIJ *)b->A->data, *o = (Mat_SeqBAIJ *)b->B->data;
2767:   PetscInt        M = B->rmap->n / B->rmap->bs, i, *ii, *jj, cnt, j, k, rstart = B->rmap->rstart / B->rmap->bs;
2768:   const PetscInt *id = d->i, *jd = d->j, *io = o->i, *jo = o->j, *garray = b->garray;

2770:   PetscFunctionBegin;
2771:   PetscCall(PetscMalloc1(M + 1, &ii));
2772:   ii[0] = 0;
2773:   for (i = 0; i < M; i++) {
2774:     PetscCheck((id[i + 1] - id[i]) >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Indices wrong %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT, i, id[i], id[i + 1]);
2775:     PetscCheck((io[i + 1] - io[i]) >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Indices wrong %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT, i, io[i], io[i + 1]);
2776:     ii[i + 1] = ii[i] + id[i + 1] - id[i] + io[i + 1] - io[i];
2777:     /* remove one from count of matrix has diagonal */
2778:     for (j = id[i]; j < id[i + 1]; j++) {
2779:       if (jd[j] == i) {
2780:         ii[i + 1]--;
2781:         break;
2782:       }
2783:     }
2784:   }
2785:   PetscCall(PetscMalloc1(ii[M], &jj));
2786:   cnt = 0;
2787:   for (i = 0; i < M; i++) {
2788:     for (j = io[i]; j < io[i + 1]; j++) {
2789:       if (garray[jo[j]] > rstart) break;
2790:       jj[cnt++] = garray[jo[j]];
2791:     }
2792:     for (k = id[i]; k < id[i + 1]; k++) {
2793:       if (jd[k] != i) jj[cnt++] = rstart + jd[k];
2794:     }
2795:     for (; j < io[i + 1]; j++) jj[cnt++] = garray[jo[j]];
2796:   }
2797:   PetscCall(MatCreateMPIAdj(PetscObjectComm((PetscObject)B), M, B->cmap->N / B->rmap->bs, ii, jj, NULL, adj));
2798:   PetscFunctionReturn(PETSC_SUCCESS);
2799: }

2801: #include <../src/mat/impls/aij/mpi/mpiaij.h>

2803: PETSC_INTERN PetscErrorCode MatConvert_SeqBAIJ_SeqAIJ(Mat, MatType, MatReuse, Mat *);

2805: PETSC_INTERN PetscErrorCode MatConvert_MPIBAIJ_MPIAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
2806: {
2807:   Mat_MPIBAIJ *a = (Mat_MPIBAIJ *)A->data;
2808:   Mat_MPIAIJ  *b;
2809:   Mat          B;

2811:   PetscFunctionBegin;
2812:   PetscCheck(A->assembled, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Matrix must be assembled");

2814:   if (reuse == MAT_REUSE_MATRIX) {
2815:     B = *newmat;
2816:   } else {
2817:     PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
2818:     PetscCall(MatSetType(B, MATMPIAIJ));
2819:     PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
2820:     PetscCall(MatSetBlockSizes(B, A->rmap->bs, A->cmap->bs));
2821:     PetscCall(MatSeqAIJSetPreallocation(B, 0, NULL));
2822:     PetscCall(MatMPIAIJSetPreallocation(B, 0, NULL, 0, NULL));
2823:   }
2824:   b = (Mat_MPIAIJ *)B->data;

2826:   if (reuse == MAT_REUSE_MATRIX) {
2827:     PetscCall(MatConvert_SeqBAIJ_SeqAIJ(a->A, MATSEQAIJ, MAT_REUSE_MATRIX, &b->A));
2828:     PetscCall(MatConvert_SeqBAIJ_SeqAIJ(a->B, MATSEQAIJ, MAT_REUSE_MATRIX, &b->B));
2829:   } else {
2830:     PetscInt   *garray = a->garray;
2831:     Mat_SeqAIJ *bB;
2832:     PetscInt    bs, nnz;
2833:     PetscCall(MatDestroy(&b->A));
2834:     PetscCall(MatDestroy(&b->B));
2835:     /* just clear out the data structure */
2836:     PetscCall(MatDisAssemble_MPIAIJ(B, PETSC_FALSE));
2837:     PetscCall(MatConvert_SeqBAIJ_SeqAIJ(a->A, MATSEQAIJ, MAT_INITIAL_MATRIX, &b->A));
2838:     PetscCall(MatConvert_SeqBAIJ_SeqAIJ(a->B, MATSEQAIJ, MAT_INITIAL_MATRIX, &b->B));

2840:     /* Global numbering for b->B columns */
2841:     bB  = (Mat_SeqAIJ *)b->B->data;
2842:     bs  = A->rmap->bs;
2843:     nnz = bB->i[A->rmap->n];
2844:     for (PetscInt k = 0; k < nnz; k++) {
2845:       PetscInt bj = bB->j[k] / bs;
2846:       PetscInt br = bB->j[k] % bs;
2847:       bB->j[k]    = garray[bj] * bs + br;
2848:     }
2849:   }
2850:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
2851:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));

2853:   if (reuse == MAT_INPLACE_MATRIX) {
2854:     PetscCall(MatHeaderReplace(A, &B));
2855:   } else {
2856:     *newmat = B;
2857:   }
2858:   PetscFunctionReturn(PETSC_SUCCESS);
2859: }

2861: /*MC
2862:    MATMPIBAIJ - MATMPIBAIJ = "mpibaij" - A matrix type to be used for distributed block sparse matrices.

2864:    Options Database Keys:
2865: + -mat_type mpibaij - sets the matrix type to `MATMPIBAIJ` during a call to `MatSetFromOptions()`
2866: . -mat_block_size <bs> - set the blocksize used to store the matrix
2867: . -mat_baij_mult_version version - indicate the version of the matrix-vector product to use  (0 often indicates using BLAS)
2868: - -mat_use_hash_table <fact> - set hash table factor

2870:    Level: beginner

2872:    Note:
2873:     `MatSetOption(A, MAT_STRUCTURE_ONLY, PETSC_TRUE)` may be called for this matrix type. In this no
2874:     space is allocated for the nonzero entries and any entries passed with `MatSetValues()` are ignored

2876: .seealso: `Mat`, `MATBAIJ`, `MATSEQBAIJ`, `MatCreateBAIJ`
2877: M*/

2879: PETSC_INTERN PetscErrorCode MatConvert_MPIBAIJ_MPIBSTRM(Mat, MatType, MatReuse, Mat *);

2881: PETSC_EXTERN PetscErrorCode MatCreate_MPIBAIJ(Mat B)
2882: {
2883:   Mat_MPIBAIJ *b;
2884:   PetscBool    flg = PETSC_FALSE;

2886:   PetscFunctionBegin;
2887:   PetscCall(PetscNew(&b));
2888:   B->data      = (void *)b;
2889:   B->ops[0]    = MatOps_Values;
2890:   B->assembled = PETSC_FALSE;

2892:   B->insertmode = NOT_SET_VALUES;
2893:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)B), &b->rank));
2894:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)B), &b->size));

2896:   /* build local table of row and column ownerships */
2897:   PetscCall(PetscMalloc1(b->size + 1, &b->rangebs));

2899:   /* build cache for off array entries formed */
2900:   PetscCall(MatStashCreate_Private(PetscObjectComm((PetscObject)B), 1, &B->stash));

2902:   b->donotstash  = PETSC_FALSE;
2903:   b->colmap      = NULL;
2904:   b->garray      = NULL;
2905:   b->roworiented = PETSC_TRUE;

2907:   /* stuff used in block assembly */
2908:   b->barray = NULL;

2910:   /* stuff used for matrix vector multiply */
2911:   b->lvec  = NULL;
2912:   b->Mvctx = NULL;

2914:   /* stuff for MatGetRow() */
2915:   b->rowindices   = NULL;
2916:   b->rowvalues    = NULL;
2917:   b->getrowactive = PETSC_FALSE;

2919:   /* hash table stuff */
2920:   b->ht           = NULL;
2921:   b->hd           = NULL;
2922:   b->ht_size      = 0;
2923:   b->ht_flag      = PETSC_FALSE;
2924:   b->ht_fact      = 0;
2925:   b->ht_total_ct  = 0;
2926:   b->ht_insert_ct = 0;

2928:   /* stuff for MatCreateSubMatrices_MPIBAIJ_local() */
2929:   b->ijonly = PETSC_FALSE;

2931:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_mpibaij_mpiadj_C", MatConvert_MPIBAIJ_MPIAdj));
2932:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_mpibaij_mpiaij_C", MatConvert_MPIBAIJ_MPIAIJ));
2933:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_mpibaij_mpisbaij_C", MatConvert_MPIBAIJ_MPISBAIJ));
2934: #if defined(PETSC_HAVE_HYPRE)
2935:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_mpibaij_hypre_C", MatConvert_AIJ_HYPRE));
2936: #endif
2937:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatStoreValues_C", MatStoreValues_MPIBAIJ));
2938:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatRetrieveValues_C", MatRetrieveValues_MPIBAIJ));
2939:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatMPIBAIJSetPreallocation_C", MatMPIBAIJSetPreallocation_MPIBAIJ));
2940:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatMPIBAIJSetPreallocationCSR_C", MatMPIBAIJSetPreallocationCSR_MPIBAIJ));
2941:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDiagonalScaleLocal_C", MatDiagonalScaleLocal_MPIBAIJ));
2942:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSetHashTableFactor_C", MatSetHashTableFactor_MPIBAIJ));
2943:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_mpibaij_is_C", MatConvert_XAIJ_IS));
2944:   PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATMPIBAIJ));

2946:   PetscOptionsBegin(PetscObjectComm((PetscObject)B), NULL, "Options for loading MPIBAIJ matrix 1", "Mat");
2947:   PetscCall(PetscOptionsName("-mat_use_hash_table", "Use hash table to save time in constructing matrix", "MatSetOption", &flg));
2948:   if (flg) {
2949:     PetscReal fact = 1.39;
2950:     PetscCall(MatSetOption(B, MAT_USE_HASH_TABLE, PETSC_TRUE));
2951:     PetscCall(PetscOptionsReal("-mat_use_hash_table", "Use hash table factor", "MatMPIBAIJSetHashTableFactor", fact, &fact, NULL));
2952:     if (fact <= 1.0) fact = 1.39;
2953:     PetscCall(MatMPIBAIJSetHashTableFactor(B, fact));
2954:     PetscCall(PetscInfo(B, "Hash table Factor used %5.2g\n", (double)fact));
2955:   }
2956:   PetscOptionsEnd();
2957:   PetscFunctionReturn(PETSC_SUCCESS);
2958: }

2960: // PetscClangLinter pragma disable: -fdoc-section-header-unknown
2961: /*MC
2962:    MATBAIJ - MATBAIJ = "baij" - A matrix type to be used for block sparse matrices.

2964:    This matrix type is identical to `MATSEQBAIJ` when constructed with a single process communicator,
2965:    and `MATMPIBAIJ` otherwise.

2967:    Options Database Keys:
2968: . -mat_type baij - sets the matrix type to `MATBAIJ` during a call to `MatSetFromOptions()`

2970:   Level: beginner

2972: .seealso: `Mat`, `MatCreateBAIJ()`, `MATSEQBAIJ`, `MATMPIBAIJ`, `MatMPIBAIJSetPreallocation()`, `MatMPIBAIJSetPreallocationCSR()`
2973: M*/

2975: /*@
2976:   MatMPIBAIJSetPreallocation - Allocates memory for a sparse parallel matrix in `MATMPIBAIJ` format
2977:   (block compressed row).

2979:   Collective

2981:   Input Parameters:
2982: + B     - the matrix
2983: . bs    - size of block, the blocks are ALWAYS square. One can use `MatSetBlockSizes()` to set a different row and column blocksize but the row
2984:           blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with `MatCreateVecs()`
2985: . d_nz  - number of block nonzeros per block row in diagonal portion of local
2986:            submatrix  (same for all local rows)
2987: . d_nnz - array containing the number of block nonzeros in the various block rows
2988:            of the in diagonal portion of the local (possibly different for each block
2989:            row) or `NULL`.  If you plan to factor the matrix you must leave room for the diagonal entry and
2990:            set it even if it is zero.
2991: . o_nz  - number of block nonzeros per block row in the off-diagonal portion of local
2992:            submatrix (same for all local rows).
2993: - o_nnz - array containing the number of nonzeros in the various block rows of the
2994:            off-diagonal portion of the local submatrix (possibly different for
2995:            each block row) or `NULL`.

2997:    If the *_nnz parameter is given then the *_nz parameter is ignored

2999:   Options Database Keys:
3000: + -mat_block_size            - size of the blocks to use
3001: - -mat_use_hash_table <fact> - set hash table factor

3003:   Level: intermediate

3005:   Notes:
3006:   For good matrix assembly performance
3007:   the user should preallocate the matrix storage by setting the parameters
3008:   `d_nz` (or `d_nnz`) and `o_nz` (or `o_nnz`).  By setting these parameters accurately,
3009:   performance can be increased by more than a factor of 50.

3011:   If `PETSC_DECIDE` or  `PETSC_DETERMINE` is used for a particular argument on one processor
3012:   than it must be used on all processors that share the object for that argument.

3014:   Storage Information:
3015:   For a square global matrix we define each processor's diagonal portion
3016:   to be its local rows and the corresponding columns (a square submatrix);
3017:   each processor's off-diagonal portion encompasses the remainder of the
3018:   local matrix (a rectangular submatrix).

3020:   The user can specify preallocated storage for the diagonal part of
3021:   the local submatrix with either `d_nz` or `d_nnz` (not both).  Set
3022:   `d_nz` = `PETSC_DEFAULT` and `d_nnz` = `NULL` for PETSc to control dynamic
3023:   memory allocation.  Likewise, specify preallocated storage for the
3024:   off-diagonal part of the local submatrix with `o_nz` or `o_nnz` (not both).

3026:   Consider a processor that owns rows 3, 4 and 5 of a parallel matrix. In
3027:   the figure below we depict these three local rows and all columns (0-11).

3029: .vb
3030:            0 1 2 3 4 5 6 7 8 9 10 11
3031:           --------------------------
3032:    row 3  |o o o d d d o o o o  o  o
3033:    row 4  |o o o d d d o o o o  o  o
3034:    row 5  |o o o d d d o o o o  o  o
3035:           --------------------------
3036: .ve

3038:   Thus, any entries in the d locations are stored in the d (diagonal)
3039:   submatrix, and any entries in the o locations are stored in the
3040:   o (off-diagonal) submatrix.  Note that the d and the o submatrices are
3041:   stored simply in the `MATSEQBAIJ` format for compressed row storage.

3043:   Now `d_nz` should indicate the number of block nonzeros per row in the d matrix,
3044:   and `o_nz` should indicate the number of block nonzeros per row in the o matrix.
3045:   In general, for PDE problems in which most nonzeros are near the diagonal,
3046:   one expects `d_nz` >> `o_nz`.

3048:   You can call `MatGetInfo()` to get information on how effective the preallocation was;
3049:   for example the fields mallocs,nz_allocated,nz_used,nz_unneeded;
3050:   You can also run with the option `-info` and look for messages with the string
3051:   malloc in them to see if additional memory allocation was needed.

3053: .seealso: `Mat`, `MATMPIBAIJ`, `MatCreate()`, `MatCreateSeqBAIJ()`, `MatSetValues()`, `MatCreateBAIJ()`, `MatMPIBAIJSetPreallocationCSR()`, `PetscSplitOwnership()`
3054: @*/
3055: PetscErrorCode MatMPIBAIJSetPreallocation(Mat B, PetscInt bs, PetscInt d_nz, const PetscInt d_nnz[], PetscInt o_nz, const PetscInt o_nnz[])
3056: {
3057:   PetscFunctionBegin;
3061:   PetscTryMethod(B, "MatMPIBAIJSetPreallocation_C", (Mat, PetscInt, PetscInt, const PetscInt[], PetscInt, const PetscInt[]), (B, bs, d_nz, d_nnz, o_nz, o_nnz));
3062:   PetscFunctionReturn(PETSC_SUCCESS);
3063: }

3065: // PetscClangLinter pragma disable: -fdoc-section-header-unknown
3066: /*@
3067:   MatCreateBAIJ - Creates a sparse parallel matrix in `MATBAIJ` format
3068:   (block compressed row).

3070:   Collective

3072:   Input Parameters:
3073: + comm  - MPI communicator
3074: . bs    - size of block, the blocks are ALWAYS square. One can use `MatSetBlockSizes()` to set a different row and column blocksize but the row
3075:           blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with `MatCreateVecs()`
3076: . m     - number of local rows (or `PETSC_DECIDE` to have calculated if M is given)
3077:           This value should be the same as the local size used in creating the
3078:           y vector for the matrix-vector product y = Ax.
3079: . n     - number of local columns (or `PETSC_DECIDE` to have calculated if N is given)
3080:           This value should be the same as the local size used in creating the
3081:           x vector for the matrix-vector product y = Ax.
3082: . M     - number of global rows (or `PETSC_DETERMINE` to have calculated if m is given)
3083: . N     - number of global columns (or `PETSC_DETERMINE` to have calculated if n is given)
3084: . d_nz  - number of nonzero blocks per block row in diagonal portion of local
3085:           submatrix  (same for all local rows)
3086: . d_nnz - array containing the number of nonzero blocks in the various block rows
3087:           of the in diagonal portion of the local (possibly different for each block
3088:           row) or NULL.  If you plan to factor the matrix you must leave room for the diagonal entry
3089:           and set it even if it is zero.
3090: . o_nz  - number of nonzero blocks per block row in the off-diagonal portion of local
3091:           submatrix (same for all local rows).
3092: - o_nnz - array containing the number of nonzero blocks in the various block rows of the
3093:           off-diagonal portion of the local submatrix (possibly different for
3094:           each block row) or NULL.

3096:   Output Parameter:
3097: . A - the matrix

3099:   Options Database Keys:
3100: + -mat_block_size            - size of the blocks to use
3101: - -mat_use_hash_table <fact> - set hash table factor

3103:   Level: intermediate

3105:   Notes:
3106:   It is recommended that one use `MatCreateFromOptions()` or the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()`,
3107:   MatXXXXSetPreallocation() paradigm instead of this routine directly.
3108:   [MatXXXXSetPreallocation() is, for example, `MatSeqBAIJSetPreallocation()`]

3110:   For good matrix assembly performance
3111:   the user should preallocate the matrix storage by setting the parameters
3112:   `d_nz` (or `d_nnz`) and `o_nz` (or `o_nnz`).  By setting these parameters accurately,
3113:   performance can be increased by more than a factor of 50.

3115:   If the *_nnz parameter is given then the *_nz parameter is ignored

3117:   A nonzero block is any block that as 1 or more nonzeros in it

3119:   The user MUST specify either the local or global matrix dimensions
3120:   (possibly both).

3122:   If `PETSC_DECIDE` or  `PETSC_DETERMINE` is used for a particular argument on one processor
3123:   than it must be used on all processors that share the object for that argument.

3125:   If `m` and `n` are not `PETSC_DECIDE`, then the values determine the `PetscLayout` of the matrix and the ranges returned by
3126:   `MatGetOwnershipRange()`,  `MatGetOwnershipRanges()`, `MatGetOwnershipRangeColumn()`, and `MatGetOwnershipRangesColumn()`.

3128:   Storage Information:
3129:   For a square global matrix we define each processor's diagonal portion
3130:   to be its local rows and the corresponding columns (a square submatrix);
3131:   each processor's off-diagonal portion encompasses the remainder of the
3132:   local matrix (a rectangular submatrix).

3134:   The user can specify preallocated storage for the diagonal part of
3135:   the local submatrix with either d_nz or d_nnz (not both).  Set
3136:   `d_nz` = `PETSC_DEFAULT` and `d_nnz` = `NULL` for PETSc to control dynamic
3137:   memory allocation.  Likewise, specify preallocated storage for the
3138:   off-diagonal part of the local submatrix with `o_nz` or `o_nnz` (not both).

3140:   Consider a processor that owns rows 3, 4 and 5 of a parallel matrix. In
3141:   the figure below we depict these three local rows and all columns (0-11).

3143: .vb
3144:            0 1 2 3 4 5 6 7 8 9 10 11
3145:           --------------------------
3146:    row 3  |o o o d d d o o o o  o  o
3147:    row 4  |o o o d d d o o o o  o  o
3148:    row 5  |o o o d d d o o o o  o  o
3149:           --------------------------
3150: .ve

3152:   Thus, any entries in the d locations are stored in the d (diagonal)
3153:   submatrix, and any entries in the o locations are stored in the
3154:   o (off-diagonal) submatrix.  Note that the d and the o submatrices are
3155:   stored simply in the `MATSEQBAIJ` format for compressed row storage.

3157:   Now `d_nz` should indicate the number of block nonzeros per row in the d matrix,
3158:   and `o_nz` should indicate the number of block nonzeros per row in the o matrix.
3159:   In general, for PDE problems in which most nonzeros are near the diagonal,
3160:   one expects `d_nz` >> `o_nz`.

3162: .seealso: `Mat`, `MatCreate()`, `MatCreateSeqBAIJ()`, `MatSetValues()`, `MatMPIBAIJSetPreallocation()`, `MatMPIBAIJSetPreallocationCSR()`,
3163:           `MatGetOwnershipRange()`,  `MatGetOwnershipRanges()`, `MatGetOwnershipRangeColumn()`, `MatGetOwnershipRangesColumn()`, `PetscLayout`
3164: @*/
3165: PetscErrorCode MatCreateBAIJ(MPI_Comm comm, PetscInt bs, PetscInt m, PetscInt n, PetscInt M, PetscInt N, PetscInt d_nz, const PetscInt d_nnz[], PetscInt o_nz, const PetscInt o_nnz[], Mat *A)
3166: {
3167:   PetscMPIInt size;

3169:   PetscFunctionBegin;
3170:   PetscCall(MatCreate(comm, A));
3171:   PetscCall(MatSetSizes(*A, m, n, M, N));
3172:   PetscCallMPI(MPI_Comm_size(comm, &size));
3173:   if (size > 1) {
3174:     PetscCall(MatSetType(*A, MATMPIBAIJ));
3175:     PetscCall(MatMPIBAIJSetPreallocation(*A, bs, d_nz, d_nnz, o_nz, o_nnz));
3176:   } else {
3177:     PetscCall(MatSetType(*A, MATSEQBAIJ));
3178:     PetscCall(MatSeqBAIJSetPreallocation(*A, bs, d_nz, d_nnz));
3179:   }
3180:   PetscFunctionReturn(PETSC_SUCCESS);
3181: }

3183: static PetscErrorCode MatDuplicate_MPIBAIJ(Mat matin, MatDuplicateOption cpvalues, Mat *newmat)
3184: {
3185:   Mat          mat;
3186:   Mat_MPIBAIJ *a, *oldmat = (Mat_MPIBAIJ *)matin->data;
3187:   PetscInt     len = 0;

3189:   PetscFunctionBegin;
3190:   *newmat = NULL;
3191:   PetscCall(MatCreate(PetscObjectComm((PetscObject)matin), &mat));
3192:   PetscCall(MatSetSizes(mat, matin->rmap->n, matin->cmap->n, matin->rmap->N, matin->cmap->N));
3193:   PetscCall(MatSetType(mat, ((PetscObject)matin)->type_name));

3195:   PetscCall(PetscLayoutReference(matin->rmap, &mat->rmap));
3196:   PetscCall(PetscLayoutReference(matin->cmap, &mat->cmap));
3197:   if (matin->hash_active) {
3198:     PetscCall(MatSetUp(mat));
3199:   } else {
3200:     mat->factortype   = matin->factortype;
3201:     mat->preallocated = PETSC_TRUE;
3202:     mat->assembled    = PETSC_TRUE;
3203:     mat->insertmode   = NOT_SET_VALUES;

3205:     a             = (Mat_MPIBAIJ *)mat->data;
3206:     mat->rmap->bs = matin->rmap->bs;
3207:     a->bs2        = oldmat->bs2;
3208:     a->mbs        = oldmat->mbs;
3209:     a->nbs        = oldmat->nbs;
3210:     a->Mbs        = oldmat->Mbs;
3211:     a->Nbs        = oldmat->Nbs;

3213:     a->size         = oldmat->size;
3214:     a->rank         = oldmat->rank;
3215:     a->donotstash   = oldmat->donotstash;
3216:     a->roworiented  = oldmat->roworiented;
3217:     a->rowindices   = NULL;
3218:     a->rowvalues    = NULL;
3219:     a->getrowactive = PETSC_FALSE;
3220:     a->barray       = NULL;
3221:     a->rstartbs     = oldmat->rstartbs;
3222:     a->rendbs       = oldmat->rendbs;
3223:     a->cstartbs     = oldmat->cstartbs;
3224:     a->cendbs       = oldmat->cendbs;

3226:     /* hash table stuff */
3227:     a->ht           = NULL;
3228:     a->hd           = NULL;
3229:     a->ht_size      = 0;
3230:     a->ht_flag      = oldmat->ht_flag;
3231:     a->ht_fact      = oldmat->ht_fact;
3232:     a->ht_total_ct  = 0;
3233:     a->ht_insert_ct = 0;

3235:     PetscCall(PetscArraycpy(a->rangebs, oldmat->rangebs, a->size + 1));
3236:     if (oldmat->colmap) {
3237: #if defined(PETSC_USE_CTABLE)
3238:       PetscCall(PetscHMapIDuplicate(oldmat->colmap, &a->colmap));
3239: #else
3240:       PetscCall(PetscMalloc1(a->Nbs, &a->colmap));
3241:       PetscCall(PetscArraycpy(a->colmap, oldmat->colmap, a->Nbs));
3242: #endif
3243:     } else a->colmap = NULL;

3245:     if (oldmat->garray && (len = ((Mat_SeqBAIJ *)oldmat->B->data)->nbs)) {
3246:       PetscCall(PetscMalloc1(len, &a->garray));
3247:       PetscCall(PetscArraycpy(a->garray, oldmat->garray, len));
3248:     } else a->garray = NULL;

3250:     PetscCall(MatStashCreate_Private(PetscObjectComm((PetscObject)matin), matin->rmap->bs, &mat->bstash));
3251:     PetscCall(VecDuplicate(oldmat->lvec, &a->lvec));
3252:     PetscCall(VecScatterCopy(oldmat->Mvctx, &a->Mvctx));

3254:     PetscCall(MatDuplicate(oldmat->A, cpvalues, &a->A));
3255:     PetscCall(MatDuplicate(oldmat->B, cpvalues, &a->B));
3256:   }
3257:   PetscCall(PetscFunctionListDuplicate(((PetscObject)matin)->qlist, &((PetscObject)mat)->qlist));
3258:   *newmat = mat;
3259:   PetscFunctionReturn(PETSC_SUCCESS);
3260: }

3262: /* Used for both MPIBAIJ and MPISBAIJ matrices */
3263: PetscErrorCode MatLoad_MPIBAIJ_Binary(Mat mat, PetscViewer viewer)
3264: {
3265:   PetscInt     header[4], M, N, nz, bs, m, n, mbs, nbs, rows, cols, sum, i, j, k;
3266:   PetscInt    *rowidxs, *colidxs, rs, cs, ce;
3267:   PetscScalar *matvals;

3269:   PetscFunctionBegin;
3270:   PetscCall(PetscViewerSetUp(viewer));

3272:   /* read in matrix header */
3273:   PetscCall(PetscViewerBinaryRead(viewer, header, 4, NULL, PETSC_INT));
3274:   PetscCheck(header[0] == MAT_FILE_CLASSID, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Not a matrix object in file");
3275:   M  = header[1];
3276:   N  = header[2];
3277:   nz = header[3];
3278:   PetscCheck(M >= 0, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Matrix row size (%" PetscInt_FMT ") in file is negative", M);
3279:   PetscCheck(N >= 0, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Matrix column size (%" PetscInt_FMT ") in file is negative", N);
3280:   PetscCheck(nz >= 0, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Matrix stored in special format on disk, cannot load as MPIBAIJ");

3282:   /* set block sizes from the viewer's .info file */
3283:   PetscCall(MatLoad_Binary_BlockSizes(mat, viewer));
3284:   /* set local sizes if not set already */
3285:   if (mat->rmap->n < 0 && M == N) mat->rmap->n = mat->cmap->n;
3286:   if (mat->cmap->n < 0 && M == N) mat->cmap->n = mat->rmap->n;
3287:   /* set global sizes if not set already */
3288:   if (mat->rmap->N < 0) mat->rmap->N = M;
3289:   if (mat->cmap->N < 0) mat->cmap->N = N;
3290:   PetscCall(PetscLayoutSetUp(mat->rmap));
3291:   PetscCall(PetscLayoutSetUp(mat->cmap));

3293:   /* check if the matrix sizes are correct */
3294:   PetscCall(MatGetSize(mat, &rows, &cols));
3295:   PetscCheck(M == rows && N == cols, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Matrix in file of different sizes (%" PetscInt_FMT ", %" PetscInt_FMT ") than the input matrix (%" PetscInt_FMT ", %" PetscInt_FMT ")", M, N, rows, cols);
3296:   PetscCall(MatGetBlockSize(mat, &bs));
3297:   PetscCall(MatGetLocalSize(mat, &m, &n));
3298:   PetscCall(PetscLayoutGetRange(mat->rmap, &rs, NULL));
3299:   PetscCall(PetscLayoutGetRange(mat->cmap, &cs, &ce));
3300:   mbs = m / bs;
3301:   nbs = n / bs;

3303:   /* read in row lengths and build row indices */
3304:   PetscCall(PetscMalloc1(m + 1, &rowidxs));
3305:   PetscCall(PetscViewerBinaryReadAll(viewer, rowidxs + 1, m, PETSC_DECIDE, M, PETSC_INT));
3306:   rowidxs[0] = 0;
3307:   for (i = 0; i < m; i++) rowidxs[i + 1] += rowidxs[i];
3308:   PetscCallMPI(MPIU_Allreduce(&rowidxs[m], &sum, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)viewer)));
3309:   PetscCheck(sum == nz, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Inconsistent matrix data in file: nonzeros = %" PetscInt_FMT ", sum-row-lengths = %" PetscInt_FMT, nz, sum);

3311:   /* read in column indices and matrix values */
3312:   PetscCall(PetscMalloc2(rowidxs[m], &colidxs, rowidxs[m], &matvals));
3313:   PetscCall(PetscViewerBinaryReadAll(viewer, colidxs, rowidxs[m], PETSC_DETERMINE, PETSC_DETERMINE, PETSC_INT));
3314:   PetscCall(PetscViewerBinaryReadAll(viewer, matvals, rowidxs[m], PETSC_DETERMINE, PETSC_DETERMINE, PETSC_SCALAR));

3316:   {                /* preallocate matrix storage */
3317:     PetscBT    bt; /* helper bit set to count diagonal nonzeros */
3318:     PetscHSetI ht; /* helper hash set to count off-diagonal nonzeros */
3319:     PetscBool  sbaij, done;
3320:     PetscInt  *d_nnz, *o_nnz;

3322:     PetscCall(PetscBTCreate(nbs, &bt));
3323:     PetscCall(PetscHSetICreate(&ht));
3324:     PetscCall(PetscCalloc2(mbs, &d_nnz, mbs, &o_nnz));
3325:     PetscCall(PetscObjectTypeCompare((PetscObject)mat, MATMPISBAIJ, &sbaij));
3326:     for (i = 0; i < mbs; i++) {
3327:       PetscCall(PetscBTMemzero(nbs, bt));
3328:       PetscCall(PetscHSetIClear(ht));
3329:       for (k = 0; k < bs; k++) {
3330:         PetscInt row = bs * i + k;
3331:         for (j = rowidxs[row]; j < rowidxs[row + 1]; j++) {
3332:           PetscInt col = colidxs[j];
3333:           if (!sbaij || col >= row) {
3334:             if (col >= cs && col < ce) {
3335:               if (!PetscBTLookupSet(bt, (col - cs) / bs)) d_nnz[i]++;
3336:             } else {
3337:               PetscCall(PetscHSetIQueryAdd(ht, col / bs, &done));
3338:               if (done) o_nnz[i]++;
3339:             }
3340:           }
3341:         }
3342:       }
3343:     }
3344:     PetscCall(PetscBTDestroy(&bt));
3345:     PetscCall(PetscHSetIDestroy(&ht));
3346:     PetscCall(MatMPIBAIJSetPreallocation(mat, bs, 0, d_nnz, 0, o_nnz));
3347:     PetscCall(MatMPISBAIJSetPreallocation(mat, bs, 0, d_nnz, 0, o_nnz));
3348:     PetscCall(PetscFree2(d_nnz, o_nnz));
3349:   }

3351:   /* store matrix values */
3352:   for (i = 0; i < m; i++) {
3353:     PetscInt row = rs + i, s = rowidxs[i], e = rowidxs[i + 1];
3354:     PetscUseTypeMethod(mat, setvalues, 1, &row, e - s, colidxs + s, matvals + s, INSERT_VALUES);
3355:   }

3357:   PetscCall(PetscFree(rowidxs));
3358:   PetscCall(PetscFree2(colidxs, matvals));
3359:   PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
3360:   PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
3361:   PetscFunctionReturn(PETSC_SUCCESS);
3362: }

3364: PetscErrorCode MatLoad_MPIBAIJ(Mat mat, PetscViewer viewer)
3365: {
3366:   PetscBool isbinary;

3368:   PetscFunctionBegin;
3369:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
3370:   PetscCheck(isbinary, PetscObjectComm((PetscObject)viewer), PETSC_ERR_SUP, "Viewer type %s not yet supported for reading %s matrices", ((PetscObject)viewer)->type_name, ((PetscObject)mat)->type_name);
3371:   PetscCall(MatLoad_MPIBAIJ_Binary(mat, viewer));
3372:   PetscFunctionReturn(PETSC_SUCCESS);
3373: }

3375: /*@
3376:   MatMPIBAIJSetHashTableFactor - Sets the factor required to compute the size of the matrices hash table

3378:   Input Parameters:
3379: + mat  - the matrix
3380: - fact - factor

3382:   Options Database Key:
3383: . -mat_use_hash_table <fact> - provide the factor

3385:   Level: advanced

3387: .seealso: `Mat`, `MATMPIBAIJ`, `MatSetOption()`
3388: @*/
3389: PetscErrorCode MatMPIBAIJSetHashTableFactor(Mat mat, PetscReal fact)
3390: {
3391:   PetscFunctionBegin;
3392:   PetscTryMethod(mat, "MatSetHashTableFactor_C", (Mat, PetscReal), (mat, fact));
3393:   PetscFunctionReturn(PETSC_SUCCESS);
3394: }

3396: PetscErrorCode MatSetHashTableFactor_MPIBAIJ(Mat mat, PetscReal fact)
3397: {
3398:   Mat_MPIBAIJ *baij;

3400:   PetscFunctionBegin;
3401:   baij          = (Mat_MPIBAIJ *)mat->data;
3402:   baij->ht_fact = fact;
3403:   PetscFunctionReturn(PETSC_SUCCESS);
3404: }

3406: PetscErrorCode MatMPIBAIJGetSeqBAIJ(Mat A, Mat *Ad, Mat *Ao, const PetscInt *colmap[])
3407: {
3408:   Mat_MPIBAIJ *a = (Mat_MPIBAIJ *)A->data;
3409:   PetscBool    flg;

3411:   PetscFunctionBegin;
3412:   PetscCall(PetscObjectTypeCompare((PetscObject)A, MATMPIBAIJ, &flg));
3413:   PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "This function requires a MATMPIBAIJ matrix as input");
3414:   if (Ad) *Ad = a->A;
3415:   if (Ao) *Ao = a->B;
3416:   if (colmap) *colmap = a->garray;
3417:   PetscFunctionReturn(PETSC_SUCCESS);
3418: }

3420: /*
3421:     Special version for direct calls from Fortran (to eliminate two function call overheads
3422: */
3423: #if defined(PETSC_HAVE_FORTRAN_CAPS)
3424:   #define matmpibaijsetvaluesblocked_ MATMPIBAIJSETVALUESBLOCKED
3425: #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE)
3426:   #define matmpibaijsetvaluesblocked_ matmpibaijsetvaluesblocked
3427: #endif

3429: // PetscClangLinter pragma disable: -fdoc-synopsis-matching-symbol-name
3430: /*@C
3431:   MatMPIBAIJSetValuesBlocked - Direct Fortran call to replace call to `MatSetValuesBlocked()`

3433:   Collective

3435:   Input Parameters:
3436: + matin  - the matrix
3437: . min    - number of input rows
3438: . im     - input rows
3439: . nin    - number of input columns
3440: . in     - input columns
3441: . v      - numerical values input
3442: - addvin - `INSERT_VALUES` or `ADD_VALUES`

3444:   Level: advanced

3446:   Developer Notes:
3447:   This has a complete copy of `MatSetValuesBlocked_MPIBAIJ()` which is terrible code un-reuse.

3449: .seealso: `Mat`, `MatSetValuesBlocked()`
3450: @*/
3451: PETSC_EXTERN PetscErrorCode matmpibaijsetvaluesblocked_(Mat *matin, PetscInt *min, const PetscInt im[], PetscInt *nin, const PetscInt in[], const MatScalar v[], InsertMode *addvin)
3452: {
3453:   /* convert input arguments to C version */
3454:   Mat        mat = *matin;
3455:   PetscInt   m = *min, n = *nin;
3456:   InsertMode addv = *addvin;

3458:   Mat_MPIBAIJ     *baij = (Mat_MPIBAIJ *)mat->data;
3459:   const MatScalar *value;
3460:   MatScalar       *barray      = baij->barray;
3461:   PetscBool        roworiented = baij->roworiented;
3462:   PetscInt         i, j, ii, jj, row, col, rstart = baij->rstartbs;
3463:   PetscInt         rend = baij->rendbs, cstart = baij->cstartbs, stepval;
3464:   PetscInt         cend = baij->cendbs, bs = mat->rmap->bs, bs2 = baij->bs2;

3466:   PetscFunctionBegin;
3467:   /* tasks normally handled by MatSetValuesBlocked() */
3468:   if (mat->insertmode == NOT_SET_VALUES) mat->insertmode = addv;
3469:   else PetscCheck(mat->insertmode == addv, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Cannot mix add values and insert values");
3470:   PetscCheck(!mat->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
3471:   if (mat->assembled) {
3472:     mat->was_assembled = PETSC_TRUE;
3473:     mat->assembled     = PETSC_FALSE;
3474:   }
3475:   PetscCall(PetscLogEventBegin(MAT_SetValues, mat, 0, 0, 0));

3477:   if (!barray) {
3478:     PetscCall(PetscMalloc1(bs2, &barray));
3479:     baij->barray = barray;
3480:   }

3482:   if (roworiented) stepval = (n - 1) * bs;
3483:   else stepval = (m - 1) * bs;

3485:   for (i = 0; i < m; i++) {
3486:     if (im[i] < 0) continue;
3487:     PetscCheck(im[i] < baij->Mbs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row too large, row %" PetscInt_FMT " max %" PetscInt_FMT, im[i], baij->Mbs - 1);
3488:     if (im[i] >= rstart && im[i] < rend) {
3489:       row = im[i] - rstart;
3490:       for (j = 0; j < n; j++) {
3491:         /* If NumCol = 1 then a copy is not required */
3492:         if ((roworiented) && (n == 1)) {
3493:           barray = (MatScalar *)v + i * bs2;
3494:         } else if ((!roworiented) && (m == 1)) {
3495:           barray = (MatScalar *)v + j * bs2;
3496:         } else { /* Here a copy is required */
3497:           if (roworiented) {
3498:             value = v + i * (stepval + bs) * bs + j * bs;
3499:           } else {
3500:             value = v + j * (stepval + bs) * bs + i * bs;
3501:           }
3502:           for (ii = 0; ii < bs; ii++, value += stepval) {
3503:             for (jj = 0; jj < bs; jj++) *barray++ = *value++;
3504:           }
3505:           barray -= bs2;
3506:         }

3508:         if (in[j] >= cstart && in[j] < cend) {
3509:           col = in[j] - cstart;
3510:           PetscCall(MatSetValuesBlocked_SeqBAIJ_Inlined(baij->A, row, col, barray, addv, im[i], in[j]));
3511:         } else if (in[j] < 0) {
3512:           continue;
3513:         } else {
3514:           PetscCheck(in[j] < baij->Nbs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Column too large, col %" PetscInt_FMT " max %" PetscInt_FMT, in[j], baij->Nbs - 1);
3515:           if (mat->was_assembled) {
3516:             if (!baij->colmap) PetscCall(MatCreateColmap_MPIBAIJ_Private(mat));

3518: #if defined(PETSC_USE_DEBUG)
3519:   #if defined(PETSC_USE_CTABLE)
3520:             {
3521:               PetscInt data;
3522:               PetscCall(PetscHMapIGetWithDefault(baij->colmap, in[j] + 1, 0, &data));
3523:               PetscCheck((data - 1) % bs == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Incorrect colmap");
3524:             }
3525:   #else
3526:             PetscCheck((baij->colmap[in[j]] - 1) % bs == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Incorrect colmap");
3527:   #endif
3528: #endif
3529: #if defined(PETSC_USE_CTABLE)
3530:             PetscCall(PetscHMapIGetWithDefault(baij->colmap, in[j] + 1, 0, &col));
3531:             col = (col - 1) / bs;
3532: #else
3533:             col = (baij->colmap[in[j]] - 1) / bs;
3534: #endif
3535:             if (col < 0 && !((Mat_SeqBAIJ *)baij->A->data)->nonew) {
3536:               PetscCall(MatDisAssemble_MPIBAIJ(mat));
3537:               col = in[j];
3538:             }
3539:           } else col = in[j];
3540:           PetscCall(MatSetValuesBlocked_SeqBAIJ_Inlined(baij->B, row, col, barray, addv, im[i], in[j]));
3541:         }
3542:       }
3543:     } else {
3544:       if (!baij->donotstash) {
3545:         if (roworiented) {
3546:           PetscCall(MatStashValuesRowBlocked_Private(&mat->bstash, im[i], n, in, v, m, n, i));
3547:         } else {
3548:           PetscCall(MatStashValuesColBlocked_Private(&mat->bstash, im[i], n, in, v, m, n, i));
3549:         }
3550:       }
3551:     }
3552:   }

3554:   /* task normally handled by MatSetValuesBlocked() */
3555:   PetscCall(PetscLogEventEnd(MAT_SetValues, mat, 0, 0, 0));
3556:   PetscFunctionReturn(PETSC_SUCCESS);
3557: }

3559: /*@
3560:   MatCreateMPIBAIJWithArrays - creates a `MATMPIBAIJ` matrix using arrays that contain in standard block CSR format for the local rows.

3562:   Collective

3564:   Input Parameters:
3565: + comm - MPI communicator
3566: . bs   - the block size, only a block size of 1 is supported
3567: . m    - number of local rows (Cannot be `PETSC_DECIDE`)
3568: . n    - This value should be the same as the local size used in creating the
3569:          x vector for the matrix-vector product $ y = Ax $. (or `PETSC_DECIDE` to have
3570:          calculated if `N` is given) For square matrices `n` is almost always `m`.
3571: . M    - number of global rows (or `PETSC_DETERMINE` to have calculated if `m` is given)
3572: . N    - number of global columns (or `PETSC_DETERMINE` to have calculated if `n` is given)
3573: . i    - row indices; that is i[0] = 0, i[row] = i[row-1] + number of block elements in that rowth block row of the matrix
3574: . j    - column indices
3575: - a    - matrix values

3577:   Output Parameter:
3578: . mat - the matrix

3580:   Level: intermediate

3582:   Notes:
3583:   The `i`, `j`, and `a` arrays ARE copied by this routine into the internal format used by PETSc;
3584:   thus you CANNOT change the matrix entries by changing the values of a[] after you have
3585:   called this routine. Use `MatCreateMPIAIJWithSplitArrays()` to avoid needing to copy the arrays.

3587:   The order of the entries in values is the same as the block compressed sparse row storage format; that is, it is
3588:   the same as a three dimensional array in Fortran values(bs,bs,nnz) that contains the first column of the first
3589:   block, followed by the second column of the first block etc etc.  That is, the blocks are contiguous in memory
3590:   with column-major ordering within blocks.

3592:   The `i` and `j` indices are 0 based, and `i` indices are indices corresponding to the local `j` array.

3594: .seealso: `Mat`, `MatCreate()`, `MatCreateSeqAIJ()`, `MatSetValues()`, `MatMPIAIJSetPreallocation()`, `MatMPIAIJSetPreallocationCSR()`,
3595:           `MATMPIAIJ`, `MatCreateAIJ()`, `MatCreateMPIAIJWithSplitArrays()`
3596: @*/
3597: PetscErrorCode MatCreateMPIBAIJWithArrays(MPI_Comm comm, PetscInt bs, PetscInt m, PetscInt n, PetscInt M, PetscInt N, const PetscInt i[], const PetscInt j[], const PetscScalar a[], Mat *mat)
3598: {
3599:   PetscFunctionBegin;
3600:   PetscCheck(!i[0], PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "i (row indices) must start with 0");
3601:   PetscCheck(m >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "local number of rows (m) cannot be PETSC_DECIDE, or negative");
3602:   PetscCall(MatCreate(comm, mat));
3603:   PetscCall(MatSetSizes(*mat, m, n, M, N));
3604:   PetscCall(MatSetType(*mat, MATMPIBAIJ));
3605:   PetscCall(MatSetBlockSize(*mat, bs));
3606:   PetscCall(MatSetUp(*mat));
3607:   PetscCall(MatSetOption(*mat, MAT_ROW_ORIENTED, PETSC_FALSE));
3608:   PetscCall(MatMPIBAIJSetPreallocationCSR(*mat, bs, i, j, a));
3609:   PetscCall(MatSetOption(*mat, MAT_ROW_ORIENTED, PETSC_TRUE));
3610:   PetscFunctionReturn(PETSC_SUCCESS);
3611: }

3613: PetscErrorCode MatCreateMPIMatConcatenateSeqMat_MPIBAIJ(MPI_Comm comm, Mat inmat, PetscInt n, MatReuse scall, Mat *outmat)
3614: {
3615:   PetscInt     m, N, i, rstart, nnz, Ii, bs, cbs;
3616:   PetscInt    *indx;
3617:   PetscScalar *values;

3619:   PetscFunctionBegin;
3620:   PetscCall(MatGetSize(inmat, &m, &N));
3621:   if (scall == MAT_INITIAL_MATRIX) { /* symbolic phase */
3622:     Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)inmat->data;
3623:     PetscInt    *dnz, *onz, mbs, Nbs, nbs;
3624:     PetscInt    *bindx, rmax = a->rmax, j;
3625:     PetscMPIInt  rank, size;

3627:     PetscCall(MatGetBlockSizes(inmat, &bs, &cbs));
3628:     mbs = m / bs;
3629:     Nbs = N / cbs;
3630:     if (n == PETSC_DECIDE) PetscCall(PetscSplitOwnershipBlock(comm, cbs, &n, &N));
3631:     nbs = n / cbs;

3633:     PetscCall(PetscMalloc1(rmax, &bindx));
3634:     MatPreallocateBegin(comm, mbs, nbs, dnz, onz); /* inline function, output __end and __rstart are used below */

3636:     PetscCallMPI(MPI_Comm_rank(comm, &rank));
3637:     PetscCallMPI(MPI_Comm_rank(comm, &size));
3638:     if (rank == size - 1) {
3639:       /* Check sum(nbs) = Nbs */
3640:       PetscCheck(__end == Nbs, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Sum of local block columns %" PetscInt_FMT " != global block columns %" PetscInt_FMT, __end, Nbs);
3641:     }

3643:     rstart = __rstart; /* block rstart of *outmat; see inline function MatPreallocateBegin */
3644:     for (i = 0; i < mbs; i++) {
3645:       PetscCall(MatGetRow_SeqBAIJ(inmat, i * bs, &nnz, &indx, NULL)); /* non-blocked nnz and indx */
3646:       nnz = nnz / bs;
3647:       for (j = 0; j < nnz; j++) bindx[j] = indx[j * bs] / bs;
3648:       PetscCall(MatPreallocateSet(i + rstart, nnz, bindx, dnz, onz));
3649:       PetscCall(MatRestoreRow_SeqBAIJ(inmat, i * bs, &nnz, &indx, NULL));
3650:     }
3651:     PetscCall(PetscFree(bindx));

3653:     PetscCall(MatCreate(comm, outmat));
3654:     PetscCall(MatSetSizes(*outmat, m, n, PETSC_DETERMINE, PETSC_DETERMINE));
3655:     PetscCall(MatSetBlockSizes(*outmat, bs, cbs));
3656:     PetscCall(MatSetType(*outmat, MATBAIJ));
3657:     PetscCall(MatSeqBAIJSetPreallocation(*outmat, bs, 0, dnz));
3658:     PetscCall(MatMPIBAIJSetPreallocation(*outmat, bs, 0, dnz, 0, onz));
3659:     MatPreallocateEnd(dnz, onz);
3660:     PetscCall(MatSetOption(*outmat, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
3661:   }

3663:   /* numeric phase */
3664:   PetscCall(MatGetBlockSizes(inmat, &bs, &cbs));
3665:   PetscCall(MatGetOwnershipRange(*outmat, &rstart, NULL));

3667:   for (i = 0; i < m; i++) {
3668:     PetscCall(MatGetRow_SeqBAIJ(inmat, i, &nnz, &indx, &values));
3669:     Ii = i + rstart;
3670:     PetscCall(MatSetValues(*outmat, 1, &Ii, nnz, indx, values, INSERT_VALUES));
3671:     PetscCall(MatRestoreRow_SeqBAIJ(inmat, i, &nnz, &indx, &values));
3672:   }
3673:   PetscCall(MatAssemblyBegin(*outmat, MAT_FINAL_ASSEMBLY));
3674:   PetscCall(MatAssemblyEnd(*outmat, MAT_FINAL_ASSEMBLY));
3675:   PetscFunctionReturn(PETSC_SUCCESS);
3676: }