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