Actual source code: shell.c

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

176:          On input Y already contains A*x

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

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

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

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

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

222:   Not Collective

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

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

230:   Level: advanced

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

823:   Logically Collective; No Fortran Support

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

834:   Level: advanced

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1172: static PetscErrorCode MatShift_Shell(Mat Y, PetscScalar a)
1173: {
1174:   Mat_Shell *shell = (Mat_Shell *)Y->data;
1175:   PetscBool  flg;

1177:   PetscFunctionBegin;
1178:   PetscCall(MatHasCongruentLayouts(Y, &flg));
1179:   PetscCheck(flg, PetscObjectComm((PetscObject)Y), PETSC_ERR_SUP, "Cannot shift shell matrix if it is not congruent");
1180:   if (shell->left || shell->right) {
1181:     if (!shell->dshift) {
1182:       PetscCall(VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift));
1183:       PetscCall(VecSet(shell->dshift, a));
1184:     } else {
1185:       if (shell->left) PetscCall(VecPointwiseMult(shell->dshift, shell->dshift, shell->left));
1186:       if (shell->right) PetscCall(VecPointwiseMult(shell->dshift, shell->dshift, shell->right));
1187:       PetscCall(VecShift(shell->dshift, a));
1188:     }
1189:     if (shell->left) PetscCall(VecPointwiseDivide(shell->dshift, shell->dshift, shell->left));
1190:     if (shell->right) PetscCall(VecPointwiseDivide(shell->dshift, shell->dshift, shell->right));
1191:   } else shell->vshift += a;
1192:   if (shell->zrows) PetscCall(VecShift(shell->zvals, a));
1193:   PetscFunctionReturn(PETSC_SUCCESS);
1194: }

1196: static PetscErrorCode MatDiagonalSet_Shell_Private(Mat A, Vec D, PetscScalar s)
1197: {
1198:   Mat_Shell *shell = (Mat_Shell *)A->data;

1200:   PetscFunctionBegin;
1201:   if (!shell->dshift) PetscCall(VecDuplicate(D, &shell->dshift));
1202:   if (shell->left || shell->right) {
1203:     if (!shell->right_work) PetscCall(VecDuplicate(shell->left ? shell->left : shell->right, &shell->right_work));
1204:     if (shell->left && shell->right) {
1205:       PetscCall(VecPointwiseDivide(shell->right_work, D, shell->left));
1206:       PetscCall(VecPointwiseDivide(shell->right_work, shell->right_work, shell->right));
1207:     } else if (shell->left) {
1208:       PetscCall(VecPointwiseDivide(shell->right_work, D, shell->left));
1209:     } else {
1210:       PetscCall(VecPointwiseDivide(shell->right_work, D, shell->right));
1211:     }
1212:     PetscCall(VecAXPY(shell->dshift, s, shell->right_work));
1213:   } else {
1214:     PetscCall(VecAXPY(shell->dshift, s, D));
1215:   }
1216:   PetscFunctionReturn(PETSC_SUCCESS);
1217: }

1219: static PetscErrorCode MatDiagonalSet_Shell(Mat A, Vec D, InsertMode ins)
1220: {
1221:   Mat_Shell *shell = (Mat_Shell *)A->data;
1222:   Vec        d;
1223:   PetscBool  flg;

1225:   PetscFunctionBegin;
1226:   PetscCall(MatHasCongruentLayouts(A, &flg));
1227:   PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Cannot diagonal set or shift shell matrix if it is not congruent");
1228:   if (ins == INSERT_VALUES) {
1229:     PetscCall(VecDuplicate(D, &d));
1230:     PetscCall(MatGetDiagonal(A, d));
1231:     PetscCall(MatDiagonalSet_Shell_Private(A, d, -1.));
1232:     PetscCall(MatDiagonalSet_Shell_Private(A, D, 1.));
1233:     PetscCall(VecDestroy(&d));
1234:     if (shell->zrows) PetscCall(VecCopy(D, shell->zvals));
1235:   } else {
1236:     PetscCall(MatDiagonalSet_Shell_Private(A, D, 1.));
1237:     if (shell->zrows) PetscCall(VecAXPY(shell->zvals, 1.0, D));
1238:   }
1239:   PetscFunctionReturn(PETSC_SUCCESS);
1240: }

1242: static PetscErrorCode MatScale_Shell(Mat Y, PetscScalar a)
1243: {
1244:   Mat_Shell *shell = (Mat_Shell *)Y->data;

1246:   PetscFunctionBegin;
1247:   shell->vscale *= a;
1248:   shell->vshift *= a;
1249:   if (shell->dshift) PetscCall(VecScale(shell->dshift, a));
1250:   shell->axpy_vscale *= a;
1251:   if (shell->zrows) PetscCall(VecScale(shell->zvals, a));
1252:   PetscFunctionReturn(PETSC_SUCCESS);
1253: }

1255: static PetscErrorCode MatDiagonalScale_Shell(Mat Y, Vec left, Vec right)
1256: {
1257:   Mat_Shell *shell = (Mat_Shell *)Y->data;

1259:   PetscFunctionBegin;
1260:   if (left) {
1261:     if (!shell->left) {
1262:       PetscCall(VecDuplicate(left, &shell->left));
1263:       PetscCall(VecCopy(left, shell->left));
1264:     } else {
1265:       PetscCall(VecPointwiseMult(shell->left, shell->left, left));
1266:     }
1267:     if (shell->zrows) PetscCall(VecPointwiseMult(shell->zvals, shell->zvals, left));
1268:   }
1269:   if (right) {
1270:     if (!shell->right) {
1271:       PetscCall(VecDuplicate(right, &shell->right));
1272:       PetscCall(VecCopy(right, shell->right));
1273:     } else {
1274:       PetscCall(VecPointwiseMult(shell->right, shell->right, right));
1275:     }
1276:     if (shell->zrows) {
1277:       if (!shell->left_work) PetscCall(MatCreateVecs(Y, NULL, &shell->left_work));
1278:       PetscCall(VecSet(shell->zvals_w, 1.0));
1279:       PetscCall(VecScatterBegin(shell->zvals_sct_c, right, shell->zvals_w, INSERT_VALUES, SCATTER_FORWARD));
1280:       PetscCall(VecScatterEnd(shell->zvals_sct_c, right, shell->zvals_w, INSERT_VALUES, SCATTER_FORWARD));
1281:       PetscCall(VecPointwiseMult(shell->zvals, shell->zvals, shell->zvals_w));
1282:     }
1283:   }
1284:   if (shell->axpy) PetscCall(MatDiagonalScale(shell->axpy, left, right));
1285:   PetscFunctionReturn(PETSC_SUCCESS);
1286: }

1288: PETSC_INTERN PetscErrorCode MatAssemblyEnd_Shell(Mat Y, MatAssemblyType t)
1289: {
1290:   Mat_Shell *shell = (Mat_Shell *)Y->data;

1292:   PetscFunctionBegin;
1293:   if (t == MAT_FINAL_ASSEMBLY) {
1294:     shell->vshift      = 0.0;
1295:     shell->vscale      = 1.0;
1296:     shell->axpy_vscale = 0.0;
1297:     shell->axpy_state  = 0;
1298:     PetscCall(VecDestroy(&shell->dshift));
1299:     PetscCall(VecDestroy(&shell->left));
1300:     PetscCall(VecDestroy(&shell->right));
1301:     PetscCall(MatDestroy(&shell->axpy));
1302:     PetscCall(VecDestroy(&shell->axpy_left));
1303:     PetscCall(VecDestroy(&shell->axpy_right));
1304:     PetscCall(VecScatterDestroy(&shell->zvals_sct_c));
1305:     PetscCall(VecScatterDestroy(&shell->zvals_sct_r));
1306:     PetscCall(ISDestroy(&shell->zrows));
1307:     PetscCall(ISDestroy(&shell->zcols));
1308:   }
1309:   PetscFunctionReturn(PETSC_SUCCESS);
1310: }

1312: static PetscErrorCode MatMissingDiagonal_Shell(Mat A, PetscBool *missing, PetscInt *d)
1313: {
1314:   PetscFunctionBegin;
1315:   *missing = PETSC_FALSE;
1316:   PetscFunctionReturn(PETSC_SUCCESS);
1317: }

1319: static PetscErrorCode MatAXPY_Shell(Mat Y, PetscScalar a, Mat X, MatStructure str)
1320: {
1321:   Mat_Shell *shell = (Mat_Shell *)Y->data;

1323:   PetscFunctionBegin;
1324:   if (X == Y) {
1325:     PetscCall(MatScale(Y, 1.0 + a));
1326:     PetscFunctionReturn(PETSC_SUCCESS);
1327:   }
1328:   if (!shell->axpy) {
1329:     PetscCall(MatConvertFrom_Shell(X, MATSHELL, MAT_INITIAL_MATRIX, &shell->axpy));
1330:     shell->axpy_vscale = a;
1331:     PetscCall(PetscObjectStateGet((PetscObject)X, &shell->axpy_state));
1332:   } else {
1333:     PetscCall(MatAXPY(shell->axpy, a / shell->axpy_vscale, X, str));
1334:   }
1335:   PetscFunctionReturn(PETSC_SUCCESS);
1336: }

1338: static struct _MatOps MatOps_Values = {NULL,
1339:                                        NULL,
1340:                                        NULL,
1341:                                        NULL,
1342:                                        /* 4*/ MatMultAdd_Shell,
1343:                                        NULL,
1344:                                        MatMultTransposeAdd_Shell,
1345:                                        NULL,
1346:                                        NULL,
1347:                                        NULL,
1348:                                        /*10*/ NULL,
1349:                                        NULL,
1350:                                        NULL,
1351:                                        NULL,
1352:                                        NULL,
1353:                                        /*15*/ NULL,
1354:                                        NULL,
1355:                                        NULL,
1356:                                        MatDiagonalScale_Shell,
1357:                                        NULL,
1358:                                        /*20*/ NULL,
1359:                                        MatAssemblyEnd_Shell,
1360:                                        NULL,
1361:                                        NULL,
1362:                                        /*24*/ MatZeroRows_Shell,
1363:                                        NULL,
1364:                                        NULL,
1365:                                        NULL,
1366:                                        NULL,
1367:                                        /*29*/ NULL,
1368:                                        NULL,
1369:                                        NULL,
1370:                                        NULL,
1371:                                        NULL,
1372:                                        /*34*/ MatDuplicate_Shell,
1373:                                        NULL,
1374:                                        NULL,
1375:                                        NULL,
1376:                                        NULL,
1377:                                        /*39*/ MatAXPY_Shell,
1378:                                        NULL,
1379:                                        NULL,
1380:                                        NULL,
1381:                                        MatCopy_Shell,
1382:                                        /*44*/ NULL,
1383:                                        MatScale_Shell,
1384:                                        MatShift_Shell,
1385:                                        MatDiagonalSet_Shell,
1386:                                        MatZeroRowsColumns_Shell,
1387:                                        /*49*/ NULL,
1388:                                        NULL,
1389:                                        NULL,
1390:                                        NULL,
1391:                                        NULL,
1392:                                        /*54*/ NULL,
1393:                                        NULL,
1394:                                        NULL,
1395:                                        NULL,
1396:                                        NULL,
1397:                                        /*59*/ NULL,
1398:                                        MatDestroy_Shell,
1399:                                        NULL,
1400:                                        MatConvertFrom_Shell,
1401:                                        NULL,
1402:                                        /*64*/ NULL,
1403:                                        NULL,
1404:                                        NULL,
1405:                                        NULL,
1406:                                        NULL,
1407:                                        /*69*/ NULL,
1408:                                        NULL,
1409:                                        MatConvert_Shell,
1410:                                        NULL,
1411:                                        NULL,
1412:                                        /*74*/ NULL,
1413:                                        NULL,
1414:                                        NULL,
1415:                                        NULL,
1416:                                        NULL,
1417:                                        /*79*/ NULL,
1418:                                        NULL,
1419:                                        NULL,
1420:                                        NULL,
1421:                                        NULL,
1422:                                        /*84*/ NULL,
1423:                                        NULL,
1424:                                        NULL,
1425:                                        NULL,
1426:                                        NULL,
1427:                                        /*89*/ NULL,
1428:                                        NULL,
1429:                                        NULL,
1430:                                        NULL,
1431:                                        NULL,
1432:                                        /*94*/ NULL,
1433:                                        NULL,
1434:                                        NULL,
1435:                                        NULL,
1436:                                        NULL,
1437:                                        /*99*/ NULL,
1438:                                        NULL,
1439:                                        NULL,
1440:                                        NULL,
1441:                                        NULL,
1442:                                        /*104*/ NULL,
1443:                                        NULL,
1444:                                        NULL,
1445:                                        NULL,
1446:                                        NULL,
1447:                                        /*109*/ NULL,
1448:                                        NULL,
1449:                                        NULL,
1450:                                        NULL,
1451:                                        MatMissingDiagonal_Shell,
1452:                                        /*114*/ NULL,
1453:                                        NULL,
1454:                                        NULL,
1455:                                        NULL,
1456:                                        NULL,
1457:                                        /*119*/ NULL,
1458:                                        NULL,
1459:                                        NULL,
1460:                                        MatMultHermitianTransposeAdd_Shell,
1461:                                        NULL,
1462:                                        /*124*/ NULL,
1463:                                        NULL,
1464:                                        NULL,
1465:                                        NULL,
1466:                                        NULL,
1467:                                        /*129*/ NULL,
1468:                                        NULL,
1469:                                        NULL,
1470:                                        NULL,
1471:                                        NULL,
1472:                                        /*134*/ NULL,
1473:                                        NULL,
1474:                                        NULL,
1475:                                        NULL,
1476:                                        NULL,
1477:                                        /*139*/ NULL,
1478:                                        NULL,
1479:                                        NULL,
1480:                                        NULL,
1481:                                        NULL,
1482:                                        /*144*/ NULL,
1483:                                        NULL,
1484:                                        NULL,
1485:                                        NULL,
1486:                                        NULL,
1487:                                        NULL,
1488:                                        /*150*/ NULL,
1489:                                        NULL,
1490:                                        NULL};

1492: static PetscErrorCode MatShellSetContext_Shell(Mat mat, void *ctx)
1493: {
1494:   Mat_Shell *shell = (Mat_Shell *)mat->data;

1496:   PetscFunctionBegin;
1497:   if (ctx) {
1498:     PetscContainer ctxcontainer;
1499:     PetscCall(PetscContainerCreate(PetscObjectComm((PetscObject)mat), &ctxcontainer));
1500:     PetscCall(PetscContainerSetPointer(ctxcontainer, ctx));
1501:     PetscCall(PetscObjectCompose((PetscObject)mat, "MatShell ctx", (PetscObject)ctxcontainer));
1502:     shell->ctxcontainer = ctxcontainer;
1503:     PetscCall(PetscContainerDestroy(&ctxcontainer));
1504:   } else {
1505:     PetscCall(PetscObjectCompose((PetscObject)mat, "MatShell ctx", NULL));
1506:     shell->ctxcontainer = NULL;
1507:   }
1508:   PetscFunctionReturn(PETSC_SUCCESS);
1509: }

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

1515:   PetscFunctionBegin;
1516:   if (shell->ctxcontainer) PetscCall(PetscContainerSetUserDestroy(shell->ctxcontainer, f));
1517:   PetscFunctionReturn(PETSC_SUCCESS);
1518: }

1520: PetscErrorCode MatShellSetContext_Immutable(Mat mat, void *ctx)
1521: {
1522:   PetscFunctionBegin;
1523:   SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "Cannot call MatShellSetContext() for a %s, it is used internally by the structure", ((PetscObject)mat)->type_name);
1524:   PetscFunctionReturn(PETSC_SUCCESS);
1525: }

1527: PetscErrorCode MatShellSetContextDestroy_Immutable(Mat mat, PetscErrorCode (*f)(void *))
1528: {
1529:   PetscFunctionBegin;
1530:   SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "Cannot call MatShellSetContextDestroy() for a %s, it is used internally by the structure", ((PetscObject)mat)->type_name);
1531:   PetscFunctionReturn(PETSC_SUCCESS);
1532: }

1534: PetscErrorCode MatShellSetManageScalingShifts_Immutable(Mat mat)
1535: {
1536:   PetscFunctionBegin;
1537:   SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "Cannot call MatShellSetManageScalingShifts() for a %s, it is used internally by the structure", ((PetscObject)mat)->type_name);
1538:   PetscFunctionReturn(PETSC_SUCCESS);
1539: }

1541: static PetscErrorCode MatShellSetVecType_Shell(Mat mat, VecType vtype)
1542: {
1543:   PetscFunctionBegin;
1544:   PetscCall(PetscFree(mat->defaultvectype));
1545:   PetscCall(PetscStrallocpy(vtype, (char **)&mat->defaultvectype));
1546:   PetscFunctionReturn(PETSC_SUCCESS);
1547: }

1549: static PetscErrorCode MatShellSetManageScalingShifts_Shell(Mat A)
1550: {
1551:   Mat_Shell *shell = (Mat_Shell *)A->data;

1553:   PetscFunctionBegin;
1554:   shell->managescalingshifts = PETSC_FALSE;
1555:   A->ops->diagonalset        = NULL;
1556:   A->ops->diagonalscale      = NULL;
1557:   A->ops->scale              = NULL;
1558:   A->ops->shift              = NULL;
1559:   A->ops->axpy               = NULL;
1560:   PetscFunctionReturn(PETSC_SUCCESS);
1561: }

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

1567:   PetscFunctionBegin;
1568:   PetscCheck(shell->managescalingshifts, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Cannot get scaling and shifts if MatShellSetManageScalingShifts() has been called");
1569:   if (vshift == MAT_SHELL_NOT_ALLOWED) PetscCheck(shell->vshift == 0.0, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Non-trivial member in parent MatShell: vshift != 0.0, set via MatShift()");
1570:   else if (vshift) *vshift = shell->vshift;
1571:   if (vscale == MAT_SHELL_NOT_ALLOWED) PetscCheck(shell->vscale == 1.0, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Non-trivial member in parent MatShell: vscale != 1.0, set via MatScale()");
1572:   else if (vscale) *vscale = shell->vscale;
1573:   if (dshift == MAT_SHELL_NOT_ALLOWED) PetscCheck(!shell->dshift, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Non-trivial member in parent MatShell: dshift, set via MatDiagonalSet()");
1574:   else if (dshift) *dshift = shell->dshift;
1575:   if (left == MAT_SHELL_NOT_ALLOWED) PetscCheck(!shell->left, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Non-trivial member in parent MatShell: left, set via MatDiagonalScale()");
1576:   else if (left) *left = shell->left;
1577:   if (right == MAT_SHELL_NOT_ALLOWED) PetscCheck(!shell->right, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Non-trivial member in parent MatShell: right, set via MatDiagonalScale()");
1578:   else if (right) *right = shell->right;
1579:   if (axpy == MAT_SHELL_NOT_ALLOWED) PetscCheck(!shell->axpy, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Non-trivial member in parent MatShell: axpy, set via MatAXPY()");
1580:   else if (axpy) *axpy = shell->axpy;
1581:   if (zrows == MAT_SHELL_NOT_ALLOWED) PetscCheck(!shell->zrows, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Non-trivial member in parent MatShell: zrows, set via MatZeroRows()");
1582:   else if (zrows) *zrows = shell->zrows;
1583:   if (zcols == MAT_SHELL_NOT_ALLOWED) PetscCheck(!shell->zcols, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Non-trivial member in parent MatShell: zcols, set via MatZeroRowsColumns()");
1584:   else if (zcols) *zcols = shell->zcols;
1585:   PetscFunctionReturn(PETSC_SUCCESS);
1586: }

1588: static PetscErrorCode MatShellSetOperation_Shell(Mat mat, MatOperation op, void (*f)(void))
1589: {
1590:   Mat_Shell *shell = (Mat_Shell *)mat->data;

1592:   PetscFunctionBegin;
1593:   switch (op) {
1594:   case MATOP_DESTROY:
1595:     shell->ops->destroy = (PetscErrorCode(*)(Mat))f;
1596:     break;
1597:   case MATOP_VIEW:
1598:     if (!mat->ops->viewnative) mat->ops->viewnative = mat->ops->view;
1599:     mat->ops->view = (PetscErrorCode(*)(Mat, PetscViewer))f;
1600:     break;
1601:   case MATOP_COPY:
1602:     shell->ops->copy = (PetscErrorCode(*)(Mat, Mat, MatStructure))f;
1603:     break;
1604:   case MATOP_DIAGONAL_SET:
1605:   case MATOP_DIAGONAL_SCALE:
1606:   case MATOP_SHIFT:
1607:   case MATOP_SCALE:
1608:   case MATOP_AXPY:
1609:   case MATOP_ZERO_ROWS:
1610:   case MATOP_ZERO_ROWS_COLUMNS:
1611:     PetscCheck(!shell->managescalingshifts, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "MATSHELL is managing scalings and shifts, see MatShellSetManageScalingShifts()");
1612:     (((void (**)(void))mat->ops)[op]) = f;
1613:     break;
1614:   case MATOP_GET_DIAGONAL:
1615:     if (shell->managescalingshifts) {
1616:       shell->ops->getdiagonal = (PetscErrorCode(*)(Mat, Vec))f;
1617:       mat->ops->getdiagonal   = MatGetDiagonal_Shell;
1618:     } else {
1619:       shell->ops->getdiagonal = NULL;
1620:       mat->ops->getdiagonal   = (PetscErrorCode(*)(Mat, Vec))f;
1621:     }
1622:     break;
1623:   case MATOP_MULT:
1624:     if (shell->managescalingshifts) {
1625:       shell->ops->mult = (PetscErrorCode(*)(Mat, Vec, Vec))f;
1626:       mat->ops->mult   = MatMult_Shell;
1627:     } else {
1628:       shell->ops->mult = NULL;
1629:       mat->ops->mult   = (PetscErrorCode(*)(Mat, Vec, Vec))f;
1630:     }
1631:     break;
1632:   case MATOP_MULT_TRANSPOSE:
1633:     if (shell->managescalingshifts) {
1634:       shell->ops->multtranspose = (PetscErrorCode(*)(Mat, Vec, Vec))f;
1635:       mat->ops->multtranspose   = MatMultTranspose_Shell;
1636:     } else {
1637:       shell->ops->multtranspose = NULL;
1638:       mat->ops->multtranspose   = (PetscErrorCode(*)(Mat, Vec, Vec))f;
1639:     }
1640:     break;
1641:   case MATOP_MULT_HERMITIAN_TRANSPOSE:
1642:     if (shell->managescalingshifts) {
1643:       shell->ops->multhermitiantranspose = (PetscErrorCode(*)(Mat, Vec, Vec))f;
1644:       mat->ops->multhermitiantranspose   = MatMultHermitianTranspose_Shell;
1645:     } else {
1646:       shell->ops->multhermitiantranspose = NULL;
1647:       mat->ops->multhermitiantranspose   = (PetscErrorCode(*)(Mat, Vec, Vec))f;
1648:     }
1649:     break;
1650:   default:
1651:     (((void (**)(void))mat->ops)[op]) = f;
1652:     break;
1653:   }
1654:   PetscFunctionReturn(PETSC_SUCCESS);
1655: }

1657: static PetscErrorCode MatShellGetOperation_Shell(Mat mat, MatOperation op, void (**f)(void))
1658: {
1659:   Mat_Shell *shell = (Mat_Shell *)mat->data;

1661:   PetscFunctionBegin;
1662:   switch (op) {
1663:   case MATOP_DESTROY:
1664:     *f = (void (*)(void))shell->ops->destroy;
1665:     break;
1666:   case MATOP_VIEW:
1667:     *f = (void (*)(void))mat->ops->view;
1668:     break;
1669:   case MATOP_COPY:
1670:     *f = (void (*)(void))shell->ops->copy;
1671:     break;
1672:   case MATOP_DIAGONAL_SET:
1673:   case MATOP_DIAGONAL_SCALE:
1674:   case MATOP_SHIFT:
1675:   case MATOP_SCALE:
1676:   case MATOP_AXPY:
1677:   case MATOP_ZERO_ROWS:
1678:   case MATOP_ZERO_ROWS_COLUMNS:
1679:     *f = (((void (**)(void))mat->ops)[op]);
1680:     break;
1681:   case MATOP_GET_DIAGONAL:
1682:     if (shell->ops->getdiagonal) *f = (void (*)(void))shell->ops->getdiagonal;
1683:     else *f = (((void (**)(void))mat->ops)[op]);
1684:     break;
1685:   case MATOP_MULT:
1686:     if (shell->ops->mult) *f = (void (*)(void))shell->ops->mult;
1687:     else *f = (((void (**)(void))mat->ops)[op]);
1688:     break;
1689:   case MATOP_MULT_TRANSPOSE:
1690:     if (shell->ops->multtranspose) *f = (void (*)(void))shell->ops->multtranspose;
1691:     else *f = (((void (**)(void))mat->ops)[op]);
1692:     break;
1693:   case MATOP_MULT_HERMITIAN_TRANSPOSE:
1694:     if (shell->ops->multhermitiantranspose) *f = (void (*)(void))shell->ops->multhermitiantranspose;
1695:     else *f = (((void (**)(void))mat->ops)[op]);
1696:     break;
1697:   default:
1698:     *f = (((void (**)(void))mat->ops)[op]);
1699:   }
1700:   PetscFunctionReturn(PETSC_SUCCESS);
1701: }

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

1706:   Level: advanced

1708: .seealso: [](ch_matrices), `Mat`, `MatCreateShell()`
1709: M*/

1711: PETSC_EXTERN PetscErrorCode MatCreate_Shell(Mat A)
1712: {
1713:   Mat_Shell *b;

1715:   PetscFunctionBegin;
1716:   PetscCall(PetscNew(&b));
1717:   A->data   = (void *)b;
1718:   A->ops[0] = MatOps_Values;

1720:   b->ctxcontainer        = NULL;
1721:   b->vshift              = 0.0;
1722:   b->vscale              = 1.0;
1723:   b->managescalingshifts = PETSC_TRUE;
1724:   A->assembled           = PETSC_TRUE;
1725:   A->preallocated        = PETSC_FALSE;

1727:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellGetContext_C", MatShellGetContext_Shell));
1728:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetContext_C", MatShellSetContext_Shell));
1729:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetContextDestroy_C", MatShellSetContextDestroy_Shell));
1730:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetVecType_C", MatShellSetVecType_Shell));
1731:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetManageScalingShifts_C", MatShellSetManageScalingShifts_Shell));
1732:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellGetScalingShifts_C", MatShellGetScalingShifts_Shell));
1733:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetOperation_C", MatShellSetOperation_Shell));
1734:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellGetOperation_C", MatShellGetOperation_Shell));
1735:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetMatProductOperation_C", MatShellSetMatProductOperation_Shell));
1736:   PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATSHELL));
1737:   PetscFunctionReturn(PETSC_SUCCESS);
1738: }

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

1744:   Collective

1746:   Input Parameters:
1747: + comm - MPI communicator
1748: . m    - number of local rows (or `PETSC_DECIDE` to have calculated if `M` is given)
1749: . n    - number of local columns (or `PETSC_DECIDE` to have calculated if `N` is given)
1750: . M    - number of global rows (may be `PETSC_DETERMINE` to have calculated if `m` is given)
1751: . N    - number of global columns (may be `PETSC_DETERMINE` to have calculated if `n` is given)
1752: - ctx  - pointer to data needed by the shell matrix routines

1754:   Output Parameter:
1755: . A - the matrix

1757:   Level: advanced

1759:   Example Usage:
1760: .vb
1761:   extern PetscErrorCode mult(Mat, Vec, Vec);

1763:   MatCreateShell(comm, m, n, M, N, ctx, &mat);
1764:   MatShellSetOperation(mat, MATOP_MULT, (void(*)(void))mult);
1765:   // Use matrix for operations that have been set
1766:   MatDestroy(mat);
1767: .ve

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

1774:   PETSc requires that matrices and vectors being used for certain
1775:   operations are partitioned accordingly.  For example, when
1776:   creating a shell matrix, `A`, that supports parallel matrix-vector
1777:   products using `MatMult`(A,x,y) the user should set the number
1778:   of local matrix rows to be the number of local elements of the
1779:   corresponding result vector, y. Note that this is information is
1780:   required for use of the matrix interface routines, even though
1781:   the shell matrix may not actually be physically partitioned.
1782:   For example,

1784: .vb
1785:      Vec x, y
1786:      extern PetscErrorCode mult(Mat,Vec,Vec);
1787:      Mat A

1789:      VecCreateMPI(comm,PETSC_DECIDE,M,&y);
1790:      VecCreateMPI(comm,PETSC_DECIDE,N,&x);
1791:      VecGetLocalSize(y,&m);
1792:      VecGetLocalSize(x,&n);
1793:      MatCreateShell(comm,m,n,M,N,ctx,&A);
1794:      MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
1795:      MatMult(A,x,y);
1796:      MatDestroy(&A);
1797:      VecDestroy(&y);
1798:      VecDestroy(&x);
1799: .ve

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

1804:   Developer Notes:
1805:   For rectangular matrices do all the scalings and shifts make sense?

1807:   Regarding shifting and scaling. The general form is

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

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

1815:   A is the user provided function.

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

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

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

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

1832: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatShellSetOperation()`, `MatHasOperation()`, `MatShellGetContext()`, `MatShellSetContext()`, `MatShellSetManageScalingShifts()`, `MatShellSetMatProductOperation()`
1833: @*/
1834: PetscErrorCode MatCreateShell(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, PetscInt N, void *ctx, Mat *A)
1835: {
1836:   PetscFunctionBegin;
1837:   PetscCall(MatCreate(comm, A));
1838:   PetscCall(MatSetSizes(*A, m, n, M, N));
1839:   PetscCall(MatSetType(*A, MATSHELL));
1840:   PetscCall(MatShellSetContext(*A, ctx));
1841:   PetscCall(MatSetUp(*A));
1842:   PetscFunctionReturn(PETSC_SUCCESS);
1843: }

1845: /*@
1846:   MatShellSetContext - sets the context for a `MATSHELL` shell matrix

1848:   Logically Collective

1850:   Input Parameters:
1851: + mat - the `MATSHELL` shell matrix
1852: - ctx - the context

1854:   Level: advanced

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

1860: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`
1861: @*/
1862: PetscErrorCode MatShellSetContext(Mat mat, void *ctx)
1863: {
1864:   PetscFunctionBegin;
1866:   PetscTryMethod(mat, "MatShellSetContext_C", (Mat, void *), (mat, ctx));
1867:   PetscFunctionReturn(PETSC_SUCCESS);
1868: }

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

1873:   Logically Collective

1875:   Input Parameters:
1876: + mat - the shell matrix
1877: - f   - the context destroy function

1879:   Level: advanced

1881:   Note:
1882:   If the `MatShell` is never duplicated, the behavior of this function is equivalent
1883:   to `MatShellSetOperation`(`Mat`,`MATOP_DESTROY`,f). However, `MatShellSetContextDestroy()`
1884:   ensures proper reference counting for the user provided context data in the case that
1885:   the `MATSHELL` is duplicated.

1887: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellSetContext()`
1888: @*/
1889: PetscErrorCode MatShellSetContextDestroy(Mat mat, PetscErrorCode (*f)(void *))
1890: {
1891:   PetscFunctionBegin;
1893:   PetscTryMethod(mat, "MatShellSetContextDestroy_C", (Mat, PetscErrorCode(*)(void *)), (mat, f));
1894:   PetscFunctionReturn(PETSC_SUCCESS);
1895: }

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

1900:   Logically Collective

1902:   Input Parameters:
1903: + mat   - the `MATSHELL` shell matrix
1904: - vtype - type to use for creating vectors

1906:   Level: advanced

1908: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateVecs()`
1909: @*/
1910: PetscErrorCode MatShellSetVecType(Mat mat, VecType vtype)
1911: {
1912:   PetscFunctionBegin;
1913:   PetscTryMethod(mat, "MatShellSetVecType_C", (Mat, VecType), (mat, vtype));
1914:   PetscFunctionReturn(PETSC_SUCCESS);
1915: }

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

1921:   Logically Collective

1923:   Input Parameter:
1924: . A - the `MATSHELL` shell matrix

1926:   Level: advanced

1928: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellSetContext()`, `MatShellSetOperation()`
1929: @*/
1930: PetscErrorCode MatShellSetManageScalingShifts(Mat A)
1931: {
1932:   PetscFunctionBegin;
1934:   PetscTryMethod(A, "MatShellSetManageScalingShifts_C", (Mat), (A));
1935:   PetscFunctionReturn(PETSC_SUCCESS);
1936: }

1938: // PetscClangLinter pragma disable: -fdoc-internal-linkage
1939: /*@C
1940:   MatShellGetScalingShifts - Gets members of a `MATSHELL` used internally for scaling and
1941:   shifting the `Mat` or calling `MatAXPY()`, `MatZeroRows()`, or `MatZeroRowsColumns()` with it

1943:   Logically Collective

1945:   Input Parameter:
1946: . A - the `MATSHELL` shell matrix

1948:   Output Parameters:
1949: + vshift - `PetscScalar` pointer (can be `NULL`), see `MatShift()`, or `MAT_SHELL_NOT_ALLOWED` if the internal shift should be 0
1950: . vscale - `PetscScalar` pointer (can be `NULL`), see `MatScale()`, or `MAT_SHELL_NOT_ALLOWED` if the internal scaling should be 1
1951: . dshift - `Vec` pointer (can be `NULL`), see `MatDiagonalSet()`, or `MAT_SHELL_NOT_ALLOWED` if the internal shift should be `NULL`
1952: . left   - `Vec` pointer (can be `NULL`), see `MatDiagonalScale()`, or `MAT_SHELL_NOT_ALLOWED` if the internal scaling should be `NULL`
1953: . right  - `Vec` pointer (can be `NULL`), see `MatDiagonalScale()`, or `MAT_SHELL_NOT_ALLOWED` if the internal scaling should be `NULL`
1954: . axpy   - `Mat` pointer (can be `NULL`), or `MAT_SHELL_NOT_ALLOWED` if `MatAXPY()` should have not been called on `A`
1955: . zrows  - `Vec` pointer (can be `NULL`), or `MAT_SHELL_NOT_ALLOWED` if `MatZeroRows()` should have not been called on `A`
1956: - zcols  - `Vec` pointer (can be `NULL`), or `MAT_SHELL_NOT_ALLOWED` if `MatZeroRowsColumns()` should have not been called on `A`

1958:   Level: advanced

1960:   Developer Notes:
1961:   This is mostly useful to check for corner-cases in `MatType` deriving from
1962:   `MATSHELL`, e.g, `MATCOMPOSITE` or `MATVIRTUALTRANSPOSE`, since scaling and
1963:   shifts often require extra work which is not always implemented.

1965: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShift()`, `MatScale()`, `MatDiagonalSet()`, `MatDiagonalScale()`, `MatAXPY()`, `MatZeroRows()`, `MatZeroRowsColumns()`, `MatShellSetManageScalingShifts()`
1966: @*/
1967: PetscErrorCode MatShellGetScalingShifts(Mat A, PetscScalar *vshift, PetscScalar *vscale, Vec *dshift, Vec *left, Vec *right, Mat *axpy, IS *zrows, IS *zcols)
1968: {
1969:   PetscFunctionBegin;
1971:   PetscTryMethod(A, "MatShellGetScalingShifts_C", (Mat, PetscScalar *, PetscScalar *, Vec *, Vec *, Vec *, Mat *, IS *, IS *), (A, vshift, vscale, dshift, left, right, axpy, zrows, zcols));
1972:   PetscFunctionReturn(PETSC_SUCCESS);
1973: }

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

1978:   Logically Collective; No Fortran Support

1980:   Input Parameters:
1981: + mat  - the `MATSHELL` shell matrix
1982: . f    - the function
1983: . base - differences are computed around this vector, see `MatMFFDSetBase()`, for Jacobians this is the point at which the Jacobian is being evaluated
1984: - ctx  - an optional context for the function

1986:   Output Parameter:
1987: . flg - `PETSC_TRUE` if the multiply is likely correct

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

1992:   Level: advanced

1994: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMultTranspose()`
1995: @*/
1996: PetscErrorCode MatShellTestMult(Mat mat, PetscErrorCode (*f)(void *, Vec, Vec), Vec base, void *ctx, PetscBool *flg)
1997: {
1998:   PetscInt  m, n;
1999:   Mat       mf, Dmf, Dmat, Ddiff;
2000:   PetscReal Diffnorm, Dmfnorm;
2001:   PetscBool v = PETSC_FALSE, flag = PETSC_TRUE;

2003:   PetscFunctionBegin;
2005:   PetscCall(PetscOptionsHasName(NULL, ((PetscObject)mat)->prefix, "-mat_shell_test_mult_view", &v));
2006:   PetscCall(MatGetLocalSize(mat, &m, &n));
2007:   PetscCall(MatCreateMFFD(PetscObjectComm((PetscObject)mat), m, n, PETSC_DECIDE, PETSC_DECIDE, &mf));
2008:   PetscCall(MatMFFDSetFunction(mf, f, ctx));
2009:   PetscCall(MatMFFDSetBase(mf, base, NULL));

2011:   PetscCall(MatComputeOperator(mf, MATAIJ, &Dmf));
2012:   PetscCall(MatComputeOperator(mat, MATAIJ, &Dmat));

2014:   PetscCall(MatDuplicate(Dmat, MAT_COPY_VALUES, &Ddiff));
2015:   PetscCall(MatAXPY(Ddiff, -1.0, Dmf, DIFFERENT_NONZERO_PATTERN));
2016:   PetscCall(MatNorm(Ddiff, NORM_FROBENIUS, &Diffnorm));
2017:   PetscCall(MatNorm(Dmf, NORM_FROBENIUS, &Dmfnorm));
2018:   if (Diffnorm / Dmfnorm > 10 * PETSC_SQRT_MACHINE_EPSILON) {
2019:     flag = PETSC_FALSE;
2020:     if (v) {
2021:       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)));
2022:       PetscCall(MatViewFromOptions(Ddiff, (PetscObject)mat, "-mat_shell_test_mult_view"));
2023:       PetscCall(MatViewFromOptions(Dmf, (PetscObject)mat, "-mat_shell_test_mult_view"));
2024:       PetscCall(MatViewFromOptions(Dmat, (PetscObject)mat, "-mat_shell_test_mult_view"));
2025:     }
2026:   } else if (v) {
2027:     PetscCall(PetscPrintf(PetscObjectComm((PetscObject)mat), "MATSHELL and matrix-free multiple appear to produce the same results\n"));
2028:   }
2029:   if (flg) *flg = flag;
2030:   PetscCall(MatDestroy(&Ddiff));
2031:   PetscCall(MatDestroy(&mf));
2032:   PetscCall(MatDestroy(&Dmf));
2033:   PetscCall(MatDestroy(&Dmat));
2034:   PetscFunctionReturn(PETSC_SUCCESS);
2035: }

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

2040:   Logically Collective; No Fortran Support

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

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

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

2054:   Level: advanced

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

2066:   PetscFunctionBegin;
2068:   PetscCall(PetscOptionsHasName(NULL, ((PetscObject)mat)->prefix, "-mat_shell_test_mult_transpose_view", &v));
2069:   PetscCall(MatCreateVecs(mat, &x, &y));
2070:   PetscCall(VecDuplicate(y, &z));
2071:   PetscCall(MatGetLocalSize(mat, &m, &n));
2072:   PetscCall(MatGetSize(mat, &M, &N));
2073:   PetscCall(MatCreateMFFD(PetscObjectComm((PetscObject)mat), m, n, M, N, &mf));
2074:   PetscCall(MatMFFDSetFunction(mf, f, ctx));
2075:   PetscCall(MatMFFDSetBase(mf, base, NULL));
2076:   PetscCall(MatComputeOperator(mf, MATAIJ, &Dmf));
2077:   PetscCall(MatTranspose(Dmf, MAT_INPLACE_MATRIX, &Dmf));
2078:   PetscCall(MatComputeOperatorTranspose(mat, MATAIJ, &Dmat));

2080:   PetscCall(MatDuplicate(Dmat, MAT_COPY_VALUES, &Ddiff));
2081:   PetscCall(MatAXPY(Ddiff, -1.0, Dmf, DIFFERENT_NONZERO_PATTERN));
2082:   PetscCall(MatNorm(Ddiff, NORM_FROBENIUS, &Diffnorm));
2083:   PetscCall(MatNorm(Dmf, NORM_FROBENIUS, &Dmfnorm));
2084:   if (Diffnorm / Dmfnorm > 10 * PETSC_SQRT_MACHINE_EPSILON) {
2085:     flag = PETSC_FALSE;
2086:     if (v) {
2087:       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)));
2088:       PetscCall(MatViewFromOptions(Ddiff, (PetscObject)mat, "-mat_shell_test_mult_transpose_view"));
2089:       PetscCall(MatViewFromOptions(Dmf, (PetscObject)mat, "-mat_shell_test_mult_transpose_view"));
2090:       PetscCall(MatViewFromOptions(Dmat, (PetscObject)mat, "-mat_shell_test_mult_transpose_view"));
2091:     }
2092:   } else if (v) {
2093:     PetscCall(PetscPrintf(PetscObjectComm((PetscObject)mat), "MATSHELL transpose and matrix-free multiple appear to produce the same results\n"));
2094:   }
2095:   if (flg) *flg = flag;
2096:   PetscCall(MatDestroy(&mf));
2097:   PetscCall(MatDestroy(&Dmat));
2098:   PetscCall(MatDestroy(&Ddiff));
2099:   PetscCall(MatDestroy(&Dmf));
2100:   PetscCall(VecDestroy(&x));
2101:   PetscCall(VecDestroy(&y));
2102:   PetscCall(VecDestroy(&z));
2103:   PetscFunctionReturn(PETSC_SUCCESS);
2104: }

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

2109:   Logically Collective

2111:   Input Parameters:
2112: + mat - the `MATSHELL` shell matrix
2113: . op  - the name of the operation
2114: - g   - the function that provides the operation.

2116:   Level: advanced

2118:   Example Usage:
2119: .vb
2120:   extern PetscErrorCode usermult(Mat, Vec, Vec);

2122:   MatCreateShell(comm, m, n, M, N, ctx, &A);
2123:   MatShellSetOperation(A, MATOP_MULT, (void(*)(void))usermult);
2124: .ve

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

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

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

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

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

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

2152: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellSetContext()`, `MatSetOperation()`, `MatShellSetManageScalingShifts()`, `MatShellSetMatProductOperation()`
2153: @*/
2154: PetscErrorCode MatShellSetOperation(Mat mat, MatOperation op, void (*g)(void))
2155: {
2156:   PetscFunctionBegin;
2158:   PetscTryMethod(mat, "MatShellSetOperation_C", (Mat, MatOperation, void (*)(void)), (mat, op, g));
2159:   PetscFunctionReturn(PETSC_SUCCESS);
2160: }

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

2165:   Not Collective

2167:   Input Parameters:
2168: + mat - the `MATSHELL` shell matrix
2169: - op  - the name of the operation

2171:   Output Parameter:
2172: . g - the function that provides the operation.

2174:   Level: advanced

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

2182:   All user-provided functions have the same calling
2183:   sequence as the usual matrix interface routines, since they
2184:   are intended to be accessed via the usual matrix interface
2185:   routines, e.g.,
2186: $       MatMult(Mat, Vec, Vec) -> usermult(Mat, Vec, Vec)

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

2192: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellSetOperation()`, `MatShellSetContext()`
2193: @*/
2194: PetscErrorCode MatShellGetOperation(Mat mat, MatOperation op, void (**g)(void))
2195: {
2196:   PetscFunctionBegin;
2198:   PetscUseMethod(mat, "MatShellGetOperation_C", (Mat, MatOperation, void (**)(void)), (mat, op, g));
2199:   PetscFunctionReturn(PETSC_SUCCESS);
2200: }

2202: /*@
2203:   MatIsShell - Inquires if a matrix is derived from `MATSHELL`

2205:   Input Parameter:
2206: . mat - the matrix

2208:   Output Parameter:
2209: . flg - the Boolean value

2211:   Level: developer

2213: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MATMFFD`, `MatCreateShell()`, `MATTRANSPOSEVIRTUAL`, `MATSCHURCOMPLEMENT`
2214: @*/
2215: PetscErrorCode MatIsShell(Mat mat, PetscBool *flg)
2216: {
2217:   PetscFunctionBegin;
2219:   PetscAssertPointer(flg, 2);
2220:   *flg = (PetscBool)(mat->ops->destroy == MatDestroy_Shell);
2221:   PetscFunctionReturn(PETSC_SUCCESS);
2222: }