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 <petsc/private/matimpl.h>
  8: #include <petsc/private/vecimpl.h>

 10: struct _MatShellOps {
 11:   /*  3 */ PetscErrorCode (*mult)(Mat, Vec, Vec);
 12:   /*  5 */ PetscErrorCode (*multtranspose)(Mat, Vec, Vec);
 13:   /* 17 */ PetscErrorCode (*getdiagonal)(Mat, Vec);
 14:   /* 43 */ PetscErrorCode (*copy)(Mat, Mat, MatStructure);
 15:   /* 60 */ PetscErrorCode (*destroy)(Mat);
 16: };

 18: struct _n_MatShellMatFunctionList {
 19:   PetscErrorCode (*symbolic)(Mat, Mat, Mat, void **);
 20:   PetscErrorCode (*numeric)(Mat, Mat, Mat, void *);
 21:   PetscErrorCode (*destroy)(void *);
 22:   MatProductType ptype;
 23:   char          *composedname; /* string to identify routine with double dispatch */
 24:   char          *resultname;   /* result matrix type */

 26:   struct _n_MatShellMatFunctionList *next;
 27: };
 28: typedef struct _n_MatShellMatFunctionList *MatShellMatFunctionList;

 30: typedef struct {
 31:   struct _MatShellOps ops[1];

 33:   /* The user will manage the scaling and shifts for the MATSHELL, not the default */
 34:   PetscBool managescalingshifts;

 36:   /* support for MatScale, MatShift and MatMultAdd */
 37:   PetscScalar vscale, vshift;
 38:   Vec         dshift;
 39:   Vec         left, right;
 40:   Vec         left_work, right_work;
 41:   Vec         left_add_work, right_add_work;

 43:   /* support for MatAXPY */
 44:   Mat              axpy;
 45:   PetscScalar      axpy_vscale;
 46:   Vec              axpy_left, axpy_right;
 47:   PetscObjectState axpy_state;

 49:   /* support for ZeroRows/Columns operations */
 50:   IS         zrows;
 51:   IS         zcols;
 52:   Vec        zvals;
 53:   Vec        zvals_w;
 54:   VecScatter zvals_sct_r;
 55:   VecScatter zvals_sct_c;

 57:   /* MatMat operations */
 58:   MatShellMatFunctionList matmat;

 60:   /* user defined context */
 61:   PetscContainer ctxcontainer;
 62: } Mat_Shell;

 64: /*
 65:      Store and scale values on zeroed rows
 66:      xx = [x_1, 0], 0 on zeroed columns
 67: */
 68: static PetscErrorCode MatShellPreZeroRight(Mat A, Vec x, Vec *xx)
 69: {
 70:   Mat_Shell *shell = (Mat_Shell *)A->data;

 72:   PetscFunctionBegin;
 73:   *xx = x;
 74:   if (shell->zrows) {
 75:     PetscCall(VecSet(shell->zvals_w, 0.0));
 76:     PetscCall(VecScatterBegin(shell->zvals_sct_c, x, shell->zvals_w, INSERT_VALUES, SCATTER_FORWARD));
 77:     PetscCall(VecScatterEnd(shell->zvals_sct_c, x, shell->zvals_w, INSERT_VALUES, SCATTER_FORWARD));
 78:     PetscCall(VecPointwiseMult(shell->zvals_w, shell->zvals_w, shell->zvals));
 79:   }
 80:   if (shell->zcols) {
 81:     if (!shell->right_work) PetscCall(MatCreateVecs(A, &shell->right_work, NULL));
 82:     PetscCall(VecCopy(x, shell->right_work));
 83:     PetscCall(VecISSet(shell->right_work, shell->zcols, 0.0));
 84:     *xx = shell->right_work;
 85:   }
 86:   PetscFunctionReturn(PETSC_SUCCESS);
 87: }

 89: /* Insert properly diagonally scaled values stored in MatShellPreZeroRight */
 90: static PetscErrorCode MatShellPostZeroLeft(Mat A, Vec x)
 91: {
 92:   Mat_Shell *shell = (Mat_Shell *)A->data;

 94:   PetscFunctionBegin;
 95:   if (shell->zrows) {
 96:     PetscCall(VecScatterBegin(shell->zvals_sct_r, shell->zvals_w, x, INSERT_VALUES, SCATTER_REVERSE));
 97:     PetscCall(VecScatterEnd(shell->zvals_sct_r, shell->zvals_w, x, INSERT_VALUES, SCATTER_REVERSE));
 98:   }
 99:   PetscFunctionReturn(PETSC_SUCCESS);
100: }

102: /*
103:      Store and scale values on zeroed rows
104:      xx = [x_1, 0], 0 on zeroed rows
105: */
106: static PetscErrorCode MatShellPreZeroLeft(Mat A, Vec x, Vec *xx)
107: {
108:   Mat_Shell *shell = (Mat_Shell *)A->data;

110:   PetscFunctionBegin;
111:   *xx = NULL;
112:   if (!shell->zrows) {
113:     *xx = x;
114:   } else {
115:     if (!shell->left_work) PetscCall(MatCreateVecs(A, NULL, &shell->left_work));
116:     PetscCall(VecCopy(x, shell->left_work));
117:     PetscCall(VecSet(shell->zvals_w, 0.0));
118:     PetscCall(VecScatterBegin(shell->zvals_sct_r, shell->zvals_w, shell->left_work, INSERT_VALUES, SCATTER_REVERSE));
119:     PetscCall(VecScatterEnd(shell->zvals_sct_r, shell->zvals_w, shell->left_work, INSERT_VALUES, SCATTER_REVERSE));
120:     PetscCall(VecScatterBegin(shell->zvals_sct_r, x, shell->zvals_w, INSERT_VALUES, SCATTER_FORWARD));
121:     PetscCall(VecScatterEnd(shell->zvals_sct_r, x, shell->zvals_w, INSERT_VALUES, SCATTER_FORWARD));
122:     PetscCall(VecPointwiseMult(shell->zvals_w, shell->zvals_w, shell->zvals));
123:     *xx = shell->left_work;
124:   }
125:   PetscFunctionReturn(PETSC_SUCCESS);
126: }

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

133:   PetscFunctionBegin;
134:   if (shell->zcols) PetscCall(VecISSet(x, shell->zcols, 0.0));
135:   if (shell->zrows) {
136:     PetscCall(VecScatterBegin(shell->zvals_sct_c, shell->zvals_w, x, ADD_VALUES, SCATTER_REVERSE));
137:     PetscCall(VecScatterEnd(shell->zvals_sct_c, shell->zvals_w, x, ADD_VALUES, SCATTER_REVERSE));
138:   }
139:   PetscFunctionReturn(PETSC_SUCCESS);
140: }

142: /*
143:       xx = diag(left)*x
144: */
145: static PetscErrorCode MatShellPreScaleLeft(Mat A, Vec x, Vec *xx)
146: {
147:   Mat_Shell *shell = (Mat_Shell *)A->data;

149:   PetscFunctionBegin;
150:   *xx = NULL;
151:   if (!shell->left) {
152:     *xx = x;
153:   } else {
154:     if (!shell->left_work) PetscCall(VecDuplicate(shell->left, &shell->left_work));
155:     PetscCall(VecPointwiseMult(shell->left_work, x, shell->left));
156:     *xx = shell->left_work;
157:   }
158:   PetscFunctionReturn(PETSC_SUCCESS);
159: }

161: /*
162:      xx = diag(right)*x
163: */
164: static PetscErrorCode MatShellPreScaleRight(Mat A, Vec x, Vec *xx)
165: {
166:   Mat_Shell *shell = (Mat_Shell *)A->data;

168:   PetscFunctionBegin;
169:   *xx = NULL;
170:   if (!shell->right) {
171:     *xx = x;
172:   } else {
173:     if (!shell->right_work) PetscCall(VecDuplicate(shell->right, &shell->right_work));
174:     PetscCall(VecPointwiseMult(shell->right_work, x, shell->right));
175:     *xx = shell->right_work;
176:   }
177:   PetscFunctionReturn(PETSC_SUCCESS);
178: }

180: /*
181:     x = diag(left)*x
182: */
183: static PetscErrorCode MatShellPostScaleLeft(Mat A, Vec x)
184: {
185:   Mat_Shell *shell = (Mat_Shell *)A->data;

187:   PetscFunctionBegin;
188:   if (shell->left) PetscCall(VecPointwiseMult(x, x, shell->left));
189:   PetscFunctionReturn(PETSC_SUCCESS);
190: }

192: /*
193:     x = diag(right)*x
194: */
195: static PetscErrorCode MatShellPostScaleRight(Mat A, Vec x)
196: {
197:   Mat_Shell *shell = (Mat_Shell *)A->data;

199:   PetscFunctionBegin;
200:   if (shell->right) PetscCall(VecPointwiseMult(x, x, shell->right));
201:   PetscFunctionReturn(PETSC_SUCCESS);
202: }

204: /*
205:          Y = vscale*Y + diag(dshift)*X + vshift*X

207:          On input Y already contains A*x
208: */
209: static PetscErrorCode MatShellShiftAndScale(Mat A, Vec X, Vec Y)
210: {
211:   Mat_Shell *shell = (Mat_Shell *)A->data;

213:   PetscFunctionBegin;
214:   if (shell->dshift) { /* get arrays because there is no VecPointwiseMultAdd() */
215:     PetscInt           i, m;
216:     const PetscScalar *x, *d;
217:     PetscScalar       *y;
218:     PetscCall(VecGetLocalSize(X, &m));
219:     PetscCall(VecGetArrayRead(shell->dshift, &d));
220:     PetscCall(VecGetArrayRead(X, &x));
221:     PetscCall(VecGetArray(Y, &y));
222:     for (i = 0; i < m; i++) y[i] = shell->vscale * y[i] + d[i] * x[i];
223:     PetscCall(VecRestoreArrayRead(shell->dshift, &d));
224:     PetscCall(VecRestoreArrayRead(X, &x));
225:     PetscCall(VecRestoreArray(Y, &y));
226:   } else {
227:     PetscCall(VecScale(Y, shell->vscale));
228:   }
229:   if (shell->vshift != 0.0) PetscCall(VecAXPY(Y, shell->vshift, X)); /* if test is for non-square matrices */
230:   PetscFunctionReturn(PETSC_SUCCESS);
231: }

233: static PetscErrorCode MatShellGetContext_Shell(Mat mat, void *ctx)
234: {
235:   Mat_Shell *shell = (Mat_Shell *)mat->data;

237:   PetscFunctionBegin;
238:   if (shell->ctxcontainer) PetscCall(PetscContainerGetPointer(shell->ctxcontainer, (void **)ctx));
239:   else *(void **)ctx = NULL;
240:   PetscFunctionReturn(PETSC_SUCCESS);
241: }

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

246:   Not Collective

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

251:   Output Parameter:
252: . ctx - the user provided context

254:   Level: advanced

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

260: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellSetOperation()`, `MatShellSetContext()`
261: @*/
262: PetscErrorCode MatShellGetContext(Mat mat, void *ctx)
263: {
264:   PetscFunctionBegin;
266:   PetscAssertPointer(ctx, 2);
267:   PetscUseMethod(mat, "MatShellGetContext_C", (Mat, void *), (mat, ctx));
268:   PetscFunctionReturn(PETSC_SUCCESS);
269: }

271: static PetscErrorCode MatZeroRowsColumns_Local_Shell(Mat mat, PetscInt nr, PetscInt rows[], PetscInt nc, PetscInt cols[], PetscScalar diag, PetscBool rc)
272: {
273:   Mat_Shell      *shell = (Mat_Shell *)mat->data;
274:   Vec             x = NULL, b = NULL;
275:   IS              is1, is2;
276:   const PetscInt *ridxs;
277:   PetscInt       *idxs, *gidxs;
278:   PetscInt        cum, rst, cst, i;

280:   PetscFunctionBegin;
281:   if (!shell->zvals) PetscCall(MatCreateVecs(mat, NULL, &shell->zvals));
282:   if (!shell->zvals_w) PetscCall(VecDuplicate(shell->zvals, &shell->zvals_w));
283:   PetscCall(MatGetOwnershipRange(mat, &rst, NULL));
284:   PetscCall(MatGetOwnershipRangeColumn(mat, &cst, NULL));

286:   /* Expand/create index set of zeroed rows */
287:   PetscCall(PetscMalloc1(nr, &idxs));
288:   for (i = 0; i < nr; i++) idxs[i] = rows[i] + rst;
289:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nr, idxs, PETSC_OWN_POINTER, &is1));
290:   PetscCall(ISSort(is1));
291:   PetscCall(VecISSet(shell->zvals, is1, diag));
292:   if (shell->zrows) {
293:     PetscCall(ISSum(shell->zrows, is1, &is2));
294:     PetscCall(ISDestroy(&shell->zrows));
295:     PetscCall(ISDestroy(&is1));
296:     shell->zrows = is2;
297:   } else shell->zrows = is1;

299:   /* Create scatters for diagonal values communications */
300:   PetscCall(VecScatterDestroy(&shell->zvals_sct_c));
301:   PetscCall(VecScatterDestroy(&shell->zvals_sct_r));

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

307:   /* col scatter: from right vector to left vector */
308:   PetscCall(ISGetIndices(shell->zrows, &ridxs));
309:   PetscCall(ISGetLocalSize(shell->zrows, &nr));
310:   PetscCall(PetscMalloc1(nr, &gidxs));
311:   for (i = 0, cum = 0; i < nr; i++) {
312:     if (ridxs[i] >= mat->cmap->N) continue;
313:     gidxs[cum] = ridxs[i];
314:     cum++;
315:   }
316:   PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)mat), cum, gidxs, PETSC_OWN_POINTER, &is1));
317:   PetscCall(VecScatterCreate(x, is1, shell->zvals_w, is1, &shell->zvals_sct_c));
318:   PetscCall(ISDestroy(&is1));
319:   PetscCall(VecDestroy(&x));
320:   PetscCall(VecDestroy(&b));

322:   /* Expand/create index set of zeroed columns */
323:   if (rc) {
324:     PetscCall(PetscMalloc1(nc, &idxs));
325:     for (i = 0; i < nc; i++) idxs[i] = cols[i] + cst;
326:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nc, idxs, PETSC_OWN_POINTER, &is1));
327:     PetscCall(ISSort(is1));
328:     if (shell->zcols) {
329:       PetscCall(ISSum(shell->zcols, is1, &is2));
330:       PetscCall(ISDestroy(&shell->zcols));
331:       PetscCall(ISDestroy(&is1));
332:       shell->zcols = is2;
333:     } else shell->zcols = is1;
334:   }
335:   PetscFunctionReturn(PETSC_SUCCESS);
336: }

338: static PetscErrorCode MatZeroRows_Shell(Mat mat, PetscInt n, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
339: {
340:   Mat_Shell *shell = (Mat_Shell *)mat->data;
341:   PetscInt   nr, *lrows;

343:   PetscFunctionBegin;
344:   if (x && b) {
345:     Vec          xt;
346:     PetscScalar *vals;
347:     PetscInt    *gcols, i, st, nl, nc;

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

353:     PetscCall(MatCreateVecs(mat, &xt, NULL));
354:     PetscCall(VecCopy(x, xt));
355:     PetscCall(PetscCalloc1(nc, &vals));
356:     PetscCall(VecSetValues(xt, nc, gcols, vals, INSERT_VALUES)); /* xt = [x1, 0] */
357:     PetscCall(PetscFree(vals));
358:     PetscCall(VecAssemblyBegin(xt));
359:     PetscCall(VecAssemblyEnd(xt));
360:     PetscCall(VecAYPX(xt, -1.0, x)); /* xt = [0, x2] */

362:     PetscCall(VecGetOwnershipRange(xt, &st, NULL));
363:     PetscCall(VecGetLocalSize(xt, &nl));
364:     PetscCall(VecGetArray(xt, &vals));
365:     for (i = 0; i < nl; i++) {
366:       PetscInt g = i + st;
367:       if (g > mat->rmap->N) continue;
368:       if (PetscAbsScalar(vals[i]) == 0.0) continue;
369:       PetscCall(VecSetValue(b, g, diag * vals[i], INSERT_VALUES));
370:     }
371:     PetscCall(VecRestoreArray(xt, &vals));
372:     PetscCall(VecAssemblyBegin(b));
373:     PetscCall(VecAssemblyEnd(b)); /* b  = [b1, x2 * diag] */
374:     PetscCall(VecDestroy(&xt));
375:     PetscCall(PetscFree(gcols));
376:   }
377:   PetscCall(PetscLayoutMapLocal(mat->rmap, n, rows, &nr, &lrows, NULL));
378:   PetscCall(MatZeroRowsColumns_Local_Shell(mat, nr, lrows, 0, NULL, diag, PETSC_FALSE));
379:   if (shell->axpy) PetscCall(MatZeroRows(shell->axpy, n, rows, 0.0, NULL, NULL));
380:   PetscCall(PetscFree(lrows));
381:   PetscFunctionReturn(PETSC_SUCCESS);
382: }

384: static PetscErrorCode MatZeroRowsColumns_Shell(Mat mat, PetscInt n, const PetscInt rowscols[], PetscScalar diag, Vec x, Vec b)
385: {
386:   Mat_Shell *shell = (Mat_Shell *)mat->data;
387:   PetscInt  *lrows, *lcols;
388:   PetscInt   nr, nc;
389:   PetscBool  congruent;

391:   PetscFunctionBegin;
392:   if (x && b) {
393:     Vec          xt, bt;
394:     PetscScalar *vals;
395:     PetscInt    *grows, *gcols, i, st, nl;

397:     PetscCall(PetscMalloc2(n, &grows, n, &gcols));
398:     for (i = 0, nr = 0; i < n; i++)
399:       if (rowscols[i] < mat->rmap->N) grows[nr++] = rowscols[i];
400:     for (i = 0, nc = 0; i < n; i++)
401:       if (rowscols[i] < mat->cmap->N) gcols[nc++] = rowscols[i];
402:     PetscCall(PetscCalloc1(n, &vals));

404:     PetscCall(MatCreateVecs(mat, &xt, &bt));
405:     PetscCall(VecCopy(x, xt));
406:     PetscCall(VecSetValues(xt, nc, gcols, vals, INSERT_VALUES)); /* xt = [x1, 0] */
407:     PetscCall(VecAssemblyBegin(xt));
408:     PetscCall(VecAssemblyEnd(xt));
409:     PetscCall(VecAXPY(xt, -1.0, x));                             /* xt = [0, -x2] */
410:     PetscCall(MatMult(mat, xt, bt));                             /* bt = [-A12*x2,-A22*x2] */
411:     PetscCall(VecSetValues(bt, nr, grows, vals, INSERT_VALUES)); /* bt = [-A12*x2,0] */
412:     PetscCall(VecAssemblyBegin(bt));
413:     PetscCall(VecAssemblyEnd(bt));
414:     PetscCall(VecAXPY(b, 1.0, bt));                              /* b  = [b1 - A12*x2, b2] */
415:     PetscCall(VecSetValues(bt, nr, grows, vals, INSERT_VALUES)); /* b  = [b1 - A12*x2, 0] */
416:     PetscCall(VecAssemblyBegin(bt));
417:     PetscCall(VecAssemblyEnd(bt));
418:     PetscCall(PetscFree(vals));

420:     PetscCall(VecGetOwnershipRange(xt, &st, NULL));
421:     PetscCall(VecGetLocalSize(xt, &nl));
422:     PetscCall(VecGetArray(xt, &vals));
423:     for (i = 0; i < nl; i++) {
424:       PetscInt g = i + st;
425:       if (g > mat->rmap->N) continue;
426:       if (PetscAbsScalar(vals[i]) == 0.0) continue;
427:       PetscCall(VecSetValue(b, g, -diag * vals[i], INSERT_VALUES));
428:     }
429:     PetscCall(VecRestoreArray(xt, &vals));
430:     PetscCall(VecAssemblyBegin(b));
431:     PetscCall(VecAssemblyEnd(b)); /* b  = [b1 - A12*x2, x2 * diag] */
432:     PetscCall(VecDestroy(&xt));
433:     PetscCall(VecDestroy(&bt));
434:     PetscCall(PetscFree2(grows, gcols));
435:   }
436:   PetscCall(PetscLayoutMapLocal(mat->rmap, n, rowscols, &nr, &lrows, NULL));
437:   PetscCall(MatHasCongruentLayouts(mat, &congruent));
438:   if (congruent) {
439:     nc    = nr;
440:     lcols = lrows;
441:   } else { /* MatZeroRowsColumns implicitly assumes the rowscols indices are for a square matrix, here we handle a more general case */
442:     PetscInt i, nt, *t;

444:     PetscCall(PetscMalloc1(n, &t));
445:     for (i = 0, nt = 0; i < n; i++)
446:       if (rowscols[i] < mat->cmap->N) t[nt++] = rowscols[i];
447:     PetscCall(PetscLayoutMapLocal(mat->cmap, nt, t, &nc, &lcols, NULL));
448:     PetscCall(PetscFree(t));
449:   }
450:   PetscCall(MatZeroRowsColumns_Local_Shell(mat, nr, lrows, nc, lcols, diag, PETSC_TRUE));
451:   if (!congruent) PetscCall(PetscFree(lcols));
452:   PetscCall(PetscFree(lrows));
453:   if (shell->axpy) PetscCall(MatZeroRowsColumns(shell->axpy, n, rowscols, 0.0, NULL, NULL));
454:   PetscFunctionReturn(PETSC_SUCCESS);
455: }

457: static PetscErrorCode MatDestroy_Shell(Mat mat)
458: {
459:   Mat_Shell              *shell = (Mat_Shell *)mat->data;
460:   MatShellMatFunctionList matmat;

462:   PetscFunctionBegin;
463:   if (shell->ops->destroy) PetscCall((*shell->ops->destroy)(mat));
464:   PetscCall(PetscMemzero(shell->ops, sizeof(struct _MatShellOps)));
465:   PetscCall(VecDestroy(&shell->left));
466:   PetscCall(VecDestroy(&shell->right));
467:   PetscCall(VecDestroy(&shell->dshift));
468:   PetscCall(VecDestroy(&shell->left_work));
469:   PetscCall(VecDestroy(&shell->right_work));
470:   PetscCall(VecDestroy(&shell->left_add_work));
471:   PetscCall(VecDestroy(&shell->right_add_work));
472:   PetscCall(VecDestroy(&shell->axpy_left));
473:   PetscCall(VecDestroy(&shell->axpy_right));
474:   PetscCall(MatDestroy(&shell->axpy));
475:   PetscCall(VecDestroy(&shell->zvals_w));
476:   PetscCall(VecDestroy(&shell->zvals));
477:   PetscCall(VecScatterDestroy(&shell->zvals_sct_c));
478:   PetscCall(VecScatterDestroy(&shell->zvals_sct_r));
479:   PetscCall(ISDestroy(&shell->zrows));
480:   PetscCall(ISDestroy(&shell->zcols));

482:   matmat = shell->matmat;
483:   while (matmat) {
484:     MatShellMatFunctionList next = matmat->next;

486:     PetscCall(PetscObjectComposeFunction((PetscObject)mat, matmat->composedname, NULL));
487:     PetscCall(PetscFree(matmat->composedname));
488:     PetscCall(PetscFree(matmat->resultname));
489:     PetscCall(PetscFree(matmat));
490:     matmat = next;
491:   }
492:   PetscCall(MatShellSetContext(mat, NULL));
493:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatShellGetContext_C", NULL));
494:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatShellSetContext_C", NULL));
495:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatShellSetContextDestroy_C", NULL));
496:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatShellSetVecType_C", NULL));
497:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatShellSetManageScalingShifts_C", NULL));
498:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatShellSetOperation_C", NULL));
499:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatShellGetOperation_C", NULL));
500:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatShellSetMatProductOperation_C", NULL));
501:   PetscCall(PetscFree(mat->data));
502:   PetscFunctionReturn(PETSC_SUCCESS);
503: }

505: typedef struct {
506:   PetscErrorCode (*numeric)(Mat, Mat, Mat, void *);
507:   PetscErrorCode (*destroy)(void *);
508:   void *userdata;
509:   Mat   B;
510:   Mat   Bt;
511:   Mat   axpy;
512: } MatMatDataShell;

514: static PetscErrorCode DestroyMatMatDataShell(void *data)
515: {
516:   MatMatDataShell *mmdata = (MatMatDataShell *)data;

518:   PetscFunctionBegin;
519:   if (mmdata->destroy) PetscCall((*mmdata->destroy)(mmdata->userdata));
520:   PetscCall(MatDestroy(&mmdata->B));
521:   PetscCall(MatDestroy(&mmdata->Bt));
522:   PetscCall(MatDestroy(&mmdata->axpy));
523:   PetscCall(PetscFree(mmdata));
524:   PetscFunctionReturn(PETSC_SUCCESS);
525: }

527: static PetscErrorCode MatProductNumeric_Shell_X(Mat D)
528: {
529:   Mat_Product     *product;
530:   Mat              A, B;
531:   MatMatDataShell *mdata;
532:   PetscScalar      zero = 0.0;

534:   PetscFunctionBegin;
535:   MatCheckProduct(D, 1);
536:   product = D->product;
537:   PetscCheck(product->data, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "Product data empty");
538:   A     = product->A;
539:   B     = product->B;
540:   mdata = (MatMatDataShell *)product->data;
541:   if (mdata->numeric) {
542:     Mat_Shell *shell                = (Mat_Shell *)A->data;
543:     PetscErrorCode (*stashsym)(Mat) = D->ops->productsymbolic;
544:     PetscErrorCode (*stashnum)(Mat) = D->ops->productnumeric;
545:     PetscBool useBmdata = PETSC_FALSE, newB = PETSC_TRUE;

547:     if (shell->managescalingshifts) {
548:       PetscCheck(!shell->zcols && !shell->zrows, PetscObjectComm((PetscObject)D), PETSC_ERR_SUP, "MatProduct not supported with zeroed rows/columns");
549:       if (shell->right || shell->left) {
550:         useBmdata = PETSC_TRUE;
551:         if (!mdata->B) {
552:           PetscCall(MatDuplicate(B, MAT_SHARE_NONZERO_PATTERN, &mdata->B));
553:         } else {
554:           newB = PETSC_FALSE;
555:         }
556:         PetscCall(MatCopy(B, mdata->B, SAME_NONZERO_PATTERN));
557:       }
558:       switch (product->type) {
559:       case MATPRODUCT_AB: /* s L A R B + v L R B + L D R B */
560:         if (shell->right) PetscCall(MatDiagonalScale(mdata->B, shell->right, NULL));
561:         break;
562:       case MATPRODUCT_AtB: /* s R A^t L B + v R L B + R D L B */
563:         if (shell->left) PetscCall(MatDiagonalScale(mdata->B, shell->left, NULL));
564:         break;
565:       case MATPRODUCT_ABt: /* s L A R B^t + v L R B^t + L D R B^t */
566:         if (shell->right) PetscCall(MatDiagonalScale(mdata->B, NULL, shell->right));
567:         break;
568:       case MATPRODUCT_RARt: /* s B L A R B^t + v B L R B^t + B L D R B^t */
569:         if (shell->right && shell->left) {
570:           PetscBool flg;

572:           PetscCall(VecEqual(shell->right, shell->left, &flg));
573:           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,
574:                      ((PetscObject)B)->type_name);
575:         }
576:         if (shell->right) PetscCall(MatDiagonalScale(mdata->B, NULL, shell->right));
577:         break;
578:       case MATPRODUCT_PtAP: /* s B^t L A R B + v B^t L R B + B^t L D R B */
579:         if (shell->right && shell->left) {
580:           PetscBool flg;

582:           PetscCall(VecEqual(shell->right, shell->left, &flg));
583:           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,
584:                      ((PetscObject)B)->type_name);
585:         }
586:         if (shell->right) PetscCall(MatDiagonalScale(mdata->B, shell->right, NULL));
587:         break;
588:       default:
589:         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);
590:       }
591:     }
592:     /* allow the user to call MatMat operations on D */
593:     D->product              = NULL;
594:     D->ops->productsymbolic = NULL;
595:     D->ops->productnumeric  = NULL;

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

599:     /* clear any leftover user data and restore D pointers */
600:     PetscCall(MatProductClear(D));
601:     D->ops->productsymbolic = stashsym;
602:     D->ops->productnumeric  = stashnum;
603:     D->product              = product;

605:     if (shell->managescalingshifts) {
606:       PetscCall(MatScale(D, shell->vscale));
607:       switch (product->type) {
608:       case MATPRODUCT_AB:  /* s L A R B + v L R B + L D R B */
609:       case MATPRODUCT_ABt: /* s L A R B^t + v L R B^t + L D R B^t */
610:         if (shell->left) {
611:           PetscCall(MatDiagonalScale(D, shell->left, NULL));
612:           if (shell->dshift || shell->vshift != zero) {
613:             if (!shell->left_work) PetscCall(MatCreateVecs(A, NULL, &shell->left_work));
614:             if (shell->dshift) {
615:               PetscCall(VecCopy(shell->dshift, shell->left_work));
616:               PetscCall(VecShift(shell->left_work, shell->vshift));
617:               PetscCall(VecPointwiseMult(shell->left_work, shell->left_work, shell->left));
618:             } else {
619:               PetscCall(VecSet(shell->left_work, shell->vshift));
620:             }
621:             if (product->type == MATPRODUCT_ABt) {
622:               MatReuse     reuse = mdata->Bt ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX;
623:               MatStructure str   = mdata->Bt ? SUBSET_NONZERO_PATTERN : DIFFERENT_NONZERO_PATTERN;

625:               PetscCall(MatTranspose(mdata->B, reuse, &mdata->Bt));
626:               PetscCall(MatDiagonalScale(mdata->Bt, shell->left_work, NULL));
627:               PetscCall(MatAXPY(D, 1.0, mdata->Bt, str));
628:             } else {
629:               MatStructure str = newB ? DIFFERENT_NONZERO_PATTERN : SUBSET_NONZERO_PATTERN;

631:               PetscCall(MatDiagonalScale(mdata->B, shell->left_work, NULL));
632:               PetscCall(MatAXPY(D, 1.0, mdata->B, str));
633:             }
634:           }
635:         }
636:         break;
637:       case MATPRODUCT_AtB: /* s R A^t L B + v R L B + R D L B */
638:         if (shell->right) {
639:           PetscCall(MatDiagonalScale(D, shell->right, NULL));
640:           if (shell->dshift || shell->vshift != zero) {
641:             MatStructure str = newB ? DIFFERENT_NONZERO_PATTERN : SUBSET_NONZERO_PATTERN;

643:             if (!shell->right_work) PetscCall(MatCreateVecs(A, &shell->right_work, NULL));
644:             if (shell->dshift) {
645:               PetscCall(VecCopy(shell->dshift, shell->right_work));
646:               PetscCall(VecShift(shell->right_work, shell->vshift));
647:               PetscCall(VecPointwiseMult(shell->right_work, shell->right_work, shell->right));
648:             } else {
649:               PetscCall(VecSet(shell->right_work, shell->vshift));
650:             }
651:             PetscCall(MatDiagonalScale(mdata->B, shell->right_work, NULL));
652:             PetscCall(MatAXPY(D, 1.0, mdata->B, str));
653:           }
654:         }
655:         break;
656:       case MATPRODUCT_PtAP: /* s B^t L A R B + v B^t L R B + B^t L D R B */
657:       case MATPRODUCT_RARt: /* s B L A R B^t + v B L R B^t + B L D R B^t */
658:         PetscCheck(!shell->dshift && shell->vshift == zero, PetscObjectComm((PetscObject)D), PETSC_ERR_SUP, "MatProductSymbolic type %s not supported for %s and %s matrices with diagonal shift", MatProductTypes[product->type], ((PetscObject)A)->type_name,
659:                    ((PetscObject)B)->type_name);
660:         break;
661:       default:
662:         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);
663:       }
664:       if (shell->axpy && shell->axpy_vscale != zero) {
665:         Mat              X;
666:         PetscObjectState axpy_state;
667:         MatStructure     str = DIFFERENT_NONZERO_PATTERN; /* not sure it is safe to ever use SUBSET_NONZERO_PATTERN */

669:         PetscCall(MatShellGetContext(shell->axpy, &X));
670:         PetscCall(PetscObjectStateGet((PetscObject)X, &axpy_state));
671:         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,...)");
672:         if (!mdata->axpy) {
673:           str = DIFFERENT_NONZERO_PATTERN;
674:           PetscCall(MatProductCreate(shell->axpy, B, NULL, &mdata->axpy));
675:           PetscCall(MatProductSetType(mdata->axpy, product->type));
676:           PetscCall(MatProductSetFromOptions(mdata->axpy));
677:           PetscCall(MatProductSymbolic(mdata->axpy));
678:         } else { /* May be that shell->axpy has changed */
679:           PetscBool flg;

681:           PetscCall(MatProductReplaceMats(shell->axpy, B, NULL, mdata->axpy));
682:           PetscCall(MatHasOperation(mdata->axpy, MATOP_PRODUCTSYMBOLIC, &flg));
683:           if (!flg) {
684:             str = DIFFERENT_NONZERO_PATTERN;
685:             PetscCall(MatProductSetFromOptions(mdata->axpy));
686:             PetscCall(MatProductSymbolic(mdata->axpy));
687:           }
688:         }
689:         PetscCall(MatProductNumeric(mdata->axpy));
690:         PetscCall(MatAXPY(D, shell->axpy_vscale, mdata->axpy, str));
691:       }
692:     }
693:   } else SETERRQ(PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "Missing numeric operation");
694:   PetscFunctionReturn(PETSC_SUCCESS);
695: }

697: static PetscErrorCode MatProductSymbolic_Shell_X(Mat D)
698: {
699:   Mat_Product            *product;
700:   Mat                     A, B;
701:   MatShellMatFunctionList matmat;
702:   Mat_Shell              *shell;
703:   PetscBool               flg = PETSC_FALSE;
704:   char                    composedname[256];
705:   MatMatDataShell        *mdata;

707:   PetscFunctionBegin;
708:   MatCheckProduct(D, 1);
709:   product = D->product;
710:   PetscCheck(!product->data, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "Product data not empty");
711:   A      = product->A;
712:   B      = product->B;
713:   shell  = (Mat_Shell *)A->data;
714:   matmat = shell->matmat;
715:   PetscCall(PetscSNPrintf(composedname, sizeof(composedname), "MatProductSetFromOptions_%s_%s_C", ((PetscObject)A)->type_name, ((PetscObject)B)->type_name));
716:   while (matmat) {
717:     PetscCall(PetscStrcmp(composedname, matmat->composedname, &flg));
718:     flg = (PetscBool)(flg && (matmat->ptype == product->type));
719:     if (flg) break;
720:     matmat = matmat->next;
721:   }
722:   PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Composedname \"%s\" for product type %s not found", composedname, MatProductTypes[product->type]);
723:   switch (product->type) {
724:   case MATPRODUCT_AB:
725:     PetscCall(MatSetSizes(D, A->rmap->n, B->cmap->n, A->rmap->N, B->cmap->N));
726:     break;
727:   case MATPRODUCT_AtB:
728:     PetscCall(MatSetSizes(D, A->cmap->n, B->cmap->n, A->cmap->N, B->cmap->N));
729:     break;
730:   case MATPRODUCT_ABt:
731:     PetscCall(MatSetSizes(D, A->rmap->n, B->rmap->n, A->rmap->N, B->rmap->N));
732:     break;
733:   case MATPRODUCT_RARt:
734:     PetscCall(MatSetSizes(D, B->rmap->n, B->rmap->n, B->rmap->N, B->rmap->N));
735:     break;
736:   case MATPRODUCT_PtAP:
737:     PetscCall(MatSetSizes(D, B->cmap->n, B->cmap->n, B->cmap->N, B->cmap->N));
738:     break;
739:   default:
740:     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);
741:   }
742:   /* respect users who passed in a matrix for which resultname is the base type */
743:   if (matmat->resultname) {
744:     PetscCall(PetscObjectBaseTypeCompare((PetscObject)D, matmat->resultname, &flg));
745:     if (!flg) PetscCall(MatSetType(D, matmat->resultname));
746:   }
747:   /* If matrix type was not set or different, we need to reset this pointers */
748:   D->ops->productsymbolic = MatProductSymbolic_Shell_X;
749:   D->ops->productnumeric  = MatProductNumeric_Shell_X;
750:   /* attach product data */
751:   PetscCall(PetscNew(&mdata));
752:   mdata->numeric = matmat->numeric;
753:   mdata->destroy = matmat->destroy;
754:   if (matmat->symbolic) {
755:     PetscCall((*matmat->symbolic)(A, B, D, &mdata->userdata));
756:   } else { /* call general setup if symbolic operation not provided */
757:     PetscCall(MatSetUp(D));
758:   }
759:   PetscCheck(D->product, PetscObjectComm((PetscObject)D), PETSC_ERR_COR, "Product disappeared after user symbolic phase");
760:   PetscCheck(!D->product->data, PetscObjectComm((PetscObject)D), PETSC_ERR_COR, "Product data not empty after user symbolic phase");
761:   D->product->data    = mdata;
762:   D->product->destroy = DestroyMatMatDataShell;
763:   /* Be sure to reset these pointers if the user did something unexpected */
764:   D->ops->productsymbolic = MatProductSymbolic_Shell_X;
765:   D->ops->productnumeric  = MatProductNumeric_Shell_X;
766:   PetscFunctionReturn(PETSC_SUCCESS);
767: }

769: static PetscErrorCode MatProductSetFromOptions_Shell_X(Mat D)
770: {
771:   Mat_Product            *product;
772:   Mat                     A, B;
773:   MatShellMatFunctionList matmat;
774:   Mat_Shell              *shell;
775:   PetscBool               flg;
776:   char                    composedname[256];

778:   PetscFunctionBegin;
779:   MatCheckProduct(D, 1);
780:   product = D->product;
781:   A       = product->A;
782:   B       = product->B;
783:   PetscCall(MatIsShell(A, &flg));
784:   if (!flg) PetscFunctionReturn(PETSC_SUCCESS);
785:   shell  = (Mat_Shell *)A->data;
786:   matmat = shell->matmat;
787:   PetscCall(PetscSNPrintf(composedname, sizeof(composedname), "MatProductSetFromOptions_%s_%s_C", ((PetscObject)A)->type_name, ((PetscObject)B)->type_name));
788:   while (matmat) {
789:     PetscCall(PetscStrcmp(composedname, matmat->composedname, &flg));
790:     flg = (PetscBool)(flg && (matmat->ptype == product->type));
791:     if (flg) break;
792:     matmat = matmat->next;
793:   }
794:   if (flg) {
795:     D->ops->productsymbolic = MatProductSymbolic_Shell_X;
796:   } else PetscCall(PetscInfo(D, "  symbolic product %s not registered for product type %s\n", composedname, MatProductTypes[product->type]));
797:   PetscFunctionReturn(PETSC_SUCCESS);
798: }

800: 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)
801: {
802:   PetscBool               flg;
803:   Mat_Shell              *shell;
804:   MatShellMatFunctionList matmat;

806:   PetscFunctionBegin;
807:   PetscCheck(numeric, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing numeric routine");
808:   PetscCheck(composedname, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing composed name");

810:   /* add product callback */
811:   shell  = (Mat_Shell *)A->data;
812:   matmat = shell->matmat;
813:   if (!matmat) {
814:     PetscCall(PetscNew(&shell->matmat));
815:     matmat = shell->matmat;
816:   } else {
817:     MatShellMatFunctionList entry = matmat;
818:     while (entry) {
819:       PetscCall(PetscStrcmp(composedname, entry->composedname, &flg));
820:       flg    = (PetscBool)(flg && (entry->ptype == ptype));
821:       matmat = entry;
822:       if (flg) goto set;
823:       entry = entry->next;
824:     }
825:     PetscCall(PetscNew(&matmat->next));
826:     matmat = matmat->next;
827:   }

829: set:
830:   matmat->symbolic = symbolic;
831:   matmat->numeric  = numeric;
832:   matmat->destroy  = destroy;
833:   matmat->ptype    = ptype;
834:   PetscCall(PetscFree(matmat->composedname));
835:   PetscCall(PetscFree(matmat->resultname));
836:   PetscCall(PetscStrallocpy(composedname, &matmat->composedname));
837:   PetscCall(PetscStrallocpy(resultname, &matmat->resultname));
838:   PetscCall(PetscInfo(A, "Composing %s for product type %s with result %s\n", matmat->composedname, MatProductTypes[matmat->ptype], matmat->resultname ? matmat->resultname : "not specified"));
839:   PetscCall(PetscObjectComposeFunction((PetscObject)A, matmat->composedname, MatProductSetFromOptions_Shell_X));
840:   PetscFunctionReturn(PETSC_SUCCESS);
841: }

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

846:   Logically Collective; No Fortran Support

848:   Input Parameters:
849: + A        - the `MATSHELL` shell matrix
850: . ptype    - the product type
851: . symbolic - the function for the symbolic phase (can be `NULL`)
852: . numeric  - the function for the numerical phase
853: . destroy  - the function for the destruction of the needed data generated during the symbolic phase (can be `NULL`)
854: . Btype    - the matrix type for the matrix to be multiplied against
855: - Ctype    - the matrix type for the result (can be `NULL`)

857:   Level: advanced

859:   Example Usage:
860: .vb
861:   extern PetscErrorCode usersymbolic(Mat, Mat, Mat, void**);
862:   extern PetscErrorCode usernumeric(Mat, Mat, Mat, void*);
863:   extern PetscErrorCode userdestroy(void*);

865:   MatCreateShell(comm, m, n, M, N, ctx, &A);
866:   MatShellSetMatProductOperation(
867:     A, MATPRODUCT_AB, usersymbolic, usernumeric, userdestroy,MATSEQAIJ, MATDENSE
868:   );
869:   // create B of type SEQAIJ etc..
870:   MatProductCreate(A, B, PETSC_NULLPTR, &C);
871:   MatProductSetType(C, MATPRODUCT_AB);
872:   MatProductSetFromOptions(C);
873:   MatProductSymbolic(C); // actually runs the user defined symbolic operation
874:   MatProductNumeric(C); // actually runs the user defined numeric operation
875:   // use C = A * B
876: .ve

878:   Notes:
879:   `MATPRODUCT_ABC` is not supported yet.

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

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

888: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellSetContext()`, `MatSetOperation()`, `MatProductType`, `MatType`, `MatSetUp()`
889: @*/
890: 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)
891: {
892:   PetscFunctionBegin;
895:   PetscCheck(ptype != MATPRODUCT_ABC, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not for product type %s", MatProductTypes[ptype]);
896:   PetscCheck(numeric, PetscObjectComm((PetscObject)A), PETSC_ERR_USER, "Missing numeric routine, argument 4");
897:   PetscAssertPointer(Btype, 6);
898:   if (Ctype) PetscAssertPointer(Ctype, 7);
899:   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));
900:   PetscFunctionReturn(PETSC_SUCCESS);
901: }

903: 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)
904: {
905:   PetscBool   flg;
906:   char        composedname[256];
907:   MatRootName Bnames = MatRootNameList, Cnames = MatRootNameList;
908:   PetscMPIInt size;

910:   PetscFunctionBegin;
912:   while (Bnames) { /* user passed in the root name */
913:     PetscCall(PetscStrcmp(Btype, Bnames->rname, &flg));
914:     if (flg) break;
915:     Bnames = Bnames->next;
916:   }
917:   while (Cnames) { /* user passed in the root name */
918:     PetscCall(PetscStrcmp(Ctype, Cnames->rname, &flg));
919:     if (flg) break;
920:     Cnames = Cnames->next;
921:   }
922:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
923:   Btype = Bnames ? (size > 1 ? Bnames->mname : Bnames->sname) : Btype;
924:   Ctype = Cnames ? (size > 1 ? Cnames->mname : Cnames->sname) : Ctype;
925:   PetscCall(PetscSNPrintf(composedname, sizeof(composedname), "MatProductSetFromOptions_%s_%s_C", ((PetscObject)A)->type_name, Btype));
926:   PetscCall(MatShellSetMatProductOperation_Private(A, ptype, symbolic, numeric, destroy, composedname, Ctype));
927:   PetscFunctionReturn(PETSC_SUCCESS);
928: }

930: static PetscErrorCode MatCopy_Shell(Mat A, Mat B, MatStructure str)
931: {
932:   Mat_Shell              *shellA = (Mat_Shell *)A->data, *shellB = (Mat_Shell *)B->data;
933:   PetscBool               matflg;
934:   MatShellMatFunctionList matmatA;

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

940:   B->ops[0]      = A->ops[0];
941:   shellB->ops[0] = shellA->ops[0];

943:   if (shellA->ops->copy) PetscCall((*shellA->ops->copy)(A, B, str));
944:   shellB->vscale = shellA->vscale;
945:   shellB->vshift = shellA->vshift;
946:   if (shellA->dshift) {
947:     if (!shellB->dshift) PetscCall(VecDuplicate(shellA->dshift, &shellB->dshift));
948:     PetscCall(VecCopy(shellA->dshift, shellB->dshift));
949:   } else {
950:     PetscCall(VecDestroy(&shellB->dshift));
951:   }
952:   if (shellA->left) {
953:     if (!shellB->left) PetscCall(VecDuplicate(shellA->left, &shellB->left));
954:     PetscCall(VecCopy(shellA->left, shellB->left));
955:   } else {
956:     PetscCall(VecDestroy(&shellB->left));
957:   }
958:   if (shellA->right) {
959:     if (!shellB->right) PetscCall(VecDuplicate(shellA->right, &shellB->right));
960:     PetscCall(VecCopy(shellA->right, shellB->right));
961:   } else {
962:     PetscCall(VecDestroy(&shellB->right));
963:   }
964:   PetscCall(MatDestroy(&shellB->axpy));
965:   shellB->axpy_vscale = 0.0;
966:   shellB->axpy_state  = 0;
967:   if (shellA->axpy) {
968:     PetscCall(PetscObjectReference((PetscObject)shellA->axpy));
969:     shellB->axpy        = shellA->axpy;
970:     shellB->axpy_vscale = shellA->axpy_vscale;
971:     shellB->axpy_state  = shellA->axpy_state;
972:   }
973:   if (shellA->zrows) {
974:     PetscCall(ISDuplicate(shellA->zrows, &shellB->zrows));
975:     if (shellA->zcols) PetscCall(ISDuplicate(shellA->zcols, &shellB->zcols));
976:     PetscCall(VecDuplicate(shellA->zvals, &shellB->zvals));
977:     PetscCall(VecCopy(shellA->zvals, shellB->zvals));
978:     PetscCall(VecDuplicate(shellA->zvals_w, &shellB->zvals_w));
979:     PetscCall(PetscObjectReference((PetscObject)shellA->zvals_sct_r));
980:     PetscCall(PetscObjectReference((PetscObject)shellA->zvals_sct_c));
981:     shellB->zvals_sct_r = shellA->zvals_sct_r;
982:     shellB->zvals_sct_c = shellA->zvals_sct_c;
983:   }

985:   matmatA = shellA->matmat;
986:   if (matmatA) {
987:     while (matmatA->next) {
988:       PetscCall(MatShellSetMatProductOperation_Private(B, matmatA->ptype, matmatA->symbolic, matmatA->numeric, matmatA->destroy, matmatA->composedname, matmatA->resultname));
989:       matmatA = matmatA->next;
990:     }
991:   }
992:   PetscFunctionReturn(PETSC_SUCCESS);
993: }

995: static PetscErrorCode MatDuplicate_Shell(Mat mat, MatDuplicateOption op, Mat *M)
996: {
997:   PetscFunctionBegin;
998:   PetscCall(MatCreateShell(PetscObjectComm((PetscObject)mat), mat->rmap->n, mat->cmap->n, mat->rmap->N, mat->cmap->N, NULL, M));
999:   ((Mat_Shell *)(*M)->data)->ctxcontainer = ((Mat_Shell *)mat->data)->ctxcontainer;
1000:   PetscCall(PetscObjectCompose((PetscObject)(*M), "MatShell ctx", (PetscObject)((Mat_Shell *)(*M)->data)->ctxcontainer));
1001:   PetscCall(PetscObjectChangeTypeName((PetscObject)(*M), ((PetscObject)mat)->type_name));
1002:   if (op != MAT_DO_NOT_COPY_VALUES) PetscCall(MatCopy(mat, *M, SAME_NONZERO_PATTERN));
1003:   PetscFunctionReturn(PETSC_SUCCESS);
1004: }

1006: static PetscErrorCode MatMult_Shell(Mat A, Vec x, Vec y)
1007: {
1008:   Mat_Shell       *shell = (Mat_Shell *)A->data;
1009:   Vec              xx;
1010:   PetscObjectState instate, outstate;

1012:   PetscFunctionBegin;
1013:   PetscCall(MatShellPreZeroRight(A, x, &xx));
1014:   PetscCall(MatShellPreScaleRight(A, xx, &xx));
1015:   PetscCall(PetscObjectStateGet((PetscObject)y, &instate));
1016:   PetscCall((*shell->ops->mult)(A, xx, y));
1017:   PetscCall(PetscObjectStateGet((PetscObject)y, &outstate));
1018:   if (instate == outstate) {
1019:     /* increase the state of the output vector since the user did not update its state themself as should have been done */
1020:     PetscCall(PetscObjectStateIncrease((PetscObject)y));
1021:   }
1022:   PetscCall(MatShellShiftAndScale(A, xx, y));
1023:   PetscCall(MatShellPostScaleLeft(A, y));
1024:   PetscCall(MatShellPostZeroLeft(A, y));

1026:   if (shell->axpy) {
1027:     Mat              X;
1028:     PetscObjectState axpy_state;

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

1034:     PetscCall(MatCreateVecs(shell->axpy, shell->axpy_right ? NULL : &shell->axpy_right, shell->axpy_left ? NULL : &shell->axpy_left));
1035:     PetscCall(VecCopy(x, shell->axpy_right));
1036:     PetscCall(MatMult(shell->axpy, shell->axpy_right, shell->axpy_left));
1037:     PetscCall(VecAXPY(y, shell->axpy_vscale, shell->axpy_left));
1038:   }
1039:   PetscFunctionReturn(PETSC_SUCCESS);
1040: }

1042: static PetscErrorCode MatMultAdd_Shell(Mat A, Vec x, Vec y, Vec z)
1043: {
1044:   Mat_Shell *shell = (Mat_Shell *)A->data;

1046:   PetscFunctionBegin;
1047:   if (y == z) {
1048:     if (!shell->right_add_work) PetscCall(VecDuplicate(z, &shell->right_add_work));
1049:     PetscCall(MatMult(A, x, shell->right_add_work));
1050:     PetscCall(VecAXPY(z, 1.0, shell->right_add_work));
1051:   } else {
1052:     PetscCall(MatMult(A, x, z));
1053:     PetscCall(VecAXPY(z, 1.0, y));
1054:   }
1055:   PetscFunctionReturn(PETSC_SUCCESS);
1056: }

1058: static PetscErrorCode MatMultTranspose_Shell(Mat A, Vec x, Vec y)
1059: {
1060:   Mat_Shell       *shell = (Mat_Shell *)A->data;
1061:   Vec              xx;
1062:   PetscObjectState instate, outstate;

1064:   PetscFunctionBegin;
1065:   PetscCall(MatShellPreZeroLeft(A, x, &xx));
1066:   PetscCall(MatShellPreScaleLeft(A, xx, &xx));
1067:   PetscCall(PetscObjectStateGet((PetscObject)y, &instate));
1068:   PetscCall((*shell->ops->multtranspose)(A, xx, y));
1069:   PetscCall(PetscObjectStateGet((PetscObject)y, &outstate));
1070:   if (instate == outstate) {
1071:     /* increase the state of the output vector since the user did not update its state themself as should have been done */
1072:     PetscCall(PetscObjectStateIncrease((PetscObject)y));
1073:   }
1074:   PetscCall(MatShellShiftAndScale(A, xx, y));
1075:   PetscCall(MatShellPostScaleRight(A, y));
1076:   PetscCall(MatShellPostZeroRight(A, y));

1078:   if (shell->axpy) {
1079:     Mat              X;
1080:     PetscObjectState axpy_state;

1082:     PetscCall(MatShellGetContext(shell->axpy, &X));
1083:     PetscCall(PetscObjectStateGet((PetscObject)X, &axpy_state));
1084:     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,...)");
1085:     PetscCall(MatCreateVecs(shell->axpy, shell->axpy_right ? NULL : &shell->axpy_right, shell->axpy_left ? NULL : &shell->axpy_left));
1086:     PetscCall(VecCopy(x, shell->axpy_left));
1087:     PetscCall(MatMultTranspose(shell->axpy, shell->axpy_left, shell->axpy_right));
1088:     PetscCall(VecAXPY(y, shell->axpy_vscale, shell->axpy_right));
1089:   }
1090:   PetscFunctionReturn(PETSC_SUCCESS);
1091: }

1093: static PetscErrorCode MatMultTransposeAdd_Shell(Mat A, Vec x, Vec y, Vec z)
1094: {
1095:   Mat_Shell *shell = (Mat_Shell *)A->data;

1097:   PetscFunctionBegin;
1098:   if (y == z) {
1099:     if (!shell->left_add_work) PetscCall(VecDuplicate(z, &shell->left_add_work));
1100:     PetscCall(MatMultTranspose(A, x, shell->left_add_work));
1101:     PetscCall(VecAXPY(z, 1.0, shell->left_add_work));
1102:   } else {
1103:     PetscCall(MatMultTranspose(A, x, z));
1104:     PetscCall(VecAXPY(z, 1.0, y));
1105:   }
1106:   PetscFunctionReturn(PETSC_SUCCESS);
1107: }

1109: /*
1110:           diag(left)(vscale*A + diag(dshift) + vshift I)diag(right)
1111: */
1112: static PetscErrorCode MatGetDiagonal_Shell(Mat A, Vec v)
1113: {
1114:   Mat_Shell *shell = (Mat_Shell *)A->data;

1116:   PetscFunctionBegin;
1117:   if (shell->ops->getdiagonal) {
1118:     PetscCall((*shell->ops->getdiagonal)(A, v));
1119:   } else SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Must provide shell matrix with routine to return diagonal using\nMatShellSetOperation(S,MATOP_GET_DIAGONAL,...)");
1120:   PetscCall(VecScale(v, shell->vscale));
1121:   if (shell->dshift) PetscCall(VecAXPY(v, 1.0, shell->dshift));
1122:   PetscCall(VecShift(v, shell->vshift));
1123:   if (shell->left) PetscCall(VecPointwiseMult(v, v, shell->left));
1124:   if (shell->right) PetscCall(VecPointwiseMult(v, v, shell->right));
1125:   if (shell->zrows) {
1126:     PetscCall(VecScatterBegin(shell->zvals_sct_r, shell->zvals, v, INSERT_VALUES, SCATTER_REVERSE));
1127:     PetscCall(VecScatterEnd(shell->zvals_sct_r, shell->zvals, v, INSERT_VALUES, SCATTER_REVERSE));
1128:   }
1129:   if (shell->axpy) {
1130:     Mat              X;
1131:     PetscObjectState axpy_state;

1133:     PetscCall(MatShellGetContext(shell->axpy, &X));
1134:     PetscCall(PetscObjectStateGet((PetscObject)X, &axpy_state));
1135:     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,...)");
1136:     PetscCall(MatCreateVecs(shell->axpy, NULL, shell->axpy_left ? NULL : &shell->axpy_left));
1137:     PetscCall(MatGetDiagonal(shell->axpy, shell->axpy_left));
1138:     PetscCall(VecAXPY(v, shell->axpy_vscale, shell->axpy_left));
1139:   }
1140:   PetscFunctionReturn(PETSC_SUCCESS);
1141: }

1143: static PetscErrorCode MatShift_Shell(Mat Y, PetscScalar a)
1144: {
1145:   Mat_Shell *shell = (Mat_Shell *)Y->data;
1146:   PetscBool  flg;

1148:   PetscFunctionBegin;
1149:   PetscCall(MatHasCongruentLayouts(Y, &flg));
1150:   PetscCheck(flg, PetscObjectComm((PetscObject)Y), PETSC_ERR_SUP, "Cannot shift shell matrix if it is not congruent");
1151:   if (shell->left || shell->right) {
1152:     if (!shell->dshift) {
1153:       PetscCall(VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift));
1154:       PetscCall(VecSet(shell->dshift, a));
1155:     } else {
1156:       if (shell->left) PetscCall(VecPointwiseMult(shell->dshift, shell->dshift, shell->left));
1157:       if (shell->right) PetscCall(VecPointwiseMult(shell->dshift, shell->dshift, shell->right));
1158:       PetscCall(VecShift(shell->dshift, a));
1159:     }
1160:     if (shell->left) PetscCall(VecPointwiseDivide(shell->dshift, shell->dshift, shell->left));
1161:     if (shell->right) PetscCall(VecPointwiseDivide(shell->dshift, shell->dshift, shell->right));
1162:   } else shell->vshift += a;
1163:   if (shell->zrows) PetscCall(VecShift(shell->zvals, a));
1164:   PetscFunctionReturn(PETSC_SUCCESS);
1165: }

1167: static PetscErrorCode MatDiagonalSet_Shell_Private(Mat A, Vec D, PetscScalar s)
1168: {
1169:   Mat_Shell *shell = (Mat_Shell *)A->data;

1171:   PetscFunctionBegin;
1172:   if (!shell->dshift) PetscCall(VecDuplicate(D, &shell->dshift));
1173:   if (shell->left || shell->right) {
1174:     if (!shell->right_work) PetscCall(VecDuplicate(shell->left ? shell->left : shell->right, &shell->right_work));
1175:     if (shell->left && shell->right) {
1176:       PetscCall(VecPointwiseDivide(shell->right_work, D, shell->left));
1177:       PetscCall(VecPointwiseDivide(shell->right_work, shell->right_work, shell->right));
1178:     } else if (shell->left) {
1179:       PetscCall(VecPointwiseDivide(shell->right_work, D, shell->left));
1180:     } else {
1181:       PetscCall(VecPointwiseDivide(shell->right_work, D, shell->right));
1182:     }
1183:     PetscCall(VecAXPY(shell->dshift, s, shell->right_work));
1184:   } else {
1185:     PetscCall(VecAXPY(shell->dshift, s, D));
1186:   }
1187:   PetscFunctionReturn(PETSC_SUCCESS);
1188: }

1190: static PetscErrorCode MatDiagonalSet_Shell(Mat A, Vec D, InsertMode ins)
1191: {
1192:   Mat_Shell *shell = (Mat_Shell *)A->data;
1193:   Vec        d;
1194:   PetscBool  flg;

1196:   PetscFunctionBegin;
1197:   PetscCall(MatHasCongruentLayouts(A, &flg));
1198:   PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Cannot diagonal set or shift shell matrix if it is not congruent");
1199:   if (ins == INSERT_VALUES) {
1200:     PetscCall(VecDuplicate(D, &d));
1201:     PetscCall(MatGetDiagonal(A, d));
1202:     PetscCall(MatDiagonalSet_Shell_Private(A, d, -1.));
1203:     PetscCall(MatDiagonalSet_Shell_Private(A, D, 1.));
1204:     PetscCall(VecDestroy(&d));
1205:     if (shell->zrows) PetscCall(VecCopy(D, shell->zvals));
1206:   } else {
1207:     PetscCall(MatDiagonalSet_Shell_Private(A, D, 1.));
1208:     if (shell->zrows) PetscCall(VecAXPY(shell->zvals, 1.0, D));
1209:   }
1210:   PetscFunctionReturn(PETSC_SUCCESS);
1211: }

1213: static PetscErrorCode MatScale_Shell(Mat Y, PetscScalar a)
1214: {
1215:   Mat_Shell *shell = (Mat_Shell *)Y->data;

1217:   PetscFunctionBegin;
1218:   shell->vscale *= a;
1219:   shell->vshift *= a;
1220:   if (shell->dshift) PetscCall(VecScale(shell->dshift, a));
1221:   shell->axpy_vscale *= a;
1222:   if (shell->zrows) PetscCall(VecScale(shell->zvals, a));
1223:   PetscFunctionReturn(PETSC_SUCCESS);
1224: }

1226: static PetscErrorCode MatDiagonalScale_Shell(Mat Y, Vec left, Vec right)
1227: {
1228:   Mat_Shell *shell = (Mat_Shell *)Y->data;

1230:   PetscFunctionBegin;
1231:   if (left) {
1232:     if (!shell->left) {
1233:       PetscCall(VecDuplicate(left, &shell->left));
1234:       PetscCall(VecCopy(left, shell->left));
1235:     } else {
1236:       PetscCall(VecPointwiseMult(shell->left, shell->left, left));
1237:     }
1238:     if (shell->zrows) PetscCall(VecPointwiseMult(shell->zvals, shell->zvals, left));
1239:   }
1240:   if (right) {
1241:     if (!shell->right) {
1242:       PetscCall(VecDuplicate(right, &shell->right));
1243:       PetscCall(VecCopy(right, shell->right));
1244:     } else {
1245:       PetscCall(VecPointwiseMult(shell->right, shell->right, right));
1246:     }
1247:     if (shell->zrows) {
1248:       if (!shell->left_work) PetscCall(MatCreateVecs(Y, NULL, &shell->left_work));
1249:       PetscCall(VecSet(shell->zvals_w, 1.0));
1250:       PetscCall(VecScatterBegin(shell->zvals_sct_c, right, shell->zvals_w, INSERT_VALUES, SCATTER_FORWARD));
1251:       PetscCall(VecScatterEnd(shell->zvals_sct_c, right, shell->zvals_w, INSERT_VALUES, SCATTER_FORWARD));
1252:       PetscCall(VecPointwiseMult(shell->zvals, shell->zvals, shell->zvals_w));
1253:     }
1254:   }
1255:   if (shell->axpy) PetscCall(MatDiagonalScale(shell->axpy, left, right));
1256:   PetscFunctionReturn(PETSC_SUCCESS);
1257: }

1259: static PetscErrorCode MatAssemblyEnd_Shell(Mat Y, MatAssemblyType t)
1260: {
1261:   Mat_Shell *shell = (Mat_Shell *)Y->data;

1263:   PetscFunctionBegin;
1264:   if (t == MAT_FINAL_ASSEMBLY) {
1265:     shell->vshift      = 0.0;
1266:     shell->vscale      = 1.0;
1267:     shell->axpy_vscale = 0.0;
1268:     shell->axpy_state  = 0;
1269:     PetscCall(VecDestroy(&shell->dshift));
1270:     PetscCall(VecDestroy(&shell->left));
1271:     PetscCall(VecDestroy(&shell->right));
1272:     PetscCall(MatDestroy(&shell->axpy));
1273:     PetscCall(VecDestroy(&shell->axpy_left));
1274:     PetscCall(VecDestroy(&shell->axpy_right));
1275:     PetscCall(VecScatterDestroy(&shell->zvals_sct_c));
1276:     PetscCall(VecScatterDestroy(&shell->zvals_sct_r));
1277:     PetscCall(ISDestroy(&shell->zrows));
1278:     PetscCall(ISDestroy(&shell->zcols));
1279:   }
1280:   PetscFunctionReturn(PETSC_SUCCESS);
1281: }

1283: static PetscErrorCode MatMissingDiagonal_Shell(Mat A, PetscBool *missing, PetscInt *d)
1284: {
1285:   PetscFunctionBegin;
1286:   *missing = PETSC_FALSE;
1287:   PetscFunctionReturn(PETSC_SUCCESS);
1288: }

1290: static PetscErrorCode MatAXPY_Shell(Mat Y, PetscScalar a, Mat X, MatStructure str)
1291: {
1292:   Mat_Shell *shell = (Mat_Shell *)Y->data;

1294:   PetscFunctionBegin;
1295:   if (X == Y) {
1296:     PetscCall(MatScale(Y, 1.0 + a));
1297:     PetscFunctionReturn(PETSC_SUCCESS);
1298:   }
1299:   if (!shell->axpy) {
1300:     PetscCall(MatConvertFrom_Shell(X, MATSHELL, MAT_INITIAL_MATRIX, &shell->axpy));
1301:     shell->axpy_vscale = a;
1302:     PetscCall(PetscObjectStateGet((PetscObject)X, &shell->axpy_state));
1303:   } else {
1304:     PetscCall(MatAXPY(shell->axpy, a / shell->axpy_vscale, X, str));
1305:   }
1306:   PetscFunctionReturn(PETSC_SUCCESS);
1307: }

1309: static struct _MatOps MatOps_Values = {NULL,
1310:                                        NULL,
1311:                                        NULL,
1312:                                        NULL,
1313:                                        /* 4*/ MatMultAdd_Shell,
1314:                                        NULL,
1315:                                        MatMultTransposeAdd_Shell,
1316:                                        NULL,
1317:                                        NULL,
1318:                                        NULL,
1319:                                        /*10*/ NULL,
1320:                                        NULL,
1321:                                        NULL,
1322:                                        NULL,
1323:                                        NULL,
1324:                                        /*15*/ NULL,
1325:                                        NULL,
1326:                                        NULL,
1327:                                        MatDiagonalScale_Shell,
1328:                                        NULL,
1329:                                        /*20*/ NULL,
1330:                                        MatAssemblyEnd_Shell,
1331:                                        NULL,
1332:                                        NULL,
1333:                                        /*24*/ MatZeroRows_Shell,
1334:                                        NULL,
1335:                                        NULL,
1336:                                        NULL,
1337:                                        NULL,
1338:                                        /*29*/ NULL,
1339:                                        NULL,
1340:                                        NULL,
1341:                                        NULL,
1342:                                        NULL,
1343:                                        /*34*/ MatDuplicate_Shell,
1344:                                        NULL,
1345:                                        NULL,
1346:                                        NULL,
1347:                                        NULL,
1348:                                        /*39*/ MatAXPY_Shell,
1349:                                        NULL,
1350:                                        NULL,
1351:                                        NULL,
1352:                                        MatCopy_Shell,
1353:                                        /*44*/ NULL,
1354:                                        MatScale_Shell,
1355:                                        MatShift_Shell,
1356:                                        MatDiagonalSet_Shell,
1357:                                        MatZeroRowsColumns_Shell,
1358:                                        /*49*/ NULL,
1359:                                        NULL,
1360:                                        NULL,
1361:                                        NULL,
1362:                                        NULL,
1363:                                        /*54*/ NULL,
1364:                                        NULL,
1365:                                        NULL,
1366:                                        NULL,
1367:                                        NULL,
1368:                                        /*59*/ NULL,
1369:                                        MatDestroy_Shell,
1370:                                        NULL,
1371:                                        MatConvertFrom_Shell,
1372:                                        NULL,
1373:                                        /*64*/ NULL,
1374:                                        NULL,
1375:                                        NULL,
1376:                                        NULL,
1377:                                        NULL,
1378:                                        /*69*/ NULL,
1379:                                        NULL,
1380:                                        MatConvert_Shell,
1381:                                        NULL,
1382:                                        NULL,
1383:                                        /*74*/ NULL,
1384:                                        NULL,
1385:                                        NULL,
1386:                                        NULL,
1387:                                        NULL,
1388:                                        /*79*/ NULL,
1389:                                        NULL,
1390:                                        NULL,
1391:                                        NULL,
1392:                                        NULL,
1393:                                        /*84*/ NULL,
1394:                                        NULL,
1395:                                        NULL,
1396:                                        NULL,
1397:                                        NULL,
1398:                                        /*89*/ NULL,
1399:                                        NULL,
1400:                                        NULL,
1401:                                        NULL,
1402:                                        NULL,
1403:                                        /*94*/ NULL,
1404:                                        NULL,
1405:                                        NULL,
1406:                                        NULL,
1407:                                        NULL,
1408:                                        /*99*/ NULL,
1409:                                        NULL,
1410:                                        NULL,
1411:                                        NULL,
1412:                                        NULL,
1413:                                        /*104*/ NULL,
1414:                                        NULL,
1415:                                        NULL,
1416:                                        NULL,
1417:                                        NULL,
1418:                                        /*109*/ NULL,
1419:                                        NULL,
1420:                                        NULL,
1421:                                        NULL,
1422:                                        MatMissingDiagonal_Shell,
1423:                                        /*114*/ NULL,
1424:                                        NULL,
1425:                                        NULL,
1426:                                        NULL,
1427:                                        NULL,
1428:                                        /*119*/ NULL,
1429:                                        NULL,
1430:                                        NULL,
1431:                                        NULL,
1432:                                        NULL,
1433:                                        /*124*/ NULL,
1434:                                        NULL,
1435:                                        NULL,
1436:                                        NULL,
1437:                                        NULL,
1438:                                        /*129*/ NULL,
1439:                                        NULL,
1440:                                        NULL,
1441:                                        NULL,
1442:                                        NULL,
1443:                                        /*134*/ NULL,
1444:                                        NULL,
1445:                                        NULL,
1446:                                        NULL,
1447:                                        NULL,
1448:                                        /*139*/ NULL,
1449:                                        NULL,
1450:                                        NULL,
1451:                                        NULL,
1452:                                        NULL,
1453:                                        /*144*/ NULL,
1454:                                        NULL,
1455:                                        NULL,
1456:                                        NULL,
1457:                                        NULL,
1458:                                        NULL,
1459:                                        /*150*/ NULL,
1460:                                        NULL};

1462: static PetscErrorCode MatShellSetContext_Shell(Mat mat, void *ctx)
1463: {
1464:   Mat_Shell *shell = (Mat_Shell *)mat->data;

1466:   PetscFunctionBegin;
1467:   if (ctx) {
1468:     PetscContainer ctxcontainer;
1469:     PetscCall(PetscContainerCreate(PetscObjectComm((PetscObject)mat), &ctxcontainer));
1470:     PetscCall(PetscContainerSetPointer(ctxcontainer, ctx));
1471:     PetscCall(PetscObjectCompose((PetscObject)mat, "MatShell ctx", (PetscObject)ctxcontainer));
1472:     shell->ctxcontainer = ctxcontainer;
1473:     PetscCall(PetscContainerDestroy(&ctxcontainer));
1474:   } else {
1475:     PetscCall(PetscObjectCompose((PetscObject)mat, "MatShell ctx", NULL));
1476:     shell->ctxcontainer = NULL;
1477:   }

1479:   PetscFunctionReturn(PETSC_SUCCESS);
1480: }

1482: static PetscErrorCode MatShellSetContextDestroy_Shell(Mat mat, PetscErrorCode (*f)(void *))
1483: {
1484:   Mat_Shell *shell = (Mat_Shell *)mat->data;

1486:   PetscFunctionBegin;
1487:   if (shell->ctxcontainer) PetscCall(PetscContainerSetUserDestroy(shell->ctxcontainer, f));
1488:   PetscFunctionReturn(PETSC_SUCCESS);
1489: }

1491: static PetscErrorCode MatShellSetVecType_Shell(Mat mat, VecType vtype)
1492: {
1493:   PetscFunctionBegin;
1494:   PetscCall(PetscFree(mat->defaultvectype));
1495:   PetscCall(PetscStrallocpy(vtype, (char **)&mat->defaultvectype));
1496:   PetscFunctionReturn(PETSC_SUCCESS);
1497: }

1499: static PetscErrorCode MatShellSetManageScalingShifts_Shell(Mat A)
1500: {
1501:   Mat_Shell *shell = (Mat_Shell *)A->data;

1503:   PetscFunctionBegin;
1504:   shell->managescalingshifts = PETSC_FALSE;
1505:   A->ops->diagonalset        = NULL;
1506:   A->ops->diagonalscale      = NULL;
1507:   A->ops->scale              = NULL;
1508:   A->ops->shift              = NULL;
1509:   A->ops->axpy               = NULL;
1510:   PetscFunctionReturn(PETSC_SUCCESS);
1511: }

1513: static PetscErrorCode MatShellSetOperation_Shell(Mat mat, MatOperation op, void (*f)(void))
1514: {
1515:   Mat_Shell *shell = (Mat_Shell *)mat->data;

1517:   PetscFunctionBegin;
1518:   switch (op) {
1519:   case MATOP_DESTROY:
1520:     shell->ops->destroy = (PetscErrorCode(*)(Mat))f;
1521:     break;
1522:   case MATOP_VIEW:
1523:     if (!mat->ops->viewnative) mat->ops->viewnative = mat->ops->view;
1524:     mat->ops->view = (PetscErrorCode(*)(Mat, PetscViewer))f;
1525:     break;
1526:   case MATOP_COPY:
1527:     shell->ops->copy = (PetscErrorCode(*)(Mat, Mat, MatStructure))f;
1528:     break;
1529:   case MATOP_DIAGONAL_SET:
1530:   case MATOP_DIAGONAL_SCALE:
1531:   case MATOP_SHIFT:
1532:   case MATOP_SCALE:
1533:   case MATOP_AXPY:
1534:   case MATOP_ZERO_ROWS:
1535:   case MATOP_ZERO_ROWS_COLUMNS:
1536:     PetscCheck(!shell->managescalingshifts, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "MATSHELL is managing scalings and shifts, see MatShellSetManageScalingShifts()");
1537:     (((void (**)(void))mat->ops)[op]) = f;
1538:     break;
1539:   case MATOP_GET_DIAGONAL:
1540:     if (shell->managescalingshifts) {
1541:       shell->ops->getdiagonal = (PetscErrorCode(*)(Mat, Vec))f;
1542:       mat->ops->getdiagonal   = MatGetDiagonal_Shell;
1543:     } else {
1544:       shell->ops->getdiagonal = NULL;
1545:       mat->ops->getdiagonal   = (PetscErrorCode(*)(Mat, Vec))f;
1546:     }
1547:     break;
1548:   case MATOP_MULT:
1549:     if (shell->managescalingshifts) {
1550:       shell->ops->mult = (PetscErrorCode(*)(Mat, Vec, Vec))f;
1551:       mat->ops->mult   = MatMult_Shell;
1552:     } else {
1553:       shell->ops->mult = NULL;
1554:       mat->ops->mult   = (PetscErrorCode(*)(Mat, Vec, Vec))f;
1555:     }
1556:     break;
1557:   case MATOP_MULT_TRANSPOSE:
1558:     if (shell->managescalingshifts) {
1559:       shell->ops->multtranspose = (PetscErrorCode(*)(Mat, Vec, Vec))f;
1560:       mat->ops->multtranspose   = MatMultTranspose_Shell;
1561:     } else {
1562:       shell->ops->multtranspose = NULL;
1563:       mat->ops->multtranspose   = (PetscErrorCode(*)(Mat, Vec, Vec))f;
1564:     }
1565:     break;
1566:   default:
1567:     (((void (**)(void))mat->ops)[op]) = f;
1568:     break;
1569:   }
1570:   PetscFunctionReturn(PETSC_SUCCESS);
1571: }

1573: static PetscErrorCode MatShellGetOperation_Shell(Mat mat, MatOperation op, void (**f)(void))
1574: {
1575:   Mat_Shell *shell = (Mat_Shell *)mat->data;

1577:   PetscFunctionBegin;
1578:   switch (op) {
1579:   case MATOP_DESTROY:
1580:     *f = (void (*)(void))shell->ops->destroy;
1581:     break;
1582:   case MATOP_VIEW:
1583:     *f = (void (*)(void))mat->ops->view;
1584:     break;
1585:   case MATOP_COPY:
1586:     *f = (void (*)(void))shell->ops->copy;
1587:     break;
1588:   case MATOP_DIAGONAL_SET:
1589:   case MATOP_DIAGONAL_SCALE:
1590:   case MATOP_SHIFT:
1591:   case MATOP_SCALE:
1592:   case MATOP_AXPY:
1593:   case MATOP_ZERO_ROWS:
1594:   case MATOP_ZERO_ROWS_COLUMNS:
1595:     *f = (((void (**)(void))mat->ops)[op]);
1596:     break;
1597:   case MATOP_GET_DIAGONAL:
1598:     if (shell->ops->getdiagonal) *f = (void (*)(void))shell->ops->getdiagonal;
1599:     else *f = (((void (**)(void))mat->ops)[op]);
1600:     break;
1601:   case MATOP_MULT:
1602:     if (shell->ops->mult) *f = (void (*)(void))shell->ops->mult;
1603:     else *f = (((void (**)(void))mat->ops)[op]);
1604:     break;
1605:   case MATOP_MULT_TRANSPOSE:
1606:     if (shell->ops->multtranspose) *f = (void (*)(void))shell->ops->multtranspose;
1607:     else *f = (((void (**)(void))mat->ops)[op]);
1608:     break;
1609:   default:
1610:     *f = (((void (**)(void))mat->ops)[op]);
1611:   }
1612:   PetscFunctionReturn(PETSC_SUCCESS);
1613: }

1615: /*MC
1616:    MATSHELL - MATSHELL = "shell" - A matrix type to be used to define your own matrix type -- perhaps matrix-free.

1618:   Level: advanced

1620: .seealso: [](ch_matrices), `Mat`, `MatCreateShell()`
1621: M*/

1623: PETSC_EXTERN PetscErrorCode MatCreate_Shell(Mat A)
1624: {
1625:   Mat_Shell *b;

1627:   PetscFunctionBegin;
1628:   PetscCall(PetscNew(&b));
1629:   A->data   = (void *)b;
1630:   A->ops[0] = MatOps_Values;

1632:   b->ctxcontainer        = NULL;
1633:   b->vshift              = 0.0;
1634:   b->vscale              = 1.0;
1635:   b->managescalingshifts = PETSC_TRUE;
1636:   A->assembled           = PETSC_TRUE;
1637:   A->preallocated        = PETSC_FALSE;

1639:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellGetContext_C", MatShellGetContext_Shell));
1640:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetContext_C", MatShellSetContext_Shell));
1641:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetContextDestroy_C", MatShellSetContextDestroy_Shell));
1642:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetVecType_C", MatShellSetVecType_Shell));
1643:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetManageScalingShifts_C", MatShellSetManageScalingShifts_Shell));
1644:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetOperation_C", MatShellSetOperation_Shell));
1645:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellGetOperation_C", MatShellGetOperation_Shell));
1646:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetMatProductOperation_C", MatShellSetMatProductOperation_Shell));
1647:   PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATSHELL));
1648:   PetscFunctionReturn(PETSC_SUCCESS);
1649: }

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

1655:   Collective

1657:   Input Parameters:
1658: + comm - MPI communicator
1659: . m    - number of local rows (must be given)
1660: . n    - number of local columns (must be given)
1661: . M    - number of global rows (may be `PETSC_DETERMINE`)
1662: . N    - number of global columns (may be `PETSC_DETERMINE`)
1663: - ctx  - pointer to data needed by the shell matrix routines

1665:   Output Parameter:
1666: . A - the matrix

1668:   Level: advanced

1670:   Example Usage:
1671: .vb
1672:   extern PetscErrorCode mult(Mat, Vec, Vec);

1674:   MatCreateShell(comm, m, n, M, N, ctx, &mat);
1675:   MatShellSetOperation(mat, MATOP_MULT, (void(*)(void))mult);
1676:   // Use matrix for operations that have been set
1677:   MatDestroy(mat);
1678: .ve

1680:   Notes:
1681:   The shell matrix type is intended to provide a simple class to use
1682:   with `KSP` (such as, for use with matrix-free methods). You should not
1683:   use the shell type if you plan to define a complete matrix class.

1685:   PETSc requires that matrices and vectors being used for certain
1686:   operations are partitioned accordingly.  For example, when
1687:   creating a shell matrix, `A`, that supports parallel matrix-vector
1688:   products using `MatMult`(A,x,y) the user should set the number
1689:   of local matrix rows to be the number of local elements of the
1690:   corresponding result vector, y. Note that this is information is
1691:   required for use of the matrix interface routines, even though
1692:   the shell matrix may not actually be physically partitioned.
1693:   For example,

1695: .vb
1696:      Vec x, y
1697:      extern PetscErrorCode mult(Mat,Vec,Vec);
1698:      Mat A

1700:      VecCreateMPI(comm,PETSC_DECIDE,M,&y);
1701:      VecCreateMPI(comm,PETSC_DECIDE,N,&x);
1702:      VecGetLocalSize(y,&m);
1703:      VecGetLocalSize(x,&n);
1704:      MatCreateShell(comm,m,n,M,N,ctx,&A);
1705:      MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
1706:      MatMult(A,x,y);
1707:      MatDestroy(&A);
1708:      VecDestroy(&y);
1709:      VecDestroy(&x);
1710: .ve

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

1715:   Developer Notes:
1716:   For rectangular matrices do all the scalings and shifts make sense?

1718:   Regarding shifting and scaling. The general form is

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

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

1726:   A is the user provided function.

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

1733:   Matrix product operations (i.e. `MatMat()`, `MatTransposeMat()` etc) can be specified using `MatShellSetMatProductOperation()`

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

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

1743: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatShellSetOperation()`, `MatHasOperation()`, `MatShellGetContext()`, `MatShellSetContext()`, `MatShellSetManageScalingShifts()`, `MatShellSetMatProductOperation()`
1744: @*/
1745: PetscErrorCode MatCreateShell(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, PetscInt N, void *ctx, Mat *A)
1746: {
1747:   PetscFunctionBegin;
1748:   PetscCall(MatCreate(comm, A));
1749:   PetscCall(MatSetSizes(*A, m, n, M, N));
1750:   PetscCall(MatSetType(*A, MATSHELL));
1751:   PetscCall(MatShellSetContext(*A, ctx));
1752:   PetscCall(MatSetUp(*A));
1753:   PetscFunctionReturn(PETSC_SUCCESS);
1754: }

1756: /*@
1757:   MatShellSetContext - sets the context for a `MATSHELL` shell matrix

1759:   Logically Collective

1761:   Input Parameters:
1762: + mat - the `MATSHELL` shell matrix
1763: - ctx - the context

1765:   Level: advanced

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

1771: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`
1772: @*/
1773: PetscErrorCode MatShellSetContext(Mat mat, void *ctx)
1774: {
1775:   PetscFunctionBegin;
1777:   PetscTryMethod(mat, "MatShellSetContext_C", (Mat, void *), (mat, ctx));
1778:   PetscFunctionReturn(PETSC_SUCCESS);
1779: }

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

1784:   Logically Collective

1786:   Input Parameters:
1787: + mat - the shell matrix
1788: - f   - the context destroy function

1790:   Level: advanced

1792:   Note:
1793:   If the `MatShell` is never duplicated, the behavior of this function is equivalent
1794:   to `MatShellSetOperation`(`Mat`,`MATOP_DESTROY`,f). However, `MatShellSetContextDestroy()`
1795:   ensures proper reference counting for the user provided context data in the case that
1796:   the `MATSHELL` is duplicated.

1798: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellSetContext()`
1799: @*/
1800: PetscErrorCode MatShellSetContextDestroy(Mat mat, PetscErrorCode (*f)(void *))
1801: {
1802:   PetscFunctionBegin;
1804:   PetscTryMethod(mat, "MatShellSetContextDestroy_C", (Mat, PetscErrorCode(*)(void *)), (mat, f));
1805:   PetscFunctionReturn(PETSC_SUCCESS);
1806: }

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

1811:   Logically Collective

1813:   Input Parameters:
1814: + mat   - the `MATSHELL` shell matrix
1815: - vtype - type to use for creating vectors

1817:   Level: advanced

1819: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateVecs()`
1820: @*/
1821: PetscErrorCode MatShellSetVecType(Mat mat, VecType vtype)
1822: {
1823:   PetscFunctionBegin;
1824:   PetscTryMethod(mat, "MatShellSetVecType_C", (Mat, VecType), (mat, vtype));
1825:   PetscFunctionReturn(PETSC_SUCCESS);
1826: }

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

1832:   Logically Collective

1834:   Input Parameter:
1835: . A - the `MATSHELL` shell matrix

1837:   Level: advanced

1839: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellSetContext()`, `MatShellSetOperation()`
1840: @*/
1841: PetscErrorCode MatShellSetManageScalingShifts(Mat A)
1842: {
1843:   PetscFunctionBegin;
1845:   PetscTryMethod(A, "MatShellSetManageScalingShifts_C", (Mat), (A));
1846:   PetscFunctionReturn(PETSC_SUCCESS);
1847: }

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

1852:   Logically Collective; No Fortran Support

1854:   Input Parameters:
1855: + mat  - the `MATSHELL` shell matrix
1856: . f    - the function
1857: . base - differences are computed around this vector, see `MatMFFDSetBase()`, for Jacobians this is the point at which the Jacobian is being evaluated
1858: - ctx  - an optional context for the function

1860:   Output Parameter:
1861: . flg - `PETSC_TRUE` if the multiply is likely correct

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

1866:   Level: advanced

1868: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMultTranspose()`
1869: @*/
1870: PetscErrorCode MatShellTestMult(Mat mat, PetscErrorCode (*f)(void *, Vec, Vec), Vec base, void *ctx, PetscBool *flg)
1871: {
1872:   PetscInt  m, n;
1873:   Mat       mf, Dmf, Dmat, Ddiff;
1874:   PetscReal Diffnorm, Dmfnorm;
1875:   PetscBool v = PETSC_FALSE, flag = PETSC_TRUE;

1877:   PetscFunctionBegin;
1879:   PetscCall(PetscOptionsHasName(NULL, ((PetscObject)mat)->prefix, "-mat_shell_test_mult_view", &v));
1880:   PetscCall(MatGetLocalSize(mat, &m, &n));
1881:   PetscCall(MatCreateMFFD(PetscObjectComm((PetscObject)mat), m, n, PETSC_DECIDE, PETSC_DECIDE, &mf));
1882:   PetscCall(MatMFFDSetFunction(mf, f, ctx));
1883:   PetscCall(MatMFFDSetBase(mf, base, NULL));

1885:   PetscCall(MatComputeOperator(mf, MATAIJ, &Dmf));
1886:   PetscCall(MatComputeOperator(mat, MATAIJ, &Dmat));

1888:   PetscCall(MatDuplicate(Dmat, MAT_COPY_VALUES, &Ddiff));
1889:   PetscCall(MatAXPY(Ddiff, -1.0, Dmf, DIFFERENT_NONZERO_PATTERN));
1890:   PetscCall(MatNorm(Ddiff, NORM_FROBENIUS, &Diffnorm));
1891:   PetscCall(MatNorm(Dmf, NORM_FROBENIUS, &Dmfnorm));
1892:   if (Diffnorm / Dmfnorm > 10 * PETSC_SQRT_MACHINE_EPSILON) {
1893:     flag = PETSC_FALSE;
1894:     if (v) {
1895:       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)));
1896:       PetscCall(MatViewFromOptions(Ddiff, (PetscObject)mat, "-mat_shell_test_mult_view"));
1897:       PetscCall(MatViewFromOptions(Dmf, (PetscObject)mat, "-mat_shell_test_mult_view"));
1898:       PetscCall(MatViewFromOptions(Dmat, (PetscObject)mat, "-mat_shell_test_mult_view"));
1899:     }
1900:   } else if (v) {
1901:     PetscCall(PetscPrintf(PetscObjectComm((PetscObject)mat), "MATSHELL and matrix-free multiple appear to produce the same results\n"));
1902:   }
1903:   if (flg) *flg = flag;
1904:   PetscCall(MatDestroy(&Ddiff));
1905:   PetscCall(MatDestroy(&mf));
1906:   PetscCall(MatDestroy(&Dmf));
1907:   PetscCall(MatDestroy(&Dmat));
1908:   PetscFunctionReturn(PETSC_SUCCESS);
1909: }

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

1914:   Logically Collective; No Fortran Support

1916:   Input Parameters:
1917: + mat  - the `MATSHELL` shell matrix
1918: . f    - the function
1919: . base - differences are computed around this vector, see `MatMFFDSetBase()`, for Jacobians this is the point at which the Jacobian is being evaluated
1920: - ctx  - an optional context for the function

1922:   Output Parameter:
1923: . flg - `PETSC_TRUE` if the multiply is likely correct

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

1928:   Level: advanced

1930: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMult()`
1931: @*/
1932: PetscErrorCode MatShellTestMultTranspose(Mat mat, PetscErrorCode (*f)(void *, Vec, Vec), Vec base, void *ctx, PetscBool *flg)
1933: {
1934:   Vec       x, y, z;
1935:   PetscInt  m, n, M, N;
1936:   Mat       mf, Dmf, Dmat, Ddiff;
1937:   PetscReal Diffnorm, Dmfnorm;
1938:   PetscBool v = PETSC_FALSE, flag = PETSC_TRUE;

1940:   PetscFunctionBegin;
1942:   PetscCall(PetscOptionsHasName(NULL, ((PetscObject)mat)->prefix, "-mat_shell_test_mult_transpose_view", &v));
1943:   PetscCall(MatCreateVecs(mat, &x, &y));
1944:   PetscCall(VecDuplicate(y, &z));
1945:   PetscCall(MatGetLocalSize(mat, &m, &n));
1946:   PetscCall(MatGetSize(mat, &M, &N));
1947:   PetscCall(MatCreateMFFD(PetscObjectComm((PetscObject)mat), m, n, M, N, &mf));
1948:   PetscCall(MatMFFDSetFunction(mf, f, ctx));
1949:   PetscCall(MatMFFDSetBase(mf, base, NULL));
1950:   PetscCall(MatComputeOperator(mf, MATAIJ, &Dmf));
1951:   PetscCall(MatTranspose(Dmf, MAT_INPLACE_MATRIX, &Dmf));
1952:   PetscCall(MatComputeOperatorTranspose(mat, MATAIJ, &Dmat));

1954:   PetscCall(MatDuplicate(Dmat, MAT_COPY_VALUES, &Ddiff));
1955:   PetscCall(MatAXPY(Ddiff, -1.0, Dmf, DIFFERENT_NONZERO_PATTERN));
1956:   PetscCall(MatNorm(Ddiff, NORM_FROBENIUS, &Diffnorm));
1957:   PetscCall(MatNorm(Dmf, NORM_FROBENIUS, &Dmfnorm));
1958:   if (Diffnorm / Dmfnorm > 10 * PETSC_SQRT_MACHINE_EPSILON) {
1959:     flag = PETSC_FALSE;
1960:     if (v) {
1961:       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)));
1962:       PetscCall(MatViewFromOptions(Ddiff, (PetscObject)mat, "-mat_shell_test_mult_transpose_view"));
1963:       PetscCall(MatViewFromOptions(Dmf, (PetscObject)mat, "-mat_shell_test_mult_transpose_view"));
1964:       PetscCall(MatViewFromOptions(Dmat, (PetscObject)mat, "-mat_shell_test_mult_transpose_view"));
1965:     }
1966:   } else if (v) {
1967:     PetscCall(PetscPrintf(PetscObjectComm((PetscObject)mat), "MATSHELL transpose and matrix-free multiple appear to produce the same results\n"));
1968:   }
1969:   if (flg) *flg = flag;
1970:   PetscCall(MatDestroy(&mf));
1971:   PetscCall(MatDestroy(&Dmat));
1972:   PetscCall(MatDestroy(&Ddiff));
1973:   PetscCall(MatDestroy(&Dmf));
1974:   PetscCall(VecDestroy(&x));
1975:   PetscCall(VecDestroy(&y));
1976:   PetscCall(VecDestroy(&z));
1977:   PetscFunctionReturn(PETSC_SUCCESS);
1978: }

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

1983:   Logically Collective

1985:   Input Parameters:
1986: + mat - the `MATSHELL` shell matrix
1987: . op  - the name of the operation
1988: - g   - the function that provides the operation.

1990:   Level: advanced

1992:   Example Usage:
1993: .vb
1994:   extern PetscErrorCode usermult(Mat, Vec, Vec);

1996:   MatCreateShell(comm, m, n, M, N, ctx, &A);
1997:   MatShellSetOperation(A, MATOP_MULT, (void(*)(void))usermult);
1998: .ve

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

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

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

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

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

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

2026: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellSetContext()`, `MatSetOperation()`, `MatShellSetManageScalingShifts()`, `MatShellSetMatProductOperation()`
2027: @*/
2028: PetscErrorCode MatShellSetOperation(Mat mat, MatOperation op, void (*g)(void))
2029: {
2030:   PetscFunctionBegin;
2032:   PetscTryMethod(mat, "MatShellSetOperation_C", (Mat, MatOperation, void (*)(void)), (mat, op, g));
2033:   PetscFunctionReturn(PETSC_SUCCESS);
2034: }

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

2039:   Not Collective

2041:   Input Parameters:
2042: + mat - the `MATSHELL` shell matrix
2043: - op  - the name of the operation

2045:   Output Parameter:
2046: . g - the function that provides the operation.

2048:   Level: advanced

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

2056:   All user-provided functions have the same calling
2057:   sequence as the usual matrix interface routines, since they
2058:   are intended to be accessed via the usual matrix interface
2059:   routines, e.g.,
2060: $       MatMult(Mat, Vec, Vec) -> usermult(Mat, Vec, Vec)

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

2066: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellSetOperation()`, `MatShellSetContext()`
2067: @*/
2068: PetscErrorCode MatShellGetOperation(Mat mat, MatOperation op, void (**g)(void))
2069: {
2070:   PetscFunctionBegin;
2072:   PetscUseMethod(mat, "MatShellGetOperation_C", (Mat, MatOperation, void (**)(void)), (mat, op, g));
2073:   PetscFunctionReturn(PETSC_SUCCESS);
2074: }

2076: /*@
2077:   MatIsShell - Inquires if a matrix is derived from `MATSHELL`

2079:   Input Parameter:
2080: . mat - the matrix

2082:   Output Parameter:
2083: . flg - the boolean value

2085:   Level: developer

2087:   Developer Notes:
2088:   In the future, we should allow the object type name to be changed still using the `MATSHELL` data structure for other matrices
2089:   (i.e. `MATTRANSPOSEVIRTUAL`, `MATSCHURCOMPLEMENT` etc)

2091: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MATMFFD`, `MatCreateShell()`, `MATTRANSPOSEVIRTUAL`, `MATSCHURCOMPLEMENT`
2092: @*/
2093: PetscErrorCode MatIsShell(Mat mat, PetscBool *flg)
2094: {
2095:   PetscFunctionBegin;
2097:   PetscAssertPointer(flg, 2);
2098:   *flg = (PetscBool)(mat->ops->destroy == MatDestroy_Shell);
2099:   PetscFunctionReturn(PETSC_SUCCESS);
2100: }