Actual source code: shell.c

      1: /*
      2:    This provides a simple shell for Fortran (and C programmers) to
      3:   create a very simple matrix class for use with KSP without coding
      4:   much of anything.
      5: */

      7: #include <../src/mat/impls/shell/shell.h>

      9: /*
     10:      Store and scale values on zeroed rows
     11:      xx = [x_1, 0], 0 on zeroed columns
     12: */
     13: static PetscErrorCode MatShellPreZeroRight(Mat A, Vec x, Vec *xx)
     14: {
     15:   Mat_Shell *shell = (Mat_Shell *)A->data;

     17:   PetscFunctionBegin;
     18:   *xx = x;
     19:   if (shell->zrows) {
     20:     PetscCall(VecSet(shell->zvals_w, 0.0));
     21:     PetscCall(VecScatterBegin(shell->zvals_sct_c, x, shell->zvals_w, INSERT_VALUES, SCATTER_FORWARD));
     22:     PetscCall(VecScatterEnd(shell->zvals_sct_c, x, shell->zvals_w, INSERT_VALUES, SCATTER_FORWARD));
     23:     PetscCall(VecPointwiseMult(shell->zvals_w, shell->zvals_w, shell->zvals));
     24:   }
     25:   if (shell->zcols) {
     26:     if (!shell->right_work) PetscCall(MatCreateVecs(A, &shell->right_work, NULL));
     27:     PetscCall(VecCopy(x, shell->right_work));
     28:     PetscCall(VecISSet(shell->right_work, shell->zcols, 0.0));
     29:     *xx = shell->right_work;
     30:   }
     31:   PetscFunctionReturn(PETSC_SUCCESS);
     32: }

     34: /* Insert properly diagonally scaled values stored in MatShellPreZeroRight */
     35: static PetscErrorCode MatShellPostZeroLeft(Mat A, Vec x)
     36: {
     37:   Mat_Shell *shell = (Mat_Shell *)A->data;

     39:   PetscFunctionBegin;
     40:   if (shell->zrows) {
     41:     PetscCall(VecScatterBegin(shell->zvals_sct_r, shell->zvals_w, x, INSERT_VALUES, SCATTER_REVERSE));
     42:     PetscCall(VecScatterEnd(shell->zvals_sct_r, shell->zvals_w, x, INSERT_VALUES, SCATTER_REVERSE));
     43:   }
     44:   PetscFunctionReturn(PETSC_SUCCESS);
     45: }

     47: /*
     48:      Store and scale values on zeroed rows
     49:      xx = [x_1, 0], 0 on zeroed rows
     50: */
     51: static PetscErrorCode MatShellPreZeroLeft(Mat A, Vec x, Vec *xx)
     52: {
     53:   Mat_Shell *shell = (Mat_Shell *)A->data;

     55:   PetscFunctionBegin;
     56:   *xx = NULL;
     57:   if (!shell->zrows) {
     58:     *xx = x;
     59:   } else {
     60:     if (!shell->left_work) PetscCall(MatCreateVecs(A, NULL, &shell->left_work));
     61:     PetscCall(VecCopy(x, shell->left_work));
     62:     PetscCall(VecSet(shell->zvals_w, 0.0));
     63:     PetscCall(VecScatterBegin(shell->zvals_sct_r, shell->zvals_w, shell->left_work, INSERT_VALUES, SCATTER_REVERSE));
     64:     PetscCall(VecScatterEnd(shell->zvals_sct_r, shell->zvals_w, shell->left_work, INSERT_VALUES, SCATTER_REVERSE));
     65:     PetscCall(VecScatterBegin(shell->zvals_sct_r, x, shell->zvals_w, INSERT_VALUES, SCATTER_FORWARD));
     66:     PetscCall(VecScatterEnd(shell->zvals_sct_r, x, shell->zvals_w, INSERT_VALUES, SCATTER_FORWARD));
     67:     PetscCall(VecPointwiseMult(shell->zvals_w, shell->zvals_w, shell->zvals));
     68:     *xx = shell->left_work;
     69:   }
     70:   PetscFunctionReturn(PETSC_SUCCESS);
     71: }

     73: /* Zero zero-columns contributions, sum contributions from properly scaled values stored in MatShellPreZeroLeft */
     74: static PetscErrorCode MatShellPostZeroRight(Mat A, Vec x)
     75: {
     76:   Mat_Shell *shell = (Mat_Shell *)A->data;

     78:   PetscFunctionBegin;
     79:   if (shell->zcols) PetscCall(VecISSet(x, shell->zcols, 0.0));
     80:   if (shell->zrows) {
     81:     PetscCall(VecScatterBegin(shell->zvals_sct_c, shell->zvals_w, x, ADD_VALUES, SCATTER_REVERSE));
     82:     PetscCall(VecScatterEnd(shell->zvals_sct_c, shell->zvals_w, x, ADD_VALUES, SCATTER_REVERSE));
     83:   }
     84:   PetscFunctionReturn(PETSC_SUCCESS);
     85: }

     87: /*
     88:       xx = diag(left)*x
     89: */
     90: static PetscErrorCode MatShellPreScaleLeft(Mat A, Vec x, Vec *xx, PetscBool conjugate)
     91: {
     92:   Mat_Shell *shell = (Mat_Shell *)A->data;

     94:   PetscFunctionBegin;
     95:   *xx = NULL;
     96:   if (!shell->left) {
     97:     *xx = x;
     98:   } else {
     99:     if (!shell->left_work) PetscCall(VecDuplicate(shell->left, &shell->left_work));
    100:     if (conjugate) { /* get arrays because there is no VecPointwiseMultConj() */
    101:       PetscInt           i, m;
    102:       const PetscScalar *d, *xarray;
    103:       PetscScalar       *w;
    104:       PetscCall(VecGetLocalSize(x, &m));
    105:       PetscCall(VecGetArrayRead(shell->left, &d));
    106:       PetscCall(VecGetArrayRead(x, &xarray));
    107:       PetscCall(VecGetArrayWrite(shell->left_work, &w));
    108:       for (i = 0; i < m; i++) w[i] = PetscConj(d[i]) * xarray[i];
    109:       PetscCall(VecRestoreArrayRead(shell->dshift, &d));
    110:       PetscCall(VecRestoreArrayRead(x, &xarray));
    111:       PetscCall(VecRestoreArrayWrite(shell->left_work, &w));
    112:     } else PetscCall(VecPointwiseMult(shell->left_work, x, shell->left));
    113:     *xx = shell->left_work;
    114:   }
    115:   PetscFunctionReturn(PETSC_SUCCESS);
    116: }

    118: /*
    119:      xx = diag(right)*x
    120: */
    121: static PetscErrorCode MatShellPreScaleRight(Mat A, Vec x, Vec *xx)
    122: {
    123:   Mat_Shell *shell = (Mat_Shell *)A->data;

    125:   PetscFunctionBegin;
    126:   *xx = NULL;
    127:   if (!shell->right) {
    128:     *xx = x;
    129:   } else {
    130:     if (!shell->right_work) PetscCall(VecDuplicate(shell->right, &shell->right_work));
    131:     PetscCall(VecPointwiseMult(shell->right_work, x, shell->right));
    132:     *xx = shell->right_work;
    133:   }
    134:   PetscFunctionReturn(PETSC_SUCCESS);
    135: }

    137: /*
    138:     x = diag(left)*x
    139: */
    140: static PetscErrorCode MatShellPostScaleLeft(Mat A, Vec x)
    141: {
    142:   Mat_Shell *shell = (Mat_Shell *)A->data;

    144:   PetscFunctionBegin;
    145:   if (shell->left) PetscCall(VecPointwiseMult(x, x, shell->left));
    146:   PetscFunctionReturn(PETSC_SUCCESS);
    147: }

    149: /*
    150:     x = diag(right)*x
    151: */
    152: static PetscErrorCode MatShellPostScaleRight(Mat A, Vec x, PetscBool conjugate)
    153: {
    154:   Mat_Shell *shell = (Mat_Shell *)A->data;

    156:   PetscFunctionBegin;
    157:   if (shell->right) {
    158:     if (conjugate) { /* get arrays because there is no VecPointwiseMultConj() */
    159:       PetscInt           i, m;
    160:       const PetscScalar *d;
    161:       PetscScalar       *xarray;
    162:       PetscCall(VecGetLocalSize(x, &m));
    163:       PetscCall(VecGetArrayRead(shell->right, &d));
    164:       PetscCall(VecGetArray(x, &xarray));
    165:       for (i = 0; i < m; i++) xarray[i] = PetscConj(d[i]) * xarray[i];
    166:       PetscCall(VecRestoreArrayRead(shell->dshift, &d));
    167:       PetscCall(VecRestoreArray(x, &xarray));
    168:     } else PetscCall(VecPointwiseMult(x, x, shell->right));
    169:   }
    170:   PetscFunctionReturn(PETSC_SUCCESS);
    171: }

    173: /*
    174:          Y = vscale*Y + diag(dshift)*X + vshift*X

    176:          On input Y already contains A*x

    178:          If conjugate=PETSC_TRUE then vscale, dshift, and vshift are conjugated
    179: */
    180: static PetscErrorCode MatShellShiftAndScale(Mat A, Vec X, Vec Y, PetscBool conjugate)
    181: {
    182:   Mat_Shell  *shell  = (Mat_Shell *)A->data;
    183:   PetscScalar vscale = conjugate ? PetscConj(shell->vscale) : shell->vscale;
    184:   PetscScalar vshift = conjugate ? PetscConj(shell->vshift) : shell->vshift;

    186:   PetscFunctionBegin;
    187:   if (shell->dshift) { /* get arrays because there is no VecPointwiseMultAdd() */
    188:     PetscInt           i, m;
    189:     const PetscScalar *x, *d;
    190:     PetscScalar       *y;
    191:     PetscCall(VecGetLocalSize(X, &m));
    192:     PetscCall(VecGetArrayRead(shell->dshift, &d));
    193:     PetscCall(VecGetArrayRead(X, &x));
    194:     PetscCall(VecGetArray(Y, &y));
    195:     if (conjugate)
    196:       for (i = 0; i < m; i++) y[i] = vscale * y[i] + PetscConj(d[i]) * x[i];
    197:     else
    198:       for (i = 0; i < m; i++) y[i] = vscale * y[i] + d[i] * x[i];
    199:     PetscCall(VecRestoreArrayRead(shell->dshift, &d));
    200:     PetscCall(VecRestoreArrayRead(X, &x));
    201:     PetscCall(VecRestoreArray(Y, &y));
    202:   } else {
    203:     PetscCall(VecScale(Y, vscale));
    204:   }
    205:   if (vshift != 0.0) PetscCall(VecAXPY(Y, vshift, X)); /* if test is for non-square matrices */
    206:   PetscFunctionReturn(PETSC_SUCCESS);
    207: }

    209: static PetscErrorCode MatShellGetContext_Shell(Mat mat, void *ctx)
    210: {
    211:   Mat_Shell *shell = (Mat_Shell *)mat->data;

    213:   PetscFunctionBegin;
    214:   if (shell->ctxcontainer) PetscCall(PetscContainerGetPointer(shell->ctxcontainer, (void **)ctx));
    215:   else *(void **)ctx = NULL;
    216:   PetscFunctionReturn(PETSC_SUCCESS);
    217: }

    219: /*@
    220:   MatShellGetContext - Returns the user-provided context associated with a `MATSHELL` shell matrix.

    222:   Not Collective

    224:   Input Parameter:
    225: . mat - the matrix, should have been created with `MatCreateShell()`

    227:   Output Parameter:
    228: . ctx - the user provided context

    230:   Level: advanced

    232:   Fortran Notes:
    233:   You must write a Fortran interface definition for this
    234:   function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.

    236: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellSetOperation()`, `MatShellSetContext()`
    237: @*/
    238: PetscErrorCode MatShellGetContext(Mat mat, void *ctx)
    239: {
    240:   PetscFunctionBegin;
    242:   PetscAssertPointer(ctx, 2);
    243:   PetscUseMethod(mat, "MatShellGetContext_C", (Mat, void *), (mat, ctx));
    244:   PetscFunctionReturn(PETSC_SUCCESS);
    245: }

    247: static PetscErrorCode MatZeroRowsColumns_Local_Shell(Mat mat, PetscInt nr, PetscInt rows[], PetscInt nc, PetscInt cols[], PetscScalar diag, PetscBool rc)
    248: {
    249:   Mat_Shell      *shell = (Mat_Shell *)mat->data;
    250:   Vec             x = NULL, b = NULL;
    251:   IS              is1, is2;
    252:   const PetscInt *ridxs;
    253:   PetscInt       *idxs, *gidxs;
    254:   PetscInt        cum, rst, cst, i;

    256:   PetscFunctionBegin;
    257:   if (!shell->zvals) PetscCall(MatCreateVecs(mat, NULL, &shell->zvals));
    258:   if (!shell->zvals_w) PetscCall(VecDuplicate(shell->zvals, &shell->zvals_w));
    259:   PetscCall(MatGetOwnershipRange(mat, &rst, NULL));
    260:   PetscCall(MatGetOwnershipRangeColumn(mat, &cst, NULL));

    262:   /* Expand/create index set of zeroed rows */
    263:   PetscCall(PetscMalloc1(nr, &idxs));
    264:   for (i = 0; i < nr; i++) idxs[i] = rows[i] + rst;
    265:   PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)mat), nr, idxs, PETSC_OWN_POINTER, &is1));
    266:   PetscCall(ISSort(is1));
    267:   PetscCall(VecISSet(shell->zvals, is1, diag));
    268:   if (shell->zrows) {
    269:     PetscCall(ISSum(shell->zrows, is1, &is2));
    270:     PetscCall(ISDestroy(&shell->zrows));
    271:     PetscCall(ISDestroy(&is1));
    272:     shell->zrows = is2;
    273:   } else shell->zrows = is1;

    275:   /* Create scatters for diagonal values communications */
    276:   PetscCall(VecScatterDestroy(&shell->zvals_sct_c));
    277:   PetscCall(VecScatterDestroy(&shell->zvals_sct_r));

    279:   /* row scatter: from/to left vector */
    280:   PetscCall(MatCreateVecs(mat, &x, &b));
    281:   PetscCall(VecScatterCreate(b, shell->zrows, shell->zvals_w, shell->zrows, &shell->zvals_sct_r));

    283:   /* col scatter: from right vector to left vector */
    284:   PetscCall(ISGetIndices(shell->zrows, &ridxs));
    285:   PetscCall(ISGetLocalSize(shell->zrows, &nr));
    286:   PetscCall(PetscMalloc1(nr, &gidxs));
    287:   for (i = 0, cum = 0; i < nr; i++) {
    288:     if (ridxs[i] >= mat->cmap->N) continue;
    289:     gidxs[cum] = ridxs[i];
    290:     cum++;
    291:   }
    292:   PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)mat), cum, gidxs, PETSC_OWN_POINTER, &is1));
    293:   PetscCall(VecScatterCreate(x, is1, shell->zvals_w, is1, &shell->zvals_sct_c));
    294:   PetscCall(ISDestroy(&is1));
    295:   PetscCall(VecDestroy(&x));
    296:   PetscCall(VecDestroy(&b));

    298:   /* Expand/create index set of zeroed columns */
    299:   if (rc) {
    300:     PetscCall(PetscMalloc1(nc, &idxs));
    301:     for (i = 0; i < nc; i++) idxs[i] = cols[i] + cst;
    302:     PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)mat), nc, idxs, PETSC_OWN_POINTER, &is1));
    303:     PetscCall(ISSort(is1));
    304:     if (shell->zcols) {
    305:       PetscCall(ISSum(shell->zcols, is1, &is2));
    306:       PetscCall(ISDestroy(&shell->zcols));
    307:       PetscCall(ISDestroy(&is1));
    308:       shell->zcols = is2;
    309:     } else shell->zcols = is1;
    310:   }
    311:   PetscFunctionReturn(PETSC_SUCCESS);
    312: }

    314: static PetscErrorCode MatZeroRows_Shell(Mat mat, PetscInt n, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
    315: {
    316:   Mat_Shell *shell = (Mat_Shell *)mat->data;
    317:   PetscInt   nr, *lrows;

    319:   PetscFunctionBegin;
    320:   if (x && b) {
    321:     Vec          xt;
    322:     PetscScalar *vals;
    323:     PetscInt    *gcols, i, st, nl, nc;

    325:     PetscCall(PetscMalloc1(n, &gcols));
    326:     for (i = 0, nc = 0; i < n; i++)
    327:       if (rows[i] < mat->cmap->N) gcols[nc++] = rows[i];

    329:     PetscCall(MatCreateVecs(mat, &xt, NULL));
    330:     PetscCall(VecCopy(x, xt));
    331:     PetscCall(PetscCalloc1(nc, &vals));
    332:     PetscCall(VecSetValues(xt, nc, gcols, vals, INSERT_VALUES)); /* xt = [x1, 0] */
    333:     PetscCall(PetscFree(vals));
    334:     PetscCall(VecAssemblyBegin(xt));
    335:     PetscCall(VecAssemblyEnd(xt));
    336:     PetscCall(VecAYPX(xt, -1.0, x)); /* xt = [0, x2] */

    338:     PetscCall(VecGetOwnershipRange(xt, &st, NULL));
    339:     PetscCall(VecGetLocalSize(xt, &nl));
    340:     PetscCall(VecGetArray(xt, &vals));
    341:     for (i = 0; i < nl; i++) {
    342:       PetscInt g = i + st;
    343:       if (g > mat->rmap->N) continue;
    344:       if (PetscAbsScalar(vals[i]) == 0.0) continue;
    345:       PetscCall(VecSetValue(b, g, diag * vals[i], INSERT_VALUES));
    346:     }
    347:     PetscCall(VecRestoreArray(xt, &vals));
    348:     PetscCall(VecAssemblyBegin(b));
    349:     PetscCall(VecAssemblyEnd(b)); /* b  = [b1, x2 * diag] */
    350:     PetscCall(VecDestroy(&xt));
    351:     PetscCall(PetscFree(gcols));
    352:   }
    353:   PetscCall(PetscLayoutMapLocal(mat->rmap, n, rows, &nr, &lrows, NULL));
    354:   PetscCall(MatZeroRowsColumns_Local_Shell(mat, nr, lrows, 0, NULL, diag, PETSC_FALSE));
    355:   if (shell->axpy) PetscCall(MatZeroRows(shell->axpy, n, rows, 0.0, NULL, NULL));
    356:   PetscCall(PetscFree(lrows));
    357:   PetscFunctionReturn(PETSC_SUCCESS);
    358: }

    360: static PetscErrorCode MatZeroRowsColumns_Shell(Mat mat, PetscInt n, const PetscInt rowscols[], PetscScalar diag, Vec x, Vec b)
    361: {
    362:   Mat_Shell *shell = (Mat_Shell *)mat->data;
    363:   PetscInt  *lrows, *lcols;
    364:   PetscInt   nr, nc;
    365:   PetscBool  congruent;

    367:   PetscFunctionBegin;
    368:   if (x && b) {
    369:     Vec          xt, bt;
    370:     PetscScalar *vals;
    371:     PetscInt    *grows, *gcols, i, st, nl;

    373:     PetscCall(PetscMalloc2(n, &grows, n, &gcols));
    374:     for (i = 0, nr = 0; i < n; i++)
    375:       if (rowscols[i] < mat->rmap->N) grows[nr++] = rowscols[i];
    376:     for (i = 0, nc = 0; i < n; i++)
    377:       if (rowscols[i] < mat->cmap->N) gcols[nc++] = rowscols[i];
    378:     PetscCall(PetscCalloc1(n, &vals));

    380:     PetscCall(MatCreateVecs(mat, &xt, &bt));
    381:     PetscCall(VecCopy(x, xt));
    382:     PetscCall(VecSetValues(xt, nc, gcols, vals, INSERT_VALUES)); /* xt = [x1, 0] */
    383:     PetscCall(VecAssemblyBegin(xt));
    384:     PetscCall(VecAssemblyEnd(xt));
    385:     PetscCall(VecAXPY(xt, -1.0, x));                             /* xt = [0, -x2] */
    386:     PetscCall(MatMult(mat, xt, bt));                             /* bt = [-A12*x2,-A22*x2] */
    387:     PetscCall(VecSetValues(bt, nr, grows, vals, INSERT_VALUES)); /* bt = [-A12*x2,0] */
    388:     PetscCall(VecAssemblyBegin(bt));
    389:     PetscCall(VecAssemblyEnd(bt));
    390:     PetscCall(VecAXPY(b, 1.0, bt));                              /* b  = [b1 - A12*x2, b2] */
    391:     PetscCall(VecSetValues(bt, nr, grows, vals, INSERT_VALUES)); /* b  = [b1 - A12*x2, 0] */
    392:     PetscCall(VecAssemblyBegin(bt));
    393:     PetscCall(VecAssemblyEnd(bt));
    394:     PetscCall(PetscFree(vals));

    396:     PetscCall(VecGetOwnershipRange(xt, &st, NULL));
    397:     PetscCall(VecGetLocalSize(xt, &nl));
    398:     PetscCall(VecGetArray(xt, &vals));
    399:     for (i = 0; i < nl; i++) {
    400:       PetscInt g = i + st;
    401:       if (g > mat->rmap->N) continue;
    402:       if (PetscAbsScalar(vals[i]) == 0.0) continue;
    403:       PetscCall(VecSetValue(b, g, -diag * vals[i], INSERT_VALUES));
    404:     }
    405:     PetscCall(VecRestoreArray(xt, &vals));
    406:     PetscCall(VecAssemblyBegin(b));
    407:     PetscCall(VecAssemblyEnd(b)); /* b  = [b1 - A12*x2, x2 * diag] */
    408:     PetscCall(VecDestroy(&xt));
    409:     PetscCall(VecDestroy(&bt));
    410:     PetscCall(PetscFree2(grows, gcols));
    411:   }
    412:   PetscCall(PetscLayoutMapLocal(mat->rmap, n, rowscols, &nr, &lrows, NULL));
    413:   PetscCall(MatHasCongruentLayouts(mat, &congruent));
    414:   if (congruent) {
    415:     nc    = nr;
    416:     lcols = lrows;
    417:   } else { /* MatZeroRowsColumns implicitly assumes the rowscols indices are for a square matrix, here we handle a more general case */
    418:     PetscInt i, nt, *t;

    420:     PetscCall(PetscMalloc1(n, &t));
    421:     for (i = 0, nt = 0; i < n; i++)
    422:       if (rowscols[i] < mat->cmap->N) t[nt++] = rowscols[i];
    423:     PetscCall(PetscLayoutMapLocal(mat->cmap, nt, t, &nc, &lcols, NULL));
    424:     PetscCall(PetscFree(t));
    425:   }
    426:   PetscCall(MatZeroRowsColumns_Local_Shell(mat, nr, lrows, nc, lcols, diag, PETSC_TRUE));
    427:   if (!congruent) PetscCall(PetscFree(lcols));
    428:   PetscCall(PetscFree(lrows));
    429:   if (shell->axpy) PetscCall(MatZeroRowsColumns(shell->axpy, n, rowscols, 0.0, NULL, NULL));
    430:   PetscFunctionReturn(PETSC_SUCCESS);
    431: }

    433: static PetscErrorCode MatDestroy_Shell(Mat mat)
    434: {
    435:   Mat_Shell              *shell = (Mat_Shell *)mat->data;
    436:   MatShellMatFunctionList matmat;

    438:   PetscFunctionBegin;
    439:   if (shell->ops->destroy) PetscCall((*shell->ops->destroy)(mat));
    440:   PetscCall(PetscMemzero(shell->ops, sizeof(struct _MatShellOps)));
    441:   PetscCall(VecDestroy(&shell->left));
    442:   PetscCall(VecDestroy(&shell->right));
    443:   PetscCall(VecDestroy(&shell->dshift));
    444:   PetscCall(VecDestroy(&shell->left_work));
    445:   PetscCall(VecDestroy(&shell->right_work));
    446:   PetscCall(VecDestroy(&shell->left_add_work));
    447:   PetscCall(VecDestroy(&shell->right_add_work));
    448:   PetscCall(VecDestroy(&shell->axpy_left));
    449:   PetscCall(VecDestroy(&shell->axpy_right));
    450:   PetscCall(MatDestroy(&shell->axpy));
    451:   PetscCall(VecDestroy(&shell->zvals_w));
    452:   PetscCall(VecDestroy(&shell->zvals));
    453:   PetscCall(VecScatterDestroy(&shell->zvals_sct_c));
    454:   PetscCall(VecScatterDestroy(&shell->zvals_sct_r));
    455:   PetscCall(ISDestroy(&shell->zrows));
    456:   PetscCall(ISDestroy(&shell->zcols));

    458:   matmat = shell->matmat;
    459:   while (matmat) {
    460:     MatShellMatFunctionList next = matmat->next;

    462:     PetscCall(PetscObjectComposeFunction((PetscObject)mat, matmat->composedname, NULL));
    463:     PetscCall(PetscFree(matmat->composedname));
    464:     PetscCall(PetscFree(matmat->resultname));
    465:     PetscCall(PetscFree(matmat));
    466:     matmat = next;
    467:   }
    468:   PetscCall(MatShellSetContext(mat, NULL));
    469:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatShellGetContext_C", NULL));
    470:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatShellSetContext_C", NULL));
    471:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatShellSetContextDestroy_C", NULL));
    472:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatShellSetVecType_C", NULL));
    473:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatShellSetManageScalingShifts_C", NULL));
    474:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatShellGetScalingShifts_C", NULL));
    475:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatShellSetOperation_C", NULL));
    476:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatShellGetOperation_C", NULL));
    477:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatShellSetMatProductOperation_C", NULL));
    478:   PetscCall(PetscFree(mat->data));
    479:   PetscFunctionReturn(PETSC_SUCCESS);
    480: }

    482: typedef struct {
    483:   PetscErrorCode (*numeric)(Mat, Mat, Mat, void *);
    484:   PetscErrorCode (*destroy)(void *);
    485:   void *userdata;
    486:   Mat   B;
    487:   Mat   Bt;
    488:   Mat   axpy;
    489: } MatMatDataShell;

    491: static PetscErrorCode DestroyMatMatDataShell(void *data)
    492: {
    493:   MatMatDataShell *mmdata = (MatMatDataShell *)data;

    495:   PetscFunctionBegin;
    496:   if (mmdata->destroy) PetscCall((*mmdata->destroy)(mmdata->userdata));
    497:   PetscCall(MatDestroy(&mmdata->B));
    498:   PetscCall(MatDestroy(&mmdata->Bt));
    499:   PetscCall(MatDestroy(&mmdata->axpy));
    500:   PetscCall(PetscFree(mmdata));
    501:   PetscFunctionReturn(PETSC_SUCCESS);
    502: }

    504: static PetscErrorCode MatProductNumeric_Shell_X(Mat D)
    505: {
    506:   Mat_Product     *product;
    507:   Mat              A, B;
    508:   MatMatDataShell *mdata;
    509:   PetscScalar      zero = 0.0;

    511:   PetscFunctionBegin;
    512:   MatCheckProduct(D, 1);
    513:   product = D->product;
    514:   PetscCheck(product->data, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "Product data empty");
    515:   A     = product->A;
    516:   B     = product->B;
    517:   mdata = (MatMatDataShell *)product->data;
    518:   if (mdata->numeric) {
    519:     Mat_Shell *shell                = (Mat_Shell *)A->data;
    520:     PetscErrorCode (*stashsym)(Mat) = D->ops->productsymbolic;
    521:     PetscErrorCode (*stashnum)(Mat) = D->ops->productnumeric;
    522:     PetscBool useBmdata = PETSC_FALSE, newB = PETSC_TRUE;

    524:     if (shell->managescalingshifts) {
    525:       PetscCheck(!shell->zcols && !shell->zrows, PetscObjectComm((PetscObject)D), PETSC_ERR_SUP, "MatProduct not supported with zeroed rows/columns");
    526:       if (shell->right || shell->left) {
    527:         useBmdata = PETSC_TRUE;
    528:         if (!mdata->B) {
    529:           PetscCall(MatDuplicate(B, MAT_SHARE_NONZERO_PATTERN, &mdata->B));
    530:         } else {
    531:           newB = PETSC_FALSE;
    532:         }
    533:         PetscCall(MatCopy(B, mdata->B, SAME_NONZERO_PATTERN));
    534:       }
    535:       switch (product->type) {
    536:       case MATPRODUCT_AB: /* s L A R B + v L R B + L D R B */
    537:         if (shell->right) PetscCall(MatDiagonalScale(mdata->B, shell->right, NULL));
    538:         break;
    539:       case MATPRODUCT_AtB: /* s R A^t L B + v R L B + R D L B */
    540:         if (shell->left) PetscCall(MatDiagonalScale(mdata->B, shell->left, NULL));
    541:         break;
    542:       case MATPRODUCT_ABt: /* s L A R B^t + v L R B^t + L D R B^t */
    543:         if (shell->right) PetscCall(MatDiagonalScale(mdata->B, NULL, shell->right));
    544:         break;
    545:       case MATPRODUCT_RARt: /* s B L A R B^t + v B L R B^t + B L D R B^t */
    546:         if (shell->right && shell->left) {
    547:           PetscBool flg;

    549:           PetscCall(VecEqual(shell->right, shell->left, &flg));
    550:           PetscCheck(flg, PetscObjectComm((PetscObject)D), PETSC_ERR_SUP, "MatProductSymbolic type %s not supported for %s and %s matrices because left scaling != from right scaling", MatProductTypes[product->type], ((PetscObject)A)->type_name,
    551:                      ((PetscObject)B)->type_name);
    552:         }
    553:         if (shell->right) PetscCall(MatDiagonalScale(mdata->B, NULL, shell->right));
    554:         break;
    555:       case MATPRODUCT_PtAP: /* s B^t L A R B + v B^t L R B + B^t L D R B */
    556:         if (shell->right && shell->left) {
    557:           PetscBool flg;

    559:           PetscCall(VecEqual(shell->right, shell->left, &flg));
    560:           PetscCheck(flg, PetscObjectComm((PetscObject)D), PETSC_ERR_SUP, "MatProductSymbolic type %s not supported for %s and %s matrices because left scaling != from right scaling", MatProductTypes[product->type], ((PetscObject)A)->type_name,
    561:                      ((PetscObject)B)->type_name);
    562:         }
    563:         if (shell->right) PetscCall(MatDiagonalScale(mdata->B, shell->right, NULL));
    564:         break;
    565:       default:
    566:         SETERRQ(PetscObjectComm((PetscObject)D), PETSC_ERR_SUP, "MatProductSymbolic type %s not supported for %s and %s matrices", MatProductTypes[product->type], ((PetscObject)A)->type_name, ((PetscObject)B)->type_name);
    567:       }
    568:     }
    569:     /* allow the user to call MatMat operations on D */
    570:     D->product              = NULL;
    571:     D->ops->productsymbolic = NULL;
    572:     D->ops->productnumeric  = NULL;

    574:     PetscCall((*mdata->numeric)(A, useBmdata ? mdata->B : B, D, mdata->userdata));

    576:     /* clear any leftover user data and restore D pointers */
    577:     PetscCall(MatProductClear(D));
    578:     D->ops->productsymbolic = stashsym;
    579:     D->ops->productnumeric  = stashnum;
    580:     D->product              = product;

    582:     if (shell->managescalingshifts) {
    583:       PetscCall(MatScale(D, shell->vscale));
    584:       switch (product->type) {
    585:       case MATPRODUCT_AB:  /* s L A R B + v L R B + L D R B */
    586:       case MATPRODUCT_ABt: /* s L A R B^t + v L R B^t + L D R B^t */
    587:         if (shell->left) {
    588:           PetscCall(MatDiagonalScale(D, shell->left, NULL));
    589:           if (shell->dshift || shell->vshift != zero) {
    590:             if (!shell->left_work) PetscCall(MatCreateVecs(A, NULL, &shell->left_work));
    591:             if (shell->dshift) {
    592:               PetscCall(VecCopy(shell->dshift, shell->left_work));
    593:               PetscCall(VecShift(shell->left_work, shell->vshift));
    594:               PetscCall(VecPointwiseMult(shell->left_work, shell->left_work, shell->left));
    595:             } else {
    596:               PetscCall(VecSet(shell->left_work, shell->vshift));
    597:             }
    598:             if (product->type == MATPRODUCT_ABt) {
    599:               MatReuse     reuse = mdata->Bt ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX;
    600:               MatStructure str   = mdata->Bt ? SUBSET_NONZERO_PATTERN : DIFFERENT_NONZERO_PATTERN;

    602:               PetscCall(MatTranspose(mdata->B, reuse, &mdata->Bt));
    603:               PetscCall(MatDiagonalScale(mdata->Bt, shell->left_work, NULL));
    604:               PetscCall(MatAXPY(D, 1.0, mdata->Bt, str));
    605:             } else {
    606:               MatStructure str = newB ? DIFFERENT_NONZERO_PATTERN : SUBSET_NONZERO_PATTERN;

    608:               PetscCall(MatDiagonalScale(mdata->B, shell->left_work, NULL));
    609:               PetscCall(MatAXPY(D, 1.0, mdata->B, str));
    610:             }
    611:           }
    612:         }
    613:         break;
    614:       case MATPRODUCT_AtB: /* s R A^t L B + v R L B + R D L B */
    615:         if (shell->right) {
    616:           PetscCall(MatDiagonalScale(D, shell->right, NULL));
    617:           if (shell->dshift || shell->vshift != zero) {
    618:             MatStructure str = newB ? DIFFERENT_NONZERO_PATTERN : SUBSET_NONZERO_PATTERN;

    620:             if (!shell->right_work) PetscCall(MatCreateVecs(A, &shell->right_work, NULL));
    621:             if (shell->dshift) {
    622:               PetscCall(VecCopy(shell->dshift, shell->right_work));
    623:               PetscCall(VecShift(shell->right_work, shell->vshift));
    624:               PetscCall(VecPointwiseMult(shell->right_work, shell->right_work, shell->right));
    625:             } else {
    626:               PetscCall(VecSet(shell->right_work, shell->vshift));
    627:             }
    628:             PetscCall(MatDiagonalScale(mdata->B, shell->right_work, NULL));
    629:             PetscCall(MatAXPY(D, 1.0, mdata->B, str));
    630:           }
    631:         }
    632:         break;
    633:       case MATPRODUCT_PtAP: /* s B^t L A R B + v B^t L R B + B^t L D R B */
    634:       case MATPRODUCT_RARt: /* s B L A R B^t + v B L R B^t + B L D R B^t */
    635:         PetscCheck(!shell->dshift && shell->vshift == zero, PetscObjectComm((PetscObject)D), PETSC_ERR_SUP, "MatProductSymbolic type %s not supported for %s and %s matrices with diagonal shift", MatProductTypes[product->type], ((PetscObject)A)->type_name,
    636:                    ((PetscObject)B)->type_name);
    637:         break;
    638:       default:
    639:         SETERRQ(PetscObjectComm((PetscObject)D), PETSC_ERR_SUP, "MatProductSymbolic type %s not supported for %s and %s matrices", MatProductTypes[product->type], ((PetscObject)A)->type_name, ((PetscObject)B)->type_name);
    640:       }
    641:       if (shell->axpy && shell->axpy_vscale != zero) {
    642:         Mat              X;
    643:         PetscObjectState axpy_state;
    644:         MatStructure     str = DIFFERENT_NONZERO_PATTERN; /* not sure it is safe to ever use SUBSET_NONZERO_PATTERN */

    646:         PetscCall(MatShellGetContext(shell->axpy, &X));
    647:         PetscCall(PetscObjectStateGet((PetscObject)X, &axpy_state));
    648:         PetscCheck(shell->axpy_state == axpy_state, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Invalid AXPY state: cannot modify the X matrix passed to MatAXPY(Y,a,X,...)");
    649:         if (!mdata->axpy) {
    650:           str = DIFFERENT_NONZERO_PATTERN;
    651:           PetscCall(MatProductCreate(shell->axpy, B, NULL, &mdata->axpy));
    652:           PetscCall(MatProductSetType(mdata->axpy, product->type));
    653:           PetscCall(MatProductSetFromOptions(mdata->axpy));
    654:           PetscCall(MatProductSymbolic(mdata->axpy));
    655:         } else { /* May be that shell->axpy has changed */
    656:           PetscBool flg;

    658:           PetscCall(MatProductReplaceMats(shell->axpy, B, NULL, mdata->axpy));
    659:           PetscCall(MatHasOperation(mdata->axpy, MATOP_PRODUCTSYMBOLIC, &flg));
    660:           if (!flg) {
    661:             str = DIFFERENT_NONZERO_PATTERN;
    662:             PetscCall(MatProductSetFromOptions(mdata->axpy));
    663:             PetscCall(MatProductSymbolic(mdata->axpy));
    664:           }
    665:         }
    666:         PetscCall(MatProductNumeric(mdata->axpy));
    667:         PetscCall(MatAXPY(D, shell->axpy_vscale, mdata->axpy, str));
    668:       }
    669:     }
    670:   } else SETERRQ(PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "Missing numeric operation");
    671:   PetscFunctionReturn(PETSC_SUCCESS);
    672: }

    674: static PetscErrorCode MatProductSymbolic_Shell_X(Mat D)
    675: {
    676:   Mat_Product            *product;
    677:   Mat                     A, B;
    678:   MatShellMatFunctionList matmat;
    679:   Mat_Shell              *shell;
    680:   PetscBool               flg = PETSC_FALSE;
    681:   char                    composedname[256];
    682:   MatMatDataShell        *mdata;

    684:   PetscFunctionBegin;
    685:   MatCheckProduct(D, 1);
    686:   product = D->product;
    687:   PetscCheck(!product->data, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "Product data not empty");
    688:   A      = product->A;
    689:   B      = product->B;
    690:   shell  = (Mat_Shell *)A->data;
    691:   matmat = shell->matmat;
    692:   PetscCall(PetscSNPrintf(composedname, sizeof(composedname), "MatProductSetFromOptions_%s_%s_C", ((PetscObject)A)->type_name, ((PetscObject)B)->type_name));
    693:   while (matmat) {
    694:     PetscCall(PetscStrcmp(composedname, matmat->composedname, &flg));
    695:     flg = (PetscBool)(flg && (matmat->ptype == product->type));
    696:     if (flg) break;
    697:     matmat = matmat->next;
    698:   }
    699:   PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Composedname \"%s\" for product type %s not found", composedname, MatProductTypes[product->type]);
    700:   switch (product->type) {
    701:   case MATPRODUCT_AB:
    702:     PetscCall(MatSetSizes(D, A->rmap->n, B->cmap->n, A->rmap->N, B->cmap->N));
    703:     break;
    704:   case MATPRODUCT_AtB:
    705:     PetscCall(MatSetSizes(D, A->cmap->n, B->cmap->n, A->cmap->N, B->cmap->N));
    706:     break;
    707:   case MATPRODUCT_ABt:
    708:     PetscCall(MatSetSizes(D, A->rmap->n, B->rmap->n, A->rmap->N, B->rmap->N));
    709:     break;
    710:   case MATPRODUCT_RARt:
    711:     PetscCall(MatSetSizes(D, B->rmap->n, B->rmap->n, B->rmap->N, B->rmap->N));
    712:     break;
    713:   case MATPRODUCT_PtAP:
    714:     PetscCall(MatSetSizes(D, B->cmap->n, B->cmap->n, B->cmap->N, B->cmap->N));
    715:     break;
    716:   default:
    717:     SETERRQ(PetscObjectComm((PetscObject)D), PETSC_ERR_SUP, "MatProductSymbolic type %s not supported for %s and %s matrices", MatProductTypes[product->type], ((PetscObject)A)->type_name, ((PetscObject)B)->type_name);
    718:   }
    719:   /* respect users who passed in a matrix for which resultname is the base type */
    720:   if (matmat->resultname) {
    721:     PetscCall(PetscObjectBaseTypeCompare((PetscObject)D, matmat->resultname, &flg));
    722:     if (!flg) PetscCall(MatSetType(D, matmat->resultname));
    723:   }
    724:   /* If matrix type was not set or different, we need to reset this pointers */
    725:   D->ops->productsymbolic = MatProductSymbolic_Shell_X;
    726:   D->ops->productnumeric  = MatProductNumeric_Shell_X;
    727:   /* attach product data */
    728:   PetscCall(PetscNew(&mdata));
    729:   mdata->numeric = matmat->numeric;
    730:   mdata->destroy = matmat->destroy;
    731:   if (matmat->symbolic) {
    732:     PetscCall((*matmat->symbolic)(A, B, D, &mdata->userdata));
    733:   } else { /* call general setup if symbolic operation not provided */
    734:     PetscCall(MatSetUp(D));
    735:   }
    736:   PetscCheck(D->product, PetscObjectComm((PetscObject)D), PETSC_ERR_COR, "Product disappeared after user symbolic phase");
    737:   PetscCheck(!D->product->data, PetscObjectComm((PetscObject)D), PETSC_ERR_COR, "Product data not empty after user symbolic phase");
    738:   D->product->data    = mdata;
    739:   D->product->destroy = DestroyMatMatDataShell;
    740:   /* Be sure to reset these pointers if the user did something unexpected */
    741:   D->ops->productsymbolic = MatProductSymbolic_Shell_X;
    742:   D->ops->productnumeric  = MatProductNumeric_Shell_X;
    743:   PetscFunctionReturn(PETSC_SUCCESS);
    744: }

    746: static PetscErrorCode MatProductSetFromOptions_Shell_X(Mat D)
    747: {
    748:   Mat_Product            *product;
    749:   Mat                     A, B;
    750:   MatShellMatFunctionList matmat;
    751:   Mat_Shell              *shell;
    752:   PetscBool               flg;
    753:   char                    composedname[256];

    755:   PetscFunctionBegin;
    756:   MatCheckProduct(D, 1);
    757:   product = D->product;
    758:   A       = product->A;
    759:   B       = product->B;
    760:   PetscCall(MatIsShell(A, &flg));
    761:   if (!flg) PetscFunctionReturn(PETSC_SUCCESS);
    762:   shell  = (Mat_Shell *)A->data;
    763:   matmat = shell->matmat;
    764:   PetscCall(PetscSNPrintf(composedname, sizeof(composedname), "MatProductSetFromOptions_%s_%s_C", ((PetscObject)A)->type_name, ((PetscObject)B)->type_name));
    765:   while (matmat) {
    766:     PetscCall(PetscStrcmp(composedname, matmat->composedname, &flg));
    767:     flg = (PetscBool)(flg && (matmat->ptype == product->type));
    768:     if (flg) break;
    769:     matmat = matmat->next;
    770:   }
    771:   if (flg) {
    772:     D->ops->productsymbolic = MatProductSymbolic_Shell_X;
    773:   } else PetscCall(PetscInfo(D, "  symbolic product %s not registered for product type %s\n", composedname, MatProductTypes[product->type]));
    774:   PetscFunctionReturn(PETSC_SUCCESS);
    775: }

    777: static PetscErrorCode MatShellSetMatProductOperation_Private(Mat A, MatProductType ptype, PetscErrorCode (*symbolic)(Mat, Mat, Mat, void **), PetscErrorCode (*numeric)(Mat, Mat, Mat, void *), PetscErrorCode (*destroy)(void *), char *composedname, const char *resultname)
    778: {
    779:   PetscBool               flg;
    780:   Mat_Shell              *shell;
    781:   MatShellMatFunctionList matmat;

    783:   PetscFunctionBegin;
    784:   PetscCheck(numeric, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing numeric routine");
    785:   PetscCheck(composedname, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing composed name");

    787:   /* add product callback */
    788:   shell  = (Mat_Shell *)A->data;
    789:   matmat = shell->matmat;
    790:   if (!matmat) {
    791:     PetscCall(PetscNew(&shell->matmat));
    792:     matmat = shell->matmat;
    793:   } else {
    794:     MatShellMatFunctionList entry = matmat;
    795:     while (entry) {
    796:       PetscCall(PetscStrcmp(composedname, entry->composedname, &flg));
    797:       flg    = (PetscBool)(flg && (entry->ptype == ptype));
    798:       matmat = entry;
    799:       if (flg) goto set;
    800:       entry = entry->next;
    801:     }
    802:     PetscCall(PetscNew(&matmat->next));
    803:     matmat = matmat->next;
    804:   }

    806: set:
    807:   matmat->symbolic = symbolic;
    808:   matmat->numeric  = numeric;
    809:   matmat->destroy  = destroy;
    810:   matmat->ptype    = ptype;
    811:   PetscCall(PetscFree(matmat->composedname));
    812:   PetscCall(PetscFree(matmat->resultname));
    813:   PetscCall(PetscStrallocpy(composedname, &matmat->composedname));
    814:   PetscCall(PetscStrallocpy(resultname, &matmat->resultname));
    815:   PetscCall(PetscInfo(A, "Composing %s for product type %s with result %s\n", matmat->composedname, MatProductTypes[matmat->ptype], matmat->resultname ? matmat->resultname : "not specified"));
    816:   PetscCall(PetscObjectComposeFunction((PetscObject)A, matmat->composedname, MatProductSetFromOptions_Shell_X));
    817:   PetscFunctionReturn(PETSC_SUCCESS);
    818: }

    820: /*@C
    821:   MatShellSetMatProductOperation - Allows user to set a matrix matrix operation for a `MATSHELL` shell matrix.

    823:   Logically Collective; No Fortran Support

    825:   Input Parameters:
    826: + A        - the `MATSHELL` shell matrix
    827: . ptype    - the product type
    828: . symbolic - the function for the symbolic phase (can be `NULL`)
    829: . numeric  - the function for the numerical phase
    830: . destroy  - the function for the destruction of the needed data generated during the symbolic phase (can be `NULL`)
    831: . Btype    - the matrix type for the matrix to be multiplied against
    832: - Ctype    - the matrix type for the result (can be `NULL`)

    834:   Level: advanced

    836:   Example Usage:
    837: .vb
    838:   extern PetscErrorCode usersymbolic(Mat, Mat, Mat, void**);
    839:   extern PetscErrorCode usernumeric(Mat, Mat, Mat, void*);
    840:   extern PetscErrorCode userdestroy(void*);

    842:   MatCreateShell(comm, m, n, M, N, ctx, &A);
    843:   MatShellSetMatProductOperation(
    844:     A, MATPRODUCT_AB, usersymbolic, usernumeric, userdestroy,MATSEQAIJ, MATDENSE
    845:   );
    846:   // create B of type SEQAIJ etc..
    847:   MatProductCreate(A, B, PETSC_NULLPTR, &C);
    848:   MatProductSetType(C, MATPRODUCT_AB);
    849:   MatProductSetFromOptions(C);
    850:   MatProductSymbolic(C); // actually runs the user defined symbolic operation
    851:   MatProductNumeric(C); // actually runs the user defined numeric operation
    852:   // use C = A * B
    853: .ve

    855:   Notes:
    856:   `MATPRODUCT_ABC` is not supported yet.

    858:   If the symbolic phase is not specified, `MatSetUp()` is called on the result matrix that must have its type set if Ctype is `NULL`.

    860:   Any additional data needed by the matrix product needs to be returned during the symbolic phase and destroyed with the destroy callback.
    861:   PETSc will take care of calling the user-defined callbacks.
    862:   It is allowed to specify the same callbacks for different Btype matrix types.
    863:   The couple (Btype,ptype) uniquely identifies the operation, the last specified callbacks takes precedence.

    865: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellSetContext()`, `MatSetOperation()`, `MatProductType`, `MatType`, `MatSetUp()`
    866: @*/
    867: PetscErrorCode MatShellSetMatProductOperation(Mat A, MatProductType ptype, PetscErrorCode (*symbolic)(Mat, Mat, Mat, void **), PetscErrorCode (*numeric)(Mat, Mat, Mat, void *), PetscErrorCode (*destroy)(void *), MatType Btype, MatType Ctype)
    868: {
    869:   PetscFunctionBegin;
    872:   PetscCheck(ptype != MATPRODUCT_ABC, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not for product type %s", MatProductTypes[ptype]);
    873:   PetscCheck(numeric, PetscObjectComm((PetscObject)A), PETSC_ERR_USER, "Missing numeric routine, argument 4");
    874:   PetscAssertPointer(Btype, 6);
    875:   if (Ctype) PetscAssertPointer(Ctype, 7);
    876:   PetscTryMethod(A, "MatShellSetMatProductOperation_C", (Mat, MatProductType, PetscErrorCode (*)(Mat, Mat, Mat, void **), PetscErrorCode (*)(Mat, Mat, Mat, void *), PetscErrorCode (*)(void *), MatType, MatType), (A, ptype, symbolic, numeric, destroy, Btype, Ctype));
    877:   PetscFunctionReturn(PETSC_SUCCESS);
    878: }

    880: static PetscErrorCode MatShellSetMatProductOperation_Shell(Mat A, MatProductType ptype, PetscErrorCode (*symbolic)(Mat, Mat, Mat, void **), PetscErrorCode (*numeric)(Mat, Mat, Mat, void *), PetscErrorCode (*destroy)(void *), MatType Btype, MatType Ctype)
    881: {
    882:   PetscBool   flg;
    883:   char        composedname[256];
    884:   MatRootName Bnames = MatRootNameList, Cnames = MatRootNameList;
    885:   PetscMPIInt size;

    887:   PetscFunctionBegin;
    889:   while (Bnames) { /* user passed in the root name */
    890:     PetscCall(PetscStrcmp(Btype, Bnames->rname, &flg));
    891:     if (flg) break;
    892:     Bnames = Bnames->next;
    893:   }
    894:   while (Cnames) { /* user passed in the root name */
    895:     PetscCall(PetscStrcmp(Ctype, Cnames->rname, &flg));
    896:     if (flg) break;
    897:     Cnames = Cnames->next;
    898:   }
    899:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
    900:   Btype = Bnames ? (size > 1 ? Bnames->mname : Bnames->sname) : Btype;
    901:   Ctype = Cnames ? (size > 1 ? Cnames->mname : Cnames->sname) : Ctype;
    902:   PetscCall(PetscSNPrintf(composedname, sizeof(composedname), "MatProductSetFromOptions_%s_%s_C", ((PetscObject)A)->type_name, Btype));
    903:   PetscCall(MatShellSetMatProductOperation_Private(A, ptype, symbolic, numeric, destroy, composedname, Ctype));
    904:   PetscFunctionReturn(PETSC_SUCCESS);
    905: }

    907: static PetscErrorCode MatCopy_Shell(Mat A, Mat B, MatStructure str)
    908: {
    909:   Mat_Shell              *shellA = (Mat_Shell *)A->data, *shellB = (Mat_Shell *)B->data;
    910:   PetscBool               matflg;
    911:   MatShellMatFunctionList matmatA;

    913:   PetscFunctionBegin;
    914:   PetscCall(MatIsShell(B, &matflg));
    915:   PetscCheck(matflg, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Matrix %s not derived from MATSHELL", ((PetscObject)B)->type_name);

    917:   B->ops[0]      = A->ops[0];
    918:   shellB->ops[0] = shellA->ops[0];

    920:   if (shellA->ops->copy) PetscCall((*shellA->ops->copy)(A, B, str));
    921:   shellB->vscale = shellA->vscale;
    922:   shellB->vshift = shellA->vshift;
    923:   if (shellA->dshift) {
    924:     if (!shellB->dshift) PetscCall(VecDuplicate(shellA->dshift, &shellB->dshift));
    925:     PetscCall(VecCopy(shellA->dshift, shellB->dshift));
    926:   } else {
    927:     PetscCall(VecDestroy(&shellB->dshift));
    928:   }
    929:   if (shellA->left) {
    930:     if (!shellB->left) PetscCall(VecDuplicate(shellA->left, &shellB->left));
    931:     PetscCall(VecCopy(shellA->left, shellB->left));
    932:   } else {
    933:     PetscCall(VecDestroy(&shellB->left));
    934:   }
    935:   if (shellA->right) {
    936:     if (!shellB->right) PetscCall(VecDuplicate(shellA->right, &shellB->right));
    937:     PetscCall(VecCopy(shellA->right, shellB->right));
    938:   } else {
    939:     PetscCall(VecDestroy(&shellB->right));
    940:   }
    941:   PetscCall(MatDestroy(&shellB->axpy));
    942:   shellB->axpy_vscale = 0.0;
    943:   shellB->axpy_state  = 0;
    944:   if (shellA->axpy) {
    945:     PetscCall(PetscObjectReference((PetscObject)shellA->axpy));
    946:     shellB->axpy        = shellA->axpy;
    947:     shellB->axpy_vscale = shellA->axpy_vscale;
    948:     shellB->axpy_state  = shellA->axpy_state;
    949:   }
    950:   if (shellA->zrows) {
    951:     PetscCall(ISDuplicate(shellA->zrows, &shellB->zrows));
    952:     if (shellA->zcols) PetscCall(ISDuplicate(shellA->zcols, &shellB->zcols));
    953:     PetscCall(VecDuplicate(shellA->zvals, &shellB->zvals));
    954:     PetscCall(VecCopy(shellA->zvals, shellB->zvals));
    955:     PetscCall(VecDuplicate(shellA->zvals_w, &shellB->zvals_w));
    956:     PetscCall(PetscObjectReference((PetscObject)shellA->zvals_sct_r));
    957:     PetscCall(PetscObjectReference((PetscObject)shellA->zvals_sct_c));
    958:     shellB->zvals_sct_r = shellA->zvals_sct_r;
    959:     shellB->zvals_sct_c = shellA->zvals_sct_c;
    960:   }

    962:   matmatA = shellA->matmat;
    963:   if (matmatA) {
    964:     while (matmatA->next) {
    965:       PetscCall(MatShellSetMatProductOperation_Private(B, matmatA->ptype, matmatA->symbolic, matmatA->numeric, matmatA->destroy, matmatA->composedname, matmatA->resultname));
    966:       matmatA = matmatA->next;
    967:     }
    968:   }
    969:   PetscFunctionReturn(PETSC_SUCCESS);
    970: }

    972: static PetscErrorCode MatDuplicate_Shell(Mat mat, MatDuplicateOption op, Mat *M)
    973: {
    974:   PetscFunctionBegin;
    975:   PetscCall(MatCreateShell(PetscObjectComm((PetscObject)mat), mat->rmap->n, mat->cmap->n, mat->rmap->N, mat->cmap->N, NULL, M));
    976:   ((Mat_Shell *)(*M)->data)->ctxcontainer = ((Mat_Shell *)mat->data)->ctxcontainer;
    977:   PetscCall(PetscObjectCompose((PetscObject)*M, "MatShell ctx", (PetscObject)((Mat_Shell *)(*M)->data)->ctxcontainer));
    978:   PetscCall(PetscObjectChangeTypeName((PetscObject)*M, ((PetscObject)mat)->type_name));
    979:   if (op == MAT_COPY_VALUES) PetscCall(MatCopy(mat, *M, SAME_NONZERO_PATTERN));
    980:   PetscCall(PetscObjectCopyFortranFunctionPointers((PetscObject)mat, (PetscObject)*M));
    981:   PetscFunctionReturn(PETSC_SUCCESS);
    982: }

    984: static PetscErrorCode MatMult_Shell(Mat A, Vec x, Vec y)
    985: {
    986:   Mat_Shell       *shell = (Mat_Shell *)A->data;
    987:   Vec              xx;
    988:   PetscObjectState instate, outstate;

    990:   PetscFunctionBegin;
    991:   PetscCall(MatShellPreZeroRight(A, x, &xx));
    992:   PetscCall(MatShellPreScaleRight(A, xx, &xx));
    993:   PetscCall(PetscObjectStateGet((PetscObject)y, &instate));
    994:   PetscCall((*shell->ops->mult)(A, xx, y));
    995:   PetscCall(PetscObjectStateGet((PetscObject)y, &outstate));
    996:   if (instate == outstate) {
    997:     /* increase the state of the output vector since the user did not update its state themself as should have been done */
    998:     PetscCall(PetscObjectStateIncrease((PetscObject)y));
    999:   }
    1000:   PetscCall(MatShellShiftAndScale(A, xx, y, PETSC_FALSE));
    1001:   PetscCall(MatShellPostScaleLeft(A, y));
    1002:   PetscCall(MatShellPostZeroLeft(A, y));

    1004:   if (shell->axpy) {
    1005:     Mat              X;
    1006:     PetscObjectState axpy_state;

    1008:     PetscCall(MatShellGetContext(shell->axpy, &X));
    1009:     PetscCall(PetscObjectStateGet((PetscObject)X, &axpy_state));
    1010:     PetscCheck(shell->axpy_state == axpy_state, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Invalid AXPY state: cannot modify the X matrix passed to MatAXPY(Y,a,X,...)");

    1012:     PetscCall(MatCreateVecs(shell->axpy, shell->axpy_right ? NULL : &shell->axpy_right, shell->axpy_left ? NULL : &shell->axpy_left));
    1013:     PetscCall(VecCopy(x, shell->axpy_right));
    1014:     PetscCall(MatMult(shell->axpy, shell->axpy_right, shell->axpy_left));
    1015:     PetscCall(VecAXPY(y, shell->axpy_vscale, shell->axpy_left));
    1016:   }
    1017:   PetscFunctionReturn(PETSC_SUCCESS);
    1018: }

    1020: static PetscErrorCode MatMultAdd_Shell(Mat A, Vec x, Vec y, Vec z)
    1021: {
    1022:   Mat_Shell *shell = (Mat_Shell *)A->data;

    1024:   PetscFunctionBegin;
    1025:   if (y == z) {
    1026:     if (!shell->right_add_work) PetscCall(VecDuplicate(z, &shell->right_add_work));
    1027:     PetscCall(MatMult(A, x, shell->right_add_work));
    1028:     PetscCall(VecAXPY(z, 1.0, shell->right_add_work));
    1029:   } else {
    1030:     PetscCall(MatMult(A, x, z));
    1031:     PetscCall(VecAXPY(z, 1.0, y));
    1032:   }
    1033:   PetscFunctionReturn(PETSC_SUCCESS);
    1034: }

    1036: static PetscErrorCode MatMultTranspose_Shell(Mat A, Vec x, Vec y)
    1037: {
    1038:   Mat_Shell       *shell = (Mat_Shell *)A->data;
    1039:   Vec              xx;
    1040:   PetscObjectState instate, outstate;

    1042:   PetscFunctionBegin;
    1043:   PetscCall(MatShellPreZeroLeft(A, x, &xx));
    1044:   PetscCall(MatShellPreScaleLeft(A, xx, &xx, PETSC_FALSE));
    1045:   PetscCall(PetscObjectStateGet((PetscObject)y, &instate));
    1046:   PetscCall((*shell->ops->multtranspose)(A, xx, y));
    1047:   PetscCall(PetscObjectStateGet((PetscObject)y, &outstate));
    1048:   if (instate == outstate) {
    1049:     /* increase the state of the output vector since the user did not update its state themself as should have been done */
    1050:     PetscCall(PetscObjectStateIncrease((PetscObject)y));
    1051:   }
    1052:   PetscCall(MatShellShiftAndScale(A, xx, y, PETSC_FALSE));
    1053:   PetscCall(MatShellPostScaleRight(A, y, PETSC_FALSE));
    1054:   PetscCall(MatShellPostZeroRight(A, y));

    1056:   if (shell->axpy) {
    1057:     Mat              X;
    1058:     PetscObjectState axpy_state;

    1060:     PetscCall(MatShellGetContext(shell->axpy, &X));
    1061:     PetscCall(PetscObjectStateGet((PetscObject)X, &axpy_state));
    1062:     PetscCheck(shell->axpy_state == axpy_state, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Invalid AXPY state: cannot modify the X matrix passed to MatAXPY(Y,a,X,...)");
    1063:     PetscCall(MatCreateVecs(shell->axpy, shell->axpy_right ? NULL : &shell->axpy_right, shell->axpy_left ? NULL : &shell->axpy_left));
    1064:     PetscCall(VecCopy(x, shell->axpy_left));
    1065:     PetscCall(MatMultTranspose(shell->axpy, shell->axpy_left, shell->axpy_right));
    1066:     PetscCall(VecAXPY(y, shell->axpy_vscale, shell->axpy_right));
    1067:   }
    1068:   PetscFunctionReturn(PETSC_SUCCESS);
    1069: }

    1071: static PetscErrorCode MatMultHermitianTranspose_Shell(Mat A, Vec x, Vec y)
    1072: {
    1073:   Mat_Shell       *shell = (Mat_Shell *)A->data;
    1074:   Vec              xx;
    1075:   PetscObjectState instate, outstate;

    1077:   PetscFunctionBegin;
    1078:   PetscCall(MatShellPreZeroLeft(A, x, &xx));
    1079:   PetscCall(MatShellPreScaleLeft(A, xx, &xx, PETSC_TRUE));
    1080:   PetscCall(PetscObjectStateGet((PetscObject)y, &instate));
    1081:   PetscCall((*shell->ops->multhermitiantranspose)(A, xx, y));
    1082:   PetscCall(PetscObjectStateGet((PetscObject)y, &outstate));
    1083:   if (instate == outstate) {
    1084:     /* increase the state of the output vector since the user did not update its state themself as should have been done */
    1085:     PetscCall(PetscObjectStateIncrease((PetscObject)y));
    1086:   }
    1087:   PetscCall(MatShellShiftAndScale(A, xx, y, PETSC_TRUE));
    1088:   PetscCall(MatShellPostScaleRight(A, y, PETSC_TRUE));
    1089:   PetscCall(MatShellPostZeroRight(A, y));

    1091:   if (shell->axpy) {
    1092:     Mat              X;
    1093:     PetscObjectState axpy_state;

    1095:     PetscCall(MatShellGetContext(shell->axpy, &X));
    1096:     PetscCall(PetscObjectStateGet((PetscObject)X, &axpy_state));
    1097:     PetscCheck(shell->axpy_state == axpy_state, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Invalid AXPY state: cannot modify the X matrix passed to MatAXPY(Y,a,X,...)");
    1098:     PetscCall(MatCreateVecs(shell->axpy, shell->axpy_right ? NULL : &shell->axpy_right, shell->axpy_left ? NULL : &shell->axpy_left));
    1099:     PetscCall(VecCopy(x, shell->axpy_left));
    1100:     PetscCall(MatMultHermitianTranspose(shell->axpy, shell->axpy_left, shell->axpy_right));
    1101:     PetscCall(VecAXPY(y, PetscConj(shell->axpy_vscale), shell->axpy_right));
    1102:   }
    1103:   PetscFunctionReturn(PETSC_SUCCESS);
    1104: }

    1106: static PetscErrorCode MatMultTransposeAdd_Shell(Mat A, Vec x, Vec y, Vec z)
    1107: {
    1108:   Mat_Shell *shell = (Mat_Shell *)A->data;

    1110:   PetscFunctionBegin;
    1111:   if (y == z) {
    1112:     if (!shell->left_add_work) PetscCall(VecDuplicate(z, &shell->left_add_work));
    1113:     PetscCall(MatMultTranspose(A, x, shell->left_add_work));
    1114:     PetscCall(VecAXPY(z, 1.0, shell->left_add_work));
    1115:   } else {
    1116:     PetscCall(MatMultTranspose(A, x, z));
    1117:     PetscCall(VecAXPY(z, 1.0, y));
    1118:   }
    1119:   PetscFunctionReturn(PETSC_SUCCESS);
    1120: }

    1122: static PetscErrorCode MatMultHermitianTransposeAdd_Shell(Mat A, Vec x, Vec y, Vec z)
    1123: {
    1124:   Mat_Shell *shell = (Mat_Shell *)A->data;

    1126:   PetscFunctionBegin;
    1127:   if (y == z) {
    1128:     if (!shell->left_add_work) PetscCall(VecDuplicate(z, &shell->left_add_work));
    1129:     PetscCall(MatMultHermitianTranspose(A, x, shell->left_add_work));
    1130:     PetscCall(VecAXPY(z, 1.0, shell->left_add_work));
    1131:   } else {
    1132:     PetscCall(MatMultHermitianTranspose(A, x, z));
    1133:     PetscCall(VecAXPY(z, 1.0, y));
    1134:   }
    1135:   PetscFunctionReturn(PETSC_SUCCESS);
    1136: }

    1138: /*
    1139:           diag(left)(vscale*A + diag(dshift) + vshift I)diag(right)
    1140: */
    1141: static PetscErrorCode MatGetDiagonal_Shell(Mat A, Vec v)
    1142: {
    1143:   Mat_Shell *shell = (Mat_Shell *)A->data;

    1145:   PetscFunctionBegin;
    1146:   if (shell->ops->getdiagonal) {
    1147:     PetscCall((*shell->ops->getdiagonal)(A, v));
    1148:   } else SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Must provide shell matrix with routine to return diagonal using\nMatShellSetOperation(S,MATOP_GET_DIAGONAL,...)");
    1149:   PetscCall(VecScale(v, shell->vscale));
    1150:   if (shell->dshift) PetscCall(VecAXPY(v, 1.0, shell->dshift));
    1151:   PetscCall(VecShift(v, shell->vshift));
    1152:   if (shell->left) PetscCall(VecPointwiseMult(v, v, shell->left));
    1153:   if (shell->right) PetscCall(VecPointwiseMult(v, v, shell->right));
    1154:   if (shell->zrows) {
    1155:     PetscCall(VecScatterBegin(shell->zvals_sct_r, shell->zvals, v, INSERT_VALUES, SCATTER_REVERSE));
    1156:     PetscCall(VecScatterEnd(shell->zvals_sct_r, shell->zvals, v, INSERT_VALUES, SCATTER_REVERSE));
    1157:   }
    1158:   if (shell->axpy) {
    1159:     Mat              X;
    1160:     PetscObjectState axpy_state;

    1162:     PetscCall(MatShellGetContext(shell->axpy, &X));
    1163:     PetscCall(PetscObjectStateGet((PetscObject)X, &axpy_state));
    1164:     PetscCheck(shell->axpy_state == axpy_state, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Invalid AXPY state: cannot modify the X matrix passed to MatAXPY(Y,a,X,...)");
    1165:     PetscCall(MatCreateVecs(shell->axpy, NULL, shell->axpy_left ? NULL : &shell->axpy_left));
    1166:     PetscCall(MatGetDiagonal(shell->axpy, shell->axpy_left));
    1167:     PetscCall(VecAXPY(v, shell->axpy_vscale, shell->axpy_left));
    1168:   }
    1169:   PetscFunctionReturn(PETSC_SUCCESS);
    1170: }

    1172: static PetscErrorCode MatGetDiagonalBlock_Shell(Mat A, Mat *b)
    1173: {
    1174:   Mat_Shell *shell = (Mat_Shell *)A->data;
    1175:   Vec        left = NULL, right = NULL;

    1177:   PetscFunctionBegin;
    1178:   PetscCheck(!shell->zrows && !shell->zcols, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Cannot call MatGetDiagonalBlock() if MatZeroRows() or MatZeroRowsColumns() has been called on the input Mat"); // TODO FIXME
    1179:   PetscCheck(!shell->axpy, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Cannot call MatGetDiagonalBlock() if MatAXPY() has been called on the input Mat");                                               // TODO FIXME
    1180:   PetscCheck(!shell->dshift, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Cannot call MatGetDiagonalBlock() if MatDiagonalSet() has been called on the input Mat");                                      // TODO FIXME
    1181:   if (shell->ops->getdiagonalblock) {
    1182:     PetscCall((*shell->ops->getdiagonalblock)(A, b));
    1183:   } else SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Must provide shell matrix with routine to return diagonal block using\nMatShellSetOperation(S,MATOP_GET_DIAGONAL_BLOCK,...)");
    1184:   PetscCall(MatScale(*b, shell->vscale));
    1185:   PetscCall(MatShift(*b, shell->vshift));
    1186:   if (shell->left) {
    1187:     PetscCall(VecCreateLocalVector(shell->left, &left));
    1188:     PetscCall(VecGetLocalVectorRead(shell->left, left));
    1189:   }
    1190:   if (shell->right) {
    1191:     PetscCall(VecCreateLocalVector(shell->right, &right));
    1192:     PetscCall(VecGetLocalVectorRead(shell->right, right));
    1193:   }
    1194:   PetscCall(MatDiagonalScale(*b, left, right));
    1195:   if (shell->left) {
    1196:     PetscCall(VecRestoreLocalVectorRead(shell->left, left));
    1197:     PetscCall(VecDestroy(&left));
    1198:   }
    1199:   if (shell->right) {
    1200:     PetscCall(VecRestoreLocalVectorRead(shell->right, right));
    1201:     PetscCall(VecDestroy(&right));
    1202:   }
    1203:   PetscFunctionReturn(PETSC_SUCCESS);
    1204: }

    1206: static PetscErrorCode MatShift_Shell(Mat Y, PetscScalar a)
    1207: {
    1208:   Mat_Shell *shell = (Mat_Shell *)Y->data;
    1209:   PetscBool  flg;

    1211:   PetscFunctionBegin;
    1212:   PetscCall(MatHasCongruentLayouts(Y, &flg));
    1213:   PetscCheck(flg, PetscObjectComm((PetscObject)Y), PETSC_ERR_SUP, "Cannot shift shell matrix if it is not congruent");
    1214:   if (shell->left || shell->right) {
    1215:     if (!shell->dshift) {
    1216:       PetscCall(VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift));
    1217:       PetscCall(VecSet(shell->dshift, a));
    1218:     } else {
    1219:       if (shell->left) PetscCall(VecPointwiseMult(shell->dshift, shell->dshift, shell->left));
    1220:       if (shell->right) PetscCall(VecPointwiseMult(shell->dshift, shell->dshift, shell->right));
    1221:       PetscCall(VecShift(shell->dshift, a));
    1222:     }
    1223:     if (shell->left) PetscCall(VecPointwiseDivide(shell->dshift, shell->dshift, shell->left));
    1224:     if (shell->right) PetscCall(VecPointwiseDivide(shell->dshift, shell->dshift, shell->right));
    1225:   } else shell->vshift += a;
    1226:   if (shell->zrows) PetscCall(VecShift(shell->zvals, a));
    1227:   PetscFunctionReturn(PETSC_SUCCESS);
    1228: }

    1230: static PetscErrorCode MatDiagonalSet_Shell_Private(Mat A, Vec D, PetscScalar s)
    1231: {
    1232:   Mat_Shell *shell = (Mat_Shell *)A->data;

    1234:   PetscFunctionBegin;
    1235:   if (!shell->dshift) PetscCall(VecDuplicate(D, &shell->dshift));
    1236:   if (shell->left || shell->right) {
    1237:     if (!shell->right_work) PetscCall(VecDuplicate(shell->left ? shell->left : shell->right, &shell->right_work));
    1238:     if (shell->left && shell->right) {
    1239:       PetscCall(VecPointwiseDivide(shell->right_work, D, shell->left));
    1240:       PetscCall(VecPointwiseDivide(shell->right_work, shell->right_work, shell->right));
    1241:     } else if (shell->left) {
    1242:       PetscCall(VecPointwiseDivide(shell->right_work, D, shell->left));
    1243:     } else {
    1244:       PetscCall(VecPointwiseDivide(shell->right_work, D, shell->right));
    1245:     }
    1246:     PetscCall(VecAXPY(shell->dshift, s, shell->right_work));
    1247:   } else {
    1248:     PetscCall(VecAXPY(shell->dshift, s, D));
    1249:   }
    1250:   PetscFunctionReturn(PETSC_SUCCESS);
    1251: }

    1253: static PetscErrorCode MatDiagonalSet_Shell(Mat A, Vec D, InsertMode ins)
    1254: {
    1255:   Mat_Shell *shell = (Mat_Shell *)A->data;
    1256:   Vec        d;
    1257:   PetscBool  flg;

    1259:   PetscFunctionBegin;
    1260:   PetscCall(MatHasCongruentLayouts(A, &flg));
    1261:   PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Cannot diagonal set or shift shell matrix if it is not congruent");
    1262:   if (ins == INSERT_VALUES) {
    1263:     PetscCall(VecDuplicate(D, &d));
    1264:     PetscCall(MatGetDiagonal(A, d));
    1265:     PetscCall(MatDiagonalSet_Shell_Private(A, d, -1.));
    1266:     PetscCall(MatDiagonalSet_Shell_Private(A, D, 1.));
    1267:     PetscCall(VecDestroy(&d));
    1268:     if (shell->zrows) PetscCall(VecCopy(D, shell->zvals));
    1269:   } else {
    1270:     PetscCall(MatDiagonalSet_Shell_Private(A, D, 1.));
    1271:     if (shell->zrows) PetscCall(VecAXPY(shell->zvals, 1.0, D));
    1272:   }
    1273:   PetscFunctionReturn(PETSC_SUCCESS);
    1274: }

    1276: static PetscErrorCode MatScale_Shell(Mat Y, PetscScalar a)
    1277: {
    1278:   Mat_Shell *shell = (Mat_Shell *)Y->data;

    1280:   PetscFunctionBegin;
    1281:   shell->vscale *= a;
    1282:   shell->vshift *= a;
    1283:   if (shell->dshift) PetscCall(VecScale(shell->dshift, a));
    1284:   shell->axpy_vscale *= a;
    1285:   if (shell->zrows) PetscCall(VecScale(shell->zvals, a));
    1286:   PetscFunctionReturn(PETSC_SUCCESS);
    1287: }

    1289: static PetscErrorCode MatDiagonalScale_Shell(Mat Y, Vec left, Vec right)
    1290: {
    1291:   Mat_Shell *shell = (Mat_Shell *)Y->data;

    1293:   PetscFunctionBegin;
    1294:   if (left) {
    1295:     if (!shell->left) {
    1296:       PetscCall(VecDuplicate(left, &shell->left));
    1297:       PetscCall(VecCopy(left, shell->left));
    1298:     } else {
    1299:       PetscCall(VecPointwiseMult(shell->left, shell->left, left));
    1300:     }
    1301:     if (shell->zrows) PetscCall(VecPointwiseMult(shell->zvals, shell->zvals, left));
    1302:   }
    1303:   if (right) {
    1304:     if (!shell->right) {
    1305:       PetscCall(VecDuplicate(right, &shell->right));
    1306:       PetscCall(VecCopy(right, shell->right));
    1307:     } else {
    1308:       PetscCall(VecPointwiseMult(shell->right, shell->right, right));
    1309:     }
    1310:     if (shell->zrows) {
    1311:       if (!shell->left_work) PetscCall(MatCreateVecs(Y, NULL, &shell->left_work));
    1312:       PetscCall(VecSet(shell->zvals_w, 1.0));
    1313:       PetscCall(VecScatterBegin(shell->zvals_sct_c, right, shell->zvals_w, INSERT_VALUES, SCATTER_FORWARD));
    1314:       PetscCall(VecScatterEnd(shell->zvals_sct_c, right, shell->zvals_w, INSERT_VALUES, SCATTER_FORWARD));
    1315:       PetscCall(VecPointwiseMult(shell->zvals, shell->zvals, shell->zvals_w));
    1316:     }
    1317:   }
    1318:   if (shell->axpy) PetscCall(MatDiagonalScale(shell->axpy, left, right));
    1319:   PetscFunctionReturn(PETSC_SUCCESS);
    1320: }

    1322: PETSC_INTERN PetscErrorCode MatAssemblyEnd_Shell(Mat Y, MatAssemblyType t)
    1323: {
    1324:   Mat_Shell *shell = (Mat_Shell *)Y->data;

    1326:   PetscFunctionBegin;
    1327:   if (t == MAT_FINAL_ASSEMBLY) {
    1328:     shell->vshift      = 0.0;
    1329:     shell->vscale      = 1.0;
    1330:     shell->axpy_vscale = 0.0;
    1331:     shell->axpy_state  = 0;
    1332:     PetscCall(VecDestroy(&shell->dshift));
    1333:     PetscCall(VecDestroy(&shell->left));
    1334:     PetscCall(VecDestroy(&shell->right));
    1335:     PetscCall(MatDestroy(&shell->axpy));
    1336:     PetscCall(VecDestroy(&shell->axpy_left));
    1337:     PetscCall(VecDestroy(&shell->axpy_right));
    1338:     PetscCall(VecScatterDestroy(&shell->zvals_sct_c));
    1339:     PetscCall(VecScatterDestroy(&shell->zvals_sct_r));
    1340:     PetscCall(ISDestroy(&shell->zrows));
    1341:     PetscCall(ISDestroy(&shell->zcols));
    1342:   }
    1343:   PetscFunctionReturn(PETSC_SUCCESS);
    1344: }

    1346: static PetscErrorCode MatMissingDiagonal_Shell(Mat A, PetscBool *missing, PetscInt *d)
    1347: {
    1348:   PetscFunctionBegin;
    1349:   *missing = PETSC_FALSE;
    1350:   PetscFunctionReturn(PETSC_SUCCESS);
    1351: }

    1353: static PetscErrorCode MatAXPY_Shell(Mat Y, PetscScalar a, Mat X, MatStructure str)
    1354: {
    1355:   Mat_Shell *shell = (Mat_Shell *)Y->data;

    1357:   PetscFunctionBegin;
    1358:   if (X == Y) {
    1359:     PetscCall(MatScale(Y, 1.0 + a));
    1360:     PetscFunctionReturn(PETSC_SUCCESS);
    1361:   }
    1362:   if (!shell->axpy) {
    1363:     PetscCall(MatConvertFrom_Shell(X, MATSHELL, MAT_INITIAL_MATRIX, &shell->axpy));
    1364:     shell->axpy_vscale = a;
    1365:     PetscCall(PetscObjectStateGet((PetscObject)X, &shell->axpy_state));
    1366:   } else {
    1367:     PetscCall(MatAXPY(shell->axpy, a / shell->axpy_vscale, X, str));
    1368:   }
    1369:   PetscFunctionReturn(PETSC_SUCCESS);
    1370: }

    1372: static struct _MatOps MatOps_Values = {NULL,
    1373:                                        NULL,
    1374:                                        NULL,
    1375:                                        NULL,
    1376:                                        /* 4*/ MatMultAdd_Shell,
    1377:                                        NULL,
    1378:                                        MatMultTransposeAdd_Shell,
    1379:                                        NULL,
    1380:                                        NULL,
    1381:                                        NULL,
    1382:                                        /*10*/ NULL,
    1383:                                        NULL,
    1384:                                        NULL,
    1385:                                        NULL,
    1386:                                        NULL,
    1387:                                        /*15*/ NULL,
    1388:                                        NULL,
    1389:                                        NULL,
    1390:                                        MatDiagonalScale_Shell,
    1391:                                        NULL,
    1392:                                        /*20*/ NULL,
    1393:                                        MatAssemblyEnd_Shell,
    1394:                                        NULL,
    1395:                                        NULL,
    1396:                                        /*24*/ MatZeroRows_Shell,
    1397:                                        NULL,
    1398:                                        NULL,
    1399:                                        NULL,
    1400:                                        NULL,
    1401:                                        /*29*/ NULL,
    1402:                                        NULL,
    1403:                                        NULL,
    1404:                                        /*32*/ NULL,
    1405:                                        NULL,
    1406:                                        /*34*/ MatDuplicate_Shell,
    1407:                                        NULL,
    1408:                                        NULL,
    1409:                                        NULL,
    1410:                                        NULL,
    1411:                                        /*39*/ MatAXPY_Shell,
    1412:                                        NULL,
    1413:                                        NULL,
    1414:                                        NULL,
    1415:                                        MatCopy_Shell,
    1416:                                        /*44*/ NULL,
    1417:                                        MatScale_Shell,
    1418:                                        MatShift_Shell,
    1419:                                        MatDiagonalSet_Shell,
    1420:                                        MatZeroRowsColumns_Shell,
    1421:                                        /*49*/ NULL,
    1422:                                        NULL,
    1423:                                        NULL,
    1424:                                        NULL,
    1425:                                        NULL,
    1426:                                        /*54*/ NULL,
    1427:                                        NULL,
    1428:                                        NULL,
    1429:                                        NULL,
    1430:                                        NULL,
    1431:                                        /*59*/ NULL,
    1432:                                        MatDestroy_Shell,
    1433:                                        NULL,
    1434:                                        MatConvertFrom_Shell,
    1435:                                        NULL,
    1436:                                        /*64*/ NULL,
    1437:                                        NULL,
    1438:                                        NULL,
    1439:                                        NULL,
    1440:                                        NULL,
    1441:                                        /*69*/ NULL,
    1442:                                        NULL,
    1443:                                        MatConvert_Shell,
    1444:                                        NULL,
    1445:                                        NULL,
    1446:                                        /*74*/ NULL,
    1447:                                        NULL,
    1448:                                        NULL,
    1449:                                        NULL,
    1450:                                        NULL,
    1451:                                        /*79*/ NULL,
    1452:                                        NULL,
    1453:                                        NULL,
    1454:                                        NULL,
    1455:                                        NULL,
    1456:                                        /*84*/ NULL,
    1457:                                        NULL,
    1458:                                        NULL,
    1459:                                        NULL,
    1460:                                        NULL,
    1461:                                        /*89*/ NULL,
    1462:                                        NULL,
    1463:                                        NULL,
    1464:                                        NULL,
    1465:                                        NULL,
    1466:                                        /*94*/ NULL,
    1467:                                        NULL,
    1468:                                        NULL,
    1469:                                        NULL,
    1470:                                        NULL,
    1471:                                        /*99*/ NULL,
    1472:                                        NULL,
    1473:                                        NULL,
    1474:                                        NULL,
    1475:                                        NULL,
    1476:                                        /*104*/ NULL,
    1477:                                        NULL,
    1478:                                        NULL,
    1479:                                        NULL,
    1480:                                        NULL,
    1481:                                        /*109*/ NULL,
    1482:                                        NULL,
    1483:                                        NULL,
    1484:                                        NULL,
    1485:                                        MatMissingDiagonal_Shell,
    1486:                                        /*114*/ NULL,
    1487:                                        NULL,
    1488:                                        NULL,
    1489:                                        NULL,
    1490:                                        NULL,
    1491:                                        /*119*/ NULL,
    1492:                                        NULL,
    1493:                                        NULL,
    1494:                                        MatMultHermitianTransposeAdd_Shell,
    1495:                                        NULL,
    1496:                                        /*124*/ NULL,
    1497:                                        NULL,
    1498:                                        NULL,
    1499:                                        NULL,
    1500:                                        NULL,
    1501:                                        /*129*/ NULL,
    1502:                                        NULL,
    1503:                                        NULL,
    1504:                                        NULL,
    1505:                                        NULL,
    1506:                                        /*134*/ NULL,
    1507:                                        NULL,
    1508:                                        NULL,
    1509:                                        NULL,
    1510:                                        NULL,
    1511:                                        /*139*/ NULL,
    1512:                                        NULL,
    1513:                                        NULL,
    1514:                                        NULL,
    1515:                                        NULL,
    1516:                                        /*144*/ NULL,
    1517:                                        NULL,
    1518:                                        NULL,
    1519:                                        NULL,
    1520:                                        NULL,
    1521:                                        NULL,
    1522:                                        /*150*/ NULL,
    1523:                                        NULL,
    1524:                                        NULL,
    1525:                                        NULL,
    1526:                                        NULL,
    1527:                                        NULL};

    1529: static PetscErrorCode MatShellSetContext_Shell(Mat mat, void *ctx)
    1530: {
    1531:   Mat_Shell *shell = (Mat_Shell *)mat->data;

    1533:   PetscFunctionBegin;
    1534:   if (ctx) {
    1535:     PetscContainer ctxcontainer;
    1536:     PetscCall(PetscContainerCreate(PetscObjectComm((PetscObject)mat), &ctxcontainer));
    1537:     PetscCall(PetscContainerSetPointer(ctxcontainer, ctx));
    1538:     PetscCall(PetscObjectCompose((PetscObject)mat, "MatShell ctx", (PetscObject)ctxcontainer));
    1539:     shell->ctxcontainer = ctxcontainer;
    1540:     PetscCall(PetscContainerDestroy(&ctxcontainer));
    1541:   } else {
    1542:     PetscCall(PetscObjectCompose((PetscObject)mat, "MatShell ctx", NULL));
    1543:     shell->ctxcontainer = NULL;
    1544:   }
    1545:   PetscFunctionReturn(PETSC_SUCCESS);
    1546: }

    1548: static PetscErrorCode MatShellSetContextDestroy_Shell(Mat mat, PetscErrorCode (*f)(void *))
    1549: {
    1550:   Mat_Shell *shell = (Mat_Shell *)mat->data;

    1552:   PetscFunctionBegin;
    1553:   if (shell->ctxcontainer) PetscCall(PetscContainerSetUserDestroy(shell->ctxcontainer, f));
    1554:   PetscFunctionReturn(PETSC_SUCCESS);
    1555: }

    1557: PetscErrorCode MatShellSetContext_Immutable(Mat mat, void *ctx)
    1558: {
    1559:   PetscFunctionBegin;
    1560:   SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "Cannot call MatShellSetContext() for a %s, it is used internally by the structure", ((PetscObject)mat)->type_name);
    1561:   PetscFunctionReturn(PETSC_SUCCESS);
    1562: }

    1564: PetscErrorCode MatShellSetContextDestroy_Immutable(Mat mat, PetscErrorCode (*f)(void *))
    1565: {
    1566:   PetscFunctionBegin;
    1567:   SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "Cannot call MatShellSetContextDestroy() for a %s, it is used internally by the structure", ((PetscObject)mat)->type_name);
    1568:   PetscFunctionReturn(PETSC_SUCCESS);
    1569: }

    1571: PetscErrorCode MatShellSetManageScalingShifts_Immutable(Mat mat)
    1572: {
    1573:   PetscFunctionBegin;
    1574:   SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "Cannot call MatShellSetManageScalingShifts() for a %s, it is used internally by the structure", ((PetscObject)mat)->type_name);
    1575:   PetscFunctionReturn(PETSC_SUCCESS);
    1576: }

    1578: static PetscErrorCode MatShellSetVecType_Shell(Mat mat, VecType vtype)
    1579: {
    1580:   PetscFunctionBegin;
    1581:   PetscCall(PetscFree(mat->defaultvectype));
    1582:   PetscCall(PetscStrallocpy(vtype, (char **)&mat->defaultvectype));
    1583:   PetscFunctionReturn(PETSC_SUCCESS);
    1584: }

    1586: static PetscErrorCode MatShellSetManageScalingShifts_Shell(Mat A)
    1587: {
    1588:   Mat_Shell *shell = (Mat_Shell *)A->data;

    1590:   PetscFunctionBegin;
    1591:   shell->managescalingshifts = PETSC_FALSE;
    1592:   A->ops->diagonalset        = NULL;
    1593:   A->ops->diagonalscale      = NULL;
    1594:   A->ops->scale              = NULL;
    1595:   A->ops->shift              = NULL;
    1596:   A->ops->axpy               = NULL;
    1597:   PetscFunctionReturn(PETSC_SUCCESS);
    1598: }

    1600: static PetscErrorCode MatShellGetScalingShifts_Shell(Mat A, PetscScalar *vshift, PetscScalar *vscale, Vec *dshift, Vec *left, Vec *right, Mat *axpy, IS *zrows, IS *zcols)
    1601: {
    1602:   Mat_Shell *shell = (Mat_Shell *)A->data;

    1604:   PetscFunctionBegin;
    1605:   PetscCheck(shell->managescalingshifts, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Cannot get scaling and shifts if MatShellSetManageScalingShifts() has been called");
    1606:   if (vshift == MAT_SHELL_NOT_ALLOWED) PetscCheck(shell->vshift == 0.0, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Non-trivial member in parent MatShell: vshift != 0.0, set via MatShift()");
    1607:   else if (vshift) *vshift = shell->vshift;
    1608:   if (vscale == MAT_SHELL_NOT_ALLOWED) PetscCheck(shell->vscale == 1.0, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Non-trivial member in parent MatShell: vscale != 1.0, set via MatScale()");
    1609:   else if (vscale) *vscale = shell->vscale;
    1610:   if (dshift == MAT_SHELL_NOT_ALLOWED) PetscCheck(!shell->dshift, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Non-trivial member in parent MatShell: dshift, set via MatDiagonalSet()");
    1611:   else if (dshift) *dshift = shell->dshift;
    1612:   if (left == MAT_SHELL_NOT_ALLOWED) PetscCheck(!shell->left, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Non-trivial member in parent MatShell: left, set via MatDiagonalScale()");
    1613:   else if (left) *left = shell->left;
    1614:   if (right == MAT_SHELL_NOT_ALLOWED) PetscCheck(!shell->right, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Non-trivial member in parent MatShell: right, set via MatDiagonalScale()");
    1615:   else if (right) *right = shell->right;
    1616:   if (axpy == MAT_SHELL_NOT_ALLOWED) PetscCheck(!shell->axpy, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Non-trivial member in parent MatShell: axpy, set via MatAXPY()");
    1617:   else if (axpy) *axpy = shell->axpy;
    1618:   if (zrows == MAT_SHELL_NOT_ALLOWED) PetscCheck(!shell->zrows, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Non-trivial member in parent MatShell: zrows, set via MatZeroRows()");
    1619:   else if (zrows) *zrows = shell->zrows;
    1620:   if (zcols == MAT_SHELL_NOT_ALLOWED) PetscCheck(!shell->zcols, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Non-trivial member in parent MatShell: zcols, set via MatZeroRowsColumns()");
    1621:   else if (zcols) *zcols = shell->zcols;
    1622:   PetscFunctionReturn(PETSC_SUCCESS);
    1623: }

    1625: static PetscErrorCode MatShellSetOperation_Shell(Mat mat, MatOperation op, void (*f)(void))
    1626: {
    1627:   Mat_Shell *shell = (Mat_Shell *)mat->data;

    1629:   PetscFunctionBegin;
    1630:   switch (op) {
    1631:   case MATOP_DESTROY:
    1632:     shell->ops->destroy = (PetscErrorCode (*)(Mat))f;
    1633:     break;
    1634:   case MATOP_VIEW:
    1635:     if (!mat->ops->viewnative) mat->ops->viewnative = mat->ops->view;
    1636:     mat->ops->view = (PetscErrorCode (*)(Mat, PetscViewer))f;
    1637:     break;
    1638:   case MATOP_COPY:
    1639:     shell->ops->copy = (PetscErrorCode (*)(Mat, Mat, MatStructure))f;
    1640:     break;
    1641:   case MATOP_DIAGONAL_SET:
    1642:   case MATOP_DIAGONAL_SCALE:
    1643:   case MATOP_SHIFT:
    1644:   case MATOP_SCALE:
    1645:   case MATOP_AXPY:
    1646:   case MATOP_ZERO_ROWS:
    1647:   case MATOP_ZERO_ROWS_COLUMNS:
    1648:     PetscCheck(!shell->managescalingshifts, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "MATSHELL is managing scalings and shifts, see MatShellSetManageScalingShifts()");
    1649:     (((void (**)(void))mat->ops)[op]) = f;
    1650:     break;
    1651:   case MATOP_GET_DIAGONAL:
    1652:     if (shell->managescalingshifts) {
    1653:       shell->ops->getdiagonal = (PetscErrorCode (*)(Mat, Vec))f;
    1654:       mat->ops->getdiagonal   = MatGetDiagonal_Shell;
    1655:     } else {
    1656:       shell->ops->getdiagonal = NULL;
    1657:       mat->ops->getdiagonal   = (PetscErrorCode (*)(Mat, Vec))f;
    1658:     }
    1659:     break;
    1660:   case MATOP_GET_DIAGONAL_BLOCK:
    1661:     if (shell->managescalingshifts) {
    1662:       shell->ops->getdiagonalblock = (PetscErrorCode (*)(Mat, Mat *))f;
    1663:       mat->ops->getdiagonalblock   = MatGetDiagonalBlock_Shell;
    1664:     } else {
    1665:       shell->ops->getdiagonalblock = NULL;
    1666:       mat->ops->getdiagonalblock   = (PetscErrorCode (*)(Mat, Mat *))f;
    1667:     }
    1668:     break;
    1669:   case MATOP_MULT:
    1670:     if (shell->managescalingshifts) {
    1671:       shell->ops->mult = (PetscErrorCode (*)(Mat, Vec, Vec))f;
    1672:       mat->ops->mult   = MatMult_Shell;
    1673:     } else {
    1674:       shell->ops->mult = NULL;
    1675:       mat->ops->mult   = (PetscErrorCode (*)(Mat, Vec, Vec))f;
    1676:     }
    1677:     break;
    1678:   case MATOP_MULT_TRANSPOSE:
    1679:     if (shell->managescalingshifts) {
    1680:       shell->ops->multtranspose = (PetscErrorCode (*)(Mat, Vec, Vec))f;
    1681:       mat->ops->multtranspose   = MatMultTranspose_Shell;
    1682:     } else {
    1683:       shell->ops->multtranspose = NULL;
    1684:       mat->ops->multtranspose   = (PetscErrorCode (*)(Mat, Vec, Vec))f;
    1685:     }
    1686:     break;
    1687:   case MATOP_MULT_HERMITIAN_TRANSPOSE:
    1688:     if (shell->managescalingshifts) {
    1689:       shell->ops->multhermitiantranspose = (PetscErrorCode (*)(Mat, Vec, Vec))f;
    1690:       mat->ops->multhermitiantranspose   = MatMultHermitianTranspose_Shell;
    1691:     } else {
    1692:       shell->ops->multhermitiantranspose = NULL;
    1693:       mat->ops->multhermitiantranspose   = (PetscErrorCode (*)(Mat, Vec, Vec))f;
    1694:     }
    1695:     break;
    1696:   default:
    1697:     (((void (**)(void))mat->ops)[op]) = f;
    1698:     break;
    1699:   }
    1700:   PetscFunctionReturn(PETSC_SUCCESS);
    1701: }

    1703: static PetscErrorCode MatShellGetOperation_Shell(Mat mat, MatOperation op, void (**f)(void))
    1704: {
    1705:   Mat_Shell *shell = (Mat_Shell *)mat->data;

    1707:   PetscFunctionBegin;
    1708:   switch (op) {
    1709:   case MATOP_DESTROY:
    1710:     *f = (void (*)(void))shell->ops->destroy;
    1711:     break;
    1712:   case MATOP_VIEW:
    1713:     *f = (void (*)(void))mat->ops->view;
    1714:     break;
    1715:   case MATOP_COPY:
    1716:     *f = (void (*)(void))shell->ops->copy;
    1717:     break;
    1718:   case MATOP_DIAGONAL_SET:
    1719:   case MATOP_DIAGONAL_SCALE:
    1720:   case MATOP_SHIFT:
    1721:   case MATOP_SCALE:
    1722:   case MATOP_AXPY:
    1723:   case MATOP_ZERO_ROWS:
    1724:   case MATOP_ZERO_ROWS_COLUMNS:
    1725:     *f = (((void (**)(void))mat->ops)[op]);
    1726:     break;
    1727:   case MATOP_GET_DIAGONAL:
    1728:     if (shell->ops->getdiagonal) *f = (void (*)(void))shell->ops->getdiagonal;
    1729:     else *f = (((void (**)(void))mat->ops)[op]);
    1730:     break;
    1731:   case MATOP_GET_DIAGONAL_BLOCK:
    1732:     if (shell->ops->getdiagonalblock) *f = (void (*)(void))shell->ops->getdiagonalblock;
    1733:     else *f = (((void (**)(void))mat->ops)[op]);
    1734:     break;
    1735:   case MATOP_MULT:
    1736:     if (shell->ops->mult) *f = (void (*)(void))shell->ops->mult;
    1737:     else *f = (((void (**)(void))mat->ops)[op]);
    1738:     break;
    1739:   case MATOP_MULT_TRANSPOSE:
    1740:     if (shell->ops->multtranspose) *f = (void (*)(void))shell->ops->multtranspose;
    1741:     else *f = (((void (**)(void))mat->ops)[op]);
    1742:     break;
    1743:   case MATOP_MULT_HERMITIAN_TRANSPOSE:
    1744:     if (shell->ops->multhermitiantranspose) *f = (void (*)(void))shell->ops->multhermitiantranspose;
    1745:     else *f = (((void (**)(void))mat->ops)[op]);
    1746:     break;
    1747:   default:
    1748:     *f = (((void (**)(void))mat->ops)[op]);
    1749:   }
    1750:   PetscFunctionReturn(PETSC_SUCCESS);
    1751: }

    1753: /*MC
    1754:   MATSHELL - MATSHELL = "shell" - A matrix type to be used to define your own matrix type with its own data structure -- perhaps matrix-free.

    1756:   Level: advanced

    1758:   Notes:
    1759:   See `MatCreateShell()` for details on the usage of `MATSHELL`

    1761:   `PCSHELL` can be used in conjunction with `MATSHELL` to provide a custom preconditioner appropriate for your `MATSHELL`. Since
    1762:   many standard preconditioners such as `PCILU` depend on having an explicit representation of the matrix entries they cannot be used
    1763:   directly with `MATSHELL`.

    1765: .seealso: [](ch_matrices), `Mat`, `MatCreateShell()`, `PCSHELL`
    1766: M*/

    1768: PETSC_EXTERN PetscErrorCode MatCreate_Shell(Mat A)
    1769: {
    1770:   Mat_Shell *b;

    1772:   PetscFunctionBegin;
    1773:   PetscCall(PetscNew(&b));
    1774:   A->data   = (void *)b;
    1775:   A->ops[0] = MatOps_Values;

    1777:   b->ctxcontainer        = NULL;
    1778:   b->vshift              = 0.0;
    1779:   b->vscale              = 1.0;
    1780:   b->managescalingshifts = PETSC_TRUE;
    1781:   A->assembled           = PETSC_TRUE;
    1782:   A->preallocated        = PETSC_FALSE;

    1784:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellGetContext_C", MatShellGetContext_Shell));
    1785:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetContext_C", MatShellSetContext_Shell));
    1786:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetContextDestroy_C", MatShellSetContextDestroy_Shell));
    1787:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetVecType_C", MatShellSetVecType_Shell));
    1788:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetManageScalingShifts_C", MatShellSetManageScalingShifts_Shell));
    1789:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellGetScalingShifts_C", MatShellGetScalingShifts_Shell));
    1790:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetOperation_C", MatShellSetOperation_Shell));
    1791:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellGetOperation_C", MatShellGetOperation_Shell));
    1792:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetMatProductOperation_C", MatShellSetMatProductOperation_Shell));
    1793:   PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATSHELL));
    1794:   PetscFunctionReturn(PETSC_SUCCESS);
    1795: }

    1797: /*@C
    1798:   MatCreateShell - Creates a new matrix of `MatType` `MATSHELL` for use with a user-defined
    1799:   private matrix data storage format.

    1801:   Collective

    1803:   Input Parameters:
    1804: + comm - MPI communicator
    1805: . m    - number of local rows (or `PETSC_DECIDE` to have calculated if `M` is given)
    1806: . n    - number of local columns (or `PETSC_DECIDE` to have calculated if `N` is given)
    1807: . M    - number of global rows (may be `PETSC_DETERMINE` to have calculated if `m` is given)
    1808: . N    - number of global columns (may be `PETSC_DETERMINE` to have calculated if `n` is given)
    1809: - ctx  - pointer to data needed by the shell matrix routines

    1811:   Output Parameter:
    1812: . A - the matrix

    1814:   Level: advanced

    1816:   Example Usage:
    1817: .vb
    1818:   extern PetscErrorCode mult(Mat, Vec, Vec);

    1820:   MatCreateShell(comm, m, n, M, N, ctx, &mat);
    1821:   MatShellSetOperation(mat, MATOP_MULT, (void(*)(void))mult);
    1822:   MatShellSetContext(mat,ctx);
    1823:   // Use matrix for operations that have been set
    1824:   MatDestroy(mat);
    1825: .ve

    1827:   Notes:
    1828:   The shell matrix type is intended to provide a simple way for users to write a custom matrix specifically for their application.

    1830:   `MatCreateShell()` is used in conjunction with `MatShellSetContext()` and `MatShellSetOperation()`.

    1832:   PETSc requires that matrices and vectors being used for certain
    1833:   operations are partitioned accordingly.  For example, when
    1834:   creating a shell matrix, `A`, that supports parallel matrix-vector
    1835:   products using `MatMult`(A,x,y) the user should set the number
    1836:   of local matrix rows to be the number of local elements of the
    1837:   corresponding result vector, y. Note that this is information is
    1838:   required for use of the matrix interface routines, even though
    1839:   the shell matrix may not actually be physically partitioned.
    1840:   For example,

    1842: .vb
    1843:      Vec x, y
    1844:      extern PetscErrorCode mult(Mat,Vec,Vec);
    1845:      Mat A

    1847:      VecCreateMPI(comm,PETSC_DECIDE,M,&y);
    1848:      VecCreateMPI(comm,PETSC_DECIDE,N,&x);
    1849:      VecGetLocalSize(y,&m);
    1850:      VecGetLocalSize(x,&n);
    1851:      MatCreateShell(comm,m,n,M,N,ctx,&A);
    1852:      MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
    1853:      MatMult(A,x,y);
    1854:      MatDestroy(&A);
    1855:      VecDestroy(&y);
    1856:      VecDestroy(&x);
    1857: .ve

    1859:   `MATSHELL` handles `MatShift()`, `MatDiagonalSet()`, `MatDiagonalScale()`, `MatAXPY()`, `MatScale()`, `MatZeroRows()` and `MatZeroRowsColumns()`
    1860:   internally, so these  operations cannot be overwritten unless `MatShellSetManageScalingShifts()` is called.

    1862:   Developer Notes:
    1863:   For rectangular matrices do all the scalings and shifts make sense?

    1865:   Regarding shifting and scaling. The general form is

    1867:   diag(left)(vscale*A + diag(dshift) + vshift I)diag(right)

    1869:   The order you apply the operations is important. For example if you have a dshift then
    1870:   apply a MatScale(s) you get s*vscale*A + s*diag(shift). But if you first scale and then shift
    1871:   you get s*vscale*A + diag(shift)

    1873:   A is the user provided function.

    1875:   `KSP`/`PC` uses changes in the `Mat`'s "state" to decide if preconditioners need to be rebuilt `PCSetUp()` only calls the setup() for
    1876:   for the `PC` implementation if the `Mat` state has increased from the previous call. Thus to get changes in a `MATSHELL` to trigger
    1877:   an update in the preconditioner you must call `MatAssemblyBegin()` and `MatAssemblyEnd()` or `PetscObjectStateIncrease`((`PetscObject`)mat);
    1878:   each time the `MATSHELL` matrix has changed.

    1880:   Matrix-matrix product operations can be specified using `MatShellSetMatProductOperation()`

    1882:   Calling `MatAssemblyBegin()`/`MatAssemblyEnd()` on a `MATSHELL` removes any previously supplied shift and scales that were provided
    1883:   with `MatDiagonalSet()`, `MatShift()`, `MatScale()`, or `MatDiagonalScale()`.

    1885:   Fortran Notes:
    1886:   To use this from Fortran with a `ctx` you must write an interface definition for this
    1887:   function and for `MatShellGetContext()` that tells Fortran the Fortran derived data type you are passing
    1888:   in as the `ctx` argument.

    1890: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatShellSetOperation()`, `MatHasOperation()`, `MatShellGetContext()`, `MatShellSetContext()`, `MatShellSetManageScalingShifts()`, `MatShellSetMatProductOperation()`
    1891: @*/
    1892: PetscErrorCode MatCreateShell(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, PetscInt N, void *ctx, Mat *A)
    1893: {
    1894:   PetscFunctionBegin;
    1895:   PetscCall(MatCreate(comm, A));
    1896:   PetscCall(MatSetSizes(*A, m, n, M, N));
    1897:   PetscCall(MatSetType(*A, MATSHELL));
    1898:   PetscCall(MatShellSetContext(*A, ctx));
    1899:   PetscCall(MatSetUp(*A));
    1900:   PetscFunctionReturn(PETSC_SUCCESS);
    1901: }

    1903: /*@
    1904:   MatShellSetContext - sets the context for a `MATSHELL` shell matrix

    1906:   Logically Collective

    1908:   Input Parameters:
    1909: + mat - the `MATSHELL` shell matrix
    1910: - ctx - the context

    1912:   Level: advanced

    1914:   Note:
    1915:   This provides an easy way, along with `MatCreateShell()` and `MatShellSetOperation()` to provide a custom matrix format
    1916:   specifically for your application.

    1918:   Fortran Notes:
    1919:   You must write a Fortran interface definition for this
    1920:   function that tells Fortran the Fortran derived data type that you are passing in as the `ctx` argument.

    1922: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`
    1923: @*/
    1924: PetscErrorCode MatShellSetContext(Mat mat, void *ctx)
    1925: {
    1926:   PetscFunctionBegin;
    1928:   PetscTryMethod(mat, "MatShellSetContext_C", (Mat, void *), (mat, ctx));
    1929:   PetscFunctionReturn(PETSC_SUCCESS);
    1930: }

    1932: /*@C
    1933:   MatShellSetContextDestroy - sets the destroy function for a `MATSHELL` shell matrix context

    1935:   Logically Collective

    1937:   Input Parameters:
    1938: + mat - the shell matrix
    1939: - f   - the context destroy function

    1941:   Level: advanced

    1943:   Note:
    1944:   If the `MatShell` is never duplicated, the behavior of this function is equivalent
    1945:   to `MatShellSetOperation`(`Mat`,`MATOP_DESTROY`,f). However, `MatShellSetContextDestroy()`
    1946:   ensures proper reference counting for the user provided context data in the case that
    1947:   the `MATSHELL` is duplicated.

    1949: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellSetContext()`
    1950: @*/
    1951: PetscErrorCode MatShellSetContextDestroy(Mat mat, PetscErrorCode (*f)(void *))
    1952: {
    1953:   PetscFunctionBegin;
    1955:   PetscTryMethod(mat, "MatShellSetContextDestroy_C", (Mat, PetscErrorCode (*)(void *)), (mat, f));
    1956:   PetscFunctionReturn(PETSC_SUCCESS);
    1957: }

    1959: /*@C
    1960:   MatShellSetVecType - Sets the `VecType` of `Vec` returned by `MatCreateVecs()`

    1962:   Logically Collective

    1964:   Input Parameters:
    1965: + mat   - the `MATSHELL` shell matrix
    1966: - vtype - type to use for creating vectors

    1968:   Level: advanced

    1970: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateVecs()`
    1971: @*/
    1972: PetscErrorCode MatShellSetVecType(Mat mat, VecType vtype)
    1973: {
    1974:   PetscFunctionBegin;
    1975:   PetscTryMethod(mat, "MatShellSetVecType_C", (Mat, VecType), (mat, vtype));
    1976:   PetscFunctionReturn(PETSC_SUCCESS);
    1977: }

    1979: /*@
    1980:   MatShellSetManageScalingShifts - Allows the user to control the scaling and shift operations of the `MATSHELL`. Must be called immediately
    1981:   after `MatCreateShell()`

    1983:   Logically Collective

    1985:   Input Parameter:
    1986: . A - the `MATSHELL` shell matrix

    1988:   Level: advanced

    1990: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellSetContext()`, `MatShellSetOperation()`
    1991: @*/
    1992: PetscErrorCode MatShellSetManageScalingShifts(Mat A)
    1993: {
    1994:   PetscFunctionBegin;
    1996:   PetscTryMethod(A, "MatShellSetManageScalingShifts_C", (Mat), (A));
    1997:   PetscFunctionReturn(PETSC_SUCCESS);
    1998: }

    2000: /*@C
    2001:   MatShellGetScalingShifts - Gets members of a `MATSHELL` used internally for scaling and
    2002:   shifting the `Mat` or calling `MatAXPY()`, `MatZeroRows()`, or `MatZeroRowsColumns()` with it

    2004:   Logically Collective

    2006:   Input Parameter:
    2007: . A - the `MATSHELL` shell matrix

    2009:   Output Parameters:
    2010: + vshift - `PetscScalar` pointer (can be `NULL`), see `MatShift()`, or `MAT_SHELL_NOT_ALLOWED` if the internal shift should be 0
    2011: . vscale - `PetscScalar` pointer (can be `NULL`), see `MatScale()`, or `MAT_SHELL_NOT_ALLOWED` if the internal scaling should be 1
    2012: . dshift - `Vec` pointer (can be `NULL`), see `MatDiagonalSet()`, or `MAT_SHELL_NOT_ALLOWED` if the internal shift should be `NULL`
    2013: . left   - `Vec` pointer (can be `NULL`), see `MatDiagonalScale()`, or `MAT_SHELL_NOT_ALLOWED` if the internal scaling should be `NULL`
    2014: . right  - `Vec` pointer (can be `NULL`), see `MatDiagonalScale()`, or `MAT_SHELL_NOT_ALLOWED` if the internal scaling should be `NULL`
    2015: . axpy   - `Mat` pointer (can be `NULL`), or `MAT_SHELL_NOT_ALLOWED` if `MatAXPY()` should have not been called on `A`
    2016: . zrows  - `Vec` pointer (can be `NULL`), or `MAT_SHELL_NOT_ALLOWED` if `MatZeroRows()` should have not been called on `A`
    2017: - zcols  - `Vec` pointer (can be `NULL`), or `MAT_SHELL_NOT_ALLOWED` if `MatZeroRowsColumns()` should have not been called on `A`

    2019:   Level: advanced

    2021:   Developer Notes:
    2022:   This is mostly useful to check for corner-cases in `MatType` deriving from
    2023:   `MATSHELL`, e.g, `MATCOMPOSITE` or `MATTRANSPOSEVIRTUAL`, since scaling and
    2024:   shifts often require extra work which is not always implemented.

    2026: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShift()`, `MatScale()`, `MatDiagonalSet()`, `MatDiagonalScale()`, `MatAXPY()`, `MatZeroRows()`, `MatZeroRowsColumns()`, `MatShellSetManageScalingShifts()`
    2027: @*/
    2028: PetscErrorCode MatShellGetScalingShifts(Mat A, PetscScalar *vshift, PetscScalar *vscale, Vec *dshift, Vec *left, Vec *right, Mat *axpy, IS *zrows, IS *zcols)
    2029: {
    2030:   PetscFunctionBegin;
    2032:   PetscTryMethod(A, "MatShellGetScalingShifts_C", (Mat, PetscScalar *, PetscScalar *, Vec *, Vec *, Vec *, Mat *, IS *, IS *), (A, vshift, vscale, dshift, left, right, axpy, zrows, zcols));
    2033:   PetscFunctionReturn(PETSC_SUCCESS);
    2034: }

    2036: /*@C
    2037:   MatShellTestMult - Compares the multiply routine provided to the `MATSHELL` with differencing on a given function.

    2039:   Logically Collective; No Fortran Support

    2041:   Input Parameters:
    2042: + mat  - the `MATSHELL` shell matrix
    2043: . f    - the function
    2044: . base - differences are computed around this vector, see `MatMFFDSetBase()`, for Jacobians this is the point at which the Jacobian is being evaluated
    2045: - ctx  - an optional context for the function

    2047:   Output Parameter:
    2048: . flg - `PETSC_TRUE` if the multiply is likely correct

    2050:   Options Database Key:
    2051: . -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference

    2053:   Level: advanced

    2055: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMultTranspose()`
    2056: @*/
    2057: PetscErrorCode MatShellTestMult(Mat mat, PetscErrorCode (*f)(void *, Vec, Vec), Vec base, void *ctx, PetscBool *flg)
    2058: {
    2059:   PetscInt  m, n;
    2060:   Mat       mf, Dmf, Dmat, Ddiff;
    2061:   PetscReal Diffnorm, Dmfnorm;
    2062:   PetscBool v = PETSC_FALSE, flag = PETSC_TRUE;

    2064:   PetscFunctionBegin;
    2066:   PetscCall(PetscOptionsHasName(NULL, ((PetscObject)mat)->prefix, "-mat_shell_test_mult_view", &v));
    2067:   PetscCall(MatGetLocalSize(mat, &m, &n));
    2068:   PetscCall(MatCreateMFFD(PetscObjectComm((PetscObject)mat), m, n, PETSC_DECIDE, PETSC_DECIDE, &mf));
    2069:   PetscCall(MatMFFDSetFunction(mf, f, ctx));
    2070:   PetscCall(MatMFFDSetBase(mf, base, NULL));

    2072:   PetscCall(MatComputeOperator(mf, MATAIJ, &Dmf));
    2073:   PetscCall(MatComputeOperator(mat, MATAIJ, &Dmat));

    2075:   PetscCall(MatDuplicate(Dmat, MAT_COPY_VALUES, &Ddiff));
    2076:   PetscCall(MatAXPY(Ddiff, -1.0, Dmf, DIFFERENT_NONZERO_PATTERN));
    2077:   PetscCall(MatNorm(Ddiff, NORM_FROBENIUS, &Diffnorm));
    2078:   PetscCall(MatNorm(Dmf, NORM_FROBENIUS, &Dmfnorm));
    2079:   if (Diffnorm / Dmfnorm > 10 * PETSC_SQRT_MACHINE_EPSILON) {
    2080:     flag = PETSC_FALSE;
    2081:     if (v) {
    2082:       PetscCall(PetscPrintf(PetscObjectComm((PetscObject)mat), "MATSHELL and matrix-free multiple appear to produce different results.\n Norm Ratio %g Difference results followed by finite difference one\n", (double)(Diffnorm / Dmfnorm)));
    2083:       PetscCall(MatViewFromOptions(Ddiff, (PetscObject)mat, "-mat_shell_test_mult_view"));
    2084:       PetscCall(MatViewFromOptions(Dmf, (PetscObject)mat, "-mat_shell_test_mult_view"));
    2085:       PetscCall(MatViewFromOptions(Dmat, (PetscObject)mat, "-mat_shell_test_mult_view"));
    2086:     }
    2087:   } else if (v) {
    2088:     PetscCall(PetscPrintf(PetscObjectComm((PetscObject)mat), "MATSHELL and matrix-free multiple appear to produce the same results\n"));
    2089:   }
    2090:   if (flg) *flg = flag;
    2091:   PetscCall(MatDestroy(&Ddiff));
    2092:   PetscCall(MatDestroy(&mf));
    2093:   PetscCall(MatDestroy(&Dmf));
    2094:   PetscCall(MatDestroy(&Dmat));
    2095:   PetscFunctionReturn(PETSC_SUCCESS);
    2096: }

    2098: /*@C
    2099:   MatShellTestMultTranspose - Compares the multiply transpose routine provided to the `MATSHELL` with differencing on a given function.

    2101:   Logically Collective; No Fortran Support

    2103:   Input Parameters:
    2104: + mat  - the `MATSHELL` shell matrix
    2105: . f    - the function
    2106: . base - differences are computed around this vector, see `MatMFFDSetBase()`, for Jacobians this is the point at which the Jacobian is being evaluated
    2107: - ctx  - an optional context for the function

    2109:   Output Parameter:
    2110: . flg - `PETSC_TRUE` if the multiply is likely correct

    2112:   Options Database Key:
    2113: . -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference

    2115:   Level: advanced

    2117: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMult()`
    2118: @*/
    2119: PetscErrorCode MatShellTestMultTranspose(Mat mat, PetscErrorCode (*f)(void *, Vec, Vec), Vec base, void *ctx, PetscBool *flg)
    2120: {
    2121:   Vec       x, y, z;
    2122:   PetscInt  m, n, M, N;
    2123:   Mat       mf, Dmf, Dmat, Ddiff;
    2124:   PetscReal Diffnorm, Dmfnorm;
    2125:   PetscBool v = PETSC_FALSE, flag = PETSC_TRUE;

    2127:   PetscFunctionBegin;
    2129:   PetscCall(PetscOptionsHasName(NULL, ((PetscObject)mat)->prefix, "-mat_shell_test_mult_transpose_view", &v));
    2130:   PetscCall(MatCreateVecs(mat, &x, &y));
    2131:   PetscCall(VecDuplicate(y, &z));
    2132:   PetscCall(MatGetLocalSize(mat, &m, &n));
    2133:   PetscCall(MatGetSize(mat, &M, &N));
    2134:   PetscCall(MatCreateMFFD(PetscObjectComm((PetscObject)mat), m, n, M, N, &mf));
    2135:   PetscCall(MatMFFDSetFunction(mf, f, ctx));
    2136:   PetscCall(MatMFFDSetBase(mf, base, NULL));
    2137:   PetscCall(MatComputeOperator(mf, MATAIJ, &Dmf));
    2138:   PetscCall(MatTranspose(Dmf, MAT_INPLACE_MATRIX, &Dmf));
    2139:   PetscCall(MatComputeOperatorTranspose(mat, MATAIJ, &Dmat));

    2141:   PetscCall(MatDuplicate(Dmat, MAT_COPY_VALUES, &Ddiff));
    2142:   PetscCall(MatAXPY(Ddiff, -1.0, Dmf, DIFFERENT_NONZERO_PATTERN));
    2143:   PetscCall(MatNorm(Ddiff, NORM_FROBENIUS, &Diffnorm));
    2144:   PetscCall(MatNorm(Dmf, NORM_FROBENIUS, &Dmfnorm));
    2145:   if (Diffnorm / Dmfnorm > 10 * PETSC_SQRT_MACHINE_EPSILON) {
    2146:     flag = PETSC_FALSE;
    2147:     if (v) {
    2148:       PetscCall(PetscPrintf(PetscObjectComm((PetscObject)mat), "MATSHELL and matrix-free multiple appear to produce different results.\n Norm Ratio %g Difference results followed by finite difference one\n", (double)(Diffnorm / Dmfnorm)));
    2149:       PetscCall(MatViewFromOptions(Ddiff, (PetscObject)mat, "-mat_shell_test_mult_transpose_view"));
    2150:       PetscCall(MatViewFromOptions(Dmf, (PetscObject)mat, "-mat_shell_test_mult_transpose_view"));
    2151:       PetscCall(MatViewFromOptions(Dmat, (PetscObject)mat, "-mat_shell_test_mult_transpose_view"));
    2152:     }
    2153:   } else if (v) {
    2154:     PetscCall(PetscPrintf(PetscObjectComm((PetscObject)mat), "MATSHELL transpose and matrix-free multiple appear to produce the same results\n"));
    2155:   }
    2156:   if (flg) *flg = flag;
    2157:   PetscCall(MatDestroy(&mf));
    2158:   PetscCall(MatDestroy(&Dmat));
    2159:   PetscCall(MatDestroy(&Ddiff));
    2160:   PetscCall(MatDestroy(&Dmf));
    2161:   PetscCall(VecDestroy(&x));
    2162:   PetscCall(VecDestroy(&y));
    2163:   PetscCall(VecDestroy(&z));
    2164:   PetscFunctionReturn(PETSC_SUCCESS);
    2165: }

    2167: /*@C
    2168:   MatShellSetOperation - Allows user to set a matrix operation for a `MATSHELL` shell matrix.

    2170:   Logically Collective

    2172:   Input Parameters:
    2173: + mat - the `MATSHELL` shell matrix
    2174: . op  - the name of the operation
    2175: - g   - the function that provides the operation.

    2177:   Level: advanced

    2179:   Example Usage:
    2180: .vb
    2181:   extern PetscErrorCode usermult(Mat, Vec, Vec);

    2183:   MatCreateShell(comm, m, n, M, N, ctx, &A);
    2184:   MatShellSetOperation(A, MATOP_MULT, (void(*)(void))usermult);
    2185: .ve

    2187:   Notes:
    2188:   See the file include/petscmat.h for a complete list of matrix
    2189:   operations, which all have the form MATOP_<OPERATION>, where
    2190:   <OPERATION> is the name (in all capital letters) of the
    2191:   user interface routine (e.g., `MatMult()` -> `MATOP_MULT`).

    2193:   All user-provided functions (except for `MATOP_DESTROY`) should have the same calling
    2194:   sequence as the usual matrix interface routines, since they
    2195:   are intended to be accessed via the usual matrix interface
    2196:   routines, e.g.,
    2197: $       MatMult(Mat, Vec, Vec) -> usermult(Mat, Vec, Vec)

    2199:   In particular each function MUST return an error code of 0 on success and
    2200:   nonzero on failure.

    2202:   Within each user-defined routine, the user should call
    2203:   `MatShellGetContext()` to obtain the user-defined context that was
    2204:   set by `MatCreateShell()`.

    2206:   Use `MatSetOperation()` to set an operation for any matrix type. For matrix product operations (i.e. `MatMatXXX()`, `MatTransposeMatXXX()` etc)
    2207:   use `MatShellSetMatProductOperation()`

    2209:   Fortran Notes:
    2210:   For `MatCreateVecs()` the user code should check if the input left or right matrix is -1 and in that case not
    2211:   generate a matrix. See src/mat/tests/ex120f.F

    2213: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellSetContext()`, `MatSetOperation()`, `MatShellSetManageScalingShifts()`, `MatShellSetMatProductOperation()`
    2214: @*/
    2215: PetscErrorCode MatShellSetOperation(Mat mat, MatOperation op, void (*g)(void))
    2216: {
    2217:   PetscFunctionBegin;
    2219:   PetscTryMethod(mat, "MatShellSetOperation_C", (Mat, MatOperation, void (*)(void)), (mat, op, g));
    2220:   PetscFunctionReturn(PETSC_SUCCESS);
    2221: }

    2223: /*@C
    2224:   MatShellGetOperation - Gets a matrix function for a `MATSHELL` shell matrix.

    2226:   Not Collective

    2228:   Input Parameters:
    2229: + mat - the `MATSHELL` shell matrix
    2230: - op  - the name of the operation

    2232:   Output Parameter:
    2233: . g - the function that provides the operation.

    2235:   Level: advanced

    2237:   Notes:
    2238:   See the file include/petscmat.h for a complete list of matrix
    2239:   operations, which all have the form MATOP_<OPERATION>, where
    2240:   <OPERATION> is the name (in all capital letters) of the
    2241:   user interface routine (e.g., `MatMult()` -> `MATOP_MULT`).

    2243:   All user-provided functions have the same calling
    2244:   sequence as the usual matrix interface routines, since they
    2245:   are intended to be accessed via the usual matrix interface
    2246:   routines, e.g.,
    2247: $       MatMult(Mat, Vec, Vec) -> usermult(Mat, Vec, Vec)

    2249:   Within each user-defined routine, the user should call
    2250:   `MatShellGetContext()` to obtain the user-defined context that was
    2251:   set by `MatCreateShell()`.

    2253: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellSetOperation()`, `MatShellSetContext()`
    2254: @*/
    2255: PetscErrorCode MatShellGetOperation(Mat mat, MatOperation op, void (**g)(void))
    2256: {
    2257:   PetscFunctionBegin;
    2259:   PetscUseMethod(mat, "MatShellGetOperation_C", (Mat, MatOperation, void (**)(void)), (mat, op, g));
    2260:   PetscFunctionReturn(PETSC_SUCCESS);
    2261: }

    2263: /*@
    2264:   MatIsShell - Inquires if a matrix is derived from `MATSHELL`

    2266:   Input Parameter:
    2267: . mat - the matrix

    2269:   Output Parameter:
    2270: . flg - the Boolean value

    2272:   Level: developer

    2274: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MATMFFD`, `MatCreateShell()`, `MATTRANSPOSEVIRTUAL`, `MATSCHURCOMPLEMENT`
    2275: @*/
    2276: PetscErrorCode MatIsShell(Mat mat, PetscBool *flg)
    2277: {
    2278:   PetscFunctionBegin;
    2280:   PetscAssertPointer(flg, 2);
    2281:   *flg = (PetscBool)(mat->ops->destroy == MatDestroy_Shell);
    2282:   PetscFunctionReturn(PETSC_SUCCESS);
    2283: }