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, PetscCtxRt ctx)
210: {
211: Mat_Shell *shell = (Mat_Shell *)mat->data;
213: PetscFunctionBegin;
214: if (shell->ctxcontainer) PetscCall(PetscContainerGetPointer(shell->ctxcontainer, 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: This only works when the context is a Fortran derived type or a `PetscObject`. Declare `ctx` with
234: .vb
235: type(tUsertype), pointer :: ctx
236: .ve
238: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellSetOperation()`, `MatShellSetContext()`
239: @*/
240: PetscErrorCode MatShellGetContext(Mat mat, PetscCtxRt ctx)
241: {
242: PetscFunctionBegin;
244: PetscAssertPointer(ctx, 2);
245: PetscUseMethod(mat, "MatShellGetContext_C", (Mat, void *), (mat, ctx));
246: PetscFunctionReturn(PETSC_SUCCESS);
247: }
249: static PetscErrorCode MatZeroRowsColumns_Local_Shell(Mat mat, PetscInt nr, PetscInt rows[], PetscInt nc, PetscInt cols[], PetscScalar diag, PetscBool rc)
250: {
251: Mat_Shell *shell = (Mat_Shell *)mat->data;
252: Vec x = NULL, b = NULL;
253: IS is1, is2;
254: const PetscInt *ridxs;
255: PetscInt *idxs, *gidxs;
256: PetscInt cum, rst, cst, i;
258: PetscFunctionBegin;
259: if (!shell->zvals) PetscCall(MatCreateVecs(mat, NULL, &shell->zvals));
260: if (!shell->zvals_w) PetscCall(VecDuplicate(shell->zvals, &shell->zvals_w));
261: PetscCall(MatGetOwnershipRange(mat, &rst, NULL));
262: PetscCall(MatGetOwnershipRangeColumn(mat, &cst, NULL));
264: /* Expand/create index set of zeroed rows */
265: PetscCall(PetscMalloc1(nr, &idxs));
266: for (i = 0; i < nr; i++) idxs[i] = rows[i] + rst;
267: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)mat), nr, idxs, PETSC_OWN_POINTER, &is1));
268: PetscCall(ISSort(is1));
269: PetscCall(VecISSet(shell->zvals, is1, diag));
270: if (shell->zrows) {
271: PetscCall(ISSum(shell->zrows, is1, &is2));
272: PetscCall(ISDestroy(&shell->zrows));
273: PetscCall(ISDestroy(&is1));
274: shell->zrows = is2;
275: } else shell->zrows = is1;
277: /* Create scatters for diagonal values communications */
278: PetscCall(VecScatterDestroy(&shell->zvals_sct_c));
279: PetscCall(VecScatterDestroy(&shell->zvals_sct_r));
281: /* row scatter: from/to left vector */
282: PetscCall(MatCreateVecs(mat, &x, &b));
283: PetscCall(VecScatterCreate(b, shell->zrows, shell->zvals_w, shell->zrows, &shell->zvals_sct_r));
285: /* col scatter: from right vector to left vector */
286: PetscCall(ISGetIndices(shell->zrows, &ridxs));
287: PetscCall(ISGetLocalSize(shell->zrows, &nr));
288: PetscCall(PetscMalloc1(nr, &gidxs));
289: for (i = 0, cum = 0; i < nr; i++) {
290: if (ridxs[i] >= mat->cmap->N) continue;
291: gidxs[cum] = ridxs[i];
292: cum++;
293: }
294: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)mat), cum, gidxs, PETSC_OWN_POINTER, &is1));
295: PetscCall(VecScatterCreate(x, is1, shell->zvals_w, is1, &shell->zvals_sct_c));
296: PetscCall(ISDestroy(&is1));
297: PetscCall(VecDestroy(&x));
298: PetscCall(VecDestroy(&b));
300: /* Expand/create index set of zeroed columns */
301: if (rc) {
302: PetscCall(PetscMalloc1(nc, &idxs));
303: for (i = 0; i < nc; i++) idxs[i] = cols[i] + cst;
304: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)mat), nc, idxs, PETSC_OWN_POINTER, &is1));
305: PetscCall(ISSort(is1));
306: if (shell->zcols) {
307: PetscCall(ISSum(shell->zcols, is1, &is2));
308: PetscCall(ISDestroy(&shell->zcols));
309: PetscCall(ISDestroy(&is1));
310: shell->zcols = is2;
311: } else shell->zcols = is1;
312: }
313: PetscFunctionReturn(PETSC_SUCCESS);
314: }
316: static PetscErrorCode MatZeroRows_Shell(Mat mat, PetscInt n, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
317: {
318: Mat_Shell *shell = (Mat_Shell *)mat->data;
319: PetscInt nr, *lrows;
321: PetscFunctionBegin;
322: if (x && b) {
323: Vec xt;
324: PetscScalar *vals;
325: PetscInt *gcols, i, st, nl, nc;
327: PetscCall(PetscMalloc1(n, &gcols));
328: for (i = 0, nc = 0; i < n; i++)
329: if (rows[i] < mat->cmap->N) gcols[nc++] = rows[i];
331: PetscCall(MatCreateVecs(mat, &xt, NULL));
332: PetscCall(VecCopy(x, xt));
333: PetscCall(PetscCalloc1(nc, &vals));
334: PetscCall(VecSetValues(xt, nc, gcols, vals, INSERT_VALUES)); /* xt = [x1, 0] */
335: PetscCall(PetscFree(vals));
336: PetscCall(VecAssemblyBegin(xt));
337: PetscCall(VecAssemblyEnd(xt));
338: PetscCall(VecAYPX(xt, -1.0, x)); /* xt = [0, x2] */
340: PetscCall(VecGetOwnershipRange(xt, &st, NULL));
341: PetscCall(VecGetLocalSize(xt, &nl));
342: PetscCall(VecGetArray(xt, &vals));
343: for (i = 0; i < nl; i++) {
344: PetscInt g = i + st;
345: if (g > mat->rmap->N) continue;
346: if (PetscAbsScalar(vals[i]) == 0.0) continue;
347: PetscCall(VecSetValue(b, g, diag * vals[i], INSERT_VALUES));
348: }
349: PetscCall(VecRestoreArray(xt, &vals));
350: PetscCall(VecAssemblyBegin(b));
351: PetscCall(VecAssemblyEnd(b)); /* b = [b1, x2 * diag] */
352: PetscCall(VecDestroy(&xt));
353: PetscCall(PetscFree(gcols));
354: }
355: PetscCall(PetscLayoutMapLocal(mat->rmap, n, rows, &nr, &lrows, NULL));
356: PetscCall(MatZeroRowsColumns_Local_Shell(mat, nr, lrows, 0, NULL, diag, PETSC_FALSE));
357: if (shell->axpy) PetscCall(MatZeroRows(shell->axpy, n, rows, 0.0, NULL, NULL));
358: PetscCall(PetscFree(lrows));
359: PetscFunctionReturn(PETSC_SUCCESS);
360: }
362: static PetscErrorCode MatZeroRowsColumns_Shell(Mat mat, PetscInt n, const PetscInt rowscols[], PetscScalar diag, Vec x, Vec b)
363: {
364: Mat_Shell *shell = (Mat_Shell *)mat->data;
365: PetscInt *lrows, *lcols;
366: PetscInt nr, nc;
367: PetscBool congruent;
369: PetscFunctionBegin;
370: if (x && b) {
371: Vec xt, bt;
372: PetscScalar *vals;
373: PetscInt *grows, *gcols, i, st, nl;
375: PetscCall(PetscMalloc2(n, &grows, n, &gcols));
376: for (i = 0, nr = 0; i < n; i++)
377: if (rowscols[i] < mat->rmap->N) grows[nr++] = rowscols[i];
378: for (i = 0, nc = 0; i < n; i++)
379: if (rowscols[i] < mat->cmap->N) gcols[nc++] = rowscols[i];
380: PetscCall(PetscCalloc1(n, &vals));
382: PetscCall(MatCreateVecs(mat, &xt, &bt));
383: PetscCall(VecCopy(x, xt));
384: PetscCall(VecSetValues(xt, nc, gcols, vals, INSERT_VALUES)); /* xt = [x1, 0] */
385: PetscCall(VecAssemblyBegin(xt));
386: PetscCall(VecAssemblyEnd(xt));
387: PetscCall(VecAXPY(xt, -1.0, x)); /* xt = [0, -x2] */
388: PetscCall(MatMult(mat, xt, bt)); /* bt = [-A12*x2,-A22*x2] */
389: PetscCall(VecSetValues(bt, nr, grows, vals, INSERT_VALUES)); /* bt = [-A12*x2,0] */
390: PetscCall(VecAssemblyBegin(bt));
391: PetscCall(VecAssemblyEnd(bt));
392: PetscCall(VecAXPY(b, 1.0, bt)); /* b = [b1 - A12*x2, b2] */
393: PetscCall(VecSetValues(bt, nr, grows, vals, INSERT_VALUES)); /* b = [b1 - A12*x2, 0] */
394: PetscCall(VecAssemblyBegin(bt));
395: PetscCall(VecAssemblyEnd(bt));
396: PetscCall(PetscFree(vals));
398: PetscCall(VecGetOwnershipRange(xt, &st, NULL));
399: PetscCall(VecGetLocalSize(xt, &nl));
400: PetscCall(VecGetArray(xt, &vals));
401: for (i = 0; i < nl; i++) {
402: PetscInt g = i + st;
403: if (g > mat->rmap->N) continue;
404: if (PetscAbsScalar(vals[i]) == 0.0) continue;
405: PetscCall(VecSetValue(b, g, -diag * vals[i], INSERT_VALUES));
406: }
407: PetscCall(VecRestoreArray(xt, &vals));
408: PetscCall(VecAssemblyBegin(b));
409: PetscCall(VecAssemblyEnd(b)); /* b = [b1 - A12*x2, x2 * diag] */
410: PetscCall(VecDestroy(&xt));
411: PetscCall(VecDestroy(&bt));
412: PetscCall(PetscFree2(grows, gcols));
413: }
414: PetscCall(PetscLayoutMapLocal(mat->rmap, n, rowscols, &nr, &lrows, NULL));
415: PetscCall(MatHasCongruentLayouts(mat, &congruent));
416: if (congruent) {
417: nc = nr;
418: lcols = lrows;
419: } else { /* MatZeroRowsColumns implicitly assumes the rowscols indices are for a square matrix, here we handle a more general case */
420: PetscInt i, nt, *t;
422: PetscCall(PetscMalloc1(n, &t));
423: for (i = 0, nt = 0; i < n; i++)
424: if (rowscols[i] < mat->cmap->N) t[nt++] = rowscols[i];
425: PetscCall(PetscLayoutMapLocal(mat->cmap, nt, t, &nc, &lcols, NULL));
426: PetscCall(PetscFree(t));
427: }
428: PetscCall(MatZeroRowsColumns_Local_Shell(mat, nr, lrows, nc, lcols, diag, PETSC_TRUE));
429: if (!congruent) PetscCall(PetscFree(lcols));
430: PetscCall(PetscFree(lrows));
431: if (shell->axpy) PetscCall(MatZeroRowsColumns(shell->axpy, n, rowscols, 0.0, NULL, NULL));
432: PetscFunctionReturn(PETSC_SUCCESS);
433: }
435: static PetscErrorCode MatDestroy_Shell(Mat mat)
436: {
437: Mat_Shell *shell = (Mat_Shell *)mat->data;
438: MatShellMatFunctionList matmat;
440: PetscFunctionBegin;
441: if (shell->ops->destroy) PetscCall((*shell->ops->destroy)(mat));
442: PetscCall(PetscMemzero(shell->ops, sizeof(struct _MatShellOps)));
443: PetscCall(VecDestroy(&shell->left));
444: PetscCall(VecDestroy(&shell->right));
445: PetscCall(VecDestroy(&shell->dshift));
446: PetscCall(VecDestroy(&shell->left_work));
447: PetscCall(VecDestroy(&shell->right_work));
448: PetscCall(VecDestroy(&shell->left_add_work));
449: PetscCall(VecDestroy(&shell->right_add_work));
450: PetscCall(VecDestroy(&shell->axpy_left));
451: PetscCall(VecDestroy(&shell->axpy_right));
452: PetscCall(MatDestroy(&shell->axpy));
453: PetscCall(VecDestroy(&shell->zvals_w));
454: PetscCall(VecDestroy(&shell->zvals));
455: PetscCall(VecScatterDestroy(&shell->zvals_sct_c));
456: PetscCall(VecScatterDestroy(&shell->zvals_sct_r));
457: PetscCall(ISDestroy(&shell->zrows));
458: PetscCall(ISDestroy(&shell->zcols));
460: matmat = shell->matmat;
461: while (matmat) {
462: MatShellMatFunctionList next = matmat->next;
464: PetscCall(PetscObjectComposeFunction((PetscObject)mat, matmat->composedname, NULL));
465: PetscCall(PetscFree(matmat->composedname));
466: PetscCall(PetscFree(matmat->resultname));
467: PetscCall(PetscFree(matmat));
468: matmat = next;
469: }
470: PetscCall(MatShellSetContext(mat, NULL));
471: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatShellGetContext_C", NULL));
472: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatShellSetContext_C", NULL));
473: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatShellSetContextDestroy_C", NULL));
474: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatShellSetVecType_C", NULL));
475: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatShellSetManageScalingShifts_C", NULL));
476: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatShellGetScalingShifts_C", NULL));
477: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatShellSetOperation_C", NULL));
478: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatShellGetOperation_C", NULL));
479: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatShellSetMatProductOperation_C", NULL));
480: PetscCall(PetscFree(mat->data));
481: PetscFunctionReturn(PETSC_SUCCESS);
482: }
484: typedef struct {
485: PetscErrorCode (*numeric)(Mat, Mat, Mat, void *);
486: PetscCtxDestroyFn *destroy;
487: void *ctx;
488: Mat B;
489: Mat Bt;
490: Mat axpy;
491: } MatProductCtx_MatMatShell;
493: static PetscErrorCode MatProductCtxDestroy_MatMatShell(PetscCtxRt data)
494: {
495: MatProductCtx_MatMatShell *mmdata = *(MatProductCtx_MatMatShell **)data;
497: PetscFunctionBegin;
498: if (mmdata->destroy) PetscCall((*mmdata->destroy)(&mmdata->ctx));
499: PetscCall(MatDestroy(&mmdata->B));
500: PetscCall(MatDestroy(&mmdata->Bt));
501: PetscCall(MatDestroy(&mmdata->axpy));
502: PetscCall(PetscFree(mmdata));
503: PetscFunctionReturn(PETSC_SUCCESS);
504: }
506: static PetscErrorCode MatProductNumeric_Shell_X(Mat D)
507: {
508: Mat_Product *product;
509: Mat A, B;
510: MatProductCtx_MatMatShell *mdata;
511: Mat_Shell *shell;
512: PetscBool useBmdata = PETSC_FALSE, newB = PETSC_TRUE;
513: PetscErrorCode (*stashsym)(Mat), (*stashnum)(Mat);
515: PetscFunctionBegin;
516: MatCheckProduct(D, 1);
517: product = D->product;
518: PetscCheck(product->data, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "Product data empty");
519: A = product->A;
520: B = product->B;
521: mdata = (MatProductCtx_MatMatShell *)product->data;
522: PetscCheck(mdata->numeric, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "Missing numeric operation");
523: shell = (Mat_Shell *)A->data;
524: stashsym = D->ops->productsymbolic;
525: stashnum = D->ops->productnumeric;
526: if (shell->managescalingshifts) {
527: PetscCheck(!shell->zcols && !shell->zrows, PetscObjectComm((PetscObject)D), PETSC_ERR_SUP, "MatProduct not supported with zeroed rows/columns");
528: if (shell->right || shell->left) {
529: useBmdata = PETSC_TRUE;
530: if (!mdata->B) PetscCall(MatDuplicate(B, MAT_SHARE_NONZERO_PATTERN, &mdata->B));
531: else newB = PETSC_FALSE;
532: PetscCall(MatCopy(B, mdata->B, SAME_NONZERO_PATTERN));
533: }
534: switch (product->type) {
535: case MATPRODUCT_AB: /* s L A R B + v L R B + L D R B */
536: if (shell->right) PetscCall(MatDiagonalScale(mdata->B, shell->right, NULL));
537: break;
538: case MATPRODUCT_AtB: /* s R A^t L B + v R L B + R D L B */
539: if (shell->left) PetscCall(MatDiagonalScale(mdata->B, shell->left, NULL));
540: break;
541: case MATPRODUCT_ABt: /* s L A R B^t + v L R B^t + L D R B^t */
542: if (shell->right) PetscCall(MatDiagonalScale(mdata->B, NULL, shell->right));
543: break;
544: case MATPRODUCT_RARt: /* s B L A R B^t + v B L R B^t + B L D R B^t */
545: if (shell->right && shell->left) {
546: PetscBool flg;
548: PetscCall(VecEqual(shell->right, shell->left, &flg));
549: 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,
550: ((PetscObject)B)->type_name);
551: }
552: if (shell->right) PetscCall(MatDiagonalScale(mdata->B, NULL, shell->right));
553: break;
554: case MATPRODUCT_PtAP: /* s B^t L A R B + v B^t L R B + B^t L D R B */
555: if (shell->right && shell->left) {
556: PetscBool flg;
558: PetscCall(VecEqual(shell->right, shell->left, &flg));
559: 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,
560: ((PetscObject)B)->type_name);
561: }
562: if (shell->right) PetscCall(MatDiagonalScale(mdata->B, shell->right, NULL));
563: break;
564: default:
565: 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);
566: }
567: }
568: /* allow the user to call MatMat operations on D */
569: D->product = NULL;
570: D->ops->productsymbolic = NULL;
571: D->ops->productnumeric = NULL;
573: PetscCall((*mdata->numeric)(A, useBmdata ? mdata->B : B, D, mdata->ctx));
575: /* clear any leftover user data and restore D pointers */
576: PetscCall(MatProductClear(D));
577: D->ops->productsymbolic = stashsym;
578: D->ops->productnumeric = stashnum;
579: D->product = product;
581: if (shell->managescalingshifts) {
582: PetscCall(MatScale(D, shell->vscale));
583: switch (product->type) {
584: case MATPRODUCT_AB: /* s L A R B + v L R B + L D R B */
585: case MATPRODUCT_ABt: /* s L A R B^t + v L R B^t + L D R B^t */
586: if (shell->left) {
587: PetscCall(MatDiagonalScale(D, shell->left, NULL));
588: if (shell->dshift || shell->vshift != 0.0) {
589: if (!shell->left_work) PetscCall(MatCreateVecs(A, NULL, &shell->left_work));
590: if (shell->dshift) {
591: PetscCall(VecCopy(shell->dshift, shell->left_work));
592: PetscCall(VecShift(shell->left_work, shell->vshift));
593: PetscCall(VecPointwiseMult(shell->left_work, shell->left_work, shell->left));
594: } else {
595: PetscCall(VecSet(shell->left_work, shell->vshift));
596: }
597: if (product->type == MATPRODUCT_ABt) {
598: MatReuse reuse = mdata->Bt ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX;
599: MatStructure str = mdata->Bt ? SUBSET_NONZERO_PATTERN : DIFFERENT_NONZERO_PATTERN;
601: PetscCall(MatTranspose(mdata->B, reuse, &mdata->Bt));
602: PetscCall(MatDiagonalScale(mdata->Bt, shell->left_work, NULL));
603: PetscCall(MatAXPY(D, 1.0, mdata->Bt, str));
604: } else {
605: MatStructure str = newB ? DIFFERENT_NONZERO_PATTERN : SUBSET_NONZERO_PATTERN;
607: PetscCall(MatDiagonalScale(mdata->B, shell->left_work, NULL));
608: PetscCall(MatAXPY(D, 1.0, mdata->B, str));
609: }
610: }
611: }
612: break;
613: case MATPRODUCT_AtB: /* s R A^t L B + v R L B + R D L B */
614: if (shell->right) {
615: PetscCall(MatDiagonalScale(D, shell->right, NULL));
616: if (shell->dshift || shell->vshift != 0.0) {
617: MatStructure str = newB ? DIFFERENT_NONZERO_PATTERN : SUBSET_NONZERO_PATTERN;
619: if (!shell->right_work) PetscCall(MatCreateVecs(A, &shell->right_work, NULL));
620: if (shell->dshift) {
621: PetscCall(VecCopy(shell->dshift, shell->right_work));
622: PetscCall(VecShift(shell->right_work, shell->vshift));
623: PetscCall(VecPointwiseMult(shell->right_work, shell->right_work, shell->right));
624: } else {
625: PetscCall(VecSet(shell->right_work, shell->vshift));
626: }
627: PetscCall(MatDiagonalScale(mdata->B, shell->right_work, NULL));
628: PetscCall(MatAXPY(D, 1.0, mdata->B, str));
629: }
630: }
631: break;
632: case MATPRODUCT_PtAP: /* s B^t L A R B + v B^t L R B + B^t L D R B */
633: case MATPRODUCT_RARt: /* s B L A R B^t + v B L R B^t + B L D R B^t */
634: PetscCheck(!shell->dshift && shell->vshift == 0.0, 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,
635: ((PetscObject)B)->type_name);
636: break;
637: default:
638: 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);
639: }
640: if (shell->axpy && shell->axpy_vscale != 0.0) {
641: Mat X;
642: PetscObjectState axpy_state;
643: MatStructure str = DIFFERENT_NONZERO_PATTERN; /* not sure it is safe to ever use SUBSET_NONZERO_PATTERN */
645: PetscCall(MatShellGetContext(shell->axpy, &X));
646: PetscCall(PetscObjectStateGet((PetscObject)X, &axpy_state));
647: 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,...)");
648: if (!mdata->axpy) {
649: str = DIFFERENT_NONZERO_PATTERN;
650: PetscCall(MatProductCreate(shell->axpy, B, NULL, &mdata->axpy));
651: PetscCall(MatProductSetType(mdata->axpy, product->type));
652: PetscCall(MatProductSetFromOptions(mdata->axpy));
653: PetscCall(MatProductSymbolic(mdata->axpy));
654: } else { /* May be that shell->axpy has changed */
655: PetscBool flg;
657: PetscCall(MatProductReplaceMats(shell->axpy, B, NULL, mdata->axpy));
658: PetscCall(MatHasOperation(mdata->axpy, MATOP_PRODUCTSYMBOLIC, &flg));
659: if (!flg) {
660: str = DIFFERENT_NONZERO_PATTERN;
661: PetscCall(MatProductSetFromOptions(mdata->axpy));
662: PetscCall(MatProductSymbolic(mdata->axpy));
663: }
664: }
665: PetscCall(MatProductNumeric(mdata->axpy));
666: PetscCall(MatAXPY(D, shell->axpy_vscale, mdata->axpy, str));
667: }
668: }
669: PetscFunctionReturn(PETSC_SUCCESS);
670: }
672: static PetscErrorCode MatProductSymbolic_Shell_X(Mat D)
673: {
674: Mat_Product *product;
675: Mat A, B;
676: MatShellMatFunctionList matmat;
677: Mat_Shell *shell;
678: PetscBool flg = PETSC_FALSE;
679: char composedname[256];
680: MatProductCtx_MatMatShell *mdata;
682: PetscFunctionBegin;
683: MatCheckProduct(D, 1);
684: product = D->product;
685: PetscCheck(!product->data, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "Product data not empty");
686: A = product->A;
687: B = product->B;
688: shell = (Mat_Shell *)A->data;
689: matmat = shell->matmat;
690: PetscCall(PetscSNPrintf(composedname, sizeof(composedname), "MatProductSetFromOptions_%s_%s_C", ((PetscObject)A)->type_name, ((PetscObject)B)->type_name));
691: while (matmat) {
692: PetscCall(PetscStrcmp(composedname, matmat->composedname, &flg));
693: flg = (PetscBool)(flg && (matmat->ptype == product->type));
694: if (flg) break;
695: matmat = matmat->next;
696: }
697: PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Composedname \"%s\" for product type %s not found", composedname, MatProductTypes[product->type]);
698: switch (product->type) {
699: case MATPRODUCT_AB:
700: PetscCall(MatSetSizes(D, A->rmap->n, B->cmap->n, A->rmap->N, B->cmap->N));
701: break;
702: case MATPRODUCT_AtB:
703: PetscCall(MatSetSizes(D, A->cmap->n, B->cmap->n, A->cmap->N, B->cmap->N));
704: break;
705: case MATPRODUCT_ABt:
706: PetscCall(MatSetSizes(D, A->rmap->n, B->rmap->n, A->rmap->N, B->rmap->N));
707: break;
708: case MATPRODUCT_RARt:
709: PetscCall(MatSetSizes(D, B->rmap->n, B->rmap->n, B->rmap->N, B->rmap->N));
710: break;
711: case MATPRODUCT_PtAP:
712: PetscCall(MatSetSizes(D, B->cmap->n, B->cmap->n, B->cmap->N, B->cmap->N));
713: break;
714: default:
715: 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);
716: }
717: /* respect users who passed in a matrix for which resultname is the base type */
718: if (matmat->resultname) {
719: PetscCall(PetscObjectBaseTypeCompare((PetscObject)D, matmat->resultname, &flg));
720: if (!flg) PetscCall(MatSetType(D, matmat->resultname));
721: }
722: /* If matrix type was not set or different, we need to reset this pointers */
723: D->ops->productsymbolic = MatProductSymbolic_Shell_X;
724: D->ops->productnumeric = MatProductNumeric_Shell_X;
725: /* attach product data */
726: PetscCall(PetscNew(&mdata));
727: mdata->numeric = matmat->numeric;
728: mdata->destroy = matmat->destroy;
729: if (matmat->symbolic) PetscCall((*matmat->symbolic)(A, B, D, &mdata->ctx));
730: else PetscCall(MatSetUp(D)); /* call general setup if symbolic operation not provided */
731: PetscCheck(D->product, PetscObjectComm((PetscObject)D), PETSC_ERR_COR, "Product disappeared after user symbolic phase");
732: PetscCheck(!D->product->data, PetscObjectComm((PetscObject)D), PETSC_ERR_COR, "Product data not empty after user symbolic phase");
733: D->product->data = mdata;
734: D->product->destroy = MatProductCtxDestroy_MatMatShell;
735: /* Be sure to reset these pointers if the user did something unexpected */
736: D->ops->productsymbolic = MatProductSymbolic_Shell_X;
737: D->ops->productnumeric = MatProductNumeric_Shell_X;
738: PetscFunctionReturn(PETSC_SUCCESS);
739: }
741: static PetscErrorCode MatProductSetFromOptions_Shell_X(Mat D)
742: {
743: Mat_Product *product;
744: Mat A, B;
745: MatShellMatFunctionList matmat;
746: Mat_Shell *shell;
747: PetscBool flg;
748: char composedname[256];
750: PetscFunctionBegin;
751: MatCheckProduct(D, 1);
752: product = D->product;
753: A = product->A;
754: B = product->B;
755: PetscCall(MatIsShell(A, &flg));
756: if (!flg) PetscFunctionReturn(PETSC_SUCCESS);
757: shell = (Mat_Shell *)A->data;
758: matmat = shell->matmat;
759: PetscCall(PetscSNPrintf(composedname, sizeof(composedname), "MatProductSetFromOptions_%s_%s_C", ((PetscObject)A)->type_name, ((PetscObject)B)->type_name));
760: while (matmat) {
761: PetscCall(PetscStrcmp(composedname, matmat->composedname, &flg));
762: flg = (PetscBool)(flg && (matmat->ptype == product->type));
763: if (flg) break;
764: matmat = matmat->next;
765: }
766: if (flg) {
767: D->ops->productsymbolic = MatProductSymbolic_Shell_X;
768: } else PetscCall(PetscInfo(D, " symbolic product %s not registered for product type %s\n", composedname, MatProductTypes[product->type]));
769: PetscFunctionReturn(PETSC_SUCCESS);
770: }
772: static PetscErrorCode MatShellSetMatProductOperation_Private(Mat A, MatProductType ptype, PetscErrorCode (*symbolic)(Mat, Mat, Mat, void **), PetscErrorCode (*numeric)(Mat, Mat, Mat, void *), PetscCtxDestroyFn *destroy, char *composedname, const char *resultname)
773: {
774: PetscBool flg;
775: Mat_Shell *shell;
776: MatShellMatFunctionList matmat;
778: PetscFunctionBegin;
779: PetscCheck(numeric, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing numeric routine");
780: PetscCheck(composedname, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing composed name");
782: /* add product callback */
783: shell = (Mat_Shell *)A->data;
784: matmat = shell->matmat;
785: if (!matmat) {
786: PetscCall(PetscNew(&shell->matmat));
787: matmat = shell->matmat;
788: } else {
789: MatShellMatFunctionList entry = matmat;
790: while (entry) {
791: PetscCall(PetscStrcmp(composedname, entry->composedname, &flg));
792: flg = (PetscBool)(flg && (entry->ptype == ptype));
793: matmat = entry;
794: if (flg) goto set;
795: entry = entry->next;
796: }
797: PetscCall(PetscNew(&matmat->next));
798: matmat = matmat->next;
799: }
801: set:
802: matmat->symbolic = symbolic;
803: matmat->numeric = numeric;
804: matmat->destroy = destroy;
805: matmat->ptype = ptype;
806: PetscCall(PetscFree(matmat->composedname));
807: PetscCall(PetscFree(matmat->resultname));
808: PetscCall(PetscStrallocpy(composedname, &matmat->composedname));
809: PetscCall(PetscStrallocpy(resultname, &matmat->resultname));
810: PetscCall(PetscInfo(A, "Composing %s for product type %s with result %s\n", matmat->composedname, MatProductTypes[matmat->ptype], matmat->resultname ? matmat->resultname : "not specified"));
811: PetscCall(PetscObjectComposeFunction((PetscObject)A, matmat->composedname, MatProductSetFromOptions_Shell_X));
812: PetscFunctionReturn(PETSC_SUCCESS);
813: }
815: /*@C
816: MatShellSetMatProductOperation - Allows user to set a matrix matrix operation for a `MATSHELL` shell matrix.
818: Logically Collective; No Fortran Support
820: Input Parameters:
821: + A - the `MATSHELL` shell matrix
822: . ptype - the product type
823: . symbolic - the function for the symbolic phase (can be `NULL`)
824: . numeric - the function for the numerical phase
825: . destroy - the function for the destruction of the needed data generated during the symbolic phase (can be `NULL`)
826: . Btype - the matrix type for the matrix to be multiplied against
827: - Ctype - the matrix type for the result (can be `NULL`)
829: Level: advanced
831: Example Usage:
832: .vb
833: extern PetscErrorCode usersymbolic(Mat, Mat, Mat, void**);
834: extern PetscErrorCode usernumeric(Mat, Mat, Mat, void*);
835: PetscCtxDestroyFn *ctxdestroy;
837: MatCreateShell(comm, m, n, M, N, ctx, &A);
838: MatShellSetMatProductOperation(
839: A, MATPRODUCT_AB, usersymbolic, usernumeric, ctxdestroy,MATSEQAIJ, MATDENSE
840: );
841: // create B of type SEQAIJ etc..
842: MatProductCreate(A, B, PETSC_NULLPTR, &C);
843: MatProductSetType(C, MATPRODUCT_AB);
844: MatProductSetFromOptions(C);
845: MatProductSymbolic(C); // actually runs the user defined symbolic operation
846: MatProductNumeric(C); // actually runs the user defined numeric operation
847: // use C = A * B
848: .ve
850: Notes:
851: `MATPRODUCT_ABC` is not supported yet.
853: If the symbolic phase is not specified, `MatSetUp()` is called on the result matrix that must have its type set if `Ctype` is `NULL`.
855: Any additional data needed by the matrix product needs to be returned during the symbolic phase and destroyed with the destroy callback.
856: PETSc will take care of calling the user-defined callbacks.
857: It is allowed to specify the same callbacks for different `Btype` matrix types.
858: The couple (`Btype`,`ptype`) uniquely identifies the operation, the last specified callbacks takes precedence.
860: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellSetContext()`, `MatSetOperation()`, `MatProductType`, `MatType`, `MatSetUp()`
861: @*/
862: PetscErrorCode MatShellSetMatProductOperation(Mat A, MatProductType ptype, PetscErrorCode (*symbolic)(Mat, Mat, Mat, void **), PetscErrorCode (*numeric)(Mat, Mat, Mat, void *), PetscCtxDestroyFn *destroy, MatType Btype, MatType Ctype)
863: {
864: PetscFunctionBegin;
867: PetscCheck(ptype != MATPRODUCT_ABC, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not for product type %s", MatProductTypes[ptype]);
868: PetscCheck(numeric, PetscObjectComm((PetscObject)A), PETSC_ERR_USER, "Missing numeric routine, argument 4");
869: PetscAssertPointer(Btype, 6);
870: if (Ctype) PetscAssertPointer(Ctype, 7);
871: PetscTryMethod(A, "MatShellSetMatProductOperation_C", (Mat, MatProductType, PetscErrorCode (*)(Mat, Mat, Mat, void **), PetscErrorCode (*)(Mat, Mat, Mat, void *), PetscCtxDestroyFn *, MatType, MatType), (A, ptype, symbolic, numeric, destroy, Btype, Ctype));
872: PetscFunctionReturn(PETSC_SUCCESS);
873: }
875: static PetscErrorCode MatShellSetMatProductOperation_Shell(Mat A, MatProductType ptype, PetscErrorCode (*symbolic)(Mat, Mat, Mat, void **), PetscErrorCode (*numeric)(Mat, Mat, Mat, void *), PetscCtxDestroyFn *destroy, MatType Btype, MatType Ctype)
876: {
877: PetscBool flg;
878: char composedname[256];
879: MatRootName Bnames = MatRootNameList, Cnames = MatRootNameList;
880: PetscMPIInt size;
882: PetscFunctionBegin;
884: while (Bnames) { /* user passed in the root name */
885: PetscCall(PetscStrcmp(Btype, Bnames->rname, &flg));
886: if (flg) break;
887: Bnames = Bnames->next;
888: }
889: while (Cnames) { /* user passed in the root name */
890: PetscCall(PetscStrcmp(Ctype, Cnames->rname, &flg));
891: if (flg) break;
892: Cnames = Cnames->next;
893: }
894: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
895: Btype = Bnames ? (size > 1 ? Bnames->mname : Bnames->sname) : Btype;
896: Ctype = Cnames ? (size > 1 ? Cnames->mname : Cnames->sname) : Ctype;
897: PetscCall(PetscSNPrintf(composedname, sizeof(composedname), "MatProductSetFromOptions_%s_%s_C", ((PetscObject)A)->type_name, Btype));
898: PetscCall(MatShellSetMatProductOperation_Private(A, ptype, symbolic, numeric, destroy, composedname, Ctype));
899: PetscFunctionReturn(PETSC_SUCCESS);
900: }
902: static PetscErrorCode MatCopy_Shell(Mat A, Mat B, MatStructure str)
903: {
904: Mat_Shell *shellA = (Mat_Shell *)A->data, *shellB = (Mat_Shell *)B->data;
905: PetscBool matflg;
906: MatShellMatFunctionList matmatA;
908: PetscFunctionBegin;
909: PetscCall(MatIsShell(B, &matflg));
910: PetscCheck(matflg, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Matrix %s not derived from MATSHELL", ((PetscObject)B)->type_name);
912: B->ops[0] = A->ops[0];
913: shellB->ops[0] = shellA->ops[0];
915: if (shellA->ops->copy) PetscCall((*shellA->ops->copy)(A, B, str));
916: shellB->vscale = shellA->vscale;
917: shellB->vshift = shellA->vshift;
918: if (shellA->dshift) {
919: if (!shellB->dshift) PetscCall(VecDuplicate(shellA->dshift, &shellB->dshift));
920: PetscCall(VecCopy(shellA->dshift, shellB->dshift));
921: } else {
922: PetscCall(VecDestroy(&shellB->dshift));
923: }
924: if (shellA->left) {
925: if (!shellB->left) PetscCall(VecDuplicate(shellA->left, &shellB->left));
926: PetscCall(VecCopy(shellA->left, shellB->left));
927: } else {
928: PetscCall(VecDestroy(&shellB->left));
929: }
930: if (shellA->right) {
931: if (!shellB->right) PetscCall(VecDuplicate(shellA->right, &shellB->right));
932: PetscCall(VecCopy(shellA->right, shellB->right));
933: } else {
934: PetscCall(VecDestroy(&shellB->right));
935: }
936: PetscCall(MatDestroy(&shellB->axpy));
937: shellB->axpy_vscale = 0.0;
938: shellB->axpy_state = 0;
939: if (shellA->axpy) {
940: PetscCall(PetscObjectReference((PetscObject)shellA->axpy));
941: shellB->axpy = shellA->axpy;
942: shellB->axpy_vscale = shellA->axpy_vscale;
943: shellB->axpy_state = shellA->axpy_state;
944: }
945: if (shellA->zrows) {
946: PetscCall(ISDuplicate(shellA->zrows, &shellB->zrows));
947: if (shellA->zcols) PetscCall(ISDuplicate(shellA->zcols, &shellB->zcols));
948: PetscCall(VecDuplicate(shellA->zvals, &shellB->zvals));
949: PetscCall(VecCopy(shellA->zvals, shellB->zvals));
950: PetscCall(VecDuplicate(shellA->zvals_w, &shellB->zvals_w));
951: PetscCall(PetscObjectReference((PetscObject)shellA->zvals_sct_r));
952: PetscCall(PetscObjectReference((PetscObject)shellA->zvals_sct_c));
953: shellB->zvals_sct_r = shellA->zvals_sct_r;
954: shellB->zvals_sct_c = shellA->zvals_sct_c;
955: }
957: matmatA = shellA->matmat;
958: if (matmatA) {
959: while (matmatA->next) {
960: PetscCall(MatShellSetMatProductOperation_Private(B, matmatA->ptype, matmatA->symbolic, matmatA->numeric, matmatA->destroy, matmatA->composedname, matmatA->resultname));
961: matmatA = matmatA->next;
962: }
963: }
964: PetscFunctionReturn(PETSC_SUCCESS);
965: }
967: static PetscErrorCode MatDuplicate_Shell(Mat mat, MatDuplicateOption op, Mat *M)
968: {
969: PetscFunctionBegin;
970: PetscCall(MatCreateShell(PetscObjectComm((PetscObject)mat), mat->rmap->n, mat->cmap->n, mat->rmap->N, mat->cmap->N, NULL, M));
971: ((Mat_Shell *)(*M)->data)->ctxcontainer = ((Mat_Shell *)mat->data)->ctxcontainer;
972: PetscCall(PetscObjectCompose((PetscObject)*M, "MatShell ctx", (PetscObject)((Mat_Shell *)(*M)->data)->ctxcontainer));
973: PetscCall(PetscObjectChangeTypeName((PetscObject)*M, ((PetscObject)mat)->type_name));
974: if (op == MAT_COPY_VALUES) PetscCall(MatCopy(mat, *M, SAME_NONZERO_PATTERN));
975: PetscCall(PetscObjectCopyFortranFunctionPointers((PetscObject)mat, (PetscObject)*M));
976: PetscFunctionReturn(PETSC_SUCCESS);
977: }
979: static PetscErrorCode MatMult_Shell(Mat A, Vec x, Vec y)
980: {
981: Mat_Shell *shell = (Mat_Shell *)A->data;
982: Vec xx;
983: PetscObjectState instate, outstate;
985: PetscFunctionBegin;
986: PetscCall(MatShellPreZeroRight(A, x, &xx));
987: PetscCall(MatShellPreScaleRight(A, xx, &xx));
988: PetscCall(PetscObjectStateGet((PetscObject)y, &instate));
989: PetscCall((*shell->ops->mult)(A, xx, y));
990: PetscCall(PetscObjectStateGet((PetscObject)y, &outstate));
991: if (instate == outstate) {
992: /* increase the state of the output vector since the user did not update its state themself as should have been done */
993: PetscCall(PetscObjectStateIncrease((PetscObject)y));
994: }
995: PetscCall(MatShellShiftAndScale(A, xx, y, PETSC_FALSE));
996: PetscCall(MatShellPostScaleLeft(A, y));
997: PetscCall(MatShellPostZeroLeft(A, y));
999: if (shell->axpy) {
1000: Mat X;
1001: PetscObjectState axpy_state;
1003: PetscCall(MatShellGetContext(shell->axpy, &X));
1004: PetscCall(PetscObjectStateGet((PetscObject)X, &axpy_state));
1005: 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,...)");
1007: PetscCall(MatCreateVecs(shell->axpy, shell->axpy_right ? NULL : &shell->axpy_right, shell->axpy_left ? NULL : &shell->axpy_left));
1008: PetscCall(VecCopy(x, shell->axpy_right));
1009: PetscCall(MatMult(shell->axpy, shell->axpy_right, shell->axpy_left));
1010: PetscCall(VecAXPY(y, shell->axpy_vscale, shell->axpy_left));
1011: }
1012: PetscFunctionReturn(PETSC_SUCCESS);
1013: }
1015: static PetscErrorCode MatMultAdd_Shell(Mat A, Vec x, Vec y, Vec z)
1016: {
1017: Mat_Shell *shell = (Mat_Shell *)A->data;
1019: PetscFunctionBegin;
1020: if (y == z) {
1021: if (!shell->right_add_work) PetscCall(VecDuplicate(z, &shell->right_add_work));
1022: PetscCall(MatMult(A, x, shell->right_add_work));
1023: PetscCall(VecAXPY(z, 1.0, shell->right_add_work));
1024: } else {
1025: PetscCall(MatMult(A, x, z));
1026: PetscCall(VecAXPY(z, 1.0, y));
1027: }
1028: PetscFunctionReturn(PETSC_SUCCESS);
1029: }
1031: static PetscErrorCode MatMultTranspose_Shell(Mat A, Vec x, Vec y)
1032: {
1033: Mat_Shell *shell = (Mat_Shell *)A->data;
1034: Vec xx;
1035: PetscObjectState instate, outstate;
1037: PetscFunctionBegin;
1038: PetscCall(MatShellPreZeroLeft(A, x, &xx));
1039: PetscCall(MatShellPreScaleLeft(A, xx, &xx, PETSC_FALSE));
1040: PetscCall(PetscObjectStateGet((PetscObject)y, &instate));
1041: PetscCall((*shell->ops->multtranspose)(A, xx, y));
1042: PetscCall(PetscObjectStateGet((PetscObject)y, &outstate));
1043: if (instate == outstate) {
1044: /* increase the state of the output vector since the user did not update its state themself as should have been done */
1045: PetscCall(PetscObjectStateIncrease((PetscObject)y));
1046: }
1047: PetscCall(MatShellShiftAndScale(A, xx, y, PETSC_FALSE));
1048: PetscCall(MatShellPostScaleRight(A, y, PETSC_FALSE));
1049: PetscCall(MatShellPostZeroRight(A, y));
1051: if (shell->axpy) {
1052: Mat X;
1053: PetscObjectState axpy_state;
1055: PetscCall(MatShellGetContext(shell->axpy, &X));
1056: PetscCall(PetscObjectStateGet((PetscObject)X, &axpy_state));
1057: 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,...)");
1058: PetscCall(MatCreateVecs(shell->axpy, shell->axpy_right ? NULL : &shell->axpy_right, shell->axpy_left ? NULL : &shell->axpy_left));
1059: PetscCall(VecCopy(x, shell->axpy_left));
1060: PetscCall(MatMultTranspose(shell->axpy, shell->axpy_left, shell->axpy_right));
1061: PetscCall(VecAXPY(y, shell->axpy_vscale, shell->axpy_right));
1062: }
1063: PetscFunctionReturn(PETSC_SUCCESS);
1064: }
1066: static PetscErrorCode MatMultHermitianTranspose_Shell(Mat A, Vec x, Vec y)
1067: {
1068: Mat_Shell *shell = (Mat_Shell *)A->data;
1069: Vec xx;
1070: PetscObjectState instate, outstate;
1072: PetscFunctionBegin;
1073: PetscCall(MatShellPreZeroLeft(A, x, &xx));
1074: PetscCall(MatShellPreScaleLeft(A, xx, &xx, PETSC_TRUE));
1075: PetscCall(PetscObjectStateGet((PetscObject)y, &instate));
1076: PetscCall((*shell->ops->multhermitiantranspose)(A, xx, y));
1077: PetscCall(PetscObjectStateGet((PetscObject)y, &outstate));
1078: if (instate == outstate) {
1079: /* increase the state of the output vector since the user did not update its state themself as should have been done */
1080: PetscCall(PetscObjectStateIncrease((PetscObject)y));
1081: }
1082: PetscCall(MatShellShiftAndScale(A, xx, y, PETSC_TRUE));
1083: PetscCall(MatShellPostScaleRight(A, y, PETSC_TRUE));
1084: PetscCall(MatShellPostZeroRight(A, y));
1086: if (shell->axpy) {
1087: Mat X;
1088: PetscObjectState axpy_state;
1090: PetscCall(MatShellGetContext(shell->axpy, &X));
1091: PetscCall(PetscObjectStateGet((PetscObject)X, &axpy_state));
1092: 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,...)");
1093: PetscCall(MatCreateVecs(shell->axpy, shell->axpy_right ? NULL : &shell->axpy_right, shell->axpy_left ? NULL : &shell->axpy_left));
1094: PetscCall(VecCopy(x, shell->axpy_left));
1095: PetscCall(MatMultHermitianTranspose(shell->axpy, shell->axpy_left, shell->axpy_right));
1096: PetscCall(VecAXPY(y, PetscConj(shell->axpy_vscale), shell->axpy_right));
1097: }
1098: PetscFunctionReturn(PETSC_SUCCESS);
1099: }
1101: static PetscErrorCode MatMultTransposeAdd_Shell(Mat A, Vec x, Vec y, Vec z)
1102: {
1103: Mat_Shell *shell = (Mat_Shell *)A->data;
1105: PetscFunctionBegin;
1106: if (y == z) {
1107: if (!shell->left_add_work) PetscCall(VecDuplicate(z, &shell->left_add_work));
1108: PetscCall(MatMultTranspose(A, x, shell->left_add_work));
1109: PetscCall(VecAXPY(z, 1.0, shell->left_add_work));
1110: } else {
1111: PetscCall(MatMultTranspose(A, x, z));
1112: PetscCall(VecAXPY(z, 1.0, y));
1113: }
1114: PetscFunctionReturn(PETSC_SUCCESS);
1115: }
1117: static PetscErrorCode MatMultHermitianTransposeAdd_Shell(Mat A, Vec x, Vec y, Vec z)
1118: {
1119: Mat_Shell *shell = (Mat_Shell *)A->data;
1121: PetscFunctionBegin;
1122: if (y == z) {
1123: if (!shell->left_add_work) PetscCall(VecDuplicate(z, &shell->left_add_work));
1124: PetscCall(MatMultHermitianTranspose(A, x, shell->left_add_work));
1125: PetscCall(VecAXPY(z, 1.0, shell->left_add_work));
1126: } else {
1127: PetscCall(MatMultHermitianTranspose(A, x, z));
1128: PetscCall(VecAXPY(z, 1.0, y));
1129: }
1130: PetscFunctionReturn(PETSC_SUCCESS);
1131: }
1133: /*
1134: diag(left)(vscale*A + diag(dshift) + vshift I)diag(right)
1135: */
1136: static PetscErrorCode MatGetDiagonal_Shell(Mat A, Vec v)
1137: {
1138: Mat_Shell *shell = (Mat_Shell *)A->data;
1140: PetscFunctionBegin;
1141: PetscCheck(shell->ops->getdiagonal, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Must provide shell matrix with routine to return diagonal using MatShellSetOperation(S, MATOP_GET_DIAGONAL...)");
1142: PetscCall((*shell->ops->getdiagonal)(A, v));
1143: PetscCall(VecScale(v, shell->vscale));
1144: if (shell->dshift) PetscCall(VecAXPY(v, 1.0, shell->dshift));
1145: PetscCall(VecShift(v, shell->vshift));
1146: if (shell->left) PetscCall(VecPointwiseMult(v, v, shell->left));
1147: if (shell->right) PetscCall(VecPointwiseMult(v, v, shell->right));
1148: if (shell->zrows) {
1149: PetscCall(VecScatterBegin(shell->zvals_sct_r, shell->zvals, v, INSERT_VALUES, SCATTER_REVERSE));
1150: PetscCall(VecScatterEnd(shell->zvals_sct_r, shell->zvals, v, INSERT_VALUES, SCATTER_REVERSE));
1151: }
1152: if (shell->axpy) {
1153: Mat X;
1154: PetscObjectState axpy_state;
1156: PetscCall(MatShellGetContext(shell->axpy, &X));
1157: PetscCall(PetscObjectStateGet((PetscObject)X, &axpy_state));
1158: 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,...)");
1159: PetscCall(MatCreateVecs(shell->axpy, NULL, shell->axpy_left ? NULL : &shell->axpy_left));
1160: PetscCall(MatGetDiagonal(shell->axpy, shell->axpy_left));
1161: PetscCall(VecAXPY(v, shell->axpy_vscale, shell->axpy_left));
1162: }
1163: PetscFunctionReturn(PETSC_SUCCESS);
1164: }
1166: static PetscErrorCode MatGetDiagonalBlock_Shell(Mat A, Mat *b)
1167: {
1168: Mat_Shell *shell = (Mat_Shell *)A->data;
1169: Vec left = NULL, right = NULL;
1171: PetscFunctionBegin;
1172: 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
1173: PetscCheck(!shell->axpy, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Cannot call MatGetDiagonalBlock() if MatAXPY() has been called on the input Mat"); // TODO FIXME
1174: PetscCheck(!shell->dshift, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Cannot call MatGetDiagonalBlock() if MatDiagonalSet() has been called on the input Mat"); // TODO FIXME
1175: PetscCheck(shell->ops->getdiagonalblock, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Must provide shell matrix with routine to return diagonal block using MatShellSetOperation(S, MATOP_GET_DIAGONAL_BLOCK...)");
1176: PetscCall((*shell->ops->getdiagonalblock)(A, b));
1177: PetscCall(MatScale(*b, shell->vscale));
1178: PetscCall(MatShift(*b, shell->vshift));
1179: if (shell->left) {
1180: PetscCall(VecCreateLocalVector(shell->left, &left));
1181: PetscCall(VecGetLocalVectorRead(shell->left, left));
1182: }
1183: if (shell->right) {
1184: PetscCall(VecCreateLocalVector(shell->right, &right));
1185: PetscCall(VecGetLocalVectorRead(shell->right, right));
1186: }
1187: PetscCall(MatDiagonalScale(*b, left, right));
1188: if (shell->left) {
1189: PetscCall(VecRestoreLocalVectorRead(shell->left, left));
1190: PetscCall(VecDestroy(&left));
1191: }
1192: if (shell->right) {
1193: PetscCall(VecRestoreLocalVectorRead(shell->right, right));
1194: PetscCall(VecDestroy(&right));
1195: }
1196: PetscFunctionReturn(PETSC_SUCCESS);
1197: }
1199: static PetscErrorCode MatShift_Shell(Mat Y, PetscScalar a)
1200: {
1201: Mat_Shell *shell = (Mat_Shell *)Y->data;
1202: PetscBool flg;
1204: PetscFunctionBegin;
1205: PetscCall(MatHasCongruentLayouts(Y, &flg));
1206: PetscCheck(flg, PetscObjectComm((PetscObject)Y), PETSC_ERR_SUP, "Cannot shift shell matrix if it is not congruent");
1207: if (shell->left || shell->right) {
1208: if (!shell->dshift) {
1209: PetscCall(VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift));
1210: PetscCall(VecSet(shell->dshift, a));
1211: } else {
1212: if (shell->left) PetscCall(VecPointwiseMult(shell->dshift, shell->dshift, shell->left));
1213: if (shell->right) PetscCall(VecPointwiseMult(shell->dshift, shell->dshift, shell->right));
1214: PetscCall(VecShift(shell->dshift, a));
1215: }
1216: if (shell->left) PetscCall(VecPointwiseDivide(shell->dshift, shell->dshift, shell->left));
1217: if (shell->right) PetscCall(VecPointwiseDivide(shell->dshift, shell->dshift, shell->right));
1218: } else shell->vshift += a;
1219: if (shell->zrows) PetscCall(VecShift(shell->zvals, a));
1220: PetscFunctionReturn(PETSC_SUCCESS);
1221: }
1223: static PetscErrorCode MatDiagonalSet_Shell_Private(Mat A, Vec D, PetscScalar s)
1224: {
1225: Mat_Shell *shell = (Mat_Shell *)A->data;
1227: PetscFunctionBegin;
1228: if (!shell->dshift) PetscCall(VecDuplicate(D, &shell->dshift));
1229: if (shell->left || shell->right) {
1230: if (!shell->right_work) PetscCall(VecDuplicate(shell->left ? shell->left : shell->right, &shell->right_work));
1231: if (shell->left && shell->right) {
1232: PetscCall(VecPointwiseDivide(shell->right_work, D, shell->left));
1233: PetscCall(VecPointwiseDivide(shell->right_work, shell->right_work, shell->right));
1234: } else if (shell->left) {
1235: PetscCall(VecPointwiseDivide(shell->right_work, D, shell->left));
1236: } else {
1237: PetscCall(VecPointwiseDivide(shell->right_work, D, shell->right));
1238: }
1239: PetscCall(VecAXPY(shell->dshift, s, shell->right_work));
1240: } else {
1241: PetscCall(VecAXPY(shell->dshift, s, D));
1242: }
1243: PetscFunctionReturn(PETSC_SUCCESS);
1244: }
1246: static PetscErrorCode MatDiagonalSet_Shell(Mat A, Vec D, InsertMode ins)
1247: {
1248: Mat_Shell *shell = (Mat_Shell *)A->data;
1249: Vec d;
1250: PetscBool flg;
1252: PetscFunctionBegin;
1253: PetscCall(MatHasCongruentLayouts(A, &flg));
1254: PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Cannot diagonal set or shift shell matrix if it is not congruent");
1255: if (ins == INSERT_VALUES) {
1256: PetscCall(VecDuplicate(D, &d));
1257: PetscCall(MatGetDiagonal(A, d));
1258: PetscCall(MatDiagonalSet_Shell_Private(A, d, -1.));
1259: PetscCall(MatDiagonalSet_Shell_Private(A, D, 1.));
1260: PetscCall(VecDestroy(&d));
1261: if (shell->zrows) PetscCall(VecCopy(D, shell->zvals));
1262: } else {
1263: PetscCall(MatDiagonalSet_Shell_Private(A, D, 1.));
1264: if (shell->zrows) PetscCall(VecAXPY(shell->zvals, 1.0, D));
1265: }
1266: PetscFunctionReturn(PETSC_SUCCESS);
1267: }
1269: static PetscErrorCode MatScale_Shell(Mat Y, PetscScalar a)
1270: {
1271: Mat_Shell *shell = (Mat_Shell *)Y->data;
1273: PetscFunctionBegin;
1274: shell->vscale *= a;
1275: shell->vshift *= a;
1276: if (shell->dshift) PetscCall(VecScale(shell->dshift, a));
1277: shell->axpy_vscale *= a;
1278: if (shell->zrows) PetscCall(VecScale(shell->zvals, a));
1279: PetscFunctionReturn(PETSC_SUCCESS);
1280: }
1282: static PetscErrorCode MatDiagonalScale_Shell(Mat Y, Vec left, Vec right)
1283: {
1284: Mat_Shell *shell = (Mat_Shell *)Y->data;
1286: PetscFunctionBegin;
1287: if (left) {
1288: if (!shell->left) {
1289: PetscCall(VecDuplicate(left, &shell->left));
1290: PetscCall(VecCopy(left, shell->left));
1291: } else {
1292: PetscCall(VecPointwiseMult(shell->left, shell->left, left));
1293: }
1294: if (shell->zrows) PetscCall(VecPointwiseMult(shell->zvals, shell->zvals, left));
1295: }
1296: if (right) {
1297: if (!shell->right) {
1298: PetscCall(VecDuplicate(right, &shell->right));
1299: PetscCall(VecCopy(right, shell->right));
1300: } else {
1301: PetscCall(VecPointwiseMult(shell->right, shell->right, right));
1302: }
1303: if (shell->zrows) {
1304: if (!shell->left_work) PetscCall(MatCreateVecs(Y, NULL, &shell->left_work));
1305: PetscCall(VecSet(shell->zvals_w, 1.0));
1306: PetscCall(VecScatterBegin(shell->zvals_sct_c, right, shell->zvals_w, INSERT_VALUES, SCATTER_FORWARD));
1307: PetscCall(VecScatterEnd(shell->zvals_sct_c, right, shell->zvals_w, INSERT_VALUES, SCATTER_FORWARD));
1308: PetscCall(VecPointwiseMult(shell->zvals, shell->zvals, shell->zvals_w));
1309: }
1310: }
1311: if (shell->axpy) PetscCall(MatDiagonalScale(shell->axpy, left, right));
1312: PetscFunctionReturn(PETSC_SUCCESS);
1313: }
1315: PETSC_INTERN PetscErrorCode MatAssemblyEnd_Shell(Mat Y, MatAssemblyType t)
1316: {
1317: Mat_Shell *shell = (Mat_Shell *)Y->data;
1319: PetscFunctionBegin;
1320: if (t == MAT_FINAL_ASSEMBLY) {
1321: shell->vshift = 0.0;
1322: shell->vscale = 1.0;
1323: shell->axpy_vscale = 0.0;
1324: shell->axpy_state = 0;
1325: PetscCall(VecDestroy(&shell->dshift));
1326: PetscCall(VecDestroy(&shell->left));
1327: PetscCall(VecDestroy(&shell->right));
1328: PetscCall(MatDestroy(&shell->axpy));
1329: PetscCall(VecDestroy(&shell->axpy_left));
1330: PetscCall(VecDestroy(&shell->axpy_right));
1331: PetscCall(VecScatterDestroy(&shell->zvals_sct_c));
1332: PetscCall(VecScatterDestroy(&shell->zvals_sct_r));
1333: PetscCall(ISDestroy(&shell->zrows));
1334: PetscCall(ISDestroy(&shell->zcols));
1335: }
1336: PetscFunctionReturn(PETSC_SUCCESS);
1337: }
1339: static PetscErrorCode MatAXPY_Shell(Mat Y, PetscScalar a, Mat X, MatStructure str)
1340: {
1341: Mat_Shell *shell = (Mat_Shell *)Y->data;
1343: PetscFunctionBegin;
1344: if (X == Y) {
1345: PetscCall(MatScale(Y, 1.0 + a));
1346: PetscFunctionReturn(PETSC_SUCCESS);
1347: }
1348: if (!shell->axpy) {
1349: PetscCall(MatConvertFrom_Shell(X, MATSHELL, MAT_INITIAL_MATRIX, &shell->axpy));
1350: shell->axpy_vscale = a;
1351: PetscCall(PetscObjectStateGet((PetscObject)X, &shell->axpy_state));
1352: } else {
1353: PetscCall(MatAXPY(shell->axpy, a / shell->axpy_vscale, X, str));
1354: }
1355: PetscFunctionReturn(PETSC_SUCCESS);
1356: }
1358: static struct _MatOps MatOps_Values = {NULL,
1359: NULL,
1360: NULL,
1361: NULL,
1362: /* 4*/ MatMultAdd_Shell,
1363: NULL,
1364: MatMultTransposeAdd_Shell,
1365: NULL,
1366: NULL,
1367: NULL,
1368: /*10*/ NULL,
1369: NULL,
1370: NULL,
1371: NULL,
1372: NULL,
1373: /*15*/ NULL,
1374: NULL,
1375: NULL,
1376: MatDiagonalScale_Shell,
1377: NULL,
1378: /*20*/ NULL,
1379: MatAssemblyEnd_Shell,
1380: NULL,
1381: NULL,
1382: /*24*/ MatZeroRows_Shell,
1383: NULL,
1384: NULL,
1385: NULL,
1386: NULL,
1387: /*29*/ NULL,
1388: NULL,
1389: NULL,
1390: /*32*/ NULL,
1391: NULL,
1392: /*34*/ MatDuplicate_Shell,
1393: NULL,
1394: NULL,
1395: NULL,
1396: NULL,
1397: /*39*/ MatAXPY_Shell,
1398: NULL,
1399: NULL,
1400: NULL,
1401: MatCopy_Shell,
1402: /*44*/ NULL,
1403: MatScale_Shell,
1404: MatShift_Shell,
1405: MatDiagonalSet_Shell,
1406: MatZeroRowsColumns_Shell,
1407: /*49*/ NULL,
1408: NULL,
1409: NULL,
1410: NULL,
1411: NULL,
1412: /*54*/ NULL,
1413: NULL,
1414: NULL,
1415: NULL,
1416: NULL,
1417: /*59*/ NULL,
1418: MatDestroy_Shell,
1419: NULL,
1420: MatConvertFrom_Shell,
1421: NULL,
1422: /*64*/ NULL,
1423: NULL,
1424: NULL,
1425: NULL,
1426: NULL,
1427: /*69*/ NULL,
1428: MatConvert_Shell,
1429: NULL,
1430: NULL,
1431: NULL,
1432: /*74*/ NULL,
1433: NULL,
1434: NULL,
1435: NULL,
1436: NULL,
1437: /*79*/ NULL,
1438: NULL,
1439: NULL,
1440: NULL,
1441: NULL,
1442: /*84*/ NULL,
1443: NULL,
1444: NULL,
1445: NULL,
1446: NULL,
1447: /*89*/ NULL,
1448: NULL,
1449: NULL,
1450: NULL,
1451: NULL,
1452: /*94*/ NULL,
1453: NULL,
1454: NULL,
1455: NULL,
1456: NULL,
1457: /*99*/ NULL,
1458: NULL,
1459: NULL,
1460: NULL,
1461: NULL,
1462: /*104*/ NULL,
1463: NULL,
1464: NULL,
1465: NULL,
1466: NULL,
1467: /*109*/ NULL,
1468: NULL,
1469: NULL,
1470: MatMultHermitianTransposeAdd_Shell,
1471: NULL,
1472: /*114*/ NULL,
1473: NULL,
1474: NULL,
1475: NULL,
1476: NULL,
1477: /*119*/ NULL,
1478: NULL,
1479: NULL,
1480: NULL,
1481: NULL,
1482: /*124*/ NULL,
1483: NULL,
1484: NULL,
1485: NULL,
1486: NULL,
1487: /*129*/ NULL,
1488: NULL,
1489: NULL,
1490: NULL,
1491: NULL,
1492: /*134*/ NULL,
1493: NULL,
1494: NULL,
1495: NULL,
1496: NULL,
1497: /*139*/ NULL,
1498: NULL,
1499: NULL,
1500: NULL,
1501: NULL,
1502: /*144*/ MatADot_Default,
1503: MatANorm_Default,
1504: NULL,
1505: NULL};
1507: static PetscErrorCode MatShellSetContext_Shell(Mat mat, PetscCtx ctx)
1508: {
1509: Mat_Shell *shell = (Mat_Shell *)mat->data;
1511: PetscFunctionBegin;
1512: if (ctx) {
1513: PetscContainer ctxcontainer;
1514: PetscCall(PetscContainerCreate(PetscObjectComm((PetscObject)mat), &ctxcontainer));
1515: PetscCall(PetscContainerSetPointer(ctxcontainer, ctx));
1516: PetscCall(PetscObjectCompose((PetscObject)mat, "MatShell ctx", (PetscObject)ctxcontainer));
1517: shell->ctxcontainer = ctxcontainer;
1518: PetscCall(PetscContainerDestroy(&ctxcontainer));
1519: } else {
1520: PetscCall(PetscObjectCompose((PetscObject)mat, "MatShell ctx", NULL));
1521: shell->ctxcontainer = NULL;
1522: }
1523: PetscFunctionReturn(PETSC_SUCCESS);
1524: }
1526: static PetscErrorCode MatShellSetContextDestroy_Shell(Mat mat, PetscCtxDestroyFn *f)
1527: {
1528: Mat_Shell *shell = (Mat_Shell *)mat->data;
1530: PetscFunctionBegin;
1531: if (shell->ctxcontainer) PetscCall(PetscContainerSetCtxDestroy(shell->ctxcontainer, f));
1532: PetscFunctionReturn(PETSC_SUCCESS);
1533: }
1535: PetscErrorCode MatShellSetContext_Immutable(Mat mat, PetscCtx ctx)
1536: {
1537: PetscFunctionBegin;
1538: 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);
1539: PetscFunctionReturn(PETSC_SUCCESS);
1540: }
1542: PetscErrorCode MatShellSetContextDestroy_Immutable(Mat mat, PetscCtxDestroyFn *f)
1543: {
1544: PetscFunctionBegin;
1545: 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);
1546: PetscFunctionReturn(PETSC_SUCCESS);
1547: }
1549: PetscErrorCode MatShellSetManageScalingShifts_Immutable(Mat mat)
1550: {
1551: PetscFunctionBegin;
1552: 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);
1553: PetscFunctionReturn(PETSC_SUCCESS);
1554: }
1556: static PetscErrorCode MatShellSetVecType_Shell(Mat mat, VecType vtype)
1557: {
1558: PetscFunctionBegin;
1559: PetscCall(PetscFree(mat->defaultvectype));
1560: PetscCall(PetscStrallocpy(vtype, (char **)&mat->defaultvectype));
1561: PetscFunctionReturn(PETSC_SUCCESS);
1562: }
1564: static PetscErrorCode MatShellSetManageScalingShifts_Shell(Mat A)
1565: {
1566: Mat_Shell *shell = (Mat_Shell *)A->data;
1568: PetscFunctionBegin;
1569: shell->managescalingshifts = PETSC_FALSE;
1570: A->ops->diagonalset = NULL;
1571: A->ops->diagonalscale = NULL;
1572: A->ops->scale = NULL;
1573: A->ops->shift = NULL;
1574: A->ops->axpy = NULL;
1575: PetscFunctionReturn(PETSC_SUCCESS);
1576: }
1578: static PetscErrorCode MatShellGetScalingShifts_Shell(Mat A, PetscScalar *vshift, PetscScalar *vscale, Vec *dshift, Vec *left, Vec *right, Mat *axpy, IS *zrows, IS *zcols)
1579: {
1580: Mat_Shell *shell = (Mat_Shell *)A->data;
1582: PetscFunctionBegin;
1583: PetscCheck(shell->managescalingshifts, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Cannot get scaling and shifts if MatShellSetManageScalingShifts() has been called");
1584: 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()");
1585: else if (vshift) *vshift = shell->vshift;
1586: 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()");
1587: else if (vscale) *vscale = shell->vscale;
1588: 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()");
1589: else if (dshift) *dshift = shell->dshift;
1590: 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()");
1591: else if (left) *left = shell->left;
1592: 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()");
1593: else if (right) *right = shell->right;
1594: 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()");
1595: else if (axpy) *axpy = shell->axpy;
1596: 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()");
1597: else if (zrows) *zrows = shell->zrows;
1598: 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()");
1599: else if (zcols) *zcols = shell->zcols;
1600: PetscFunctionReturn(PETSC_SUCCESS);
1601: }
1603: static PetscErrorCode MatShellSetOperation_Shell(Mat mat, MatOperation op, PetscErrorCodeFn *f)
1604: {
1605: Mat_Shell *shell = (Mat_Shell *)mat->data;
1607: PetscFunctionBegin;
1608: switch (op) {
1609: case MATOP_DESTROY:
1610: shell->ops->destroy = (PetscErrorCode (*)(Mat))f;
1611: break;
1612: case MATOP_VIEW:
1613: if (!mat->ops->viewnative) mat->ops->viewnative = mat->ops->view;
1614: mat->ops->view = (PetscErrorCode (*)(Mat, PetscViewer))f;
1615: break;
1616: case MATOP_COPY:
1617: shell->ops->copy = (PetscErrorCode (*)(Mat, Mat, MatStructure))f;
1618: break;
1619: case MATOP_DIAGONAL_SET:
1620: case MATOP_DIAGONAL_SCALE:
1621: case MATOP_SHIFT:
1622: case MATOP_SCALE:
1623: case MATOP_AXPY:
1624: case MATOP_ZERO_ROWS:
1625: case MATOP_ZERO_ROWS_COLUMNS:
1626: case MATOP_ZERO_ROWS_LOCAL:
1627: case MATOP_ZERO_ROWS_COLUMNS_LOCAL:
1628: PetscCheck(!shell->managescalingshifts, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "MATSHELL is managing scalings and shifts, see MatShellSetManageScalingShifts()");
1629: (((PetscErrorCodeFn **)mat->ops)[op]) = f;
1630: break;
1631: case MATOP_GET_DIAGONAL:
1632: if (shell->managescalingshifts) {
1633: shell->ops->getdiagonal = (PetscErrorCode (*)(Mat, Vec))f;
1634: mat->ops->getdiagonal = MatGetDiagonal_Shell;
1635: } else {
1636: shell->ops->getdiagonal = NULL;
1637: mat->ops->getdiagonal = (PetscErrorCode (*)(Mat, Vec))f;
1638: }
1639: break;
1640: case MATOP_GET_DIAGONAL_BLOCK:
1641: if (shell->managescalingshifts) {
1642: shell->ops->getdiagonalblock = (PetscErrorCode (*)(Mat, Mat *))f;
1643: mat->ops->getdiagonalblock = MatGetDiagonalBlock_Shell;
1644: } else {
1645: shell->ops->getdiagonalblock = NULL;
1646: mat->ops->getdiagonalblock = (PetscErrorCode (*)(Mat, Mat *))f;
1647: }
1648: break;
1649: case MATOP_MULT:
1650: if (shell->managescalingshifts) {
1651: shell->ops->mult = (PetscErrorCode (*)(Mat, Vec, Vec))f;
1652: mat->ops->mult = MatMult_Shell;
1653: } else {
1654: shell->ops->mult = NULL;
1655: mat->ops->mult = (PetscErrorCode (*)(Mat, Vec, Vec))f;
1656: }
1657: break;
1658: case MATOP_MULT_TRANSPOSE:
1659: if (shell->managescalingshifts) {
1660: shell->ops->multtranspose = (PetscErrorCode (*)(Mat, Vec, Vec))f;
1661: mat->ops->multtranspose = MatMultTranspose_Shell;
1662: } else {
1663: shell->ops->multtranspose = NULL;
1664: mat->ops->multtranspose = (PetscErrorCode (*)(Mat, Vec, Vec))f;
1665: }
1666: break;
1667: case MATOP_MULT_HERMITIAN_TRANSPOSE:
1668: if (shell->managescalingshifts) {
1669: shell->ops->multhermitiantranspose = (PetscErrorCode (*)(Mat, Vec, Vec))f;
1670: mat->ops->multhermitiantranspose = MatMultHermitianTranspose_Shell;
1671: } else {
1672: shell->ops->multhermitiantranspose = NULL;
1673: mat->ops->multhermitiantranspose = (PetscErrorCode (*)(Mat, Vec, Vec))f;
1674: }
1675: break;
1676: default:
1677: (((PetscErrorCodeFn **)mat->ops)[op]) = f;
1678: break;
1679: }
1680: PetscFunctionReturn(PETSC_SUCCESS);
1681: }
1683: static PetscErrorCode MatShellGetOperation_Shell(Mat mat, MatOperation op, PetscErrorCodeFn **f)
1684: {
1685: Mat_Shell *shell = (Mat_Shell *)mat->data;
1687: PetscFunctionBegin;
1688: switch (op) {
1689: case MATOP_DESTROY:
1690: *f = (PetscErrorCodeFn *)shell->ops->destroy;
1691: break;
1692: case MATOP_VIEW:
1693: *f = (PetscErrorCodeFn *)mat->ops->view;
1694: break;
1695: case MATOP_COPY:
1696: *f = (PetscErrorCodeFn *)shell->ops->copy;
1697: break;
1698: case MATOP_DIAGONAL_SET:
1699: case MATOP_DIAGONAL_SCALE:
1700: case MATOP_SHIFT:
1701: case MATOP_SCALE:
1702: case MATOP_AXPY:
1703: case MATOP_ZERO_ROWS:
1704: case MATOP_ZERO_ROWS_COLUMNS:
1705: *f = (((PetscErrorCodeFn **)mat->ops)[op]);
1706: break;
1707: case MATOP_GET_DIAGONAL:
1708: if (shell->ops->getdiagonal) *f = (PetscErrorCodeFn *)shell->ops->getdiagonal;
1709: else *f = (((PetscErrorCodeFn **)mat->ops)[op]);
1710: break;
1711: case MATOP_GET_DIAGONAL_BLOCK:
1712: if (shell->ops->getdiagonalblock) *f = (PetscErrorCodeFn *)shell->ops->getdiagonalblock;
1713: else *f = (((PetscErrorCodeFn **)mat->ops)[op]);
1714: break;
1715: case MATOP_MULT:
1716: if (shell->ops->mult) *f = (PetscErrorCodeFn *)shell->ops->mult;
1717: else *f = (((PetscErrorCodeFn **)mat->ops)[op]);
1718: break;
1719: case MATOP_MULT_TRANSPOSE:
1720: if (shell->ops->multtranspose) *f = (PetscErrorCodeFn *)shell->ops->multtranspose;
1721: else *f = (((PetscErrorCodeFn **)mat->ops)[op]);
1722: break;
1723: case MATOP_MULT_HERMITIAN_TRANSPOSE:
1724: if (shell->ops->multhermitiantranspose) *f = (PetscErrorCodeFn *)shell->ops->multhermitiantranspose;
1725: else *f = (((PetscErrorCodeFn **)mat->ops)[op]);
1726: break;
1727: default:
1728: *f = (((PetscErrorCodeFn **)mat->ops)[op]);
1729: }
1730: PetscFunctionReturn(PETSC_SUCCESS);
1731: }
1733: /*MC
1734: MATSHELL - MATSHELL = "shell" - A matrix type to be used to define your own matrix type with its own data structure -- perhaps matrix-free.
1736: Level: advanced
1738: Notes:
1739: See `MatCreateShell()` for details on the usage of `MATSHELL`
1741: `PCSHELL` can be used in conjunction with `MATSHELL` to provide a custom preconditioner appropriate for your `MATSHELL`. Since
1742: many standard preconditioners such as `PCILU` depend on having an explicit representation of the matrix entries they cannot be used
1743: directly with `MATSHELL`.
1745: .seealso: [](ch_matrices), `Mat`, `MatCreateShell()`, `PCSHELL`
1746: M*/
1748: PETSC_EXTERN PetscErrorCode MatCreate_Shell(Mat A)
1749: {
1750: Mat_Shell *b;
1752: PetscFunctionBegin;
1753: PetscCall(PetscNew(&b));
1754: A->data = (void *)b;
1755: A->ops[0] = MatOps_Values;
1757: b->ctxcontainer = NULL;
1758: b->vshift = 0.0;
1759: b->vscale = 1.0;
1760: b->managescalingshifts = PETSC_TRUE;
1761: A->assembled = PETSC_TRUE;
1762: A->preallocated = PETSC_FALSE;
1764: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellGetContext_C", MatShellGetContext_Shell));
1765: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetContext_C", MatShellSetContext_Shell));
1766: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetContextDestroy_C", MatShellSetContextDestroy_Shell));
1767: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetVecType_C", MatShellSetVecType_Shell));
1768: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetManageScalingShifts_C", MatShellSetManageScalingShifts_Shell));
1769: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellGetScalingShifts_C", MatShellGetScalingShifts_Shell));
1770: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetOperation_C", MatShellSetOperation_Shell));
1771: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellGetOperation_C", MatShellGetOperation_Shell));
1772: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetMatProductOperation_C", MatShellSetMatProductOperation_Shell));
1773: PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATSHELL));
1774: PetscFunctionReturn(PETSC_SUCCESS);
1775: }
1777: /*@C
1778: MatCreateShell - Creates a new matrix of `MatType` `MATSHELL` for use with a user-defined
1779: private matrix data storage format.
1781: Collective
1783: Input Parameters:
1784: + comm - MPI communicator
1785: . m - number of local rows (or `PETSC_DECIDE` to have calculated if `M` is given)
1786: . n - number of local columns (or `PETSC_DECIDE` to have calculated if `N` is given)
1787: . M - number of global rows (may be `PETSC_DETERMINE` to have calculated if `m` is given)
1788: . N - number of global columns (may be `PETSC_DETERMINE` to have calculated if `n` is given)
1789: - ctx - pointer to data needed by the shell matrix routines
1791: Output Parameter:
1792: . A - the matrix
1794: Level: advanced
1796: Example Usage:
1797: .vb
1798: extern PetscErrorCode mult(Mat, Vec, Vec);
1800: MatCreateShell(comm, m, n, M, N, ctx, &mat);
1801: MatShellSetOperation(mat, MATOP_MULT, (void(*)(void))mult);
1802: MatShellSetContext(mat,ctx);
1803: // Use matrix for operations that have been set
1804: MatDestroy(mat);
1805: .ve
1807: Notes:
1808: The shell matrix type is intended to provide a simple way for users to write a custom matrix specifically for their application.
1810: `MatCreateShell()` is used in conjunction with `MatShellSetContext()` and `MatShellSetOperation()`.
1812: PETSc requires that matrices and vectors being used for certain
1813: operations are partitioned accordingly. For example, when
1814: creating a shell matrix, `A`, that supports parallel matrix-vector
1815: products using `MatMult`(A,x,y) the user should set the number
1816: of local matrix rows to be the number of local elements of the
1817: corresponding result vector, y. Note that this is information is
1818: required for use of the matrix interface routines, even though
1819: the shell matrix may not actually be physically partitioned.
1820: For example,
1822: .vb
1823: Vec x, y
1824: extern PetscErrorCode mult(Mat,Vec,Vec);
1825: Mat A
1827: VecCreateMPI(comm,PETSC_DECIDE,M,&y);
1828: VecCreateMPI(comm,PETSC_DECIDE,N,&x);
1829: VecGetLocalSize(y,&m);
1830: VecGetLocalSize(x,&n);
1831: MatCreateShell(comm,m,n,M,N,ctx,&A);
1832: MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
1833: MatMult(A,x,y);
1834: MatDestroy(&A);
1835: VecDestroy(&y);
1836: VecDestroy(&x);
1837: .ve
1839: `MATSHELL` handles `MatShift()`, `MatDiagonalSet()`, `MatDiagonalScale()`, `MatAXPY()`, `MatScale()`, `MatZeroRows()` and `MatZeroRowsColumns()`
1840: internally, so these operations cannot be overwritten unless `MatShellSetManageScalingShifts()` is called.
1842: Developer Notes:
1843: For rectangular matrices do all the scalings and shifts make sense?
1845: Regarding shifting and scaling. The general form is
1847: diag(left)(vscale*A + diag(dshift) + vshift I)diag(right)
1849: The order you apply the operations is important. For example if you have a dshift then
1850: apply a MatScale(s) you get s*vscale*A + s*diag(shift). But if you first scale and then shift
1851: you get s*vscale*A + diag(shift)
1853: A is the user provided function.
1855: `KSP`/`PC` uses changes in the `Mat`'s "state" to decide if preconditioners need to be rebuilt `PCSetUp()` only calls the setup() for
1856: for the `PC` implementation if the `Mat` state has increased from the previous call. Thus to get changes in a `MATSHELL` to trigger
1857: an update in the preconditioner you must call `MatAssemblyBegin()` and `MatAssemblyEnd()` or `PetscObjectStateIncrease`((`PetscObject`)mat);
1858: each time the `MATSHELL` matrix has changed.
1860: Matrix-matrix product operations can be specified using `MatShellSetMatProductOperation()`
1862: Calling `MatAssemblyBegin()`/`MatAssemblyEnd()` on a `MATSHELL` removes any previously supplied shift and scales that were provided
1863: with `MatDiagonalSet()`, `MatShift()`, `MatScale()`, or `MatDiagonalScale()`.
1865: Fortran Notes:
1866: To use this from Fortran with a `ctx` you must write an interface definition for this
1867: function and for `MatShellGetContext()` that tells Fortran the Fortran derived data type you are passing
1868: in as the `ctx` argument.
1870: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatShellSetOperation()`, `MatHasOperation()`, `MatShellGetContext()`, `MatShellSetContext()`, `MatShellSetManageScalingShifts()`, `MatShellSetMatProductOperation()`
1871: @*/
1872: PetscErrorCode MatCreateShell(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, PetscInt N, PetscCtx ctx, Mat *A)
1873: {
1874: PetscFunctionBegin;
1875: PetscCall(MatCreate(comm, A));
1876: PetscCall(MatSetSizes(*A, m, n, M, N));
1877: PetscCall(MatSetType(*A, MATSHELL));
1878: PetscCall(MatShellSetContext(*A, ctx));
1879: PetscCall(MatSetUp(*A));
1880: PetscFunctionReturn(PETSC_SUCCESS);
1881: }
1883: /*@
1884: MatShellSetContext - sets the context for a `MATSHELL` shell matrix
1886: Logically Collective
1888: Input Parameters:
1889: + mat - the `MATSHELL` shell matrix
1890: - ctx - the context
1892: Level: advanced
1894: Note:
1895: This provides an easy way, along with `MatCreateShell()` and `MatShellSetOperation()` to provide a custom matrix format
1896: specifically for your application.
1898: Fortran Notes:
1899: You must write a Fortran interface definition for this
1900: function that tells Fortran the Fortran derived data type that you are passing in as the `ctx` argument.
1902: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`
1903: @*/
1904: PetscErrorCode MatShellSetContext(Mat mat, PetscCtx ctx)
1905: {
1906: PetscFunctionBegin;
1908: PetscTryMethod(mat, "MatShellSetContext_C", (Mat, void *), (mat, ctx));
1909: PetscFunctionReturn(PETSC_SUCCESS);
1910: }
1912: /*@C
1913: MatShellSetContextDestroy - sets the destroy function for a `MATSHELL` shell matrix context
1915: Logically Collective
1917: Input Parameters:
1918: + mat - the shell matrix
1919: - f - the context destroy function, see `PetscCtxDestroyFn` for calling sequence
1921: Level: advanced
1923: Note:
1924: If the `MatShell` is never duplicated, the behavior of this function is equivalent
1925: to `MatShellSetOperation`(`Mat`,`MATOP_DESTROY`,f). However, `MatShellSetContextDestroy()`
1926: ensures proper reference counting for the user provided context data in the case that
1927: the `MATSHELL` is duplicated.
1929: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellSetContext()`,
1930: `PetscCtxDestroyFn`
1931: @*/
1932: PetscErrorCode MatShellSetContextDestroy(Mat mat, PetscCtxDestroyFn *f)
1933: {
1934: PetscFunctionBegin;
1936: PetscTryMethod(mat, "MatShellSetContextDestroy_C", (Mat, PetscCtxDestroyFn *), (mat, f));
1937: PetscFunctionReturn(PETSC_SUCCESS);
1938: }
1940: /*@C
1941: MatShellSetVecType - Sets the `VecType` of `Vec` returned by `MatCreateVecs()`
1943: Logically Collective
1945: Input Parameters:
1946: + mat - the `MATSHELL` shell matrix
1947: - vtype - type to use for creating vectors
1949: Level: advanced
1951: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateVecs()`
1952: @*/
1953: PetscErrorCode MatShellSetVecType(Mat mat, VecType vtype)
1954: {
1955: PetscFunctionBegin;
1956: PetscTryMethod(mat, "MatShellSetVecType_C", (Mat, VecType), (mat, vtype));
1957: PetscFunctionReturn(PETSC_SUCCESS);
1958: }
1960: /*@
1961: MatShellSetManageScalingShifts - Allows the user to control the scaling and shift operations of the `MATSHELL`. Must be called immediately
1962: after `MatCreateShell()`
1964: Logically Collective
1966: Input Parameter:
1967: . A - the `MATSHELL` shell matrix
1969: Level: advanced
1971: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellSetContext()`, `MatShellSetOperation()`
1972: @*/
1973: PetscErrorCode MatShellSetManageScalingShifts(Mat A)
1974: {
1975: PetscFunctionBegin;
1977: PetscTryMethod(A, "MatShellSetManageScalingShifts_C", (Mat), (A));
1978: PetscFunctionReturn(PETSC_SUCCESS);
1979: }
1981: /*@C
1982: MatShellGetScalingShifts - Gets members of a `MATSHELL` used internally for scaling and
1983: shifting the `Mat` or calling `MatAXPY()`, `MatZeroRows()`, or `MatZeroRowsColumns()` with it
1985: Logically Collective
1987: Input Parameter:
1988: . A - the `MATSHELL` shell matrix
1990: Output Parameters:
1991: + vshift - `PetscScalar` pointer (can be `NULL`), see `MatShift()`, or `MAT_SHELL_NOT_ALLOWED` if the internal shift should be 0
1992: . vscale - `PetscScalar` pointer (can be `NULL`), see `MatScale()`, or `MAT_SHELL_NOT_ALLOWED` if the internal scaling should be 1
1993: . dshift - `Vec` pointer (can be `NULL`), see `MatDiagonalSet()`, or `MAT_SHELL_NOT_ALLOWED` if the internal shift should be `NULL`
1994: . left - `Vec` pointer (can be `NULL`), see `MatDiagonalScale()`, or `MAT_SHELL_NOT_ALLOWED` if the internal scaling should be `NULL`
1995: . right - `Vec` pointer (can be `NULL`), see `MatDiagonalScale()`, or `MAT_SHELL_NOT_ALLOWED` if the internal scaling should be `NULL`
1996: . axpy - `Mat` pointer (can be `NULL`), or `MAT_SHELL_NOT_ALLOWED` if `MatAXPY()` should have not been called on `A`
1997: . zrows - `Vec` pointer (can be `NULL`), or `MAT_SHELL_NOT_ALLOWED` if `MatZeroRows()` should have not been called on `A`
1998: - zcols - `Vec` pointer (can be `NULL`), or `MAT_SHELL_NOT_ALLOWED` if `MatZeroRowsColumns()` should have not been called on `A`
2000: Level: advanced
2002: Developer Notes:
2003: This is mostly useful to check for corner-cases in `MatType` deriving from
2004: `MATSHELL`, e.g, `MATCOMPOSITE` or `MATTRANSPOSEVIRTUAL`, since scaling and
2005: shifts often require extra work which is not always implemented.
2007: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShift()`, `MatScale()`, `MatDiagonalSet()`, `MatDiagonalScale()`, `MatAXPY()`, `MatZeroRows()`, `MatZeroRowsColumns()`, `MatShellSetManageScalingShifts()`
2008: @*/
2009: PetscErrorCode MatShellGetScalingShifts(Mat A, PetscScalar *vshift, PetscScalar *vscale, Vec *dshift, Vec *left, Vec *right, Mat *axpy, IS *zrows, IS *zcols)
2010: {
2011: PetscFunctionBegin;
2013: PetscTryMethod(A, "MatShellGetScalingShifts_C", (Mat, PetscScalar *, PetscScalar *, Vec *, Vec *, Vec *, Mat *, IS *, IS *), (A, vshift, vscale, dshift, left, right, axpy, zrows, zcols));
2014: PetscFunctionReturn(PETSC_SUCCESS);
2015: }
2017: /*@C
2018: MatShellTestMult - Compares the multiply routine provided to the `MATSHELL` with differencing on a given function.
2020: Logically Collective; No Fortran Support
2022: Input Parameters:
2023: + mat - the `MATSHELL` shell matrix
2024: . f - the function
2025: . base - differences are computed around this vector, see `MatMFFDSetBase()`, for Jacobians this is the point at which the Jacobian is being evaluated
2026: - ctx - an optional context for the function
2028: Output Parameter:
2029: . flg - `PETSC_TRUE` if the multiply is likely correct
2031: Options Database Key:
2032: . -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference
2034: Level: advanced
2036: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMultTranspose()`
2037: @*/
2038: PetscErrorCode MatShellTestMult(Mat mat, PetscErrorCode (*f)(void *, Vec, Vec), Vec base, PetscCtx ctx, PetscBool *flg)
2039: {
2040: PetscInt m, n;
2041: Mat mf, Dmf, Dmat, Ddiff;
2042: PetscReal Diffnorm, Dmfnorm;
2043: PetscBool v = PETSC_FALSE, flag = PETSC_TRUE;
2045: PetscFunctionBegin;
2047: PetscCall(PetscOptionsHasName(NULL, ((PetscObject)mat)->prefix, "-mat_shell_test_mult_view", &v));
2048: PetscCall(MatGetLocalSize(mat, &m, &n));
2049: PetscCall(MatCreateMFFD(PetscObjectComm((PetscObject)mat), m, n, PETSC_DECIDE, PETSC_DECIDE, &mf));
2050: PetscCall(MatMFFDSetFunction(mf, f, ctx));
2051: PetscCall(MatMFFDSetBase(mf, base, NULL));
2053: PetscCall(MatComputeOperator(mf, MATAIJ, &Dmf));
2054: PetscCall(MatComputeOperator(mat, MATAIJ, &Dmat));
2056: PetscCall(MatDuplicate(Dmat, MAT_COPY_VALUES, &Ddiff));
2057: PetscCall(MatAXPY(Ddiff, -1.0, Dmf, DIFFERENT_NONZERO_PATTERN));
2058: PetscCall(MatNorm(Ddiff, NORM_FROBENIUS, &Diffnorm));
2059: PetscCall(MatNorm(Dmf, NORM_FROBENIUS, &Dmfnorm));
2060: if (Diffnorm / Dmfnorm > 10 * PETSC_SQRT_MACHINE_EPSILON) {
2061: flag = PETSC_FALSE;
2062: if (v) {
2063: 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)));
2064: PetscCall(MatViewFromOptions(Ddiff, (PetscObject)mat, "-mat_shell_test_mult_view"));
2065: PetscCall(MatViewFromOptions(Dmf, (PetscObject)mat, "-mat_shell_test_mult_view"));
2066: PetscCall(MatViewFromOptions(Dmat, (PetscObject)mat, "-mat_shell_test_mult_view"));
2067: }
2068: } else if (v) {
2069: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)mat), "MATSHELL and matrix-free multiple appear to produce the same results\n"));
2070: }
2071: if (flg) *flg = flag;
2072: PetscCall(MatDestroy(&Ddiff));
2073: PetscCall(MatDestroy(&mf));
2074: PetscCall(MatDestroy(&Dmf));
2075: PetscCall(MatDestroy(&Dmat));
2076: PetscFunctionReturn(PETSC_SUCCESS);
2077: }
2079: /*@C
2080: MatShellTestMultTranspose - Compares the multiply transpose routine provided to the `MATSHELL` with differencing on a given function.
2082: Logically Collective; No Fortran Support
2084: Input Parameters:
2085: + mat - the `MATSHELL` shell matrix
2086: . f - the function
2087: . base - differences are computed around this vector, see `MatMFFDSetBase()`, for Jacobians this is the point at which the Jacobian is being evaluated
2088: - ctx - an optional context for the function
2090: Output Parameter:
2091: . flg - `PETSC_TRUE` if the multiply is likely correct
2093: Options Database Key:
2094: . -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference
2096: Level: advanced
2098: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMult()`
2099: @*/
2100: PetscErrorCode MatShellTestMultTranspose(Mat mat, PetscErrorCode (*f)(void *, Vec, Vec), Vec base, PetscCtx ctx, PetscBool *flg)
2101: {
2102: Vec x, y, z;
2103: PetscInt m, n, M, N;
2104: Mat mf, Dmf, Dmat, Ddiff;
2105: PetscReal Diffnorm, Dmfnorm;
2106: PetscBool v = PETSC_FALSE, flag = PETSC_TRUE;
2108: PetscFunctionBegin;
2110: PetscCall(PetscOptionsHasName(NULL, ((PetscObject)mat)->prefix, "-mat_shell_test_mult_transpose_view", &v));
2111: PetscCall(MatCreateVecs(mat, &x, &y));
2112: PetscCall(VecDuplicate(y, &z));
2113: PetscCall(MatGetLocalSize(mat, &m, &n));
2114: PetscCall(MatGetSize(mat, &M, &N));
2115: PetscCall(MatCreateMFFD(PetscObjectComm((PetscObject)mat), m, n, M, N, &mf));
2116: PetscCall(MatMFFDSetFunction(mf, f, ctx));
2117: PetscCall(MatMFFDSetBase(mf, base, NULL));
2118: PetscCall(MatComputeOperator(mf, MATAIJ, &Dmf));
2119: PetscCall(MatTranspose(Dmf, MAT_INPLACE_MATRIX, &Dmf));
2120: PetscCall(MatComputeOperatorTranspose(mat, MATAIJ, &Dmat));
2122: PetscCall(MatDuplicate(Dmat, MAT_COPY_VALUES, &Ddiff));
2123: PetscCall(MatAXPY(Ddiff, -1.0, Dmf, DIFFERENT_NONZERO_PATTERN));
2124: PetscCall(MatNorm(Ddiff, NORM_FROBENIUS, &Diffnorm));
2125: PetscCall(MatNorm(Dmf, NORM_FROBENIUS, &Dmfnorm));
2126: if (Diffnorm / Dmfnorm > 10 * PETSC_SQRT_MACHINE_EPSILON) {
2127: flag = PETSC_FALSE;
2128: if (v) {
2129: 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)));
2130: PetscCall(MatViewFromOptions(Ddiff, (PetscObject)mat, "-mat_shell_test_mult_transpose_view"));
2131: PetscCall(MatViewFromOptions(Dmf, (PetscObject)mat, "-mat_shell_test_mult_transpose_view"));
2132: PetscCall(MatViewFromOptions(Dmat, (PetscObject)mat, "-mat_shell_test_mult_transpose_view"));
2133: }
2134: } else if (v) {
2135: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)mat), "MATSHELL transpose and matrix-free multiple appear to produce the same results\n"));
2136: }
2137: if (flg) *flg = flag;
2138: PetscCall(MatDestroy(&mf));
2139: PetscCall(MatDestroy(&Dmat));
2140: PetscCall(MatDestroy(&Ddiff));
2141: PetscCall(MatDestroy(&Dmf));
2142: PetscCall(VecDestroy(&x));
2143: PetscCall(VecDestroy(&y));
2144: PetscCall(VecDestroy(&z));
2145: PetscFunctionReturn(PETSC_SUCCESS);
2146: }
2148: /*@C
2149: MatShellSetOperation - Allows user to set a matrix operation for a `MATSHELL` shell matrix.
2151: Logically Collective
2153: Input Parameters:
2154: + mat - the `MATSHELL` shell matrix
2155: . op - the name of the operation
2156: - g - the function that provides the operation.
2158: Level: advanced
2160: Example Usage:
2161: .vb
2162: extern PetscErrorCode usermult(Mat, Vec, Vec);
2164: MatCreateShell(comm, m, n, M, N, ctx, &A);
2165: MatShellSetOperation(A, MATOP_MULT, (void(*)(void))usermult);
2166: .ve
2168: Notes:
2169: See the file include/petscmat.h for a complete list of matrix
2170: operations, which all have the form MATOP_<OPERATION>, where
2171: <OPERATION> is the name (in all capital letters) of the
2172: user interface routine (e.g., `MatMult()` -> `MATOP_MULT`).
2174: All user-provided functions (except for `MATOP_DESTROY`) should have the same calling
2175: sequence as the usual matrix interface routines, since they
2176: are intended to be accessed via the usual matrix interface
2177: routines, e.g.,
2178: .vb
2179: MatMult(Mat, Vec, Vec) -> usermult(Mat, Vec, Vec)
2180: .ve
2182: In particular each function MUST return an error code of `PETSC_SUCCESS` on success and
2183: nonzero on failure.
2185: Within each user-defined routine, the user should call
2186: `MatShellGetContext()` to obtain the user-defined context that was
2187: set by `MatCreateShell()`.
2189: Use `MatSetOperation()` to set an operation for any matrix type. For matrix product operations (i.e. `MatMatXXX()`, `MatTransposeMatXXX()` etc)
2190: use `MatShellSetMatProductOperation()`
2192: Fortran Note:
2193: For `MatCreateVecs()` the user code should check if the input left or right vector is -1 and in that case not
2194: generate a vector. See `src/mat/tests/ex120f.F`
2196: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellSetContext()`, `MatSetOperation()`, `MatShellSetManageScalingShifts()`, `MatShellSetMatProductOperation()`
2197: @*/
2198: PetscErrorCode MatShellSetOperation(Mat mat, MatOperation op, PetscErrorCodeFn *g)
2199: {
2200: PetscFunctionBegin;
2202: PetscTryMethod(mat, "MatShellSetOperation_C", (Mat, MatOperation, PetscErrorCodeFn *), (mat, op, g));
2203: PetscFunctionReturn(PETSC_SUCCESS);
2204: }
2206: /*@C
2207: MatShellGetOperation - Gets a matrix function for a `MATSHELL` shell matrix.
2209: Not Collective
2211: Input Parameters:
2212: + mat - the `MATSHELL` shell matrix
2213: - op - the name of the operation
2215: Output Parameter:
2216: . g - the function that provides the operation.
2218: Level: advanced
2220: Notes:
2221: See the file include/petscmat.h for a complete list of matrix
2222: operations, which all have the form MATOP_<OPERATION>, where
2223: <OPERATION> is the name (in all capital letters) of the
2224: user interface routine (e.g., `MatMult()` -> `MATOP_MULT`).
2226: All user-provided functions have the same calling
2227: sequence as the usual matrix interface routines, since they
2228: are intended to be accessed via the usual matrix interface
2229: routines, e.g.,
2230: .vb
2231: MatMult(Mat, Vec, Vec) -> usermult(Mat, Vec, Vec)
2232: .ve
2234: Within each user-defined routine, the user should call
2235: `MatShellGetContext()` to obtain the user-defined context that was
2236: set by `MatCreateShell()`.
2238: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellSetOperation()`, `MatShellSetContext()`
2239: @*/
2240: PetscErrorCode MatShellGetOperation(Mat mat, MatOperation op, PetscErrorCodeFn **g)
2241: {
2242: PetscFunctionBegin;
2244: PetscUseMethod(mat, "MatShellGetOperation_C", (Mat, MatOperation, PetscErrorCodeFn **), (mat, op, g));
2245: PetscFunctionReturn(PETSC_SUCCESS);
2246: }
2248: /*@
2249: MatIsShell - Inquires if a matrix is derived from `MATSHELL`
2251: Input Parameter:
2252: . mat - the matrix
2254: Output Parameter:
2255: . flg - the Boolean value
2257: Level: developer
2259: .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MATMFFD`, `MatCreateShell()`, `MATTRANSPOSEVIRTUAL`, `MATSCHURCOMPLEMENT`
2260: @*/
2261: PetscErrorCode MatIsShell(Mat mat, PetscBool *flg)
2262: {
2263: PetscFunctionBegin;
2265: PetscAssertPointer(flg, 2);
2266: *flg = (PetscBool)(mat->ops->destroy == MatDestroy_Shell);
2267: PetscFunctionReturn(PETSC_SUCCESS);
2268: }