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