Actual source code: mpisbaij.c

      1: #include <../src/mat/impls/baij/mpi/mpibaij.h>
      2: #include <../src/mat/impls/sbaij/mpi/mpisbaij.h>
      3: #include <../src/mat/impls/sbaij/seq/sbaij.h>
      4: #include <petscblaslapack.h>
      5: #include <petscsf.h>

      7: static PetscErrorCode MatDestroy_MPISBAIJ(Mat mat)
      8: {
      9:   Mat_MPISBAIJ *baij = (Mat_MPISBAIJ *)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(VecDestroy(&baij->slvec0));
     26:   PetscCall(VecDestroy(&baij->slvec0b));
     27:   PetscCall(VecDestroy(&baij->slvec1));
     28:   PetscCall(VecDestroy(&baij->slvec1a));
     29:   PetscCall(VecDestroy(&baij->slvec1b));
     30:   PetscCall(VecScatterDestroy(&baij->sMvctx));
     31:   PetscCall(PetscFree2(baij->rowvalues, baij->rowindices));
     32:   PetscCall(PetscFree(baij->barray));
     33:   PetscCall(PetscFree(baij->hd));
     34:   PetscCall(VecDestroy(&baij->diag));
     35:   PetscCall(VecDestroy(&baij->bb1));
     36:   PetscCall(VecDestroy(&baij->xx1));
     37: #if defined(PETSC_USE_REAL_MAT_SINGLE)
     38:   PetscCall(PetscFree(baij->setvaluescopy));
     39: #endif
     40:   PetscCall(PetscFree(baij->in_loc));
     41:   PetscCall(PetscFree(baij->v_loc));
     42:   PetscCall(PetscFree(baij->rangebs));
     43:   PetscCall(PetscFree(mat->data));

     45:   PetscCall(PetscObjectChangeTypeName((PetscObject)mat, NULL));
     46:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatStoreValues_C", NULL));
     47:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatRetrieveValues_C", NULL));
     48:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMPISBAIJSetPreallocation_C", NULL));
     49:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMPISBAIJSetPreallocationCSR_C", NULL));
     50: #if defined(PETSC_HAVE_ELEMENTAL)
     51:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_mpisbaij_elemental_C", NULL));
     52: #endif
     53: #if defined(PETSC_HAVE_SCALAPACK) && (defined(PETSC_USE_REAL_SINGLE) || defined(PETSC_USE_REAL_DOUBLE))
     54:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_mpisbaij_scalapack_C", NULL));
     55: #endif
     56:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_mpisbaij_mpiaij_C", NULL));
     57:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_mpisbaij_mpibaij_C", NULL));
     58:   PetscFunctionReturn(PETSC_SUCCESS);
     59: }

     61: /* defines MatSetValues_MPI_Hash(), MatAssemblyBegin_MPI_Hash(), MatAssemblyEnd_MPI_Hash(), MatSetUp_MPI_Hash() */
     62: #define TYPE SBAIJ
     63: #define TYPE_SBAIJ
     64: #include "../src/mat/impls/aij/mpi/mpihashmat.h"
     65: #undef TYPE
     66: #undef TYPE_SBAIJ

     68: #if defined(PETSC_HAVE_ELEMENTAL)
     69: PETSC_INTERN PetscErrorCode MatConvert_MPISBAIJ_Elemental(Mat, MatType, MatReuse, Mat *);
     70: #endif
     71: #if defined(PETSC_HAVE_SCALAPACK) && (defined(PETSC_USE_REAL_SINGLE) || defined(PETSC_USE_REAL_DOUBLE))
     72: PETSC_INTERN PetscErrorCode MatConvert_SBAIJ_ScaLAPACK(Mat, MatType, MatReuse, Mat *);
     73: #endif

     75: /* This could be moved to matimpl.h */
     76: static PetscErrorCode MatPreallocateWithMats_Private(Mat B, PetscInt nm, Mat X[], PetscBool symm[], PetscBool fill)
     77: {
     78:   Mat       preallocator;
     79:   PetscInt  r, rstart, rend;
     80:   PetscInt  bs, i, m, n, M, N;
     81:   PetscBool cong = PETSC_TRUE;

     83:   PetscFunctionBegin;
     86:   for (i = 0; i < nm; i++) {
     88:     PetscCall(PetscLayoutCompare(B->rmap, X[i]->rmap, &cong));
     89:     PetscCheck(cong, PetscObjectComm((PetscObject)B), PETSC_ERR_SUP, "Not for different layouts");
     90:   }
     92:   PetscCall(MatGetBlockSize(B, &bs));
     93:   PetscCall(MatGetSize(B, &M, &N));
     94:   PetscCall(MatGetLocalSize(B, &m, &n));
     95:   PetscCall(MatCreate(PetscObjectComm((PetscObject)B), &preallocator));
     96:   PetscCall(MatSetType(preallocator, MATPREALLOCATOR));
     97:   PetscCall(MatSetBlockSize(preallocator, bs));
     98:   PetscCall(MatSetSizes(preallocator, m, n, M, N));
     99:   PetscCall(MatSetUp(preallocator));
    100:   PetscCall(MatGetOwnershipRange(preallocator, &rstart, &rend));
    101:   for (r = rstart; r < rend; ++r) {
    102:     PetscInt           ncols;
    103:     const PetscInt    *row;
    104:     const PetscScalar *vals;

    106:     for (i = 0; i < nm; i++) {
    107:       PetscCall(MatGetRow(X[i], r, &ncols, &row, &vals));
    108:       PetscCall(MatSetValues(preallocator, 1, &r, ncols, row, vals, INSERT_VALUES));
    109:       if (symm && symm[i]) PetscCall(MatSetValues(preallocator, ncols, row, 1, &r, vals, INSERT_VALUES));
    110:       PetscCall(MatRestoreRow(X[i], r, &ncols, &row, &vals));
    111:     }
    112:   }
    113:   PetscCall(MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY));
    114:   PetscCall(MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY));
    115:   PetscCall(MatPreallocatorPreallocate(preallocator, fill, B));
    116:   PetscCall(MatDestroy(&preallocator));
    117:   PetscFunctionReturn(PETSC_SUCCESS);
    118: }

    120: PETSC_INTERN PetscErrorCode MatConvert_MPISBAIJ_Basic(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
    121: {
    122:   Mat      B;
    123:   PetscInt r;

    125:   PetscFunctionBegin;
    126:   if (reuse != MAT_REUSE_MATRIX) {
    127:     PetscBool symm = PETSC_TRUE, isdense;
    128:     PetscInt  bs;

    130:     PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
    131:     PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
    132:     PetscCall(MatSetType(B, newtype));
    133:     PetscCall(MatGetBlockSize(A, &bs));
    134:     PetscCall(MatSetBlockSize(B, bs));
    135:     PetscCall(PetscLayoutSetUp(B->rmap));
    136:     PetscCall(PetscLayoutSetUp(B->cmap));
    137:     PetscCall(PetscObjectTypeCompareAny((PetscObject)B, &isdense, MATSEQDENSE, MATMPIDENSE, MATSEQDENSECUDA, ""));
    138:     if (!isdense) {
    139:       PetscCall(MatGetRowUpperTriangular(A));
    140:       PetscCall(MatPreallocateWithMats_Private(B, 1, &A, &symm, PETSC_TRUE));
    141:       PetscCall(MatRestoreRowUpperTriangular(A));
    142:     } else {
    143:       PetscCall(MatSetUp(B));
    144:     }
    145:   } else {
    146:     B = *newmat;
    147:     PetscCall(MatZeroEntries(B));
    148:   }

    150:   PetscCall(MatGetRowUpperTriangular(A));
    151:   for (r = A->rmap->rstart; r < A->rmap->rend; r++) {
    152:     PetscInt           ncols;
    153:     const PetscInt    *row;
    154:     const PetscScalar *vals;

    156:     PetscCall(MatGetRow(A, r, &ncols, &row, &vals));
    157:     PetscCall(MatSetValues(B, 1, &r, ncols, row, vals, INSERT_VALUES));
    158: #if defined(PETSC_USE_COMPLEX)
    159:     if (A->hermitian == PETSC_BOOL3_TRUE) {
    160:       PetscInt i;
    161:       for (i = 0; i < ncols; i++) PetscCall(MatSetValue(B, row[i], r, PetscConj(vals[i]), INSERT_VALUES));
    162:     } else {
    163:       PetscCall(MatSetValues(B, ncols, row, 1, &r, vals, INSERT_VALUES));
    164:     }
    165: #else
    166:     PetscCall(MatSetValues(B, ncols, row, 1, &r, vals, INSERT_VALUES));
    167: #endif
    168:     PetscCall(MatRestoreRow(A, r, &ncols, &row, &vals));
    169:   }
    170:   PetscCall(MatRestoreRowUpperTriangular(A));
    171:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
    172:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));

    174:   if (reuse == MAT_INPLACE_MATRIX) {
    175:     PetscCall(MatHeaderReplace(A, &B));
    176:   } else {
    177:     *newmat = B;
    178:   }
    179:   PetscFunctionReturn(PETSC_SUCCESS);
    180: }

    182: static PetscErrorCode MatStoreValues_MPISBAIJ(Mat mat)
    183: {
    184:   Mat_MPISBAIJ *aij = (Mat_MPISBAIJ *)mat->data;

    186:   PetscFunctionBegin;
    187:   PetscCall(MatStoreValues(aij->A));
    188:   PetscCall(MatStoreValues(aij->B));
    189:   PetscFunctionReturn(PETSC_SUCCESS);
    190: }

    192: static PetscErrorCode MatRetrieveValues_MPISBAIJ(Mat mat)
    193: {
    194:   Mat_MPISBAIJ *aij = (Mat_MPISBAIJ *)mat->data;

    196:   PetscFunctionBegin;
    197:   PetscCall(MatRetrieveValues(aij->A));
    198:   PetscCall(MatRetrieveValues(aij->B));
    199:   PetscFunctionReturn(PETSC_SUCCESS);
    200: }

    202: #define MatSetValues_SeqSBAIJ_A_Private(row, col, value, addv, orow, ocol) \
    203:   do { \
    204:     brow = row / bs; \
    205:     rp   = aj + ai[brow]; \
    206:     ap   = aa + bs2 * ai[brow]; \
    207:     rmax = aimax[brow]; \
    208:     nrow = ailen[brow]; \
    209:     bcol = col / bs; \
    210:     ridx = row % bs; \
    211:     cidx = col % bs; \
    212:     low  = 0; \
    213:     high = nrow; \
    214:     while (high - low > 3) { \
    215:       t = (low + high) / 2; \
    216:       if (rp[t] > bcol) high = t; \
    217:       else low = t; \
    218:     } \
    219:     for (_i = low; _i < high; _i++) { \
    220:       if (rp[_i] > bcol) break; \
    221:       if (rp[_i] == bcol) { \
    222:         bap = ap + bs2 * _i + bs * cidx + ridx; \
    223:         if (addv == ADD_VALUES) *bap += value; \
    224:         else *bap = value; \
    225:         goto a_noinsert; \
    226:       } \
    227:     } \
    228:     if (a->nonew == 1) goto a_noinsert; \
    229:     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); \
    230:     MatSeqXAIJReallocateAIJ(A, a->mbs, bs2, nrow, brow, bcol, rmax, aa, ai, aj, rp, ap, aimax, a->nonew, MatScalar); \
    231:     N = nrow++ - 1; \
    232:     /* shift up all the later entries in this row */ \
    233:     PetscCall(PetscArraymove(rp + _i + 1, rp + _i, N - _i + 1)); \
    234:     PetscCall(PetscArraymove(ap + bs2 * (_i + 1), ap + bs2 * _i, bs2 * (N - _i + 1))); \
    235:     PetscCall(PetscArrayzero(ap + bs2 * _i, bs2)); \
    236:     rp[_i]                          = bcol; \
    237:     ap[bs2 * _i + bs * cidx + ridx] = value; \
    238:   a_noinsert:; \
    239:     ailen[brow] = nrow; \
    240:   } while (0)

    242: #define MatSetValues_SeqSBAIJ_B_Private(row, col, value, addv, orow, ocol) \
    243:   do { \
    244:     brow = row / bs; \
    245:     rp   = bj + bi[brow]; \
    246:     ap   = ba + bs2 * bi[brow]; \
    247:     rmax = bimax[brow]; \
    248:     nrow = bilen[brow]; \
    249:     bcol = col / bs; \
    250:     ridx = row % bs; \
    251:     cidx = col % bs; \
    252:     low  = 0; \
    253:     high = nrow; \
    254:     while (high - low > 3) { \
    255:       t = (low + high) / 2; \
    256:       if (rp[t] > bcol) high = t; \
    257:       else low = t; \
    258:     } \
    259:     for (_i = low; _i < high; _i++) { \
    260:       if (rp[_i] > bcol) break; \
    261:       if (rp[_i] == bcol) { \
    262:         bap = ap + bs2 * _i + bs * cidx + ridx; \
    263:         if (addv == ADD_VALUES) *bap += value; \
    264:         else *bap = value; \
    265:         goto b_noinsert; \
    266:       } \
    267:     } \
    268:     if (b->nonew == 1) goto b_noinsert; \
    269:     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); \
    270:     MatSeqXAIJReallocateAIJ(B, b->mbs, bs2, nrow, brow, bcol, rmax, ba, bi, bj, rp, ap, bimax, b->nonew, MatScalar); \
    271:     N = nrow++ - 1; \
    272:     /* shift up all the later entries in this row */ \
    273:     PetscCall(PetscArraymove(rp + _i + 1, rp + _i, N - _i + 1)); \
    274:     PetscCall(PetscArraymove(ap + bs2 * (_i + 1), ap + bs2 * _i, bs2 * (N - _i + 1))); \
    275:     PetscCall(PetscArrayzero(ap + bs2 * _i, bs2)); \
    276:     rp[_i]                          = bcol; \
    277:     ap[bs2 * _i + bs * cidx + ridx] = value; \
    278:   b_noinsert:; \
    279:     bilen[brow] = nrow; \
    280:   } while (0)

    282: /* Only add/insert a(i,j) with i<=j (blocks).
    283:    Any a(i,j) with i>j input by user is ignored or generates an error
    284: */
    285: static PetscErrorCode MatSetValues_MPISBAIJ(Mat mat, PetscInt m, const PetscInt im[], PetscInt n, const PetscInt in[], const PetscScalar v[], InsertMode addv)
    286: {
    287:   Mat_MPISBAIJ *baij = (Mat_MPISBAIJ *)mat->data;
    288:   MatScalar     value;
    289:   PetscBool     roworiented = baij->roworiented;
    290:   PetscInt      i, j, row, col;
    291:   PetscInt      rstart_orig = mat->rmap->rstart;
    292:   PetscInt      rend_orig = mat->rmap->rend, cstart_orig = mat->cmap->rstart;
    293:   PetscInt      cend_orig = mat->cmap->rend, bs = mat->rmap->bs;

    295:   /* Some Variables required in the macro */
    296:   Mat           A     = baij->A;
    297:   Mat_SeqSBAIJ *a     = (Mat_SeqSBAIJ *)A->data;
    298:   PetscInt     *aimax = a->imax, *ai = a->i, *ailen = a->ilen, *aj = a->j;
    299:   MatScalar    *aa = a->a;

    301:   Mat          B     = baij->B;
    302:   Mat_SeqBAIJ *b     = (Mat_SeqBAIJ *)B->data;
    303:   PetscInt    *bimax = b->imax, *bi = b->i, *bilen = b->ilen, *bj = b->j;
    304:   MatScalar   *ba = b->a;

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

    310:   /* for stash */
    311:   PetscInt   n_loc, *in_loc = NULL;
    312:   MatScalar *v_loc = NULL;

    314:   PetscFunctionBegin;
    315:   if (!baij->donotstash) {
    316:     if (n > baij->n_loc) {
    317:       PetscCall(PetscFree(baij->in_loc));
    318:       PetscCall(PetscFree(baij->v_loc));
    319:       PetscCall(PetscMalloc1(n, &baij->in_loc));
    320:       PetscCall(PetscMalloc1(n, &baij->v_loc));

    322:       baij->n_loc = n;
    323:     }
    324:     in_loc = baij->in_loc;
    325:     v_loc  = baij->v_loc;
    326:   }

    328:   for (i = 0; i < m; i++) {
    329:     if (im[i] < 0) continue;
    330:     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);
    331:     if (im[i] >= rstart_orig && im[i] < rend_orig) { /* this processor entry */
    332:       row = im[i] - rstart_orig;                     /* local row index */
    333:       for (j = 0; j < n; j++) {
    334:         if (im[i] / bs > in[j] / bs) {
    335:           PetscCheck(a->ignore_ltriangular, PETSC_COMM_SELF, PETSC_ERR_USER, "Lower triangular value cannot be set for sbaij format. Ignoring these values, run with -mat_ignore_lower_triangular or call MatSetOption(mat,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE)");
    336:           continue; /* ignore lower triangular blocks */
    337:         }
    338:         if (in[j] >= cstart_orig && in[j] < cend_orig) { /* diag entry (A) */
    339:           col  = in[j] - cstart_orig;                    /* local col index */
    340:           brow = row / bs;
    341:           bcol = col / bs;
    342:           if (brow > bcol) continue; /* ignore lower triangular blocks of A */
    343:           if (roworiented) value = v[i * n + j];
    344:           else value = v[i + j * m];
    345:           MatSetValues_SeqSBAIJ_A_Private(row, col, value, addv, im[i], in[j]);
    346:           /* PetscCall(MatSetValues_SeqBAIJ(baij->A,1,&row,1,&col,&value,addv)); */
    347:         } else if (in[j] < 0) {
    348:           continue;
    349:         } else {
    350:           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);
    351:           /* off-diag entry (B) */
    352:           if (mat->was_assembled) {
    353:             if (!baij->colmap) PetscCall(MatCreateColmap_MPIBAIJ_Private(mat));
    354: #if defined(PETSC_USE_CTABLE)
    355:             PetscCall(PetscHMapIGetWithDefault(baij->colmap, in[j] / bs + 1, 0, &col));
    356:             col = col - 1;
    357: #else
    358:             col = baij->colmap[in[j] / bs] - 1;
    359: #endif
    360:             if (col < 0 && !((Mat_SeqSBAIJ *)baij->A->data)->nonew) {
    361:               PetscCall(MatDisAssemble_MPISBAIJ(mat));
    362:               col = in[j];
    363:               /* Reinitialize the variables required by MatSetValues_SeqBAIJ_B_Private() */
    364:               B     = baij->B;
    365:               b     = (Mat_SeqBAIJ *)B->data;
    366:               bimax = b->imax;
    367:               bi    = b->i;
    368:               bilen = b->ilen;
    369:               bj    = b->j;
    370:               ba    = b->a;
    371:             } else col += in[j] % bs;
    372:           } else col = in[j];
    373:           if (roworiented) value = v[i * n + j];
    374:           else value = v[i + j * m];
    375:           MatSetValues_SeqSBAIJ_B_Private(row, col, value, addv, im[i], in[j]);
    376:           /* PetscCall(MatSetValues_SeqBAIJ(baij->B,1,&row,1,&col,&value,addv)); */
    377:         }
    378:       }
    379:     } else { /* off processor entry */
    380:       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]);
    381:       if (!baij->donotstash) {
    382:         mat->assembled = PETSC_FALSE;
    383:         n_loc          = 0;
    384:         for (j = 0; j < n; j++) {
    385:           if (im[i] / bs > in[j] / bs) continue; /* ignore lower triangular blocks */
    386:           in_loc[n_loc] = in[j];
    387:           if (roworiented) {
    388:             v_loc[n_loc] = v[i * n + j];
    389:           } else {
    390:             v_loc[n_loc] = v[j * m + i];
    391:           }
    392:           n_loc++;
    393:         }
    394:         PetscCall(MatStashValuesRow_Private(&mat->stash, im[i], n_loc, in_loc, v_loc, PETSC_FALSE));
    395:       }
    396:     }
    397:   }
    398:   PetscFunctionReturn(PETSC_SUCCESS);
    399: }

    401: static inline PetscErrorCode MatSetValuesBlocked_SeqSBAIJ_Inlined(Mat A, PetscInt row, PetscInt col, const PetscScalar v[], InsertMode is, PetscInt orow, PetscInt ocol)
    402: {
    403:   Mat_SeqSBAIJ      *a = (Mat_SeqSBAIJ *)A->data;
    404:   PetscInt          *rp, low, high, t, ii, jj, nrow, i, rmax, N;
    405:   PetscInt          *imax = a->imax, *ai = a->i, *ailen = a->ilen;
    406:   PetscInt          *aj = a->j, nonew = a->nonew, bs2 = a->bs2, bs = A->rmap->bs;
    407:   PetscBool          roworiented = a->roworiented;
    408:   const PetscScalar *value       = v;
    409:   MatScalar         *ap, *aa = a->a, *bap;

    411:   PetscFunctionBegin;
    412:   if (col < row) {
    413:     PetscCheck(a->ignore_ltriangular, PETSC_COMM_SELF, PETSC_ERR_USER, "Lower triangular value cannot be set for sbaij format. Ignoring these values, run with -mat_ignore_lower_triangular or call MatSetOption(mat,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE)");
    414:     PetscFunctionReturn(PETSC_SUCCESS); /* ignore lower triangular block */
    415:   }
    416:   rp    = aj + ai[row];
    417:   ap    = aa + bs2 * ai[row];
    418:   rmax  = imax[row];
    419:   nrow  = ailen[row];
    420:   value = v;
    421:   low   = 0;
    422:   high  = nrow;

    424:   while (high - low > 7) {
    425:     t = (low + high) / 2;
    426:     if (rp[t] > col) high = t;
    427:     else low = t;
    428:   }
    429:   for (i = low; i < high; i++) {
    430:     if (rp[i] > col) break;
    431:     if (rp[i] == col) {
    432:       bap = ap + bs2 * i;
    433:       if (roworiented) {
    434:         if (is == ADD_VALUES) {
    435:           for (ii = 0; ii < bs; ii++) {
    436:             for (jj = ii; jj < bs2; jj += bs) bap[jj] += *value++;
    437:           }
    438:         } else {
    439:           for (ii = 0; ii < bs; ii++) {
    440:             for (jj = ii; jj < bs2; jj += bs) bap[jj] = *value++;
    441:           }
    442:         }
    443:       } else {
    444:         if (is == ADD_VALUES) {
    445:           for (ii = 0; ii < bs; ii++) {
    446:             for (jj = 0; jj < bs; jj++) *bap++ += *value++;
    447:           }
    448:         } else {
    449:           for (ii = 0; ii < bs; ii++) {
    450:             for (jj = 0; jj < bs; jj++) *bap++ = *value++;
    451:           }
    452:         }
    453:       }
    454:       goto noinsert2;
    455:     }
    456:   }
    457:   if (nonew == 1) goto noinsert2;
    458:   PetscCheck(nonew != -1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Inserting a new block index nonzero block (%" PetscInt_FMT ", %" PetscInt_FMT ") in the matrix", orow, ocol);
    459:   MatSeqXAIJReallocateAIJ(A, a->mbs, bs2, nrow, row, col, rmax, aa, ai, aj, rp, ap, imax, nonew, MatScalar);
    460:   N = nrow++ - 1;
    461:   high++;
    462:   /* shift up all the later entries in this row */
    463:   PetscCall(PetscArraymove(rp + i + 1, rp + i, N - i + 1));
    464:   PetscCall(PetscArraymove(ap + bs2 * (i + 1), ap + bs2 * i, bs2 * (N - i + 1)));
    465:   rp[i] = col;
    466:   bap   = ap + bs2 * i;
    467:   if (roworiented) {
    468:     for (ii = 0; ii < bs; ii++) {
    469:       for (jj = ii; jj < bs2; jj += bs) bap[jj] = *value++;
    470:     }
    471:   } else {
    472:     for (ii = 0; ii < bs; ii++) {
    473:       for (jj = 0; jj < bs; jj++) *bap++ = *value++;
    474:     }
    475:   }
    476: noinsert2:;
    477:   ailen[row] = nrow;
    478:   PetscFunctionReturn(PETSC_SUCCESS);
    479: }

    481: /*
    482:    This routine is exactly duplicated in mpibaij.c
    483: */
    484: static inline PetscErrorCode MatSetValuesBlocked_SeqBAIJ_Inlined(Mat A, PetscInt row, PetscInt col, const PetscScalar v[], InsertMode is, PetscInt orow, PetscInt ocol)
    485: {
    486:   Mat_SeqBAIJ       *a = (Mat_SeqBAIJ *)A->data;
    487:   PetscInt          *rp, low, high, t, ii, jj, nrow, i, rmax, N;
    488:   PetscInt          *imax = a->imax, *ai = a->i, *ailen = a->ilen;
    489:   PetscInt          *aj = a->j, nonew = a->nonew, bs2 = a->bs2, bs = A->rmap->bs;
    490:   PetscBool          roworiented = a->roworiented;
    491:   const PetscScalar *value       = v;
    492:   MatScalar         *ap, *aa = a->a, *bap;

    494:   PetscFunctionBegin;
    495:   rp    = aj + ai[row];
    496:   ap    = aa + bs2 * ai[row];
    497:   rmax  = imax[row];
    498:   nrow  = ailen[row];
    499:   low   = 0;
    500:   high  = nrow;
    501:   value = v;
    502:   while (high - low > 7) {
    503:     t = (low + high) / 2;
    504:     if (rp[t] > col) high = t;
    505:     else low = t;
    506:   }
    507:   for (i = low; i < high; i++) {
    508:     if (rp[i] > col) break;
    509:     if (rp[i] == col) {
    510:       bap = ap + bs2 * i;
    511:       if (roworiented) {
    512:         if (is == ADD_VALUES) {
    513:           for (ii = 0; ii < bs; ii++) {
    514:             for (jj = ii; jj < bs2; jj += bs) bap[jj] += *value++;
    515:           }
    516:         } else {
    517:           for (ii = 0; ii < bs; ii++) {
    518:             for (jj = ii; jj < bs2; jj += bs) bap[jj] = *value++;
    519:           }
    520:         }
    521:       } else {
    522:         if (is == ADD_VALUES) {
    523:           for (ii = 0; ii < bs; ii++, value += bs) {
    524:             for (jj = 0; jj < bs; jj++) bap[jj] += value[jj];
    525:             bap += bs;
    526:           }
    527:         } else {
    528:           for (ii = 0; ii < bs; ii++, value += bs) {
    529:             for (jj = 0; jj < bs; jj++) bap[jj] = value[jj];
    530:             bap += bs;
    531:           }
    532:         }
    533:       }
    534:       goto noinsert2;
    535:     }
    536:   }
    537:   if (nonew == 1) goto noinsert2;
    538:   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);
    539:   MatSeqXAIJReallocateAIJ(A, a->mbs, bs2, nrow, row, col, rmax, aa, ai, aj, rp, ap, imax, nonew, MatScalar);
    540:   N = nrow++ - 1;
    541:   high++;
    542:   /* shift up all the later entries in this row */
    543:   PetscCall(PetscArraymove(rp + i + 1, rp + i, N - i + 1));
    544:   PetscCall(PetscArraymove(ap + bs2 * (i + 1), ap + bs2 * i, bs2 * (N - i + 1)));
    545:   rp[i] = col;
    546:   bap   = ap + bs2 * i;
    547:   if (roworiented) {
    548:     for (ii = 0; ii < bs; ii++) {
    549:       for (jj = ii; jj < bs2; jj += bs) bap[jj] = *value++;
    550:     }
    551:   } else {
    552:     for (ii = 0; ii < bs; ii++) {
    553:       for (jj = 0; jj < bs; jj++) *bap++ = *value++;
    554:     }
    555:   }
    556: noinsert2:;
    557:   ailen[row] = nrow;
    558:   PetscFunctionReturn(PETSC_SUCCESS);
    559: }

    561: /*
    562:     This routine could be optimized by removing the need for the block copy below and passing stride information
    563:   to the above inline routines; similarly in MatSetValuesBlocked_MPIBAIJ()
    564: */
    565: static PetscErrorCode MatSetValuesBlocked_MPISBAIJ(Mat mat, PetscInt m, const PetscInt im[], PetscInt n, const PetscInt in[], const MatScalar v[], InsertMode addv)
    566: {
    567:   Mat_MPISBAIJ    *baij = (Mat_MPISBAIJ *)mat->data;
    568:   const MatScalar *value;
    569:   MatScalar       *barray      = baij->barray;
    570:   PetscBool        roworiented = baij->roworiented, ignore_ltriangular = ((Mat_SeqSBAIJ *)baij->A->data)->ignore_ltriangular;
    571:   PetscInt         i, j, ii, jj, row, col, rstart = baij->rstartbs;
    572:   PetscInt         rend = baij->rendbs, cstart = baij->cstartbs, stepval;
    573:   PetscInt         cend = baij->cendbs, bs = mat->rmap->bs, bs2 = baij->bs2;

    575:   PetscFunctionBegin;
    576:   if (!barray) {
    577:     PetscCall(PetscMalloc1(bs2, &barray));
    578:     baij->barray = barray;
    579:   }

    581:   if (roworiented) {
    582:     stepval = (n - 1) * bs;
    583:   } else {
    584:     stepval = (m - 1) * bs;
    585:   }
    586:   for (i = 0; i < m; i++) {
    587:     if (im[i] < 0) continue;
    588:     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);
    589:     if (im[i] >= rstart && im[i] < rend) {
    590:       row = im[i] - rstart;
    591:       for (j = 0; j < n; j++) {
    592:         if (im[i] > in[j]) {
    593:           PetscCheck(ignore_ltriangular, PETSC_COMM_SELF, PETSC_ERR_USER, "Lower triangular value cannot be set for sbaij format. Ignoring these values, run with -mat_ignore_lower_triangular or call MatSetOption(mat,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE)");
    594:           continue; /* ignore lower triangular blocks */
    595:         }
    596:         /* If NumCol = 1 then a copy is not required */
    597:         if (roworiented && n == 1) {
    598:           barray = (MatScalar *)v + i * bs2;
    599:         } else if ((!roworiented) && (m == 1)) {
    600:           barray = (MatScalar *)v + j * bs2;
    601:         } else { /* Here a copy is required */
    602:           if (roworiented) {
    603:             value = v + i * (stepval + bs) * bs + j * bs;
    604:           } else {
    605:             value = v + j * (stepval + bs) * bs + i * bs;
    606:           }
    607:           for (ii = 0; ii < bs; ii++, value += stepval) {
    608:             for (jj = 0; jj < bs; jj++) *barray++ = *value++;
    609:           }
    610:           barray -= bs2;
    611:         }

    613:         if (in[j] >= cstart && in[j] < cend) {
    614:           col = in[j] - cstart;
    615:           PetscCall(MatSetValuesBlocked_SeqSBAIJ_Inlined(baij->A, row, col, barray, addv, im[i], in[j]));
    616:         } else if (in[j] < 0) {
    617:           continue;
    618:         } else {
    619:           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);
    620:           if (mat->was_assembled) {
    621:             if (!baij->colmap) PetscCall(MatCreateColmap_MPIBAIJ_Private(mat));

    623: #if defined(PETSC_USE_CTABLE)
    624:             PetscCall(PetscHMapIGetWithDefault(baij->colmap, in[j] + 1, 0, &col));
    625:             col = col < 1 ? -1 : (col - 1) / bs;
    626: #else
    627:             col = baij->colmap[in[j]] < 1 ? -1 : (baij->colmap[in[j]] - 1) / bs;
    628: #endif
    629:             if (col < 0 && !((Mat_SeqBAIJ *)baij->A->data)->nonew) {
    630:               PetscCall(MatDisAssemble_MPISBAIJ(mat));
    631:               col = in[j];
    632:             }
    633:           } else col = in[j];
    634:           PetscCall(MatSetValuesBlocked_SeqBAIJ_Inlined(baij->B, row, col, barray, addv, im[i], in[j]));
    635:         }
    636:       }
    637:     } else {
    638:       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]);
    639:       if (!baij->donotstash) {
    640:         if (roworiented) {
    641:           PetscCall(MatStashValuesRowBlocked_Private(&mat->bstash, im[i], n, in, v, m, n, i));
    642:         } else {
    643:           PetscCall(MatStashValuesColBlocked_Private(&mat->bstash, im[i], n, in, v, m, n, i));
    644:         }
    645:       }
    646:     }
    647:   }
    648:   PetscFunctionReturn(PETSC_SUCCESS);
    649: }

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

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

    688: static PetscErrorCode MatNorm_MPISBAIJ(Mat mat, NormType type, PetscReal *norm)
    689: {
    690:   Mat_MPISBAIJ *baij = (Mat_MPISBAIJ *)mat->data;
    691:   PetscReal     sum[2], *lnorm2;

    693:   PetscFunctionBegin;
    694:   if (baij->size == 1) {
    695:     PetscCall(MatNorm(baij->A, type, norm));
    696:   } else {
    697:     if (type == NORM_FROBENIUS) {
    698:       PetscCall(PetscMalloc1(2, &lnorm2));
    699:       PetscCall(MatNorm(baij->A, type, lnorm2));
    700:       *lnorm2 = (*lnorm2) * (*lnorm2);
    701:       lnorm2++; /* square power of norm(A) */
    702:       PetscCall(MatNorm(baij->B, type, lnorm2));
    703:       *lnorm2 = (*lnorm2) * (*lnorm2);
    704:       lnorm2--; /* square power of norm(B) */
    705:       PetscCallMPI(MPIU_Allreduce(lnorm2, sum, 2, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)mat)));
    706:       *norm = PetscSqrtReal(sum[0] + 2 * sum[1]);
    707:       PetscCall(PetscFree(lnorm2));
    708:     } else if (type == NORM_INFINITY || type == NORM_1) { /* max row/column sum */
    709:       Mat_SeqSBAIJ *amat = (Mat_SeqSBAIJ *)baij->A->data;
    710:       Mat_SeqBAIJ  *bmat = (Mat_SeqBAIJ *)baij->B->data;
    711:       PetscReal    *rsum, vabs;
    712:       PetscInt     *jj, *garray = baij->garray, rstart = baij->rstartbs, nz;
    713:       PetscInt      brow, bcol, col, bs = baij->A->rmap->bs, row, grow, gcol, mbs = amat->mbs;
    714:       MatScalar    *v;

    716:       PetscCall(PetscCalloc1(mat->cmap->N, &rsum));
    717:       /* Amat */
    718:       v  = amat->a;
    719:       jj = amat->j;
    720:       for (brow = 0; brow < mbs; brow++) {
    721:         grow = bs * (rstart + brow);
    722:         nz   = amat->i[brow + 1] - amat->i[brow];
    723:         for (bcol = 0; bcol < nz; bcol++) {
    724:           gcol = bs * (rstart + *jj);
    725:           jj++;
    726:           for (col = 0; col < bs; col++) {
    727:             for (row = 0; row < bs; row++) {
    728:               vabs = PetscAbsScalar(*v);
    729:               v++;
    730:               rsum[gcol + col] += vabs;
    731:               /* non-diagonal block */
    732:               if (bcol > 0 && vabs > 0.0) rsum[grow + row] += vabs;
    733:             }
    734:           }
    735:         }
    736:         PetscCall(PetscLogFlops(nz * bs * bs));
    737:       }
    738:       /* Bmat */
    739:       v  = bmat->a;
    740:       jj = bmat->j;
    741:       for (brow = 0; brow < mbs; brow++) {
    742:         grow = bs * (rstart + brow);
    743:         nz   = bmat->i[brow + 1] - bmat->i[brow];
    744:         for (bcol = 0; bcol < nz; bcol++) {
    745:           gcol = bs * garray[*jj];
    746:           jj++;
    747:           for (col = 0; col < bs; col++) {
    748:             for (row = 0; row < bs; row++) {
    749:               vabs = PetscAbsScalar(*v);
    750:               v++;
    751:               rsum[gcol + col] += vabs;
    752:               rsum[grow + row] += vabs;
    753:             }
    754:           }
    755:         }
    756:         PetscCall(PetscLogFlops(nz * bs * bs));
    757:       }
    758:       PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, rsum, mat->cmap->N, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)mat)));
    759:       *norm = 0.0;
    760:       for (col = 0; col < mat->cmap->N; col++) {
    761:         if (rsum[col] > *norm) *norm = rsum[col];
    762:       }
    763:       PetscCall(PetscFree(rsum));
    764:     } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for this norm yet");
    765:   }
    766:   PetscFunctionReturn(PETSC_SUCCESS);
    767: }

    769: static PetscErrorCode MatAssemblyBegin_MPISBAIJ(Mat mat, MatAssemblyType mode)
    770: {
    771:   Mat_MPISBAIJ *baij = (Mat_MPISBAIJ *)mat->data;
    772:   PetscInt      nstash, reallocs;

    774:   PetscFunctionBegin;
    775:   if (baij->donotstash || mat->nooffprocentries) PetscFunctionReturn(PETSC_SUCCESS);

    777:   PetscCall(MatStashScatterBegin_Private(mat, &mat->stash, mat->rmap->range));
    778:   PetscCall(MatStashScatterBegin_Private(mat, &mat->bstash, baij->rangebs));
    779:   PetscCall(MatStashGetInfo_Private(&mat->stash, &nstash, &reallocs));
    780:   PetscCall(PetscInfo(mat, "Stash has %" PetscInt_FMT " entries,uses %" PetscInt_FMT " mallocs.\n", nstash, reallocs));
    781:   PetscCall(MatStashGetInfo_Private(&mat->stash, &nstash, &reallocs));
    782:   PetscCall(PetscInfo(mat, "Block-Stash has %" PetscInt_FMT " entries, uses %" PetscInt_FMT " mallocs.\n", nstash, reallocs));
    783:   PetscFunctionReturn(PETSC_SUCCESS);
    784: }

    786: static PetscErrorCode MatAssemblyEnd_MPISBAIJ(Mat mat, MatAssemblyType mode)
    787: {
    788:   Mat_MPISBAIJ *baij = (Mat_MPISBAIJ *)mat->data;
    789:   Mat_SeqSBAIJ *a    = (Mat_SeqSBAIJ *)baij->A->data;
    790:   PetscInt      i, j, rstart, ncols, flg, bs2 = baij->bs2;
    791:   PetscInt     *row, *col;
    792:   PetscBool     all_assembled;
    793:   PetscMPIInt   n;
    794:   PetscBool     r1, r2, r3;
    795:   MatScalar    *val;

    797:   /* do not use 'b=(Mat_SeqBAIJ*)baij->B->data' as B can be reset in disassembly */
    798:   PetscFunctionBegin;
    799:   if (!baij->donotstash && !mat->nooffprocentries) {
    800:     while (1) {
    801:       PetscCall(MatStashScatterGetMesg_Private(&mat->stash, &n, &row, &col, &val, &flg));
    802:       if (!flg) break;

    804:       for (i = 0; i < n;) {
    805:         /* Now identify the consecutive vals belonging to the same row */
    806:         for (j = i, rstart = row[j]; j < n; j++) {
    807:           if (row[j] != rstart) break;
    808:         }
    809:         if (j < n) ncols = j - i;
    810:         else ncols = n - i;
    811:         /* Now assemble all these values with a single function call */
    812:         PetscCall(MatSetValues_MPISBAIJ(mat, 1, row + i, ncols, col + i, val + i, mat->insertmode));
    813:         i = j;
    814:       }
    815:     }
    816:     PetscCall(MatStashScatterEnd_Private(&mat->stash));
    817:     /* Now process the block-stash. Since the values are stashed column-oriented,
    818:        set the row-oriented flag to column-oriented, and after MatSetValues()
    819:        restore the original flags */
    820:     r1 = baij->roworiented;
    821:     r2 = a->roworiented;
    822:     r3 = ((Mat_SeqBAIJ *)baij->B->data)->roworiented;

    824:     baij->roworiented = PETSC_FALSE;
    825:     a->roworiented    = PETSC_FALSE;

    827:     ((Mat_SeqBAIJ *)baij->B->data)->roworiented = PETSC_FALSE; /* b->roworiented */
    828:     while (1) {
    829:       PetscCall(MatStashScatterGetMesg_Private(&mat->bstash, &n, &row, &col, &val, &flg));
    830:       if (!flg) break;

    832:       for (i = 0; i < n;) {
    833:         /* Now identify the consecutive vals belonging to the same row */
    834:         for (j = i, rstart = row[j]; j < n; j++) {
    835:           if (row[j] != rstart) break;
    836:         }
    837:         if (j < n) ncols = j - i;
    838:         else ncols = n - i;
    839:         PetscCall(MatSetValuesBlocked_MPISBAIJ(mat, 1, row + i, ncols, col + i, val + i * bs2, mat->insertmode));
    840:         i = j;
    841:       }
    842:     }
    843:     PetscCall(MatStashScatterEnd_Private(&mat->bstash));

    845:     baij->roworiented = r1;
    846:     a->roworiented    = r2;

    848:     ((Mat_SeqBAIJ *)baij->B->data)->roworiented = r3; /* b->roworiented */
    849:   }

    851:   PetscCall(MatAssemblyBegin(baij->A, mode));
    852:   PetscCall(MatAssemblyEnd(baij->A, mode));

    854:   /* determine if any process has disassembled, if so we must
    855:      also disassemble ourselves, in order that we may reassemble. */
    856:   /*
    857:      if nonzero structure of submatrix B cannot change then we know that
    858:      no process disassembled thus we can skip this stuff
    859:   */
    860:   if (!((Mat_SeqBAIJ *)baij->B->data)->nonew) {
    861:     PetscCallMPI(MPIU_Allreduce(&mat->was_assembled, &all_assembled, 1, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)mat)));
    862:     if (mat->was_assembled && !all_assembled) PetscCall(MatDisAssemble_MPISBAIJ(mat));
    863:   }

    865:   if (!mat->was_assembled && mode == MAT_FINAL_ASSEMBLY) PetscCall(MatSetUpMultiply_MPISBAIJ(mat)); /* setup Mvctx and sMvctx */
    866:   PetscCall(MatAssemblyBegin(baij->B, mode));
    867:   PetscCall(MatAssemblyEnd(baij->B, mode));

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

    871:   baij->rowvalues = NULL;

    873:   /* if no new nonzero locations are allowed in matrix then only set the matrix state the first time through */
    874:   if ((!mat->was_assembled && mode == MAT_FINAL_ASSEMBLY) || !((Mat_SeqBAIJ *)baij->A->data)->nonew) {
    875:     PetscObjectState state = baij->A->nonzerostate + baij->B->nonzerostate;
    876:     PetscCallMPI(MPIU_Allreduce(&state, &mat->nonzerostate, 1, MPIU_INT64, MPI_SUM, PetscObjectComm((PetscObject)mat)));
    877:   }
    878:   PetscFunctionReturn(PETSC_SUCCESS);
    879: }

    881: extern PetscErrorCode MatSetValues_MPIBAIJ(Mat, PetscInt, const PetscInt[], PetscInt, const PetscInt[], const PetscScalar[], InsertMode);
    882: #include <petscdraw.h>
    883: static PetscErrorCode MatView_MPISBAIJ_ASCIIorDraworSocket(Mat mat, PetscViewer viewer)
    884: {
    885:   Mat_MPISBAIJ     *baij = (Mat_MPISBAIJ *)mat->data;
    886:   PetscInt          bs   = mat->rmap->bs;
    887:   PetscMPIInt       rank = baij->rank;
    888:   PetscBool         isascii, isdraw;
    889:   PetscViewer       sviewer;
    890:   PetscViewerFormat format;

    892:   PetscFunctionBegin;
    893:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    894:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
    895:   if (isascii) {
    896:     PetscCall(PetscViewerGetFormat(viewer, &format));
    897:     if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
    898:       MatInfo info;
    899:       PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)mat), &rank));
    900:       PetscCall(MatGetInfo(mat, MAT_LOCAL, &info));
    901:       PetscCall(PetscViewerASCIIPushSynchronized(viewer));
    902:       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,
    903:                                                    mat->rmap->bs, info.memory));
    904:       PetscCall(MatGetInfo(baij->A, MAT_LOCAL, &info));
    905:       PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] on-diagonal part: nz %" PetscInt_FMT " \n", rank, (PetscInt)info.nz_used));
    906:       PetscCall(MatGetInfo(baij->B, MAT_LOCAL, &info));
    907:       PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] off-diagonal part: nz %" PetscInt_FMT " \n", rank, (PetscInt)info.nz_used));
    908:       PetscCall(PetscViewerFlush(viewer));
    909:       PetscCall(PetscViewerASCIIPopSynchronized(viewer));
    910:       PetscCall(PetscViewerASCIIPrintf(viewer, "Information on VecScatter used in matrix-vector product: \n"));
    911:       PetscCall(VecScatterView(baij->Mvctx, viewer));
    912:       PetscFunctionReturn(PETSC_SUCCESS);
    913:     } else if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_FACTOR_INFO) PetscFunctionReturn(PETSC_SUCCESS);
    914:   }

    916:   if (isdraw) {
    917:     PetscDraw draw;
    918:     PetscBool isnull;
    919:     PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
    920:     PetscCall(PetscDrawIsNull(draw, &isnull));
    921:     if (isnull) PetscFunctionReturn(PETSC_SUCCESS);
    922:   }

    924:   {
    925:     /* assemble the entire matrix onto first processor. */
    926:     Mat           A;
    927:     Mat_SeqSBAIJ *Aloc;
    928:     Mat_SeqBAIJ  *Bloc;
    929:     PetscInt      M = mat->rmap->N, N = mat->cmap->N, *ai, *aj, col, i, j, k, *rvals, mbs = baij->mbs;
    930:     MatScalar    *a;
    931:     const char   *matname;

    933:     /* Should this be the same type as mat? */
    934:     PetscCall(MatCreate(PetscObjectComm((PetscObject)mat), &A));
    935:     if (rank == 0) {
    936:       PetscCall(MatSetSizes(A, M, N, M, N));
    937:     } else {
    938:       PetscCall(MatSetSizes(A, 0, 0, M, N));
    939:     }
    940:     PetscCall(MatSetType(A, MATMPISBAIJ));
    941:     PetscCall(MatMPISBAIJSetPreallocation(A, mat->rmap->bs, 0, NULL, 0, NULL));
    942:     PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_FALSE));

    944:     /* copy over the A part */
    945:     Aloc = (Mat_SeqSBAIJ *)baij->A->data;
    946:     ai   = Aloc->i;
    947:     aj   = Aloc->j;
    948:     a    = Aloc->a;
    949:     PetscCall(PetscMalloc1(bs, &rvals));

    951:     for (i = 0; i < mbs; i++) {
    952:       rvals[0] = bs * (baij->rstartbs + i);
    953:       for (j = 1; j < bs; j++) rvals[j] = rvals[j - 1] + 1;
    954:       for (j = ai[i]; j < ai[i + 1]; j++) {
    955:         col = (baij->cstartbs + aj[j]) * bs;
    956:         for (k = 0; k < bs; k++) {
    957:           PetscCall(MatSetValues_MPISBAIJ(A, bs, rvals, 1, &col, a, INSERT_VALUES));
    958:           col++;
    959:           a += bs;
    960:         }
    961:       }
    962:     }
    963:     /* copy over the B part */
    964:     Bloc = (Mat_SeqBAIJ *)baij->B->data;
    965:     ai   = Bloc->i;
    966:     aj   = Bloc->j;
    967:     a    = Bloc->a;
    968:     for (i = 0; i < mbs; i++) {
    969:       rvals[0] = bs * (baij->rstartbs + i);
    970:       for (j = 1; j < bs; j++) rvals[j] = rvals[j - 1] + 1;
    971:       for (j = ai[i]; j < ai[i + 1]; j++) {
    972:         col = baij->garray[aj[j]] * bs;
    973:         for (k = 0; k < bs; k++) {
    974:           PetscCall(MatSetValues_MPIBAIJ(A, bs, rvals, 1, &col, a, INSERT_VALUES));
    975:           col++;
    976:           a += bs;
    977:         }
    978:       }
    979:     }
    980:     PetscCall(PetscFree(rvals));
    981:     PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    982:     PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    983:     /*
    984:        Everyone has to call to draw the matrix since the graphics waits are
    985:        synchronized across all processors that share the PetscDraw object
    986:     */
    987:     PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
    988:     if (((PetscObject)mat)->name) PetscCall(PetscObjectGetName((PetscObject)mat, &matname));
    989:     if (rank == 0) {
    990:       if (((PetscObject)mat)->name) PetscCall(PetscObjectSetName((PetscObject)((Mat_MPISBAIJ *)A->data)->A, matname));
    991:       PetscCall(MatView_SeqSBAIJ(((Mat_MPISBAIJ *)A->data)->A, sviewer));
    992:     }
    993:     PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
    994:     PetscCall(MatDestroy(&A));
    995:   }
    996:   PetscFunctionReturn(PETSC_SUCCESS);
    997: }

    999: /* Used for both MPIBAIJ and MPISBAIJ matrices */
    1000: #define MatView_MPISBAIJ_Binary MatView_MPIBAIJ_Binary

    1002: static PetscErrorCode MatView_MPISBAIJ(Mat mat, PetscViewer viewer)
    1003: {
    1004:   PetscBool isascii, isdraw, issocket, isbinary;

    1006:   PetscFunctionBegin;
    1007:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    1008:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
    1009:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSOCKET, &issocket));
    1010:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
    1011:   if (isascii || isdraw || issocket) PetscCall(MatView_MPISBAIJ_ASCIIorDraworSocket(mat, viewer));
    1012:   else if (isbinary) PetscCall(MatView_MPISBAIJ_Binary(mat, viewer));
    1013:   PetscFunctionReturn(PETSC_SUCCESS);
    1014: }

    1016: #if defined(PETSC_USE_COMPLEX)
    1017: static PetscErrorCode MatMult_MPISBAIJ_Hermitian(Mat A, Vec xx, Vec yy)
    1018: {
    1019:   Mat_MPISBAIJ      *a   = (Mat_MPISBAIJ *)A->data;
    1020:   PetscInt           mbs = a->mbs, bs = A->rmap->bs;
    1021:   PetscScalar       *from;
    1022:   const PetscScalar *x;

    1024:   PetscFunctionBegin;
    1025:   /* diagonal part */
    1026:   PetscCall((*a->A->ops->mult)(a->A, xx, a->slvec1a));
    1027:   /* since a->slvec1b shares memory (dangerously) with a->slec1 changes to a->slec1 will affect it */
    1028:   PetscCall(PetscObjectStateIncrease((PetscObject)a->slvec1b));
    1029:   PetscCall(VecZeroEntries(a->slvec1b));

    1031:   /* subdiagonal part */
    1032:   PetscCheck(a->B->ops->multhermitiantranspose, PetscObjectComm((PetscObject)a->B), PETSC_ERR_SUP, "Not for type %s", ((PetscObject)a->B)->type_name);
    1033:   PetscCall((*a->B->ops->multhermitiantranspose)(a->B, xx, a->slvec0b));

    1035:   /* copy x into the vec slvec0 */
    1036:   PetscCall(VecGetArray(a->slvec0, &from));
    1037:   PetscCall(VecGetArrayRead(xx, &x));

    1039:   PetscCall(PetscArraycpy(from, x, bs * mbs));
    1040:   PetscCall(VecRestoreArray(a->slvec0, &from));
    1041:   PetscCall(VecRestoreArrayRead(xx, &x));

    1043:   PetscCall(VecScatterBegin(a->sMvctx, a->slvec0, a->slvec1, ADD_VALUES, SCATTER_FORWARD));
    1044:   PetscCall(VecScatterEnd(a->sMvctx, a->slvec0, a->slvec1, ADD_VALUES, SCATTER_FORWARD));
    1045:   /* supperdiagonal part */
    1046:   PetscCall((*a->B->ops->multadd)(a->B, a->slvec1b, a->slvec1a, yy));
    1047:   PetscFunctionReturn(PETSC_SUCCESS);
    1048: }
    1049: #endif

    1051: static PetscErrorCode MatMult_MPISBAIJ(Mat A, Vec xx, Vec yy)
    1052: {
    1053:   Mat_MPISBAIJ      *a   = (Mat_MPISBAIJ *)A->data;
    1054:   PetscInt           mbs = a->mbs, bs = A->rmap->bs;
    1055:   PetscScalar       *from;
    1056:   const PetscScalar *x;

    1058:   PetscFunctionBegin;
    1059:   /* diagonal part */
    1060:   PetscCall((*a->A->ops->mult)(a->A, xx, a->slvec1a));
    1061:   /* since a->slvec1b shares memory (dangerously) with a->slec1 changes to a->slec1 will affect it */
    1062:   PetscCall(PetscObjectStateIncrease((PetscObject)a->slvec1b));
    1063:   PetscCall(VecZeroEntries(a->slvec1b));

    1065:   /* subdiagonal part */
    1066:   PetscCall((*a->B->ops->multtranspose)(a->B, xx, a->slvec0b));

    1068:   /* copy x into the vec slvec0 */
    1069:   PetscCall(VecGetArray(a->slvec0, &from));
    1070:   PetscCall(VecGetArrayRead(xx, &x));

    1072:   PetscCall(PetscArraycpy(from, x, bs * mbs));
    1073:   PetscCall(VecRestoreArray(a->slvec0, &from));
    1074:   PetscCall(VecRestoreArrayRead(xx, &x));

    1076:   PetscCall(VecScatterBegin(a->sMvctx, a->slvec0, a->slvec1, ADD_VALUES, SCATTER_FORWARD));
    1077:   PetscCall(VecScatterEnd(a->sMvctx, a->slvec0, a->slvec1, ADD_VALUES, SCATTER_FORWARD));
    1078:   /* supperdiagonal part */
    1079:   PetscCall((*a->B->ops->multadd)(a->B, a->slvec1b, a->slvec1a, yy));
    1080:   PetscFunctionReturn(PETSC_SUCCESS);
    1081: }

    1083: #if PetscDefined(USE_COMPLEX)
    1084: static PetscErrorCode MatMultAdd_MPISBAIJ_Hermitian(Mat A, Vec xx, Vec yy, Vec zz)
    1085: {
    1086:   Mat_MPISBAIJ      *a   = (Mat_MPISBAIJ *)A->data;
    1087:   PetscInt           mbs = a->mbs, bs = A->rmap->bs;
    1088:   PetscScalar       *from;
    1089:   const PetscScalar *x;

    1091:   PetscFunctionBegin;
    1092:   /* diagonal part */
    1093:   PetscCall((*a->A->ops->multadd)(a->A, xx, yy, a->slvec1a));
    1094:   PetscCall(PetscObjectStateIncrease((PetscObject)a->slvec1b));
    1095:   PetscCall(VecZeroEntries(a->slvec1b));

    1097:   /* subdiagonal part */
    1098:   PetscCheck(a->B->ops->multhermitiantranspose, PetscObjectComm((PetscObject)a->B), PETSC_ERR_SUP, "Not for type %s", ((PetscObject)a->B)->type_name);
    1099:   PetscCall((*a->B->ops->multhermitiantranspose)(a->B, xx, a->slvec0b));

    1101:   /* copy x into the vec slvec0 */
    1102:   PetscCall(VecGetArray(a->slvec0, &from));
    1103:   PetscCall(VecGetArrayRead(xx, &x));
    1104:   PetscCall(PetscArraycpy(from, x, bs * mbs));
    1105:   PetscCall(VecRestoreArray(a->slvec0, &from));

    1107:   PetscCall(VecScatterBegin(a->sMvctx, a->slvec0, a->slvec1, ADD_VALUES, SCATTER_FORWARD));
    1108:   PetscCall(VecRestoreArrayRead(xx, &x));
    1109:   PetscCall(VecScatterEnd(a->sMvctx, a->slvec0, a->slvec1, ADD_VALUES, SCATTER_FORWARD));

    1111:   /* supperdiagonal part */
    1112:   PetscCall((*a->B->ops->multadd)(a->B, a->slvec1b, a->slvec1a, zz));
    1113:   PetscFunctionReturn(PETSC_SUCCESS);
    1114: }
    1115: #endif

    1117: static PetscErrorCode MatMultAdd_MPISBAIJ(Mat A, Vec xx, Vec yy, Vec zz)
    1118: {
    1119:   Mat_MPISBAIJ      *a   = (Mat_MPISBAIJ *)A->data;
    1120:   PetscInt           mbs = a->mbs, bs = A->rmap->bs;
    1121:   PetscScalar       *from;
    1122:   const PetscScalar *x;

    1124:   PetscFunctionBegin;
    1125:   /* diagonal part */
    1126:   PetscCall((*a->A->ops->multadd)(a->A, xx, yy, a->slvec1a));
    1127:   PetscCall(PetscObjectStateIncrease((PetscObject)a->slvec1b));
    1128:   PetscCall(VecZeroEntries(a->slvec1b));

    1130:   /* subdiagonal part */
    1131:   PetscCall((*a->B->ops->multtranspose)(a->B, xx, a->slvec0b));

    1133:   /* copy x into the vec slvec0 */
    1134:   PetscCall(VecGetArray(a->slvec0, &from));
    1135:   PetscCall(VecGetArrayRead(xx, &x));
    1136:   PetscCall(PetscArraycpy(from, x, bs * mbs));
    1137:   PetscCall(VecRestoreArray(a->slvec0, &from));

    1139:   PetscCall(VecScatterBegin(a->sMvctx, a->slvec0, a->slvec1, ADD_VALUES, SCATTER_FORWARD));
    1140:   PetscCall(VecRestoreArrayRead(xx, &x));
    1141:   PetscCall(VecScatterEnd(a->sMvctx, a->slvec0, a->slvec1, ADD_VALUES, SCATTER_FORWARD));

    1143:   /* supperdiagonal part */
    1144:   PetscCall((*a->B->ops->multadd)(a->B, a->slvec1b, a->slvec1a, zz));
    1145:   PetscFunctionReturn(PETSC_SUCCESS);
    1146: }

    1148: /*
    1149:   This only works correctly for square matrices where the subblock A->A is the
    1150:    diagonal block
    1151: */
    1152: static PetscErrorCode MatGetDiagonal_MPISBAIJ(Mat A, Vec v)
    1153: {
    1154:   Mat_MPISBAIJ *a = (Mat_MPISBAIJ *)A->data;

    1156:   PetscFunctionBegin;
    1157:   /* PetscCheck(a->rmap->N == a->cmap->N,PETSC_COMM_SELF,PETSC_ERR_SUP,"Supports only square matrix where A->A is diag block"); */
    1158:   PetscCall(MatGetDiagonal(a->A, v));
    1159:   PetscFunctionReturn(PETSC_SUCCESS);
    1160: }

    1162: static PetscErrorCode MatScale_MPISBAIJ(Mat A, PetscScalar aa)
    1163: {
    1164:   Mat_MPISBAIJ *a = (Mat_MPISBAIJ *)A->data;

    1166:   PetscFunctionBegin;
    1167:   PetscCall(MatScale(a->A, aa));
    1168:   PetscCall(MatScale(a->B, aa));
    1169:   PetscFunctionReturn(PETSC_SUCCESS);
    1170: }

    1172: static PetscErrorCode MatGetRow_MPISBAIJ(Mat matin, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v)
    1173: {
    1174:   Mat_MPISBAIJ *mat = (Mat_MPISBAIJ *)matin->data;
    1175:   PetscScalar  *vworkA, *vworkB, **pvA, **pvB, *v_p;
    1176:   PetscInt      bs = matin->rmap->bs, bs2 = mat->bs2, i, *cworkA, *cworkB, **pcA, **pcB;
    1177:   PetscInt      nztot, nzA, nzB, lrow, brstart = matin->rmap->rstart, brend = matin->rmap->rend;
    1178:   PetscInt     *cmap, *idx_p, cstart = mat->rstartbs;

    1180:   PetscFunctionBegin;
    1181:   PetscCheck(!mat->getrowactive, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Already active");
    1182:   mat->getrowactive = PETSC_TRUE;

    1184:   if (!mat->rowvalues && (idx || v)) {
    1185:     /*
    1186:         allocate enough space to hold information from the longest row.
    1187:     */
    1188:     Mat_SeqSBAIJ *Aa  = (Mat_SeqSBAIJ *)mat->A->data;
    1189:     Mat_SeqBAIJ  *Ba  = (Mat_SeqBAIJ *)mat->B->data;
    1190:     PetscInt      max = 1, mbs = mat->mbs, tmp;
    1191:     for (i = 0; i < mbs; i++) {
    1192:       tmp = Aa->i[i + 1] - Aa->i[i] + Ba->i[i + 1] - Ba->i[i]; /* row length */
    1193:       if (max < tmp) max = tmp;
    1194:     }
    1195:     PetscCall(PetscMalloc2(max * bs2, &mat->rowvalues, max * bs2, &mat->rowindices));
    1196:   }

    1198:   PetscCheck(row >= brstart && row < brend, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only local rows");
    1199:   lrow = row - brstart; /* local row index */

    1201:   pvA = &vworkA;
    1202:   pcA = &cworkA;
    1203:   pvB = &vworkB;
    1204:   pcB = &cworkB;
    1205:   if (!v) {
    1206:     pvA = NULL;
    1207:     pvB = NULL;
    1208:   }
    1209:   if (!idx) {
    1210:     pcA = NULL;
    1211:     if (!v) pcB = NULL;
    1212:   }
    1213:   PetscCall((*mat->A->ops->getrow)(mat->A, lrow, &nzA, pcA, pvA));
    1214:   PetscCall((*mat->B->ops->getrow)(mat->B, lrow, &nzB, pcB, pvB));
    1215:   nztot = nzA + nzB;

    1217:   cmap = mat->garray;
    1218:   if (v || idx) {
    1219:     if (nztot) {
    1220:       /* Sort by increasing column numbers, assuming A and B already sorted */
    1221:       PetscInt imark = -1;
    1222:       if (v) {
    1223:         *v = v_p = mat->rowvalues;
    1224:         for (i = 0; i < nzB; i++) {
    1225:           if (cmap[cworkB[i] / bs] < cstart) v_p[i] = vworkB[i];
    1226:           else break;
    1227:         }
    1228:         imark = i;
    1229:         for (i = 0; i < nzA; i++) v_p[imark + i] = vworkA[i];
    1230:         for (i = imark; i < nzB; i++) v_p[nzA + i] = vworkB[i];
    1231:       }
    1232:       if (idx) {
    1233:         *idx = idx_p = mat->rowindices;
    1234:         if (imark > -1) {
    1235:           for (i = 0; i < imark; i++) idx_p[i] = cmap[cworkB[i] / bs] * bs + cworkB[i] % bs;
    1236:         } else {
    1237:           for (i = 0; i < nzB; i++) {
    1238:             if (cmap[cworkB[i] / bs] < cstart) idx_p[i] = cmap[cworkB[i] / bs] * bs + cworkB[i] % bs;
    1239:             else break;
    1240:           }
    1241:           imark = i;
    1242:         }
    1243:         for (i = 0; i < nzA; i++) idx_p[imark + i] = cstart * bs + cworkA[i];
    1244:         for (i = imark; i < nzB; i++) idx_p[nzA + i] = cmap[cworkB[i] / bs] * bs + cworkB[i] % bs;
    1245:       }
    1246:     } else {
    1247:       if (idx) *idx = NULL;
    1248:       if (v) *v = NULL;
    1249:     }
    1250:   }
    1251:   *nz = nztot;
    1252:   PetscCall((*mat->A->ops->restorerow)(mat->A, lrow, &nzA, pcA, pvA));
    1253:   PetscCall((*mat->B->ops->restorerow)(mat->B, lrow, &nzB, pcB, pvB));
    1254:   PetscFunctionReturn(PETSC_SUCCESS);
    1255: }

    1257: static PetscErrorCode MatRestoreRow_MPISBAIJ(Mat mat, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v)
    1258: {
    1259:   Mat_MPISBAIJ *baij = (Mat_MPISBAIJ *)mat->data;

    1261:   PetscFunctionBegin;
    1262:   PetscCheck(baij->getrowactive, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "MatGetRow() must be called first");
    1263:   baij->getrowactive = PETSC_FALSE;
    1264:   PetscFunctionReturn(PETSC_SUCCESS);
    1265: }

    1267: static PetscErrorCode MatGetRowUpperTriangular_MPISBAIJ(Mat A)
    1268: {
    1269:   Mat_MPISBAIJ *a  = (Mat_MPISBAIJ *)A->data;
    1270:   Mat_SeqSBAIJ *aA = (Mat_SeqSBAIJ *)a->A->data;

    1272:   PetscFunctionBegin;
    1273:   aA->getrow_utriangular = PETSC_TRUE;
    1274:   PetscFunctionReturn(PETSC_SUCCESS);
    1275: }
    1276: static PetscErrorCode MatRestoreRowUpperTriangular_MPISBAIJ(Mat A)
    1277: {
    1278:   Mat_MPISBAIJ *a  = (Mat_MPISBAIJ *)A->data;
    1279:   Mat_SeqSBAIJ *aA = (Mat_SeqSBAIJ *)a->A->data;

    1281:   PetscFunctionBegin;
    1282:   aA->getrow_utriangular = PETSC_FALSE;
    1283:   PetscFunctionReturn(PETSC_SUCCESS);
    1284: }

    1286: static PetscErrorCode MatConjugate_MPISBAIJ(Mat mat)
    1287: {
    1288:   Mat_MPISBAIJ *a = (Mat_MPISBAIJ *)mat->data;

    1290:   PetscFunctionBegin;
    1291:   PetscCall(MatConjugate(a->A));
    1292:   PetscCall(MatConjugate(a->B));
    1293:   PetscFunctionReturn(PETSC_SUCCESS);
    1294: }

    1296: static PetscErrorCode MatRealPart_MPISBAIJ(Mat A)
    1297: {
    1298:   Mat_MPISBAIJ *a = (Mat_MPISBAIJ *)A->data;

    1300:   PetscFunctionBegin;
    1301:   PetscCall(MatRealPart(a->A));
    1302:   PetscCall(MatRealPart(a->B));
    1303:   PetscFunctionReturn(PETSC_SUCCESS);
    1304: }

    1306: static PetscErrorCode MatImaginaryPart_MPISBAIJ(Mat A)
    1307: {
    1308:   Mat_MPISBAIJ *a = (Mat_MPISBAIJ *)A->data;

    1310:   PetscFunctionBegin;
    1311:   PetscCall(MatImaginaryPart(a->A));
    1312:   PetscCall(MatImaginaryPart(a->B));
    1313:   PetscFunctionReturn(PETSC_SUCCESS);
    1314: }

    1316: /* Check if isrow is a subset of iscol_local, called by MatCreateSubMatrix_MPISBAIJ()
    1317:    Input: isrow       - distributed(parallel),
    1318:           iscol_local - locally owned (seq)
    1319: */
    1320: static PetscErrorCode ISEqual_private(IS isrow, IS iscol_local, PetscBool *flg)
    1321: {
    1322:   PetscInt        sz1, sz2, *a1, *a2, i, j, k, nmatch;
    1323:   const PetscInt *ptr1, *ptr2;

    1325:   PetscFunctionBegin;
    1326:   *flg = PETSC_FALSE;
    1327:   PetscCall(ISGetLocalSize(isrow, &sz1));
    1328:   PetscCall(ISGetLocalSize(iscol_local, &sz2));
    1329:   if (sz1 > sz2) PetscFunctionReturn(PETSC_SUCCESS);

    1331:   PetscCall(ISGetIndices(isrow, &ptr1));
    1332:   PetscCall(ISGetIndices(iscol_local, &ptr2));

    1334:   PetscCall(PetscMalloc1(sz1, &a1));
    1335:   PetscCall(PetscMalloc1(sz2, &a2));
    1336:   PetscCall(PetscArraycpy(a1, ptr1, sz1));
    1337:   PetscCall(PetscArraycpy(a2, ptr2, sz2));
    1338:   PetscCall(PetscSortInt(sz1, a1));
    1339:   PetscCall(PetscSortInt(sz2, a2));

    1341:   nmatch = 0;
    1342:   k      = 0;
    1343:   for (i = 0; i < sz1; i++) {
    1344:     for (j = k; j < sz2; j++) {
    1345:       if (a1[i] == a2[j]) {
    1346:         k = j;
    1347:         nmatch++;
    1348:         break;
    1349:       }
    1350:     }
    1351:   }
    1352:   PetscCall(ISRestoreIndices(isrow, &ptr1));
    1353:   PetscCall(ISRestoreIndices(iscol_local, &ptr2));
    1354:   PetscCall(PetscFree(a1));
    1355:   PetscCall(PetscFree(a2));
    1356:   if (nmatch < sz1) {
    1357:     *flg = PETSC_FALSE;
    1358:   } else {
    1359:     *flg = PETSC_TRUE;
    1360:   }
    1361:   PetscFunctionReturn(PETSC_SUCCESS);
    1362: }

    1364: static PetscErrorCode MatCreateSubMatrix_MPISBAIJ(Mat mat, IS isrow, IS iscol, MatReuse call, Mat *newmat)
    1365: {
    1366:   Mat       C[2];
    1367:   IS        iscol_local, isrow_local;
    1368:   PetscInt  csize, csize_local, rsize;
    1369:   PetscBool isequal, issorted, isidentity = PETSC_FALSE;

    1371:   PetscFunctionBegin;
    1372:   PetscCall(ISGetLocalSize(iscol, &csize));
    1373:   PetscCall(ISGetLocalSize(isrow, &rsize));
    1374:   if (call == MAT_REUSE_MATRIX) {
    1375:     PetscCall(PetscObjectQuery((PetscObject)*newmat, "ISAllGather", (PetscObject *)&iscol_local));
    1376:     PetscCheck(iscol_local, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Submatrix passed in was not used before, cannot reuse");
    1377:   } else {
    1378:     PetscCall(ISAllGather(iscol, &iscol_local));
    1379:     PetscCall(ISSorted(iscol_local, &issorted));
    1380:     PetscCheck(issorted, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "For symmetric format, iscol must be sorted");
    1381:   }
    1382:   PetscCall(ISEqual_private(isrow, iscol_local, &isequal));
    1383:   if (!isequal) {
    1384:     PetscCall(ISGetLocalSize(iscol_local, &csize_local));
    1385:     isidentity = (PetscBool)(mat->cmap->N == csize_local);
    1386:     if (!isidentity) {
    1387:       if (call == MAT_REUSE_MATRIX) {
    1388:         PetscCall(PetscObjectQuery((PetscObject)*newmat, "ISAllGather_other", (PetscObject *)&isrow_local));
    1389:         PetscCheck(isrow_local, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Submatrix passed in was not used before, cannot reuse");
    1390:       } else {
    1391:         PetscCall(ISAllGather(isrow, &isrow_local));
    1392:         PetscCall(ISSorted(isrow_local, &issorted));
    1393:         PetscCheck(issorted, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "For symmetric format, isrow must be sorted");
    1394:       }
    1395:     }
    1396:   }
    1397:   /* now call MatCreateSubMatrix_MPIBAIJ() */
    1398:   PetscCall(MatCreateSubMatrix_MPIBAIJ_Private(mat, isrow, iscol_local, csize, isequal || isidentity ? call : MAT_INITIAL_MATRIX, isequal || isidentity ? newmat : C, (PetscBool)(isequal || isidentity)));
    1399:   if (!isequal && !isidentity) {
    1400:     if (call == MAT_INITIAL_MATRIX) {
    1401:       IS       intersect;
    1402:       PetscInt ni;

    1404:       PetscCall(ISIntersect(isrow_local, iscol_local, &intersect));
    1405:       PetscCall(ISGetLocalSize(intersect, &ni));
    1406:       PetscCall(ISDestroy(&intersect));
    1407:       PetscCheck(ni == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot create such a submatrix: for symmetric format, when requesting an off-diagonal submatrix, isrow and iscol should have an empty intersection (number of common indices is %" PetscInt_FMT ")", ni);
    1408:     }
    1409:     PetscCall(MatCreateSubMatrix_MPIBAIJ_Private(mat, iscol, isrow_local, rsize, MAT_INITIAL_MATRIX, C + 1, PETSC_FALSE));
    1410:     PetscCall(MatTranspose(C[1], MAT_INPLACE_MATRIX, C + 1));
    1411:     PetscCall(MatAXPY(C[0], 1.0, C[1], DIFFERENT_NONZERO_PATTERN));
    1412:     if (call == MAT_REUSE_MATRIX) PetscCall(MatCopy(C[0], *newmat, SAME_NONZERO_PATTERN));
    1413:     else if (mat->rmap->bs == 1) PetscCall(MatConvert(C[0], MATAIJ, MAT_INITIAL_MATRIX, newmat));
    1414:     else PetscCall(MatCopy(C[0], *newmat, SAME_NONZERO_PATTERN));
    1415:     PetscCall(MatDestroy(C));
    1416:     PetscCall(MatDestroy(C + 1));
    1417:   }
    1418:   if (call == MAT_INITIAL_MATRIX) {
    1419:     if (!isequal && !isidentity) {
    1420:       PetscCall(PetscObjectCompose((PetscObject)*newmat, "ISAllGather_other", (PetscObject)isrow_local));
    1421:       PetscCall(ISDestroy(&isrow_local));
    1422:     }
    1423:     PetscCall(PetscObjectCompose((PetscObject)*newmat, "ISAllGather", (PetscObject)iscol_local));
    1424:     PetscCall(ISDestroy(&iscol_local));
    1425:   }
    1426:   PetscFunctionReturn(PETSC_SUCCESS);
    1427: }

    1429: static PetscErrorCode MatZeroEntries_MPISBAIJ(Mat A)
    1430: {
    1431:   Mat_MPISBAIJ *l = (Mat_MPISBAIJ *)A->data;

    1433:   PetscFunctionBegin;
    1434:   PetscCall(MatZeroEntries(l->A));
    1435:   PetscCall(MatZeroEntries(l->B));
    1436:   PetscFunctionReturn(PETSC_SUCCESS);
    1437: }

    1439: static PetscErrorCode MatGetInfo_MPISBAIJ(Mat matin, MatInfoType flag, MatInfo *info)
    1440: {
    1441:   Mat_MPISBAIJ  *a = (Mat_MPISBAIJ *)matin->data;
    1442:   Mat            A = a->A, B = a->B;
    1443:   PetscLogDouble isend[5], irecv[5];

    1445:   PetscFunctionBegin;
    1446:   info->block_size = (PetscReal)matin->rmap->bs;

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

    1450:   isend[0] = info->nz_used;
    1451:   isend[1] = info->nz_allocated;
    1452:   isend[2] = info->nz_unneeded;
    1453:   isend[3] = info->memory;
    1454:   isend[4] = info->mallocs;

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

    1458:   isend[0] += info->nz_used;
    1459:   isend[1] += info->nz_allocated;
    1460:   isend[2] += info->nz_unneeded;
    1461:   isend[3] += info->memory;
    1462:   isend[4] += info->mallocs;
    1463:   if (flag == MAT_LOCAL) {
    1464:     info->nz_used      = isend[0];
    1465:     info->nz_allocated = isend[1];
    1466:     info->nz_unneeded  = isend[2];
    1467:     info->memory       = isend[3];
    1468:     info->mallocs      = isend[4];
    1469:   } else if (flag == MAT_GLOBAL_MAX) {
    1470:     PetscCallMPI(MPIU_Allreduce(isend, irecv, 5, MPIU_PETSCLOGDOUBLE, MPI_MAX, PetscObjectComm((PetscObject)matin)));

    1472:     info->nz_used      = irecv[0];
    1473:     info->nz_allocated = irecv[1];
    1474:     info->nz_unneeded  = irecv[2];
    1475:     info->memory       = irecv[3];
    1476:     info->mallocs      = irecv[4];
    1477:   } else if (flag == MAT_GLOBAL_SUM) {
    1478:     PetscCallMPI(MPIU_Allreduce(isend, irecv, 5, MPIU_PETSCLOGDOUBLE, MPI_SUM, PetscObjectComm((PetscObject)matin)));

    1480:     info->nz_used      = irecv[0];
    1481:     info->nz_allocated = irecv[1];
    1482:     info->nz_unneeded  = irecv[2];
    1483:     info->memory       = irecv[3];
    1484:     info->mallocs      = irecv[4];
    1485:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown MatInfoType argument %d", (int)flag);
    1486:   info->fill_ratio_given  = 0; /* no parallel LU/ILU/Cholesky */
    1487:   info->fill_ratio_needed = 0;
    1488:   info->factor_mallocs    = 0;
    1489:   PetscFunctionReturn(PETSC_SUCCESS);
    1490: }

    1492: static PetscErrorCode MatSetOption_MPISBAIJ(Mat A, MatOption op, PetscBool flg)
    1493: {
    1494:   Mat_MPISBAIJ *a  = (Mat_MPISBAIJ *)A->data;
    1495:   Mat_SeqSBAIJ *aA = (Mat_SeqSBAIJ *)a->A->data;

    1497:   PetscFunctionBegin;
    1498:   switch (op) {
    1499:   case MAT_NEW_NONZERO_LOCATIONS:
    1500:   case MAT_NEW_NONZERO_ALLOCATION_ERR:
    1501:   case MAT_UNUSED_NONZERO_LOCATION_ERR:
    1502:   case MAT_KEEP_NONZERO_PATTERN:
    1503:   case MAT_NEW_NONZERO_LOCATION_ERR:
    1504:     MatCheckPreallocated(A, 1);
    1505:     PetscCall(MatSetOption(a->A, op, flg));
    1506:     PetscCall(MatSetOption(a->B, op, flg));
    1507:     break;
    1508:   case MAT_ROW_ORIENTED:
    1509:     MatCheckPreallocated(A, 1);
    1510:     a->roworiented = flg;

    1512:     PetscCall(MatSetOption(a->A, op, flg));
    1513:     PetscCall(MatSetOption(a->B, op, flg));
    1514:     break;
    1515:   case MAT_IGNORE_OFF_PROC_ENTRIES:
    1516:     a->donotstash = flg;
    1517:     break;
    1518:   case MAT_USE_HASH_TABLE:
    1519:     a->ht_flag = flg;
    1520:     break;
    1521:   case MAT_HERMITIAN:
    1522:     if (a->A && A->rmap->n == A->cmap->n) PetscCall(MatSetOption(a->A, op, flg));
    1523: #if defined(PETSC_USE_COMPLEX)
    1524:     if (flg) { /* need different mat-vec ops */
    1525:       A->ops->mult             = MatMult_MPISBAIJ_Hermitian;
    1526:       A->ops->multadd          = MatMultAdd_MPISBAIJ_Hermitian;
    1527:       A->ops->multtranspose    = NULL;
    1528:       A->ops->multtransposeadd = NULL;
    1529:     }
    1530: #endif
    1531:     break;
    1532:   case MAT_SPD:
    1533:   case MAT_SYMMETRIC:
    1534:     if (a->A && A->rmap->n == A->cmap->n) PetscCall(MatSetOption(a->A, op, flg));
    1535: #if defined(PETSC_USE_COMPLEX)
    1536:     if (flg) { /* restore to use default mat-vec ops */
    1537:       A->ops->mult             = MatMult_MPISBAIJ;
    1538:       A->ops->multadd          = MatMultAdd_MPISBAIJ;
    1539:       A->ops->multtranspose    = MatMult_MPISBAIJ;
    1540:       A->ops->multtransposeadd = MatMultAdd_MPISBAIJ;
    1541:     }
    1542: #endif
    1543:     break;
    1544:   case MAT_STRUCTURALLY_SYMMETRIC:
    1545:     if (a->A && A->rmap->n == A->cmap->n) PetscCall(MatSetOption(a->A, op, flg));
    1546:     break;
    1547:   case MAT_IGNORE_LOWER_TRIANGULAR:
    1548:   case MAT_ERROR_LOWER_TRIANGULAR:
    1549:     aA->ignore_ltriangular = flg;
    1550:     break;
    1551:   case MAT_GETROW_UPPERTRIANGULAR:
    1552:     aA->getrow_utriangular = flg;
    1553:     break;
    1554:   default:
    1555:     break;
    1556:   }
    1557:   PetscFunctionReturn(PETSC_SUCCESS);
    1558: }

    1560: static PetscErrorCode MatTranspose_MPISBAIJ(Mat A, MatReuse reuse, Mat *B)
    1561: {
    1562:   PetscFunctionBegin;
    1563:   if (reuse == MAT_REUSE_MATRIX) PetscCall(MatTransposeCheckNonzeroState_Private(A, *B));
    1564:   if (reuse == MAT_INITIAL_MATRIX) {
    1565:     PetscCall(MatDuplicate(A, MAT_COPY_VALUES, B));
    1566:   } else if (reuse == MAT_REUSE_MATRIX) {
    1567:     PetscCall(MatCopy(A, *B, SAME_NONZERO_PATTERN));
    1568:   }
    1569:   PetscFunctionReturn(PETSC_SUCCESS);
    1570: }

    1572: static PetscErrorCode MatDiagonalScale_MPISBAIJ(Mat mat, Vec ll, Vec rr)
    1573: {
    1574:   Mat_MPISBAIJ *baij = (Mat_MPISBAIJ *)mat->data;
    1575:   Mat           a = baij->A, b = baij->B;
    1576:   PetscInt      nv, m, n;

    1578:   PetscFunctionBegin;
    1579:   if (!ll) PetscFunctionReturn(PETSC_SUCCESS);

    1581:   PetscCall(MatGetLocalSize(mat, &m, &n));
    1582:   PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "For symmetric format, local size %" PetscInt_FMT " %" PetscInt_FMT " must be same", m, n);

    1584:   PetscCall(VecGetLocalSize(rr, &nv));
    1585:   PetscCheck(nv == n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Left and right vector non-conforming local size");

    1587:   PetscCall(VecScatterBegin(baij->Mvctx, rr, baij->lvec, INSERT_VALUES, SCATTER_FORWARD));

    1589:   /* left diagonalscale the off-diagonal part */
    1590:   PetscUseTypeMethod(b, diagonalscale, ll, NULL);

    1592:   /* scale the diagonal part */
    1593:   PetscUseTypeMethod(a, diagonalscale, ll, rr);

    1595:   /* right diagonalscale the off-diagonal part */
    1596:   PetscCall(VecScatterEnd(baij->Mvctx, rr, baij->lvec, INSERT_VALUES, SCATTER_FORWARD));
    1597:   PetscUseTypeMethod(b, diagonalscale, NULL, baij->lvec);
    1598:   PetscFunctionReturn(PETSC_SUCCESS);
    1599: }

    1601: static PetscErrorCode MatSetUnfactored_MPISBAIJ(Mat A)
    1602: {
    1603:   Mat_MPISBAIJ *a = (Mat_MPISBAIJ *)A->data;

    1605:   PetscFunctionBegin;
    1606:   PetscCall(MatSetUnfactored(a->A));
    1607:   PetscFunctionReturn(PETSC_SUCCESS);
    1608: }

    1610: static PetscErrorCode MatDuplicate_MPISBAIJ(Mat, MatDuplicateOption, Mat *);

    1612: static PetscErrorCode MatEqual_MPISBAIJ(Mat A, Mat B, PetscBool *flag)
    1613: {
    1614:   Mat_MPISBAIJ *matB = (Mat_MPISBAIJ *)B->data, *matA = (Mat_MPISBAIJ *)A->data;
    1615:   Mat           a, b, c, d;
    1616:   PetscBool     flg;

    1618:   PetscFunctionBegin;
    1619:   a = matA->A;
    1620:   b = matA->B;
    1621:   c = matB->A;
    1622:   d = matB->B;

    1624:   PetscCall(MatEqual(a, c, &flg));
    1625:   if (flg) PetscCall(MatEqual(b, d, &flg));
    1626:   PetscCallMPI(MPIU_Allreduce(&flg, flag, 1, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A)));
    1627:   PetscFunctionReturn(PETSC_SUCCESS);
    1628: }

    1630: static PetscErrorCode MatCopy_MPISBAIJ(Mat A, Mat B, MatStructure str)
    1631: {
    1632:   PetscBool isbaij;

    1634:   PetscFunctionBegin;
    1635:   PetscCall(PetscObjectTypeCompareAny((PetscObject)B, &isbaij, MATSEQSBAIJ, MATMPISBAIJ, ""));
    1636:   PetscCheck(isbaij, PetscObjectComm((PetscObject)B), PETSC_ERR_SUP, "Not for matrix type %s", ((PetscObject)B)->type_name);
    1637:   /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */
    1638:   if ((str != SAME_NONZERO_PATTERN) || (A->ops->copy != B->ops->copy)) {
    1639:     PetscCall(MatGetRowUpperTriangular(A));
    1640:     PetscCall(MatCopy_Basic(A, B, str));
    1641:     PetscCall(MatRestoreRowUpperTriangular(A));
    1642:   } else {
    1643:     Mat_MPISBAIJ *a = (Mat_MPISBAIJ *)A->data;
    1644:     Mat_MPISBAIJ *b = (Mat_MPISBAIJ *)B->data;

    1646:     PetscCall(MatCopy(a->A, b->A, str));
    1647:     PetscCall(MatCopy(a->B, b->B, str));
    1648:   }
    1649:   PetscCall(PetscObjectStateIncrease((PetscObject)B));
    1650:   PetscFunctionReturn(PETSC_SUCCESS);
    1651: }

    1653: static PetscErrorCode MatAXPY_MPISBAIJ(Mat Y, PetscScalar a, Mat X, MatStructure str)
    1654: {
    1655:   Mat_MPISBAIJ *xx = (Mat_MPISBAIJ *)X->data, *yy = (Mat_MPISBAIJ *)Y->data;
    1656:   PetscBLASInt  bnz, one                          = 1;
    1657:   Mat_SeqSBAIJ *xa, *ya;
    1658:   Mat_SeqBAIJ  *xb, *yb;

    1660:   PetscFunctionBegin;
    1661:   if (str == SAME_NONZERO_PATTERN) {
    1662:     PetscScalar alpha = a;
    1663:     xa                = (Mat_SeqSBAIJ *)xx->A->data;
    1664:     ya                = (Mat_SeqSBAIJ *)yy->A->data;
    1665:     PetscCall(PetscBLASIntCast(xa->nz, &bnz));
    1666:     PetscCallBLAS("BLASaxpy", BLASaxpy_(&bnz, &alpha, xa->a, &one, ya->a, &one));
    1667:     xb = (Mat_SeqBAIJ *)xx->B->data;
    1668:     yb = (Mat_SeqBAIJ *)yy->B->data;
    1669:     PetscCall(PetscBLASIntCast(xb->nz, &bnz));
    1670:     PetscCallBLAS("BLASaxpy", BLASaxpy_(&bnz, &alpha, xb->a, &one, yb->a, &one));
    1671:     PetscCall(PetscObjectStateIncrease((PetscObject)Y));
    1672:   } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */
    1673:     PetscCall(MatSetOption(X, MAT_GETROW_UPPERTRIANGULAR, PETSC_TRUE));
    1674:     PetscCall(MatAXPY_Basic(Y, a, X, str));
    1675:     PetscCall(MatSetOption(X, MAT_GETROW_UPPERTRIANGULAR, PETSC_FALSE));
    1676:   } else {
    1677:     Mat       B;
    1678:     PetscInt *nnz_d, *nnz_o, bs = Y->rmap->bs;
    1679:     PetscCheck(bs == X->rmap->bs, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Matrices must have same block size");
    1680:     PetscCall(MatGetRowUpperTriangular(X));
    1681:     PetscCall(MatGetRowUpperTriangular(Y));
    1682:     PetscCall(PetscMalloc1(yy->A->rmap->N, &nnz_d));
    1683:     PetscCall(PetscMalloc1(yy->B->rmap->N, &nnz_o));
    1684:     PetscCall(MatCreate(PetscObjectComm((PetscObject)Y), &B));
    1685:     PetscCall(PetscObjectSetName((PetscObject)B, ((PetscObject)Y)->name));
    1686:     PetscCall(MatSetSizes(B, Y->rmap->n, Y->cmap->n, Y->rmap->N, Y->cmap->N));
    1687:     PetscCall(MatSetBlockSizesFromMats(B, Y, Y));
    1688:     PetscCall(MatSetType(B, MATMPISBAIJ));
    1689:     PetscCall(MatAXPYGetPreallocation_SeqSBAIJ(yy->A, xx->A, nnz_d));
    1690:     PetscCall(MatAXPYGetPreallocation_MPIBAIJ(yy->B, yy->garray, xx->B, xx->garray, nnz_o));
    1691:     PetscCall(MatMPISBAIJSetPreallocation(B, bs, 0, nnz_d, 0, nnz_o));
    1692:     PetscCall(MatAXPY_BasicWithPreallocation(B, Y, a, X, str));
    1693:     PetscCall(MatHeaderMerge(Y, &B));
    1694:     PetscCall(PetscFree(nnz_d));
    1695:     PetscCall(PetscFree(nnz_o));
    1696:     PetscCall(MatRestoreRowUpperTriangular(X));
    1697:     PetscCall(MatRestoreRowUpperTriangular(Y));
    1698:   }
    1699:   PetscFunctionReturn(PETSC_SUCCESS);
    1700: }

    1702: static PetscErrorCode MatCreateSubMatrices_MPISBAIJ(Mat A, PetscInt n, const IS irow[], const IS icol[], MatReuse scall, Mat *B[])
    1703: {
    1704:   PetscInt  i;
    1705:   PetscBool flg;

    1707:   PetscFunctionBegin;
    1708:   PetscCall(MatCreateSubMatrices_MPIBAIJ(A, n, irow, icol, scall, B)); /* B[] are sbaij matrices */
    1709:   for (i = 0; i < n; i++) {
    1710:     PetscCall(ISEqual(irow[i], icol[i], &flg));
    1711:     if (!flg) PetscCall(MatSeqSBAIJZeroOps_Private(*B[i]));
    1712:   }
    1713:   PetscFunctionReturn(PETSC_SUCCESS);
    1714: }

    1716: static PetscErrorCode MatShift_MPISBAIJ(Mat Y, PetscScalar a)
    1717: {
    1718:   Mat_MPISBAIJ *maij = (Mat_MPISBAIJ *)Y->data;
    1719:   Mat_SeqSBAIJ *aij  = (Mat_SeqSBAIJ *)maij->A->data;

    1721:   PetscFunctionBegin;
    1722:   if (!Y->preallocated) {
    1723:     PetscCall(MatMPISBAIJSetPreallocation(Y, Y->rmap->bs, 1, NULL, 0, NULL));
    1724:   } else if (!aij->nz) {
    1725:     PetscInt nonew = aij->nonew;
    1726:     PetscCall(MatSeqSBAIJSetPreallocation(maij->A, Y->rmap->bs, 1, NULL));
    1727:     aij->nonew = nonew;
    1728:   }
    1729:   PetscCall(MatShift_Basic(Y, a));
    1730:   PetscFunctionReturn(PETSC_SUCCESS);
    1731: }

    1733: static PetscErrorCode MatZeroRowsColumns_MPISBAIJ(Mat A, PetscInt N, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
    1734: {
    1735:   Mat_MPISBAIJ      *l = (Mat_MPISBAIJ *)A->data;
    1736:   PetscMPIInt        n, p = 0;
    1737:   PetscInt           i, j, k, r, len = 0, row, col, count;
    1738:   PetscInt          *lrows, *owners = A->rmap->range;
    1739:   PetscSFNode       *rrows;
    1740:   PetscSF            sf;
    1741:   const PetscScalar *xx;
    1742:   PetscScalar       *bb, *mask;
    1743:   Vec                xmask, lmask, lvec_contrib = NULL;
    1744:   Mat_SeqBAIJ       *baij = (Mat_SeqBAIJ *)l->B->data;
    1745:   PetscInt           bs = A->rmap->bs, bs2 = baij->bs2;
    1746:   PetscScalar       *aa;

    1748:   PetscFunctionBegin;
    1749:   PetscCall(PetscMPIIntCast(A->rmap->n, &n));
    1750:   /* create PetscSF where leaves are input rows and roots are owned rows */
    1751:   PetscCall(PetscMalloc1(n, &lrows));
    1752:   for (r = 0; r < n; ++r) lrows[r] = -1;
    1753:   PetscCall(PetscMalloc1(N, &rrows));
    1754:   for (r = 0; r < N; ++r) {
    1755:     const PetscInt idx = rows[r];
    1756:     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);
    1757:     if (idx < owners[p] || owners[p + 1] <= idx) { /* short-circuit the search if the last p owns this row too */
    1758:       PetscCall(PetscLayoutFindOwner(A->rmap, idx, &p));
    1759:     }
    1760:     rrows[r].rank  = p;
    1761:     rrows[r].index = rows[r] - owners[p];
    1762:   }
    1763:   PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)A), &sf));
    1764:   PetscCall(PetscSFSetGraph(sf, n, N, NULL, PETSC_OWN_POINTER, rrows, PETSC_OWN_POINTER));
    1765:   /* collect flags for rows to be zeroed */
    1766:   PetscCall(PetscSFReduceBegin(sf, MPIU_INT, (PetscInt *)rows, lrows, MPI_LOR));
    1767:   PetscCall(PetscSFReduceEnd(sf, MPIU_INT, (PetscInt *)rows, lrows, MPI_LOR));
    1768:   PetscCall(PetscSFDestroy(&sf));
    1769:   /* compress and put in row numbers */
    1770:   for (r = 0; r < n; ++r) {
    1771:     if (lrows[r] >= 0) lrows[len++] = r;
    1772:   }
    1773:   /* zero diagonal part of matrix */
    1774:   PetscCall(MatZeroRowsColumns(l->A, len, lrows, diag, x, b));
    1775:   /* handle off-diagonal part of matrix */
    1776:   PetscCall(MatCreateVecs(A, &xmask, NULL));
    1777:   PetscCall(VecDuplicate(l->lvec, &lmask));
    1778:   PetscCall(VecGetArray(xmask, &bb));
    1779:   for (i = 0; i < len; i++) bb[lrows[i]] = 1;
    1780:   PetscCall(VecRestoreArray(xmask, &bb));
    1781:   PetscCall(VecScatterBegin(l->Mvctx, xmask, lmask, ADD_VALUES, SCATTER_FORWARD));
    1782:   PetscCall(VecScatterEnd(l->Mvctx, xmask, lmask, ADD_VALUES, SCATTER_FORWARD));
    1783:   PetscCall(VecDestroy(&xmask));
    1784:   if (x) {
    1785:     PetscCall(VecScatterBegin(l->Mvctx, x, l->lvec, INSERT_VALUES, SCATTER_FORWARD));
    1786:     PetscCall(VecScatterEnd(l->Mvctx, x, l->lvec, INSERT_VALUES, SCATTER_FORWARD));
    1787:     PetscCall(VecGetArrayRead(l->lvec, &xx));
    1788:     PetscCall(VecGetArray(b, &bb));
    1789:   }
    1790:   PetscCall(VecGetArray(lmask, &mask));
    1791:   /* MPISBAIJ stores only the upper off-diagonal in l->B; for each zeroed local row r and
    1792:      non-zeroed off-process column c in that row, accumulate -A[r,c] * x[r] into lvec_contrib.
    1793:      A SCATTER_REVERSE below sends these contributions to b[c] on the owning (higher-rank)
    1794:      process, the missing symmetric lower-triangular update. We skip entries where c is
    1795:      also a zeroed row (mask[col] != 0) since b[c] = diag * x[c] is handled separately. */
    1796:   if (x) {
    1797:     const PetscScalar *x_vals;
    1798:     PetscScalar       *c_vals;

    1800:     PetscCall(VecDuplicate(l->lvec, &lvec_contrib));
    1801:     PetscCall(VecGetArray(lvec_contrib, &c_vals));
    1802:     PetscCall(VecGetArrayRead(x, &x_vals));
    1803:     /* Only accumulate b[c] -= A[r,c] * x[r] when off-process col c is not also a zeroed row
    1804:        (mask[c] non-zero means col c is zeroed, so b[c] = diag * x[c] is already set).
    1805:        This mirrors the MatSeqSBAIJ pattern: if (zeroed[r] && !zeroed[c]) bb[c] -= A[r,c] * x[r].
    1806:        c_vals is indexed by the local B column index. */
    1807:     for (i = 0; i < len; ++i) {
    1808:       row = lrows[i];
    1809:       for (j = baij->i[row / bs]; j < baij->i[row / bs + 1]; ++j) {
    1810:         for (k = 0; k < bs; ++k) {
    1811:           col = baij->j[j] * bs + k;
    1812:           if (!PetscAbsScalar(mask[col])) {
    1813:             aa = baij->a + j * bs2 + (row % bs) + bs * k;
    1814:             c_vals[col] -= aa[0] * x_vals[row];
    1815:           }
    1816:         }
    1817:       }
    1818:     }
    1819:     PetscCall(VecRestoreArrayRead(x, &x_vals));
    1820:     PetscCall(VecRestoreArray(lvec_contrib, &c_vals));
    1821:   }
    1822:   /* remove zeroed rows of off-diagonal matrix */
    1823:   for (i = 0; i < len; ++i) {
    1824:     row   = lrows[i];
    1825:     count = (baij->i[row / bs + 1] - baij->i[row / bs]) * bs;
    1826:     aa    = baij->a + baij->i[row / bs] * bs2 + (row % bs);
    1827:     for (k = 0; k < count; ++k) {
    1828:       aa[0] = 0.0;
    1829:       aa += bs;
    1830:     }
    1831:   }
    1832:   /* loop over all elements of off process part of matrix zeroing removed columns */
    1833:   for (i = 0; i < l->B->rmap->N; ++i) {
    1834:     row = i / bs;
    1835:     for (j = baij->i[row]; j < baij->i[row + 1]; ++j) {
    1836:       for (k = 0; k < bs; ++k) {
    1837:         col = bs * baij->j[j] + k;
    1838:         if (PetscAbsScalar(mask[col])) {
    1839:           aa = baij->a + j * bs2 + (i % bs) + bs * k;
    1840:           if (x) bb[i] -= aa[0] * xx[col];
    1841:           aa[0] = 0.0;
    1842:         }
    1843:       }
    1844:     }
    1845:   }
    1846:   if (x) {
    1847:     PetscCall(VecRestoreArray(b, &bb));
    1848:     PetscCall(VecRestoreArrayRead(l->lvec, &xx));
    1849:     /* scatter the accumulated contributions to b[c] on higher-rank processes owning column c */
    1850:     PetscCall(VecScatterBegin(l->Mvctx, lvec_contrib, b, ADD_VALUES, SCATTER_REVERSE));
    1851:     PetscCall(VecScatterEnd(l->Mvctx, lvec_contrib, b, ADD_VALUES, SCATTER_REVERSE));
    1852:     PetscCall(VecDestroy(&lvec_contrib));
    1853:   }
    1854:   PetscCall(VecRestoreArray(lmask, &mask));
    1855:   PetscCall(VecDestroy(&lmask));
    1856:   PetscCall(PetscFree(lrows));

    1858:   /* only change matrix nonzero state if pattern was allowed to be changed */
    1859:   if (!((Mat_SeqSBAIJ *)l->A->data)->nonew) {
    1860:     PetscObjectState state = l->A->nonzerostate + l->B->nonzerostate;
    1861:     PetscCallMPI(MPIU_Allreduce(&state, &A->nonzerostate, 1, MPIU_INT64, MPI_SUM, PetscObjectComm((PetscObject)A)));
    1862:   }
    1863:   PetscFunctionReturn(PETSC_SUCCESS);
    1864: }

    1866: static PetscErrorCode MatGetDiagonalBlock_MPISBAIJ(Mat A, Mat *a)
    1867: {
    1868:   PetscFunctionBegin;
    1869:   *a = ((Mat_MPISBAIJ *)A->data)->A;
    1870:   PetscFunctionReturn(PETSC_SUCCESS);
    1871: }

    1873: static PetscErrorCode MatEliminateZeros_MPISBAIJ(Mat A, PetscBool keep)
    1874: {
    1875:   Mat_MPISBAIJ *a = (Mat_MPISBAIJ *)A->data;

    1877:   PetscFunctionBegin;
    1878:   PetscCall(MatEliminateZeros_SeqSBAIJ(a->A, keep));       // possibly keep zero diagonal coefficients
    1879:   PetscCall(MatEliminateZeros_SeqBAIJ(a->B, PETSC_FALSE)); // never keep zero diagonal coefficients
    1880:   PetscFunctionReturn(PETSC_SUCCESS);
    1881: }

    1883: static PetscErrorCode MatLoad_MPISBAIJ(Mat, PetscViewer);
    1884: static PetscErrorCode MatGetRowMaxAbs_MPISBAIJ(Mat, Vec, PetscInt[]);
    1885: static PetscErrorCode MatSOR_MPISBAIJ(Mat, Vec, PetscReal, MatSORType, PetscReal, PetscInt, PetscInt, Vec);

    1887: static struct _MatOps MatOps_Values = {MatSetValues_MPISBAIJ,
    1888:                                        MatGetRow_MPISBAIJ,
    1889:                                        MatRestoreRow_MPISBAIJ,
    1890:                                        MatMult_MPISBAIJ,
    1891:                                        /*  4*/ MatMultAdd_MPISBAIJ,
    1892:                                        MatMult_MPISBAIJ, /* transpose versions are same as non-transpose */
    1893:                                        MatMultAdd_MPISBAIJ,
    1894:                                        NULL,
    1895:                                        NULL,
    1896:                                        NULL,
    1897:                                        /* 10*/ NULL,
    1898:                                        NULL,
    1899:                                        NULL,
    1900:                                        MatSOR_MPISBAIJ,
    1901:                                        MatTranspose_MPISBAIJ,
    1902:                                        /* 15*/ MatGetInfo_MPISBAIJ,
    1903:                                        MatEqual_MPISBAIJ,
    1904:                                        MatGetDiagonal_MPISBAIJ,
    1905:                                        MatDiagonalScale_MPISBAIJ,
    1906:                                        MatNorm_MPISBAIJ,
    1907:                                        /* 20*/ MatAssemblyBegin_MPISBAIJ,
    1908:                                        MatAssemblyEnd_MPISBAIJ,
    1909:                                        MatSetOption_MPISBAIJ,
    1910:                                        MatZeroEntries_MPISBAIJ,
    1911:                                        /* 24*/ NULL,
    1912:                                        NULL,
    1913:                                        NULL,
    1914:                                        NULL,
    1915:                                        NULL,
    1916:                                        /* 29*/ MatSetUp_MPI_Hash,
    1917:                                        NULL,
    1918:                                        NULL,
    1919:                                        MatGetDiagonalBlock_MPISBAIJ,
    1920:                                        NULL,
    1921:                                        /* 34*/ MatDuplicate_MPISBAIJ,
    1922:                                        NULL,
    1923:                                        NULL,
    1924:                                        NULL,
    1925:                                        NULL,
    1926:                                        /* 39*/ MatAXPY_MPISBAIJ,
    1927:                                        MatCreateSubMatrices_MPISBAIJ,
    1928:                                        MatIncreaseOverlap_MPISBAIJ,
    1929:                                        MatGetValues_MPISBAIJ,
    1930:                                        MatCopy_MPISBAIJ,
    1931:                                        /* 44*/ NULL,
    1932:                                        MatScale_MPISBAIJ,
    1933:                                        MatShift_MPISBAIJ,
    1934:                                        NULL,
    1935:                                        MatZeroRowsColumns_MPISBAIJ,
    1936:                                        /* 49*/ NULL,
    1937:                                        NULL,
    1938:                                        NULL,
    1939:                                        NULL,
    1940:                                        NULL,
    1941:                                        /* 54*/ NULL,
    1942:                                        NULL,
    1943:                                        MatSetUnfactored_MPISBAIJ,
    1944:                                        NULL,
    1945:                                        MatSetValuesBlocked_MPISBAIJ,
    1946:                                        /* 59*/ MatCreateSubMatrix_MPISBAIJ,
    1947:                                        NULL,
    1948:                                        NULL,
    1949:                                        NULL,
    1950:                                        NULL,
    1951:                                        /* 64*/ NULL,
    1952:                                        NULL,
    1953:                                        NULL,
    1954:                                        NULL,
    1955:                                        MatGetRowMaxAbs_MPISBAIJ,
    1956:                                        /* 69*/ NULL,
    1957:                                        MatConvert_MPISBAIJ_Basic,
    1958:                                        NULL,
    1959:                                        NULL,
    1960:                                        NULL,
    1961:                                        NULL,
    1962:                                        NULL,
    1963:                                        NULL,
    1964:                                        NULL,
    1965:                                        MatLoad_MPISBAIJ,
    1966:                                        /* 79*/ NULL,
    1967:                                        NULL,
    1968:                                        NULL,
    1969:                                        NULL,
    1970:                                        NULL,
    1971:                                        /* 84*/ NULL,
    1972:                                        NULL,
    1973:                                        NULL,
    1974:                                        NULL,
    1975:                                        NULL,
    1976:                                        /* 89*/ NULL,
    1977:                                        NULL,
    1978:                                        NULL,
    1979:                                        NULL,
    1980:                                        MatConjugate_MPISBAIJ,
    1981:                                        /* 94*/ NULL,
    1982:                                        NULL,
    1983:                                        MatRealPart_MPISBAIJ,
    1984:                                        MatImaginaryPart_MPISBAIJ,
    1985:                                        MatGetRowUpperTriangular_MPISBAIJ,
    1986:                                        /* 99*/ MatRestoreRowUpperTriangular_MPISBAIJ,
    1987:                                        NULL,
    1988:                                        NULL,
    1989:                                        NULL,
    1990:                                        NULL,
    1991:                                        /*104*/ NULL,
    1992:                                        NULL,
    1993:                                        NULL,
    1994:                                        NULL,
    1995:                                        NULL,
    1996:                                        /*109*/ NULL,
    1997:                                        NULL,
    1998:                                        NULL,
    1999:                                        NULL,
    2000:                                        NULL,
    2001:                                        /*114*/ NULL,
    2002:                                        NULL,
    2003:                                        NULL,
    2004:                                        NULL,
    2005:                                        NULL,
    2006:                                        /*119*/ NULL,
    2007:                                        NULL,
    2008:                                        NULL,
    2009:                                        NULL,
    2010:                                        NULL,
    2011:                                        /*124*/ NULL,
    2012:                                        NULL,
    2013:                                        MatSetBlockSizes_Default,
    2014:                                        NULL,
    2015:                                        NULL,
    2016:                                        /*129*/ NULL,
    2017:                                        MatCreateMPIMatConcatenateSeqMat_MPISBAIJ,
    2018:                                        NULL,
    2019:                                        NULL,
    2020:                                        NULL,
    2021:                                        /*134*/ NULL,
    2022:                                        NULL,
    2023:                                        MatEliminateZeros_MPISBAIJ,
    2024:                                        NULL,
    2025:                                        NULL,
    2026:                                        /*139*/ NULL,
    2027:                                        NULL,
    2028:                                        MatCopyHashToXAIJ_MPI_Hash,
    2029:                                        NULL,
    2030:                                        NULL,
    2031:                                        /*144*/ MatADot_Default,
    2032:                                        MatANorm_Default,
    2033:                                        NULL,
    2034:                                        NULL};

    2036: static PetscErrorCode MatMPISBAIJSetPreallocation_MPISBAIJ(Mat B, PetscInt bs, PetscInt d_nz, const PetscInt *d_nnz, PetscInt o_nz, const PetscInt *o_nnz)
    2037: {
    2038:   Mat_MPISBAIJ *b = (Mat_MPISBAIJ *)B->data;
    2039:   PetscInt      i, mbs, Mbs;
    2040:   PetscMPIInt   size;

    2042:   PetscFunctionBegin;
    2043:   if (B->hash_active) {
    2044:     B->ops[0]      = b->cops;
    2045:     B->hash_active = PETSC_FALSE;
    2046:   }
    2047:   if (!B->preallocated) PetscCall(MatStashCreate_Private(PetscObjectComm((PetscObject)B), bs, &B->bstash));
    2048:   PetscCall(MatSetBlockSize(B, bs));
    2049:   PetscCall(PetscLayoutSetUp(B->rmap));
    2050:   PetscCall(PetscLayoutSetUp(B->cmap));
    2051:   PetscCall(PetscLayoutGetBlockSize(B->rmap, &bs));
    2052:   PetscCheck(B->rmap->N <= B->cmap->N, PetscObjectComm((PetscObject)B), PETSC_ERR_SUP, "MPISBAIJ matrix cannot have more rows %" PetscInt_FMT " than columns %" PetscInt_FMT, B->rmap->N, B->cmap->N);
    2053:   PetscCheck(B->rmap->n <= B->cmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "MPISBAIJ matrix cannot have more local rows %" PetscInt_FMT " than columns %" PetscInt_FMT, B->rmap->n, B->cmap->n);

    2055:   mbs = B->rmap->n / bs;
    2056:   Mbs = B->rmap->N / bs;
    2057:   PetscCheck(mbs * bs == B->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "No of local rows %" PetscInt_FMT " must be divisible by blocksize %" PetscInt_FMT, B->rmap->N, bs);

    2059:   B->rmap->bs = bs;
    2060:   b->bs2      = bs * bs;
    2061:   b->mbs      = mbs;
    2062:   b->Mbs      = Mbs;
    2063:   b->nbs      = B->cmap->n / bs;
    2064:   b->Nbs      = B->cmap->N / bs;

    2066:   for (i = 0; i <= b->size; i++) b->rangebs[i] = B->rmap->range[i] / bs;
    2067:   b->rstartbs = B->rmap->rstart / bs;
    2068:   b->rendbs   = B->rmap->rend / bs;

    2070:   b->cstartbs = B->cmap->rstart / bs;
    2071:   b->cendbs   = B->cmap->rend / bs;

    2073: #if defined(PETSC_USE_CTABLE)
    2074:   PetscCall(PetscHMapIDestroy(&b->colmap));
    2075: #else
    2076:   PetscCall(PetscFree(b->colmap));
    2077: #endif
    2078:   PetscCall(PetscFree(b->garray));
    2079:   PetscCall(VecDestroy(&b->lvec));
    2080:   PetscCall(VecScatterDestroy(&b->Mvctx));
    2081:   PetscCall(VecDestroy(&b->slvec0));
    2082:   PetscCall(VecDestroy(&b->slvec0b));
    2083:   PetscCall(VecDestroy(&b->slvec1));
    2084:   PetscCall(VecDestroy(&b->slvec1a));
    2085:   PetscCall(VecDestroy(&b->slvec1b));
    2086:   PetscCall(VecScatterDestroy(&b->sMvctx));

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

    2090:   MatSeqXAIJGetOptions_Private(b->B);
    2091:   PetscCall(MatDestroy(&b->B));
    2092:   PetscCall(MatCreate(PETSC_COMM_SELF, &b->B));
    2093:   PetscCall(MatSetSizes(b->B, B->rmap->n, size > 1 ? B->cmap->N : 0, B->rmap->n, size > 1 ? B->cmap->N : 0));
    2094:   PetscCall(MatSetType(b->B, MATSEQBAIJ));
    2095:   MatSeqXAIJRestoreOptions_Private(b->B);

    2097:   MatSeqXAIJGetOptions_Private(b->A);
    2098:   PetscCall(MatDestroy(&b->A));
    2099:   PetscCall(MatCreate(PETSC_COMM_SELF, &b->A));
    2100:   PetscCall(MatSetSizes(b->A, B->rmap->n, B->cmap->n, B->rmap->n, B->cmap->n));
    2101:   PetscCall(MatSetType(b->A, MATSEQSBAIJ));
    2102:   MatSeqXAIJRestoreOptions_Private(b->A);

    2104:   PetscCall(MatSeqSBAIJSetPreallocation(b->A, bs, d_nz, d_nnz));
    2105:   PetscCall(MatSeqBAIJSetPreallocation(b->B, bs, o_nz, o_nnz));

    2107:   B->preallocated  = PETSC_TRUE;
    2108:   B->was_assembled = PETSC_FALSE;
    2109:   B->assembled     = PETSC_FALSE;
    2110:   PetscFunctionReturn(PETSC_SUCCESS);
    2111: }

    2113: static PetscErrorCode MatMPISBAIJSetPreallocationCSR_MPISBAIJ(Mat B, PetscInt bs, const PetscInt ii[], const PetscInt jj[], const PetscScalar V[])
    2114: {
    2115:   PetscInt        m, rstart, cend;
    2116:   PetscInt        i, j, d, nz, bd, nz_max = 0, *d_nnz = NULL, *o_nnz = NULL;
    2117:   const PetscInt *JJ          = NULL;
    2118:   PetscScalar    *values      = NULL;
    2119:   PetscBool       roworiented = ((Mat_MPISBAIJ *)B->data)->roworiented;
    2120:   PetscBool       nooffprocentries;

    2122:   PetscFunctionBegin;
    2123:   PetscCheck(bs >= 1, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "Invalid block size specified, must be positive but it is %" PetscInt_FMT, bs);
    2124:   PetscCall(PetscLayoutSetBlockSize(B->rmap, bs));
    2125:   PetscCall(PetscLayoutSetBlockSize(B->cmap, bs));
    2126:   PetscCall(PetscLayoutSetUp(B->rmap));
    2127:   PetscCall(PetscLayoutSetUp(B->cmap));
    2128:   PetscCall(PetscLayoutGetBlockSize(B->rmap, &bs));
    2129:   m      = B->rmap->n / bs;
    2130:   rstart = B->rmap->rstart / bs;
    2131:   cend   = B->cmap->rend / bs;

    2133:   PetscCheck(!ii[0], PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "ii[0] must be 0 but it is %" PetscInt_FMT, ii[0]);
    2134:   PetscCall(PetscMalloc2(m, &d_nnz, m, &o_nnz));
    2135:   for (i = 0; i < m; i++) {
    2136:     nz = ii[i + 1] - ii[i];
    2137:     PetscCheck(nz >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local row %" PetscInt_FMT " has a negative number of columns %" PetscInt_FMT, i, nz);
    2138:     /* count the ones on the diagonal and above, split into diagonal and off-diagonal portions. */
    2139:     JJ = jj + ii[i];
    2140:     bd = 0;
    2141:     for (j = 0; j < nz; j++) {
    2142:       if (*JJ >= i + rstart) break;
    2143:       JJ++;
    2144:       bd++;
    2145:     }
    2146:     d = 0;
    2147:     for (; j < nz; j++) {
    2148:       if (*JJ++ >= cend) break;
    2149:       d++;
    2150:     }
    2151:     d_nnz[i] = d;
    2152:     o_nnz[i] = nz - d - bd;
    2153:     nz       = nz - bd;
    2154:     nz_max   = PetscMax(nz_max, nz);
    2155:   }
    2156:   PetscCall(MatMPISBAIJSetPreallocation(B, bs, 0, d_nnz, 0, o_nnz));
    2157:   PetscCall(MatSetOption(B, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE));
    2158:   PetscCall(PetscFree2(d_nnz, o_nnz));

    2160:   values = (PetscScalar *)V;
    2161:   if (!values) PetscCall(PetscCalloc1(bs * bs * nz_max, &values));
    2162:   for (i = 0; i < m; i++) {
    2163:     PetscInt        row   = i + rstart;
    2164:     PetscInt        ncols = ii[i + 1] - ii[i];
    2165:     const PetscInt *icols = jj + ii[i];
    2166:     if (bs == 1 || !roworiented) { /* block ordering matches the non-nested layout of MatSetValues so we can insert entire rows */
    2167:       const PetscScalar *svals = values + (V ? (bs * bs * ii[i]) : 0);
    2168:       PetscCall(MatSetValuesBlocked_MPISBAIJ(B, 1, &row, ncols, icols, svals, INSERT_VALUES));
    2169:     } else { /* block ordering does not match so we can only insert one block at a time. */
    2170:       PetscInt j;
    2171:       for (j = 0; j < ncols; j++) {
    2172:         const PetscScalar *svals = values + (V ? (bs * bs * (ii[i] + j)) : 0);
    2173:         PetscCall(MatSetValuesBlocked_MPISBAIJ(B, 1, &row, 1, &icols[j], svals, INSERT_VALUES));
    2174:       }
    2175:     }
    2176:   }

    2178:   if (!V) PetscCall(PetscFree(values));
    2179:   nooffprocentries    = B->nooffprocentries;
    2180:   B->nooffprocentries = PETSC_TRUE;
    2181:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
    2182:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
    2183:   B->nooffprocentries = nooffprocentries;

    2185:   PetscCall(MatSetOption(B, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
    2186:   PetscFunctionReturn(PETSC_SUCCESS);
    2187: }

    2189: /*MC
    2190:    MATMPISBAIJ - MATMPISBAIJ = "mpisbaij" - A matrix type to be used for distributed symmetric sparse block matrices,
    2191:    based on block compressed sparse row format.  Only the upper triangular portion of the "diagonal" portion of
    2192:    the matrix is stored.

    2194:    For complex numbers by default this matrix is symmetric, NOT Hermitian symmetric. To make it Hermitian symmetric you
    2195:    can call `MatSetOption`(`Mat`, `MAT_HERMITIAN`);

    2197:    Options Database Key:
    2198: . -mat_type mpisbaij - sets the matrix type to "mpisbaij" during a call to `MatSetFromOptions()`

    2200:    Level: beginner

    2202:    Note:
    2203:      The number of rows in the matrix must be less than or equal to the number of columns. Similarly the number of rows in the
    2204:      diagonal portion of the matrix of each process has to less than or equal the number of columns.

    2206: .seealso: [](ch_matrices), `Mat`, `MATSBAIJ`, `MATBAIJ`, `MatCreateBAIJ()`, `MATSEQSBAIJ`, `MatType`
    2207: M*/

    2209: PETSC_EXTERN PetscErrorCode MatCreate_MPISBAIJ(Mat B)
    2210: {
    2211:   Mat_MPISBAIJ *b;
    2212:   PetscBool     flg = PETSC_FALSE;

    2214:   PetscFunctionBegin;
    2215:   PetscCall(PetscNew(&b));
    2216:   B->data   = (void *)b;
    2217:   B->ops[0] = MatOps_Values;

    2219:   B->ops->destroy = MatDestroy_MPISBAIJ;
    2220:   B->ops->view    = MatView_MPISBAIJ;
    2221:   B->assembled    = PETSC_FALSE;
    2222:   B->insertmode   = NOT_SET_VALUES;

    2224:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)B), &b->rank));
    2225:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)B), &b->size));

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

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

    2233:   b->donotstash  = PETSC_FALSE;
    2234:   b->colmap      = NULL;
    2235:   b->garray      = NULL;
    2236:   b->roworiented = PETSC_TRUE;

    2238:   /* stuff used in block assembly */
    2239:   b->barray = NULL;

    2241:   /* stuff used for matrix vector multiply */
    2242:   b->lvec    = NULL;
    2243:   b->Mvctx   = NULL;
    2244:   b->slvec0  = NULL;
    2245:   b->slvec0b = NULL;
    2246:   b->slvec1  = NULL;
    2247:   b->slvec1a = NULL;
    2248:   b->slvec1b = NULL;
    2249:   b->sMvctx  = NULL;

    2251:   /* stuff for MatGetRow() */
    2252:   b->rowindices   = NULL;
    2253:   b->rowvalues    = NULL;
    2254:   b->getrowactive = PETSC_FALSE;

    2256:   /* hash table stuff */
    2257:   b->ht           = NULL;
    2258:   b->hd           = NULL;
    2259:   b->ht_size      = 0;
    2260:   b->ht_flag      = PETSC_FALSE;
    2261:   b->ht_fact      = 0;
    2262:   b->ht_total_ct  = 0;
    2263:   b->ht_insert_ct = 0;

    2265:   /* stuff for MatCreateSubMatrices_MPIBAIJ_local() */
    2266:   b->ijonly = PETSC_FALSE;

    2268:   b->in_loc = NULL;
    2269:   b->v_loc  = NULL;
    2270:   b->n_loc  = 0;

    2272:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatStoreValues_C", MatStoreValues_MPISBAIJ));
    2273:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatRetrieveValues_C", MatRetrieveValues_MPISBAIJ));
    2274:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatMPISBAIJSetPreallocation_C", MatMPISBAIJSetPreallocation_MPISBAIJ));
    2275:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatMPISBAIJSetPreallocationCSR_C", MatMPISBAIJSetPreallocationCSR_MPISBAIJ));
    2276: #if defined(PETSC_HAVE_ELEMENTAL)
    2277:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_mpisbaij_elemental_C", MatConvert_MPISBAIJ_Elemental));
    2278: #endif
    2279: #if defined(PETSC_HAVE_SCALAPACK) && (defined(PETSC_USE_REAL_SINGLE) || defined(PETSC_USE_REAL_DOUBLE))
    2280:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_mpisbaij_scalapack_C", MatConvert_SBAIJ_ScaLAPACK));
    2281: #endif
    2282:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_mpisbaij_mpiaij_C", MatConvert_MPISBAIJ_Basic));
    2283:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_mpisbaij_mpibaij_C", MatConvert_MPISBAIJ_Basic));

    2285:   B->symmetric                   = PETSC_BOOL3_TRUE;
    2286:   B->structurally_symmetric      = PETSC_BOOL3_TRUE;
    2287:   B->symmetry_eternal            = PETSC_TRUE;
    2288:   B->structural_symmetry_eternal = PETSC_TRUE;
    2289: #if !defined(PETSC_USE_COMPLEX)
    2290:   B->hermitian = PETSC_BOOL3_TRUE;
    2291: #endif

    2293:   PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATMPISBAIJ));
    2294:   PetscOptionsBegin(PetscObjectComm((PetscObject)B), NULL, "Options for loading MPISBAIJ matrix 1", "Mat");
    2295:   PetscCall(PetscOptionsBool("-mat_use_hash_table", "Use hash table to save memory in constructing matrix", "MatSetOption", flg, &flg, NULL));
    2296:   if (flg) {
    2297:     PetscReal fact = 1.39;
    2298:     PetscCall(MatSetOption(B, MAT_USE_HASH_TABLE, PETSC_TRUE));
    2299:     PetscCall(PetscOptionsReal("-mat_use_hash_table", "Use hash table factor", "MatMPIBAIJSetHashTableFactor", fact, &fact, NULL));
    2300:     if (fact <= 1.0) fact = 1.39;
    2301:     PetscCall(MatMPIBAIJSetHashTableFactor(B, fact));
    2302:     PetscCall(PetscInfo(B, "Hash table Factor used %5.2g\n", (double)fact));
    2303:   }
    2304:   PetscOptionsEnd();
    2305:   PetscFunctionReturn(PETSC_SUCCESS);
    2306: }

    2308: // PetscClangLinter pragma disable: -fdoc-section-header-unknown
    2309: /*MC
    2310:    MATSBAIJ - MATSBAIJ = "sbaij" - A matrix type to be used for symmetric block sparse matrices.

    2312:    This matrix type is identical to `MATSEQSBAIJ` when constructed with a single process communicator,
    2313:    and `MATMPISBAIJ` otherwise.

    2315:    Options Database Key:
    2316: . -mat_type sbaij - sets the matrix type to `MATSBAIJ` during a call to `MatSetFromOptions()`

    2318:   Level: beginner

    2320: .seealso: [](ch_matrices), `Mat`, `MATSEQSBAIJ`, `MATMPISBAIJ`, `MatCreateSBAIJ()`, `MATSEQBAIJ`, `MATMPIBAIJ`
    2321: M*/

    2323: /*@
    2324:   MatMPISBAIJSetPreallocation - For good matrix assembly performance
    2325:   the user should preallocate the matrix storage by setting the parameters
    2326:   d_nz (or d_nnz) and o_nz (or o_nnz).  By setting these parameters accurately,
    2327:   performance can be increased by more than a factor of 50.

    2329:   Collective

    2331:   Input Parameters:
    2332: + B     - the matrix
    2333: . bs    - size of block, the blocks are ALWAYS square. One can use MatSetBlockSizes() to set a different row and column blocksize but the row
    2334:           blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs()
    2335: . d_nz  - number of block nonzeros per block row in diagonal portion of local
    2336:           submatrix  (same for all local rows)
    2337: . d_nnz - array containing the number of block nonzeros in the various block rows
    2338:           in the upper triangular and diagonal part of the in diagonal portion of the local
    2339:           (possibly different for each block row) or `NULL`.  If you plan to factor the matrix you must leave room
    2340:           for the diagonal entry and set a value even if it is zero.
    2341: . o_nz  - number of block nonzeros per block row in the off-diagonal portion of local
    2342:           submatrix (same for all local rows).
    2343: - o_nnz - array containing the number of nonzeros in the various block rows of the
    2344:           off-diagonal portion of the local submatrix that is right of the diagonal
    2345:           (possibly different for each block row) or `NULL`.

    2347:   Options Database Keys:
    2348: + -mat_no_unroll  - uses code that does not unroll the loops in the
    2349:                     block calculations (much slower)
    2350: - -mat_block_size - size of the blocks to use

    2352:   Level: intermediate

    2354:   Notes:

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

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

    2361:   Storage Information:
    2362:   For a square global matrix we define each processor's diagonal portion
    2363:   to be its local rows and the corresponding columns (a square submatrix);
    2364:   each processor's off-diagonal portion encompasses the remainder of the
    2365:   local matrix (a rectangular submatrix).

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

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

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

    2381: .vb
    2382:            0 1 2 3 4 5 6 7 8 9 10 11
    2383:           --------------------------
    2384:    row 3  |. . . d d d o o o o  o  o
    2385:    row 4  |. . . d d d o o o o  o  o
    2386:    row 5  |. . . d d d o o o o  o  o
    2387:           --------------------------
    2388: .ve

    2390:   Thus, any entries in the d locations are stored in the d (diagonal)
    2391:   submatrix, and any entries in the o locations are stored in the
    2392:   o (off-diagonal) submatrix.  Note that the d matrix is stored in
    2393:   `MATSEQSBAIJ` format and the o submatrix in `MATSEQBAIJ` format.

    2395:   Now `d_nz` should indicate the number of block nonzeros per row in the upper triangular
    2396:   plus the diagonal part of the d matrix,
    2397:   and `o_nz` should indicate the number of block nonzeros per row in the o matrix

    2399:   In general, for PDE problems in which most nonzeros are near the diagonal,
    2400:   one expects `d_nz` >> `o_nz`.

    2402: .seealso: [](ch_matrices), `Mat`, `MATMPISBAIJ`, `MATSBAIJ`, `MatCreate()`, `MatCreateSeqSBAIJ()`, `MatSetValues()`, `MatCreateBAIJ()`, `PetscSplitOwnership()`
    2403: @*/
    2404: PetscErrorCode MatMPISBAIJSetPreallocation(Mat B, PetscInt bs, PetscInt d_nz, const PetscInt d_nnz[], PetscInt o_nz, const PetscInt o_nnz[])
    2405: {
    2406:   PetscFunctionBegin;
    2410:   PetscTryMethod(B, "MatMPISBAIJSetPreallocation_C", (Mat, PetscInt, PetscInt, const PetscInt[], PetscInt, const PetscInt[]), (B, bs, d_nz, d_nnz, o_nz, o_nnz));
    2411:   PetscFunctionReturn(PETSC_SUCCESS);
    2412: }

    2414: // PetscClangLinter pragma disable: -fdoc-section-header-unknown
    2415: /*@
    2416:   MatCreateSBAIJ - Creates a sparse parallel matrix in symmetric block AIJ format, `MATSBAIJ`,
    2417:   (block compressed row).  For good matrix assembly performance
    2418:   the user should preallocate the matrix storage by setting the parameters
    2419:   `d_nz` (or `d_nnz`) and `o_nz` (or `o_nnz`).

    2421:   Collective

    2423:   Input Parameters:
    2424: + comm  - MPI communicator
    2425: . bs    - size of block, the blocks are ALWAYS square. One can use `MatSetBlockSizes()` to set a different row and column blocksize but the row
    2426:           blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with `MatCreateVecs()`
    2427: . m     - number of local rows (or `PETSC_DECIDE` to have calculated if `M` is given)
    2428:           This value should be the same as the local size used in creating the
    2429:           y vector for the matrix-vector product y = Ax.
    2430: . n     - number of local columns (or `PETSC_DECIDE` to have calculated if `N` is given)
    2431:           This value should be the same as the local size used in creating the
    2432:           x vector for the matrix-vector product y = Ax.
    2433: . M     - number of global rows (or `PETSC_DETERMINE` to have calculated if `m` is given)
    2434: . N     - number of global columns (or `PETSC_DETERMINE` to have calculated if `n` is given)
    2435: . d_nz  - number of block nonzeros per block row in diagonal portion of local
    2436:           submatrix (same for all local rows)
    2437: . d_nnz - array containing the number of block nonzeros in the various block rows
    2438:           in the upper triangular portion of the in diagonal portion of the local
    2439:           (possibly different for each block block row) or `NULL`.
    2440:           If you plan to factor the matrix you must leave room for the diagonal entry and
    2441:           set its value even if it is zero.
    2442: . o_nz  - number of block nonzeros per block row in the off-diagonal portion of local
    2443:           submatrix (same for all local rows).
    2444: - o_nnz - array containing the number of nonzeros in the various block rows of the
    2445:           off-diagonal portion of the local submatrix (possibly different for
    2446:           each block row) or `NULL`.

    2448:   Output Parameter:
    2449: . A - the matrix

    2451:   Options Database Keys:
    2452: + -mat_no_unroll  - uses code that does not unroll the loops in the
    2453:                     block calculations (much slower)
    2454: . -mat_block_size - size of the blocks to use
    2455: - -mat_mpi        - use the parallel matrix data structures even on one processor
    2456:                     (defaults to using SeqBAIJ format on one processor)

    2458:   Level: intermediate

    2460:   Notes:
    2461:   It is recommended that one use `MatCreateFromOptions()` or the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()`,
    2462:   MatXXXXSetPreallocation() paradigm instead of this routine directly.
    2463:   [MatXXXXSetPreallocation() is, for example, `MatSeqAIJSetPreallocation()`]

    2465:   The number of rows and columns must be divisible by blocksize.
    2466:   This matrix type does not support complex Hermitian operation.

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

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

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

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

    2479:   Storage Information:
    2480:   For a square global matrix we define each processor's diagonal portion
    2481:   to be its local rows and the corresponding columns (a square submatrix);
    2482:   each processor's off-diagonal portion encompasses the remainder of the
    2483:   local matrix (a rectangular submatrix).

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

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

    2494: .vb
    2495:            0 1 2 3 4 5 6 7 8 9 10 11
    2496:           --------------------------
    2497:    row 3  |. . . d d d o o o o  o  o
    2498:    row 4  |. . . d d d o o o o  o  o
    2499:    row 5  |. . . d d d o o o o  o  o
    2500:           --------------------------
    2501: .ve

    2503:   Thus, any entries in the d locations are stored in the d (diagonal)
    2504:   submatrix, and any entries in the o locations are stored in the
    2505:   o (off-diagonal) submatrix. Note that the d matrix is stored in
    2506:   `MATSEQSBAIJ` format and the o submatrix in `MATSEQBAIJ` format.

    2508:   Now `d_nz` should indicate the number of block nonzeros per row in the upper triangular
    2509:   plus the diagonal part of the d matrix,
    2510:   and `o_nz` should indicate the number of block nonzeros per row in the o matrix.
    2511:   In general, for PDE problems in which most nonzeros are near the diagonal,
    2512:   one expects `d_nz` >> `o_nz`.

    2514: .seealso: [](ch_matrices), `Mat`, `MATSBAIJ`, `MatCreate()`, `MatCreateSeqSBAIJ()`, `MatSetValues()`, `MatCreateBAIJ()`,
    2515:           `MatGetOwnershipRange()`, `MatGetOwnershipRanges()`, `MatGetOwnershipRangeColumn()`, `MatGetOwnershipRangesColumn()`, `PetscLayout`
    2516: @*/
    2517: PetscErrorCode MatCreateSBAIJ(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)
    2518: {
    2519:   PetscMPIInt size;

    2521:   PetscFunctionBegin;
    2522:   PetscCall(MatCreate(comm, A));
    2523:   PetscCall(MatSetSizes(*A, m, n, M, N));
    2524:   PetscCallMPI(MPI_Comm_size(comm, &size));
    2525:   if (size > 1) {
    2526:     PetscCall(MatSetType(*A, MATMPISBAIJ));
    2527:     PetscCall(MatMPISBAIJSetPreallocation(*A, bs, d_nz, d_nnz, o_nz, o_nnz));
    2528:   } else {
    2529:     PetscCall(MatSetType(*A, MATSEQSBAIJ));
    2530:     PetscCall(MatSeqSBAIJSetPreallocation(*A, bs, d_nz, d_nnz));
    2531:   }
    2532:   PetscFunctionReturn(PETSC_SUCCESS);
    2533: }

    2535: static PetscErrorCode MatDuplicate_MPISBAIJ(Mat matin, MatDuplicateOption cpvalues, Mat *newmat)
    2536: {
    2537:   Mat           mat;
    2538:   Mat_MPISBAIJ *a, *oldmat = (Mat_MPISBAIJ *)matin->data;
    2539:   PetscInt      len = 0, nt, bs = matin->rmap->bs, mbs = oldmat->mbs;
    2540:   PetscScalar  *array;

    2542:   PetscFunctionBegin;
    2543:   *newmat = NULL;

    2545:   PetscCall(MatCreate(PetscObjectComm((PetscObject)matin), &mat));
    2546:   PetscCall(MatSetSizes(mat, matin->rmap->n, matin->cmap->n, matin->rmap->N, matin->cmap->N));
    2547:   PetscCall(MatSetType(mat, ((PetscObject)matin)->type_name));
    2548:   PetscCall(PetscLayoutReference(matin->rmap, &mat->rmap));
    2549:   PetscCall(PetscLayoutReference(matin->cmap, &mat->cmap));

    2551:   if (matin->hash_active) PetscCall(MatSetUp(mat));
    2552:   else {
    2553:     mat->factortype   = matin->factortype;
    2554:     mat->preallocated = PETSC_TRUE;
    2555:     mat->assembled    = PETSC_TRUE;
    2556:     mat->insertmode   = NOT_SET_VALUES;

    2558:     a      = (Mat_MPISBAIJ *)mat->data;
    2559:     a->bs2 = oldmat->bs2;
    2560:     a->mbs = oldmat->mbs;
    2561:     a->nbs = oldmat->nbs;
    2562:     a->Mbs = oldmat->Mbs;
    2563:     a->Nbs = oldmat->Nbs;

    2565:     a->size         = oldmat->size;
    2566:     a->rank         = oldmat->rank;
    2567:     a->donotstash   = oldmat->donotstash;
    2568:     a->roworiented  = oldmat->roworiented;
    2569:     a->rowindices   = NULL;
    2570:     a->rowvalues    = NULL;
    2571:     a->getrowactive = PETSC_FALSE;
    2572:     a->barray       = NULL;
    2573:     a->rstartbs     = oldmat->rstartbs;
    2574:     a->rendbs       = oldmat->rendbs;
    2575:     a->cstartbs     = oldmat->cstartbs;
    2576:     a->cendbs       = oldmat->cendbs;

    2578:     /* hash table stuff */
    2579:     a->ht           = NULL;
    2580:     a->hd           = NULL;
    2581:     a->ht_size      = 0;
    2582:     a->ht_flag      = oldmat->ht_flag;
    2583:     a->ht_fact      = oldmat->ht_fact;
    2584:     a->ht_total_ct  = 0;
    2585:     a->ht_insert_ct = 0;

    2587:     PetscCall(PetscArraycpy(a->rangebs, oldmat->rangebs, a->size + 2));
    2588:     if (oldmat->colmap) {
    2589: #if defined(PETSC_USE_CTABLE)
    2590:       PetscCall(PetscHMapIDuplicate(oldmat->colmap, &a->colmap));
    2591: #else
    2592:       PetscCall(PetscMalloc1(a->Nbs, &a->colmap));
    2593:       PetscCall(PetscArraycpy(a->colmap, oldmat->colmap, a->Nbs));
    2594: #endif
    2595:     } else a->colmap = NULL;

    2597:     if (oldmat->garray && (len = ((Mat_SeqBAIJ *)oldmat->B->data)->nbs)) {
    2598:       PetscCall(PetscMalloc1(len, &a->garray));
    2599:       PetscCall(PetscArraycpy(a->garray, oldmat->garray, len));
    2600:     } else a->garray = NULL;

    2602:     PetscCall(MatStashCreate_Private(PetscObjectComm((PetscObject)matin), matin->rmap->bs, &mat->bstash));
    2603:     PetscCall(VecDuplicate(oldmat->lvec, &a->lvec));
    2604:     PetscCall(VecScatterCopy(oldmat->Mvctx, &a->Mvctx));

    2606:     PetscCall(VecDuplicate(oldmat->slvec0, &a->slvec0));
    2607:     PetscCall(VecDuplicate(oldmat->slvec1, &a->slvec1));

    2609:     PetscCall(VecGetLocalSize(a->slvec1, &nt));
    2610:     PetscCall(VecGetArray(a->slvec1, &array));
    2611:     PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, bs * mbs, array, &a->slvec1a));
    2612:     PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, nt - bs * mbs, array + bs * mbs, &a->slvec1b));
    2613:     PetscCall(VecRestoreArray(a->slvec1, &array));
    2614:     PetscCall(VecGetArray(a->slvec0, &array));
    2615:     PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, nt - bs * mbs, array + bs * mbs, &a->slvec0b));
    2616:     PetscCall(VecRestoreArray(a->slvec0, &array));

    2618:     /* ierr =  VecScatterCopy(oldmat->sMvctx,&a->sMvctx); - not written yet, replaced by the lazy trick: */
    2619:     PetscCall(PetscObjectReference((PetscObject)oldmat->sMvctx));
    2620:     a->sMvctx = oldmat->sMvctx;

    2622:     PetscCall(MatDuplicate(oldmat->A, cpvalues, &a->A));
    2623:     PetscCall(MatDuplicate(oldmat->B, cpvalues, &a->B));
    2624:   }
    2625:   PetscCall(PetscFunctionListDuplicate(((PetscObject)matin)->qlist, &((PetscObject)mat)->qlist));
    2626:   *newmat = mat;
    2627:   PetscFunctionReturn(PETSC_SUCCESS);
    2628: }

    2630: /* Used for both MPIBAIJ and MPISBAIJ matrices */
    2631: #define MatLoad_MPISBAIJ_Binary MatLoad_MPIBAIJ_Binary

    2633: static PetscErrorCode MatLoad_MPISBAIJ(Mat mat, PetscViewer viewer)
    2634: {
    2635:   PetscBool isbinary;

    2637:   PetscFunctionBegin;
    2638:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
    2639:   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);
    2640:   PetscCall(MatLoad_MPISBAIJ_Binary(mat, viewer));
    2641:   PetscFunctionReturn(PETSC_SUCCESS);
    2642: }

    2644: static PetscErrorCode MatGetRowMaxAbs_MPISBAIJ(Mat A, Vec v, PetscInt idx[])
    2645: {
    2646:   Mat_MPISBAIJ *a = (Mat_MPISBAIJ *)A->data;
    2647:   Mat_SeqBAIJ  *b = (Mat_SeqBAIJ *)a->B->data;
    2648:   PetscReal     atmp;
    2649:   PetscReal    *work, *svalues, *rvalues;
    2650:   PetscInt      i, bs, mbs, *bi, *bj, brow, j, ncols, krow, kcol, col, row, Mbs, bcol;
    2651:   PetscMPIInt   rank, size;
    2652:   PetscInt     *rowners_bs, count, source;
    2653:   PetscScalar  *va;
    2654:   MatScalar    *ba;
    2655:   MPI_Status    stat;

    2657:   PetscFunctionBegin;
    2658:   PetscCheck(!idx, PETSC_COMM_SELF, PETSC_ERR_SUP, "Send email to petsc-maint@mcs.anl.gov");
    2659:   PetscCall(MatGetRowMaxAbs(a->A, v, NULL));
    2660:   PetscCall(VecGetArray(v, &va));

    2662:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
    2663:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank));

    2665:   bs  = A->rmap->bs;
    2666:   mbs = a->mbs;
    2667:   Mbs = a->Mbs;
    2668:   ba  = b->a;
    2669:   bi  = b->i;
    2670:   bj  = b->j;

    2672:   /* find ownerships */
    2673:   rowners_bs = A->rmap->range;

    2675:   /* each proc creates an array to be distributed */
    2676:   PetscCall(PetscCalloc1(bs * Mbs, &work));

    2678:   /* row_max for B */
    2679:   if (rank != size - 1) {
    2680:     for (i = 0; i < mbs; i++) {
    2681:       ncols = bi[1] - bi[0];
    2682:       bi++;
    2683:       brow = bs * i;
    2684:       for (j = 0; j < ncols; j++) {
    2685:         bcol = bs * (*bj);
    2686:         for (kcol = 0; kcol < bs; kcol++) {
    2687:           col = bcol + kcol;           /* local col index */
    2688:           col += rowners_bs[rank + 1]; /* global col index */
    2689:           for (krow = 0; krow < bs; krow++) {
    2690:             atmp = PetscAbsScalar(*ba);
    2691:             ba++;
    2692:             row = brow + krow; /* local row index */
    2693:             if (PetscRealPart(va[row]) < atmp) va[row] = atmp;
    2694:             if (work[col] < atmp) work[col] = atmp;
    2695:           }
    2696:         }
    2697:         bj++;
    2698:       }
    2699:     }

    2701:     /* send values to its owners */
    2702:     for (PetscMPIInt dest = rank + 1; dest < size; dest++) {
    2703:       svalues = work + rowners_bs[dest];
    2704:       count   = rowners_bs[dest + 1] - rowners_bs[dest];
    2705:       PetscCallMPI(MPIU_Send(svalues, count, MPIU_REAL, dest, rank, PetscObjectComm((PetscObject)A)));
    2706:     }
    2707:   }

    2709:   /* receive values */
    2710:   if (rank) {
    2711:     rvalues = work;
    2712:     count   = rowners_bs[rank + 1] - rowners_bs[rank];
    2713:     for (source = 0; source < rank; source++) {
    2714:       PetscCallMPI(MPIU_Recv(rvalues, count, MPIU_REAL, MPI_ANY_SOURCE, MPI_ANY_TAG, PetscObjectComm((PetscObject)A), &stat));
    2715:       /* process values */
    2716:       for (i = 0; i < count; i++) {
    2717:         if (PetscRealPart(va[i]) < rvalues[i]) va[i] = rvalues[i];
    2718:       }
    2719:     }
    2720:   }

    2722:   PetscCall(VecRestoreArray(v, &va));
    2723:   PetscCall(PetscFree(work));
    2724:   PetscFunctionReturn(PETSC_SUCCESS);
    2725: }

    2727: static PetscErrorCode MatSOR_MPISBAIJ(Mat matin, Vec bb, PetscReal omega, MatSORType flag, PetscReal fshift, PetscInt its, PetscInt lits, Vec xx)
    2728: {
    2729:   Mat_MPISBAIJ      *mat = (Mat_MPISBAIJ *)matin->data;
    2730:   PetscInt           mbs = mat->mbs, bs = matin->rmap->bs;
    2731:   PetscScalar       *x, *ptr, *from;
    2732:   Vec                bb1;
    2733:   const PetscScalar *b;

    2735:   PetscFunctionBegin;
    2736:   PetscCheck(its > 0 && lits > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Relaxation requires global its %" PetscInt_FMT " and local its %" PetscInt_FMT " both positive", its, lits);
    2737:   PetscCheck(bs <= 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "SSOR for block size > 1 is not yet implemented");

    2739:   if (flag == SOR_APPLY_UPPER) {
    2740:     PetscCall((*mat->A->ops->sor)(mat->A, bb, omega, flag, fshift, lits, 1, xx));
    2741:     PetscFunctionReturn(PETSC_SUCCESS);
    2742:   }

    2744:   if ((flag & SOR_LOCAL_SYMMETRIC_SWEEP) == SOR_LOCAL_SYMMETRIC_SWEEP) {
    2745:     if (flag & SOR_ZERO_INITIAL_GUESS) {
    2746:       PetscCall((*mat->A->ops->sor)(mat->A, bb, omega, flag, fshift, lits, lits, xx));
    2747:       its--;
    2748:     }

    2750:     PetscCall(VecDuplicate(bb, &bb1));
    2751:     while (its--) {
    2752:       /* lower triangular part: slvec0b = - B^T*xx */
    2753:       PetscCall((*mat->B->ops->multtranspose)(mat->B, xx, mat->slvec0b));

    2755:       /* copy xx into slvec0a */
    2756:       PetscCall(VecGetArray(mat->slvec0, &ptr));
    2757:       PetscCall(VecGetArray(xx, &x));
    2758:       PetscCall(PetscArraycpy(ptr, x, bs * mbs));
    2759:       PetscCall(VecRestoreArray(mat->slvec0, &ptr));

    2761:       PetscCall(VecScale(mat->slvec0, -1.0));

    2763:       /* copy bb into slvec1a */
    2764:       PetscCall(VecGetArray(mat->slvec1, &ptr));
    2765:       PetscCall(VecGetArrayRead(bb, &b));
    2766:       PetscCall(PetscArraycpy(ptr, b, bs * mbs));
    2767:       PetscCall(VecRestoreArray(mat->slvec1, &ptr));

    2769:       /* set slvec1b = 0 */
    2770:       PetscCall(PetscObjectStateIncrease((PetscObject)mat->slvec1b));
    2771:       PetscCall(VecZeroEntries(mat->slvec1b));

    2773:       PetscCall(VecScatterBegin(mat->sMvctx, mat->slvec0, mat->slvec1, ADD_VALUES, SCATTER_FORWARD));
    2774:       PetscCall(VecRestoreArray(xx, &x));
    2775:       PetscCall(VecRestoreArrayRead(bb, &b));
    2776:       PetscCall(VecScatterEnd(mat->sMvctx, mat->slvec0, mat->slvec1, ADD_VALUES, SCATTER_FORWARD));

    2778:       /* upper triangular part: bb1 = bb1 - B*x */
    2779:       PetscCall((*mat->B->ops->multadd)(mat->B, mat->slvec1b, mat->slvec1a, bb1));

    2781:       /* local diagonal sweep */
    2782:       PetscCall((*mat->A->ops->sor)(mat->A, bb1, omega, SOR_SYMMETRIC_SWEEP, fshift, lits, lits, xx));
    2783:     }
    2784:     PetscCall(VecDestroy(&bb1));
    2785:   } else if ((flag & SOR_LOCAL_FORWARD_SWEEP) && (its == 1) && (flag & SOR_ZERO_INITIAL_GUESS)) {
    2786:     PetscCall((*mat->A->ops->sor)(mat->A, bb, omega, flag, fshift, lits, 1, xx));
    2787:   } else if ((flag & SOR_LOCAL_BACKWARD_SWEEP) && (its == 1) && (flag & SOR_ZERO_INITIAL_GUESS)) {
    2788:     PetscCall((*mat->A->ops->sor)(mat->A, bb, omega, flag, fshift, lits, 1, xx));
    2789:   } else if (flag & SOR_EISENSTAT) {
    2790:     Vec                xx1;
    2791:     PetscBool          hasop;
    2792:     const PetscScalar *diag;
    2793:     PetscScalar       *sl, scale = (omega - 2.0) / omega;
    2794:     PetscInt           i, n;

    2796:     if (!mat->xx1) {
    2797:       PetscCall(VecDuplicate(bb, &mat->xx1));
    2798:       PetscCall(VecDuplicate(bb, &mat->bb1));
    2799:     }
    2800:     xx1 = mat->xx1;
    2801:     bb1 = mat->bb1;

    2803:     PetscCall((*mat->A->ops->sor)(mat->A, bb, omega, (MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_BACKWARD_SWEEP), fshift, lits, 1, xx));

    2805:     if (!mat->diag) {
    2806:       /* this is wrong for same matrix with new nonzero values */
    2807:       PetscCall(MatCreateVecs(matin, &mat->diag, NULL));
    2808:       PetscCall(MatGetDiagonal(matin, mat->diag));
    2809:     }
    2810:     PetscCall(MatHasOperation(matin, MATOP_MULT_DIAGONAL_BLOCK, &hasop));

    2812:     if (hasop) {
    2813:       PetscCall(MatMultDiagonalBlock(matin, xx, bb1));
    2814:       PetscCall(VecAYPX(mat->slvec1a, scale, bb));
    2815:     } else {
    2816:       /*
    2817:           These two lines are replaced by code that may be a bit faster for a good compiler
    2818:       PetscCall(VecPointwiseMult(mat->slvec1a,mat->diag,xx));
    2819:       PetscCall(VecAYPX(mat->slvec1a,scale,bb));
    2820:       */
    2821:       PetscCall(VecGetArray(mat->slvec1a, &sl));
    2822:       PetscCall(VecGetArrayRead(mat->diag, &diag));
    2823:       PetscCall(VecGetArrayRead(bb, &b));
    2824:       PetscCall(VecGetArray(xx, &x));
    2825:       PetscCall(VecGetLocalSize(xx, &n));
    2826:       if (omega == 1.0) {
    2827:         for (i = 0; i < n; i++) sl[i] = b[i] - diag[i] * x[i];
    2828:         PetscCall(PetscLogFlops(2.0 * n));
    2829:       } else {
    2830:         for (i = 0; i < n; i++) sl[i] = b[i] + scale * diag[i] * x[i];
    2831:         PetscCall(PetscLogFlops(3.0 * n));
    2832:       }
    2833:       PetscCall(VecRestoreArray(mat->slvec1a, &sl));
    2834:       PetscCall(VecRestoreArrayRead(mat->diag, &diag));
    2835:       PetscCall(VecRestoreArrayRead(bb, &b));
    2836:       PetscCall(VecRestoreArray(xx, &x));
    2837:     }

    2839:     /* multiply off-diagonal portion of matrix */
    2840:     PetscCall(PetscObjectStateIncrease((PetscObject)mat->slvec1b));
    2841:     PetscCall(VecZeroEntries(mat->slvec1b));
    2842:     PetscCall((*mat->B->ops->multtranspose)(mat->B, xx, mat->slvec0b));
    2843:     PetscCall(VecGetArray(mat->slvec0, &from));
    2844:     PetscCall(VecGetArray(xx, &x));
    2845:     PetscCall(PetscArraycpy(from, x, bs * mbs));
    2846:     PetscCall(VecRestoreArray(mat->slvec0, &from));
    2847:     PetscCall(VecRestoreArray(xx, &x));
    2848:     PetscCall(VecScatterBegin(mat->sMvctx, mat->slvec0, mat->slvec1, ADD_VALUES, SCATTER_FORWARD));
    2849:     PetscCall(VecScatterEnd(mat->sMvctx, mat->slvec0, mat->slvec1, ADD_VALUES, SCATTER_FORWARD));
    2850:     PetscCall((*mat->B->ops->multadd)(mat->B, mat->slvec1b, mat->slvec1a, mat->slvec1a));

    2852:     /* local sweep */
    2853:     PetscCall((*mat->A->ops->sor)(mat->A, mat->slvec1a, omega, (MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_FORWARD_SWEEP), fshift, lits, 1, xx1));
    2854:     PetscCall(VecAXPY(xx, 1.0, xx1));
    2855:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "MatSORType is not supported for SBAIJ matrix format");
    2856:   PetscFunctionReturn(PETSC_SUCCESS);
    2857: }

    2859: /*@
    2860:   MatCreateMPISBAIJWithArrays - creates a `MATMPISBAIJ` matrix using arrays that contain in standard CSR format for the local rows.

    2862:   Collective

    2864:   Input Parameters:
    2865: + comm - MPI communicator
    2866: . bs   - the block size, only a block size of 1 is supported
    2867: . m    - number of local rows (Cannot be `PETSC_DECIDE`)
    2868: . n    - This value should be the same as the local size used in creating the
    2869:          x vector for the matrix-vector product $ y = Ax $. (or `PETSC_DECIDE` to have
    2870:          calculated if `N` is given) For square matrices `n` is almost always `m`.
    2871: . M    - number of global rows (or `PETSC_DETERMINE` to have calculated if `m` is given)
    2872: . N    - number of global columns (or `PETSC_DETERMINE` to have calculated if `n` is given)
    2873: . i    - row indices; that is i[0] = 0, i[row] = i[row-1] + number of block elements in that row block row of the matrix
    2874: . j    - column indices
    2875: - a    - matrix values

    2877:   Output Parameter:
    2878: . mat - the matrix

    2880:   Level: intermediate

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

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

    2889: .seealso: [](ch_matrices), `Mat`, `MATMPISBAIJ`, `MatCreate()`, `MatCreateSeqAIJ()`, `MatSetValues()`, `MatMPIAIJSetPreallocation()`, `MatMPIAIJSetPreallocationCSR()`,
    2890:           `MATMPIAIJ`, `MatCreateAIJ()`, `MatCreateMPIAIJWithSplitArrays()`, `MatMPISBAIJSetPreallocationCSR()`
    2891: @*/
    2892: PetscErrorCode MatCreateMPISBAIJWithArrays(MPI_Comm comm, PetscInt bs, PetscInt m, PetscInt n, PetscInt M, PetscInt N, const PetscInt i[], const PetscInt j[], const PetscScalar a[], Mat *mat)
    2893: {
    2894:   PetscFunctionBegin;
    2895:   PetscCheck(!i[0], PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "i (row indices) must start with 0");
    2896:   PetscCheck(m >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "local number of rows (m) cannot be PETSC_DECIDE, or negative");
    2897:   PetscCall(MatCreate(comm, mat));
    2898:   PetscCall(MatSetSizes(*mat, m, n, M, N));
    2899:   PetscCall(MatSetType(*mat, MATMPISBAIJ));
    2900:   PetscCall(MatMPISBAIJSetPreallocationCSR(*mat, bs, i, j, a));
    2901:   PetscFunctionReturn(PETSC_SUCCESS);
    2902: }

    2904: /*@
    2905:   MatMPISBAIJSetPreallocationCSR - Creates a sparse parallel matrix in `MATMPISBAIJ` format using the given nonzero structure and (optional) numerical values

    2907:   Collective

    2909:   Input Parameters:
    2910: + B  - the matrix
    2911: . bs - the block size
    2912: . i  - the indices into `j` for the start of each local row (indices start with zero)
    2913: . j  - the column indices for each local row (indices start with zero) these must be sorted for each row
    2914: - v  - optional values in the matrix, pass `NULL` if not provided

    2916:   Level: advanced

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

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

    2926:   Any entries passed in that are below the diagonal are ignored

    2928: .seealso: [](ch_matrices), `Mat`, `MATMPISBAIJ`, `MatCreate()`, `MatCreateSeqAIJ()`, `MatSetValues()`, `MatMPIBAIJSetPreallocation()`, `MatCreateAIJ()`, `MATMPIAIJ`,
    2929:           `MatCreateMPISBAIJWithArrays()`
    2930: @*/
    2931: PetscErrorCode MatMPISBAIJSetPreallocationCSR(Mat B, PetscInt bs, const PetscInt i[], const PetscInt j[], const PetscScalar v[])
    2932: {
    2933:   PetscFunctionBegin;
    2934:   PetscTryMethod(B, "MatMPISBAIJSetPreallocationCSR_C", (Mat, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[]), (B, bs, i, j, v));
    2935:   PetscFunctionReturn(PETSC_SUCCESS);
    2936: }

    2938: PetscErrorCode MatCreateMPIMatConcatenateSeqMat_MPISBAIJ(MPI_Comm comm, Mat inmat, PetscInt n, MatReuse scall, Mat *outmat)
    2939: {
    2940:   PetscInt     m, N, i, rstart, nnz, Ii, bs, cbs;
    2941:   PetscInt    *indx;
    2942:   PetscScalar *values;

    2944:   PetscFunctionBegin;
    2945:   PetscCall(MatGetSize(inmat, &m, &N));
    2946:   if (scall == MAT_INITIAL_MATRIX) { /* symbolic phase */
    2947:     Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)inmat->data;
    2948:     PetscInt     *dnz, *onz, mbs, Nbs, nbs;
    2949:     PetscInt     *bindx, rmax = a->rmax, j;
    2950:     PetscMPIInt   rank, size;

    2952:     PetscCall(MatGetBlockSizes(inmat, &bs, &cbs));
    2953:     mbs = m / bs;
    2954:     Nbs = N / cbs;
    2955:     if (n == PETSC_DECIDE) PetscCall(PetscSplitOwnershipBlock(comm, cbs, &n, &N));
    2956:     nbs = n / cbs;

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

    2961:     PetscCallMPI(MPI_Comm_rank(comm, &rank));
    2962:     PetscCallMPI(MPI_Comm_size(comm, &size));
    2963:     if (rank == size - 1) {
    2964:       /* Check sum(nbs) = Nbs */
    2965:       PetscCheck(__end == Nbs, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Sum of local block columns %" PetscInt_FMT " != global block columns %" PetscInt_FMT, __end, Nbs);
    2966:     }

    2968:     rstart = __rstart; /* block rstart of *outmat; see inline function MatPreallocateBegin */
    2969:     PetscCall(MatSetOption(inmat, MAT_GETROW_UPPERTRIANGULAR, PETSC_TRUE));
    2970:     for (i = 0; i < mbs; i++) {
    2971:       PetscCall(MatGetRow_SeqSBAIJ(inmat, i * bs, &nnz, &indx, NULL)); /* non-blocked nnz and indx */
    2972:       nnz = nnz / bs;
    2973:       for (j = 0; j < nnz; j++) bindx[j] = indx[j * bs] / bs;
    2974:       PetscCall(MatPreallocateSet(i + rstart, nnz, bindx, dnz, onz));
    2975:       PetscCall(MatRestoreRow_SeqSBAIJ(inmat, i * bs, &nnz, &indx, NULL));
    2976:     }
    2977:     PetscCall(MatSetOption(inmat, MAT_GETROW_UPPERTRIANGULAR, PETSC_FALSE));
    2978:     PetscCall(PetscFree(bindx));

    2980:     PetscCall(MatCreate(comm, outmat));
    2981:     PetscCall(MatSetSizes(*outmat, m, n, PETSC_DETERMINE, PETSC_DETERMINE));
    2982:     PetscCall(MatSetBlockSizes(*outmat, bs, cbs));
    2983:     PetscCall(MatSetType(*outmat, MATSBAIJ));
    2984:     PetscCall(MatSeqSBAIJSetPreallocation(*outmat, bs, 0, dnz));
    2985:     PetscCall(MatMPISBAIJSetPreallocation(*outmat, bs, 0, dnz, 0, onz));
    2986:     MatPreallocateEnd(dnz, onz);
    2987:   }

    2989:   /* numeric phase */
    2990:   PetscCall(MatGetBlockSizes(inmat, &bs, &cbs));
    2991:   PetscCall(MatGetOwnershipRange(*outmat, &rstart, NULL));

    2993:   PetscCall(MatSetOption(inmat, MAT_GETROW_UPPERTRIANGULAR, PETSC_TRUE));
    2994:   for (i = 0; i < m; i++) {
    2995:     PetscCall(MatGetRow_SeqSBAIJ(inmat, i, &nnz, &indx, &values));
    2996:     Ii = i + rstart;
    2997:     PetscCall(MatSetValues(*outmat, 1, &Ii, nnz, indx, values, INSERT_VALUES));
    2998:     PetscCall(MatRestoreRow_SeqSBAIJ(inmat, i, &nnz, &indx, &values));
    2999:   }
    3000:   PetscCall(MatSetOption(inmat, MAT_GETROW_UPPERTRIANGULAR, PETSC_FALSE));
    3001:   PetscCall(MatAssemblyBegin(*outmat, MAT_FINAL_ASSEMBLY));
    3002:   PetscCall(MatAssemblyEnd(*outmat, MAT_FINAL_ASSEMBLY));
    3003:   PetscFunctionReturn(PETSC_SUCCESS);
    3004: }