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:   Mat_Shell       *shell;
510:   PetscBool        useBmdata = PETSC_FALSE, newB = PETSC_TRUE;
511:   PetscErrorCode (*stashsym)(Mat), (*stashnum)(Mat);

513:   PetscFunctionBegin;
514:   MatCheckProduct(D, 1);
515:   product = D->product;
516:   PetscCheck(product->data, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "Product data empty");
517:   A     = product->A;
518:   B     = product->B;
519:   mdata = (MatMatDataShell *)product->data;
520:   PetscCheck(mdata->numeric, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "Missing numeric operation");
521:   shell    = (Mat_Shell *)A->data;
522:   stashsym = D->ops->productsymbolic;
523:   stashnum = D->ops->productnumeric;
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 != 0.0) {
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 != 0.0) {
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 == 0.0, 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 != 0.0) {
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:   PetscFunctionReturn(PETSC_SUCCESS);
671: }

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

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

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

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

776: 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)
777: {
778:   PetscBool               flg;
779:   Mat_Shell              *shell;
780:   MatShellMatFunctionList matmat;

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

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

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

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

822:   Logically Collective; No Fortran Support

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

833:   Level: advanced

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

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

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

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

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

864: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellSetContext()`, `MatSetOperation()`, `MatProductType`, `MatType`, `MatSetUp()`
865: @*/
866: 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)
867: {
868:   PetscFunctionBegin;
871:   PetscCheck(ptype != MATPRODUCT_ABC, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not for product type %s", MatProductTypes[ptype]);
872:   PetscCheck(numeric, PetscObjectComm((PetscObject)A), PETSC_ERR_USER, "Missing numeric routine, argument 4");
873:   PetscAssertPointer(Btype, 6);
874:   if (Ctype) PetscAssertPointer(Ctype, 7);
875:   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));
876:   PetscFunctionReturn(PETSC_SUCCESS);
877: }

879: 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)
880: {
881:   PetscBool   flg;
882:   char        composedname[256];
883:   MatRootName Bnames = MatRootNameList, Cnames = MatRootNameList;
884:   PetscMPIInt size;

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

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

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

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

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

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

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

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

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

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

1007:     PetscCall(MatShellGetContext(shell->axpy, &X));
1008:     PetscCall(PetscObjectStateGet((PetscObject)X, &axpy_state));
1009:     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,...)");

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

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

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

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

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

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

1059:     PetscCall(MatShellGetContext(shell->axpy, &X));
1060:     PetscCall(PetscObjectStateGet((PetscObject)X, &axpy_state));
1061:     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,...)");
1062:     PetscCall(MatCreateVecs(shell->axpy, shell->axpy_right ? NULL : &shell->axpy_right, shell->axpy_left ? NULL : &shell->axpy_left));
1063:     PetscCall(VecCopy(x, shell->axpy_left));
1064:     PetscCall(MatMultTranspose(shell->axpy, shell->axpy_left, shell->axpy_right));
1065:     PetscCall(VecAXPY(y, shell->axpy_vscale, shell->axpy_right));
1066:   }
1067:   PetscFunctionReturn(PETSC_SUCCESS);
1068: }

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

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

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

1094:     PetscCall(MatShellGetContext(shell->axpy, &X));
1095:     PetscCall(PetscObjectStateGet((PetscObject)X, &axpy_state));
1096:     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,...)");
1097:     PetscCall(MatCreateVecs(shell->axpy, shell->axpy_right ? NULL : &shell->axpy_right, shell->axpy_left ? NULL : &shell->axpy_left));
1098:     PetscCall(VecCopy(x, shell->axpy_left));
1099:     PetscCall(MatMultHermitianTranspose(shell->axpy, shell->axpy_left, shell->axpy_right));
1100:     PetscCall(VecAXPY(y, PetscConj(shell->axpy_vscale), shell->axpy_right));
1101:   }
1102:   PetscFunctionReturn(PETSC_SUCCESS);
1103: }

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

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

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

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

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

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

1160:     PetscCall(MatShellGetContext(shell->axpy, &X));
1161:     PetscCall(PetscObjectStateGet((PetscObject)X, &axpy_state));
1162:     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,...)");
1163:     PetscCall(MatCreateVecs(shell->axpy, NULL, shell->axpy_left ? NULL : &shell->axpy_left));
1164:     PetscCall(MatGetDiagonal(shell->axpy, shell->axpy_left));
1165:     PetscCall(VecAXPY(v, shell->axpy_vscale, shell->axpy_left));
1166:   }
1167:   PetscFunctionReturn(PETSC_SUCCESS);
1168: }

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

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

1203: static PetscErrorCode MatShift_Shell(Mat Y, PetscScalar a)
1204: {
1205:   Mat_Shell *shell = (Mat_Shell *)Y->data;
1206:   PetscBool  flg;

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

1227: static PetscErrorCode MatDiagonalSet_Shell_Private(Mat A, Vec D, PetscScalar s)
1228: {
1229:   Mat_Shell *shell = (Mat_Shell *)A->data;

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

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

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

1273: static PetscErrorCode MatScale_Shell(Mat Y, PetscScalar a)
1274: {
1275:   Mat_Shell *shell = (Mat_Shell *)Y->data;

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

1286: static PetscErrorCode MatDiagonalScale_Shell(Mat Y, Vec left, Vec right)
1287: {
1288:   Mat_Shell *shell = (Mat_Shell *)Y->data;

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

1319: PETSC_INTERN PetscErrorCode MatAssemblyEnd_Shell(Mat Y, MatAssemblyType t)
1320: {
1321:   Mat_Shell *shell = (Mat_Shell *)Y->data;

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

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

1350: static PetscErrorCode MatAXPY_Shell(Mat Y, PetscScalar a, Mat X, MatStructure str)
1351: {
1352:   Mat_Shell *shell = (Mat_Shell *)Y->data;

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

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

1514: static PetscErrorCode MatShellSetContext_Shell(Mat mat, void *ctx)
1515: {
1516:   Mat_Shell *shell = (Mat_Shell *)mat->data;

1518:   PetscFunctionBegin;
1519:   if (ctx) {
1520:     PetscContainer ctxcontainer;
1521:     PetscCall(PetscContainerCreate(PetscObjectComm((PetscObject)mat), &ctxcontainer));
1522:     PetscCall(PetscContainerSetPointer(ctxcontainer, ctx));
1523:     PetscCall(PetscObjectCompose((PetscObject)mat, "MatShell ctx", (PetscObject)ctxcontainer));
1524:     shell->ctxcontainer = ctxcontainer;
1525:     PetscCall(PetscContainerDestroy(&ctxcontainer));
1526:   } else {
1527:     PetscCall(PetscObjectCompose((PetscObject)mat, "MatShell ctx", NULL));
1528:     shell->ctxcontainer = NULL;
1529:   }
1530:   PetscFunctionReturn(PETSC_SUCCESS);
1531: }

1533: static PetscErrorCode MatShellSetContextDestroy_Shell(Mat mat, PetscCtxDestroyFn *f)
1534: {
1535:   Mat_Shell *shell = (Mat_Shell *)mat->data;

1537:   PetscFunctionBegin;
1538:   if (shell->ctxcontainer) PetscCall(PetscContainerSetCtxDestroy(shell->ctxcontainer, f));
1539:   PetscFunctionReturn(PETSC_SUCCESS);
1540: }

1542: PetscErrorCode MatShellSetContext_Immutable(Mat mat, void *ctx)
1543: {
1544:   PetscFunctionBegin;
1545:   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);
1546:   PetscFunctionReturn(PETSC_SUCCESS);
1547: }

1549: PetscErrorCode MatShellSetContextDestroy_Immutable(Mat mat, PetscCtxDestroyFn *f)
1550: {
1551:   PetscFunctionBegin;
1552:   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);
1553:   PetscFunctionReturn(PETSC_SUCCESS);
1554: }

1556: PetscErrorCode MatShellSetManageScalingShifts_Immutable(Mat mat)
1557: {
1558:   PetscFunctionBegin;
1559:   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);
1560:   PetscFunctionReturn(PETSC_SUCCESS);
1561: }

1563: static PetscErrorCode MatShellSetVecType_Shell(Mat mat, VecType vtype)
1564: {
1565:   PetscFunctionBegin;
1566:   PetscCall(PetscFree(mat->defaultvectype));
1567:   PetscCall(PetscStrallocpy(vtype, (char **)&mat->defaultvectype));
1568:   PetscFunctionReturn(PETSC_SUCCESS);
1569: }

1571: static PetscErrorCode MatShellSetManageScalingShifts_Shell(Mat A)
1572: {
1573:   Mat_Shell *shell = (Mat_Shell *)A->data;

1575:   PetscFunctionBegin;
1576:   shell->managescalingshifts = PETSC_FALSE;
1577:   A->ops->diagonalset        = NULL;
1578:   A->ops->diagonalscale      = NULL;
1579:   A->ops->scale              = NULL;
1580:   A->ops->shift              = NULL;
1581:   A->ops->axpy               = NULL;
1582:   PetscFunctionReturn(PETSC_SUCCESS);
1583: }

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

1589:   PetscFunctionBegin;
1590:   PetscCheck(shell->managescalingshifts, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Cannot get scaling and shifts if MatShellSetManageScalingShifts() has been called");
1591:   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()");
1592:   else if (vshift) *vshift = shell->vshift;
1593:   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()");
1594:   else if (vscale) *vscale = shell->vscale;
1595:   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()");
1596:   else if (dshift) *dshift = shell->dshift;
1597:   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()");
1598:   else if (left) *left = shell->left;
1599:   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()");
1600:   else if (right) *right = shell->right;
1601:   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()");
1602:   else if (axpy) *axpy = shell->axpy;
1603:   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()");
1604:   else if (zrows) *zrows = shell->zrows;
1605:   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()");
1606:   else if (zcols) *zcols = shell->zcols;
1607:   PetscFunctionReturn(PETSC_SUCCESS);
1608: }

1610: static PetscErrorCode MatShellSetOperation_Shell(Mat mat, MatOperation op, void (*f)(void))
1611: {
1612:   Mat_Shell *shell = (Mat_Shell *)mat->data;

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

1688: static PetscErrorCode MatShellGetOperation_Shell(Mat mat, MatOperation op, void (**f)(void))
1689: {
1690:   Mat_Shell *shell = (Mat_Shell *)mat->data;

1692:   PetscFunctionBegin;
1693:   switch (op) {
1694:   case MATOP_DESTROY:
1695:     *f = (void (*)(void))shell->ops->destroy;
1696:     break;
1697:   case MATOP_VIEW:
1698:     *f = (void (*)(void))mat->ops->view;
1699:     break;
1700:   case MATOP_COPY:
1701:     *f = (void (*)(void))shell->ops->copy;
1702:     break;
1703:   case MATOP_DIAGONAL_SET:
1704:   case MATOP_DIAGONAL_SCALE:
1705:   case MATOP_SHIFT:
1706:   case MATOP_SCALE:
1707:   case MATOP_AXPY:
1708:   case MATOP_ZERO_ROWS:
1709:   case MATOP_ZERO_ROWS_COLUMNS:
1710:     *f = (((void (**)(void))mat->ops)[op]);
1711:     break;
1712:   case MATOP_GET_DIAGONAL:
1713:     if (shell->ops->getdiagonal) *f = (void (*)(void))shell->ops->getdiagonal;
1714:     else *f = (((void (**)(void))mat->ops)[op]);
1715:     break;
1716:   case MATOP_GET_DIAGONAL_BLOCK:
1717:     if (shell->ops->getdiagonalblock) *f = (void (*)(void))shell->ops->getdiagonalblock;
1718:     else *f = (((void (**)(void))mat->ops)[op]);
1719:     break;
1720:   case MATOP_MULT:
1721:     if (shell->ops->mult) *f = (void (*)(void))shell->ops->mult;
1722:     else *f = (((void (**)(void))mat->ops)[op]);
1723:     break;
1724:   case MATOP_MULT_TRANSPOSE:
1725:     if (shell->ops->multtranspose) *f = (void (*)(void))shell->ops->multtranspose;
1726:     else *f = (((void (**)(void))mat->ops)[op]);
1727:     break;
1728:   case MATOP_MULT_HERMITIAN_TRANSPOSE:
1729:     if (shell->ops->multhermitiantranspose) *f = (void (*)(void))shell->ops->multhermitiantranspose;
1730:     else *f = (((void (**)(void))mat->ops)[op]);
1731:     break;
1732:   default:
1733:     *f = (((void (**)(void))mat->ops)[op]);
1734:   }
1735:   PetscFunctionReturn(PETSC_SUCCESS);
1736: }

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

1741:   Level: advanced

1743:   Notes:
1744:   See `MatCreateShell()` for details on the usage of `MATSHELL`

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

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

1753: PETSC_EXTERN PetscErrorCode MatCreate_Shell(Mat A)
1754: {
1755:   Mat_Shell *b;

1757:   PetscFunctionBegin;
1758:   PetscCall(PetscNew(&b));
1759:   A->data   = (void *)b;
1760:   A->ops[0] = MatOps_Values;

1762:   b->ctxcontainer        = NULL;
1763:   b->vshift              = 0.0;
1764:   b->vscale              = 1.0;
1765:   b->managescalingshifts = PETSC_TRUE;
1766:   A->assembled           = PETSC_TRUE;
1767:   A->preallocated        = PETSC_FALSE;

1769:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellGetContext_C", MatShellGetContext_Shell));
1770:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetContext_C", MatShellSetContext_Shell));
1771:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetContextDestroy_C", MatShellSetContextDestroy_Shell));
1772:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetVecType_C", MatShellSetVecType_Shell));
1773:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetManageScalingShifts_C", MatShellSetManageScalingShifts_Shell));
1774:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellGetScalingShifts_C", MatShellGetScalingShifts_Shell));
1775:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetOperation_C", MatShellSetOperation_Shell));
1776:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellGetOperation_C", MatShellGetOperation_Shell));
1777:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetMatProductOperation_C", MatShellSetMatProductOperation_Shell));
1778:   PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATSHELL));
1779:   PetscFunctionReturn(PETSC_SUCCESS);
1780: }

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

1786:   Collective

1788:   Input Parameters:
1789: + comm - MPI communicator
1790: . m    - number of local rows (or `PETSC_DECIDE` to have calculated if `M` is given)
1791: . n    - number of local columns (or `PETSC_DECIDE` to have calculated if `N` is given)
1792: . M    - number of global rows (may be `PETSC_DETERMINE` to have calculated if `m` is given)
1793: . N    - number of global columns (may be `PETSC_DETERMINE` to have calculated if `n` is given)
1794: - ctx  - pointer to data needed by the shell matrix routines

1796:   Output Parameter:
1797: . A - the matrix

1799:   Level: advanced

1801:   Example Usage:
1802: .vb
1803:   extern PetscErrorCode mult(Mat, Vec, Vec);

1805:   MatCreateShell(comm, m, n, M, N, ctx, &mat);
1806:   MatShellSetOperation(mat, MATOP_MULT, (void(*)(void))mult);
1807:   MatShellSetContext(mat,ctx);
1808:   // Use matrix for operations that have been set
1809:   MatDestroy(mat);
1810: .ve

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

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

1817:   PETSc requires that matrices and vectors being used for certain
1818:   operations are partitioned accordingly.  For example, when
1819:   creating a shell matrix, `A`, that supports parallel matrix-vector
1820:   products using `MatMult`(A,x,y) the user should set the number
1821:   of local matrix rows to be the number of local elements of the
1822:   corresponding result vector, y. Note that this is information is
1823:   required for use of the matrix interface routines, even though
1824:   the shell matrix may not actually be physically partitioned.
1825:   For example,

1827: .vb
1828:      Vec x, y
1829:      extern PetscErrorCode mult(Mat,Vec,Vec);
1830:      Mat A

1832:      VecCreateMPI(comm,PETSC_DECIDE,M,&y);
1833:      VecCreateMPI(comm,PETSC_DECIDE,N,&x);
1834:      VecGetLocalSize(y,&m);
1835:      VecGetLocalSize(x,&n);
1836:      MatCreateShell(comm,m,n,M,N,ctx,&A);
1837:      MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
1838:      MatMult(A,x,y);
1839:      MatDestroy(&A);
1840:      VecDestroy(&y);
1841:      VecDestroy(&x);
1842: .ve

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

1847:   Developer Notes:
1848:   For rectangular matrices do all the scalings and shifts make sense?

1850:   Regarding shifting and scaling. The general form is

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

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

1858:   A is the user provided function.

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

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

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

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

1875: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatShellSetOperation()`, `MatHasOperation()`, `MatShellGetContext()`, `MatShellSetContext()`, `MatShellSetManageScalingShifts()`, `MatShellSetMatProductOperation()`
1876: @*/
1877: PetscErrorCode MatCreateShell(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, PetscInt N, void *ctx, Mat *A)
1878: {
1879:   PetscFunctionBegin;
1880:   PetscCall(MatCreate(comm, A));
1881:   PetscCall(MatSetSizes(*A, m, n, M, N));
1882:   PetscCall(MatSetType(*A, MATSHELL));
1883:   PetscCall(MatShellSetContext(*A, ctx));
1884:   PetscCall(MatSetUp(*A));
1885:   PetscFunctionReturn(PETSC_SUCCESS);
1886: }

1888: /*@
1889:   MatShellSetContext - sets the context for a `MATSHELL` shell matrix

1891:   Logically Collective

1893:   Input Parameters:
1894: + mat - the `MATSHELL` shell matrix
1895: - ctx - the context

1897:   Level: advanced

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

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

1907: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`
1908: @*/
1909: PetscErrorCode MatShellSetContext(Mat mat, void *ctx)
1910: {
1911:   PetscFunctionBegin;
1913:   PetscTryMethod(mat, "MatShellSetContext_C", (Mat, void *), (mat, ctx));
1914:   PetscFunctionReturn(PETSC_SUCCESS);
1915: }

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

1920:   Logically Collective

1922:   Input Parameters:
1923: + mat - the shell matrix
1924: - f   - the context destroy function, see `PetscCtxDestroyFn` for calling sequence

1926:   Level: advanced

1928:   Note:
1929:   If the `MatShell` is never duplicated, the behavior of this function is equivalent
1930:   to `MatShellSetOperation`(`Mat`,`MATOP_DESTROY`,f). However, `MatShellSetContextDestroy()`
1931:   ensures proper reference counting for the user provided context data in the case that
1932:   the `MATSHELL` is duplicated.

1934: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellSetContext()`,
1935:           `PetscCtxDestroyFn`
1936: @*/
1937: PetscErrorCode MatShellSetContextDestroy(Mat mat, PetscCtxDestroyFn *f)
1938: {
1939:   PetscFunctionBegin;
1941:   PetscTryMethod(mat, "MatShellSetContextDestroy_C", (Mat, PetscCtxDestroyFn *), (mat, f));
1942:   PetscFunctionReturn(PETSC_SUCCESS);
1943: }

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

1948:   Logically Collective

1950:   Input Parameters:
1951: + mat   - the `MATSHELL` shell matrix
1952: - vtype - type to use for creating vectors

1954:   Level: advanced

1956: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateVecs()`
1957: @*/
1958: PetscErrorCode MatShellSetVecType(Mat mat, VecType vtype)
1959: {
1960:   PetscFunctionBegin;
1961:   PetscTryMethod(mat, "MatShellSetVecType_C", (Mat, VecType), (mat, vtype));
1962:   PetscFunctionReturn(PETSC_SUCCESS);
1963: }

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

1969:   Logically Collective

1971:   Input Parameter:
1972: . A - the `MATSHELL` shell matrix

1974:   Level: advanced

1976: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellSetContext()`, `MatShellSetOperation()`
1977: @*/
1978: PetscErrorCode MatShellSetManageScalingShifts(Mat A)
1979: {
1980:   PetscFunctionBegin;
1982:   PetscTryMethod(A, "MatShellSetManageScalingShifts_C", (Mat), (A));
1983:   PetscFunctionReturn(PETSC_SUCCESS);
1984: }

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

1990:   Logically Collective

1992:   Input Parameter:
1993: . A - the `MATSHELL` shell matrix

1995:   Output Parameters:
1996: + vshift - `PetscScalar` pointer (can be `NULL`), see `MatShift()`, or `MAT_SHELL_NOT_ALLOWED` if the internal shift should be 0
1997: . vscale - `PetscScalar` pointer (can be `NULL`), see `MatScale()`, or `MAT_SHELL_NOT_ALLOWED` if the internal scaling should be 1
1998: . dshift - `Vec` pointer (can be `NULL`), see `MatDiagonalSet()`, or `MAT_SHELL_NOT_ALLOWED` if the internal shift should be `NULL`
1999: . left   - `Vec` pointer (can be `NULL`), see `MatDiagonalScale()`, or `MAT_SHELL_NOT_ALLOWED` if the internal scaling should be `NULL`
2000: . right  - `Vec` pointer (can be `NULL`), see `MatDiagonalScale()`, or `MAT_SHELL_NOT_ALLOWED` if the internal scaling should be `NULL`
2001: . axpy   - `Mat` pointer (can be `NULL`), or `MAT_SHELL_NOT_ALLOWED` if `MatAXPY()` should have not been called on `A`
2002: . zrows  - `Vec` pointer (can be `NULL`), or `MAT_SHELL_NOT_ALLOWED` if `MatZeroRows()` should have not been called on `A`
2003: - zcols  - `Vec` pointer (can be `NULL`), or `MAT_SHELL_NOT_ALLOWED` if `MatZeroRowsColumns()` should have not been called on `A`

2005:   Level: advanced

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

2012: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShift()`, `MatScale()`, `MatDiagonalSet()`, `MatDiagonalScale()`, `MatAXPY()`, `MatZeroRows()`, `MatZeroRowsColumns()`, `MatShellSetManageScalingShifts()`
2013: @*/
2014: PetscErrorCode MatShellGetScalingShifts(Mat A, PetscScalar *vshift, PetscScalar *vscale, Vec *dshift, Vec *left, Vec *right, Mat *axpy, IS *zrows, IS *zcols)
2015: {
2016:   PetscFunctionBegin;
2018:   PetscTryMethod(A, "MatShellGetScalingShifts_C", (Mat, PetscScalar *, PetscScalar *, Vec *, Vec *, Vec *, Mat *, IS *, IS *), (A, vshift, vscale, dshift, left, right, axpy, zrows, zcols));
2019:   PetscFunctionReturn(PETSC_SUCCESS);
2020: }

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

2025:   Logically Collective; No Fortran Support

2027:   Input Parameters:
2028: + mat  - the `MATSHELL` shell matrix
2029: . f    - the function
2030: . base - differences are computed around this vector, see `MatMFFDSetBase()`, for Jacobians this is the point at which the Jacobian is being evaluated
2031: - ctx  - an optional context for the function

2033:   Output Parameter:
2034: . flg - `PETSC_TRUE` if the multiply is likely correct

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

2039:   Level: advanced

2041: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMultTranspose()`
2042: @*/
2043: PetscErrorCode MatShellTestMult(Mat mat, PetscErrorCode (*f)(void *, Vec, Vec), Vec base, void *ctx, PetscBool *flg)
2044: {
2045:   PetscInt  m, n;
2046:   Mat       mf, Dmf, Dmat, Ddiff;
2047:   PetscReal Diffnorm, Dmfnorm;
2048:   PetscBool v = PETSC_FALSE, flag = PETSC_TRUE;

2050:   PetscFunctionBegin;
2052:   PetscCall(PetscOptionsHasName(NULL, ((PetscObject)mat)->prefix, "-mat_shell_test_mult_view", &v));
2053:   PetscCall(MatGetLocalSize(mat, &m, &n));
2054:   PetscCall(MatCreateMFFD(PetscObjectComm((PetscObject)mat), m, n, PETSC_DECIDE, PETSC_DECIDE, &mf));
2055:   PetscCall(MatMFFDSetFunction(mf, f, ctx));
2056:   PetscCall(MatMFFDSetBase(mf, base, NULL));

2058:   PetscCall(MatComputeOperator(mf, MATAIJ, &Dmf));
2059:   PetscCall(MatComputeOperator(mat, MATAIJ, &Dmat));

2061:   PetscCall(MatDuplicate(Dmat, MAT_COPY_VALUES, &Ddiff));
2062:   PetscCall(MatAXPY(Ddiff, -1.0, Dmf, DIFFERENT_NONZERO_PATTERN));
2063:   PetscCall(MatNorm(Ddiff, NORM_FROBENIUS, &Diffnorm));
2064:   PetscCall(MatNorm(Dmf, NORM_FROBENIUS, &Dmfnorm));
2065:   if (Diffnorm / Dmfnorm > 10 * PETSC_SQRT_MACHINE_EPSILON) {
2066:     flag = PETSC_FALSE;
2067:     if (v) {
2068:       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)));
2069:       PetscCall(MatViewFromOptions(Ddiff, (PetscObject)mat, "-mat_shell_test_mult_view"));
2070:       PetscCall(MatViewFromOptions(Dmf, (PetscObject)mat, "-mat_shell_test_mult_view"));
2071:       PetscCall(MatViewFromOptions(Dmat, (PetscObject)mat, "-mat_shell_test_mult_view"));
2072:     }
2073:   } else if (v) {
2074:     PetscCall(PetscPrintf(PetscObjectComm((PetscObject)mat), "MATSHELL and matrix-free multiple appear to produce the same results\n"));
2075:   }
2076:   if (flg) *flg = flag;
2077:   PetscCall(MatDestroy(&Ddiff));
2078:   PetscCall(MatDestroy(&mf));
2079:   PetscCall(MatDestroy(&Dmf));
2080:   PetscCall(MatDestroy(&Dmat));
2081:   PetscFunctionReturn(PETSC_SUCCESS);
2082: }

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

2087:   Logically Collective; No Fortran Support

2089:   Input Parameters:
2090: + mat  - the `MATSHELL` shell matrix
2091: . f    - the function
2092: . base - differences are computed around this vector, see `MatMFFDSetBase()`, for Jacobians this is the point at which the Jacobian is being evaluated
2093: - ctx  - an optional context for the function

2095:   Output Parameter:
2096: . flg - `PETSC_TRUE` if the multiply is likely correct

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

2101:   Level: advanced

2103: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMult()`
2104: @*/
2105: PetscErrorCode MatShellTestMultTranspose(Mat mat, PetscErrorCode (*f)(void *, Vec, Vec), Vec base, void *ctx, PetscBool *flg)
2106: {
2107:   Vec       x, y, z;
2108:   PetscInt  m, n, M, N;
2109:   Mat       mf, Dmf, Dmat, Ddiff;
2110:   PetscReal Diffnorm, Dmfnorm;
2111:   PetscBool v = PETSC_FALSE, flag = PETSC_TRUE;

2113:   PetscFunctionBegin;
2115:   PetscCall(PetscOptionsHasName(NULL, ((PetscObject)mat)->prefix, "-mat_shell_test_mult_transpose_view", &v));
2116:   PetscCall(MatCreateVecs(mat, &x, &y));
2117:   PetscCall(VecDuplicate(y, &z));
2118:   PetscCall(MatGetLocalSize(mat, &m, &n));
2119:   PetscCall(MatGetSize(mat, &M, &N));
2120:   PetscCall(MatCreateMFFD(PetscObjectComm((PetscObject)mat), m, n, M, N, &mf));
2121:   PetscCall(MatMFFDSetFunction(mf, f, ctx));
2122:   PetscCall(MatMFFDSetBase(mf, base, NULL));
2123:   PetscCall(MatComputeOperator(mf, MATAIJ, &Dmf));
2124:   PetscCall(MatTranspose(Dmf, MAT_INPLACE_MATRIX, &Dmf));
2125:   PetscCall(MatComputeOperatorTranspose(mat, MATAIJ, &Dmat));

2127:   PetscCall(MatDuplicate(Dmat, MAT_COPY_VALUES, &Ddiff));
2128:   PetscCall(MatAXPY(Ddiff, -1.0, Dmf, DIFFERENT_NONZERO_PATTERN));
2129:   PetscCall(MatNorm(Ddiff, NORM_FROBENIUS, &Diffnorm));
2130:   PetscCall(MatNorm(Dmf, NORM_FROBENIUS, &Dmfnorm));
2131:   if (Diffnorm / Dmfnorm > 10 * PETSC_SQRT_MACHINE_EPSILON) {
2132:     flag = PETSC_FALSE;
2133:     if (v) {
2134:       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)));
2135:       PetscCall(MatViewFromOptions(Ddiff, (PetscObject)mat, "-mat_shell_test_mult_transpose_view"));
2136:       PetscCall(MatViewFromOptions(Dmf, (PetscObject)mat, "-mat_shell_test_mult_transpose_view"));
2137:       PetscCall(MatViewFromOptions(Dmat, (PetscObject)mat, "-mat_shell_test_mult_transpose_view"));
2138:     }
2139:   } else if (v) {
2140:     PetscCall(PetscPrintf(PetscObjectComm((PetscObject)mat), "MATSHELL transpose and matrix-free multiple appear to produce the same results\n"));
2141:   }
2142:   if (flg) *flg = flag;
2143:   PetscCall(MatDestroy(&mf));
2144:   PetscCall(MatDestroy(&Dmat));
2145:   PetscCall(MatDestroy(&Ddiff));
2146:   PetscCall(MatDestroy(&Dmf));
2147:   PetscCall(VecDestroy(&x));
2148:   PetscCall(VecDestroy(&y));
2149:   PetscCall(VecDestroy(&z));
2150:   PetscFunctionReturn(PETSC_SUCCESS);
2151: }

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

2156:   Logically Collective

2158:   Input Parameters:
2159: + mat - the `MATSHELL` shell matrix
2160: . op  - the name of the operation
2161: - g   - the function that provides the operation.

2163:   Level: advanced

2165:   Example Usage:
2166: .vb
2167:   extern PetscErrorCode usermult(Mat, Vec, Vec);

2169:   MatCreateShell(comm, m, n, M, N, ctx, &A);
2170:   MatShellSetOperation(A, MATOP_MULT, (void(*)(void))usermult);
2171: .ve

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

2179:   All user-provided functions (except for `MATOP_DESTROY`) should have the same calling
2180:   sequence as the usual matrix interface routines, since they
2181:   are intended to be accessed via the usual matrix interface
2182:   routines, e.g.,
2183: .vb
2184:   MatMult(Mat, Vec, Vec) -> usermult(Mat, Vec, Vec)
2185: .ve

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

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

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

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

2201: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellSetContext()`, `MatSetOperation()`, `MatShellSetManageScalingShifts()`, `MatShellSetMatProductOperation()`
2202: @*/
2203: PetscErrorCode MatShellSetOperation(Mat mat, MatOperation op, void (*g)(void))
2204: {
2205:   PetscFunctionBegin;
2207:   PetscTryMethod(mat, "MatShellSetOperation_C", (Mat, MatOperation, void (*)(void)), (mat, op, g));
2208:   PetscFunctionReturn(PETSC_SUCCESS);
2209: }

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

2214:   Not Collective

2216:   Input Parameters:
2217: + mat - the `MATSHELL` shell matrix
2218: - op  - the name of the operation

2220:   Output Parameter:
2221: . g - the function that provides the operation.

2223:   Level: advanced

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

2231:   All user-provided functions have the same calling
2232:   sequence as the usual matrix interface routines, since they
2233:   are intended to be accessed via the usual matrix interface
2234:   routines, e.g.,
2235: .vb
2236:   MatMult(Mat, Vec, Vec) -> usermult(Mat, Vec, Vec)
2237: .ve

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

2243: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellSetOperation()`, `MatShellSetContext()`
2244: @*/
2245: PetscErrorCode MatShellGetOperation(Mat mat, MatOperation op, void (**g)(void))
2246: {
2247:   PetscFunctionBegin;
2249:   PetscUseMethod(mat, "MatShellGetOperation_C", (Mat, MatOperation, void (**)(void)), (mat, op, g));
2250:   PetscFunctionReturn(PETSC_SUCCESS);
2251: }

2253: /*@
2254:   MatIsShell - Inquires if a matrix is derived from `MATSHELL`

2256:   Input Parameter:
2257: . mat - the matrix

2259:   Output Parameter:
2260: . flg - the Boolean value

2262:   Level: developer

2264: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MATMFFD`, `MatCreateShell()`, `MATTRANSPOSEVIRTUAL`, `MATSCHURCOMPLEMENT`
2265: @*/
2266: PetscErrorCode MatIsShell(Mat mat, PetscBool *flg)
2267: {
2268:   PetscFunctionBegin;
2270:   PetscAssertPointer(flg, 2);
2271:   *flg = (PetscBool)(mat->ops->destroy == MatDestroy_Shell);
2272:   PetscFunctionReturn(PETSC_SUCCESS);
2273: }