Actual source code: shell.c
petsc-3.8.4 2018-03-24
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: /* 0 */
12: PetscErrorCode (*mult)(Mat,Vec,Vec);
13: /* 5 */
14: PetscErrorCode (*multtranspose)(Mat,Vec,Vec);
15: /* 10 */
16: /* 15 */
17: PetscErrorCode (*getdiagonal)(Mat,Vec);
18: /* 20 */
19: /* 24 */
20: /* 29 */
21: /* 34 */
22: /* 39 */
23: PetscErrorCode (*copy)(Mat,Mat,MatStructure);
24: /* 44 */
25: PetscErrorCode (*diagonalset)(Mat,Vec,InsertMode);
26: /* 49 */
27: /* 54 */
28: /* 59 */
29: PetscErrorCode (*destroy)(Mat);
30: /* 64 */
31: /* 69 */
32: /* 74 */
33: /* 79 */
34: /* 84 */
35: /* 89 */
36: /* 94 */
37: /* 99 */
38: /* 104 */
39: /* 109 */
40: /* 114 */
41: /* 119 */
42: /* 124 */
43: /* 129 */
44: /* 134 */
45: /* 139 */
46: /* 144 */
47: };
49: typedef struct {
50: struct _MatShellOps ops[1];
52: PetscScalar vscale,vshift;
53: Vec dshift;
54: Vec left,right;
55: Vec dshift_owned,left_owned,right_owned;
56: Vec left_work,right_work;
57: Vec left_add_work,right_add_work;
58: PetscBool usingscaled;
59: void *ctx;
60: } Mat_Shell;
61: /*
62: The most general expression for the matrix is
64: S = L (a A + B) R
66: where
67: A is the matrix defined by the user's function
68: a is a scalar multiple
69: L is left scaling
70: R is right scaling
71: B is a diagonal shift defined by
72: diag(dshift) if the vector dshift is non-NULL
73: vscale*identity otherwise
75: The following identities apply:
77: Scale by c:
78: c [L (a A + B) R] = L [(a c) A + c B] R
80: Shift by c:
81: [L (a A + B) R] + c = L [a A + (B + c Linv Rinv)] R
83: Diagonal scaling is achieved by simply multiplying with existing L and R vectors
85: In the data structure:
87: vscale=1.0 means no special scaling will be applied
88: dshift=NULL means a constant diagonal shift (fall back to vshift)
89: vshift=0.0 means no constant diagonal shift, note that vshift is only used if dshift is NULL
90: */
92: static PetscErrorCode MatMult_Shell(Mat,Vec,Vec);
93: static PetscErrorCode MatMultTranspose_Shell(Mat,Vec,Vec);
94: static PetscErrorCode MatGetDiagonal_Shell(Mat,Vec);
95: static PetscErrorCode MatCopy_Shell(Mat,Mat,MatStructure);
97: static PetscErrorCode MatShellUseScaledMethods(Mat Y)
98: {
99: Mat_Shell *shell = (Mat_Shell*)Y->data;
102: if (shell->usingscaled) return(0);
103: shell->ops->mult = Y->ops->mult;
104: Y->ops->mult = MatMult_Shell;
105: if (Y->ops->multtranspose) {
106: shell->ops->multtranspose = Y->ops->multtranspose;
107: Y->ops->multtranspose = MatMultTranspose_Shell;
108: }
109: if (Y->ops->getdiagonal) {
110: shell->ops->getdiagonal = Y->ops->getdiagonal;
111: Y->ops->getdiagonal = MatGetDiagonal_Shell;
112: }
113: if (Y->ops->copy) {
114: shell->ops->copy = Y->ops->copy;
115: Y->ops->copy = MatCopy_Shell;
116: }
117: shell->usingscaled = PETSC_TRUE;
118: return(0);
119: }
121: static PetscErrorCode MatShellPreScaleLeft(Mat A,Vec x,Vec *xx)
122: {
123: Mat_Shell *shell = (Mat_Shell*)A->data;
127: *xx = NULL;
128: if (!shell->left) {
129: *xx = x;
130: } else {
131: if (!shell->left_work) {VecDuplicate(shell->left,&shell->left_work);}
132: VecPointwiseMult(shell->left_work,x,shell->left);
133: *xx = shell->left_work;
134: }
135: return(0);
136: }
138: static PetscErrorCode MatShellPreScaleRight(Mat A,Vec x,Vec *xx)
139: {
140: Mat_Shell *shell = (Mat_Shell*)A->data;
144: *xx = NULL;
145: if (!shell->right) {
146: *xx = x;
147: } else {
148: if (!shell->right_work) {VecDuplicate(shell->right,&shell->right_work);}
149: VecPointwiseMult(shell->right_work,x,shell->right);
150: *xx = shell->right_work;
151: }
152: return(0);
153: }
155: static PetscErrorCode MatShellPostScaleLeft(Mat A,Vec x)
156: {
157: Mat_Shell *shell = (Mat_Shell*)A->data;
161: if (shell->left) {VecPointwiseMult(x,x,shell->left);}
162: return(0);
163: }
165: static PetscErrorCode MatShellPostScaleRight(Mat A,Vec x)
166: {
167: Mat_Shell *shell = (Mat_Shell*)A->data;
171: if (shell->right) {VecPointwiseMult(x,x,shell->right);}
172: return(0);
173: }
175: static PetscErrorCode MatShellShiftAndScale(Mat A,Vec X,Vec Y)
176: {
177: Mat_Shell *shell = (Mat_Shell*)A->data;
181: if (shell->dshift) { /* get arrays because there is no VecPointwiseMultAdd() */
182: PetscInt i,m;
183: const PetscScalar *x,*d;
184: PetscScalar *y;
185: VecGetLocalSize(X,&m);
186: VecGetArrayRead(shell->dshift,&d);
187: VecGetArrayRead(X,&x);
188: VecGetArray(Y,&y);
189: for (i=0; i<m; i++) y[i] = shell->vscale*y[i] + d[i]*x[i];
190: VecRestoreArrayRead(shell->dshift,&d);
191: VecRestoreArrayRead(X,&x);
192: VecRestoreArray(Y,&y);
193: } else if (PetscAbsScalar(shell->vshift) != 0) {
194: VecAXPBY(Y,shell->vshift,shell->vscale,X);
195: } else if (shell->vscale != 1.0) {
196: VecScale(Y,shell->vscale);
197: }
198: return(0);
199: }
201: /*@
202: MatShellGetContext - Returns the user-provided context associated with a shell matrix.
204: Not Collective
206: Input Parameter:
207: . mat - the matrix, should have been created with MatCreateShell()
209: Output Parameter:
210: . ctx - the user provided context
212: Level: advanced
214: Fortran Notes: To use this from Fortran you must write a Fortran interface definition for this
215: function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.
217: .keywords: matrix, shell, get, context
219: .seealso: MatCreateShell(), MatShellSetOperation(), MatShellSetContext()
220: @*/
221: PetscErrorCode MatShellGetContext(Mat mat,void *ctx)
222: {
224: PetscBool flg;
229: PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
230: if (flg) *(void**)ctx = ((Mat_Shell*)(mat->data))->ctx;
231: else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot get context from non-shell matrix");
232: return(0);
233: }
235: PetscErrorCode MatDestroy_Shell(Mat mat)
236: {
238: Mat_Shell *shell = (Mat_Shell*)mat->data;
241: if (shell->ops->destroy) {
242: (*shell->ops->destroy)(mat);
243: }
244: VecDestroy(&shell->left_owned);
245: VecDestroy(&shell->right_owned);
246: VecDestroy(&shell->dshift_owned);
247: VecDestroy(&shell->left_work);
248: VecDestroy(&shell->right_work);
249: VecDestroy(&shell->left_add_work);
250: VecDestroy(&shell->right_add_work);
251: PetscFree(mat->data);
252: return(0);
253: }
255: PetscErrorCode MatCopy_Shell(Mat A,Mat B,MatStructure str)
256: {
257: Mat_Shell *shellA = (Mat_Shell*)A->data,*shellB = (Mat_Shell*)B->data;
258: PetscBool matflg,shellflg;
259: PetscErrorCode ierr;
262: PetscObjectTypeCompare((PetscObject)B,MATSHELL,&matflg);
263: if(!matflg) { SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_NOTSAMETYPE,"Output matrix must be a MATSHELL"); }
265: MatShellUseScaledMethods(B);
266: PetscMemcmp(A->ops,B->ops,sizeof(struct _MatOps),&matflg);
267: PetscMemcmp(shellA->ops,shellB->ops,sizeof(struct _MatShellOps),&shellflg);
268: if (!matflg || !shellflg) { SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Matrices cannot be copied because they have different operations"); }
270: (*shellA->ops->copy)(A,B,str);
271: shellB->vscale = shellA->vscale;
272: shellB->vshift = shellA->vshift;
273: if (shellA->dshift_owned) {
274: if (!shellB->dshift_owned) {
275: VecDuplicate(shellA->dshift_owned,&shellB->dshift_owned);
276: }
277: VecCopy(shellA->dshift_owned,shellB->dshift_owned);
278: shellB->dshift = shellB->dshift_owned;
279: } else {
280: shellB->dshift = NULL;
281: }
282: if (shellA->left_owned) {
283: if (!shellB->left_owned) {
284: VecDuplicate(shellA->left_owned,&shellB->left_owned);
285: }
286: VecCopy(shellA->left_owned,shellB->left_owned);
287: shellB->left = shellB->left_owned;
288: } else {
289: shellB->left = NULL;
290: }
291: if (shellA->right_owned) {
292: if (!shellB->right_owned) {
293: VecDuplicate(shellA->right_owned,&shellB->right_owned);
294: }
295: VecCopy(shellA->right_owned,shellB->right_owned);
296: shellB->right = shellB->right_owned;
297: } else {
298: shellB->right = NULL;
299: }
300: return(0);
301: }
303: PetscErrorCode MatMult_Shell(Mat A,Vec x,Vec y)
304: {
305: Mat_Shell *shell = (Mat_Shell*)A->data;
306: PetscErrorCode ierr;
307: Vec xx;
308: PetscObjectState instate,outstate;
311: MatShellPreScaleRight(A,x,&xx);
312: PetscObjectStateGet((PetscObject)y, &instate);
313: (*shell->ops->mult)(A,xx,y);
314: PetscObjectStateGet((PetscObject)y, &outstate);
315: if (instate == outstate) {
316: /* increase the state of the output vector since the user did not update its state themself as should have been done */
317: PetscObjectStateIncrease((PetscObject)y);
318: }
319: MatShellShiftAndScale(A,xx,y);
320: MatShellPostScaleLeft(A,y);
321: return(0);
322: }
324: PetscErrorCode MatMultAdd_Shell(Mat A,Vec x,Vec y,Vec z)
325: {
326: Mat_Shell *shell = (Mat_Shell*)A->data;
330: if (y == z) {
331: if (!shell->right_add_work) {VecDuplicate(z,&shell->right_add_work);}
332: MatMult(A,x,shell->right_add_work);
333: VecAXPY(z,1.0,shell->right_add_work);
334: } else {
335: MatMult(A,x,z);
336: VecAXPY(z,1.0,y);
337: }
338: return(0);
339: }
341: PetscErrorCode MatMultTranspose_Shell(Mat A,Vec x,Vec y)
342: {
343: Mat_Shell *shell = (Mat_Shell*)A->data;
344: PetscErrorCode ierr;
345: Vec xx;
346: PetscObjectState instate,outstate;
349: MatShellPreScaleLeft(A,x,&xx);
350: PetscObjectStateGet((PetscObject)y, &instate);
351: (*shell->ops->multtranspose)(A,xx,y);
352: PetscObjectStateGet((PetscObject)y, &outstate);
353: if (instate == outstate) {
354: /* increase the state of the output vector since the user did not update its state themself as should have been done */
355: PetscObjectStateIncrease((PetscObject)y);
356: }
357: MatShellShiftAndScale(A,xx,y);
358: MatShellPostScaleRight(A,y);
359: return(0);
360: }
362: PetscErrorCode MatMultTransposeAdd_Shell(Mat A,Vec x,Vec y,Vec z)
363: {
364: Mat_Shell *shell = (Mat_Shell*)A->data;
368: if (y == z) {
369: if (!shell->left_add_work) {VecDuplicate(z,&shell->left_add_work);}
370: MatMultTranspose(A,x,shell->left_add_work);
371: VecWAXPY(z,1.0,shell->left_add_work,y);
372: } else {
373: MatMultTranspose(A,x,z);
374: VecAXPY(z,1.0,y);
375: }
376: return(0);
377: }
379: PetscErrorCode MatGetDiagonal_Shell(Mat A,Vec v)
380: {
381: Mat_Shell *shell = (Mat_Shell*)A->data;
385: (*shell->ops->getdiagonal)(A,v);
386: VecScale(v,shell->vscale);
387: if (shell->dshift) {
388: VecPointwiseMult(v,v,shell->dshift);
389: } else {
390: VecShift(v,shell->vshift);
391: }
392: if (shell->left) {VecPointwiseMult(v,v,shell->left);}
393: if (shell->right) {VecPointwiseMult(v,v,shell->right);}
394: return(0);
395: }
397: PetscErrorCode MatShift_Shell(Mat Y,PetscScalar a)
398: {
399: Mat_Shell *shell = (Mat_Shell*)Y->data;
403: if (shell->left || shell->right || shell->dshift) {
404: if (!shell->dshift) {
405: if (!shell->dshift_owned) {VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift_owned);}
406: shell->dshift = shell->dshift_owned;
407: VecSet(shell->dshift,shell->vshift+a);
408: } else {VecScale(shell->dshift,a);}
409: if (shell->left) {VecPointwiseDivide(shell->dshift,shell->dshift,shell->left);}
410: if (shell->right) {VecPointwiseDivide(shell->dshift,shell->dshift,shell->right);}
411: } else shell->vshift += a;
412: MatShellUseScaledMethods(Y);
413: return(0);
414: }
416: PetscErrorCode MatDiagonalSet_Shell(Mat A,Vec D,InsertMode ins)
417: {
418: Mat_Shell *shell = (Mat_Shell*)A->data;
422: if (shell->vscale != 1.0 || shell->left || shell->right) {
423: SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE, "Operation not supported on scaled matrix");
424: }
426: if (shell->ops->diagonalset) {
427: (*shell->ops->diagonalset)(A,D,ins);
428: }
429: shell->vshift = 0.0;
430: if (shell->dshift) { VecZeroEntries(shell->dshift); }
431: return(0);
432: }
434: PetscErrorCode MatScale_Shell(Mat Y,PetscScalar a)
435: {
436: Mat_Shell *shell = (Mat_Shell*)Y->data;
440: shell->vscale *= a;
441: if (shell->dshift) {
442: VecScale(shell->dshift,a);
443: } else shell->vshift *= a;
444: MatShellUseScaledMethods(Y);
445: return(0);
446: }
448: static PetscErrorCode MatDiagonalScale_Shell(Mat Y,Vec left,Vec right)
449: {
450: Mat_Shell *shell = (Mat_Shell*)Y->data;
454: if (left) {
455: if (!shell->left_owned) {VecDuplicate(left,&shell->left_owned);}
456: if (shell->left) {
457: VecPointwiseMult(shell->left,shell->left,left);
458: } else {
459: shell->left = shell->left_owned;
460: VecCopy(left,shell->left);
461: }
462: }
463: if (right) {
464: if (!shell->right_owned) {VecDuplicate(right,&shell->right_owned);}
465: if (shell->right) {
466: VecPointwiseMult(shell->right,shell->right,right);
467: } else {
468: shell->right = shell->right_owned;
469: VecCopy(right,shell->right);
470: }
471: }
472: MatShellUseScaledMethods(Y);
473: return(0);
474: }
476: PetscErrorCode MatAssemblyEnd_Shell(Mat Y,MatAssemblyType t)
477: {
478: Mat_Shell *shell = (Mat_Shell*)Y->data;
481: if (t == MAT_FINAL_ASSEMBLY) {
482: shell->vshift = 0.0;
483: shell->vscale = 1.0;
484: shell->dshift = NULL;
485: shell->left = NULL;
486: shell->right = NULL;
487: if (shell->ops->mult) {
488: Y->ops->mult = shell->ops->mult;
489: shell->ops->mult = NULL;
490: }
491: if (shell->ops->multtranspose) {
492: Y->ops->multtranspose = shell->ops->multtranspose;
493: shell->ops->multtranspose = NULL;
494: }
495: if (shell->ops->getdiagonal) {
496: Y->ops->getdiagonal = shell->ops->getdiagonal;
497: shell->ops->getdiagonal = NULL;
498: }
499: if (shell->ops->copy) {
500: Y->ops->copy = shell->ops->copy;
501: shell->ops->copy = NULL;
502: }
503: shell->usingscaled = PETSC_FALSE;
504: }
505: return(0);
506: }
508: PETSC_INTERN PetscErrorCode MatConvert_Shell(Mat, MatType,MatReuse,Mat*);
510: static PetscErrorCode MatMissingDiagonal_Shell(Mat A,PetscBool *missing,PetscInt *d)
511: {
513: *missing = PETSC_FALSE;
514: return(0);
515: }
517: static struct _MatOps MatOps_Values = {0,
518: 0,
519: 0,
520: 0,
521: /* 4*/ 0,
522: 0,
523: 0,
524: 0,
525: 0,
526: 0,
527: /*10*/ 0,
528: 0,
529: 0,
530: 0,
531: 0,
532: /*15*/ 0,
533: 0,
534: 0,
535: MatDiagonalScale_Shell,
536: 0,
537: /*20*/ 0,
538: MatAssemblyEnd_Shell,
539: 0,
540: 0,
541: /*24*/ 0,
542: 0,
543: 0,
544: 0,
545: 0,
546: /*29*/ 0,
547: 0,
548: 0,
549: 0,
550: 0,
551: /*34*/ 0,
552: 0,
553: 0,
554: 0,
555: 0,
556: /*39*/ 0,
557: 0,
558: 0,
559: 0,
560: 0,
561: /*44*/ 0,
562: MatScale_Shell,
563: MatShift_Shell,
564: MatDiagonalSet_Shell,
565: 0,
566: /*49*/ 0,
567: 0,
568: 0,
569: 0,
570: 0,
571: /*54*/ 0,
572: 0,
573: 0,
574: 0,
575: 0,
576: /*59*/ 0,
577: MatDestroy_Shell,
578: 0,
579: 0,
580: 0,
581: /*64*/ 0,
582: 0,
583: 0,
584: 0,
585: 0,
586: /*69*/ 0,
587: 0,
588: MatConvert_Shell,
589: 0,
590: 0,
591: /*74*/ 0,
592: 0,
593: 0,
594: 0,
595: 0,
596: /*79*/ 0,
597: 0,
598: 0,
599: 0,
600: 0,
601: /*84*/ 0,
602: 0,
603: 0,
604: 0,
605: 0,
606: /*89*/ 0,
607: 0,
608: 0,
609: 0,
610: 0,
611: /*94*/ 0,
612: 0,
613: 0,
614: 0,
615: 0,
616: /*99*/ 0,
617: 0,
618: 0,
619: 0,
620: 0,
621: /*104*/ 0,
622: 0,
623: 0,
624: 0,
625: 0,
626: /*109*/ 0,
627: 0,
628: 0,
629: 0,
630: MatMissingDiagonal_Shell,
631: /*114*/ 0,
632: 0,
633: 0,
634: 0,
635: 0,
636: /*119*/ 0,
637: 0,
638: 0,
639: 0,
640: 0,
641: /*124*/ 0,
642: 0,
643: 0,
644: 0,
645: 0,
646: /*129*/ 0,
647: 0,
648: 0,
649: 0,
650: 0,
651: /*134*/ 0,
652: 0,
653: 0,
654: 0,
655: 0,
656: /*139*/ 0,
657: 0,
658: 0
659: };
661: /*MC
662: MATSHELL - MATSHELL = "shell" - A matrix type to be used to define your own matrix type -- perhaps matrix free.
664: Level: advanced
666: .seealso: MatCreateShell
667: M*/
669: PETSC_EXTERN PetscErrorCode MatCreate_Shell(Mat A)
670: {
671: Mat_Shell *b;
675: PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));
677: PetscNewLog(A,&b);
678: A->data = (void*)b;
680: PetscLayoutSetUp(A->rmap);
681: PetscLayoutSetUp(A->cmap);
683: b->ctx = 0;
684: b->vshift = 0.0;
685: b->vscale = 1.0;
686: A->assembled = PETSC_TRUE;
687: A->preallocated = PETSC_FALSE;
689: PetscObjectChangeTypeName((PetscObject)A,MATSHELL);
690: return(0);
691: }
693: /*@C
694: MatCreateShell - Creates a new matrix class for use with a user-defined
695: private data storage format.
697: Collective on MPI_Comm
699: Input Parameters:
700: + comm - MPI communicator
701: . m - number of local rows (must be given)
702: . n - number of local columns (must be given)
703: . M - number of global rows (may be PETSC_DETERMINE)
704: . N - number of global columns (may be PETSC_DETERMINE)
705: - ctx - pointer to data needed by the shell matrix routines
707: Output Parameter:
708: . A - the matrix
710: Level: advanced
712: Usage:
713: $ extern int mult(Mat,Vec,Vec);
714: $ MatCreateShell(comm,m,n,M,N,ctx,&mat);
715: $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
716: $ [ Use matrix for operations that have been set ]
717: $ MatDestroy(mat);
719: Notes:
720: The shell matrix type is intended to provide a simple class to use
721: with KSP (such as, for use with matrix-free methods). You should not
722: use the shell type if you plan to define a complete matrix class.
724: Fortran Notes: To use this from Fortran with a ctx you must write an interface definition for this
725: function and for MatShellGetContext() that tells Fortran the Fortran derived data type you are passing
726: in as the ctx argument.
728: PETSc requires that matrices and vectors being used for certain
729: operations are partitioned accordingly. For example, when
730: creating a shell matrix, A, that supports parallel matrix-vector
731: products using MatMult(A,x,y) the user should set the number
732: of local matrix rows to be the number of local elements of the
733: corresponding result vector, y. Note that this is information is
734: required for use of the matrix interface routines, even though
735: the shell matrix may not actually be physically partitioned.
736: For example,
738: $
739: $ Vec x, y
740: $ extern int mult(Mat,Vec,Vec);
741: $ Mat A
742: $
743: $ VecCreateMPI(comm,PETSC_DECIDE,M,&y);
744: $ VecCreateMPI(comm,PETSC_DECIDE,N,&x);
745: $ VecGetLocalSize(y,&m);
746: $ VecGetLocalSize(x,&n);
747: $ MatCreateShell(comm,m,n,M,N,ctx,&A);
748: $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
749: $ MatMult(A,x,y);
750: $ MatDestroy(A);
751: $ VecDestroy(y); VecDestroy(x);
752: $
754: .keywords: matrix, shell, create
756: .seealso: MatShellSetOperation(), MatHasOperation(), MatShellGetContext(), MatShellSetContext()
757: @*/
758: PetscErrorCode MatCreateShell(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,void *ctx,Mat *A)
759: {
763: MatCreate(comm,A);
764: MatSetSizes(*A,m,n,M,N);
765: MatSetType(*A,MATSHELL);
766: MatShellSetContext(*A,ctx);
767: MatSetUp(*A);
768: return(0);
769: }
771: /*@
772: MatShellSetContext - sets the context for a shell matrix
774: Logically Collective on Mat
776: Input Parameters:
777: + mat - the shell matrix
778: - ctx - the context
780: Level: advanced
782: Fortran Notes: To use this from Fortran you must write a Fortran interface definition for this
783: function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.
785: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation()
786: @*/
787: PetscErrorCode MatShellSetContext(Mat mat,void *ctx)
788: {
789: Mat_Shell *shell = (Mat_Shell*)mat->data;
791: PetscBool flg;
795: PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
796: if (flg) {
797: shell->ctx = ctx;
798: } else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot attach context to non-shell matrix");
799: return(0);
800: }
802: /*@C
803: MatShellSetOperation - Allows user to set a matrix operation for
804: a shell matrix.
806: Logically Collective on Mat
808: Input Parameters:
809: + mat - the shell matrix
810: . op - the name of the operation
811: - f - the function that provides the operation.
813: Level: advanced
815: Usage:
816: $ extern PetscErrorCode usermult(Mat,Vec,Vec);
817: $ MatCreateShell(comm,m,n,M,N,ctx,&A);
818: $ MatShellSetOperation(A,MATOP_MULT,(void(*)(void))usermult);
820: Notes:
821: See the file include/petscmat.h for a complete list of matrix
822: operations, which all have the form MATOP_<OPERATION>, where
823: <OPERATION> is the name (in all capital letters) of the
824: user interface routine (e.g., MatMult() -> MATOP_MULT).
826: All user-provided functions (execept for MATOP_DESTROY) should have the same calling
827: sequence as the usual matrix interface routines, since they
828: are intended to be accessed via the usual matrix interface
829: routines, e.g.,
830: $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec)
832: In particular each function MUST return an error code of 0 on success and
833: nonzero on failure.
835: Within each user-defined routine, the user should call
836: MatShellGetContext() to obtain the user-defined context that was
837: set by MatCreateShell().
839: Fortran Notes: For MatCreateVecs() the user code should check if the input left or right matrix is -1 and in that case not
840: generate a matrix. See src/mat/examples/tests/ex120f.F
842: .keywords: matrix, shell, set, operation
844: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext()
845: @*/
846: PetscErrorCode MatShellSetOperation(Mat mat,MatOperation op,void (*f)(void))
847: {
849: PetscBool flg;
853: switch (op) {
854: case MATOP_DESTROY:
855: PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
856: if (flg) {
857: Mat_Shell *shell = (Mat_Shell*)mat->data;
858: shell->ops->destroy = (PetscErrorCode (*)(Mat))f;
859: } else mat->ops->destroy = (PetscErrorCode (*)(Mat))f;
860: break;
861: case MATOP_DIAGONAL_SET:
862: PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
863: if (flg) {
864: Mat_Shell *shell = (Mat_Shell*)mat->data;
865: shell->ops->diagonalset = (PetscErrorCode (*)(Mat,Vec,InsertMode))f;
866: } else {
867: mat->ops->diagonalset = (PetscErrorCode (*)(Mat,Vec,InsertMode))f;
868: }
869: break;
870: case MATOP_VIEW:
871: if (!mat->ops->viewnative) {
872: mat->ops->viewnative = mat->ops->view;
873: }
874: mat->ops->view = (PetscErrorCode (*)(Mat,PetscViewer))f;
875: break;
876: case MATOP_MULT:
877: mat->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f;
878: if (!mat->ops->multadd) {
879: PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
880: if (flg) mat->ops->multadd = MatMultAdd_Shell;
881: }
882: break;
883: case MATOP_MULT_TRANSPOSE:
884: mat->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f;
885: if (!mat->ops->multtransposeadd) {
886: PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
887: if (flg) mat->ops->multtransposeadd = MatMultTransposeAdd_Shell;
888: }
889: break;
890: default:
891: (((void(**)(void))mat->ops)[op]) = f;
892: }
893: return(0);
894: }
896: /*@C
897: MatShellGetOperation - Gets a matrix function for a shell matrix.
899: Not Collective
901: Input Parameters:
902: + mat - the shell matrix
903: - op - the name of the operation
905: Output Parameter:
906: . f - the function that provides the operation.
908: Level: advanced
910: Notes:
911: See the file include/petscmat.h for a complete list of matrix
912: operations, which all have the form MATOP_<OPERATION>, where
913: <OPERATION> is the name (in all capital letters) of the
914: user interface routine (e.g., MatMult() -> MATOP_MULT).
916: All user-provided functions have the same calling
917: sequence as the usual matrix interface routines, since they
918: are intended to be accessed via the usual matrix interface
919: routines, e.g.,
920: $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec)
922: Within each user-defined routine, the user should call
923: MatShellGetContext() to obtain the user-defined context that was
924: set by MatCreateShell().
926: .keywords: matrix, shell, set, operation
928: .seealso: MatCreateShell(), MatShellGetContext(), MatShellSetOperation(), MatShellSetContext()
929: @*/
930: PetscErrorCode MatShellGetOperation(Mat mat,MatOperation op,void(**f)(void))
931: {
933: PetscBool flg;
937: switch (op) {
938: case MATOP_DESTROY:
939: PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
940: if (flg) {
941: Mat_Shell *shell = (Mat_Shell*)mat->data;
942: *f = (void (*)(void))shell->ops->destroy;
943: } else {
944: *f = (void (*)(void))mat->ops->destroy;
945: }
946: break;
947: case MATOP_DIAGONAL_SET:
948: PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
949: if (flg) {
950: Mat_Shell *shell = (Mat_Shell*)mat->data;
951: *f = (void (*)(void))shell->ops->diagonalset;
952: } else {
953: *f = (void (*)(void))mat->ops->diagonalset;
954: }
955: break;
956: case MATOP_VIEW:
957: *f = (void (*)(void))mat->ops->view;
958: break;
959: default:
960: *f = (((void (**)(void))mat->ops)[op]);
961: }
962: return(0);
963: }