Actual source code: shell.c
1: /*
2: This provides a simple shell for Fortran (and C programmers) to
3: create a very simple matrix class for use with KSP without coding
4: much of anything.
5: */
7: #include <petsc/private/matimpl.h>
8: #include <petsc/private/vecimpl.h>
10: struct _MatShellOps {
11: /* 3 */ PetscErrorCode (*mult)(Mat, Vec, Vec);
12: /* 5 */ PetscErrorCode (*multtranspose)(Mat, Vec, Vec);
13: /* 17 */ PetscErrorCode (*getdiagonal)(Mat, Vec);
14: /* 43 */ PetscErrorCode (*copy)(Mat, Mat, MatStructure);
15: /* 60 */ PetscErrorCode (*destroy)(Mat);
16: };
18: struct _n_MatShellMatFunctionList {
19: PetscErrorCode (*symbolic)(Mat, Mat, Mat, void **);
20: PetscErrorCode (*numeric)(Mat, Mat, Mat, void *);
21: PetscErrorCode (*destroy)(void *);
22: MatProductType ptype;
23: char *composedname; /* string to identify routine with double dispatch */
24: char *resultname; /* result matrix type */
26: struct _n_MatShellMatFunctionList *next;
27: };
28: typedef struct _n_MatShellMatFunctionList *MatShellMatFunctionList;
30: typedef struct {
31: struct _MatShellOps ops[1];
33: /* The user will manage the scaling and shifts for the MATSHELL, not the default */
34: PetscBool managescalingshifts;
36: /* support for MatScale, MatShift and MatMultAdd */
37: PetscScalar vscale, vshift;
38: Vec dshift;
39: Vec left, right;
40: Vec left_work, right_work;
41: Vec left_add_work, right_add_work;
43: /* support for MatAXPY */
44: Mat axpy;
45: PetscScalar axpy_vscale;
46: Vec axpy_left, axpy_right;
47: PetscObjectState axpy_state;
49: /* support for ZeroRows/Columns operations */
50: IS zrows;
51: IS zcols;
52: Vec zvals;
53: Vec zvals_w;
54: VecScatter zvals_sct_r;
55: VecScatter zvals_sct_c;
57: /* MatMat operations */
58: MatShellMatFunctionList matmat;
60: /* user defined context */
61: PetscContainer ctxcontainer;
62: } Mat_Shell;
64: /*
65: Store and scale values on zeroed rows
66: xx = [x_1, 0], 0 on zeroed columns
67: */
68: static PetscErrorCode MatShellPreZeroRight(Mat A, Vec x, Vec *xx)
69: {
70: Mat_Shell *shell = (Mat_Shell *)A->data;
72: PetscFunctionBegin;
73: *xx = x;
74: if (shell->zrows) {
75: PetscCall(VecSet(shell->zvals_w, 0.0));
76: PetscCall(VecScatterBegin(shell->zvals_sct_c, x, shell->zvals_w, INSERT_VALUES, SCATTER_FORWARD));
77: PetscCall(VecScatterEnd(shell->zvals_sct_c, x, shell->zvals_w, INSERT_VALUES, SCATTER_FORWARD));
78: PetscCall(VecPointwiseMult(shell->zvals_w, shell->zvals_w, shell->zvals));
79: }
80: if (shell->zcols) {
81: if (!shell->right_work) PetscCall(MatCreateVecs(A, &shell->right_work, NULL));
82: PetscCall(VecCopy(x, shell->right_work));
83: PetscCall(VecISSet(shell->right_work, shell->zcols, 0.0));
84: *xx = shell->right_work;
85: }
86: PetscFunctionReturn(PETSC_SUCCESS);
87: }
89: /* Insert properly diagonally scaled values stored in MatShellPreZeroRight */
90: static PetscErrorCode MatShellPostZeroLeft(Mat A, Vec x)
91: {
92: Mat_Shell *shell = (Mat_Shell *)A->data;
94: PetscFunctionBegin;
95: if (shell->zrows) {
96: PetscCall(VecScatterBegin(shell->zvals_sct_r, shell->zvals_w, x, INSERT_VALUES, SCATTER_REVERSE));
97: PetscCall(VecScatterEnd(shell->zvals_sct_r, shell->zvals_w, x, INSERT_VALUES, SCATTER_REVERSE));
98: }
99: PetscFunctionReturn(PETSC_SUCCESS);
100: }
102: /*
103: Store and scale values on zeroed rows
104: xx = [x_1, 0], 0 on zeroed rows
105: */
106: static PetscErrorCode MatShellPreZeroLeft(Mat A, Vec x, Vec *xx)
107: {
108: Mat_Shell *shell = (Mat_Shell *)A->data;
110: PetscFunctionBegin;
111: *xx = NULL;
112: if (!shell->zrows) {
113: *xx = x;
114: } else {
115: if (!shell->left_work) PetscCall(MatCreateVecs(A, NULL, &shell->left_work));
116: PetscCall(VecCopy(x, shell->left_work));
117: PetscCall(VecSet(shell->zvals_w, 0.0));
118: PetscCall(VecScatterBegin(shell->zvals_sct_r, shell->zvals_w, shell->left_work, INSERT_VALUES, SCATTER_REVERSE));
119: PetscCall(VecScatterEnd(shell->zvals_sct_r, shell->zvals_w, shell->left_work, INSERT_VALUES, SCATTER_REVERSE));
120: PetscCall(VecScatterBegin(shell->zvals_sct_r, x, shell->zvals_w, INSERT_VALUES, SCATTER_FORWARD));
121: PetscCall(VecScatterEnd(shell->zvals_sct_r, x, shell->zvals_w, INSERT_VALUES, SCATTER_FORWARD));
122: PetscCall(VecPointwiseMult(shell->zvals_w, shell->zvals_w, shell->zvals));
123: *xx = shell->left_work;
124: }
125: PetscFunctionReturn(PETSC_SUCCESS);
126: }
128: /* Zero zero-columns contributions, sum contributions from properly scaled values stored in MatShellPreZeroLeft */
129: static PetscErrorCode MatShellPostZeroRight(Mat A, Vec x)
130: {
131: Mat_Shell *shell = (Mat_Shell *)A->data;
133: PetscFunctionBegin;
134: if (shell->zcols) PetscCall(VecISSet(x, shell->zcols, 0.0));
135: if (shell->zrows) {
136: PetscCall(VecScatterBegin(shell->zvals_sct_c, shell->zvals_w, x, ADD_VALUES, SCATTER_REVERSE));
137: PetscCall(VecScatterEnd(shell->zvals_sct_c, shell->zvals_w, x, ADD_VALUES, SCATTER_REVERSE));
138: }
139: PetscFunctionReturn(PETSC_SUCCESS);
140: }
142: /*
143: xx = diag(left)*x
144: */
145: static PetscErrorCode MatShellPreScaleLeft(Mat A, Vec x, Vec *xx)
146: {
147: Mat_Shell *shell = (Mat_Shell *)A->data;
149: PetscFunctionBegin;
150: *xx = NULL;
151: if (!shell->left) {
152: *xx = x;
153: } else {
154: if (!shell->left_work) PetscCall(VecDuplicate(shell->left, &shell->left_work));
155: PetscCall(VecPointwiseMult(shell->left_work, x, shell->left));
156: *xx = shell->left_work;
157: }
158: PetscFunctionReturn(PETSC_SUCCESS);
159: }
161: /*
162: xx = diag(right)*x
163: */
164: static PetscErrorCode MatShellPreScaleRight(Mat A, Vec x, Vec *xx)
165: {
166: Mat_Shell *shell = (Mat_Shell *)A->data;
168: PetscFunctionBegin;
169: *xx = NULL;
170: if (!shell->right) {
171: *xx = x;
172: } else {
173: if (!shell->right_work) PetscCall(VecDuplicate(shell->right, &shell->right_work));
174: PetscCall(VecPointwiseMult(shell->right_work, x, shell->right));
175: *xx = shell->right_work;
176: }
177: PetscFunctionReturn(PETSC_SUCCESS);
178: }
180: /*
181: x = diag(left)*x
182: */
183: static PetscErrorCode MatShellPostScaleLeft(Mat A, Vec x)
184: {
185: Mat_Shell *shell = (Mat_Shell *)A->data;
187: PetscFunctionBegin;
188: if (shell->left) PetscCall(VecPointwiseMult(x, x, shell->left));
189: PetscFunctionReturn(PETSC_SUCCESS);
190: }
192: /*
193: x = diag(right)*x
194: */
195: static PetscErrorCode MatShellPostScaleRight(Mat A, Vec x)
196: {
197: Mat_Shell *shell = (Mat_Shell *)A->data;
199: PetscFunctionBegin;
200: if (shell->right) PetscCall(VecPointwiseMult(x, x, shell->right));
201: PetscFunctionReturn(PETSC_SUCCESS);
202: }
204: /*
205: Y = vscale*Y + diag(dshift)*X + vshift*X
207: On input Y already contains A*x
208: */
209: static PetscErrorCode MatShellShiftAndScale(Mat A, Vec X, Vec Y)
210: {
211: Mat_Shell *shell = (Mat_Shell *)A->data;
213: PetscFunctionBegin;
214: if (shell->dshift) { /* get arrays because there is no VecPointwiseMultAdd() */
215: PetscInt i, m;
216: const PetscScalar *x, *d;
217: PetscScalar *y;
218: PetscCall(VecGetLocalSize(X, &m));
219: PetscCall(VecGetArrayRead(shell->dshift, &d));
220: PetscCall(VecGetArrayRead(X, &x));
221: PetscCall(VecGetArray(Y, &y));
222: for (i = 0; i < m; i++) y[i] = shell->vscale * y[i] + d[i] * x[i];
223: PetscCall(VecRestoreArrayRead(shell->dshift, &d));
224: PetscCall(VecRestoreArrayRead(X, &x));
225: PetscCall(VecRestoreArray(Y, &y));
226: } else {
227: PetscCall(VecScale(Y, shell->vscale));
228: }
229: if (shell->vshift != 0.0) PetscCall(VecAXPY(Y, shell->vshift, X)); /* if test is for non-square matrices */
230: PetscFunctionReturn(PETSC_SUCCESS);
231: }
233: static PetscErrorCode MatShellGetContext_Shell(Mat mat, void *ctx)
234: {
235: Mat_Shell *shell = (Mat_Shell *)mat->data;
237: PetscFunctionBegin;
238: if (shell->ctxcontainer) PetscCall(PetscContainerGetPointer(shell->ctxcontainer, (void **)ctx));
239: else *(void **)ctx = NULL;
240: PetscFunctionReturn(PETSC_SUCCESS);
241: }
243: /*@
244: MatShellGetContext - Returns the user-provided context associated with a `MATSHELL` shell matrix.
246: Not Collective
248: Input Parameter:
249: . mat - the matrix, should have been created with `MatCreateShell()`
251: Output Parameter:
252: . ctx - the user provided context
254: Level: advanced
256: Fortran Notes:
257: You must write a Fortran interface definition for this
258: function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.
260: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellSetOperation()`, `MatShellSetContext()`
261: @*/
262: PetscErrorCode MatShellGetContext(Mat mat, void *ctx)
263: {
264: PetscFunctionBegin;
266: PetscAssertPointer(ctx, 2);
267: PetscUseMethod(mat, "MatShellGetContext_C", (Mat, void *), (mat, ctx));
268: PetscFunctionReturn(PETSC_SUCCESS);
269: }
271: static PetscErrorCode MatZeroRowsColumns_Local_Shell(Mat mat, PetscInt nr, PetscInt rows[], PetscInt nc, PetscInt cols[], PetscScalar diag, PetscBool rc)
272: {
273: Mat_Shell *shell = (Mat_Shell *)mat->data;
274: Vec x = NULL, b = NULL;
275: IS is1, is2;
276: const PetscInt *ridxs;
277: PetscInt *idxs, *gidxs;
278: PetscInt cum, rst, cst, i;
280: PetscFunctionBegin;
281: if (!shell->zvals) PetscCall(MatCreateVecs(mat, NULL, &shell->zvals));
282: if (!shell->zvals_w) PetscCall(VecDuplicate(shell->zvals, &shell->zvals_w));
283: PetscCall(MatGetOwnershipRange(mat, &rst, NULL));
284: PetscCall(MatGetOwnershipRangeColumn(mat, &cst, NULL));
286: /* Expand/create index set of zeroed rows */
287: PetscCall(PetscMalloc1(nr, &idxs));
288: for (i = 0; i < nr; i++) idxs[i] = rows[i] + rst;
289: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nr, idxs, PETSC_OWN_POINTER, &is1));
290: PetscCall(ISSort(is1));
291: PetscCall(VecISSet(shell->zvals, is1, diag));
292: if (shell->zrows) {
293: PetscCall(ISSum(shell->zrows, is1, &is2));
294: PetscCall(ISDestroy(&shell->zrows));
295: PetscCall(ISDestroy(&is1));
296: shell->zrows = is2;
297: } else shell->zrows = is1;
299: /* Create scatters for diagonal values communications */
300: PetscCall(VecScatterDestroy(&shell->zvals_sct_c));
301: PetscCall(VecScatterDestroy(&shell->zvals_sct_r));
303: /* row scatter: from/to left vector */
304: PetscCall(MatCreateVecs(mat, &x, &b));
305: PetscCall(VecScatterCreate(b, shell->zrows, shell->zvals_w, shell->zrows, &shell->zvals_sct_r));
307: /* col scatter: from right vector to left vector */
308: PetscCall(ISGetIndices(shell->zrows, &ridxs));
309: PetscCall(ISGetLocalSize(shell->zrows, &nr));
310: PetscCall(PetscMalloc1(nr, &gidxs));
311: for (i = 0, cum = 0; i < nr; i++) {
312: if (ridxs[i] >= mat->cmap->N) continue;
313: gidxs[cum] = ridxs[i];
314: cum++;
315: }
316: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)mat), cum, gidxs, PETSC_OWN_POINTER, &is1));
317: PetscCall(VecScatterCreate(x, is1, shell->zvals_w, is1, &shell->zvals_sct_c));
318: PetscCall(ISDestroy(&is1));
319: PetscCall(VecDestroy(&x));
320: PetscCall(VecDestroy(&b));
322: /* Expand/create index set of zeroed columns */
323: if (rc) {
324: PetscCall(PetscMalloc1(nc, &idxs));
325: for (i = 0; i < nc; i++) idxs[i] = cols[i] + cst;
326: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nc, idxs, PETSC_OWN_POINTER, &is1));
327: PetscCall(ISSort(is1));
328: if (shell->zcols) {
329: PetscCall(ISSum(shell->zcols, is1, &is2));
330: PetscCall(ISDestroy(&shell->zcols));
331: PetscCall(ISDestroy(&is1));
332: shell->zcols = is2;
333: } else shell->zcols = is1;
334: }
335: PetscFunctionReturn(PETSC_SUCCESS);
336: }
338: static PetscErrorCode MatZeroRows_Shell(Mat mat, PetscInt n, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
339: {
340: Mat_Shell *shell = (Mat_Shell *)mat->data;
341: PetscInt nr, *lrows;
343: PetscFunctionBegin;
344: if (x && b) {
345: Vec xt;
346: PetscScalar *vals;
347: PetscInt *gcols, i, st, nl, nc;
349: PetscCall(PetscMalloc1(n, &gcols));
350: for (i = 0, nc = 0; i < n; i++)
351: if (rows[i] < mat->cmap->N) gcols[nc++] = rows[i];
353: PetscCall(MatCreateVecs(mat, &xt, NULL));
354: PetscCall(VecCopy(x, xt));
355: PetscCall(PetscCalloc1(nc, &vals));
356: PetscCall(VecSetValues(xt, nc, gcols, vals, INSERT_VALUES)); /* xt = [x1, 0] */
357: PetscCall(PetscFree(vals));
358: PetscCall(VecAssemblyBegin(xt));
359: PetscCall(VecAssemblyEnd(xt));
360: PetscCall(VecAYPX(xt, -1.0, x)); /* xt = [0, x2] */
362: PetscCall(VecGetOwnershipRange(xt, &st, NULL));
363: PetscCall(VecGetLocalSize(xt, &nl));
364: PetscCall(VecGetArray(xt, &vals));
365: for (i = 0; i < nl; i++) {
366: PetscInt g = i + st;
367: if (g > mat->rmap->N) continue;
368: if (PetscAbsScalar(vals[i]) == 0.0) continue;
369: PetscCall(VecSetValue(b, g, diag * vals[i], INSERT_VALUES));
370: }
371: PetscCall(VecRestoreArray(xt, &vals));
372: PetscCall(VecAssemblyBegin(b));
373: PetscCall(VecAssemblyEnd(b)); /* b = [b1, x2 * diag] */
374: PetscCall(VecDestroy(&xt));
375: PetscCall(PetscFree(gcols));
376: }
377: PetscCall(PetscLayoutMapLocal(mat->rmap, n, rows, &nr, &lrows, NULL));
378: PetscCall(MatZeroRowsColumns_Local_Shell(mat, nr, lrows, 0, NULL, diag, PETSC_FALSE));
379: if (shell->axpy) PetscCall(MatZeroRows(shell->axpy, n, rows, 0.0, NULL, NULL));
380: PetscCall(PetscFree(lrows));
381: PetscFunctionReturn(PETSC_SUCCESS);
382: }
384: static PetscErrorCode MatZeroRowsColumns_Shell(Mat mat, PetscInt n, const PetscInt rowscols[], PetscScalar diag, Vec x, Vec b)
385: {
386: Mat_Shell *shell = (Mat_Shell *)mat->data;
387: PetscInt *lrows, *lcols;
388: PetscInt nr, nc;
389: PetscBool congruent;
391: PetscFunctionBegin;
392: if (x && b) {
393: Vec xt, bt;
394: PetscScalar *vals;
395: PetscInt *grows, *gcols, i, st, nl;
397: PetscCall(PetscMalloc2(n, &grows, n, &gcols));
398: for (i = 0, nr = 0; i < n; i++)
399: if (rowscols[i] < mat->rmap->N) grows[nr++] = rowscols[i];
400: for (i = 0, nc = 0; i < n; i++)
401: if (rowscols[i] < mat->cmap->N) gcols[nc++] = rowscols[i];
402: PetscCall(PetscCalloc1(n, &vals));
404: PetscCall(MatCreateVecs(mat, &xt, &bt));
405: PetscCall(VecCopy(x, xt));
406: PetscCall(VecSetValues(xt, nc, gcols, vals, INSERT_VALUES)); /* xt = [x1, 0] */
407: PetscCall(VecAssemblyBegin(xt));
408: PetscCall(VecAssemblyEnd(xt));
409: PetscCall(VecAXPY(xt, -1.0, x)); /* xt = [0, -x2] */
410: PetscCall(MatMult(mat, xt, bt)); /* bt = [-A12*x2,-A22*x2] */
411: PetscCall(VecSetValues(bt, nr, grows, vals, INSERT_VALUES)); /* bt = [-A12*x2,0] */
412: PetscCall(VecAssemblyBegin(bt));
413: PetscCall(VecAssemblyEnd(bt));
414: PetscCall(VecAXPY(b, 1.0, bt)); /* b = [b1 - A12*x2, b2] */
415: PetscCall(VecSetValues(bt, nr, grows, vals, INSERT_VALUES)); /* b = [b1 - A12*x2, 0] */
416: PetscCall(VecAssemblyBegin(bt));
417: PetscCall(VecAssemblyEnd(bt));
418: PetscCall(PetscFree(vals));
420: PetscCall(VecGetOwnershipRange(xt, &st, NULL));
421: PetscCall(VecGetLocalSize(xt, &nl));
422: PetscCall(VecGetArray(xt, &vals));
423: for (i = 0; i < nl; i++) {
424: PetscInt g = i + st;
425: if (g > mat->rmap->N) continue;
426: if (PetscAbsScalar(vals[i]) == 0.0) continue;
427: PetscCall(VecSetValue(b, g, -diag * vals[i], INSERT_VALUES));
428: }
429: PetscCall(VecRestoreArray(xt, &vals));
430: PetscCall(VecAssemblyBegin(b));
431: PetscCall(VecAssemblyEnd(b)); /* b = [b1 - A12*x2, x2 * diag] */
432: PetscCall(VecDestroy(&xt));
433: PetscCall(VecDestroy(&bt));
434: PetscCall(PetscFree2(grows, gcols));
435: }
436: PetscCall(PetscLayoutMapLocal(mat->rmap, n, rowscols, &nr, &lrows, NULL));
437: PetscCall(MatHasCongruentLayouts(mat, &congruent));
438: if (congruent) {
439: nc = nr;
440: lcols = lrows;
441: } else { /* MatZeroRowsColumns implicitly assumes the rowscols indices are for a square matrix, here we handle a more general case */
442: PetscInt i, nt, *t;
444: PetscCall(PetscMalloc1(n, &t));
445: for (i = 0, nt = 0; i < n; i++)
446: if (rowscols[i] < mat->cmap->N) t[nt++] = rowscols[i];
447: PetscCall(PetscLayoutMapLocal(mat->cmap, nt, t, &nc, &lcols, NULL));
448: PetscCall(PetscFree(t));
449: }
450: PetscCall(MatZeroRowsColumns_Local_Shell(mat, nr, lrows, nc, lcols, diag, PETSC_TRUE));
451: if (!congruent) PetscCall(PetscFree(lcols));
452: PetscCall(PetscFree(lrows));
453: if (shell->axpy) PetscCall(MatZeroRowsColumns(shell->axpy, n, rowscols, 0.0, NULL, NULL));
454: PetscFunctionReturn(PETSC_SUCCESS);
455: }
457: static PetscErrorCode MatDestroy_Shell(Mat mat)
458: {
459: Mat_Shell *shell = (Mat_Shell *)mat->data;
460: MatShellMatFunctionList matmat;
462: PetscFunctionBegin;
463: if (shell->ops->destroy) PetscCall((*shell->ops->destroy)(mat));
464: PetscCall(PetscMemzero(shell->ops, sizeof(struct _MatShellOps)));
465: PetscCall(VecDestroy(&shell->left));
466: PetscCall(VecDestroy(&shell->right));
467: PetscCall(VecDestroy(&shell->dshift));
468: PetscCall(VecDestroy(&shell->left_work));
469: PetscCall(VecDestroy(&shell->right_work));
470: PetscCall(VecDestroy(&shell->left_add_work));
471: PetscCall(VecDestroy(&shell->right_add_work));
472: PetscCall(VecDestroy(&shell->axpy_left));
473: PetscCall(VecDestroy(&shell->axpy_right));
474: PetscCall(MatDestroy(&shell->axpy));
475: PetscCall(VecDestroy(&shell->zvals_w));
476: PetscCall(VecDestroy(&shell->zvals));
477: PetscCall(VecScatterDestroy(&shell->zvals_sct_c));
478: PetscCall(VecScatterDestroy(&shell->zvals_sct_r));
479: PetscCall(ISDestroy(&shell->zrows));
480: PetscCall(ISDestroy(&shell->zcols));
482: matmat = shell->matmat;
483: while (matmat) {
484: MatShellMatFunctionList next = matmat->next;
486: PetscCall(PetscObjectComposeFunction((PetscObject)mat, matmat->composedname, NULL));
487: PetscCall(PetscFree(matmat->composedname));
488: PetscCall(PetscFree(matmat->resultname));
489: PetscCall(PetscFree(matmat));
490: matmat = next;
491: }
492: PetscCall(MatShellSetContext(mat, NULL));
493: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatShellGetContext_C", NULL));
494: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatShellSetContext_C", NULL));
495: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatShellSetContextDestroy_C", NULL));
496: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatShellSetVecType_C", NULL));
497: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatShellSetManageScalingShifts_C", NULL));
498: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatShellSetOperation_C", NULL));
499: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatShellGetOperation_C", NULL));
500: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatShellSetMatProductOperation_C", NULL));
501: PetscCall(PetscFree(mat->data));
502: PetscFunctionReturn(PETSC_SUCCESS);
503: }
505: typedef struct {
506: PetscErrorCode (*numeric)(Mat, Mat, Mat, void *);
507: PetscErrorCode (*destroy)(void *);
508: void *userdata;
509: Mat B;
510: Mat Bt;
511: Mat axpy;
512: } MatMatDataShell;
514: static PetscErrorCode DestroyMatMatDataShell(void *data)
515: {
516: MatMatDataShell *mmdata = (MatMatDataShell *)data;
518: PetscFunctionBegin;
519: if (mmdata->destroy) PetscCall((*mmdata->destroy)(mmdata->userdata));
520: PetscCall(MatDestroy(&mmdata->B));
521: PetscCall(MatDestroy(&mmdata->Bt));
522: PetscCall(MatDestroy(&mmdata->axpy));
523: PetscCall(PetscFree(mmdata));
524: PetscFunctionReturn(PETSC_SUCCESS);
525: }
527: static PetscErrorCode MatProductNumeric_Shell_X(Mat D)
528: {
529: Mat_Product *product;
530: Mat A, B;
531: MatMatDataShell *mdata;
532: PetscScalar zero = 0.0;
534: PetscFunctionBegin;
535: MatCheckProduct(D, 1);
536: product = D->product;
537: PetscCheck(product->data, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "Product data empty");
538: A = product->A;
539: B = product->B;
540: mdata = (MatMatDataShell *)product->data;
541: if (mdata->numeric) {
542: Mat_Shell *shell = (Mat_Shell *)A->data;
543: PetscErrorCode (*stashsym)(Mat) = D->ops->productsymbolic;
544: PetscErrorCode (*stashnum)(Mat) = D->ops->productnumeric;
545: PetscBool useBmdata = PETSC_FALSE, newB = PETSC_TRUE;
547: if (shell->managescalingshifts) {
548: PetscCheck(!shell->zcols && !shell->zrows, PetscObjectComm((PetscObject)D), PETSC_ERR_SUP, "MatProduct not supported with zeroed rows/columns");
549: if (shell->right || shell->left) {
550: useBmdata = PETSC_TRUE;
551: if (!mdata->B) {
552: PetscCall(MatDuplicate(B, MAT_SHARE_NONZERO_PATTERN, &mdata->B));
553: } else {
554: newB = PETSC_FALSE;
555: }
556: PetscCall(MatCopy(B, mdata->B, SAME_NONZERO_PATTERN));
557: }
558: switch (product->type) {
559: case MATPRODUCT_AB: /* s L A R B + v L R B + L D R B */
560: if (shell->right) PetscCall(MatDiagonalScale(mdata->B, shell->right, NULL));
561: break;
562: case MATPRODUCT_AtB: /* s R A^t L B + v R L B + R D L B */
563: if (shell->left) PetscCall(MatDiagonalScale(mdata->B, shell->left, NULL));
564: break;
565: case MATPRODUCT_ABt: /* s L A R B^t + v L R B^t + L D R B^t */
566: if (shell->right) PetscCall(MatDiagonalScale(mdata->B, NULL, shell->right));
567: break;
568: case MATPRODUCT_RARt: /* s B L A R B^t + v B L R B^t + B L D R B^t */
569: if (shell->right && shell->left) {
570: PetscBool flg;
572: PetscCall(VecEqual(shell->right, shell->left, &flg));
573: PetscCheck(flg, PetscObjectComm((PetscObject)D), PETSC_ERR_SUP, "MatProductSymbolic type %s not supported for %s and %s matrices because left scaling != from right scaling", MatProductTypes[product->type], ((PetscObject)A)->type_name,
574: ((PetscObject)B)->type_name);
575: }
576: if (shell->right) PetscCall(MatDiagonalScale(mdata->B, NULL, shell->right));
577: break;
578: case MATPRODUCT_PtAP: /* s B^t L A R B + v B^t L R B + B^t L D R B */
579: if (shell->right && shell->left) {
580: PetscBool flg;
582: PetscCall(VecEqual(shell->right, shell->left, &flg));
583: PetscCheck(flg, PetscObjectComm((PetscObject)D), PETSC_ERR_SUP, "MatProductSymbolic type %s not supported for %s and %s matrices because left scaling != from right scaling", MatProductTypes[product->type], ((PetscObject)A)->type_name,
584: ((PetscObject)B)->type_name);
585: }
586: if (shell->right) PetscCall(MatDiagonalScale(mdata->B, shell->right, NULL));
587: break;
588: default:
589: SETERRQ(PetscObjectComm((PetscObject)D), PETSC_ERR_SUP, "MatProductSymbolic type %s not supported for %s and %s matrices", MatProductTypes[product->type], ((PetscObject)A)->type_name, ((PetscObject)B)->type_name);
590: }
591: }
592: /* allow the user to call MatMat operations on D */
593: D->product = NULL;
594: D->ops->productsymbolic = NULL;
595: D->ops->productnumeric = NULL;
597: PetscCall((*mdata->numeric)(A, useBmdata ? mdata->B : B, D, mdata->userdata));
599: /* clear any leftover user data and restore D pointers */
600: PetscCall(MatProductClear(D));
601: D->ops->productsymbolic = stashsym;
602: D->ops->productnumeric = stashnum;
603: D->product = product;
605: if (shell->managescalingshifts) {
606: PetscCall(MatScale(D, shell->vscale));
607: switch (product->type) {
608: case MATPRODUCT_AB: /* s L A R B + v L R B + L D R B */
609: case MATPRODUCT_ABt: /* s L A R B^t + v L R B^t + L D R B^t */
610: if (shell->left) {
611: PetscCall(MatDiagonalScale(D, shell->left, NULL));
612: if (shell->dshift || shell->vshift != zero) {
613: if (!shell->left_work) PetscCall(MatCreateVecs(A, NULL, &shell->left_work));
614: if (shell->dshift) {
615: PetscCall(VecCopy(shell->dshift, shell->left_work));
616: PetscCall(VecShift(shell->left_work, shell->vshift));
617: PetscCall(VecPointwiseMult(shell->left_work, shell->left_work, shell->left));
618: } else {
619: PetscCall(VecSet(shell->left_work, shell->vshift));
620: }
621: if (product->type == MATPRODUCT_ABt) {
622: MatReuse reuse = mdata->Bt ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX;
623: MatStructure str = mdata->Bt ? SUBSET_NONZERO_PATTERN : DIFFERENT_NONZERO_PATTERN;
625: PetscCall(MatTranspose(mdata->B, reuse, &mdata->Bt));
626: PetscCall(MatDiagonalScale(mdata->Bt, shell->left_work, NULL));
627: PetscCall(MatAXPY(D, 1.0, mdata->Bt, str));
628: } else {
629: MatStructure str = newB ? DIFFERENT_NONZERO_PATTERN : SUBSET_NONZERO_PATTERN;
631: PetscCall(MatDiagonalScale(mdata->B, shell->left_work, NULL));
632: PetscCall(MatAXPY(D, 1.0, mdata->B, str));
633: }
634: }
635: }
636: break;
637: case MATPRODUCT_AtB: /* s R A^t L B + v R L B + R D L B */
638: if (shell->right) {
639: PetscCall(MatDiagonalScale(D, shell->right, NULL));
640: if (shell->dshift || shell->vshift != zero) {
641: MatStructure str = newB ? DIFFERENT_NONZERO_PATTERN : SUBSET_NONZERO_PATTERN;
643: if (!shell->right_work) PetscCall(MatCreateVecs(A, &shell->right_work, NULL));
644: if (shell->dshift) {
645: PetscCall(VecCopy(shell->dshift, shell->right_work));
646: PetscCall(VecShift(shell->right_work, shell->vshift));
647: PetscCall(VecPointwiseMult(shell->right_work, shell->right_work, shell->right));
648: } else {
649: PetscCall(VecSet(shell->right_work, shell->vshift));
650: }
651: PetscCall(MatDiagonalScale(mdata->B, shell->right_work, NULL));
652: PetscCall(MatAXPY(D, 1.0, mdata->B, str));
653: }
654: }
655: break;
656: case MATPRODUCT_PtAP: /* s B^t L A R B + v B^t L R B + B^t L D R B */
657: case MATPRODUCT_RARt: /* s B L A R B^t + v B L R B^t + B L D R B^t */
658: PetscCheck(!shell->dshift && shell->vshift == zero, PetscObjectComm((PetscObject)D), PETSC_ERR_SUP, "MatProductSymbolic type %s not supported for %s and %s matrices with diagonal shift", MatProductTypes[product->type], ((PetscObject)A)->type_name,
659: ((PetscObject)B)->type_name);
660: break;
661: default:
662: SETERRQ(PetscObjectComm((PetscObject)D), PETSC_ERR_SUP, "MatProductSymbolic type %s not supported for %s and %s matrices", MatProductTypes[product->type], ((PetscObject)A)->type_name, ((PetscObject)B)->type_name);
663: }
664: if (shell->axpy && shell->axpy_vscale != zero) {
665: Mat X;
666: PetscObjectState axpy_state;
667: MatStructure str = DIFFERENT_NONZERO_PATTERN; /* not sure it is safe to ever use SUBSET_NONZERO_PATTERN */
669: PetscCall(MatShellGetContext(shell->axpy, &X));
670: PetscCall(PetscObjectStateGet((PetscObject)X, &axpy_state));
671: PetscCheck(shell->axpy_state == axpy_state, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Invalid AXPY state: cannot modify the X matrix passed to MatAXPY(Y,a,X,...)");
672: if (!mdata->axpy) {
673: str = DIFFERENT_NONZERO_PATTERN;
674: PetscCall(MatProductCreate(shell->axpy, B, NULL, &mdata->axpy));
675: PetscCall(MatProductSetType(mdata->axpy, product->type));
676: PetscCall(MatProductSetFromOptions(mdata->axpy));
677: PetscCall(MatProductSymbolic(mdata->axpy));
678: } else { /* May be that shell->axpy has changed */
679: PetscBool flg;
681: PetscCall(MatProductReplaceMats(shell->axpy, B, NULL, mdata->axpy));
682: PetscCall(MatHasOperation(mdata->axpy, MATOP_PRODUCTSYMBOLIC, &flg));
683: if (!flg) {
684: str = DIFFERENT_NONZERO_PATTERN;
685: PetscCall(MatProductSetFromOptions(mdata->axpy));
686: PetscCall(MatProductSymbolic(mdata->axpy));
687: }
688: }
689: PetscCall(MatProductNumeric(mdata->axpy));
690: PetscCall(MatAXPY(D, shell->axpy_vscale, mdata->axpy, str));
691: }
692: }
693: } else SETERRQ(PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "Missing numeric operation");
694: PetscFunctionReturn(PETSC_SUCCESS);
695: }
697: static PetscErrorCode MatProductSymbolic_Shell_X(Mat D)
698: {
699: Mat_Product *product;
700: Mat A, B;
701: MatShellMatFunctionList matmat;
702: Mat_Shell *shell;
703: PetscBool flg = PETSC_FALSE;
704: char composedname[256];
705: MatMatDataShell *mdata;
707: PetscFunctionBegin;
708: MatCheckProduct(D, 1);
709: product = D->product;
710: PetscCheck(!product->data, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "Product data not empty");
711: A = product->A;
712: B = product->B;
713: shell = (Mat_Shell *)A->data;
714: matmat = shell->matmat;
715: PetscCall(PetscSNPrintf(composedname, sizeof(composedname), "MatProductSetFromOptions_%s_%s_C", ((PetscObject)A)->type_name, ((PetscObject)B)->type_name));
716: while (matmat) {
717: PetscCall(PetscStrcmp(composedname, matmat->composedname, &flg));
718: flg = (PetscBool)(flg && (matmat->ptype == product->type));
719: if (flg) break;
720: matmat = matmat->next;
721: }
722: PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Composedname \"%s\" for product type %s not found", composedname, MatProductTypes[product->type]);
723: switch (product->type) {
724: case MATPRODUCT_AB:
725: PetscCall(MatSetSizes(D, A->rmap->n, B->cmap->n, A->rmap->N, B->cmap->N));
726: break;
727: case MATPRODUCT_AtB:
728: PetscCall(MatSetSizes(D, A->cmap->n, B->cmap->n, A->cmap->N, B->cmap->N));
729: break;
730: case MATPRODUCT_ABt:
731: PetscCall(MatSetSizes(D, A->rmap->n, B->rmap->n, A->rmap->N, B->rmap->N));
732: break;
733: case MATPRODUCT_RARt:
734: PetscCall(MatSetSizes(D, B->rmap->n, B->rmap->n, B->rmap->N, B->rmap->N));
735: break;
736: case MATPRODUCT_PtAP:
737: PetscCall(MatSetSizes(D, B->cmap->n, B->cmap->n, B->cmap->N, B->cmap->N));
738: break;
739: default:
740: SETERRQ(PetscObjectComm((PetscObject)D), PETSC_ERR_SUP, "MatProductSymbolic type %s not supported for %s and %s matrices", MatProductTypes[product->type], ((PetscObject)A)->type_name, ((PetscObject)B)->type_name);
741: }
742: /* respect users who passed in a matrix for which resultname is the base type */
743: if (matmat->resultname) {
744: PetscCall(PetscObjectBaseTypeCompare((PetscObject)D, matmat->resultname, &flg));
745: if (!flg) PetscCall(MatSetType(D, matmat->resultname));
746: }
747: /* If matrix type was not set or different, we need to reset this pointers */
748: D->ops->productsymbolic = MatProductSymbolic_Shell_X;
749: D->ops->productnumeric = MatProductNumeric_Shell_X;
750: /* attach product data */
751: PetscCall(PetscNew(&mdata));
752: mdata->numeric = matmat->numeric;
753: mdata->destroy = matmat->destroy;
754: if (matmat->symbolic) {
755: PetscCall((*matmat->symbolic)(A, B, D, &mdata->userdata));
756: } else { /* call general setup if symbolic operation not provided */
757: PetscCall(MatSetUp(D));
758: }
759: PetscCheck(D->product, PetscObjectComm((PetscObject)D), PETSC_ERR_COR, "Product disappeared after user symbolic phase");
760: PetscCheck(!D->product->data, PetscObjectComm((PetscObject)D), PETSC_ERR_COR, "Product data not empty after user symbolic phase");
761: D->product->data = mdata;
762: D->product->destroy = DestroyMatMatDataShell;
763: /* Be sure to reset these pointers if the user did something unexpected */
764: D->ops->productsymbolic = MatProductSymbolic_Shell_X;
765: D->ops->productnumeric = MatProductNumeric_Shell_X;
766: PetscFunctionReturn(PETSC_SUCCESS);
767: }
769: static PetscErrorCode MatProductSetFromOptions_Shell_X(Mat D)
770: {
771: Mat_Product *product;
772: Mat A, B;
773: MatShellMatFunctionList matmat;
774: Mat_Shell *shell;
775: PetscBool flg;
776: char composedname[256];
778: PetscFunctionBegin;
779: MatCheckProduct(D, 1);
780: product = D->product;
781: A = product->A;
782: B = product->B;
783: PetscCall(MatIsShell(A, &flg));
784: if (!flg) PetscFunctionReturn(PETSC_SUCCESS);
785: shell = (Mat_Shell *)A->data;
786: matmat = shell->matmat;
787: PetscCall(PetscSNPrintf(composedname, sizeof(composedname), "MatProductSetFromOptions_%s_%s_C", ((PetscObject)A)->type_name, ((PetscObject)B)->type_name));
788: while (matmat) {
789: PetscCall(PetscStrcmp(composedname, matmat->composedname, &flg));
790: flg = (PetscBool)(flg && (matmat->ptype == product->type));
791: if (flg) break;
792: matmat = matmat->next;
793: }
794: if (flg) {
795: D->ops->productsymbolic = MatProductSymbolic_Shell_X;
796: } else PetscCall(PetscInfo(D, " symbolic product %s not registered for product type %s\n", composedname, MatProductTypes[product->type]));
797: PetscFunctionReturn(PETSC_SUCCESS);
798: }
800: static PetscErrorCode MatShellSetMatProductOperation_Private(Mat A, MatProductType ptype, PetscErrorCode (*symbolic)(Mat, Mat, Mat, void **), PetscErrorCode (*numeric)(Mat, Mat, Mat, void *), PetscErrorCode (*destroy)(void *), char *composedname, const char *resultname)
801: {
802: PetscBool flg;
803: Mat_Shell *shell;
804: MatShellMatFunctionList matmat;
806: PetscFunctionBegin;
807: PetscCheck(numeric, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing numeric routine");
808: PetscCheck(composedname, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing composed name");
810: /* add product callback */
811: shell = (Mat_Shell *)A->data;
812: matmat = shell->matmat;
813: if (!matmat) {
814: PetscCall(PetscNew(&shell->matmat));
815: matmat = shell->matmat;
816: } else {
817: MatShellMatFunctionList entry = matmat;
818: while (entry) {
819: PetscCall(PetscStrcmp(composedname, entry->composedname, &flg));
820: flg = (PetscBool)(flg && (entry->ptype == ptype));
821: matmat = entry;
822: if (flg) goto set;
823: entry = entry->next;
824: }
825: PetscCall(PetscNew(&matmat->next));
826: matmat = matmat->next;
827: }
829: set:
830: matmat->symbolic = symbolic;
831: matmat->numeric = numeric;
832: matmat->destroy = destroy;
833: matmat->ptype = ptype;
834: PetscCall(PetscFree(matmat->composedname));
835: PetscCall(PetscFree(matmat->resultname));
836: PetscCall(PetscStrallocpy(composedname, &matmat->composedname));
837: PetscCall(PetscStrallocpy(resultname, &matmat->resultname));
838: PetscCall(PetscInfo(A, "Composing %s for product type %s with result %s\n", matmat->composedname, MatProductTypes[matmat->ptype], matmat->resultname ? matmat->resultname : "not specified"));
839: PetscCall(PetscObjectComposeFunction((PetscObject)A, matmat->composedname, MatProductSetFromOptions_Shell_X));
840: PetscFunctionReturn(PETSC_SUCCESS);
841: }
843: /*@C
844: MatShellSetMatProductOperation - Allows user to set a matrix matrix operation for a `MATSHELL` shell matrix.
846: Logically Collective; No Fortran Support
848: Input Parameters:
849: + A - the `MATSHELL` shell matrix
850: . ptype - the product type
851: . symbolic - the function for the symbolic phase (can be `NULL`)
852: . numeric - the function for the numerical phase
853: . destroy - the function for the destruction of the needed data generated during the symbolic phase (can be `NULL`)
854: . Btype - the matrix type for the matrix to be multiplied against
855: - Ctype - the matrix type for the result (can be `NULL`)
857: Level: advanced
859: Example Usage:
860: .vb
861: extern PetscErrorCode usersymbolic(Mat, Mat, Mat, void**);
862: extern PetscErrorCode usernumeric(Mat, Mat, Mat, void*);
863: extern PetscErrorCode userdestroy(void*);
865: MatCreateShell(comm, m, n, M, N, ctx, &A);
866: MatShellSetMatProductOperation(
867: A, MATPRODUCT_AB, usersymbolic, usernumeric, userdestroy,MATSEQAIJ, MATDENSE
868: );
869: // create B of type SEQAIJ etc..
870: MatProductCreate(A, B, PETSC_NULLPTR, &C);
871: MatProductSetType(C, MATPRODUCT_AB);
872: MatProductSetFromOptions(C);
873: MatProductSymbolic(C); // actually runs the user defined symbolic operation
874: MatProductNumeric(C); // actually runs the user defined numeric operation
875: // use C = A * B
876: .ve
878: Notes:
879: `MATPRODUCT_ABC` is not supported yet.
881: If the symbolic phase is not specified, `MatSetUp()` is called on the result matrix that must have its type set if Ctype is `NULL`.
883: Any additional data needed by the matrix product needs to be returned during the symbolic phase and destroyed with the destroy callback.
884: PETSc will take care of calling the user-defined callbacks.
885: It is allowed to specify the same callbacks for different Btype matrix types.
886: The couple (Btype,ptype) uniquely identifies the operation, the last specified callbacks takes precedence.
888: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellSetContext()`, `MatSetOperation()`, `MatProductType`, `MatType`, `MatSetUp()`
889: @*/
890: PetscErrorCode MatShellSetMatProductOperation(Mat A, MatProductType ptype, PetscErrorCode (*symbolic)(Mat, Mat, Mat, void **), PetscErrorCode (*numeric)(Mat, Mat, Mat, void *), PetscErrorCode (*destroy)(void *), MatType Btype, MatType Ctype)
891: {
892: PetscFunctionBegin;
895: PetscCheck(ptype != MATPRODUCT_ABC, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not for product type %s", MatProductTypes[ptype]);
896: PetscCheck(numeric, PetscObjectComm((PetscObject)A), PETSC_ERR_USER, "Missing numeric routine, argument 4");
897: PetscAssertPointer(Btype, 6);
898: if (Ctype) PetscAssertPointer(Ctype, 7);
899: PetscTryMethod(A, "MatShellSetMatProductOperation_C", (Mat, MatProductType, PetscErrorCode(*)(Mat, Mat, Mat, void **), PetscErrorCode(*)(Mat, Mat, Mat, void *), PetscErrorCode(*)(void *), MatType, MatType), (A, ptype, symbolic, numeric, destroy, Btype, Ctype));
900: PetscFunctionReturn(PETSC_SUCCESS);
901: }
903: static PetscErrorCode MatShellSetMatProductOperation_Shell(Mat A, MatProductType ptype, PetscErrorCode (*symbolic)(Mat, Mat, Mat, void **), PetscErrorCode (*numeric)(Mat, Mat, Mat, void *), PetscErrorCode (*destroy)(void *), MatType Btype, MatType Ctype)
904: {
905: PetscBool flg;
906: char composedname[256];
907: MatRootName Bnames = MatRootNameList, Cnames = MatRootNameList;
908: PetscMPIInt size;
910: PetscFunctionBegin;
912: while (Bnames) { /* user passed in the root name */
913: PetscCall(PetscStrcmp(Btype, Bnames->rname, &flg));
914: if (flg) break;
915: Bnames = Bnames->next;
916: }
917: while (Cnames) { /* user passed in the root name */
918: PetscCall(PetscStrcmp(Ctype, Cnames->rname, &flg));
919: if (flg) break;
920: Cnames = Cnames->next;
921: }
922: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
923: Btype = Bnames ? (size > 1 ? Bnames->mname : Bnames->sname) : Btype;
924: Ctype = Cnames ? (size > 1 ? Cnames->mname : Cnames->sname) : Ctype;
925: PetscCall(PetscSNPrintf(composedname, sizeof(composedname), "MatProductSetFromOptions_%s_%s_C", ((PetscObject)A)->type_name, Btype));
926: PetscCall(MatShellSetMatProductOperation_Private(A, ptype, symbolic, numeric, destroy, composedname, Ctype));
927: PetscFunctionReturn(PETSC_SUCCESS);
928: }
930: static PetscErrorCode MatCopy_Shell(Mat A, Mat B, MatStructure str)
931: {
932: Mat_Shell *shellA = (Mat_Shell *)A->data, *shellB = (Mat_Shell *)B->data;
933: PetscBool matflg;
934: MatShellMatFunctionList matmatA;
936: PetscFunctionBegin;
937: PetscCall(MatIsShell(B, &matflg));
938: PetscCheck(matflg, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Matrix %s not derived from MATSHELL", ((PetscObject)B)->type_name);
940: B->ops[0] = A->ops[0];
941: shellB->ops[0] = shellA->ops[0];
943: if (shellA->ops->copy) PetscCall((*shellA->ops->copy)(A, B, str));
944: shellB->vscale = shellA->vscale;
945: shellB->vshift = shellA->vshift;
946: if (shellA->dshift) {
947: if (!shellB->dshift) PetscCall(VecDuplicate(shellA->dshift, &shellB->dshift));
948: PetscCall(VecCopy(shellA->dshift, shellB->dshift));
949: } else {
950: PetscCall(VecDestroy(&shellB->dshift));
951: }
952: if (shellA->left) {
953: if (!shellB->left) PetscCall(VecDuplicate(shellA->left, &shellB->left));
954: PetscCall(VecCopy(shellA->left, shellB->left));
955: } else {
956: PetscCall(VecDestroy(&shellB->left));
957: }
958: if (shellA->right) {
959: if (!shellB->right) PetscCall(VecDuplicate(shellA->right, &shellB->right));
960: PetscCall(VecCopy(shellA->right, shellB->right));
961: } else {
962: PetscCall(VecDestroy(&shellB->right));
963: }
964: PetscCall(MatDestroy(&shellB->axpy));
965: shellB->axpy_vscale = 0.0;
966: shellB->axpy_state = 0;
967: if (shellA->axpy) {
968: PetscCall(PetscObjectReference((PetscObject)shellA->axpy));
969: shellB->axpy = shellA->axpy;
970: shellB->axpy_vscale = shellA->axpy_vscale;
971: shellB->axpy_state = shellA->axpy_state;
972: }
973: if (shellA->zrows) {
974: PetscCall(ISDuplicate(shellA->zrows, &shellB->zrows));
975: if (shellA->zcols) PetscCall(ISDuplicate(shellA->zcols, &shellB->zcols));
976: PetscCall(VecDuplicate(shellA->zvals, &shellB->zvals));
977: PetscCall(VecCopy(shellA->zvals, shellB->zvals));
978: PetscCall(VecDuplicate(shellA->zvals_w, &shellB->zvals_w));
979: PetscCall(PetscObjectReference((PetscObject)shellA->zvals_sct_r));
980: PetscCall(PetscObjectReference((PetscObject)shellA->zvals_sct_c));
981: shellB->zvals_sct_r = shellA->zvals_sct_r;
982: shellB->zvals_sct_c = shellA->zvals_sct_c;
983: }
985: matmatA = shellA->matmat;
986: if (matmatA) {
987: while (matmatA->next) {
988: PetscCall(MatShellSetMatProductOperation_Private(B, matmatA->ptype, matmatA->symbolic, matmatA->numeric, matmatA->destroy, matmatA->composedname, matmatA->resultname));
989: matmatA = matmatA->next;
990: }
991: }
992: PetscFunctionReturn(PETSC_SUCCESS);
993: }
995: static PetscErrorCode MatDuplicate_Shell(Mat mat, MatDuplicateOption op, Mat *M)
996: {
997: PetscFunctionBegin;
998: PetscCall(MatCreateShell(PetscObjectComm((PetscObject)mat), mat->rmap->n, mat->cmap->n, mat->rmap->N, mat->cmap->N, NULL, M));
999: ((Mat_Shell *)(*M)->data)->ctxcontainer = ((Mat_Shell *)mat->data)->ctxcontainer;
1000: PetscCall(PetscObjectCompose((PetscObject)(*M), "MatShell ctx", (PetscObject)((Mat_Shell *)(*M)->data)->ctxcontainer));
1001: PetscCall(PetscObjectChangeTypeName((PetscObject)(*M), ((PetscObject)mat)->type_name));
1002: if (op != MAT_DO_NOT_COPY_VALUES) PetscCall(MatCopy(mat, *M, SAME_NONZERO_PATTERN));
1003: PetscCall(PetscObjectCopyFortranFunctionPointers((PetscObject)mat, (PetscObject)*M));
1004: PetscFunctionReturn(PETSC_SUCCESS);
1005: }
1007: static PetscErrorCode MatMult_Shell(Mat A, Vec x, Vec y)
1008: {
1009: Mat_Shell *shell = (Mat_Shell *)A->data;
1010: Vec xx;
1011: PetscObjectState instate, outstate;
1013: PetscFunctionBegin;
1014: PetscCall(MatShellPreZeroRight(A, x, &xx));
1015: PetscCall(MatShellPreScaleRight(A, xx, &xx));
1016: PetscCall(PetscObjectStateGet((PetscObject)y, &instate));
1017: PetscCall((*shell->ops->mult)(A, xx, y));
1018: PetscCall(PetscObjectStateGet((PetscObject)y, &outstate));
1019: if (instate == outstate) {
1020: /* increase the state of the output vector since the user did not update its state themself as should have been done */
1021: PetscCall(PetscObjectStateIncrease((PetscObject)y));
1022: }
1023: PetscCall(MatShellShiftAndScale(A, xx, y));
1024: PetscCall(MatShellPostScaleLeft(A, y));
1025: PetscCall(MatShellPostZeroLeft(A, y));
1027: if (shell->axpy) {
1028: Mat X;
1029: PetscObjectState axpy_state;
1031: PetscCall(MatShellGetContext(shell->axpy, &X));
1032: PetscCall(PetscObjectStateGet((PetscObject)X, &axpy_state));
1033: 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,...)");
1035: PetscCall(MatCreateVecs(shell->axpy, shell->axpy_right ? NULL : &shell->axpy_right, shell->axpy_left ? NULL : &shell->axpy_left));
1036: PetscCall(VecCopy(x, shell->axpy_right));
1037: PetscCall(MatMult(shell->axpy, shell->axpy_right, shell->axpy_left));
1038: PetscCall(VecAXPY(y, shell->axpy_vscale, shell->axpy_left));
1039: }
1040: PetscFunctionReturn(PETSC_SUCCESS);
1041: }
1043: static PetscErrorCode MatMultAdd_Shell(Mat A, Vec x, Vec y, Vec z)
1044: {
1045: Mat_Shell *shell = (Mat_Shell *)A->data;
1047: PetscFunctionBegin;
1048: if (y == z) {
1049: if (!shell->right_add_work) PetscCall(VecDuplicate(z, &shell->right_add_work));
1050: PetscCall(MatMult(A, x, shell->right_add_work));
1051: PetscCall(VecAXPY(z, 1.0, shell->right_add_work));
1052: } else {
1053: PetscCall(MatMult(A, x, z));
1054: PetscCall(VecAXPY(z, 1.0, y));
1055: }
1056: PetscFunctionReturn(PETSC_SUCCESS);
1057: }
1059: static PetscErrorCode MatMultTranspose_Shell(Mat A, Vec x, Vec y)
1060: {
1061: Mat_Shell *shell = (Mat_Shell *)A->data;
1062: Vec xx;
1063: PetscObjectState instate, outstate;
1065: PetscFunctionBegin;
1066: PetscCall(MatShellPreZeroLeft(A, x, &xx));
1067: PetscCall(MatShellPreScaleLeft(A, xx, &xx));
1068: PetscCall(PetscObjectStateGet((PetscObject)y, &instate));
1069: PetscCall((*shell->ops->multtranspose)(A, xx, y));
1070: PetscCall(PetscObjectStateGet((PetscObject)y, &outstate));
1071: if (instate == outstate) {
1072: /* increase the state of the output vector since the user did not update its state themself as should have been done */
1073: PetscCall(PetscObjectStateIncrease((PetscObject)y));
1074: }
1075: PetscCall(MatShellShiftAndScale(A, xx, y));
1076: PetscCall(MatShellPostScaleRight(A, y));
1077: PetscCall(MatShellPostZeroRight(A, y));
1079: if (shell->axpy) {
1080: Mat X;
1081: PetscObjectState axpy_state;
1083: PetscCall(MatShellGetContext(shell->axpy, &X));
1084: PetscCall(PetscObjectStateGet((PetscObject)X, &axpy_state));
1085: 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,...)");
1086: PetscCall(MatCreateVecs(shell->axpy, shell->axpy_right ? NULL : &shell->axpy_right, shell->axpy_left ? NULL : &shell->axpy_left));
1087: PetscCall(VecCopy(x, shell->axpy_left));
1088: PetscCall(MatMultTranspose(shell->axpy, shell->axpy_left, shell->axpy_right));
1089: PetscCall(VecAXPY(y, shell->axpy_vscale, shell->axpy_right));
1090: }
1091: PetscFunctionReturn(PETSC_SUCCESS);
1092: }
1094: static PetscErrorCode MatMultTransposeAdd_Shell(Mat A, Vec x, Vec y, Vec z)
1095: {
1096: Mat_Shell *shell = (Mat_Shell *)A->data;
1098: PetscFunctionBegin;
1099: if (y == z) {
1100: if (!shell->left_add_work) PetscCall(VecDuplicate(z, &shell->left_add_work));
1101: PetscCall(MatMultTranspose(A, x, shell->left_add_work));
1102: PetscCall(VecAXPY(z, 1.0, shell->left_add_work));
1103: } else {
1104: PetscCall(MatMultTranspose(A, x, z));
1105: PetscCall(VecAXPY(z, 1.0, y));
1106: }
1107: PetscFunctionReturn(PETSC_SUCCESS);
1108: }
1110: /*
1111: diag(left)(vscale*A + diag(dshift) + vshift I)diag(right)
1112: */
1113: static PetscErrorCode MatGetDiagonal_Shell(Mat A, Vec v)
1114: {
1115: Mat_Shell *shell = (Mat_Shell *)A->data;
1117: PetscFunctionBegin;
1118: if (shell->ops->getdiagonal) {
1119: PetscCall((*shell->ops->getdiagonal)(A, v));
1120: } else SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Must provide shell matrix with routine to return diagonal using\nMatShellSetOperation(S,MATOP_GET_DIAGONAL,...)");
1121: PetscCall(VecScale(v, shell->vscale));
1122: if (shell->dshift) PetscCall(VecAXPY(v, 1.0, shell->dshift));
1123: PetscCall(VecShift(v, shell->vshift));
1124: if (shell->left) PetscCall(VecPointwiseMult(v, v, shell->left));
1125: if (shell->right) PetscCall(VecPointwiseMult(v, v, shell->right));
1126: if (shell->zrows) {
1127: PetscCall(VecScatterBegin(shell->zvals_sct_r, shell->zvals, v, INSERT_VALUES, SCATTER_REVERSE));
1128: PetscCall(VecScatterEnd(shell->zvals_sct_r, shell->zvals, v, INSERT_VALUES, SCATTER_REVERSE));
1129: }
1130: if (shell->axpy) {
1131: Mat X;
1132: PetscObjectState axpy_state;
1134: PetscCall(MatShellGetContext(shell->axpy, &X));
1135: PetscCall(PetscObjectStateGet((PetscObject)X, &axpy_state));
1136: 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,...)");
1137: PetscCall(MatCreateVecs(shell->axpy, NULL, shell->axpy_left ? NULL : &shell->axpy_left));
1138: PetscCall(MatGetDiagonal(shell->axpy, shell->axpy_left));
1139: PetscCall(VecAXPY(v, shell->axpy_vscale, shell->axpy_left));
1140: }
1141: PetscFunctionReturn(PETSC_SUCCESS);
1142: }
1144: static PetscErrorCode MatShift_Shell(Mat Y, PetscScalar a)
1145: {
1146: Mat_Shell *shell = (Mat_Shell *)Y->data;
1147: PetscBool flg;
1149: PetscFunctionBegin;
1150: PetscCall(MatHasCongruentLayouts(Y, &flg));
1151: PetscCheck(flg, PetscObjectComm((PetscObject)Y), PETSC_ERR_SUP, "Cannot shift shell matrix if it is not congruent");
1152: if (shell->left || shell->right) {
1153: if (!shell->dshift) {
1154: PetscCall(VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift));
1155: PetscCall(VecSet(shell->dshift, a));
1156: } else {
1157: if (shell->left) PetscCall(VecPointwiseMult(shell->dshift, shell->dshift, shell->left));
1158: if (shell->right) PetscCall(VecPointwiseMult(shell->dshift, shell->dshift, shell->right));
1159: PetscCall(VecShift(shell->dshift, a));
1160: }
1161: if (shell->left) PetscCall(VecPointwiseDivide(shell->dshift, shell->dshift, shell->left));
1162: if (shell->right) PetscCall(VecPointwiseDivide(shell->dshift, shell->dshift, shell->right));
1163: } else shell->vshift += a;
1164: if (shell->zrows) PetscCall(VecShift(shell->zvals, a));
1165: PetscFunctionReturn(PETSC_SUCCESS);
1166: }
1168: static PetscErrorCode MatDiagonalSet_Shell_Private(Mat A, Vec D, PetscScalar s)
1169: {
1170: Mat_Shell *shell = (Mat_Shell *)A->data;
1172: PetscFunctionBegin;
1173: if (!shell->dshift) PetscCall(VecDuplicate(D, &shell->dshift));
1174: if (shell->left || shell->right) {
1175: if (!shell->right_work) PetscCall(VecDuplicate(shell->left ? shell->left : shell->right, &shell->right_work));
1176: if (shell->left && shell->right) {
1177: PetscCall(VecPointwiseDivide(shell->right_work, D, shell->left));
1178: PetscCall(VecPointwiseDivide(shell->right_work, shell->right_work, shell->right));
1179: } else if (shell->left) {
1180: PetscCall(VecPointwiseDivide(shell->right_work, D, shell->left));
1181: } else {
1182: PetscCall(VecPointwiseDivide(shell->right_work, D, shell->right));
1183: }
1184: PetscCall(VecAXPY(shell->dshift, s, shell->right_work));
1185: } else {
1186: PetscCall(VecAXPY(shell->dshift, s, D));
1187: }
1188: PetscFunctionReturn(PETSC_SUCCESS);
1189: }
1191: static PetscErrorCode MatDiagonalSet_Shell(Mat A, Vec D, InsertMode ins)
1192: {
1193: Mat_Shell *shell = (Mat_Shell *)A->data;
1194: Vec d;
1195: PetscBool flg;
1197: PetscFunctionBegin;
1198: PetscCall(MatHasCongruentLayouts(A, &flg));
1199: PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Cannot diagonal set or shift shell matrix if it is not congruent");
1200: if (ins == INSERT_VALUES) {
1201: PetscCall(VecDuplicate(D, &d));
1202: PetscCall(MatGetDiagonal(A, d));
1203: PetscCall(MatDiagonalSet_Shell_Private(A, d, -1.));
1204: PetscCall(MatDiagonalSet_Shell_Private(A, D, 1.));
1205: PetscCall(VecDestroy(&d));
1206: if (shell->zrows) PetscCall(VecCopy(D, shell->zvals));
1207: } else {
1208: PetscCall(MatDiagonalSet_Shell_Private(A, D, 1.));
1209: if (shell->zrows) PetscCall(VecAXPY(shell->zvals, 1.0, D));
1210: }
1211: PetscFunctionReturn(PETSC_SUCCESS);
1212: }
1214: static PetscErrorCode MatScale_Shell(Mat Y, PetscScalar a)
1215: {
1216: Mat_Shell *shell = (Mat_Shell *)Y->data;
1218: PetscFunctionBegin;
1219: shell->vscale *= a;
1220: shell->vshift *= a;
1221: if (shell->dshift) PetscCall(VecScale(shell->dshift, a));
1222: shell->axpy_vscale *= a;
1223: if (shell->zrows) PetscCall(VecScale(shell->zvals, a));
1224: PetscFunctionReturn(PETSC_SUCCESS);
1225: }
1227: static PetscErrorCode MatDiagonalScale_Shell(Mat Y, Vec left, Vec right)
1228: {
1229: Mat_Shell *shell = (Mat_Shell *)Y->data;
1231: PetscFunctionBegin;
1232: if (left) {
1233: if (!shell->left) {
1234: PetscCall(VecDuplicate(left, &shell->left));
1235: PetscCall(VecCopy(left, shell->left));
1236: } else {
1237: PetscCall(VecPointwiseMult(shell->left, shell->left, left));
1238: }
1239: if (shell->zrows) PetscCall(VecPointwiseMult(shell->zvals, shell->zvals, left));
1240: }
1241: if (right) {
1242: if (!shell->right) {
1243: PetscCall(VecDuplicate(right, &shell->right));
1244: PetscCall(VecCopy(right, shell->right));
1245: } else {
1246: PetscCall(VecPointwiseMult(shell->right, shell->right, right));
1247: }
1248: if (shell->zrows) {
1249: if (!shell->left_work) PetscCall(MatCreateVecs(Y, NULL, &shell->left_work));
1250: PetscCall(VecSet(shell->zvals_w, 1.0));
1251: PetscCall(VecScatterBegin(shell->zvals_sct_c, right, shell->zvals_w, INSERT_VALUES, SCATTER_FORWARD));
1252: PetscCall(VecScatterEnd(shell->zvals_sct_c, right, shell->zvals_w, INSERT_VALUES, SCATTER_FORWARD));
1253: PetscCall(VecPointwiseMult(shell->zvals, shell->zvals, shell->zvals_w));
1254: }
1255: }
1256: if (shell->axpy) PetscCall(MatDiagonalScale(shell->axpy, left, right));
1257: PetscFunctionReturn(PETSC_SUCCESS);
1258: }
1260: static PetscErrorCode MatAssemblyEnd_Shell(Mat Y, MatAssemblyType t)
1261: {
1262: Mat_Shell *shell = (Mat_Shell *)Y->data;
1264: PetscFunctionBegin;
1265: if (t == MAT_FINAL_ASSEMBLY) {
1266: shell->vshift = 0.0;
1267: shell->vscale = 1.0;
1268: shell->axpy_vscale = 0.0;
1269: shell->axpy_state = 0;
1270: PetscCall(VecDestroy(&shell->dshift));
1271: PetscCall(VecDestroy(&shell->left));
1272: PetscCall(VecDestroy(&shell->right));
1273: PetscCall(MatDestroy(&shell->axpy));
1274: PetscCall(VecDestroy(&shell->axpy_left));
1275: PetscCall(VecDestroy(&shell->axpy_right));
1276: PetscCall(VecScatterDestroy(&shell->zvals_sct_c));
1277: PetscCall(VecScatterDestroy(&shell->zvals_sct_r));
1278: PetscCall(ISDestroy(&shell->zrows));
1279: PetscCall(ISDestroy(&shell->zcols));
1280: }
1281: PetscFunctionReturn(PETSC_SUCCESS);
1282: }
1284: static PetscErrorCode MatMissingDiagonal_Shell(Mat A, PetscBool *missing, PetscInt *d)
1285: {
1286: PetscFunctionBegin;
1287: *missing = PETSC_FALSE;
1288: PetscFunctionReturn(PETSC_SUCCESS);
1289: }
1291: static PetscErrorCode MatAXPY_Shell(Mat Y, PetscScalar a, Mat X, MatStructure str)
1292: {
1293: Mat_Shell *shell = (Mat_Shell *)Y->data;
1295: PetscFunctionBegin;
1296: if (X == Y) {
1297: PetscCall(MatScale(Y, 1.0 + a));
1298: PetscFunctionReturn(PETSC_SUCCESS);
1299: }
1300: if (!shell->axpy) {
1301: PetscCall(MatConvertFrom_Shell(X, MATSHELL, MAT_INITIAL_MATRIX, &shell->axpy));
1302: shell->axpy_vscale = a;
1303: PetscCall(PetscObjectStateGet((PetscObject)X, &shell->axpy_state));
1304: } else {
1305: PetscCall(MatAXPY(shell->axpy, a / shell->axpy_vscale, X, str));
1306: }
1307: PetscFunctionReturn(PETSC_SUCCESS);
1308: }
1310: static struct _MatOps MatOps_Values = {NULL,
1311: NULL,
1312: NULL,
1313: NULL,
1314: /* 4*/ MatMultAdd_Shell,
1315: NULL,
1316: MatMultTransposeAdd_Shell,
1317: NULL,
1318: NULL,
1319: NULL,
1320: /*10*/ NULL,
1321: NULL,
1322: NULL,
1323: NULL,
1324: NULL,
1325: /*15*/ NULL,
1326: NULL,
1327: NULL,
1328: MatDiagonalScale_Shell,
1329: NULL,
1330: /*20*/ NULL,
1331: MatAssemblyEnd_Shell,
1332: NULL,
1333: NULL,
1334: /*24*/ MatZeroRows_Shell,
1335: NULL,
1336: NULL,
1337: NULL,
1338: NULL,
1339: /*29*/ NULL,
1340: NULL,
1341: NULL,
1342: NULL,
1343: NULL,
1344: /*34*/ MatDuplicate_Shell,
1345: NULL,
1346: NULL,
1347: NULL,
1348: NULL,
1349: /*39*/ MatAXPY_Shell,
1350: NULL,
1351: NULL,
1352: NULL,
1353: MatCopy_Shell,
1354: /*44*/ NULL,
1355: MatScale_Shell,
1356: MatShift_Shell,
1357: MatDiagonalSet_Shell,
1358: MatZeroRowsColumns_Shell,
1359: /*49*/ NULL,
1360: NULL,
1361: NULL,
1362: NULL,
1363: NULL,
1364: /*54*/ NULL,
1365: NULL,
1366: NULL,
1367: NULL,
1368: NULL,
1369: /*59*/ NULL,
1370: MatDestroy_Shell,
1371: NULL,
1372: MatConvertFrom_Shell,
1373: NULL,
1374: /*64*/ NULL,
1375: NULL,
1376: NULL,
1377: NULL,
1378: NULL,
1379: /*69*/ NULL,
1380: NULL,
1381: MatConvert_Shell,
1382: NULL,
1383: NULL,
1384: /*74*/ NULL,
1385: NULL,
1386: NULL,
1387: NULL,
1388: NULL,
1389: /*79*/ NULL,
1390: NULL,
1391: NULL,
1392: NULL,
1393: NULL,
1394: /*84*/ NULL,
1395: NULL,
1396: NULL,
1397: NULL,
1398: NULL,
1399: /*89*/ NULL,
1400: NULL,
1401: NULL,
1402: NULL,
1403: NULL,
1404: /*94*/ NULL,
1405: NULL,
1406: NULL,
1407: NULL,
1408: NULL,
1409: /*99*/ NULL,
1410: NULL,
1411: NULL,
1412: NULL,
1413: NULL,
1414: /*104*/ NULL,
1415: NULL,
1416: NULL,
1417: NULL,
1418: NULL,
1419: /*109*/ NULL,
1420: NULL,
1421: NULL,
1422: NULL,
1423: MatMissingDiagonal_Shell,
1424: /*114*/ NULL,
1425: NULL,
1426: NULL,
1427: NULL,
1428: NULL,
1429: /*119*/ NULL,
1430: NULL,
1431: NULL,
1432: NULL,
1433: NULL,
1434: /*124*/ NULL,
1435: NULL,
1436: NULL,
1437: NULL,
1438: NULL,
1439: /*129*/ NULL,
1440: NULL,
1441: NULL,
1442: NULL,
1443: NULL,
1444: /*134*/ NULL,
1445: NULL,
1446: NULL,
1447: NULL,
1448: NULL,
1449: /*139*/ NULL,
1450: NULL,
1451: NULL,
1452: NULL,
1453: NULL,
1454: /*144*/ NULL,
1455: NULL,
1456: NULL,
1457: NULL,
1458: NULL,
1459: NULL,
1460: /*150*/ NULL,
1461: NULL};
1463: static PetscErrorCode MatShellSetContext_Shell(Mat mat, void *ctx)
1464: {
1465: Mat_Shell *shell = (Mat_Shell *)mat->data;
1467: PetscFunctionBegin;
1468: if (ctx) {
1469: PetscContainer ctxcontainer;
1470: PetscCall(PetscContainerCreate(PetscObjectComm((PetscObject)mat), &ctxcontainer));
1471: PetscCall(PetscContainerSetPointer(ctxcontainer, ctx));
1472: PetscCall(PetscObjectCompose((PetscObject)mat, "MatShell ctx", (PetscObject)ctxcontainer));
1473: shell->ctxcontainer = ctxcontainer;
1474: PetscCall(PetscContainerDestroy(&ctxcontainer));
1475: } else {
1476: PetscCall(PetscObjectCompose((PetscObject)mat, "MatShell ctx", NULL));
1477: shell->ctxcontainer = NULL;
1478: }
1480: PetscFunctionReturn(PETSC_SUCCESS);
1481: }
1483: static PetscErrorCode MatShellSetContextDestroy_Shell(Mat mat, PetscErrorCode (*f)(void *))
1484: {
1485: Mat_Shell *shell = (Mat_Shell *)mat->data;
1487: PetscFunctionBegin;
1488: if (shell->ctxcontainer) PetscCall(PetscContainerSetUserDestroy(shell->ctxcontainer, f));
1489: PetscFunctionReturn(PETSC_SUCCESS);
1490: }
1492: static PetscErrorCode MatShellSetVecType_Shell(Mat mat, VecType vtype)
1493: {
1494: PetscFunctionBegin;
1495: PetscCall(PetscFree(mat->defaultvectype));
1496: PetscCall(PetscStrallocpy(vtype, (char **)&mat->defaultvectype));
1497: PetscFunctionReturn(PETSC_SUCCESS);
1498: }
1500: static PetscErrorCode MatShellSetManageScalingShifts_Shell(Mat A)
1501: {
1502: Mat_Shell *shell = (Mat_Shell *)A->data;
1504: PetscFunctionBegin;
1505: shell->managescalingshifts = PETSC_FALSE;
1506: A->ops->diagonalset = NULL;
1507: A->ops->diagonalscale = NULL;
1508: A->ops->scale = NULL;
1509: A->ops->shift = NULL;
1510: A->ops->axpy = NULL;
1511: PetscFunctionReturn(PETSC_SUCCESS);
1512: }
1514: static PetscErrorCode MatShellSetOperation_Shell(Mat mat, MatOperation op, void (*f)(void))
1515: {
1516: Mat_Shell *shell = (Mat_Shell *)mat->data;
1518: PetscFunctionBegin;
1519: switch (op) {
1520: case MATOP_DESTROY:
1521: shell->ops->destroy = (PetscErrorCode(*)(Mat))f;
1522: break;
1523: case MATOP_VIEW:
1524: if (!mat->ops->viewnative) mat->ops->viewnative = mat->ops->view;
1525: mat->ops->view = (PetscErrorCode(*)(Mat, PetscViewer))f;
1526: break;
1527: case MATOP_COPY:
1528: shell->ops->copy = (PetscErrorCode(*)(Mat, Mat, MatStructure))f;
1529: break;
1530: case MATOP_DIAGONAL_SET:
1531: case MATOP_DIAGONAL_SCALE:
1532: case MATOP_SHIFT:
1533: case MATOP_SCALE:
1534: case MATOP_AXPY:
1535: case MATOP_ZERO_ROWS:
1536: case MATOP_ZERO_ROWS_COLUMNS:
1537: PetscCheck(!shell->managescalingshifts, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "MATSHELL is managing scalings and shifts, see MatShellSetManageScalingShifts()");
1538: (((void (**)(void))mat->ops)[op]) = f;
1539: break;
1540: case MATOP_GET_DIAGONAL:
1541: if (shell->managescalingshifts) {
1542: shell->ops->getdiagonal = (PetscErrorCode(*)(Mat, Vec))f;
1543: mat->ops->getdiagonal = MatGetDiagonal_Shell;
1544: } else {
1545: shell->ops->getdiagonal = NULL;
1546: mat->ops->getdiagonal = (PetscErrorCode(*)(Mat, Vec))f;
1547: }
1548: break;
1549: case MATOP_MULT:
1550: if (shell->managescalingshifts) {
1551: shell->ops->mult = (PetscErrorCode(*)(Mat, Vec, Vec))f;
1552: mat->ops->mult = MatMult_Shell;
1553: } else {
1554: shell->ops->mult = NULL;
1555: mat->ops->mult = (PetscErrorCode(*)(Mat, Vec, Vec))f;
1556: }
1557: break;
1558: case MATOP_MULT_TRANSPOSE:
1559: if (shell->managescalingshifts) {
1560: shell->ops->multtranspose = (PetscErrorCode(*)(Mat, Vec, Vec))f;
1561: mat->ops->multtranspose = MatMultTranspose_Shell;
1562: } else {
1563: shell->ops->multtranspose = NULL;
1564: mat->ops->multtranspose = (PetscErrorCode(*)(Mat, Vec, Vec))f;
1565: }
1566: break;
1567: default:
1568: (((void (**)(void))mat->ops)[op]) = f;
1569: break;
1570: }
1571: PetscFunctionReturn(PETSC_SUCCESS);
1572: }
1574: static PetscErrorCode MatShellGetOperation_Shell(Mat mat, MatOperation op, void (**f)(void))
1575: {
1576: Mat_Shell *shell = (Mat_Shell *)mat->data;
1578: PetscFunctionBegin;
1579: switch (op) {
1580: case MATOP_DESTROY:
1581: *f = (void (*)(void))shell->ops->destroy;
1582: break;
1583: case MATOP_VIEW:
1584: *f = (void (*)(void))mat->ops->view;
1585: break;
1586: case MATOP_COPY:
1587: *f = (void (*)(void))shell->ops->copy;
1588: break;
1589: case MATOP_DIAGONAL_SET:
1590: case MATOP_DIAGONAL_SCALE:
1591: case MATOP_SHIFT:
1592: case MATOP_SCALE:
1593: case MATOP_AXPY:
1594: case MATOP_ZERO_ROWS:
1595: case MATOP_ZERO_ROWS_COLUMNS:
1596: *f = (((void (**)(void))mat->ops)[op]);
1597: break;
1598: case MATOP_GET_DIAGONAL:
1599: if (shell->ops->getdiagonal) *f = (void (*)(void))shell->ops->getdiagonal;
1600: else *f = (((void (**)(void))mat->ops)[op]);
1601: break;
1602: case MATOP_MULT:
1603: if (shell->ops->mult) *f = (void (*)(void))shell->ops->mult;
1604: else *f = (((void (**)(void))mat->ops)[op]);
1605: break;
1606: case MATOP_MULT_TRANSPOSE:
1607: if (shell->ops->multtranspose) *f = (void (*)(void))shell->ops->multtranspose;
1608: else *f = (((void (**)(void))mat->ops)[op]);
1609: break;
1610: default:
1611: *f = (((void (**)(void))mat->ops)[op]);
1612: }
1613: PetscFunctionReturn(PETSC_SUCCESS);
1614: }
1616: /*MC
1617: MATSHELL - MATSHELL = "shell" - A matrix type to be used to define your own matrix type -- perhaps matrix-free.
1619: Level: advanced
1621: .seealso: [](ch_matrices), `Mat`, `MatCreateShell()`
1622: M*/
1624: PETSC_EXTERN PetscErrorCode MatCreate_Shell(Mat A)
1625: {
1626: Mat_Shell *b;
1628: PetscFunctionBegin;
1629: PetscCall(PetscNew(&b));
1630: A->data = (void *)b;
1631: A->ops[0] = MatOps_Values;
1633: b->ctxcontainer = NULL;
1634: b->vshift = 0.0;
1635: b->vscale = 1.0;
1636: b->managescalingshifts = PETSC_TRUE;
1637: A->assembled = PETSC_TRUE;
1638: A->preallocated = PETSC_FALSE;
1640: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellGetContext_C", MatShellGetContext_Shell));
1641: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetContext_C", MatShellSetContext_Shell));
1642: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetContextDestroy_C", MatShellSetContextDestroy_Shell));
1643: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetVecType_C", MatShellSetVecType_Shell));
1644: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetManageScalingShifts_C", MatShellSetManageScalingShifts_Shell));
1645: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetOperation_C", MatShellSetOperation_Shell));
1646: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellGetOperation_C", MatShellGetOperation_Shell));
1647: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetMatProductOperation_C", MatShellSetMatProductOperation_Shell));
1648: PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATSHELL));
1649: PetscFunctionReturn(PETSC_SUCCESS);
1650: }
1652: /*@C
1653: MatCreateShell - Creates a new matrix of `MatType` `MATSHELL` for use with a user-defined
1654: private data storage format.
1656: Collective
1658: Input Parameters:
1659: + comm - MPI communicator
1660: . m - number of local rows (or `PETSC_DECIDE` to have calculated if `M` is given)
1661: . n - number of local columns (or `PETSC_DECIDE` to have calculated if `N` is given)
1662: . M - number of global rows (may be `PETSC_DETERMINE` to have calculated if `m` is given)
1663: . N - number of global columns (may be `PETSC_DETERMINE` to have calculated if `n` is given)
1664: - ctx - pointer to data needed by the shell matrix routines
1666: Output Parameter:
1667: . A - the matrix
1669: Level: advanced
1671: Example Usage:
1672: .vb
1673: extern PetscErrorCode mult(Mat, Vec, Vec);
1675: MatCreateShell(comm, m, n, M, N, ctx, &mat);
1676: MatShellSetOperation(mat, MATOP_MULT, (void(*)(void))mult);
1677: // Use matrix for operations that have been set
1678: MatDestroy(mat);
1679: .ve
1681: Notes:
1682: The shell matrix type is intended to provide a simple class to use
1683: with `KSP` (such as, for use with matrix-free methods). You should not
1684: use the shell type if you plan to define a complete matrix class.
1686: PETSc requires that matrices and vectors being used for certain
1687: operations are partitioned accordingly. For example, when
1688: creating a shell matrix, `A`, that supports parallel matrix-vector
1689: products using `MatMult`(A,x,y) the user should set the number
1690: of local matrix rows to be the number of local elements of the
1691: corresponding result vector, y. Note that this is information is
1692: required for use of the matrix interface routines, even though
1693: the shell matrix may not actually be physically partitioned.
1694: For example,
1696: .vb
1697: Vec x, y
1698: extern PetscErrorCode mult(Mat,Vec,Vec);
1699: Mat A
1701: VecCreateMPI(comm,PETSC_DECIDE,M,&y);
1702: VecCreateMPI(comm,PETSC_DECIDE,N,&x);
1703: VecGetLocalSize(y,&m);
1704: VecGetLocalSize(x,&n);
1705: MatCreateShell(comm,m,n,M,N,ctx,&A);
1706: MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
1707: MatMult(A,x,y);
1708: MatDestroy(&A);
1709: VecDestroy(&y);
1710: VecDestroy(&x);
1711: .ve
1713: `MATSHELL` handles `MatShift()`, `MatDiagonalSet()`, `MatDiagonalScale()`, `MatAXPY()`, `MatScale()`, `MatZeroRows()` and `MatZeroRowsColumns()`
1714: internally, so these operations cannot be overwritten unless `MatShellSetManageScalingShifts()` is called.
1716: Developer Notes:
1717: For rectangular matrices do all the scalings and shifts make sense?
1719: Regarding shifting and scaling. The general form is
1721: diag(left)(vscale*A + diag(dshift) + vshift I)diag(right)
1723: The order you apply the operations is important. For example if you have a dshift then
1724: apply a MatScale(s) you get s*vscale*A + s*diag(shift). But if you first scale and then shift
1725: you get s*vscale*A + diag(shift)
1727: A is the user provided function.
1729: `KSP`/`PC` uses changes in the `Mat`'s "state" to decide if preconditioners need to be rebuilt `PCSetUp()` only calls the setup() for
1730: for the `PC` implementation if the `Mat` state has increased from the previous call. Thus to get changes in a `MATSHELL` to trigger
1731: an update in the preconditioner you must call `MatAssemblyBegin()` and `MatAssemblyEnd()` or `PetscObjectStateIncrease`((`PetscObject`)mat);
1732: each time the `MATSHELL` matrix has changed.
1734: Matrix product operations (i.e. `MatMat()`, `MatTransposeMat()` etc) can be specified using `MatShellSetMatProductOperation()`
1736: Calling `MatAssemblyBegin()`/`MatAssemblyEnd()` on a `MATSHELL` removes any previously supplied shift and scales that were provided
1737: with `MatDiagonalSet()`, `MatShift()`, `MatScale()`, or `MatDiagonalScale()`.
1739: Fortran Notes:
1740: To use this from Fortran with a `ctx` you must write an interface definition for this
1741: function and for `MatShellGetContext()` that tells Fortran the Fortran derived data type you are passing
1742: in as the `ctx` argument.
1744: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatShellSetOperation()`, `MatHasOperation()`, `MatShellGetContext()`, `MatShellSetContext()`, `MatShellSetManageScalingShifts()`, `MatShellSetMatProductOperation()`
1745: @*/
1746: PetscErrorCode MatCreateShell(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, PetscInt N, void *ctx, Mat *A)
1747: {
1748: PetscFunctionBegin;
1749: PetscCall(MatCreate(comm, A));
1750: PetscCall(MatSetSizes(*A, m, n, M, N));
1751: PetscCall(MatSetType(*A, MATSHELL));
1752: PetscCall(MatShellSetContext(*A, ctx));
1753: PetscCall(MatSetUp(*A));
1754: PetscFunctionReturn(PETSC_SUCCESS);
1755: }
1757: /*@
1758: MatShellSetContext - sets the context for a `MATSHELL` shell matrix
1760: Logically Collective
1762: Input Parameters:
1763: + mat - the `MATSHELL` shell matrix
1764: - ctx - the context
1766: Level: advanced
1768: Fortran Notes:
1769: You must write a Fortran interface definition for this
1770: function that tells Fortran the Fortran derived data type that you are passing in as the `ctx` argument.
1772: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`
1773: @*/
1774: PetscErrorCode MatShellSetContext(Mat mat, void *ctx)
1775: {
1776: PetscFunctionBegin;
1778: PetscTryMethod(mat, "MatShellSetContext_C", (Mat, void *), (mat, ctx));
1779: PetscFunctionReturn(PETSC_SUCCESS);
1780: }
1782: /*@C
1783: MatShellSetContextDestroy - sets the destroy function for a `MATSHELL` shell matrix context
1785: Logically Collective
1787: Input Parameters:
1788: + mat - the shell matrix
1789: - f - the context destroy function
1791: Level: advanced
1793: Note:
1794: If the `MatShell` is never duplicated, the behavior of this function is equivalent
1795: to `MatShellSetOperation`(`Mat`,`MATOP_DESTROY`,f). However, `MatShellSetContextDestroy()`
1796: ensures proper reference counting for the user provided context data in the case that
1797: the `MATSHELL` is duplicated.
1799: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellSetContext()`
1800: @*/
1801: PetscErrorCode MatShellSetContextDestroy(Mat mat, PetscErrorCode (*f)(void *))
1802: {
1803: PetscFunctionBegin;
1805: PetscTryMethod(mat, "MatShellSetContextDestroy_C", (Mat, PetscErrorCode(*)(void *)), (mat, f));
1806: PetscFunctionReturn(PETSC_SUCCESS);
1807: }
1809: /*@C
1810: MatShellSetVecType - Sets the `VecType` of `Vec` returned by `MatCreateVecs()`
1812: Logically Collective
1814: Input Parameters:
1815: + mat - the `MATSHELL` shell matrix
1816: - vtype - type to use for creating vectors
1818: Level: advanced
1820: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateVecs()`
1821: @*/
1822: PetscErrorCode MatShellSetVecType(Mat mat, VecType vtype)
1823: {
1824: PetscFunctionBegin;
1825: PetscTryMethod(mat, "MatShellSetVecType_C", (Mat, VecType), (mat, vtype));
1826: PetscFunctionReturn(PETSC_SUCCESS);
1827: }
1829: /*@
1830: MatShellSetManageScalingShifts - Allows the user to control the scaling and shift operations of the `MATSHELL`. Must be called immediately
1831: after `MatCreateShell()`
1833: Logically Collective
1835: Input Parameter:
1836: . A - the `MATSHELL` shell matrix
1838: Level: advanced
1840: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellSetContext()`, `MatShellSetOperation()`
1841: @*/
1842: PetscErrorCode MatShellSetManageScalingShifts(Mat A)
1843: {
1844: PetscFunctionBegin;
1846: PetscTryMethod(A, "MatShellSetManageScalingShifts_C", (Mat), (A));
1847: PetscFunctionReturn(PETSC_SUCCESS);
1848: }
1850: /*@C
1851: MatShellTestMult - Compares the multiply routine provided to the `MATSHELL` with differencing on a given function.
1853: Logically Collective; No Fortran Support
1855: Input Parameters:
1856: + mat - the `MATSHELL` shell matrix
1857: . f - the function
1858: . base - differences are computed around this vector, see `MatMFFDSetBase()`, for Jacobians this is the point at which the Jacobian is being evaluated
1859: - ctx - an optional context for the function
1861: Output Parameter:
1862: . flg - `PETSC_TRUE` if the multiply is likely correct
1864: Options Database Key:
1865: . -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference
1867: Level: advanced
1869: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMultTranspose()`
1870: @*/
1871: PetscErrorCode MatShellTestMult(Mat mat, PetscErrorCode (*f)(void *, Vec, Vec), Vec base, void *ctx, PetscBool *flg)
1872: {
1873: PetscInt m, n;
1874: Mat mf, Dmf, Dmat, Ddiff;
1875: PetscReal Diffnorm, Dmfnorm;
1876: PetscBool v = PETSC_FALSE, flag = PETSC_TRUE;
1878: PetscFunctionBegin;
1880: PetscCall(PetscOptionsHasName(NULL, ((PetscObject)mat)->prefix, "-mat_shell_test_mult_view", &v));
1881: PetscCall(MatGetLocalSize(mat, &m, &n));
1882: PetscCall(MatCreateMFFD(PetscObjectComm((PetscObject)mat), m, n, PETSC_DECIDE, PETSC_DECIDE, &mf));
1883: PetscCall(MatMFFDSetFunction(mf, f, ctx));
1884: PetscCall(MatMFFDSetBase(mf, base, NULL));
1886: PetscCall(MatComputeOperator(mf, MATAIJ, &Dmf));
1887: PetscCall(MatComputeOperator(mat, MATAIJ, &Dmat));
1889: PetscCall(MatDuplicate(Dmat, MAT_COPY_VALUES, &Ddiff));
1890: PetscCall(MatAXPY(Ddiff, -1.0, Dmf, DIFFERENT_NONZERO_PATTERN));
1891: PetscCall(MatNorm(Ddiff, NORM_FROBENIUS, &Diffnorm));
1892: PetscCall(MatNorm(Dmf, NORM_FROBENIUS, &Dmfnorm));
1893: if (Diffnorm / Dmfnorm > 10 * PETSC_SQRT_MACHINE_EPSILON) {
1894: flag = PETSC_FALSE;
1895: if (v) {
1896: 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)));
1897: PetscCall(MatViewFromOptions(Ddiff, (PetscObject)mat, "-mat_shell_test_mult_view"));
1898: PetscCall(MatViewFromOptions(Dmf, (PetscObject)mat, "-mat_shell_test_mult_view"));
1899: PetscCall(MatViewFromOptions(Dmat, (PetscObject)mat, "-mat_shell_test_mult_view"));
1900: }
1901: } else if (v) {
1902: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)mat), "MATSHELL and matrix-free multiple appear to produce the same results\n"));
1903: }
1904: if (flg) *flg = flag;
1905: PetscCall(MatDestroy(&Ddiff));
1906: PetscCall(MatDestroy(&mf));
1907: PetscCall(MatDestroy(&Dmf));
1908: PetscCall(MatDestroy(&Dmat));
1909: PetscFunctionReturn(PETSC_SUCCESS);
1910: }
1912: /*@C
1913: MatShellTestMultTranspose - Compares the multiply transpose routine provided to the `MATSHELL` with differencing on a given function.
1915: Logically Collective; No Fortran Support
1917: Input Parameters:
1918: + mat - the `MATSHELL` shell matrix
1919: . f - the function
1920: . base - differences are computed around this vector, see `MatMFFDSetBase()`, for Jacobians this is the point at which the Jacobian is being evaluated
1921: - ctx - an optional context for the function
1923: Output Parameter:
1924: . flg - `PETSC_TRUE` if the multiply is likely correct
1926: Options Database Key:
1927: . -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference
1929: Level: advanced
1931: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMult()`
1932: @*/
1933: PetscErrorCode MatShellTestMultTranspose(Mat mat, PetscErrorCode (*f)(void *, Vec, Vec), Vec base, void *ctx, PetscBool *flg)
1934: {
1935: Vec x, y, z;
1936: PetscInt m, n, M, N;
1937: Mat mf, Dmf, Dmat, Ddiff;
1938: PetscReal Diffnorm, Dmfnorm;
1939: PetscBool v = PETSC_FALSE, flag = PETSC_TRUE;
1941: PetscFunctionBegin;
1943: PetscCall(PetscOptionsHasName(NULL, ((PetscObject)mat)->prefix, "-mat_shell_test_mult_transpose_view", &v));
1944: PetscCall(MatCreateVecs(mat, &x, &y));
1945: PetscCall(VecDuplicate(y, &z));
1946: PetscCall(MatGetLocalSize(mat, &m, &n));
1947: PetscCall(MatGetSize(mat, &M, &N));
1948: PetscCall(MatCreateMFFD(PetscObjectComm((PetscObject)mat), m, n, M, N, &mf));
1949: PetscCall(MatMFFDSetFunction(mf, f, ctx));
1950: PetscCall(MatMFFDSetBase(mf, base, NULL));
1951: PetscCall(MatComputeOperator(mf, MATAIJ, &Dmf));
1952: PetscCall(MatTranspose(Dmf, MAT_INPLACE_MATRIX, &Dmf));
1953: PetscCall(MatComputeOperatorTranspose(mat, MATAIJ, &Dmat));
1955: PetscCall(MatDuplicate(Dmat, MAT_COPY_VALUES, &Ddiff));
1956: PetscCall(MatAXPY(Ddiff, -1.0, Dmf, DIFFERENT_NONZERO_PATTERN));
1957: PetscCall(MatNorm(Ddiff, NORM_FROBENIUS, &Diffnorm));
1958: PetscCall(MatNorm(Dmf, NORM_FROBENIUS, &Dmfnorm));
1959: if (Diffnorm / Dmfnorm > 10 * PETSC_SQRT_MACHINE_EPSILON) {
1960: flag = PETSC_FALSE;
1961: if (v) {
1962: 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)));
1963: PetscCall(MatViewFromOptions(Ddiff, (PetscObject)mat, "-mat_shell_test_mult_transpose_view"));
1964: PetscCall(MatViewFromOptions(Dmf, (PetscObject)mat, "-mat_shell_test_mult_transpose_view"));
1965: PetscCall(MatViewFromOptions(Dmat, (PetscObject)mat, "-mat_shell_test_mult_transpose_view"));
1966: }
1967: } else if (v) {
1968: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)mat), "MATSHELL transpose and matrix-free multiple appear to produce the same results\n"));
1969: }
1970: if (flg) *flg = flag;
1971: PetscCall(MatDestroy(&mf));
1972: PetscCall(MatDestroy(&Dmat));
1973: PetscCall(MatDestroy(&Ddiff));
1974: PetscCall(MatDestroy(&Dmf));
1975: PetscCall(VecDestroy(&x));
1976: PetscCall(VecDestroy(&y));
1977: PetscCall(VecDestroy(&z));
1978: PetscFunctionReturn(PETSC_SUCCESS);
1979: }
1981: /*@C
1982: MatShellSetOperation - Allows user to set a matrix operation for a `MATSHELL` shell matrix.
1984: Logically Collective
1986: Input Parameters:
1987: + mat - the `MATSHELL` shell matrix
1988: . op - the name of the operation
1989: - g - the function that provides the operation.
1991: Level: advanced
1993: Example Usage:
1994: .vb
1995: extern PetscErrorCode usermult(Mat, Vec, Vec);
1997: MatCreateShell(comm, m, n, M, N, ctx, &A);
1998: MatShellSetOperation(A, MATOP_MULT, (void(*)(void))usermult);
1999: .ve
2001: Notes:
2002: See the file include/petscmat.h for a complete list of matrix
2003: operations, which all have the form MATOP_<OPERATION>, where
2004: <OPERATION> is the name (in all capital letters) of the
2005: user interface routine (e.g., `MatMult()` -> `MATOP_MULT`).
2007: All user-provided functions (except for `MATOP_DESTROY`) should have the same calling
2008: sequence as the usual matrix interface routines, since they
2009: are intended to be accessed via the usual matrix interface
2010: routines, e.g.,
2011: $ MatMult(Mat, Vec, Vec) -> usermult(Mat, Vec, Vec)
2013: In particular each function MUST return an error code of 0 on success and
2014: nonzero on failure.
2016: Within each user-defined routine, the user should call
2017: `MatShellGetContext()` to obtain the user-defined context that was
2018: set by `MatCreateShell()`.
2020: Use `MatSetOperation()` to set an operation for any matrix type. For matrix product operations (i.e. `MatMatXXX()`, `MatTransposeMatXXX()` etc)
2021: use `MatShellSetMatProductOperation()`
2023: Fortran Notes:
2024: For `MatCreateVecs()` the user code should check if the input left or right matrix is -1 and in that case not
2025: generate a matrix. See src/mat/tests/ex120f.F
2027: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellSetContext()`, `MatSetOperation()`, `MatShellSetManageScalingShifts()`, `MatShellSetMatProductOperation()`
2028: @*/
2029: PetscErrorCode MatShellSetOperation(Mat mat, MatOperation op, void (*g)(void))
2030: {
2031: PetscFunctionBegin;
2033: PetscTryMethod(mat, "MatShellSetOperation_C", (Mat, MatOperation, void (*)(void)), (mat, op, g));
2034: PetscFunctionReturn(PETSC_SUCCESS);
2035: }
2037: /*@C
2038: MatShellGetOperation - Gets a matrix function for a `MATSHELL` shell matrix.
2040: Not Collective
2042: Input Parameters:
2043: + mat - the `MATSHELL` shell matrix
2044: - op - the name of the operation
2046: Output Parameter:
2047: . g - the function that provides the operation.
2049: Level: advanced
2051: Notes:
2052: See the file include/petscmat.h for a complete list of matrix
2053: operations, which all have the form MATOP_<OPERATION>, where
2054: <OPERATION> is the name (in all capital letters) of the
2055: user interface routine (e.g., `MatMult()` -> `MATOP_MULT`).
2057: All user-provided functions have the same calling
2058: sequence as the usual matrix interface routines, since they
2059: are intended to be accessed via the usual matrix interface
2060: routines, e.g.,
2061: $ MatMult(Mat, Vec, Vec) -> usermult(Mat, Vec, Vec)
2063: Within each user-defined routine, the user should call
2064: `MatShellGetContext()` to obtain the user-defined context that was
2065: set by `MatCreateShell()`.
2067: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellSetOperation()`, `MatShellSetContext()`
2068: @*/
2069: PetscErrorCode MatShellGetOperation(Mat mat, MatOperation op, void (**g)(void))
2070: {
2071: PetscFunctionBegin;
2073: PetscUseMethod(mat, "MatShellGetOperation_C", (Mat, MatOperation, void (**)(void)), (mat, op, g));
2074: PetscFunctionReturn(PETSC_SUCCESS);
2075: }
2077: /*@
2078: MatIsShell - Inquires if a matrix is derived from `MATSHELL`
2080: Input Parameter:
2081: . mat - the matrix
2083: Output Parameter:
2084: . flg - the boolean value
2086: Level: developer
2088: Developer Notes:
2089: In the future, we should allow the object type name to be changed still using the `MATSHELL` data structure for other matrices
2090: (i.e. `MATTRANSPOSEVIRTUAL`, `MATSCHURCOMPLEMENT` etc)
2092: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MATMFFD`, `MatCreateShell()`, `MATTRANSPOSEVIRTUAL`, `MATSCHURCOMPLEMENT`
2093: @*/
2094: PetscErrorCode MatIsShell(Mat mat, PetscBool *flg)
2095: {
2096: PetscFunctionBegin;
2098: PetscAssertPointer(flg, 2);
2099: *flg = (PetscBool)(mat->ops->destroy == MatDestroy_Shell);
2100: PetscFunctionReturn(PETSC_SUCCESS);
2101: }