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