Actual source code: shell.c


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

  8: #include <petsc/private/matimpl.h>
  9: #include <petsc/private/vecimpl.h>

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

247:     Not Collective

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

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

255:     Level: advanced

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

847:    Logically Collective; No Fortran Support

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

858:    Level: advanced

860:     Usage:
861: .vb
862:       extern PetscErrorCode usersymbolic(Mat,Mat,Mat,void**);
863:       extern PetscErrorCode usernumeric(Mat,Mat,Mat,void*);
864:       extern PetscErrorCode userdestroy(void*);
865:       MatCreateShell(comm,m,n,M,N,ctx,&A);
866:       MatShellSetMatProductOperation(A,MATPRODUCT_AB,usersymbolic,usernumeric,userdestroy,MATSEQAIJ,MATDENSE);
867:       [ create B of type SEQAIJ etc..]
868:       MatProductCreate(A,B,NULL,&C);
869:       MatProductSetType(C,MATPRODUCT_AB);
870:       MatProductSetFromOptions(C);
871:       MatProductSymbolic(C); -> actually runs the user defined symbolic operation
872:       MatProductNumeric(C); -> actually runs the user defined numeric operation
873:       [ use C = A*B ]
874: .ve

876:     Notes:
877:     `MATPRODUCT_ABC` is not supported yet.

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

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

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

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

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

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

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

938:   PetscCall(PetscMemcpy(B->ops, A->ops, sizeof(struct _MatOps)));
939:   PetscCall(PetscMemcpy(shellB->ops, shellA->ops, sizeof(struct _MatShellOps)));

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

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

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

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

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

1024:   if (shell->axpy) {
1025:     Mat              X;
1026:     PetscObjectState axpy_state;

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

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

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

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

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

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

1076:   if (shell->axpy) {
1077:     Mat              X;
1078:     PetscObjectState axpy_state;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1477:   PetscFunctionReturn(PETSC_SUCCESS);
1478: }

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

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

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

1497: PetscErrorCode MatShellSetManageScalingShifts_Shell(Mat A)
1498: {
1499:   Mat_Shell *shell = (Mat_Shell *)A->data;

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

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

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

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

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

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

1616:   Level: advanced

1618: .seealso: [](ch_matrices), `Mat`, `MatCreateShell()`
1619: M*/

1621: PETSC_EXTERN PetscErrorCode MatCreate_Shell(Mat A)
1622: {
1623:   Mat_Shell *b;

1625:   PetscFunctionBegin;
1626:   PetscCall(PetscMemcpy(A->ops, &MatOps_Values, sizeof(struct _MatOps)));

1628:   PetscCall(PetscNew(&b));
1629:   A->data = (void *)b;

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

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

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

1654:   Collective

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

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

1667:    Level: advanced

1669:   Usage:
1670: .vb
1671:     extern PetscErrorCode mult(Mat,Vec,Vec);
1672:     MatCreateShell(comm,m,n,M,N,ctx,&mat);
1673:     MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
1674:     [ Use matrix for operations that have been set ]
1675:     MatDestroy(mat);
1676: .ve

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

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

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

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

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

1713:     Developer Notes:
1714:     For rectangular matrices do all the scalings and shifts make sense?

1716:     Regarding shifting and scaling. The general form is

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

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

1724:           A is the user provided function.

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

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

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

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

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

1754: /*@
1755:     MatShellSetContext - sets the context for a `MATSHELL` shell matrix

1757:    Logically Collective

1759:     Input Parameters:
1760: +   mat - the `MATSHELL` shell matrix
1761: -   ctx - the context

1763:    Level: advanced

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

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

1779: /*@
1780:     MatShellSetContextDestroy - sets the destroy function for a `MATSHELL` shell matrix context

1782:    Logically Collective

1784:     Input Parameters:
1785: +   mat - the shell matrix
1786: -   f - the context destroy function

1788:    Level: advanced

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

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

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

1809:  Logically Collective

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

1815:  Level: advanced

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

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

1830:    Logically Collective

1832:     Input Parameter:
1833: .   mat - the `MATSHELL` shell matrix

1835:   Level: advanced

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

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

1850:    Logically Collective; No Fortran Support

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

1858:    Output Parameter:
1859: .   flg - `PETSC_TRUE` if the multiply is likely correct

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

1864:    Level: advanced

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

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

1883:   PetscCall(MatComputeOperator(mf, MATAIJ, &Dmf));
1884:   PetscCall(MatComputeOperator(mat, MATAIJ, &Dmat));

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

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

1912:    Logically Collective; No Fortran Support

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

1920:    Output Parameter:
1921: .   flg - `PETSC_TRUE` if the multiply is likely correct

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

1926:    Level: advanced

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

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

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

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

1981:    Logically Collective

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

1988:    Level: advanced

1990:     Usage:
1991: .vb
1992:       extern PetscErrorCode usermult(Mat,Vec,Vec);
1993:       MatCreateShell(comm,m,n,M,N,ctx,&A);
1994:       MatShellSetOperation(A,MATOP_MULT,(void(*)(void))usermult);
1995: .ve

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

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

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

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

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

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

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

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

2036:     Not Collective

2038:     Input Parameters:
2039: +   mat - the `MATSHELL` shell matrix
2040: -   op - the name of the operation

2042:     Output Parameter:
2043: .   g - the function that provides the operation.

2045:     Level: advanced

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

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

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

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

2073: /*@
2074:     MatIsShell - Inquires if a matrix is derived from `MATSHELL`

2076:     Input Parameter:
2077: .   mat - the matrix

2079:     Output Parameter:
2080: .   flg - the boolean value

2082:     Level: developer

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

2088: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MATMFFD`, `MatCreateShell()`, `MATTRANSPOSEVIRTUAL`, `MATSCHURCOMPLEMENT`
2089: @*/
2090: PetscErrorCode MatIsShell(Mat mat, PetscBool *flg)
2091: {
2092:   PetscFunctionBegin;
2095:   *flg = (PetscBool)(mat->ops->destroy == MatDestroy_Shell);
2096:   PetscFunctionReturn(PETSC_SUCCESS);
2097: }