Actual source code: shell.c
petsc-3.9.4 2018-09-11
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>
10: struct _MatShellOps {
11: /* 3 */ PetscErrorCode (*mult)(Mat,Vec,Vec);
12: /* 5 */ PetscErrorCode (*multtranspose)(Mat,Vec,Vec);
13: /* 17 */ PetscErrorCode (*getdiagonal)(Mat,Vec);
14: /* 43 */ PetscErrorCode (*copy)(Mat,Mat,MatStructure);
15: /* 60 */ PetscErrorCode (*destroy)(Mat);
16: };
18: typedef struct {
19: struct _MatShellOps ops[1];
21: PetscScalar vscale,vshift;
22: Vec dshift;
23: Vec left,right;
24: Vec left_work,right_work;
25: Vec left_add_work,right_add_work;
26: Mat axpy;
27: PetscScalar axpy_vscale;
28: PetscBool managescalingshifts; /* The user will manage the scaling and shifts for the MATSHELL, not the default */
29: void *ctx;
30: } Mat_Shell;
32: /*
33: xx = diag(left)*x
34: */
35: static PetscErrorCode MatShellPreScaleLeft(Mat A,Vec x,Vec *xx)
36: {
37: Mat_Shell *shell = (Mat_Shell*)A->data;
41: *xx = NULL;
42: if (!shell->left) {
43: *xx = x;
44: } else {
45: if (!shell->left_work) {VecDuplicate(shell->left,&shell->left_work);}
46: VecPointwiseMult(shell->left_work,x,shell->left);
47: *xx = shell->left_work;
48: }
49: return(0);
50: }
52: /*
53: xx = diag(right)*x
54: */
55: static PetscErrorCode MatShellPreScaleRight(Mat A,Vec x,Vec *xx)
56: {
57: Mat_Shell *shell = (Mat_Shell*)A->data;
61: *xx = NULL;
62: if (!shell->right) {
63: *xx = x;
64: } else {
65: if (!shell->right_work) {VecDuplicate(shell->right,&shell->right_work);}
66: VecPointwiseMult(shell->right_work,x,shell->right);
67: *xx = shell->right_work;
68: }
69: return(0);
70: }
72: /*
73: x = diag(left)*x
74: */
75: static PetscErrorCode MatShellPostScaleLeft(Mat A,Vec x)
76: {
77: Mat_Shell *shell = (Mat_Shell*)A->data;
81: if (shell->left) {VecPointwiseMult(x,x,shell->left);}
82: return(0);
83: }
85: /*
86: x = diag(right)*x
87: */
88: static PetscErrorCode MatShellPostScaleRight(Mat A,Vec x)
89: {
90: Mat_Shell *shell = (Mat_Shell*)A->data;
94: if (shell->right) {VecPointwiseMult(x,x,shell->right);}
95: return(0);
96: }
98: /*
99: Y = vscale*Y + diag(dshift)*X + vshift*X
101: On input Y already contains A*x
102: */
103: static PetscErrorCode MatShellShiftAndScale(Mat A,Vec X,Vec Y)
104: {
105: Mat_Shell *shell = (Mat_Shell*)A->data;
109: if (shell->dshift) { /* get arrays because there is no VecPointwiseMultAdd() */
110: PetscInt i,m;
111: const PetscScalar *x,*d;
112: PetscScalar *y;
113: VecGetLocalSize(X,&m);
114: VecGetArrayRead(shell->dshift,&d);
115: VecGetArrayRead(X,&x);
116: VecGetArray(Y,&y);
117: for (i=0; i<m; i++) y[i] = shell->vscale*y[i] + d[i]*x[i];
118: VecRestoreArrayRead(shell->dshift,&d);
119: VecRestoreArrayRead(X,&x);
120: VecRestoreArray(Y,&y);
121: } else {
122: VecScale(Y,shell->vscale);
123: }
124: if (shell->vshift != 0.0) {VecAXPY(Y,shell->vshift,X);} /* if test is for non-square matrices */
125: return(0);
126: }
128: /*@
129: MatShellGetContext - Returns the user-provided context associated with a shell matrix.
131: Not Collective
133: Input Parameter:
134: . mat - the matrix, should have been created with MatCreateShell()
136: Output Parameter:
137: . ctx - the user provided context
139: Level: advanced
141: Fortran Notes: To use this from Fortran you must write a Fortran interface definition for this
142: function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.
144: .keywords: matrix, shell, get, context
146: .seealso: MatCreateShell(), MatShellSetOperation(), MatShellSetContext()
147: @*/
148: PetscErrorCode MatShellGetContext(Mat mat,void *ctx)
149: {
151: PetscBool flg;
156: PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
157: if (flg) *(void**)ctx = ((Mat_Shell*)(mat->data))->ctx;
158: else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot get context from non-shell matrix");
159: return(0);
160: }
162: PetscErrorCode MatDestroy_Shell(Mat mat)
163: {
165: Mat_Shell *shell = (Mat_Shell*)mat->data;
168: if (shell->ops->destroy) {
169: (*shell->ops->destroy)(mat);
170: }
171: VecDestroy(&shell->left);
172: VecDestroy(&shell->right);
173: VecDestroy(&shell->dshift);
174: VecDestroy(&shell->left_work);
175: VecDestroy(&shell->right_work);
176: VecDestroy(&shell->left_add_work);
177: VecDestroy(&shell->right_add_work);
178: MatDestroy(&shell->axpy);
179: PetscFree(mat->data);
180: return(0);
181: }
183: PetscErrorCode MatCopy_Shell(Mat A,Mat B,MatStructure str)
184: {
185: Mat_Shell *shellA = (Mat_Shell*)A->data,*shellB = (Mat_Shell*)B->data;
186: PetscErrorCode ierr;
187: PetscBool matflg;
190: PetscObjectTypeCompare((PetscObject)B,MATSHELL,&matflg);
191: if (!matflg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_NOTSAMETYPE,"Output matrix must be a MATSHELL");
193: PetscMemcpy(B->ops,A->ops,sizeof(struct _MatOps));
194: PetscMemcpy(shellB->ops,shellA->ops,sizeof(struct _MatShellOps));
196: if (shellA->ops->copy) {
197: (*shellA->ops->copy)(A,B,str);
198: }
199: shellB->vscale = shellA->vscale;
200: shellB->vshift = shellA->vshift;
201: if (shellA->dshift) {
202: if (!shellB->dshift) {
203: VecDuplicate(shellA->dshift,&shellB->dshift);
204: }
205: VecCopy(shellA->dshift,shellB->dshift);
206: } else {
207: VecDestroy(&shellB->dshift);
208: }
209: if (shellA->left) {
210: if (!shellB->left) {
211: VecDuplicate(shellA->left,&shellB->left);
212: }
213: VecCopy(shellA->left,shellB->left);
214: } else {
215: VecDestroy(&shellB->left);
216: }
217: if (shellA->right) {
218: if (!shellB->right) {
219: VecDuplicate(shellA->right,&shellB->right);
220: }
221: VecCopy(shellA->right,shellB->right);
222: } else {
223: VecDestroy(&shellB->right);
224: }
225: MatDestroy(&shellB->axpy);
226: if (shellA->axpy) {
227: PetscObjectReference((PetscObject)shellA->axpy);
228: shellB->axpy = shellA->axpy;
229: shellB->axpy_vscale = shellA->axpy_vscale;
230: }
231: return(0);
232: }
234: PetscErrorCode MatDuplicate_Shell(Mat mat,MatDuplicateOption op,Mat *M)
235: {
237: void *ctx;
240: MatShellGetContext(mat,&ctx);
241: MatCreateShell(PetscObjectComm((PetscObject)mat),mat->rmap->n,mat->cmap->n,mat->rmap->N,mat->cmap->N,ctx,M);
242: MatCopy(mat,*M,SAME_NONZERO_PATTERN);
243: return(0);
244: }
246: PetscErrorCode MatMult_Shell(Mat A,Vec x,Vec y)
247: {
248: Mat_Shell *shell = (Mat_Shell*)A->data;
249: PetscErrorCode ierr;
250: Vec xx;
251: PetscObjectState instate,outstate;
254: MatShellPreScaleRight(A,x,&xx);
255: PetscObjectStateGet((PetscObject)y, &instate);
256: if (!shell->ops->mult) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Have not provided a MatMult() for this MATSHELL");
257: (*shell->ops->mult)(A,xx,y);
258: PetscObjectStateGet((PetscObject)y, &outstate);
259: if (instate == outstate) {
260: /* increase the state of the output vector since the user did not update its state themself as should have been done */
261: PetscObjectStateIncrease((PetscObject)y);
262: }
263: MatShellShiftAndScale(A,xx,y);
264: MatShellPostScaleLeft(A,y);
266: if (shell->axpy) {
267: if (!shell->left_work) {MatCreateVecs(A,&shell->left_work,NULL);}
268: MatMult(shell->axpy,x,shell->left_work);
269: VecAXPY(y,shell->axpy_vscale,shell->left_work);
270: }
271: return(0);
272: }
274: PetscErrorCode MatMultAdd_Shell(Mat A,Vec x,Vec y,Vec z)
275: {
276: Mat_Shell *shell = (Mat_Shell*)A->data;
280: if (y == z) {
281: if (!shell->right_add_work) {VecDuplicate(z,&shell->right_add_work);}
282: MatMult(A,x,shell->right_add_work);
283: VecAXPY(z,1.0,shell->right_add_work);
284: } else {
285: MatMult(A,x,z);
286: VecAXPY(z,1.0,y);
287: }
288: return(0);
289: }
291: PetscErrorCode MatMultTranspose_Shell(Mat A,Vec x,Vec y)
292: {
293: Mat_Shell *shell = (Mat_Shell*)A->data;
294: PetscErrorCode ierr;
295: Vec xx;
296: PetscObjectState instate,outstate;
299: MatShellPreScaleLeft(A,x,&xx);
300: PetscObjectStateGet((PetscObject)y, &instate);
301: if (!shell->ops->multtranspose) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Have not provided a MatMultTranspose() for this MATSHELL");
302: (*shell->ops->multtranspose)(A,xx,y);
303: PetscObjectStateGet((PetscObject)y, &outstate);
304: if (instate == outstate) {
305: /* increase the state of the output vector since the user did not update its state themself as should have been done */
306: PetscObjectStateIncrease((PetscObject)y);
307: }
308: MatShellShiftAndScale(A,xx,y);
309: MatShellPostScaleRight(A,y);
310: return(0);
311: }
313: PetscErrorCode MatMultTransposeAdd_Shell(Mat A,Vec x,Vec y,Vec z)
314: {
315: Mat_Shell *shell = (Mat_Shell*)A->data;
319: if (y == z) {
320: if (!shell->left_add_work) {VecDuplicate(z,&shell->left_add_work);}
321: MatMultTranspose(A,x,shell->left_add_work);
322: VecWAXPY(z,1.0,shell->left_add_work,y);
323: } else {
324: MatMultTranspose(A,x,z);
325: VecAXPY(z,1.0,y);
326: }
327: return(0);
328: }
330: /*
331: diag(left)(vscale*A + diag(dshift) + vshift I)diag(right)
332: */
333: PetscErrorCode MatGetDiagonal_Shell(Mat A,Vec v)
334: {
335: Mat_Shell *shell = (Mat_Shell*)A->data;
339: if (shell->ops->getdiagonal) {
340: (*shell->ops->getdiagonal)(A,v);
341: } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Must provide shell matrix with routine to return diagonal using\nMatShellSetOperation(S,MATOP_GET_DIAGONAL,...)");
342: VecScale(v,shell->vscale);
343: if (shell->dshift) {
344: VecAXPY(v,1.0,shell->dshift);
345: }
346: VecShift(v,shell->vshift);
347: if (shell->left) {VecPointwiseMult(v,v,shell->left);}
348: if (shell->right) {VecPointwiseMult(v,v,shell->right);}
349: if (shell->axpy) {
350: if (!shell->left_work) {VecDuplicate(v,&shell->left_work);}
351: MatGetDiagonal(shell->axpy,shell->left_work);
352: VecAXPY(v,shell->axpy_vscale,shell->left_work);
353: }
354: return(0);
355: }
357: PetscErrorCode MatShift_Shell(Mat Y,PetscScalar a)
358: {
359: Mat_Shell *shell = (Mat_Shell*)Y->data;
363: if (shell->left || shell->right) {
364: if (!shell->dshift) {
365: VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift);
366: VecSet(shell->dshift,a);
367: } else {
368: if (shell->left) {VecPointwiseMult(shell->dshift,shell->dshift,shell->left);}
369: if (shell->right) {VecPointwiseMult(shell->dshift,shell->dshift,shell->right);}
370: VecShift(shell->dshift,a);
371: }
372: if (shell->left) {VecPointwiseDivide(shell->dshift,shell->dshift,shell->left);}
373: if (shell->right) {VecPointwiseDivide(shell->dshift,shell->dshift,shell->right);}
374: } else shell->vshift += a;
375: return(0);
376: }
378: PetscErrorCode MatDiagonalSet_Shell(Mat A,Vec D,InsertMode ins)
379: {
380: Mat_Shell *shell = (Mat_Shell*)A->data;
384: if (ins == INSERT_VALUES) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE, "Operation not supported with INSERT_VALUES");
385: if (!shell->dshift) {VecDuplicate(D,&shell->dshift);}
386: if (shell->left || shell->right) {
387: if (!shell->right_work) {VecDuplicate(shell->left ? shell->left : shell->right, &shell->right_work);}
388: if (shell->left && shell->right) {
389: VecPointwiseDivide(shell->right_work,D,shell->left);
390: VecPointwiseDivide(shell->right_work,shell->right_work,shell->right);
391: } else if (shell->left) {
392: VecPointwiseDivide(shell->right_work,D,shell->left);
393: } else {
394: VecPointwiseDivide(shell->right_work,D,shell->right);
395: }
396: if (!shell->dshift) {
397: VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift);
398: VecCopy(shell->dshift,shell->right_work);
399: } else {
400: VecAXPY(shell->dshift,1.0,shell->right_work);
401: }
402: } else {
403: VecAXPY(shell->dshift,1.0,D);
404: }
405: return(0);
406: }
408: PetscErrorCode MatScale_Shell(Mat Y,PetscScalar a)
409: {
410: Mat_Shell *shell = (Mat_Shell*)Y->data;
414: shell->vscale *= a;
415: shell->vshift *= a;
416: if (shell->dshift) {
417: VecScale(shell->dshift,a);
418: }
419: shell->axpy_vscale *= a;
420: return(0);
421: }
423: static PetscErrorCode MatDiagonalScale_Shell(Mat Y,Vec left,Vec right)
424: {
425: Mat_Shell *shell = (Mat_Shell*)Y->data;
429: if (shell->axpy) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_SUP,"Cannot diagonal scale MATSHELL after MatAXPY operation");
430: if (left) {
431: if (!shell->left) {
432: VecDuplicate(left,&shell->left);
433: VecCopy(left,shell->left);
434: } else {
435: VecPointwiseMult(shell->left,shell->left,left);
436: }
437: }
438: if (right) {
439: if (!shell->right) {
440: VecDuplicate(right,&shell->right);
441: VecCopy(right,shell->right);
442: } else {
443: VecPointwiseMult(shell->right,shell->right,right);
444: }
445: }
446: return(0);
447: }
449: PetscErrorCode MatAssemblyEnd_Shell(Mat Y,MatAssemblyType t)
450: {
451: Mat_Shell *shell = (Mat_Shell*)Y->data;
455: if (t == MAT_FINAL_ASSEMBLY) {
456: shell->vshift = 0.0;
457: shell->vscale = 1.0;
458: VecDestroy(&shell->dshift);
459: VecDestroy(&shell->left);
460: VecDestroy(&shell->right);
461: MatDestroy(&shell->axpy);
462: }
463: return(0);
464: }
466: PETSC_INTERN PetscErrorCode MatConvert_Shell(Mat, MatType,MatReuse,Mat*);
468: static PetscErrorCode MatMissingDiagonal_Shell(Mat A,PetscBool *missing,PetscInt *d)
469: {
471: *missing = PETSC_FALSE;
472: return(0);
473: }
475: PetscErrorCode MatAXPY_Shell(Mat Y,PetscScalar a,Mat X,MatStructure str)
476: {
477: Mat_Shell *shell = (Mat_Shell*)Y->data;
481: PetscObjectReference((PetscObject)X);
482: MatDestroy(&shell->axpy);
483: shell->axpy = X;
484: shell->axpy_vscale = a;
485: return(0);
486: }
488: static struct _MatOps MatOps_Values = {0,
489: 0,
490: 0,
491: 0,
492: /* 4*/ MatMultAdd_Shell,
493: 0,
494: MatMultTransposeAdd_Shell,
495: 0,
496: 0,
497: 0,
498: /*10*/ 0,
499: 0,
500: 0,
501: 0,
502: 0,
503: /*15*/ 0,
504: 0,
505: MatGetDiagonal_Shell,
506: MatDiagonalScale_Shell,
507: 0,
508: /*20*/ 0,
509: MatAssemblyEnd_Shell,
510: 0,
511: 0,
512: /*24*/ 0,
513: 0,
514: 0,
515: 0,
516: 0,
517: /*29*/ 0,
518: 0,
519: 0,
520: 0,
521: 0,
522: /*34*/ MatDuplicate_Shell,
523: 0,
524: 0,
525: 0,
526: 0,
527: /*39*/ MatAXPY_Shell,
528: 0,
529: 0,
530: 0,
531: MatCopy_Shell,
532: /*44*/ 0,
533: MatScale_Shell,
534: MatShift_Shell,
535: MatDiagonalSet_Shell,
536: 0,
537: /*49*/ 0,
538: 0,
539: 0,
540: 0,
541: 0,
542: /*54*/ 0,
543: 0,
544: 0,
545: 0,
546: 0,
547: /*59*/ 0,
548: MatDestroy_Shell,
549: 0,
550: 0,
551: 0,
552: /*64*/ 0,
553: 0,
554: 0,
555: 0,
556: 0,
557: /*69*/ 0,
558: 0,
559: MatConvert_Shell,
560: 0,
561: 0,
562: /*74*/ 0,
563: 0,
564: 0,
565: 0,
566: 0,
567: /*79*/ 0,
568: 0,
569: 0,
570: 0,
571: 0,
572: /*84*/ 0,
573: 0,
574: 0,
575: 0,
576: 0,
577: /*89*/ 0,
578: 0,
579: 0,
580: 0,
581: 0,
582: /*94*/ 0,
583: 0,
584: 0,
585: 0,
586: 0,
587: /*99*/ 0,
588: 0,
589: 0,
590: 0,
591: 0,
592: /*104*/ 0,
593: 0,
594: 0,
595: 0,
596: 0,
597: /*109*/ 0,
598: 0,
599: 0,
600: 0,
601: MatMissingDiagonal_Shell,
602: /*114*/ 0,
603: 0,
604: 0,
605: 0,
606: 0,
607: /*119*/ 0,
608: 0,
609: 0,
610: 0,
611: 0,
612: /*124*/ 0,
613: 0,
614: 0,
615: 0,
616: 0,
617: /*129*/ 0,
618: 0,
619: 0,
620: 0,
621: 0,
622: /*134*/ 0,
623: 0,
624: 0,
625: 0,
626: 0,
627: /*139*/ 0,
628: 0,
629: 0
630: };
632: /*MC
633: MATSHELL - MATSHELL = "shell" - A matrix type to be used to define your own matrix type -- perhaps matrix free.
635: Level: advanced
637: .seealso: MatCreateShell()
638: M*/
640: PETSC_EXTERN PetscErrorCode MatCreate_Shell(Mat A)
641: {
642: Mat_Shell *b;
646: PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));
648: PetscNewLog(A,&b);
649: A->data = (void*)b;
651: PetscLayoutSetUp(A->rmap);
652: PetscLayoutSetUp(A->cmap);
654: b->ctx = 0;
655: b->vshift = 0.0;
656: b->vscale = 1.0;
657: b->managescalingshifts = PETSC_TRUE;
658: A->assembled = PETSC_TRUE;
659: A->preallocated = PETSC_FALSE;
661: PetscObjectChangeTypeName((PetscObject)A,MATSHELL);
662: return(0);
663: }
665: /*@C
666: MatCreateShell - Creates a new matrix class for use with a user-defined
667: private data storage format.
669: Collective on MPI_Comm
671: Input Parameters:
672: + comm - MPI communicator
673: . m - number of local rows (must be given)
674: . n - number of local columns (must be given)
675: . M - number of global rows (may be PETSC_DETERMINE)
676: . N - number of global columns (may be PETSC_DETERMINE)
677: - ctx - pointer to data needed by the shell matrix routines
679: Output Parameter:
680: . A - the matrix
682: Level: advanced
684: Usage:
685: $ extern int mult(Mat,Vec,Vec);
686: $ MatCreateShell(comm,m,n,M,N,ctx,&mat);
687: $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
688: $ [ Use matrix for operations that have been set ]
689: $ MatDestroy(mat);
691: Notes:
692: The shell matrix type is intended to provide a simple class to use
693: with KSP (such as, for use with matrix-free methods). You should not
694: use the shell type if you plan to define a complete matrix class.
696: Fortran Notes: To use this from Fortran with a ctx you must write an interface definition for this
697: function and for MatShellGetContext() that tells Fortran the Fortran derived data type you are passing
698: in as the ctx argument.
700: PETSc requires that matrices and vectors being used for certain
701: operations are partitioned accordingly. For example, when
702: creating a shell matrix, A, that supports parallel matrix-vector
703: products using MatMult(A,x,y) the user should set the number
704: of local matrix rows to be the number of local elements of the
705: corresponding result vector, y. Note that this is information is
706: required for use of the matrix interface routines, even though
707: the shell matrix may not actually be physically partitioned.
708: For example,
710: $
711: $ Vec x, y
712: $ extern int mult(Mat,Vec,Vec);
713: $ Mat A
714: $
715: $ VecCreateMPI(comm,PETSC_DECIDE,M,&y);
716: $ VecCreateMPI(comm,PETSC_DECIDE,N,&x);
717: $ VecGetLocalSize(y,&m);
718: $ VecGetLocalSize(x,&n);
719: $ MatCreateShell(comm,m,n,M,N,ctx,&A);
720: $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
721: $ MatMult(A,x,y);
722: $ MatDestroy(A);
723: $ VecDestroy(y); VecDestroy(x);
724: $
727: MATSHELL handles MatShift(), MatDiagonalSet(), MatDiagonalScale(), MatAXPY(), and MatScale() internally so these
728: operations cannot be overwritten unless MatShellSetManageScalingShifts() is called.
731: For rectangular matrices do all the scalings and shifts make sense?
733: Developers Notes: Regarding shifting and scaling. The general form is
735: diag(left)(vscale*A + diag(dshift) + vshift I)diag(right)
737: The order you apply the operations is important. For example if you have a dshift then
738: apply a MatScale(s) you get s*vscale*A + s*diag(shift). But if you first scale and then shift
739: you get s*vscale*A + diag(shift)
741: A is the user provided function.
743: .keywords: matrix, shell, create
745: .seealso: MatShellSetOperation(), MatHasOperation(), MatShellGetContext(), MatShellSetContext(), MATSHELL, MatShellSetManageScalingShifts()
746: @*/
747: PetscErrorCode MatCreateShell(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,void *ctx,Mat *A)
748: {
752: MatCreate(comm,A);
753: MatSetSizes(*A,m,n,M,N);
754: MatSetType(*A,MATSHELL);
755: MatShellSetContext(*A,ctx);
756: MatSetUp(*A);
757: return(0);
758: }
760: /*@
761: MatShellSetContext - sets the context for a shell matrix
763: Logically Collective on Mat
765: Input Parameters:
766: + mat - the shell matrix
767: - ctx - the context
769: Level: advanced
771: Fortran Notes: To use this from Fortran you must write a Fortran interface definition for this
772: function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.
774: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation()
775: @*/
776: PetscErrorCode MatShellSetContext(Mat mat,void *ctx)
777: {
778: Mat_Shell *shell = (Mat_Shell*)mat->data;
780: PetscBool flg;
784: PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
785: if (flg) {
786: shell->ctx = ctx;
787: } else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot attach context to non-shell matrix");
788: return(0);
789: }
791: /*@
792: MatShellSetManageScalingShifts - Allows the user to control the scaling and shift operations of the MATSHELL. Must be called immediately
793: after MatCreateShell()
795: Logically Collective on Mat
797: Input Parameter:
798: . mat - the shell matrix
800: Level: advanced
802: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatShellSetOperation()
803: @*/
804: PetscErrorCode MatShellSetManageScalingShifts(Mat A)
805: {
807: Mat_Shell *shell;
808: PetscBool flg;
812: PetscObjectTypeCompare((PetscObject)A,MATSHELL,&flg);
813: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Can only use with MATSHELL matrices");
814: shell = (Mat_Shell*)A->data;
815: shell->managescalingshifts = PETSC_FALSE;
816: A->ops->diagonalset = NULL;
817: A->ops->diagonalscale = NULL;
818: A->ops->scale = NULL;
819: A->ops->shift = NULL;
820: A->ops->axpy = NULL;
821: return(0);
822: }
824: /*@C
825: MatShellTestMult - Compares the multiply routine provided to the MATSHELL with differencing on a given function.
827: Logically Collective on Mat
829: Input Parameters:
830: + mat - the shell matrix
831: . f - the function
832: . base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated
833: - ctx - an optional context for the function
835: Output Parameter:
836: . flg - PETSC_TRUE if the multiply is likely correct
838: Options Database:
839: . -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference
841: Level: advanced
843: Fortran Notes: Not supported from Fortran
845: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellTestMultTranspose()
846: @*/
847: PetscErrorCode MatShellTestMult(Mat mat,PetscErrorCode (*f)(void*,Vec,Vec),Vec base,void *ctx,PetscBool *flg)
848: {
850: PetscInt m,n;
851: Mat mf,Dmf,Dmat,Ddiff;
852: PetscReal Diffnorm,Dmfnorm;
853: PetscBool v = PETSC_FALSE, flag = PETSC_TRUE;
857: PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_view",&v);
858: MatGetLocalSize(mat,&m,&n);
859: MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,PETSC_DECIDE,PETSC_DECIDE,&mf);
860: MatMFFDSetFunction(mf,f,ctx);
861: MatMFFDSetBase(mf,base,NULL);
863: MatComputeExplicitOperator(mf,&Dmf);
864: MatComputeExplicitOperator(mat,&Dmat);
866: MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff);
867: MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN);
868: MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm);
869: MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm);
870: if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) {
871: flag = PETSC_FALSE;
872: if (v) {
873: 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));
874: MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_view");
875: MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_view");
876: MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_view");
877: }
878: } else if (v) {
879: PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL and matrix free multiple appear to produce the same results\n");
880: }
881: if (flg) *flg = flag;
882: MatDestroy(&Ddiff);
883: MatDestroy(&mf);
884: MatDestroy(&Dmf);
885: MatDestroy(&Dmat);
886: return(0);
887: }
889: /*@C
890: MatShellTestMultTranpose - Compares the multiply transpose routine provided to the MATSHELL with differencing on a given function.
892: Logically Collective on Mat
894: Input Parameters:
895: + mat - the shell matrix
896: . f - the function
897: . base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated
898: - ctx - an optional context for the function
900: Output Parameter:
901: . flg - PETSC_TRUE if the multiply is likely correct
903: Options Database:
904: . -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference
906: Level: advanced
908: Fortran Notes: Not supported from Fortran
910: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellTestMult()
911: @*/
912: PetscErrorCode MatShellTestMultTranspose(Mat mat,PetscErrorCode (*f)(void*,Vec,Vec),Vec base,void *ctx,PetscBool *flg)
913: {
915: Vec x,y,z;
916: PetscInt m,n,M,N;
917: Mat mf,Dmf,Dmat,Ddiff;
918: PetscReal Diffnorm,Dmfnorm;
919: PetscBool v = PETSC_FALSE, flag = PETSC_TRUE;
923: PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_transpose_view",&v);
924: MatCreateVecs(mat,&x,&y);
925: VecDuplicate(y,&z);
926: MatGetLocalSize(mat,&m,&n);
927: MatGetSize(mat,&M,&N);
928: MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,M,N,&mf);
929: MatMFFDSetFunction(mf,f,ctx);
930: MatMFFDSetBase(mf,base,NULL);
931: MatComputeExplicitOperator(mf,&Dmf);
932: MatTranspose(Dmf,MAT_INPLACE_MATRIX,&Dmf);
933: MatComputeExplicitOperatorTranspose(mat,&Dmat);
935: MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff);
936: MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN);
937: MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm);
938: MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm);
939: if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) {
940: flag = PETSC_FALSE;
941: if (v) {
942: 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));
943: MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");
944: MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");
945: MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");
946: }
947: } else if (v) {
948: PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL transpose and matrix free multiple appear to produce the same results\n");
949: }
950: if (flg) *flg = flag;
951: MatDestroy(&mf);
952: MatDestroy(&Dmat);
953: MatDestroy(&Ddiff);
954: MatDestroy(&Dmf);
955: VecDestroy(&x);
956: VecDestroy(&y);
957: VecDestroy(&z);
958: return(0);
959: }
961: /*@C
962: MatShellSetOperation - Allows user to set a matrix operation for a shell matrix.
964: Logically Collective on Mat
966: Input Parameters:
967: + mat - the shell matrix
968: . op - the name of the operation
969: - f - the function that provides the operation.
971: Level: advanced
973: Usage:
974: $ extern PetscErrorCode usermult(Mat,Vec,Vec);
975: $ MatCreateShell(comm,m,n,M,N,ctx,&A);
976: $ MatShellSetOperation(A,MATOP_MULT,(void(*)(void))usermult);
978: Notes:
979: See the file include/petscmat.h for a complete list of matrix
980: operations, which all have the form MATOP_<OPERATION>, where
981: <OPERATION> is the name (in all capital letters) of the
982: user interface routine (e.g., MatMult() -> MATOP_MULT).
984: All user-provided functions (except for MATOP_DESTROY) should have the same calling
985: sequence as the usual matrix interface routines, since they
986: are intended to be accessed via the usual matrix interface
987: routines, e.g.,
988: $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec)
990: In particular each function MUST return an error code of 0 on success and
991: nonzero on failure.
993: Within each user-defined routine, the user should call
994: MatShellGetContext() to obtain the user-defined context that was
995: set by MatCreateShell().
997: Fortran Notes: For MatCreateVecs() the user code should check if the input left or right matrix is -1 and in that case not
998: generate a matrix. See src/mat/examples/tests/ex120f.F
1000: Use MatSetOperation() to set an operation for any matrix type
1002: .keywords: matrix, shell, set, operation
1004: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatSetOperation(), MatShellSetManageScalingShifts()
1005: @*/
1006: PetscErrorCode MatShellSetOperation(Mat mat,MatOperation op,void (*f)(void))
1007: {
1008: PetscBool flg;
1009: Mat_Shell *shell;
1014: PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
1015: if (!flg) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Can only use with MATSHELL matrices");
1016: shell = (Mat_Shell*)mat->data;
1018: switch (op) {
1019: case MATOP_DESTROY:
1020: shell->ops->destroy = (PetscErrorCode (*)(Mat))f;
1021: break;
1022: case MATOP_VIEW:
1023: if (!mat->ops->viewnative) {
1024: mat->ops->viewnative = mat->ops->view;
1025: }
1026: mat->ops->view = (PetscErrorCode (*)(Mat,PetscViewer))f;
1027: break;
1028: case MATOP_COPY:
1029: shell->ops->copy = (PetscErrorCode (*)(Mat,Mat,MatStructure))f;
1030: break;
1031: case MATOP_DIAGONAL_SET:
1032: case MATOP_DIAGONAL_SCALE:
1033: case MATOP_SHIFT:
1034: case MATOP_SCALE:
1035: case MATOP_AXPY:
1036: if (shell->managescalingshifts) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"MATSHELL is managing scalings and shifts, see MatShellSetManageScalingShifts()");
1037: (((void(**)(void))mat->ops)[op]) = f;
1038: break;
1039: case MATOP_GET_DIAGONAL:
1040: if (shell->managescalingshifts) {
1041: shell->ops->getdiagonal = (PetscErrorCode (*)(Mat,Vec))f;
1042: mat->ops->getdiagonal = MatGetDiagonal_Shell;
1043: } else {
1044: shell->ops->getdiagonal = NULL;
1045: mat->ops->getdiagonal = (PetscErrorCode (*)(Mat,Vec))f;
1046: }
1047: break;
1048: case MATOP_MULT:
1049: if (shell->managescalingshifts) {
1050: shell->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1051: mat->ops->mult = MatMult_Shell;
1052: } else {
1053: shell->ops->mult = NULL;
1054: mat->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1055: }
1056: break;
1057: case MATOP_MULT_TRANSPOSE:
1058: if (shell->managescalingshifts) {
1059: shell->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1060: mat->ops->multtranspose = MatMultTranspose_Shell;
1061: } else {
1062: shell->ops->multtranspose = NULL;
1063: mat->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1064: }
1065: break;
1066: default:
1067: (((void(**)(void))mat->ops)[op]) = f;
1068: break;
1069: }
1070: return(0);
1071: }
1073: /*@C
1074: MatShellGetOperation - Gets a matrix function for a shell matrix.
1076: Not Collective
1078: Input Parameters:
1079: + mat - the shell matrix
1080: - op - the name of the operation
1082: Output Parameter:
1083: . f - the function that provides the operation.
1085: Level: advanced
1087: Notes:
1088: See the file include/petscmat.h for a complete list of matrix
1089: operations, which all have the form MATOP_<OPERATION>, where
1090: <OPERATION> is the name (in all capital letters) of the
1091: user interface routine (e.g., MatMult() -> MATOP_MULT).
1093: All user-provided functions have the same calling
1094: sequence as the usual matrix interface routines, since they
1095: are intended to be accessed via the usual matrix interface
1096: routines, e.g.,
1097: $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec)
1099: Within each user-defined routine, the user should call
1100: MatShellGetContext() to obtain the user-defined context that was
1101: set by MatCreateShell().
1103: .keywords: matrix, shell, set, operation
1105: .seealso: MatCreateShell(), MatShellGetContext(), MatShellSetOperation(), MatShellSetContext()
1106: @*/
1107: PetscErrorCode MatShellGetOperation(Mat mat,MatOperation op,void(**f)(void))
1108: {
1109: PetscBool flg;
1110: Mat_Shell *shell;
1115: PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
1116: if (!flg) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Can only use with MATSHELL matrices");
1117: shell = (Mat_Shell*)mat->data;
1119: switch (op) {
1120: case MATOP_DESTROY:
1121: *f = (void (*)(void))shell->ops->destroy;
1122: break;
1123: case MATOP_VIEW:
1124: *f = (void (*)(void))mat->ops->view;
1125: break;
1126: case MATOP_COPY:
1127: *f = (void (*)(void))shell->ops->copy;
1128: break;
1129: case MATOP_DIAGONAL_SET:
1130: case MATOP_DIAGONAL_SCALE:
1131: case MATOP_SHIFT:
1132: case MATOP_SCALE:
1133: case MATOP_AXPY:
1134: *f = (((void (**)(void))mat->ops)[op]);
1135: break;
1136: case MATOP_GET_DIAGONAL:
1137: if (shell->ops->getdiagonal)
1138: *f = (void (*)(void))shell->ops->getdiagonal;
1139: else
1140: *f = (((void (**)(void))mat->ops)[op]);
1141: break;
1142: case MATOP_MULT:
1143: if (shell->ops->mult)
1144: *f = (void (*)(void))shell->ops->mult;
1145: else
1146: *f = (((void (**)(void))mat->ops)[op]);
1147: break;
1148: case MATOP_MULT_TRANSPOSE:
1149: if (shell->ops->multtranspose)
1150: *f = (void (*)(void))shell->ops->multtranspose;
1151: else
1152: *f = (((void (**)(void))mat->ops)[op]);
1153: break;
1154: default:
1155: *f = (((void (**)(void))mat->ops)[op]);
1156: }
1157: return(0);
1158: }