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