Actual source code: shell.c
petsc-3.12.5 2020-03-29
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: typedef struct {
20: struct _MatShellOps ops[1];
22: PetscScalar vscale,vshift;
23: Vec dshift;
24: Vec left,right;
25: Vec left_work,right_work;
26: Vec left_add_work,right_add_work;
27: Mat axpy;
28: PetscScalar axpy_vscale;
29: PetscBool managescalingshifts; /* The user will manage the scaling and shifts for the MATSHELL, not the default */
30: /* support for ZeroRows/Columns operations */
31: IS zrows;
32: IS zcols;
33: Vec zvals;
34: Vec zvals_w;
35: VecScatter zvals_sct_r;
36: VecScatter zvals_sct_c;
37: void *ctx;
38: } Mat_Shell;
41: /*
42: Store and scale values on zeroed rows
43: xx = [x_1, 0], 0 on zeroed columns
44: */
45: static PetscErrorCode MatShellPreZeroRight(Mat A,Vec x,Vec *xx)
46: {
47: Mat_Shell *shell = (Mat_Shell*)A->data;
51: *xx = x;
52: if (shell->zrows) {
53: VecSet(shell->zvals_w,0.0);
54: VecScatterBegin(shell->zvals_sct_c,x,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);
55: VecScatterEnd(shell->zvals_sct_c,x,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);
56: VecPointwiseMult(shell->zvals_w,shell->zvals_w,shell->zvals);
57: }
58: if (shell->zcols) {
59: if (!shell->right_work) {
60: MatCreateVecs(A,&shell->right_work,NULL);
61: }
62: VecCopy(x,shell->right_work);
63: VecISSet(shell->right_work,shell->zcols,0.0);
64: *xx = shell->right_work;
65: }
66: return(0);
67: }
69: /* Insert properly diagonally scaled values stored in MatShellPreZeroRight */
70: static PetscErrorCode MatShellPostZeroLeft(Mat A,Vec x)
71: {
72: Mat_Shell *shell = (Mat_Shell*)A->data;
76: if (shell->zrows) {
77: VecScatterBegin(shell->zvals_sct_r,shell->zvals_w,x,INSERT_VALUES,SCATTER_REVERSE);
78: VecScatterEnd(shell->zvals_sct_r,shell->zvals_w,x,INSERT_VALUES,SCATTER_REVERSE);
79: }
80: return(0);
81: }
83: /*
84: Store and scale values on zeroed rows
85: xx = [x_1, 0], 0 on zeroed rows
86: */
87: static PetscErrorCode MatShellPreZeroLeft(Mat A,Vec x,Vec *xx)
88: {
89: Mat_Shell *shell = (Mat_Shell*)A->data;
93: *xx = NULL;
94: if (!shell->zrows) {
95: *xx = x;
96: } else {
97: if (!shell->left_work) {
98: MatCreateVecs(A,NULL,&shell->left_work);
99: }
100: VecCopy(x,shell->left_work);
101: VecSet(shell->zvals_w,0.0);
102: VecScatterBegin(shell->zvals_sct_r,shell->zvals_w,shell->left_work,INSERT_VALUES,SCATTER_REVERSE);
103: VecScatterEnd(shell->zvals_sct_r,shell->zvals_w,shell->left_work,INSERT_VALUES,SCATTER_REVERSE);
104: VecScatterBegin(shell->zvals_sct_r,x,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);
105: VecScatterEnd(shell->zvals_sct_r,x,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);
106: VecPointwiseMult(shell->zvals_w,shell->zvals_w,shell->zvals);
107: *xx = shell->left_work;
108: }
109: return(0);
110: }
112: /* Zero zero-columns contributions, sum contributions from properly scaled values stored in MatShellPreZeroLeft */
113: static PetscErrorCode MatShellPostZeroRight(Mat A,Vec x)
114: {
115: Mat_Shell *shell = (Mat_Shell*)A->data;
119: if (shell->zcols) {
120: VecISSet(x,shell->zcols,0.0);
121: }
122: if (shell->zrows) {
123: VecScatterBegin(shell->zvals_sct_c,shell->zvals_w,x,ADD_VALUES,SCATTER_REVERSE);
124: VecScatterEnd(shell->zvals_sct_c,shell->zvals_w,x,ADD_VALUES,SCATTER_REVERSE);
125: }
126: return(0);
127: }
129: /*
130: xx = diag(left)*x
131: */
132: static PetscErrorCode MatShellPreScaleLeft(Mat A,Vec x,Vec *xx)
133: {
134: Mat_Shell *shell = (Mat_Shell*)A->data;
138: *xx = NULL;
139: if (!shell->left) {
140: *xx = x;
141: } else {
142: if (!shell->left_work) {VecDuplicate(shell->left,&shell->left_work);}
143: VecPointwiseMult(shell->left_work,x,shell->left);
144: *xx = shell->left_work;
145: }
146: return(0);
147: }
149: /*
150: xx = diag(right)*x
151: */
152: static PetscErrorCode MatShellPreScaleRight(Mat A,Vec x,Vec *xx)
153: {
154: Mat_Shell *shell = (Mat_Shell*)A->data;
158: *xx = NULL;
159: if (!shell->right) {
160: *xx = x;
161: } else {
162: if (!shell->right_work) {VecDuplicate(shell->right,&shell->right_work);}
163: VecPointwiseMult(shell->right_work,x,shell->right);
164: *xx = shell->right_work;
165: }
166: return(0);
167: }
169: /*
170: x = diag(left)*x
171: */
172: static PetscErrorCode MatShellPostScaleLeft(Mat A,Vec x)
173: {
174: Mat_Shell *shell = (Mat_Shell*)A->data;
178: if (shell->left) {VecPointwiseMult(x,x,shell->left);}
179: return(0);
180: }
182: /*
183: x = diag(right)*x
184: */
185: static PetscErrorCode MatShellPostScaleRight(Mat A,Vec x)
186: {
187: Mat_Shell *shell = (Mat_Shell*)A->data;
191: if (shell->right) {VecPointwiseMult(x,x,shell->right);}
192: return(0);
193: }
195: /*
196: Y = vscale*Y + diag(dshift)*X + vshift*X
198: On input Y already contains A*x
199: */
200: static PetscErrorCode MatShellShiftAndScale(Mat A,Vec X,Vec Y)
201: {
202: 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: PetscErrorCode MatShellGetContext_Shell(Mat mat,void *ctx)
226: {
228: *(void**)ctx = ((Mat_Shell*)(mat->data))->ctx;
229: return(0);
230: }
232: /*@
233: MatShellGetContext - Returns the user-provided context associated with a shell matrix.
235: Not Collective
237: Input Parameter:
238: . mat - the matrix, should have been created with MatCreateShell()
240: Output Parameter:
241: . ctx - the user provided context
243: Level: advanced
245: Fortran Notes:
246: To use this from Fortran you must write a Fortran interface definition for this
247: function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.
249: .seealso: MatCreateShell(), MatShellSetOperation(), MatShellSetContext()
250: @*/
251: PetscErrorCode MatShellGetContext(Mat mat,void *ctx)
252: {
258: PetscUseMethod(mat,"MatShellGetContext_C",(Mat,void*),(mat,ctx));
259: return(0);
260: }
262: static PetscErrorCode MatZeroRowsColumns_Local_Shell(Mat mat,PetscInt nr,PetscInt rows[],PetscInt nc,PetscInt cols[],PetscScalar diag,PetscBool rc)
263: {
265: Mat_Shell *shell = (Mat_Shell*)mat->data;
266: Vec x = NULL,b = NULL;
267: IS is1, is2;
268: const PetscInt *ridxs;
269: PetscInt *idxs,*gidxs;
270: PetscInt cum,rst,cst,i;
273: if (shell->axpy) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot diagonal scale MATSHELL after MatAXPY operation");
274: if (!shell->zvals) {
275: MatCreateVecs(mat,NULL,&shell->zvals);
276: }
277: if (!shell->zvals_w) {
278: VecDuplicate(shell->zvals,&shell->zvals_w);
279: }
280: MatGetOwnershipRange(mat,&rst,NULL);
281: MatGetOwnershipRangeColumn(mat,&cst,NULL);
283: /* Expand/create index set of zeroed rows */
284: PetscMalloc1(nr,&idxs);
285: for (i = 0; i < nr; i++) idxs[i] = rows[i] + rst;
286: ISCreateGeneral(PETSC_COMM_SELF,nr,idxs,PETSC_OWN_POINTER,&is1);
287: ISSort(is1);
288: VecISSet(shell->zvals,is1,diag);
289: if (shell->zrows) {
290: ISSum(shell->zrows,is1,&is2);
291: ISDestroy(&shell->zrows);
292: ISDestroy(&is1);
293: shell->zrows = is2;
294: } else shell->zrows = is1;
296: /* Create scatters for diagonal values communications */
297: VecScatterDestroy(&shell->zvals_sct_c);
298: VecScatterDestroy(&shell->zvals_sct_r);
300: /* row scatter: from/to left vector */
301: MatCreateVecs(mat,&x,&b);
302: VecScatterCreate(b,shell->zrows,shell->zvals_w,shell->zrows,&shell->zvals_sct_r);
304: /* col scatter: from right vector to left vector */
305: ISGetIndices(shell->zrows,&ridxs);
306: ISGetLocalSize(shell->zrows,&nr);
307: PetscMalloc1(nr,&gidxs);
308: for (i = 0, cum = 0; i < nr; i++) {
309: if (ridxs[i] >= mat->cmap->N) continue;
310: gidxs[cum] = ridxs[i];
311: cum++;
312: }
313: ISCreateGeneral(PetscObjectComm((PetscObject)mat),cum,gidxs,PETSC_OWN_POINTER,&is1);
314: VecScatterCreate(x,is1,shell->zvals_w,is1,&shell->zvals_sct_c);
315: ISDestroy(&is1);
316: VecDestroy(&x);
317: VecDestroy(&b);
319: /* Expand/create index set of zeroed columns */
320: if (rc) {
321: PetscMalloc1(nc,&idxs);
322: for (i = 0; i < nc; i++) idxs[i] = cols[i] + cst;
323: ISCreateGeneral(PETSC_COMM_SELF,nc,idxs,PETSC_OWN_POINTER,&is1);
324: ISSort(is1);
325: if (shell->zcols) {
326: ISSum(shell->zcols,is1,&is2);
327: ISDestroy(&shell->zcols);
328: ISDestroy(&is1);
329: shell->zcols = is2;
330: } else shell->zcols = is1;
331: }
332: return(0);
333: }
335: static PetscErrorCode MatZeroRows_Shell(Mat mat,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
336: {
337: PetscInt nr, *lrows;
341: if (x && b) {
342: Vec xt;
343: PetscScalar *vals;
344: PetscInt *gcols,i,st,nl,nc;
346: PetscMalloc1(n,&gcols);
347: for (i = 0, nc = 0; i < n; i++) if (rows[i] < mat->cmap->N) gcols[nc++] = rows[i];
349: MatCreateVecs(mat,&xt,NULL);
350: VecCopy(x,xt);
351: PetscCalloc1(nc,&vals);
352: VecSetValues(xt,nc,gcols,vals,INSERT_VALUES); /* xt = [x1, 0] */
353: PetscFree(vals);
354: VecAssemblyBegin(xt);
355: VecAssemblyEnd(xt);
356: VecAYPX(xt,-1.0,x); /* xt = [0, x2] */
358: VecGetOwnershipRange(xt,&st,NULL);
359: VecGetLocalSize(xt,&nl);
360: VecGetArray(xt,&vals);
361: for (i = 0; i < nl; i++) {
362: PetscInt g = i + st;
363: if (g > mat->rmap->N) continue;
364: if (PetscAbsScalar(vals[i]) == 0.0) continue;
365: VecSetValue(b,g,diag*vals[i],INSERT_VALUES);
366: }
367: VecRestoreArray(xt,&vals);
368: VecAssemblyBegin(b);
369: VecAssemblyEnd(b); /* b = [b1, x2 * diag] */
370: VecDestroy(&xt);
371: PetscFree(gcols);
372: }
373: PetscLayoutMapLocal(mat->rmap,n,rows,&nr,&lrows,NULL);
374: MatZeroRowsColumns_Local_Shell(mat,nr,lrows,0,NULL,diag,PETSC_FALSE);
375: PetscFree(lrows);
376: return(0);
377: }
379: static PetscErrorCode MatZeroRowsColumns_Shell(Mat mat,PetscInt n,const PetscInt rowscols[],PetscScalar diag,Vec x,Vec b)
380: {
381: PetscInt *lrows, *lcols;
382: PetscInt nr, nc;
383: 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: return(0);
448: }
450: PetscErrorCode MatDestroy_Shell(Mat mat)
451: {
453: Mat_Shell *shell = (Mat_Shell*)mat->data;
456: if (shell->ops->destroy) {
457: (*shell->ops->destroy)(mat);
458: }
459: VecDestroy(&shell->left);
460: VecDestroy(&shell->right);
461: VecDestroy(&shell->dshift);
462: VecDestroy(&shell->left_work);
463: VecDestroy(&shell->right_work);
464: VecDestroy(&shell->left_add_work);
465: VecDestroy(&shell->right_add_work);
466: MatDestroy(&shell->axpy);
467: VecDestroy(&shell->zvals_w);
468: VecDestroy(&shell->zvals);
469: VecScatterDestroy(&shell->zvals_sct_c);
470: VecScatterDestroy(&shell->zvals_sct_r);
471: ISDestroy(&shell->zrows);
472: ISDestroy(&shell->zcols);
473: PetscFree(mat->data);
474: PetscObjectComposeFunction((PetscObject)mat,"MatShellGetContext_C",NULL);
475: PetscObjectComposeFunction((PetscObject)mat,"MatShellSetContext_C",NULL);
476: PetscObjectComposeFunction((PetscObject)mat,"MatShellSetManageScalingShifts_C",NULL);
477: PetscObjectComposeFunction((PetscObject)mat,"MatShellSetOperation_C",NULL);
478: PetscObjectComposeFunction((PetscObject)mat,"MatShellGetOperation_C",NULL);
479: return(0);
480: }
482: PetscErrorCode MatCopy_Shell(Mat A,Mat B,MatStructure str)
483: {
484: Mat_Shell *shellA = (Mat_Shell*)A->data,*shellB = (Mat_Shell*)B->data;
485: PetscErrorCode ierr;
486: PetscBool matflg;
489: PetscObjectTypeCompare((PetscObject)B,MATSHELL,&matflg);
490: if (!matflg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_NOTSAMETYPE,"Output matrix must be a MATSHELL");
492: PetscMemcpy(B->ops,A->ops,sizeof(struct _MatOps));
493: PetscMemcpy(shellB->ops,shellA->ops,sizeof(struct _MatShellOps));
495: if (shellA->ops->copy) {
496: (*shellA->ops->copy)(A,B,str);
497: }
498: shellB->vscale = shellA->vscale;
499: shellB->vshift = shellA->vshift;
500: if (shellA->dshift) {
501: if (!shellB->dshift) {
502: VecDuplicate(shellA->dshift,&shellB->dshift);
503: }
504: VecCopy(shellA->dshift,shellB->dshift);
505: } else {
506: VecDestroy(&shellB->dshift);
507: }
508: if (shellA->left) {
509: if (!shellB->left) {
510: VecDuplicate(shellA->left,&shellB->left);
511: }
512: VecCopy(shellA->left,shellB->left);
513: } else {
514: VecDestroy(&shellB->left);
515: }
516: if (shellA->right) {
517: if (!shellB->right) {
518: VecDuplicate(shellA->right,&shellB->right);
519: }
520: VecCopy(shellA->right,shellB->right);
521: } else {
522: VecDestroy(&shellB->right);
523: }
524: MatDestroy(&shellB->axpy);
525: if (shellA->axpy) {
526: PetscObjectReference((PetscObject)shellA->axpy);
527: shellB->axpy = shellA->axpy;
528: shellB->axpy_vscale = shellA->axpy_vscale;
529: }
530: if (shellA->zrows) {
531: ISDuplicate(shellA->zrows,&shellB->zrows);
532: if (shellA->zcols) {
533: ISDuplicate(shellA->zcols,&shellB->zcols);
534: }
535: VecDuplicate(shellA->zvals,&shellB->zvals);
536: VecCopy(shellA->zvals,shellB->zvals);
537: VecDuplicate(shellA->zvals_w,&shellB->zvals_w);
538: PetscObjectReference((PetscObject)shellA->zvals_sct_r);
539: PetscObjectReference((PetscObject)shellA->zvals_sct_c);
540: shellB->zvals_sct_r = shellA->zvals_sct_r;
541: shellB->zvals_sct_c = shellA->zvals_sct_c;
542: }
543: return(0);
544: }
546: PetscErrorCode MatDuplicate_Shell(Mat mat,MatDuplicateOption op,Mat *M)
547: {
549: void *ctx;
552: MatShellGetContext(mat,&ctx);
553: MatCreateShell(PetscObjectComm((PetscObject)mat),mat->rmap->n,mat->cmap->n,mat->rmap->N,mat->cmap->N,ctx,M);
554: if (op != MAT_DO_NOT_COPY_VALUES) {
555: MatCopy(mat,*M,SAME_NONZERO_PATTERN);
556: }
557: return(0);
558: }
560: PetscErrorCode MatMult_Shell(Mat A,Vec x,Vec y)
561: {
562: Mat_Shell *shell = (Mat_Shell*)A->data;
563: PetscErrorCode ierr;
564: Vec xx;
565: PetscObjectState instate,outstate;
568: if (!shell->ops->mult) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Have not provided a MatMult() for this MATSHELL");
569: MatShellPreZeroRight(A,x,&xx);
570: MatShellPreScaleRight(A,xx,&xx);
571: PetscObjectStateGet((PetscObject)y, &instate);
572: (*shell->ops->mult)(A,xx,y);
573: PetscObjectStateGet((PetscObject)y, &outstate);
574: if (instate == outstate) {
575: /* increase the state of the output vector since the user did not update its state themself as should have been done */
576: PetscObjectStateIncrease((PetscObject)y);
577: }
578: MatShellShiftAndScale(A,xx,y);
579: MatShellPostScaleLeft(A,y);
580: MatShellPostZeroLeft(A,y);
582: if (shell->axpy) {
583: if (!shell->left_work) {MatCreateVecs(A,&shell->left_work,NULL);}
584: MatMult(shell->axpy,x,shell->left_work);
585: VecAXPY(y,shell->axpy_vscale,shell->left_work);
586: }
587: return(0);
588: }
590: PetscErrorCode MatMultAdd_Shell(Mat A,Vec x,Vec y,Vec z)
591: {
592: Mat_Shell *shell = (Mat_Shell*)A->data;
596: if (y == z) {
597: if (!shell->right_add_work) {VecDuplicate(z,&shell->right_add_work);}
598: MatMult(A,x,shell->right_add_work);
599: VecAXPY(z,1.0,shell->right_add_work);
600: } else {
601: MatMult(A,x,z);
602: VecAXPY(z,1.0,y);
603: }
604: return(0);
605: }
607: PetscErrorCode MatMultTranspose_Shell(Mat A,Vec x,Vec y)
608: {
609: Mat_Shell *shell = (Mat_Shell*)A->data;
610: PetscErrorCode ierr;
611: Vec xx;
612: PetscObjectState instate,outstate;
615: MatShellPreZeroLeft(A,x,&xx);
616: MatShellPreScaleLeft(A,xx,&xx);
617: PetscObjectStateGet((PetscObject)y, &instate);
618: if (!shell->ops->multtranspose) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Have not provided a MatMultTranspose() for this MATSHELL");
619: (*shell->ops->multtranspose)(A,xx,y);
620: PetscObjectStateGet((PetscObject)y, &outstate);
621: if (instate == outstate) {
622: /* increase the state of the output vector since the user did not update its state themself as should have been done */
623: PetscObjectStateIncrease((PetscObject)y);
624: }
625: MatShellShiftAndScale(A,xx,y);
626: MatShellPostScaleRight(A,y);
627: MatShellPostZeroRight(A,y);
629: if (shell->axpy) {
630: if (!shell->right_work) {MatCreateVecs(A,NULL,&shell->right_work);}
631: MatMultTranspose(shell->axpy,x,shell->right_work);
632: VecAXPY(y,shell->axpy_vscale,shell->right_work);
633: }
634: return(0);
635: }
637: PetscErrorCode MatMultTransposeAdd_Shell(Mat A,Vec x,Vec y,Vec z)
638: {
639: Mat_Shell *shell = (Mat_Shell*)A->data;
643: if (y == z) {
644: if (!shell->left_add_work) {VecDuplicate(z,&shell->left_add_work);}
645: MatMultTranspose(A,x,shell->left_add_work);
646: VecAXPY(z,1.0,shell->left_add_work);
647: } else {
648: MatMultTranspose(A,x,z);
649: VecAXPY(z,1.0,y);
650: }
651: return(0);
652: }
654: /*
655: diag(left)(vscale*A + diag(dshift) + vshift I)diag(right)
656: */
657: PetscErrorCode MatGetDiagonal_Shell(Mat A,Vec v)
658: {
659: Mat_Shell *shell = (Mat_Shell*)A->data;
663: if (shell->ops->getdiagonal) {
664: (*shell->ops->getdiagonal)(A,v);
665: } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Must provide shell matrix with routine to return diagonal using\nMatShellSetOperation(S,MATOP_GET_DIAGONAL,...)");
666: VecScale(v,shell->vscale);
667: if (shell->dshift) {
668: VecAXPY(v,1.0,shell->dshift);
669: }
670: VecShift(v,shell->vshift);
671: if (shell->left) {VecPointwiseMult(v,v,shell->left);}
672: if (shell->right) {VecPointwiseMult(v,v,shell->right);}
673: if (shell->zrows) {
674: VecScatterBegin(shell->zvals_sct_r,shell->zvals,v,INSERT_VALUES,SCATTER_REVERSE);
675: VecScatterEnd(shell->zvals_sct_r,shell->zvals,v,INSERT_VALUES,SCATTER_REVERSE);
676: }
677: if (shell->axpy) {
678: if (!shell->left_work) {VecDuplicate(v,&shell->left_work);}
679: MatGetDiagonal(shell->axpy,shell->left_work);
680: VecAXPY(v,shell->axpy_vscale,shell->left_work);
681: }
682: return(0);
683: }
685: PetscErrorCode MatShift_Shell(Mat Y,PetscScalar a)
686: {
687: Mat_Shell *shell = (Mat_Shell*)Y->data;
689: PetscBool flg;
692: MatHasCongruentLayouts(Y,&flg);
693: if (!flg) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_SUP,"Cannot shift shell matrix if it is not congruent");
694: if (shell->left || shell->right) {
695: if (!shell->dshift) {
696: VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift);
697: VecSet(shell->dshift,a);
698: } else {
699: if (shell->left) {VecPointwiseMult(shell->dshift,shell->dshift,shell->left);}
700: if (shell->right) {VecPointwiseMult(shell->dshift,shell->dshift,shell->right);}
701: VecShift(shell->dshift,a);
702: }
703: if (shell->left) {VecPointwiseDivide(shell->dshift,shell->dshift,shell->left);}
704: if (shell->right) {VecPointwiseDivide(shell->dshift,shell->dshift,shell->right);}
705: } else shell->vshift += a;
706: if (shell->zrows) {
707: VecShift(shell->zvals,a);
708: }
709: return(0);
710: }
712: PetscErrorCode MatDiagonalSet_Shell_Private(Mat A,Vec D,PetscScalar s)
713: {
714: Mat_Shell *shell = (Mat_Shell*)A->data;
718: if (!shell->dshift) {VecDuplicate(D,&shell->dshift);}
719: if (shell->left || shell->right) {
720: if (!shell->right_work) {VecDuplicate(shell->left ? shell->left : shell->right, &shell->right_work);}
721: if (shell->left && shell->right) {
722: VecPointwiseDivide(shell->right_work,D,shell->left);
723: VecPointwiseDivide(shell->right_work,shell->right_work,shell->right);
724: } else if (shell->left) {
725: VecPointwiseDivide(shell->right_work,D,shell->left);
726: } else {
727: VecPointwiseDivide(shell->right_work,D,shell->right);
728: }
729: VecAXPY(shell->dshift,s,shell->right_work);
730: } else {
731: VecAXPY(shell->dshift,s,D);
732: }
733: return(0);
734: }
736: PetscErrorCode MatDiagonalSet_Shell(Mat A,Vec D,InsertMode ins)
737: {
738: Mat_Shell *shell = (Mat_Shell*)A->data;
739: Vec d;
741: PetscBool flg;
744: MatHasCongruentLayouts(A,&flg);
745: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Cannot diagonal set or shift shell matrix if it is not congruent");
746: if (ins == INSERT_VALUES) {
747: if (!A->ops->getdiagonal) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Operation MATOP_GETDIAGONAL must be set first");
748: VecDuplicate(D,&d);
749: MatGetDiagonal(A,d);
750: MatDiagonalSet_Shell_Private(A,d,-1.);
751: MatDiagonalSet_Shell_Private(A,D,1.);
752: VecDestroy(&d);
753: if (shell->zrows) {
754: VecCopy(D,shell->zvals);
755: }
756: } else {
757: MatDiagonalSet_Shell_Private(A,D,1.);
758: if (shell->zrows) {
759: VecAXPY(shell->zvals,1.0,D);
760: }
761: }
762: return(0);
763: }
765: PetscErrorCode MatScale_Shell(Mat Y,PetscScalar a)
766: {
767: Mat_Shell *shell = (Mat_Shell*)Y->data;
771: shell->vscale *= a;
772: shell->vshift *= a;
773: if (shell->dshift) {
774: VecScale(shell->dshift,a);
775: }
776: shell->axpy_vscale *= a;
777: if (shell->zrows) {
778: VecScale(shell->zvals,a);
779: }
780: return(0);
781: }
783: static PetscErrorCode MatDiagonalScale_Shell(Mat Y,Vec left,Vec right)
784: {
785: Mat_Shell *shell = (Mat_Shell*)Y->data;
789: if (shell->axpy) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_SUP,"Cannot diagonal scale MATSHELL after MatAXPY operation");
790: if (left) {
791: if (!shell->left) {
792: VecDuplicate(left,&shell->left);
793: VecCopy(left,shell->left);
794: } else {
795: VecPointwiseMult(shell->left,shell->left,left);
796: }
797: if (shell->zrows) {
798: VecPointwiseMult(shell->zvals,shell->zvals,left);
799: }
800: }
801: if (right) {
802: if (!shell->right) {
803: VecDuplicate(right,&shell->right);
804: VecCopy(right,shell->right);
805: } else {
806: VecPointwiseMult(shell->right,shell->right,right);
807: }
808: if (shell->zrows) {
809: if (!shell->left_work) {
810: MatCreateVecs(Y,NULL,&shell->left_work);
811: }
812: VecSet(shell->zvals_w,1.0);
813: VecScatterBegin(shell->zvals_sct_c,right,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);
814: VecScatterEnd(shell->zvals_sct_c,right,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);
815: VecPointwiseMult(shell->zvals,shell->zvals,shell->zvals_w);
816: }
817: }
818: return(0);
819: }
821: PetscErrorCode MatAssemblyEnd_Shell(Mat Y,MatAssemblyType t)
822: {
823: Mat_Shell *shell = (Mat_Shell*)Y->data;
827: if (t == MAT_FINAL_ASSEMBLY) {
828: shell->vshift = 0.0;
829: shell->vscale = 1.0;
830: VecDestroy(&shell->dshift);
831: VecDestroy(&shell->left);
832: VecDestroy(&shell->right);
833: MatDestroy(&shell->axpy);
834: VecScatterDestroy(&shell->zvals_sct_c);
835: VecScatterDestroy(&shell->zvals_sct_r);
836: ISDestroy(&shell->zrows);
837: ISDestroy(&shell->zcols);
838: }
839: return(0);
840: }
842: static PetscErrorCode MatMissingDiagonal_Shell(Mat A,PetscBool *missing,PetscInt *d)
843: {
845: *missing = PETSC_FALSE;
846: return(0);
847: }
849: PetscErrorCode MatAXPY_Shell(Mat Y,PetscScalar a,Mat X,MatStructure str)
850: {
851: Mat_Shell *shell = (Mat_Shell*)Y->data;
855: PetscObjectReference((PetscObject)X);
856: MatDestroy(&shell->axpy);
857: shell->axpy = X;
858: shell->axpy_vscale = a;
859: return(0);
860: }
862: static struct _MatOps MatOps_Values = {0,
863: 0,
864: 0,
865: 0,
866: /* 4*/ MatMultAdd_Shell,
867: 0,
868: MatMultTransposeAdd_Shell,
869: 0,
870: 0,
871: 0,
872: /*10*/ 0,
873: 0,
874: 0,
875: 0,
876: 0,
877: /*15*/ 0,
878: 0,
879: 0,
880: MatDiagonalScale_Shell,
881: 0,
882: /*20*/ 0,
883: MatAssemblyEnd_Shell,
884: 0,
885: 0,
886: /*24*/ MatZeroRows_Shell,
887: 0,
888: 0,
889: 0,
890: 0,
891: /*29*/ 0,
892: 0,
893: 0,
894: 0,
895: 0,
896: /*34*/ MatDuplicate_Shell,
897: 0,
898: 0,
899: 0,
900: 0,
901: /*39*/ MatAXPY_Shell,
902: 0,
903: 0,
904: 0,
905: MatCopy_Shell,
906: /*44*/ 0,
907: MatScale_Shell,
908: MatShift_Shell,
909: MatDiagonalSet_Shell,
910: MatZeroRowsColumns_Shell,
911: /*49*/ 0,
912: 0,
913: 0,
914: 0,
915: 0,
916: /*54*/ 0,
917: 0,
918: 0,
919: 0,
920: 0,
921: /*59*/ 0,
922: MatDestroy_Shell,
923: 0,
924: MatConvertFrom_Shell,
925: 0,
926: /*64*/ 0,
927: 0,
928: 0,
929: 0,
930: 0,
931: /*69*/ 0,
932: 0,
933: MatConvert_Shell,
934: 0,
935: 0,
936: /*74*/ 0,
937: 0,
938: 0,
939: 0,
940: 0,
941: /*79*/ 0,
942: 0,
943: 0,
944: 0,
945: 0,
946: /*84*/ 0,
947: 0,
948: 0,
949: 0,
950: 0,
951: /*89*/ 0,
952: 0,
953: 0,
954: 0,
955: 0,
956: /*94*/ 0,
957: 0,
958: 0,
959: 0,
960: 0,
961: /*99*/ 0,
962: 0,
963: 0,
964: 0,
965: 0,
966: /*104*/ 0,
967: 0,
968: 0,
969: 0,
970: 0,
971: /*109*/ 0,
972: 0,
973: 0,
974: 0,
975: MatMissingDiagonal_Shell,
976: /*114*/ 0,
977: 0,
978: 0,
979: 0,
980: 0,
981: /*119*/ 0,
982: 0,
983: 0,
984: 0,
985: 0,
986: /*124*/ 0,
987: 0,
988: 0,
989: 0,
990: 0,
991: /*129*/ 0,
992: 0,
993: 0,
994: 0,
995: 0,
996: /*134*/ 0,
997: 0,
998: 0,
999: 0,
1000: 0,
1001: /*139*/ 0,
1002: 0,
1003: 0
1004: };
1006: PetscErrorCode MatShellSetContext_Shell(Mat mat,void *ctx)
1007: {
1008: Mat_Shell *shell = (Mat_Shell*)mat->data;
1011: shell->ctx = ctx;
1012: return(0);
1013: }
1015: PetscErrorCode MatShellSetManageScalingShifts_Shell(Mat A)
1016: {
1017: Mat_Shell *shell = (Mat_Shell*)A->data;
1020: shell->managescalingshifts = PETSC_FALSE;
1021: A->ops->diagonalset = NULL;
1022: A->ops->diagonalscale = NULL;
1023: A->ops->scale = NULL;
1024: A->ops->shift = NULL;
1025: A->ops->axpy = NULL;
1026: return(0);
1027: }
1029: PetscErrorCode MatShellSetOperation_Shell(Mat mat,MatOperation op,void (*f)(void))
1030: {
1031: Mat_Shell *shell = (Mat_Shell*)mat->data;
1034: switch (op) {
1035: case MATOP_DESTROY:
1036: shell->ops->destroy = (PetscErrorCode (*)(Mat))f;
1037: break;
1038: case MATOP_VIEW:
1039: if (!mat->ops->viewnative) {
1040: mat->ops->viewnative = mat->ops->view;
1041: }
1042: mat->ops->view = (PetscErrorCode (*)(Mat,PetscViewer))f;
1043: break;
1044: case MATOP_COPY:
1045: shell->ops->copy = (PetscErrorCode (*)(Mat,Mat,MatStructure))f;
1046: break;
1047: case MATOP_DIAGONAL_SET:
1048: case MATOP_DIAGONAL_SCALE:
1049: case MATOP_SHIFT:
1050: case MATOP_SCALE:
1051: case MATOP_AXPY:
1052: case MATOP_ZERO_ROWS:
1053: case MATOP_ZERO_ROWS_COLUMNS:
1054: if (shell->managescalingshifts) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"MATSHELL is managing scalings and shifts, see MatShellSetManageScalingShifts()");
1055: (((void(**)(void))mat->ops)[op]) = f;
1056: break;
1057: case MATOP_GET_DIAGONAL:
1058: if (shell->managescalingshifts) {
1059: shell->ops->getdiagonal = (PetscErrorCode (*)(Mat,Vec))f;
1060: mat->ops->getdiagonal = MatGetDiagonal_Shell;
1061: } else {
1062: shell->ops->getdiagonal = NULL;
1063: mat->ops->getdiagonal = (PetscErrorCode (*)(Mat,Vec))f;
1064: }
1065: break;
1066: case MATOP_MULT:
1067: if (shell->managescalingshifts) {
1068: shell->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1069: mat->ops->mult = MatMult_Shell;
1070: } else {
1071: shell->ops->mult = NULL;
1072: mat->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1073: }
1074: break;
1075: case MATOP_MULT_TRANSPOSE:
1076: if (shell->managescalingshifts) {
1077: shell->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1078: mat->ops->multtranspose = MatMultTranspose_Shell;
1079: } else {
1080: shell->ops->multtranspose = NULL;
1081: mat->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1082: }
1083: break;
1084: default:
1085: (((void(**)(void))mat->ops)[op]) = f;
1086: break;
1087: }
1088: return(0);
1089: }
1091: PetscErrorCode MatShellGetOperation_Shell(Mat mat,MatOperation op,void(**f)(void))
1092: {
1093: Mat_Shell *shell = (Mat_Shell*)mat->data;
1096: switch (op) {
1097: case MATOP_DESTROY:
1098: *f = (void (*)(void))shell->ops->destroy;
1099: break;
1100: case MATOP_VIEW:
1101: *f = (void (*)(void))mat->ops->view;
1102: break;
1103: case MATOP_COPY:
1104: *f = (void (*)(void))shell->ops->copy;
1105: break;
1106: case MATOP_DIAGONAL_SET:
1107: case MATOP_DIAGONAL_SCALE:
1108: case MATOP_SHIFT:
1109: case MATOP_SCALE:
1110: case MATOP_AXPY:
1111: case MATOP_ZERO_ROWS:
1112: case MATOP_ZERO_ROWS_COLUMNS:
1113: *f = (((void (**)(void))mat->ops)[op]);
1114: break;
1115: case MATOP_GET_DIAGONAL:
1116: if (shell->ops->getdiagonal)
1117: *f = (void (*)(void))shell->ops->getdiagonal;
1118: else
1119: *f = (((void (**)(void))mat->ops)[op]);
1120: break;
1121: case MATOP_MULT:
1122: if (shell->ops->mult)
1123: *f = (void (*)(void))shell->ops->mult;
1124: else
1125: *f = (((void (**)(void))mat->ops)[op]);
1126: break;
1127: case MATOP_MULT_TRANSPOSE:
1128: if (shell->ops->multtranspose)
1129: *f = (void (*)(void))shell->ops->multtranspose;
1130: else
1131: *f = (((void (**)(void))mat->ops)[op]);
1132: break;
1133: default:
1134: *f = (((void (**)(void))mat->ops)[op]);
1135: }
1136: return(0);
1137: }
1139: /*MC
1140: MATSHELL - MATSHELL = "shell" - A matrix type to be used to define your own matrix type -- perhaps matrix free.
1142: Level: advanced
1144: .seealso: MatCreateShell()
1145: M*/
1147: PETSC_EXTERN PetscErrorCode MatCreate_Shell(Mat A)
1148: {
1149: Mat_Shell *b;
1153: PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));
1155: PetscNewLog(A,&b);
1156: A->data = (void*)b;
1158: PetscLayoutSetUp(A->rmap);
1159: PetscLayoutSetUp(A->cmap);
1161: b->ctx = 0;
1162: b->vshift = 0.0;
1163: b->vscale = 1.0;
1164: b->managescalingshifts = PETSC_TRUE;
1165: A->assembled = PETSC_TRUE;
1166: A->preallocated = PETSC_FALSE;
1168: PetscObjectComposeFunction((PetscObject)A,"MatShellGetContext_C",MatShellGetContext_Shell);
1169: PetscObjectComposeFunction((PetscObject)A,"MatShellSetContext_C",MatShellSetContext_Shell);
1170: PetscObjectComposeFunction((PetscObject)A,"MatShellSetManageScalingShifts_C",MatShellSetManageScalingShifts_Shell);
1171: PetscObjectComposeFunction((PetscObject)A,"MatShellSetOperation_C",MatShellSetOperation_Shell);
1172: PetscObjectComposeFunction((PetscObject)A,"MatShellGetOperation_C",MatShellGetOperation_Shell);
1173: PetscObjectChangeTypeName((PetscObject)A,MATSHELL);
1174: return(0);
1175: }
1177: /*@C
1178: MatCreateShell - Creates a new matrix class for use with a user-defined
1179: private data storage format.
1181: Collective
1183: Input Parameters:
1184: + comm - MPI communicator
1185: . m - number of local rows (must be given)
1186: . n - number of local columns (must be given)
1187: . M - number of global rows (may be PETSC_DETERMINE)
1188: . N - number of global columns (may be PETSC_DETERMINE)
1189: - ctx - pointer to data needed by the shell matrix routines
1191: Output Parameter:
1192: . A - the matrix
1194: Level: advanced
1196: Usage:
1197: $ extern PetscErrorCode mult(Mat,Vec,Vec);
1198: $ MatCreateShell(comm,m,n,M,N,ctx,&mat);
1199: $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
1200: $ [ Use matrix for operations that have been set ]
1201: $ MatDestroy(mat);
1203: Notes:
1204: The shell matrix type is intended to provide a simple class to use
1205: with KSP (such as, for use with matrix-free methods). You should not
1206: use the shell type if you plan to define a complete matrix class.
1208: Fortran Notes:
1209: To use this from Fortran with a ctx you must write an interface definition for this
1210: function and for MatShellGetContext() that tells Fortran the Fortran derived data type you are passing
1211: in as the ctx argument.
1213: PETSc requires that matrices and vectors being used for certain
1214: operations are partitioned accordingly. For example, when
1215: creating a shell matrix, A, that supports parallel matrix-vector
1216: products using MatMult(A,x,y) the user should set the number
1217: of local matrix rows to be the number of local elements of the
1218: corresponding result vector, y. Note that this is information is
1219: required for use of the matrix interface routines, even though
1220: the shell matrix may not actually be physically partitioned.
1221: For example,
1223: $
1224: $ Vec x, y
1225: $ extern PetscErrorCode mult(Mat,Vec,Vec);
1226: $ Mat A
1227: $
1228: $ VecCreateMPI(comm,PETSC_DECIDE,M,&y);
1229: $ VecCreateMPI(comm,PETSC_DECIDE,N,&x);
1230: $ VecGetLocalSize(y,&m);
1231: $ VecGetLocalSize(x,&n);
1232: $ MatCreateShell(comm,m,n,M,N,ctx,&A);
1233: $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
1234: $ MatMult(A,x,y);
1235: $ MatDestroy(&A);
1236: $ VecDestroy(&y);
1237: $ VecDestroy(&x);
1238: $
1241: MATSHELL handles MatShift(), MatDiagonalSet(), MatDiagonalScale(), MatAXPY(), MatScale(), MatZeroRows() and MatZeroRowsColumns() internally, so these
1242: operations cannot be overwritten unless MatShellSetManageScalingShifts() is called.
1245: For rectangular matrices do all the scalings and shifts make sense?
1247: Developers Notes:
1248: Regarding shifting and scaling. The general form is
1250: diag(left)(vscale*A + diag(dshift) + vshift I)diag(right)
1252: The order you apply the operations is important. For example if you have a dshift then
1253: apply a MatScale(s) you get s*vscale*A + s*diag(shift). But if you first scale and then shift
1254: you get s*vscale*A + diag(shift)
1256: A is the user provided function.
1258: KSP/PC uses changes in the Mat's "state" to decide if preconditioners need to be rebuilt: PCSetUp() only calls the setup() for
1259: for the PC implementation if the Mat state has increased from the previous call. Thus to get changes in a MATSHELL to trigger
1260: an update in the preconditioner you must call MatAssemblyBegin()/MatAssemblyEnd() or PetscObjectStateIncrease((PetscObject)mat);
1261: each time the MATSHELL matrix has changed.
1263: Calling MatAssemblyBegin()/MatAssemblyEnd() on a MATSHELL removes any previously supplied shift and scales that were provided
1264: with MatDiagonalSet(), MatShift(), MatScale(), or MatDiagonalScale().
1266: .seealso: MatShellSetOperation(), MatHasOperation(), MatShellGetContext(), MatShellSetContext(), MATSHELL, MatShellSetManageScalingShifts()
1267: @*/
1268: PetscErrorCode MatCreateShell(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,void *ctx,Mat *A)
1269: {
1273: MatCreate(comm,A);
1274: MatSetSizes(*A,m,n,M,N);
1275: MatSetType(*A,MATSHELL);
1276: MatShellSetContext(*A,ctx);
1277: MatSetUp(*A);
1278: return(0);
1279: }
1282: /*@
1283: MatShellSetContext - sets the context for a shell matrix
1285: Logically Collective on Mat
1287: Input Parameters:
1288: + mat - the shell matrix
1289: - ctx - the context
1291: Level: advanced
1293: Fortran Notes:
1294: To use this from Fortran you must write a Fortran interface definition for this
1295: function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.
1297: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation()
1298: @*/
1299: PetscErrorCode MatShellSetContext(Mat mat,void *ctx)
1300: {
1305: PetscUseMethod(mat,"MatShellSetContext_C",(Mat,void*),(mat,ctx));
1306: return(0);
1307: }
1309: /*@
1310: MatShellSetManageScalingShifts - Allows the user to control the scaling and shift operations of the MATSHELL. Must be called immediately
1311: after MatCreateShell()
1313: Logically Collective on Mat
1315: Input Parameter:
1316: . mat - the shell matrix
1318: Level: advanced
1320: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatShellSetOperation()
1321: @*/
1322: PetscErrorCode MatShellSetManageScalingShifts(Mat A)
1323: {
1328: PetscUseMethod(A,"MatShellSetManageScalingShifts_C",(Mat),(A));
1329: return(0);
1330: }
1332: /*@C
1333: MatShellTestMult - Compares the multiply routine provided to the MATSHELL with differencing on a given function.
1335: Logically Collective on Mat
1337: Input Parameters:
1338: + mat - the shell matrix
1339: . f - the function
1340: . base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated
1341: - ctx - an optional context for the function
1343: Output Parameter:
1344: . flg - PETSC_TRUE if the multiply is likely correct
1346: Options Database:
1347: . -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference
1349: Level: advanced
1351: Fortran Notes:
1352: Not supported from Fortran
1354: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellTestMultTranspose()
1355: @*/
1356: PetscErrorCode MatShellTestMult(Mat mat,PetscErrorCode (*f)(void*,Vec,Vec),Vec base,void *ctx,PetscBool *flg)
1357: {
1359: PetscInt m,n;
1360: Mat mf,Dmf,Dmat,Ddiff;
1361: PetscReal Diffnorm,Dmfnorm;
1362: PetscBool v = PETSC_FALSE, flag = PETSC_TRUE;
1366: PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_view",&v);
1367: MatGetLocalSize(mat,&m,&n);
1368: MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,PETSC_DECIDE,PETSC_DECIDE,&mf);
1369: MatMFFDSetFunction(mf,f,ctx);
1370: MatMFFDSetBase(mf,base,NULL);
1372: MatComputeOperator(mf,MATAIJ,&Dmf);
1373: MatComputeOperator(mat,MATAIJ,&Dmat);
1375: MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff);
1376: MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN);
1377: MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm);
1378: MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm);
1379: if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) {
1380: flag = PETSC_FALSE;
1381: if (v) {
1382: 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));
1383: MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_view");
1384: MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_view");
1385: MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_view");
1386: }
1387: } else if (v) {
1388: PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL and matrix free multiple appear to produce the same results\n");
1389: }
1390: if (flg) *flg = flag;
1391: MatDestroy(&Ddiff);
1392: MatDestroy(&mf);
1393: MatDestroy(&Dmf);
1394: MatDestroy(&Dmat);
1395: return(0);
1396: }
1398: /*@C
1399: MatShellTestMultTranpose - Compares the multiply transpose routine provided to the MATSHELL with differencing on a given function.
1401: Logically Collective on Mat
1403: Input Parameters:
1404: + mat - the shell matrix
1405: . f - the function
1406: . base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated
1407: - ctx - an optional context for the function
1409: Output Parameter:
1410: . flg - PETSC_TRUE if the multiply is likely correct
1412: Options Database:
1413: . -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference
1415: Level: advanced
1417: Fortran Notes:
1418: Not supported from Fortran
1420: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellTestMult()
1421: @*/
1422: PetscErrorCode MatShellTestMultTranspose(Mat mat,PetscErrorCode (*f)(void*,Vec,Vec),Vec base,void *ctx,PetscBool *flg)
1423: {
1425: Vec x,y,z;
1426: PetscInt m,n,M,N;
1427: Mat mf,Dmf,Dmat,Ddiff;
1428: PetscReal Diffnorm,Dmfnorm;
1429: PetscBool v = PETSC_FALSE, flag = PETSC_TRUE;
1433: PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_transpose_view",&v);
1434: MatCreateVecs(mat,&x,&y);
1435: VecDuplicate(y,&z);
1436: MatGetLocalSize(mat,&m,&n);
1437: MatGetSize(mat,&M,&N);
1438: MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,M,N,&mf);
1439: MatMFFDSetFunction(mf,f,ctx);
1440: MatMFFDSetBase(mf,base,NULL);
1441: MatComputeOperator(mf,MATAIJ,&Dmf);
1442: MatTranspose(Dmf,MAT_INPLACE_MATRIX,&Dmf);
1443: MatComputeOperatorTranspose(mat,MATAIJ,&Dmat);
1445: MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff);
1446: MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN);
1447: MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm);
1448: MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm);
1449: if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) {
1450: flag = PETSC_FALSE;
1451: if (v) {
1452: 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));
1453: MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");
1454: MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");
1455: MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");
1456: }
1457: } else if (v) {
1458: PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL transpose and matrix free multiple appear to produce the same results\n");
1459: }
1460: if (flg) *flg = flag;
1461: MatDestroy(&mf);
1462: MatDestroy(&Dmat);
1463: MatDestroy(&Ddiff);
1464: MatDestroy(&Dmf);
1465: VecDestroy(&x);
1466: VecDestroy(&y);
1467: VecDestroy(&z);
1468: return(0);
1469: }
1471: /*@C
1472: MatShellSetOperation - Allows user to set a matrix operation for a shell matrix.
1474: Logically Collective on Mat
1476: Input Parameters:
1477: + mat - the shell matrix
1478: . op - the name of the operation
1479: - g - the function that provides the operation.
1481: Level: advanced
1483: Usage:
1484: $ extern PetscErrorCode usermult(Mat,Vec,Vec);
1485: $ MatCreateShell(comm,m,n,M,N,ctx,&A);
1486: $ MatShellSetOperation(A,MATOP_MULT,(void(*)(void))usermult);
1488: Notes:
1489: See the file include/petscmat.h for a complete list of matrix
1490: operations, which all have the form MATOP_<OPERATION>, where
1491: <OPERATION> is the name (in all capital letters) of the
1492: user interface routine (e.g., MatMult() -> MATOP_MULT).
1494: All user-provided functions (except for MATOP_DESTROY) should have the same calling
1495: sequence as the usual matrix interface routines, since they
1496: are intended to be accessed via the usual matrix interface
1497: routines, e.g.,
1498: $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec)
1500: In particular each function MUST return an error code of 0 on success and
1501: nonzero on failure.
1503: Within each user-defined routine, the user should call
1504: MatShellGetContext() to obtain the user-defined context that was
1505: set by MatCreateShell().
1507: Fortran Notes:
1508: For MatCreateVecs() the user code should check if the input left or right matrix is -1 and in that case not
1509: generate a matrix. See src/mat/examples/tests/ex120f.F
1511: Use MatSetOperation() to set an operation for any matrix type
1513: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatSetOperation(), MatShellSetManageScalingShifts()
1514: @*/
1515: PetscErrorCode MatShellSetOperation(Mat mat,MatOperation op,void (*g)(void))
1516: {
1521: PetscUseMethod(mat,"MatShellSetOperation_C",(Mat,MatOperation,void (*)(void)),(mat,op,g));
1522: return(0);
1523: }
1525: /*@C
1526: MatShellGetOperation - Gets a matrix function for a shell matrix.
1528: Not Collective
1530: Input Parameters:
1531: + mat - the shell matrix
1532: - op - the name of the operation
1534: Output Parameter:
1535: . g - the function that provides the operation.
1537: Level: advanced
1539: Notes:
1540: See the file include/petscmat.h for a complete list of matrix
1541: operations, which all have the form MATOP_<OPERATION>, where
1542: <OPERATION> is the name (in all capital letters) of the
1543: user interface routine (e.g., MatMult() -> MATOP_MULT).
1545: All user-provided functions have the same calling
1546: sequence as the usual matrix interface routines, since they
1547: are intended to be accessed via the usual matrix interface
1548: routines, e.g.,
1549: $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec)
1551: Within each user-defined routine, the user should call
1552: MatShellGetContext() to obtain the user-defined context that was
1553: set by MatCreateShell().
1555: .seealso: MatCreateShell(), MatShellGetContext(), MatShellSetOperation(), MatShellSetContext()
1556: @*/
1557: PetscErrorCode MatShellGetOperation(Mat mat,MatOperation op,void(**g)(void))
1558: {
1563: PetscUseMethod(mat,"MatShellGetOperation_C",(Mat,MatOperation,void (**)(void)),(mat,op,g));
1564: return(0);
1565: }