Actual source code: shell.c
petsc-3.13.6 2020-09-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: /* support for MatAXPY */
28: Mat axpy;
29: PetscScalar axpy_vscale;
30: PetscObjectState axpy_state;
31: PetscBool managescalingshifts; /* The user will manage the scaling and shifts for the MATSHELL, not the default */
32: /* support for ZeroRows/Columns operations */
33: IS zrows;
34: IS zcols;
35: Vec zvals;
36: Vec zvals_w;
37: VecScatter zvals_sct_r;
38: VecScatter zvals_sct_c;
39: void *ctx;
40: } Mat_Shell;
43: /*
44: Store and scale values on zeroed rows
45: xx = [x_1, 0], 0 on zeroed columns
46: */
47: static PetscErrorCode MatShellPreZeroRight(Mat A,Vec x,Vec *xx)
48: {
49: Mat_Shell *shell = (Mat_Shell*)A->data;
53: *xx = x;
54: if (shell->zrows) {
55: VecSet(shell->zvals_w,0.0);
56: VecScatterBegin(shell->zvals_sct_c,x,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);
57: VecScatterEnd(shell->zvals_sct_c,x,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);
58: VecPointwiseMult(shell->zvals_w,shell->zvals_w,shell->zvals);
59: }
60: if (shell->zcols) {
61: if (!shell->right_work) {
62: MatCreateVecs(A,&shell->right_work,NULL);
63: }
64: VecCopy(x,shell->right_work);
65: VecISSet(shell->right_work,shell->zcols,0.0);
66: *xx = shell->right_work;
67: }
68: return(0);
69: }
71: /* Insert properly diagonally scaled values stored in MatShellPreZeroRight */
72: static PetscErrorCode MatShellPostZeroLeft(Mat A,Vec x)
73: {
74: Mat_Shell *shell = (Mat_Shell*)A->data;
78: if (shell->zrows) {
79: VecScatterBegin(shell->zvals_sct_r,shell->zvals_w,x,INSERT_VALUES,SCATTER_REVERSE);
80: VecScatterEnd(shell->zvals_sct_r,shell->zvals_w,x,INSERT_VALUES,SCATTER_REVERSE);
81: }
82: return(0);
83: }
85: /*
86: Store and scale values on zeroed rows
87: xx = [x_1, 0], 0 on zeroed rows
88: */
89: static PetscErrorCode MatShellPreZeroLeft(Mat A,Vec x,Vec *xx)
90: {
91: Mat_Shell *shell = (Mat_Shell*)A->data;
95: *xx = NULL;
96: if (!shell->zrows) {
97: *xx = x;
98: } else {
99: if (!shell->left_work) {
100: MatCreateVecs(A,NULL,&shell->left_work);
101: }
102: VecCopy(x,shell->left_work);
103: VecSet(shell->zvals_w,0.0);
104: VecScatterBegin(shell->zvals_sct_r,shell->zvals_w,shell->left_work,INSERT_VALUES,SCATTER_REVERSE);
105: VecScatterEnd(shell->zvals_sct_r,shell->zvals_w,shell->left_work,INSERT_VALUES,SCATTER_REVERSE);
106: VecScatterBegin(shell->zvals_sct_r,x,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);
107: VecScatterEnd(shell->zvals_sct_r,x,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);
108: VecPointwiseMult(shell->zvals_w,shell->zvals_w,shell->zvals);
109: *xx = shell->left_work;
110: }
111: return(0);
112: }
114: /* Zero zero-columns contributions, sum contributions from properly scaled values stored in MatShellPreZeroLeft */
115: static PetscErrorCode MatShellPostZeroRight(Mat A,Vec x)
116: {
117: Mat_Shell *shell = (Mat_Shell*)A->data;
121: if (shell->zcols) {
122: VecISSet(x,shell->zcols,0.0);
123: }
124: if (shell->zrows) {
125: VecScatterBegin(shell->zvals_sct_c,shell->zvals_w,x,ADD_VALUES,SCATTER_REVERSE);
126: VecScatterEnd(shell->zvals_sct_c,shell->zvals_w,x,ADD_VALUES,SCATTER_REVERSE);
127: }
128: return(0);
129: }
131: /*
132: xx = diag(left)*x
133: */
134: static PetscErrorCode MatShellPreScaleLeft(Mat A,Vec x,Vec *xx)
135: {
136: Mat_Shell *shell = (Mat_Shell*)A->data;
140: *xx = NULL;
141: if (!shell->left) {
142: *xx = x;
143: } else {
144: if (!shell->left_work) {VecDuplicate(shell->left,&shell->left_work);}
145: VecPointwiseMult(shell->left_work,x,shell->left);
146: *xx = shell->left_work;
147: }
148: return(0);
149: }
151: /*
152: xx = diag(right)*x
153: */
154: static PetscErrorCode MatShellPreScaleRight(Mat A,Vec x,Vec *xx)
155: {
156: Mat_Shell *shell = (Mat_Shell*)A->data;
160: *xx = NULL;
161: if (!shell->right) {
162: *xx = x;
163: } else {
164: if (!shell->right_work) {VecDuplicate(shell->right,&shell->right_work);}
165: VecPointwiseMult(shell->right_work,x,shell->right);
166: *xx = shell->right_work;
167: }
168: return(0);
169: }
171: /*
172: x = diag(left)*x
173: */
174: static PetscErrorCode MatShellPostScaleLeft(Mat A,Vec x)
175: {
176: Mat_Shell *shell = (Mat_Shell*)A->data;
180: if (shell->left) {VecPointwiseMult(x,x,shell->left);}
181: return(0);
182: }
184: /*
185: x = diag(right)*x
186: */
187: static PetscErrorCode MatShellPostScaleRight(Mat A,Vec x)
188: {
189: Mat_Shell *shell = (Mat_Shell*)A->data;
193: if (shell->right) {VecPointwiseMult(x,x,shell->right);}
194: return(0);
195: }
197: /*
198: Y = vscale*Y + diag(dshift)*X + vshift*X
200: On input Y already contains A*x
201: */
202: static PetscErrorCode MatShellShiftAndScale(Mat A,Vec X,Vec Y)
203: {
204: Mat_Shell *shell = (Mat_Shell*)A->data;
208: if (shell->dshift) { /* get arrays because there is no VecPointwiseMultAdd() */
209: PetscInt i,m;
210: const PetscScalar *x,*d;
211: PetscScalar *y;
212: VecGetLocalSize(X,&m);
213: VecGetArrayRead(shell->dshift,&d);
214: VecGetArrayRead(X,&x);
215: VecGetArray(Y,&y);
216: for (i=0; i<m; i++) y[i] = shell->vscale*y[i] + d[i]*x[i];
217: VecRestoreArrayRead(shell->dshift,&d);
218: VecRestoreArrayRead(X,&x);
219: VecRestoreArray(Y,&y);
220: } else {
221: VecScale(Y,shell->vscale);
222: }
223: if (shell->vshift != 0.0) {VecAXPY(Y,shell->vshift,X);} /* if test is for non-square matrices */
224: return(0);
225: }
227: PetscErrorCode MatShellGetContext_Shell(Mat mat,void *ctx)
228: {
230: *(void**)ctx = ((Mat_Shell*)(mat->data))->ctx;
231: return(0);
232: }
234: /*@
235: MatShellGetContext - Returns the user-provided context associated with a shell matrix.
237: Not Collective
239: Input Parameter:
240: . mat - the matrix, should have been created with MatCreateShell()
242: Output Parameter:
243: . ctx - the user provided context
245: Level: advanced
247: Fortran Notes:
248: To use this from Fortran you must write a Fortran interface definition for this
249: function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.
251: .seealso: MatCreateShell(), MatShellSetOperation(), MatShellSetContext()
252: @*/
253: PetscErrorCode MatShellGetContext(Mat mat,void *ctx)
254: {
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: {
267: Mat_Shell *shell = (Mat_Shell*)mat->data;
268: Vec x = NULL,b = NULL;
269: IS is1, is2;
270: const PetscInt *ridxs;
271: PetscInt *idxs,*gidxs;
272: PetscInt cum,rst,cst,i;
275: if (!shell->zvals) {
276: MatCreateVecs(mat,NULL,&shell->zvals);
277: }
278: if (!shell->zvals_w) {
279: VecDuplicate(shell->zvals,&shell->zvals_w);
280: }
281: MatGetOwnershipRange(mat,&rst,NULL);
282: MatGetOwnershipRangeColumn(mat,&cst,NULL);
284: /* Expand/create index set of zeroed rows */
285: PetscMalloc1(nr,&idxs);
286: for (i = 0; i < nr; i++) idxs[i] = rows[i] + rst;
287: ISCreateGeneral(PETSC_COMM_SELF,nr,idxs,PETSC_OWN_POINTER,&is1);
288: ISSort(is1);
289: VecISSet(shell->zvals,is1,diag);
290: if (shell->zrows) {
291: ISSum(shell->zrows,is1,&is2);
292: ISDestroy(&shell->zrows);
293: ISDestroy(&is1);
294: shell->zrows = is2;
295: } else shell->zrows = is1;
297: /* Create scatters for diagonal values communications */
298: VecScatterDestroy(&shell->zvals_sct_c);
299: VecScatterDestroy(&shell->zvals_sct_r);
301: /* row scatter: from/to left vector */
302: MatCreateVecs(mat,&x,&b);
303: VecScatterCreate(b,shell->zrows,shell->zvals_w,shell->zrows,&shell->zvals_sct_r);
305: /* col scatter: from right vector to left vector */
306: ISGetIndices(shell->zrows,&ridxs);
307: ISGetLocalSize(shell->zrows,&nr);
308: PetscMalloc1(nr,&gidxs);
309: for (i = 0, cum = 0; i < nr; i++) {
310: if (ridxs[i] >= mat->cmap->N) continue;
311: gidxs[cum] = ridxs[i];
312: cum++;
313: }
314: ISCreateGeneral(PetscObjectComm((PetscObject)mat),cum,gidxs,PETSC_OWN_POINTER,&is1);
315: VecScatterCreate(x,is1,shell->zvals_w,is1,&shell->zvals_sct_c);
316: ISDestroy(&is1);
317: VecDestroy(&x);
318: VecDestroy(&b);
320: /* Expand/create index set of zeroed columns */
321: if (rc) {
322: PetscMalloc1(nc,&idxs);
323: for (i = 0; i < nc; i++) idxs[i] = cols[i] + cst;
324: ISCreateGeneral(PETSC_COMM_SELF,nc,idxs,PETSC_OWN_POINTER,&is1);
325: ISSort(is1);
326: if (shell->zcols) {
327: ISSum(shell->zcols,is1,&is2);
328: ISDestroy(&shell->zcols);
329: ISDestroy(&is1);
330: shell->zcols = is2;
331: } else shell->zcols = is1;
332: }
333: return(0);
334: }
336: static PetscErrorCode MatZeroRows_Shell(Mat mat,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
337: {
338: Mat_Shell *shell = (Mat_Shell*)mat->data;
339: PetscInt nr, *lrows;
343: if (x && b) {
344: Vec xt;
345: PetscScalar *vals;
346: PetscInt *gcols,i,st,nl,nc;
348: PetscMalloc1(n,&gcols);
349: for (i = 0, nc = 0; i < n; i++) if (rows[i] < mat->cmap->N) gcols[nc++] = rows[i];
351: MatCreateVecs(mat,&xt,NULL);
352: VecCopy(x,xt);
353: PetscCalloc1(nc,&vals);
354: VecSetValues(xt,nc,gcols,vals,INSERT_VALUES); /* xt = [x1, 0] */
355: PetscFree(vals);
356: VecAssemblyBegin(xt);
357: VecAssemblyEnd(xt);
358: VecAYPX(xt,-1.0,x); /* xt = [0, x2] */
360: VecGetOwnershipRange(xt,&st,NULL);
361: VecGetLocalSize(xt,&nl);
362: VecGetArray(xt,&vals);
363: for (i = 0; i < nl; i++) {
364: PetscInt g = i + st;
365: if (g > mat->rmap->N) continue;
366: if (PetscAbsScalar(vals[i]) == 0.0) continue;
367: VecSetValue(b,g,diag*vals[i],INSERT_VALUES);
368: }
369: VecRestoreArray(xt,&vals);
370: VecAssemblyBegin(b);
371: VecAssemblyEnd(b); /* b = [b1, x2 * diag] */
372: VecDestroy(&xt);
373: PetscFree(gcols);
374: }
375: PetscLayoutMapLocal(mat->rmap,n,rows,&nr,&lrows,NULL);
376: MatZeroRowsColumns_Local_Shell(mat,nr,lrows,0,NULL,diag,PETSC_FALSE);
377: if (shell->axpy) {
378: MatZeroRows(shell->axpy,n,rows,0.0,NULL,NULL);
379: }
380: PetscFree(lrows);
381: return(0);
382: }
384: static PetscErrorCode MatZeroRowsColumns_Shell(Mat mat,PetscInt n,const PetscInt rowscols[],PetscScalar diag,Vec x,Vec b)
385: {
386: Mat_Shell *shell = (Mat_Shell*)mat->data;
387: PetscInt *lrows, *lcols;
388: PetscInt nr, nc;
389: PetscBool congruent;
393: if (x && b) {
394: Vec xt, bt;
395: PetscScalar *vals;
396: PetscInt *grows,*gcols,i,st,nl;
398: PetscMalloc2(n,&grows,n,&gcols);
399: for (i = 0, nr = 0; i < n; i++) if (rowscols[i] < mat->rmap->N) grows[nr++] = rowscols[i];
400: for (i = 0, nc = 0; i < n; i++) if (rowscols[i] < mat->cmap->N) gcols[nc++] = rowscols[i];
401: PetscCalloc1(n,&vals);
403: MatCreateVecs(mat,&xt,&bt);
404: VecCopy(x,xt);
405: VecSetValues(xt,nc,gcols,vals,INSERT_VALUES); /* xt = [x1, 0] */
406: VecAssemblyBegin(xt);
407: VecAssemblyEnd(xt);
408: VecAXPY(xt,-1.0,x); /* xt = [0, -x2] */
409: MatMult(mat,xt,bt); /* bt = [-A12*x2,-A22*x2] */
410: VecSetValues(bt,nr,grows,vals,INSERT_VALUES); /* bt = [-A12*x2,0] */
411: VecAssemblyBegin(bt);
412: VecAssemblyEnd(bt);
413: VecAXPY(b,1.0,bt); /* b = [b1 - A12*x2, b2] */
414: VecSetValues(bt,nr,grows,vals,INSERT_VALUES); /* b = [b1 - A12*x2, 0] */
415: VecAssemblyBegin(bt);
416: VecAssemblyEnd(bt);
417: PetscFree(vals);
419: VecGetOwnershipRange(xt,&st,NULL);
420: VecGetLocalSize(xt,&nl);
421: VecGetArray(xt,&vals);
422: for (i = 0; i < nl; i++) {
423: PetscInt g = i + st;
424: if (g > mat->rmap->N) continue;
425: if (PetscAbsScalar(vals[i]) == 0.0) continue;
426: VecSetValue(b,g,-diag*vals[i],INSERT_VALUES);
427: }
428: VecRestoreArray(xt,&vals);
429: VecAssemblyBegin(b);
430: VecAssemblyEnd(b); /* b = [b1 - A12*x2, x2 * diag] */
431: VecDestroy(&xt);
432: VecDestroy(&bt);
433: PetscFree2(grows,gcols);
434: }
435: PetscLayoutMapLocal(mat->rmap,n,rowscols,&nr,&lrows,NULL);
436: MatHasCongruentLayouts(mat,&congruent);
437: if (congruent) {
438: nc = nr;
439: lcols = lrows;
440: } else { /* MatZeroRowsColumns implicitly assumes the rowscols indices are for a square matrix, here we handle a more general case */
441: PetscInt i,nt,*t;
443: PetscMalloc1(n,&t);
444: for (i = 0, nt = 0; i < n; i++) if (rowscols[i] < mat->cmap->N) t[nt++] = rowscols[i];
445: PetscLayoutMapLocal(mat->cmap,nt,t,&nc,&lcols,NULL);
446: PetscFree(t);
447: }
448: MatZeroRowsColumns_Local_Shell(mat,nr,lrows,nc,lcols,diag,PETSC_TRUE);
449: if (!congruent) {
450: PetscFree(lcols);
451: }
452: PetscFree(lrows);
453: if (shell->axpy) {
454: MatZeroRowsColumns(shell->axpy,n,rowscols,0.0,NULL,NULL);
455: }
456: return(0);
457: }
459: PetscErrorCode MatDestroy_Shell(Mat mat)
460: {
462: Mat_Shell *shell = (Mat_Shell*)mat->data;
465: if (shell->ops->destroy) {
466: (*shell->ops->destroy)(mat);
467: }
468: PetscMemzero(shell->ops,sizeof(struct _MatShellOps));
469: VecDestroy(&shell->left);
470: VecDestroy(&shell->right);
471: VecDestroy(&shell->dshift);
472: VecDestroy(&shell->left_work);
473: VecDestroy(&shell->right_work);
474: VecDestroy(&shell->left_add_work);
475: VecDestroy(&shell->right_add_work);
476: MatDestroy(&shell->axpy);
477: VecDestroy(&shell->zvals_w);
478: VecDestroy(&shell->zvals);
479: VecScatterDestroy(&shell->zvals_sct_c);
480: VecScatterDestroy(&shell->zvals_sct_r);
481: ISDestroy(&shell->zrows);
482: ISDestroy(&shell->zcols);
483: PetscFree(mat->data);
484: PetscObjectComposeFunction((PetscObject)mat,"MatShellGetContext_C",NULL);
485: PetscObjectComposeFunction((PetscObject)mat,"MatShellSetContext_C",NULL);
486: PetscObjectComposeFunction((PetscObject)mat,"MatShellSetVecType_C",NULL);
487: PetscObjectComposeFunction((PetscObject)mat,"MatShellSetManageScalingShifts_C",NULL);
488: PetscObjectComposeFunction((PetscObject)mat,"MatShellSetOperation_C",NULL);
489: PetscObjectComposeFunction((PetscObject)mat,"MatShellGetOperation_C",NULL);
490: return(0);
491: }
493: PetscErrorCode MatCopy_Shell(Mat A,Mat B,MatStructure str)
494: {
495: Mat_Shell *shellA = (Mat_Shell*)A->data,*shellB = (Mat_Shell*)B->data;
496: PetscErrorCode ierr;
497: PetscBool matflg;
500: PetscObjectTypeCompare((PetscObject)B,MATSHELL,&matflg);
501: if (!matflg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_NOTSAMETYPE,"Output matrix must be a MATSHELL");
503: PetscMemcpy(B->ops,A->ops,sizeof(struct _MatOps));
504: PetscMemcpy(shellB->ops,shellA->ops,sizeof(struct _MatShellOps));
506: if (shellA->ops->copy) {
507: (*shellA->ops->copy)(A,B,str);
508: }
509: shellB->vscale = shellA->vscale;
510: shellB->vshift = shellA->vshift;
511: if (shellA->dshift) {
512: if (!shellB->dshift) {
513: VecDuplicate(shellA->dshift,&shellB->dshift);
514: }
515: VecCopy(shellA->dshift,shellB->dshift);
516: } else {
517: VecDestroy(&shellB->dshift);
518: }
519: if (shellA->left) {
520: if (!shellB->left) {
521: VecDuplicate(shellA->left,&shellB->left);
522: }
523: VecCopy(shellA->left,shellB->left);
524: } else {
525: VecDestroy(&shellB->left);
526: }
527: if (shellA->right) {
528: if (!shellB->right) {
529: VecDuplicate(shellA->right,&shellB->right);
530: }
531: VecCopy(shellA->right,shellB->right);
532: } else {
533: VecDestroy(&shellB->right);
534: }
535: MatDestroy(&shellB->axpy);
536: shellB->axpy_vscale = 0.0;
537: shellB->axpy_state = 0;
538: if (shellA->axpy) {
539: PetscObjectReference((PetscObject)shellA->axpy);
540: shellB->axpy = shellA->axpy;
541: shellB->axpy_vscale = shellA->axpy_vscale;
542: shellB->axpy_state = shellA->axpy_state;
543: }
544: if (shellA->zrows) {
545: ISDuplicate(shellA->zrows,&shellB->zrows);
546: if (shellA->zcols) {
547: ISDuplicate(shellA->zcols,&shellB->zcols);
548: }
549: VecDuplicate(shellA->zvals,&shellB->zvals);
550: VecCopy(shellA->zvals,shellB->zvals);
551: VecDuplicate(shellA->zvals_w,&shellB->zvals_w);
552: PetscObjectReference((PetscObject)shellA->zvals_sct_r);
553: PetscObjectReference((PetscObject)shellA->zvals_sct_c);
554: shellB->zvals_sct_r = shellA->zvals_sct_r;
555: shellB->zvals_sct_c = shellA->zvals_sct_c;
556: }
557: return(0);
558: }
560: PetscErrorCode MatDuplicate_Shell(Mat mat,MatDuplicateOption op,Mat *M)
561: {
563: void *ctx;
566: MatShellGetContext(mat,&ctx);
567: MatCreateShell(PetscObjectComm((PetscObject)mat),mat->rmap->n,mat->cmap->n,mat->rmap->N,mat->cmap->N,ctx,M);
568: if (op != MAT_DO_NOT_COPY_VALUES) {
569: MatCopy(mat,*M,SAME_NONZERO_PATTERN);
570: }
571: return(0);
572: }
574: PetscErrorCode MatMult_Shell(Mat A,Vec x,Vec y)
575: {
576: Mat_Shell *shell = (Mat_Shell*)A->data;
577: PetscErrorCode ierr;
578: Vec xx;
579: PetscObjectState instate,outstate;
582: if (!shell->ops->mult) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Have not provided a MatMult() for this MATSHELL");
583: MatShellPreZeroRight(A,x,&xx);
584: MatShellPreScaleRight(A,xx,&xx);
585: PetscObjectStateGet((PetscObject)y, &instate);
586: (*shell->ops->mult)(A,xx,y);
587: PetscObjectStateGet((PetscObject)y, &outstate);
588: if (instate == outstate) {
589: /* increase the state of the output vector since the user did not update its state themself as should have been done */
590: PetscObjectStateIncrease((PetscObject)y);
591: }
592: MatShellShiftAndScale(A,xx,y);
593: MatShellPostScaleLeft(A,y);
594: MatShellPostZeroLeft(A,y);
596: if (shell->axpy) {
597: Mat X;
598: PetscObjectState axpy_state;
600: MatShellGetContext(shell->axpy,(void *)&X);
601: PetscObjectStateGet((PetscObject)X,&axpy_state);
602: if (shell->axpy_state != axpy_state) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Invalid AXPY state: cannot modify the X matrix passed to MatAXPY(Y,a,X,...)");
603: if (!shell->left_work) {MatCreateVecs(A,&shell->left_work,NULL);}
604: MatMult(shell->axpy,x,shell->left_work);
605: VecAXPY(y,shell->axpy_vscale,shell->left_work);
606: }
607: return(0);
608: }
610: PetscErrorCode MatMultAdd_Shell(Mat A,Vec x,Vec y,Vec z)
611: {
612: Mat_Shell *shell = (Mat_Shell*)A->data;
616: if (y == z) {
617: if (!shell->right_add_work) {VecDuplicate(z,&shell->right_add_work);}
618: MatMult(A,x,shell->right_add_work);
619: VecAXPY(z,1.0,shell->right_add_work);
620: } else {
621: MatMult(A,x,z);
622: VecAXPY(z,1.0,y);
623: }
624: return(0);
625: }
627: PetscErrorCode MatMultTranspose_Shell(Mat A,Vec x,Vec y)
628: {
629: Mat_Shell *shell = (Mat_Shell*)A->data;
630: PetscErrorCode ierr;
631: Vec xx;
632: PetscObjectState instate,outstate;
635: MatShellPreZeroLeft(A,x,&xx);
636: MatShellPreScaleLeft(A,xx,&xx);
637: PetscObjectStateGet((PetscObject)y, &instate);
638: if (!shell->ops->multtranspose) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Have not provided a MatMultTranspose() for this MATSHELL");
639: (*shell->ops->multtranspose)(A,xx,y);
640: PetscObjectStateGet((PetscObject)y, &outstate);
641: if (instate == outstate) {
642: /* increase the state of the output vector since the user did not update its state themself as should have been done */
643: PetscObjectStateIncrease((PetscObject)y);
644: }
645: MatShellShiftAndScale(A,xx,y);
646: MatShellPostScaleRight(A,y);
647: MatShellPostZeroRight(A,y);
649: if (shell->axpy) {
650: Mat X;
651: PetscObjectState axpy_state;
653: MatShellGetContext(shell->axpy,(void *)&X);
654: PetscObjectStateGet((PetscObject)X,&axpy_state);
655: if (shell->axpy_state != axpy_state) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Invalid AXPY state: cannot modify the X matrix passed to MatAXPY(Y,a,X,...)");
656: if (!shell->right_work) {MatCreateVecs(A,NULL,&shell->right_work);}
657: MatMultTranspose(shell->axpy,x,shell->right_work);
658: VecAXPY(y,shell->axpy_vscale,shell->right_work);
659: }
660: return(0);
661: }
663: PetscErrorCode MatMultTransposeAdd_Shell(Mat A,Vec x,Vec y,Vec z)
664: {
665: Mat_Shell *shell = (Mat_Shell*)A->data;
669: if (y == z) {
670: if (!shell->left_add_work) {VecDuplicate(z,&shell->left_add_work);}
671: MatMultTranspose(A,x,shell->left_add_work);
672: VecAXPY(z,1.0,shell->left_add_work);
673: } else {
674: MatMultTranspose(A,x,z);
675: VecAXPY(z,1.0,y);
676: }
677: return(0);
678: }
680: /*
681: diag(left)(vscale*A + diag(dshift) + vshift I)diag(right)
682: */
683: PetscErrorCode MatGetDiagonal_Shell(Mat A,Vec v)
684: {
685: Mat_Shell *shell = (Mat_Shell*)A->data;
689: if (shell->ops->getdiagonal) {
690: (*shell->ops->getdiagonal)(A,v);
691: } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Must provide shell matrix with routine to return diagonal using\nMatShellSetOperation(S,MATOP_GET_DIAGONAL,...)");
692: VecScale(v,shell->vscale);
693: if (shell->dshift) {
694: VecAXPY(v,1.0,shell->dshift);
695: }
696: VecShift(v,shell->vshift);
697: if (shell->left) {VecPointwiseMult(v,v,shell->left);}
698: if (shell->right) {VecPointwiseMult(v,v,shell->right);}
699: if (shell->zrows) {
700: VecScatterBegin(shell->zvals_sct_r,shell->zvals,v,INSERT_VALUES,SCATTER_REVERSE);
701: VecScatterEnd(shell->zvals_sct_r,shell->zvals,v,INSERT_VALUES,SCATTER_REVERSE);
702: }
703: if (shell->axpy) {
704: Mat X;
705: PetscObjectState axpy_state;
707: MatShellGetContext(shell->axpy,(void *)&X);
708: PetscObjectStateGet((PetscObject)X,&axpy_state);
709: if (shell->axpy_state != axpy_state) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Invalid AXPY state: cannot modify the X matrix passed to MatAXPY(Y,a,X,...)");
710: if (!shell->left_work) {VecDuplicate(v,&shell->left_work);}
711: MatGetDiagonal(shell->axpy,shell->left_work);
712: VecAXPY(v,shell->axpy_vscale,shell->left_work);
713: }
714: return(0);
715: }
717: PetscErrorCode MatShift_Shell(Mat Y,PetscScalar a)
718: {
719: Mat_Shell *shell = (Mat_Shell*)Y->data;
721: PetscBool flg;
724: MatHasCongruentLayouts(Y,&flg);
725: if (!flg) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_SUP,"Cannot shift shell matrix if it is not congruent");
726: if (shell->left || shell->right) {
727: if (!shell->dshift) {
728: VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift);
729: VecSet(shell->dshift,a);
730: } else {
731: if (shell->left) {VecPointwiseMult(shell->dshift,shell->dshift,shell->left);}
732: if (shell->right) {VecPointwiseMult(shell->dshift,shell->dshift,shell->right);}
733: VecShift(shell->dshift,a);
734: }
735: if (shell->left) {VecPointwiseDivide(shell->dshift,shell->dshift,shell->left);}
736: if (shell->right) {VecPointwiseDivide(shell->dshift,shell->dshift,shell->right);}
737: } else shell->vshift += a;
738: if (shell->zrows) {
739: VecShift(shell->zvals,a);
740: }
741: return(0);
742: }
744: PetscErrorCode MatDiagonalSet_Shell_Private(Mat A,Vec D,PetscScalar s)
745: {
746: Mat_Shell *shell = (Mat_Shell*)A->data;
750: if (!shell->dshift) {VecDuplicate(D,&shell->dshift);}
751: if (shell->left || shell->right) {
752: if (!shell->right_work) {VecDuplicate(shell->left ? shell->left : shell->right, &shell->right_work);}
753: if (shell->left && shell->right) {
754: VecPointwiseDivide(shell->right_work,D,shell->left);
755: VecPointwiseDivide(shell->right_work,shell->right_work,shell->right);
756: } else if (shell->left) {
757: VecPointwiseDivide(shell->right_work,D,shell->left);
758: } else {
759: VecPointwiseDivide(shell->right_work,D,shell->right);
760: }
761: VecAXPY(shell->dshift,s,shell->right_work);
762: } else {
763: VecAXPY(shell->dshift,s,D);
764: }
765: return(0);
766: }
768: PetscErrorCode MatDiagonalSet_Shell(Mat A,Vec D,InsertMode ins)
769: {
770: Mat_Shell *shell = (Mat_Shell*)A->data;
771: Vec d;
773: PetscBool flg;
776: MatHasCongruentLayouts(A,&flg);
777: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Cannot diagonal set or shift shell matrix if it is not congruent");
778: if (ins == INSERT_VALUES) {
779: if (!A->ops->getdiagonal) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Operation MATOP_GETDIAGONAL must be set first");
780: VecDuplicate(D,&d);
781: MatGetDiagonal(A,d);
782: MatDiagonalSet_Shell_Private(A,d,-1.);
783: MatDiagonalSet_Shell_Private(A,D,1.);
784: VecDestroy(&d);
785: if (shell->zrows) {
786: VecCopy(D,shell->zvals);
787: }
788: } else {
789: MatDiagonalSet_Shell_Private(A,D,1.);
790: if (shell->zrows) {
791: VecAXPY(shell->zvals,1.0,D);
792: }
793: }
794: return(0);
795: }
797: PetscErrorCode MatScale_Shell(Mat Y,PetscScalar a)
798: {
799: Mat_Shell *shell = (Mat_Shell*)Y->data;
803: shell->vscale *= a;
804: shell->vshift *= a;
805: if (shell->dshift) {
806: VecScale(shell->dshift,a);
807: }
808: shell->axpy_vscale *= a;
809: if (shell->zrows) {
810: VecScale(shell->zvals,a);
811: }
812: return(0);
813: }
815: static PetscErrorCode MatDiagonalScale_Shell(Mat Y,Vec left,Vec right)
816: {
817: Mat_Shell *shell = (Mat_Shell*)Y->data;
821: if (left) {
822: if (!shell->left) {
823: VecDuplicate(left,&shell->left);
824: VecCopy(left,shell->left);
825: } else {
826: VecPointwiseMult(shell->left,shell->left,left);
827: }
828: if (shell->zrows) {
829: VecPointwiseMult(shell->zvals,shell->zvals,left);
830: }
831: }
832: if (right) {
833: if (!shell->right) {
834: VecDuplicate(right,&shell->right);
835: VecCopy(right,shell->right);
836: } else {
837: VecPointwiseMult(shell->right,shell->right,right);
838: }
839: if (shell->zrows) {
840: if (!shell->left_work) {
841: MatCreateVecs(Y,NULL,&shell->left_work);
842: }
843: VecSet(shell->zvals_w,1.0);
844: VecScatterBegin(shell->zvals_sct_c,right,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);
845: VecScatterEnd(shell->zvals_sct_c,right,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);
846: VecPointwiseMult(shell->zvals,shell->zvals,shell->zvals_w);
847: }
848: }
849: if (shell->axpy) {
850: MatDiagonalScale(shell->axpy,left,right);
851: }
852: return(0);
853: }
855: PetscErrorCode MatAssemblyEnd_Shell(Mat Y,MatAssemblyType t)
856: {
857: Mat_Shell *shell = (Mat_Shell*)Y->data;
861: if (t == MAT_FINAL_ASSEMBLY) {
862: shell->vshift = 0.0;
863: shell->vscale = 1.0;
864: shell->axpy_vscale = 0.0;
865: shell->axpy_state = 0;
866: VecDestroy(&shell->dshift);
867: VecDestroy(&shell->left);
868: VecDestroy(&shell->right);
869: MatDestroy(&shell->axpy);
870: VecScatterDestroy(&shell->zvals_sct_c);
871: VecScatterDestroy(&shell->zvals_sct_r);
872: ISDestroy(&shell->zrows);
873: ISDestroy(&shell->zcols);
874: }
875: return(0);
876: }
878: static PetscErrorCode MatMissingDiagonal_Shell(Mat A,PetscBool *missing,PetscInt *d)
879: {
881: *missing = PETSC_FALSE;
882: return(0);
883: }
885: PetscErrorCode MatAXPY_Shell(Mat Y,PetscScalar a,Mat X,MatStructure str)
886: {
887: Mat_Shell *shell = (Mat_Shell*)Y->data;
891: if (X == Y) {
892: MatScale(Y,1.0 + a);
893: return(0);
894: }
895: if (!shell->axpy) {
896: MatConvertFrom_Shell(X,MATSHELL,MAT_INITIAL_MATRIX,&shell->axpy);
897: shell->axpy_vscale = a;
898: PetscObjectStateGet((PetscObject)X,&shell->axpy_state);
899: } else {
900: MatAXPY(shell->axpy,a/shell->axpy_vscale,X,str);
901: }
902: return(0);
903: }
905: static struct _MatOps MatOps_Values = {0,
906: 0,
907: 0,
908: 0,
909: /* 4*/ MatMultAdd_Shell,
910: 0,
911: MatMultTransposeAdd_Shell,
912: 0,
913: 0,
914: 0,
915: /*10*/ 0,
916: 0,
917: 0,
918: 0,
919: 0,
920: /*15*/ 0,
921: 0,
922: 0,
923: MatDiagonalScale_Shell,
924: 0,
925: /*20*/ 0,
926: MatAssemblyEnd_Shell,
927: 0,
928: 0,
929: /*24*/ MatZeroRows_Shell,
930: 0,
931: 0,
932: 0,
933: 0,
934: /*29*/ 0,
935: 0,
936: 0,
937: 0,
938: 0,
939: /*34*/ MatDuplicate_Shell,
940: 0,
941: 0,
942: 0,
943: 0,
944: /*39*/ MatAXPY_Shell,
945: 0,
946: 0,
947: 0,
948: MatCopy_Shell,
949: /*44*/ 0,
950: MatScale_Shell,
951: MatShift_Shell,
952: MatDiagonalSet_Shell,
953: MatZeroRowsColumns_Shell,
954: /*49*/ 0,
955: 0,
956: 0,
957: 0,
958: 0,
959: /*54*/ 0,
960: 0,
961: 0,
962: 0,
963: 0,
964: /*59*/ 0,
965: MatDestroy_Shell,
966: 0,
967: MatConvertFrom_Shell,
968: 0,
969: /*64*/ 0,
970: 0,
971: 0,
972: 0,
973: 0,
974: /*69*/ 0,
975: 0,
976: MatConvert_Shell,
977: 0,
978: 0,
979: /*74*/ 0,
980: 0,
981: 0,
982: 0,
983: 0,
984: /*79*/ 0,
985: 0,
986: 0,
987: 0,
988: 0,
989: /*84*/ 0,
990: 0,
991: 0,
992: 0,
993: 0,
994: /*89*/ 0,
995: 0,
996: 0,
997: 0,
998: 0,
999: /*94*/ 0,
1000: 0,
1001: 0,
1002: 0,
1003: 0,
1004: /*99*/ 0,
1005: 0,
1006: 0,
1007: 0,
1008: 0,
1009: /*104*/ 0,
1010: 0,
1011: 0,
1012: 0,
1013: 0,
1014: /*109*/ 0,
1015: 0,
1016: 0,
1017: 0,
1018: MatMissingDiagonal_Shell,
1019: /*114*/ 0,
1020: 0,
1021: 0,
1022: 0,
1023: 0,
1024: /*119*/ 0,
1025: 0,
1026: 0,
1027: 0,
1028: 0,
1029: /*124*/ 0,
1030: 0,
1031: 0,
1032: 0,
1033: 0,
1034: /*129*/ 0,
1035: 0,
1036: 0,
1037: 0,
1038: 0,
1039: /*134*/ 0,
1040: 0,
1041: 0,
1042: 0,
1043: 0,
1044: /*139*/ 0,
1045: 0,
1046: 0
1047: };
1049: PetscErrorCode MatShellSetContext_Shell(Mat mat,void *ctx)
1050: {
1051: Mat_Shell *shell = (Mat_Shell*)mat->data;
1054: shell->ctx = ctx;
1055: return(0);
1056: }
1058: static PetscErrorCode MatShellSetVecType_Shell(Mat mat,VecType vtype)
1059: {
1063: PetscFree(mat->defaultvectype);
1064: PetscStrallocpy(vtype,(char**)&mat->defaultvectype);
1065: return(0);
1066: }
1068: PetscErrorCode MatShellSetManageScalingShifts_Shell(Mat A)
1069: {
1070: Mat_Shell *shell = (Mat_Shell*)A->data;
1073: shell->managescalingshifts = PETSC_FALSE;
1074: A->ops->diagonalset = NULL;
1075: A->ops->diagonalscale = NULL;
1076: A->ops->scale = NULL;
1077: A->ops->shift = NULL;
1078: A->ops->axpy = NULL;
1079: return(0);
1080: }
1082: PetscErrorCode MatShellSetOperation_Shell(Mat mat,MatOperation op,void (*f)(void))
1083: {
1084: Mat_Shell *shell = (Mat_Shell*)mat->data;
1087: switch (op) {
1088: case MATOP_DESTROY:
1089: shell->ops->destroy = (PetscErrorCode (*)(Mat))f;
1090: break;
1091: case MATOP_VIEW:
1092: if (!mat->ops->viewnative) {
1093: mat->ops->viewnative = mat->ops->view;
1094: }
1095: mat->ops->view = (PetscErrorCode (*)(Mat,PetscViewer))f;
1096: break;
1097: case MATOP_COPY:
1098: shell->ops->copy = (PetscErrorCode (*)(Mat,Mat,MatStructure))f;
1099: break;
1100: case MATOP_DIAGONAL_SET:
1101: case MATOP_DIAGONAL_SCALE:
1102: case MATOP_SHIFT:
1103: case MATOP_SCALE:
1104: case MATOP_AXPY:
1105: case MATOP_ZERO_ROWS:
1106: case MATOP_ZERO_ROWS_COLUMNS:
1107: if (shell->managescalingshifts) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"MATSHELL is managing scalings and shifts, see MatShellSetManageScalingShifts()");
1108: (((void(**)(void))mat->ops)[op]) = f;
1109: break;
1110: case MATOP_GET_DIAGONAL:
1111: if (shell->managescalingshifts) {
1112: shell->ops->getdiagonal = (PetscErrorCode (*)(Mat,Vec))f;
1113: mat->ops->getdiagonal = MatGetDiagonal_Shell;
1114: } else {
1115: shell->ops->getdiagonal = NULL;
1116: mat->ops->getdiagonal = (PetscErrorCode (*)(Mat,Vec))f;
1117: }
1118: break;
1119: case MATOP_MULT:
1120: if (shell->managescalingshifts) {
1121: shell->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1122: mat->ops->mult = MatMult_Shell;
1123: } else {
1124: shell->ops->mult = NULL;
1125: mat->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1126: }
1127: break;
1128: case MATOP_MULT_TRANSPOSE:
1129: if (shell->managescalingshifts) {
1130: shell->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1131: mat->ops->multtranspose = MatMultTranspose_Shell;
1132: } else {
1133: shell->ops->multtranspose = NULL;
1134: mat->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1135: }
1136: break;
1137: default:
1138: (((void(**)(void))mat->ops)[op]) = f;
1139: break;
1140: }
1141: return(0);
1142: }
1144: PetscErrorCode MatShellGetOperation_Shell(Mat mat,MatOperation op,void(**f)(void))
1145: {
1146: Mat_Shell *shell = (Mat_Shell*)mat->data;
1149: switch (op) {
1150: case MATOP_DESTROY:
1151: *f = (void (*)(void))shell->ops->destroy;
1152: break;
1153: case MATOP_VIEW:
1154: *f = (void (*)(void))mat->ops->view;
1155: break;
1156: case MATOP_COPY:
1157: *f = (void (*)(void))shell->ops->copy;
1158: break;
1159: case MATOP_DIAGONAL_SET:
1160: case MATOP_DIAGONAL_SCALE:
1161: case MATOP_SHIFT:
1162: case MATOP_SCALE:
1163: case MATOP_AXPY:
1164: case MATOP_ZERO_ROWS:
1165: case MATOP_ZERO_ROWS_COLUMNS:
1166: *f = (((void (**)(void))mat->ops)[op]);
1167: break;
1168: case MATOP_GET_DIAGONAL:
1169: if (shell->ops->getdiagonal)
1170: *f = (void (*)(void))shell->ops->getdiagonal;
1171: else
1172: *f = (((void (**)(void))mat->ops)[op]);
1173: break;
1174: case MATOP_MULT:
1175: if (shell->ops->mult)
1176: *f = (void (*)(void))shell->ops->mult;
1177: else
1178: *f = (((void (**)(void))mat->ops)[op]);
1179: break;
1180: case MATOP_MULT_TRANSPOSE:
1181: if (shell->ops->multtranspose)
1182: *f = (void (*)(void))shell->ops->multtranspose;
1183: else
1184: *f = (((void (**)(void))mat->ops)[op]);
1185: break;
1186: default:
1187: *f = (((void (**)(void))mat->ops)[op]);
1188: }
1189: return(0);
1190: }
1192: /*MC
1193: MATSHELL - MATSHELL = "shell" - A matrix type to be used to define your own matrix type -- perhaps matrix free.
1195: Level: advanced
1197: .seealso: MatCreateShell()
1198: M*/
1200: PETSC_EXTERN PetscErrorCode MatCreate_Shell(Mat A)
1201: {
1202: Mat_Shell *b;
1206: PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));
1208: PetscNewLog(A,&b);
1209: A->data = (void*)b;
1211: PetscLayoutSetUp(A->rmap);
1212: PetscLayoutSetUp(A->cmap);
1214: b->ctx = 0;
1215: b->vshift = 0.0;
1216: b->vscale = 1.0;
1217: b->managescalingshifts = PETSC_TRUE;
1218: A->assembled = PETSC_TRUE;
1219: A->preallocated = PETSC_FALSE;
1221: PetscObjectComposeFunction((PetscObject)A,"MatShellGetContext_C",MatShellGetContext_Shell);
1222: PetscObjectComposeFunction((PetscObject)A,"MatShellSetContext_C",MatShellSetContext_Shell);
1223: PetscObjectComposeFunction((PetscObject)A,"MatShellSetVecType_C",MatShellSetVecType_Shell);
1224: PetscObjectComposeFunction((PetscObject)A,"MatShellSetManageScalingShifts_C",MatShellSetManageScalingShifts_Shell);
1225: PetscObjectComposeFunction((PetscObject)A,"MatShellSetOperation_C",MatShellSetOperation_Shell);
1226: PetscObjectComposeFunction((PetscObject)A,"MatShellGetOperation_C",MatShellGetOperation_Shell);
1227: PetscObjectChangeTypeName((PetscObject)A,MATSHELL);
1228: return(0);
1229: }
1231: /*@C
1232: MatCreateShell - Creates a new matrix class for use with a user-defined
1233: private data storage format.
1235: Collective
1237: Input Parameters:
1238: + comm - MPI communicator
1239: . m - number of local rows (must be given)
1240: . n - number of local columns (must be given)
1241: . M - number of global rows (may be PETSC_DETERMINE)
1242: . N - number of global columns (may be PETSC_DETERMINE)
1243: - ctx - pointer to data needed by the shell matrix routines
1245: Output Parameter:
1246: . A - the matrix
1248: Level: advanced
1250: Usage:
1251: $ extern PetscErrorCode mult(Mat,Vec,Vec);
1252: $ MatCreateShell(comm,m,n,M,N,ctx,&mat);
1253: $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
1254: $ [ Use matrix for operations that have been set ]
1255: $ MatDestroy(mat);
1257: Notes:
1258: The shell matrix type is intended to provide a simple class to use
1259: with KSP (such as, for use with matrix-free methods). You should not
1260: use the shell type if you plan to define a complete matrix class.
1262: Fortran Notes:
1263: To use this from Fortran with a ctx you must write an interface definition for this
1264: function and for MatShellGetContext() that tells Fortran the Fortran derived data type you are passing
1265: in as the ctx argument.
1267: PETSc requires that matrices and vectors being used for certain
1268: operations are partitioned accordingly. For example, when
1269: creating a shell matrix, A, that supports parallel matrix-vector
1270: products using MatMult(A,x,y) the user should set the number
1271: of local matrix rows to be the number of local elements of the
1272: corresponding result vector, y. Note that this is information is
1273: required for use of the matrix interface routines, even though
1274: the shell matrix may not actually be physically partitioned.
1275: For example,
1277: $
1278: $ Vec x, y
1279: $ extern PetscErrorCode mult(Mat,Vec,Vec);
1280: $ Mat A
1281: $
1282: $ VecCreateMPI(comm,PETSC_DECIDE,M,&y);
1283: $ VecCreateMPI(comm,PETSC_DECIDE,N,&x);
1284: $ VecGetLocalSize(y,&m);
1285: $ VecGetLocalSize(x,&n);
1286: $ MatCreateShell(comm,m,n,M,N,ctx,&A);
1287: $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
1288: $ MatMult(A,x,y);
1289: $ MatDestroy(&A);
1290: $ VecDestroy(&y);
1291: $ VecDestroy(&x);
1292: $
1295: MATSHELL handles MatShift(), MatDiagonalSet(), MatDiagonalScale(), MatAXPY(), MatScale(), MatZeroRows() and MatZeroRowsColumns() internally, so these
1296: operations cannot be overwritten unless MatShellSetManageScalingShifts() is called.
1299: For rectangular matrices do all the scalings and shifts make sense?
1301: Developers Notes:
1302: Regarding shifting and scaling. The general form is
1304: diag(left)(vscale*A + diag(dshift) + vshift I)diag(right)
1306: The order you apply the operations is important. For example if you have a dshift then
1307: apply a MatScale(s) you get s*vscale*A + s*diag(shift). But if you first scale and then shift
1308: you get s*vscale*A + diag(shift)
1310: A is the user provided function.
1312: KSP/PC uses changes in the Mat's "state" to decide if preconditioners need to be rebuilt: PCSetUp() only calls the setup() for
1313: for the PC implementation if the Mat state has increased from the previous call. Thus to get changes in a MATSHELL to trigger
1314: an update in the preconditioner you must call MatAssemblyBegin()/MatAssemblyEnd() or PetscObjectStateIncrease((PetscObject)mat);
1315: each time the MATSHELL matrix has changed.
1317: Calling MatAssemblyBegin()/MatAssemblyEnd() on a MATSHELL removes any previously supplied shift and scales that were provided
1318: with MatDiagonalSet(), MatShift(), MatScale(), or MatDiagonalScale().
1320: .seealso: MatShellSetOperation(), MatHasOperation(), MatShellGetContext(), MatShellSetContext(), MATSHELL, MatShellSetManageScalingShifts()
1321: @*/
1322: PetscErrorCode MatCreateShell(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,void *ctx,Mat *A)
1323: {
1327: MatCreate(comm,A);
1328: MatSetSizes(*A,m,n,M,N);
1329: MatSetType(*A,MATSHELL);
1330: MatShellSetContext(*A,ctx);
1331: MatSetUp(*A);
1332: return(0);
1333: }
1336: /*@
1337: MatShellSetContext - sets the context for a shell matrix
1339: Logically Collective on Mat
1341: Input Parameters:
1342: + mat - the shell matrix
1343: - ctx - the context
1345: Level: advanced
1347: Fortran Notes:
1348: To use this from Fortran you must write a Fortran interface definition for this
1349: function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.
1351: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation()
1352: @*/
1353: PetscErrorCode MatShellSetContext(Mat mat,void *ctx)
1354: {
1359: PetscTryMethod(mat,"MatShellSetContext_C",(Mat,void*),(mat,ctx));
1360: return(0);
1361: }
1363: /*@C
1364: MatShellSetVecType - Sets the type of Vec returned by MatCreateVecs()
1366: Logically collective
1368: Input Parameters:
1369: + mat - the shell matrix
1370: - vtype - type to use for creating vectors
1372: Notes:
1374: Level: advanced
1376: .seealso: MatCreateVecs()
1377: @*/
1378: PetscErrorCode MatShellSetVecType(Mat mat,VecType vtype)
1379: {
1383: PetscTryMethod(mat,"MatShellSetVecType_C",(Mat,VecType),(mat,vtype));
1384: return(0);
1385: }
1387: /*@
1388: MatShellSetManageScalingShifts - Allows the user to control the scaling and shift operations of the MATSHELL. Must be called immediately
1389: after MatCreateShell()
1391: Logically Collective on Mat
1393: Input Parameter:
1394: . mat - the shell matrix
1396: Level: advanced
1398: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatShellSetOperation()
1399: @*/
1400: PetscErrorCode MatShellSetManageScalingShifts(Mat A)
1401: {
1406: PetscTryMethod(A,"MatShellSetManageScalingShifts_C",(Mat),(A));
1407: return(0);
1408: }
1410: /*@C
1411: MatShellTestMult - Compares the multiply routine provided to the MATSHELL with differencing on a given function.
1413: Logically Collective on Mat
1415: Input Parameters:
1416: + mat - the shell matrix
1417: . f - the function
1418: . base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated
1419: - ctx - an optional context for the function
1421: Output Parameter:
1422: . flg - PETSC_TRUE if the multiply is likely correct
1424: Options Database:
1425: . -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference
1427: Level: advanced
1429: Fortran Notes:
1430: Not supported from Fortran
1432: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellTestMultTranspose()
1433: @*/
1434: PetscErrorCode MatShellTestMult(Mat mat,PetscErrorCode (*f)(void*,Vec,Vec),Vec base,void *ctx,PetscBool *flg)
1435: {
1437: PetscInt m,n;
1438: Mat mf,Dmf,Dmat,Ddiff;
1439: PetscReal Diffnorm,Dmfnorm;
1440: PetscBool v = PETSC_FALSE, flag = PETSC_TRUE;
1444: PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_view",&v);
1445: MatGetLocalSize(mat,&m,&n);
1446: MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,PETSC_DECIDE,PETSC_DECIDE,&mf);
1447: MatMFFDSetFunction(mf,f,ctx);
1448: MatMFFDSetBase(mf,base,NULL);
1450: MatComputeOperator(mf,MATAIJ,&Dmf);
1451: MatComputeOperator(mat,MATAIJ,&Dmat);
1453: MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff);
1454: MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN);
1455: MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm);
1456: MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm);
1457: if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) {
1458: flag = PETSC_FALSE;
1459: if (v) {
1460: 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));
1461: MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_view");
1462: MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_view");
1463: MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_view");
1464: }
1465: } else if (v) {
1466: PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL and matrix free multiple appear to produce the same results\n");
1467: }
1468: if (flg) *flg = flag;
1469: MatDestroy(&Ddiff);
1470: MatDestroy(&mf);
1471: MatDestroy(&Dmf);
1472: MatDestroy(&Dmat);
1473: return(0);
1474: }
1476: /*@C
1477: MatShellTestMultTranpose - Compares the multiply transpose routine provided to the MATSHELL with differencing on a given function.
1479: Logically Collective on Mat
1481: Input Parameters:
1482: + mat - the shell matrix
1483: . f - the function
1484: . base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated
1485: - ctx - an optional context for the function
1487: Output Parameter:
1488: . flg - PETSC_TRUE if the multiply is likely correct
1490: Options Database:
1491: . -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference
1493: Level: advanced
1495: Fortran Notes:
1496: Not supported from Fortran
1498: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellTestMult()
1499: @*/
1500: PetscErrorCode MatShellTestMultTranspose(Mat mat,PetscErrorCode (*f)(void*,Vec,Vec),Vec base,void *ctx,PetscBool *flg)
1501: {
1503: Vec x,y,z;
1504: PetscInt m,n,M,N;
1505: Mat mf,Dmf,Dmat,Ddiff;
1506: PetscReal Diffnorm,Dmfnorm;
1507: PetscBool v = PETSC_FALSE, flag = PETSC_TRUE;
1511: PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_transpose_view",&v);
1512: MatCreateVecs(mat,&x,&y);
1513: VecDuplicate(y,&z);
1514: MatGetLocalSize(mat,&m,&n);
1515: MatGetSize(mat,&M,&N);
1516: MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,M,N,&mf);
1517: MatMFFDSetFunction(mf,f,ctx);
1518: MatMFFDSetBase(mf,base,NULL);
1519: MatComputeOperator(mf,MATAIJ,&Dmf);
1520: MatTranspose(Dmf,MAT_INPLACE_MATRIX,&Dmf);
1521: MatComputeOperatorTranspose(mat,MATAIJ,&Dmat);
1523: MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff);
1524: MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN);
1525: MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm);
1526: MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm);
1527: if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) {
1528: flag = PETSC_FALSE;
1529: if (v) {
1530: 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));
1531: MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");
1532: MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");
1533: MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");
1534: }
1535: } else if (v) {
1536: PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL transpose and matrix free multiple appear to produce the same results\n");
1537: }
1538: if (flg) *flg = flag;
1539: MatDestroy(&mf);
1540: MatDestroy(&Dmat);
1541: MatDestroy(&Ddiff);
1542: MatDestroy(&Dmf);
1543: VecDestroy(&x);
1544: VecDestroy(&y);
1545: VecDestroy(&z);
1546: return(0);
1547: }
1549: /*@C
1550: MatShellSetOperation - Allows user to set a matrix operation for a shell matrix.
1552: Logically Collective on Mat
1554: Input Parameters:
1555: + mat - the shell matrix
1556: . op - the name of the operation
1557: - g - the function that provides the operation.
1559: Level: advanced
1561: Usage:
1562: $ extern PetscErrorCode usermult(Mat,Vec,Vec);
1563: $ MatCreateShell(comm,m,n,M,N,ctx,&A);
1564: $ MatShellSetOperation(A,MATOP_MULT,(void(*)(void))usermult);
1566: Notes:
1567: See the file include/petscmat.h for a complete list of matrix
1568: operations, which all have the form MATOP_<OPERATION>, where
1569: <OPERATION> is the name (in all capital letters) of the
1570: user interface routine (e.g., MatMult() -> MATOP_MULT).
1572: All user-provided functions (except for MATOP_DESTROY) should have the same calling
1573: sequence as the usual matrix interface routines, since they
1574: are intended to be accessed via the usual matrix interface
1575: routines, e.g.,
1576: $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec)
1578: In particular each function MUST return an error code of 0 on success and
1579: nonzero on failure.
1581: Within each user-defined routine, the user should call
1582: MatShellGetContext() to obtain the user-defined context that was
1583: set by MatCreateShell().
1585: Fortran Notes:
1586: For MatCreateVecs() the user code should check if the input left or right matrix is -1 and in that case not
1587: generate a matrix. See src/mat/tests/ex120f.F
1589: Use MatSetOperation() to set an operation for any matrix type
1591: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatSetOperation(), MatShellSetManageScalingShifts()
1592: @*/
1593: PetscErrorCode MatShellSetOperation(Mat mat,MatOperation op,void (*g)(void))
1594: {
1599: PetscTryMethod(mat,"MatShellSetOperation_C",(Mat,MatOperation,void (*)(void)),(mat,op,g));
1600: return(0);
1601: }
1603: /*@C
1604: MatShellGetOperation - Gets a matrix function for a shell matrix.
1606: Not Collective
1608: Input Parameters:
1609: + mat - the shell matrix
1610: - op - the name of the operation
1612: Output Parameter:
1613: . g - the function that provides the operation.
1615: Level: advanced
1617: Notes:
1618: See the file include/petscmat.h for a complete list of matrix
1619: operations, which all have the form MATOP_<OPERATION>, where
1620: <OPERATION> is the name (in all capital letters) of the
1621: user interface routine (e.g., MatMult() -> MATOP_MULT).
1623: All user-provided functions have the same calling
1624: sequence as the usual matrix interface routines, since they
1625: are intended to be accessed via the usual matrix interface
1626: routines, e.g.,
1627: $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec)
1629: Within each user-defined routine, the user should call
1630: MatShellGetContext() to obtain the user-defined context that was
1631: set by MatCreateShell().
1633: .seealso: MatCreateShell(), MatShellGetContext(), MatShellSetOperation(), MatShellSetContext()
1634: @*/
1635: PetscErrorCode MatShellGetOperation(Mat mat,MatOperation op,void(**g)(void))
1636: {
1641: PetscUseMethod(mat,"MatShellGetOperation_C",(Mat,MatOperation,void (**)(void)),(mat,op,g));
1642: return(0);
1643: }