Actual source code: shell.c
petsc-3.7.5 2017-01-01
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> /*I "petscmat.h" I*/
9: #include <petsc/private/vecimpl.h>
11: typedef struct {
12: PetscErrorCode (*destroy)(Mat);
13: PetscErrorCode (*mult)(Mat,Vec,Vec);
14: PetscErrorCode (*multtranspose)(Mat,Vec,Vec);
15: PetscErrorCode (*getdiagonal)(Mat,Vec);
17: PetscScalar vscale,vshift;
18: Vec dshift;
19: Vec left,right;
20: Vec dshift_owned,left_owned,right_owned;
21: Vec left_work,right_work;
22: Vec left_add_work,right_add_work;
23: PetscBool usingscaled;
24: void *ctx;
25: } Mat_Shell;
26: /*
27: The most general expression for the matrix is
29: S = L (a A + B) R
31: where
32: A is the matrix defined by the user's function
33: a is a scalar multiple
34: L is left scaling
35: R is right scaling
36: B is a diagonal shift defined by
37: diag(dshift) if the vector dshift is non-NULL
38: vscale*identity otherwise
40: The following identities apply:
42: Scale by c:
43: c [L (a A + B) R] = L [(a c) A + c B] R
45: Shift by c:
46: [L (a A + B) R] + c = L [a A + (B + c Linv Rinv)] R
48: Diagonal scaling is achieved by simply multiplying with existing L and R vectors
50: In the data structure:
52: vscale=1.0 means no special scaling will be applied
53: dshift=NULL means a constant diagonal shift (fall back to vshift)
54: vshift=0.0 means no constant diagonal shift, note that vshift is only used if dshift is NULL
55: */
57: static PetscErrorCode MatMult_Shell(Mat,Vec,Vec);
58: static PetscErrorCode MatMultTranspose_Shell(Mat,Vec,Vec);
59: static PetscErrorCode MatGetDiagonal_Shell(Mat,Vec);
63: static PetscErrorCode MatShellUseScaledMethods(Mat Y)
64: {
65: Mat_Shell *shell = (Mat_Shell*)Y->data;
68: if (shell->usingscaled) return(0);
69: shell->mult = Y->ops->mult;
70: Y->ops->mult = MatMult_Shell;
71: if (Y->ops->multtranspose) {
72: shell->multtranspose = Y->ops->multtranspose;
73: Y->ops->multtranspose = MatMultTranspose_Shell;
74: }
75: if (Y->ops->getdiagonal) {
76: shell->getdiagonal = Y->ops->getdiagonal;
77: Y->ops->getdiagonal = MatGetDiagonal_Shell;
78: }
79: shell->usingscaled = PETSC_TRUE;
80: return(0);
81: }
85: static PetscErrorCode MatShellPreScaleLeft(Mat A,Vec x,Vec *xx)
86: {
87: Mat_Shell *shell = (Mat_Shell*)A->data;
91: *xx = NULL;
92: if (!shell->left) {
93: *xx = x;
94: } else {
95: if (!shell->left_work) {VecDuplicate(shell->left,&shell->left_work);}
96: VecPointwiseMult(shell->left_work,x,shell->left);
97: *xx = shell->left_work;
98: }
99: return(0);
100: }
104: static PetscErrorCode MatShellPreScaleRight(Mat A,Vec x,Vec *xx)
105: {
106: Mat_Shell *shell = (Mat_Shell*)A->data;
110: *xx = NULL;
111: if (!shell->right) {
112: *xx = x;
113: } else {
114: if (!shell->right_work) {VecDuplicate(shell->right,&shell->right_work);}
115: VecPointwiseMult(shell->right_work,x,shell->right);
116: *xx = shell->right_work;
117: }
118: return(0);
119: }
123: static PetscErrorCode MatShellPostScaleLeft(Mat A,Vec x)
124: {
125: Mat_Shell *shell = (Mat_Shell*)A->data;
129: if (shell->left) {VecPointwiseMult(x,x,shell->left);}
130: return(0);
131: }
135: static PetscErrorCode MatShellPostScaleRight(Mat A,Vec x)
136: {
137: Mat_Shell *shell = (Mat_Shell*)A->data;
141: if (shell->right) {VecPointwiseMult(x,x,shell->right);}
142: return(0);
143: }
147: static PetscErrorCode MatShellShiftAndScale(Mat A,Vec X,Vec Y)
148: {
149: Mat_Shell *shell = (Mat_Shell*)A->data;
153: if (shell->dshift) { /* get arrays because there is no VecPointwiseMultAdd() */
154: PetscInt i,m;
155: const PetscScalar *x,*d;
156: PetscScalar *y;
157: VecGetLocalSize(X,&m);
158: VecGetArrayRead(shell->dshift,&d);
159: VecGetArrayRead(X,&x);
160: VecGetArray(Y,&y);
161: for (i=0; i<m; i++) y[i] = shell->vscale*y[i] + d[i]*x[i];
162: VecRestoreArrayRead(shell->dshift,&d);
163: VecRestoreArrayRead(X,&x);
164: VecRestoreArray(Y,&y);
165: } else if (PetscAbsScalar(shell->vshift) != 0) {
166: VecAXPBY(Y,shell->vshift,shell->vscale,X);
167: } else if (shell->vscale != 1.0) {
168: VecScale(Y,shell->vscale);
169: }
170: return(0);
171: }
175: /*@
176: MatShellGetContext - Returns the user-provided context associated with a shell matrix.
178: Not Collective
180: Input Parameter:
181: . mat - the matrix, should have been created with MatCreateShell()
183: Output Parameter:
184: . ctx - the user provided context
186: Level: advanced
188: Fortran Notes: To use this from Fortran you must write a Fortran interface definition for this
189: function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.
191: .keywords: matrix, shell, get, context
193: .seealso: MatCreateShell(), MatShellSetOperation(), MatShellSetContext()
194: @*/
195: PetscErrorCode MatShellGetContext(Mat mat,void *ctx)
196: {
198: PetscBool flg;
203: PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
204: if (flg) *(void**)ctx = ((Mat_Shell*)(mat->data))->ctx;
205: else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot get context from non-shell matrix");
206: return(0);
207: }
211: PetscErrorCode MatDestroy_Shell(Mat mat)
212: {
214: Mat_Shell *shell = (Mat_Shell*)mat->data;
217: if (shell->destroy) {
218: (*shell->destroy)(mat);
219: }
220: VecDestroy(&shell->left_owned);
221: VecDestroy(&shell->right_owned);
222: VecDestroy(&shell->dshift_owned);
223: VecDestroy(&shell->left_work);
224: VecDestroy(&shell->right_work);
225: VecDestroy(&shell->left_add_work);
226: VecDestroy(&shell->right_add_work);
227: PetscFree(mat->data);
228: return(0);
229: }
233: PetscErrorCode MatMult_Shell(Mat A,Vec x,Vec y)
234: {
235: Mat_Shell *shell = (Mat_Shell*)A->data;
236: PetscErrorCode ierr;
237: Vec xx;
238: PetscObjectState instate,outstate;
241: MatShellPreScaleRight(A,x,&xx);
242: PetscObjectStateGet((PetscObject)y, &instate);
243: (*shell->mult)(A,xx,y);
244: PetscObjectStateGet((PetscObject)y, &outstate);
245: if (instate == outstate) {
246: /* increase the state of the output vector since the user did not update its state themself as should have been done */
247: PetscObjectStateIncrease((PetscObject)y);
248: }
249: MatShellShiftAndScale(A,xx,y);
250: MatShellPostScaleLeft(A,y);
251: return(0);
252: }
256: PetscErrorCode MatMultAdd_Shell(Mat A,Vec x,Vec y,Vec z)
257: {
258: Mat_Shell *shell = (Mat_Shell*)A->data;
262: if (y == z) {
263: if (!shell->right_add_work) {VecDuplicate(z,&shell->right_add_work);}
264: MatMult(A,x,shell->right_add_work);
265: VecAXPY(z,1.0,shell->right_add_work);
266: } else {
267: MatMult(A,x,z);
268: VecAXPY(z,1.0,y);
269: }
270: return(0);
271: }
275: PetscErrorCode MatMultTranspose_Shell(Mat A,Vec x,Vec y)
276: {
277: Mat_Shell *shell = (Mat_Shell*)A->data;
278: PetscErrorCode ierr;
279: Vec xx;
280: PetscObjectState instate,outstate;
283: MatShellPreScaleLeft(A,x,&xx);
284: PetscObjectStateGet((PetscObject)y, &instate);
285: (*shell->multtranspose)(A,xx,y);
286: PetscObjectStateGet((PetscObject)y, &outstate);
287: if (instate == outstate) {
288: /* increase the state of the output vector since the user did not update its state themself as should have been done */
289: PetscObjectStateIncrease((PetscObject)y);
290: }
291: MatShellShiftAndScale(A,xx,y);
292: MatShellPostScaleRight(A,y);
293: return(0);
294: }
298: PetscErrorCode MatMultTransposeAdd_Shell(Mat A,Vec x,Vec y,Vec z)
299: {
300: Mat_Shell *shell = (Mat_Shell*)A->data;
304: if (y == z) {
305: if (!shell->left_add_work) {VecDuplicate(z,&shell->left_add_work);}
306: MatMultTranspose(A,x,shell->left_add_work);
307: VecWAXPY(z,1.0,shell->left_add_work,y);
308: } else {
309: MatMultTranspose(A,x,z);
310: VecAXPY(z,1.0,y);
311: }
312: return(0);
313: }
317: PetscErrorCode MatGetDiagonal_Shell(Mat A,Vec v)
318: {
319: Mat_Shell *shell = (Mat_Shell*)A->data;
323: (*shell->getdiagonal)(A,v);
324: VecScale(v,shell->vscale);
325: if (shell->dshift) {
326: VecPointwiseMult(v,v,shell->dshift);
327: } else {
328: VecShift(v,shell->vshift);
329: }
330: if (shell->left) {VecPointwiseMult(v,v,shell->left);}
331: if (shell->right) {VecPointwiseMult(v,v,shell->right);}
332: return(0);
333: }
337: PetscErrorCode MatShift_Shell(Mat Y,PetscScalar a)
338: {
339: Mat_Shell *shell = (Mat_Shell*)Y->data;
343: if (shell->left || shell->right || shell->dshift) {
344: if (!shell->dshift) {
345: if (!shell->dshift_owned) {VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift_owned);}
346: shell->dshift = shell->dshift_owned;
347: VecSet(shell->dshift,shell->vshift+a);
348: } else {VecScale(shell->dshift,a);}
349: if (shell->left) {VecPointwiseDivide(shell->dshift,shell->dshift,shell->left);}
350: if (shell->right) {VecPointwiseDivide(shell->dshift,shell->dshift,shell->right);}
351: } else shell->vshift += a;
352: MatShellUseScaledMethods(Y);
353: return(0);
354: }
358: PetscErrorCode MatScale_Shell(Mat Y,PetscScalar a)
359: {
360: Mat_Shell *shell = (Mat_Shell*)Y->data;
364: shell->vscale *= a;
365: if (shell->dshift) {
366: VecScale(shell->dshift,a);
367: } else shell->vshift *= a;
368: MatShellUseScaledMethods(Y);
369: return(0);
370: }
374: static PetscErrorCode MatDiagonalScale_Shell(Mat Y,Vec left,Vec right)
375: {
376: Mat_Shell *shell = (Mat_Shell*)Y->data;
380: if (left) {
381: if (!shell->left_owned) {VecDuplicate(left,&shell->left_owned);}
382: if (shell->left) {
383: VecPointwiseMult(shell->left,shell->left,left);
384: } else {
385: shell->left = shell->left_owned;
386: VecCopy(left,shell->left);
387: }
388: }
389: if (right) {
390: if (!shell->right_owned) {VecDuplicate(right,&shell->right_owned);}
391: if (shell->right) {
392: VecPointwiseMult(shell->right,shell->right,right);
393: } else {
394: shell->right = shell->right_owned;
395: VecCopy(right,shell->right);
396: }
397: }
398: MatShellUseScaledMethods(Y);
399: return(0);
400: }
404: PetscErrorCode MatAssemblyEnd_Shell(Mat Y,MatAssemblyType t)
405: {
406: Mat_Shell *shell = (Mat_Shell*)Y->data;
409: if (t == MAT_FINAL_ASSEMBLY) {
410: shell->vshift = 0.0;
411: shell->vscale = 1.0;
412: shell->dshift = NULL;
413: shell->left = NULL;
414: shell->right = NULL;
415: if (shell->mult) {
416: Y->ops->mult = shell->mult;
417: shell->mult = NULL;
418: }
419: if (shell->multtranspose) {
420: Y->ops->multtranspose = shell->multtranspose;
421: shell->multtranspose = NULL;
422: }
423: if (shell->getdiagonal) {
424: Y->ops->getdiagonal = shell->getdiagonal;
425: shell->getdiagonal = NULL;
426: }
427: shell->usingscaled = PETSC_FALSE;
428: }
429: return(0);
430: }
432: PETSC_INTERN PetscErrorCode MatConvert_Shell(Mat, MatType,MatReuse,Mat*);
436: static PetscErrorCode MatMissingDiagonal_Shell(Mat A,PetscBool *missing,PetscInt *d)
437: {
439: *missing = PETSC_FALSE;
440: return(0);
441: }
443: static struct _MatOps MatOps_Values = {0,
444: 0,
445: 0,
446: 0,
447: /* 4*/ 0,
448: 0,
449: 0,
450: 0,
451: 0,
452: 0,
453: /*10*/ 0,
454: 0,
455: 0,
456: 0,
457: 0,
458: /*15*/ 0,
459: 0,
460: 0,
461: MatDiagonalScale_Shell,
462: 0,
463: /*20*/ 0,
464: MatAssemblyEnd_Shell,
465: 0,
466: 0,
467: /*24*/ 0,
468: 0,
469: 0,
470: 0,
471: 0,
472: /*29*/ 0,
473: 0,
474: 0,
475: 0,
476: 0,
477: /*34*/ 0,
478: 0,
479: 0,
480: 0,
481: 0,
482: /*39*/ 0,
483: 0,
484: 0,
485: 0,
486: 0,
487: /*44*/ 0,
488: MatScale_Shell,
489: MatShift_Shell,
490: 0,
491: 0,
492: /*49*/ 0,
493: 0,
494: 0,
495: 0,
496: 0,
497: /*54*/ 0,
498: 0,
499: 0,
500: 0,
501: 0,
502: /*59*/ 0,
503: MatDestroy_Shell,
504: 0,
505: 0,
506: 0,
507: /*64*/ 0,
508: 0,
509: 0,
510: 0,
511: 0,
512: /*69*/ 0,
513: 0,
514: MatConvert_Shell,
515: 0,
516: 0,
517: /*74*/ 0,
518: 0,
519: 0,
520: 0,
521: 0,
522: /*79*/ 0,
523: 0,
524: 0,
525: 0,
526: 0,
527: /*84*/ 0,
528: 0,
529: 0,
530: 0,
531: 0,
532: /*89*/ 0,
533: 0,
534: 0,
535: 0,
536: 0,
537: /*94*/ 0,
538: 0,
539: 0,
540: 0,
541: 0,
542: /*99*/ 0,
543: 0,
544: 0,
545: 0,
546: 0,
547: /*104*/ 0,
548: 0,
549: 0,
550: 0,
551: 0,
552: /*109*/ 0,
553: 0,
554: 0,
555: 0,
556: MatMissingDiagonal_Shell,
557: /*114*/ 0,
558: 0,
559: 0,
560: 0,
561: 0,
562: /*119*/ 0,
563: 0,
564: 0,
565: 0,
566: 0,
567: /*124*/ 0,
568: 0,
569: 0,
570: 0,
571: 0,
572: /*129*/ 0,
573: 0,
574: 0,
575: 0,
576: 0,
577: /*134*/ 0,
578: 0,
579: 0,
580: 0,
581: 0,
582: /*139*/ 0,
583: 0,
584: 0
585: };
587: /*MC
588: MATSHELL - MATSHELL = "shell" - A matrix type to be used to define your own matrix type -- perhaps matrix free.
590: Level: advanced
592: .seealso: MatCreateShell
593: M*/
597: PETSC_EXTERN PetscErrorCode MatCreate_Shell(Mat A)
598: {
599: Mat_Shell *b;
603: PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));
605: PetscNewLog(A,&b);
606: A->data = (void*)b;
608: PetscLayoutSetUp(A->rmap);
609: PetscLayoutSetUp(A->cmap);
611: b->ctx = 0;
612: b->vshift = 0.0;
613: b->vscale = 1.0;
614: b->mult = 0;
615: b->multtranspose = 0;
616: b->getdiagonal = 0;
617: A->assembled = PETSC_TRUE;
618: A->preallocated = PETSC_FALSE;
620: PetscObjectChangeTypeName((PetscObject)A,MATSHELL);
621: return(0);
622: }
626: /*@C
627: MatCreateShell - Creates a new matrix class for use with a user-defined
628: private data storage format.
630: Collective on MPI_Comm
632: Input Parameters:
633: + comm - MPI communicator
634: . m - number of local rows (must be given)
635: . n - number of local columns (must be given)
636: . M - number of global rows (may be PETSC_DETERMINE)
637: . N - number of global columns (may be PETSC_DETERMINE)
638: - ctx - pointer to data needed by the shell matrix routines
640: Output Parameter:
641: . A - the matrix
643: Level: advanced
645: Usage:
646: $ extern int mult(Mat,Vec,Vec);
647: $ MatCreateShell(comm,m,n,M,N,ctx,&mat);
648: $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
649: $ [ Use matrix for operations that have been set ]
650: $ MatDestroy(mat);
652: Notes:
653: The shell matrix type is intended to provide a simple class to use
654: with KSP (such as, for use with matrix-free methods). You should not
655: use the shell type if you plan to define a complete matrix class.
657: Fortran Notes: To use this from Fortran with a ctx you must write an interface definition for this
658: function and for MatShellGetContext() that tells Fortran the Fortran derived data type you are passing
659: in as the ctx argument.
661: PETSc requires that matrices and vectors being used for certain
662: operations are partitioned accordingly. For example, when
663: creating a shell matrix, A, that supports parallel matrix-vector
664: products using MatMult(A,x,y) the user should set the number
665: of local matrix rows to be the number of local elements of the
666: corresponding result vector, y. Note that this is information is
667: required for use of the matrix interface routines, even though
668: the shell matrix may not actually be physically partitioned.
669: For example,
671: $
672: $ Vec x, y
673: $ extern int mult(Mat,Vec,Vec);
674: $ Mat A
675: $
676: $ VecCreateMPI(comm,PETSC_DECIDE,M,&y);
677: $ VecCreateMPI(comm,PETSC_DECIDE,N,&x);
678: $ VecGetLocalSize(y,&m);
679: $ VecGetLocalSize(x,&n);
680: $ MatCreateShell(comm,m,n,M,N,ctx,&A);
681: $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
682: $ MatMult(A,x,y);
683: $ MatDestroy(A);
684: $ VecDestroy(y); VecDestroy(x);
685: $
687: .keywords: matrix, shell, create
689: .seealso: MatShellSetOperation(), MatHasOperation(), MatShellGetContext(), MatShellSetContext()
690: @*/
691: PetscErrorCode MatCreateShell(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,void *ctx,Mat *A)
692: {
696: MatCreate(comm,A);
697: MatSetSizes(*A,m,n,M,N);
698: MatSetType(*A,MATSHELL);
699: MatShellSetContext(*A,ctx);
700: MatSetUp(*A);
701: return(0);
702: }
706: /*@
707: MatShellSetContext - sets the context for a shell matrix
709: Logically Collective on Mat
711: Input Parameters:
712: + mat - the shell matrix
713: - ctx - the context
715: Level: advanced
717: Fortran Notes: To use this from Fortran you must write a Fortran interface definition for this
718: function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.
720: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation()
721: @*/
722: PetscErrorCode MatShellSetContext(Mat mat,void *ctx)
723: {
724: Mat_Shell *shell = (Mat_Shell*)mat->data;
726: PetscBool flg;
730: PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
731: if (flg) {
732: shell->ctx = ctx;
733: } else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot attach context to non-shell matrix");
734: return(0);
735: }
739: /*@C
740: MatShellSetOperation - Allows user to set a matrix operation for
741: a shell matrix.
743: Logically Collective on Mat
745: Input Parameters:
746: + mat - the shell matrix
747: . op - the name of the operation
748: - f - the function that provides the operation.
750: Level: advanced
752: Usage:
753: $ extern PetscErrorCode usermult(Mat,Vec,Vec);
754: $ MatCreateShell(comm,m,n,M,N,ctx,&A);
755: $ MatShellSetOperation(A,MATOP_MULT,(void(*)(void))usermult);
757: Notes:
758: See the file include/petscmat.h for a complete list of matrix
759: operations, which all have the form MATOP_<OPERATION>, where
760: <OPERATION> is the name (in all capital letters) of the
761: user interface routine (e.g., MatMult() -> MATOP_MULT).
763: All user-provided functions (execept for MATOP_DESTROY) should have the same calling
764: sequence as the usual matrix interface routines, since they
765: are intended to be accessed via the usual matrix interface
766: routines, e.g.,
767: $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec)
769: In particular each function MUST return an error code of 0 on success and
770: nonzero on failure.
772: Within each user-defined routine, the user should call
773: MatShellGetContext() to obtain the user-defined context that was
774: set by MatCreateShell().
776: Fortran Notes: For MatCreateVecs() the user code should check if the input left or right matrix is -1 and in that case not
777: generate a matrix. See src/mat/examples/tests/ex120f.F
779: .keywords: matrix, shell, set, operation
781: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext()
782: @*/
783: PetscErrorCode MatShellSetOperation(Mat mat,MatOperation op,void (*f)(void))
784: {
786: PetscBool flg;
790: switch (op) {
791: case MATOP_DESTROY:
792: PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
793: if (flg) {
794: Mat_Shell *shell = (Mat_Shell*)mat->data;
795: shell->destroy = (PetscErrorCode (*)(Mat))f;
796: } else mat->ops->destroy = (PetscErrorCode (*)(Mat))f;
797: break;
798: case MATOP_VIEW:
799: mat->ops->view = (PetscErrorCode (*)(Mat,PetscViewer))f;
800: break;
801: case MATOP_MULT:
802: mat->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f;
803: if (!mat->ops->multadd) mat->ops->multadd = MatMultAdd_Shell;
804: break;
805: case MATOP_MULT_TRANSPOSE:
806: mat->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f;
807: if (!mat->ops->multtransposeadd) mat->ops->multtransposeadd = MatMultTransposeAdd_Shell;
808: break;
809: default:
810: (((void(**)(void))mat->ops)[op]) = f;
811: }
812: return(0);
813: }
817: /*@C
818: MatShellGetOperation - Gets a matrix function for a shell matrix.
820: Not Collective
822: Input Parameters:
823: + mat - the shell matrix
824: - op - the name of the operation
826: Output Parameter:
827: . f - the function that provides the operation.
829: Level: advanced
831: Notes:
832: See the file include/petscmat.h for a complete list of matrix
833: operations, which all have the form MATOP_<OPERATION>, where
834: <OPERATION> is the name (in all capital letters) of the
835: user interface routine (e.g., MatMult() -> MATOP_MULT).
837: All user-provided functions have the same calling
838: sequence as the usual matrix interface routines, since they
839: are intended to be accessed via the usual matrix interface
840: routines, e.g.,
841: $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec)
843: Within each user-defined routine, the user should call
844: MatShellGetContext() to obtain the user-defined context that was
845: set by MatCreateShell().
847: .keywords: matrix, shell, set, operation
849: .seealso: MatCreateShell(), MatShellGetContext(), MatShellSetOperation(), MatShellSetContext()
850: @*/
851: PetscErrorCode MatShellGetOperation(Mat mat,MatOperation op,void(**f)(void))
852: {
854: PetscBool flg;
858: if (op == MATOP_DESTROY) {
859: PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
860: if (flg) {
861: Mat_Shell *shell = (Mat_Shell*)mat->data;
862: *f = (void (*)(void))shell->destroy;
863: } else {
864: *f = (void (*)(void))mat->ops->destroy;
865: }
866: } else if (op == MATOP_VIEW) {
867: *f = (void (*)(void))mat->ops->view;
868: } else {
869: *f = (((void (**)(void))mat->ops)[op]);
870: }
871: return(0);
872: }