Actual source code: shell.c
petsc-3.11.4 2019-09-28
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:
142: To use this from Fortran you must write a Fortran interface definition for this
143: function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.
145: .keywords: matrix, shell, get, context
147: .seealso: MatCreateShell(), MatShellSetOperation(), MatShellSetContext()
148: @*/
149: PetscErrorCode MatShellGetContext(Mat mat,void *ctx)
150: {
152: PetscBool flg;
157: PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
158: if (flg) *(void**)ctx = ((Mat_Shell*)(mat->data))->ctx;
159: else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot get context from non-shell matrix");
160: return(0);
161: }
163: PetscErrorCode MatDestroy_Shell(Mat mat)
164: {
166: Mat_Shell *shell = (Mat_Shell*)mat->data;
169: if (shell->ops->destroy) {
170: (*shell->ops->destroy)(mat);
171: }
172: VecDestroy(&shell->left);
173: VecDestroy(&shell->right);
174: VecDestroy(&shell->dshift);
175: VecDestroy(&shell->left_work);
176: VecDestroy(&shell->right_work);
177: VecDestroy(&shell->left_add_work);
178: VecDestroy(&shell->right_add_work);
179: MatDestroy(&shell->axpy);
180: PetscFree(mat->data);
181: return(0);
182: }
184: PetscErrorCode MatCopy_Shell(Mat A,Mat B,MatStructure str)
185: {
186: Mat_Shell *shellA = (Mat_Shell*)A->data,*shellB = (Mat_Shell*)B->data;
187: PetscErrorCode ierr;
188: PetscBool matflg;
191: PetscObjectTypeCompare((PetscObject)B,MATSHELL,&matflg);
192: if (!matflg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_NOTSAMETYPE,"Output matrix must be a MATSHELL");
194: PetscMemcpy(B->ops,A->ops,sizeof(struct _MatOps));
195: PetscMemcpy(shellB->ops,shellA->ops,sizeof(struct _MatShellOps));
197: if (shellA->ops->copy) {
198: (*shellA->ops->copy)(A,B,str);
199: }
200: shellB->vscale = shellA->vscale;
201: shellB->vshift = shellA->vshift;
202: if (shellA->dshift) {
203: if (!shellB->dshift) {
204: VecDuplicate(shellA->dshift,&shellB->dshift);
205: }
206: VecCopy(shellA->dshift,shellB->dshift);
207: } else {
208: VecDestroy(&shellB->dshift);
209: }
210: if (shellA->left) {
211: if (!shellB->left) {
212: VecDuplicate(shellA->left,&shellB->left);
213: }
214: VecCopy(shellA->left,shellB->left);
215: } else {
216: VecDestroy(&shellB->left);
217: }
218: if (shellA->right) {
219: if (!shellB->right) {
220: VecDuplicate(shellA->right,&shellB->right);
221: }
222: VecCopy(shellA->right,shellB->right);
223: } else {
224: VecDestroy(&shellB->right);
225: }
226: MatDestroy(&shellB->axpy);
227: if (shellA->axpy) {
228: PetscObjectReference((PetscObject)shellA->axpy);
229: shellB->axpy = shellA->axpy;
230: shellB->axpy_vscale = shellA->axpy_vscale;
231: }
232: return(0);
233: }
235: PetscErrorCode MatDuplicate_Shell(Mat mat,MatDuplicateOption op,Mat *M)
236: {
238: void *ctx;
241: MatShellGetContext(mat,&ctx);
242: MatCreateShell(PetscObjectComm((PetscObject)mat),mat->rmap->n,mat->cmap->n,mat->rmap->N,mat->cmap->N,ctx,M);
243: if (op != MAT_DO_NOT_COPY_VALUES) {
244: MatCopy(mat,*M,SAME_NONZERO_PATTERN);
245: }
246: return(0);
247: }
249: PetscErrorCode MatMult_Shell(Mat A,Vec x,Vec y)
250: {
251: Mat_Shell *shell = (Mat_Shell*)A->data;
252: PetscErrorCode ierr;
253: Vec xx;
254: PetscObjectState instate,outstate;
257: MatShellPreScaleRight(A,x,&xx);
258: PetscObjectStateGet((PetscObject)y, &instate);
259: if (!shell->ops->mult) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Have not provided a MatMult() for this MATSHELL");
260: (*shell->ops->mult)(A,xx,y);
261: PetscObjectStateGet((PetscObject)y, &outstate);
262: if (instate == outstate) {
263: /* increase the state of the output vector since the user did not update its state themself as should have been done */
264: PetscObjectStateIncrease((PetscObject)y);
265: }
266: MatShellShiftAndScale(A,xx,y);
267: MatShellPostScaleLeft(A,y);
269: if (shell->axpy) {
270: if (!shell->left_work) {MatCreateVecs(A,&shell->left_work,NULL);}
271: MatMult(shell->axpy,x,shell->left_work);
272: VecAXPY(y,shell->axpy_vscale,shell->left_work);
273: }
274: return(0);
275: }
277: PetscErrorCode MatMultAdd_Shell(Mat A,Vec x,Vec y,Vec z)
278: {
279: Mat_Shell *shell = (Mat_Shell*)A->data;
283: if (y == z) {
284: if (!shell->right_add_work) {VecDuplicate(z,&shell->right_add_work);}
285: MatMult(A,x,shell->right_add_work);
286: VecAXPY(z,1.0,shell->right_add_work);
287: } else {
288: MatMult(A,x,z);
289: VecAXPY(z,1.0,y);
290: }
291: return(0);
292: }
294: PetscErrorCode MatMultTranspose_Shell(Mat A,Vec x,Vec y)
295: {
296: Mat_Shell *shell = (Mat_Shell*)A->data;
297: PetscErrorCode ierr;
298: Vec xx;
299: PetscObjectState instate,outstate;
302: MatShellPreScaleLeft(A,x,&xx);
303: PetscObjectStateGet((PetscObject)y, &instate);
304: if (!shell->ops->multtranspose) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Have not provided a MatMultTranspose() for this MATSHELL");
305: (*shell->ops->multtranspose)(A,xx,y);
306: PetscObjectStateGet((PetscObject)y, &outstate);
307: if (instate == outstate) {
308: /* increase the state of the output vector since the user did not update its state themself as should have been done */
309: PetscObjectStateIncrease((PetscObject)y);
310: }
311: MatShellShiftAndScale(A,xx,y);
312: MatShellPostScaleRight(A,y);
314: if (shell->axpy) {
315: if (!shell->right_work) {MatCreateVecs(A,NULL,&shell->right_work);}
316: MatMultTranspose(shell->axpy,x,shell->right_work);
317: VecAXPY(y,shell->axpy_vscale,shell->right_work);
318: }
319: return(0);
320: }
322: PetscErrorCode MatMultTransposeAdd_Shell(Mat A,Vec x,Vec y,Vec z)
323: {
324: Mat_Shell *shell = (Mat_Shell*)A->data;
328: if (y == z) {
329: if (!shell->left_add_work) {VecDuplicate(z,&shell->left_add_work);}
330: MatMultTranspose(A,x,shell->left_add_work);
331: VecAXPY(z,1.0,shell->left_add_work);
332: } else {
333: MatMultTranspose(A,x,z);
334: VecAXPY(z,1.0,y);
335: }
336: return(0);
337: }
339: /*
340: diag(left)(vscale*A + diag(dshift) + vshift I)diag(right)
341: */
342: PetscErrorCode MatGetDiagonal_Shell(Mat A,Vec v)
343: {
344: Mat_Shell *shell = (Mat_Shell*)A->data;
348: if (shell->ops->getdiagonal) {
349: (*shell->ops->getdiagonal)(A,v);
350: } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Must provide shell matrix with routine to return diagonal using\nMatShellSetOperation(S,MATOP_GET_DIAGONAL,...)");
351: VecScale(v,shell->vscale);
352: if (shell->dshift) {
353: VecAXPY(v,1.0,shell->dshift);
354: }
355: VecShift(v,shell->vshift);
356: if (shell->left) {VecPointwiseMult(v,v,shell->left);}
357: if (shell->right) {VecPointwiseMult(v,v,shell->right);}
358: if (shell->axpy) {
359: if (!shell->left_work) {VecDuplicate(v,&shell->left_work);}
360: MatGetDiagonal(shell->axpy,shell->left_work);
361: VecAXPY(v,shell->axpy_vscale,shell->left_work);
362: }
363: return(0);
364: }
366: PetscErrorCode MatShift_Shell(Mat Y,PetscScalar a)
367: {
368: Mat_Shell *shell = (Mat_Shell*)Y->data;
372: if (shell->left || shell->right) {
373: if (!shell->dshift) {
374: VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift);
375: VecSet(shell->dshift,a);
376: } else {
377: if (shell->left) {VecPointwiseMult(shell->dshift,shell->dshift,shell->left);}
378: if (shell->right) {VecPointwiseMult(shell->dshift,shell->dshift,shell->right);}
379: VecShift(shell->dshift,a);
380: }
381: if (shell->left) {VecPointwiseDivide(shell->dshift,shell->dshift,shell->left);}
382: if (shell->right) {VecPointwiseDivide(shell->dshift,shell->dshift,shell->right);}
383: } else shell->vshift += a;
384: return(0);
385: }
387: PetscErrorCode MatDiagonalSet_Shell_Private(Mat A,Vec D,PetscScalar s)
388: {
389: Mat_Shell *shell = (Mat_Shell*)A->data;
393: if (!shell->dshift) {VecDuplicate(D,&shell->dshift);}
394: if (shell->left || shell->right) {
395: if (!shell->right_work) {VecDuplicate(shell->left ? shell->left : shell->right, &shell->right_work);}
396: if (shell->left && shell->right) {
397: VecPointwiseDivide(shell->right_work,D,shell->left);
398: VecPointwiseDivide(shell->right_work,shell->right_work,shell->right);
399: } else if (shell->left) {
400: VecPointwiseDivide(shell->right_work,D,shell->left);
401: } else {
402: VecPointwiseDivide(shell->right_work,D,shell->right);
403: }
404: if (!shell->dshift) {
405: VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift);
406: VecCopy(shell->dshift,shell->right_work);
407: } else {
408: VecAXPY(shell->dshift,s,shell->right_work);
409: }
410: } else {
411: VecAXPY(shell->dshift,s,D);
412: }
413: return(0);
414: }
416: PetscErrorCode MatDiagonalSet_Shell(Mat A,Vec D,InsertMode ins)
417: {
418: Vec d;
422: if (ins == INSERT_VALUES) {
423: if (!A->ops->getdiagonal) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Operation MATOP_GETDIAGONAL must be set first");
424: VecDuplicate(D,&d);
425: MatGetDiagonal(A,d);
426: MatDiagonalSet_Shell_Private(A,d,-1.);
427: MatDiagonalSet_Shell_Private(A,D,1.);
428: VecDestroy(&d);
429: } else {
430: MatDiagonalSet_Shell_Private(A,D,1.);
431: }
432: return(0);
433: }
435: PetscErrorCode MatScale_Shell(Mat Y,PetscScalar a)
436: {
437: Mat_Shell *shell = (Mat_Shell*)Y->data;
441: shell->vscale *= a;
442: shell->vshift *= a;
443: if (shell->dshift) {
444: VecScale(shell->dshift,a);
445: }
446: shell->axpy_vscale *= a;
447: return(0);
448: }
450: static PetscErrorCode MatDiagonalScale_Shell(Mat Y,Vec left,Vec right)
451: {
452: Mat_Shell *shell = (Mat_Shell*)Y->data;
456: if (shell->axpy) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_SUP,"Cannot diagonal scale MATSHELL after MatAXPY operation");
457: if (left) {
458: if (!shell->left) {
459: VecDuplicate(left,&shell->left);
460: VecCopy(left,shell->left);
461: } else {
462: VecPointwiseMult(shell->left,shell->left,left);
463: }
464: }
465: if (right) {
466: if (!shell->right) {
467: VecDuplicate(right,&shell->right);
468: VecCopy(right,shell->right);
469: } else {
470: VecPointwiseMult(shell->right,shell->right,right);
471: }
472: }
473: return(0);
474: }
476: PetscErrorCode MatAssemblyEnd_Shell(Mat Y,MatAssemblyType t)
477: {
478: Mat_Shell *shell = (Mat_Shell*)Y->data;
482: if (t == MAT_FINAL_ASSEMBLY) {
483: shell->vshift = 0.0;
484: shell->vscale = 1.0;
485: VecDestroy(&shell->dshift);
486: VecDestroy(&shell->left);
487: VecDestroy(&shell->right);
488: MatDestroy(&shell->axpy);
489: }
490: return(0);
491: }
493: static PetscErrorCode MatMissingDiagonal_Shell(Mat A,PetscBool *missing,PetscInt *d)
494: {
496: *missing = PETSC_FALSE;
497: return(0);
498: }
500: PetscErrorCode MatAXPY_Shell(Mat Y,PetscScalar a,Mat X,MatStructure str)
501: {
502: Mat_Shell *shell = (Mat_Shell*)Y->data;
506: PetscObjectReference((PetscObject)X);
507: MatDestroy(&shell->axpy);
508: shell->axpy = X;
509: shell->axpy_vscale = a;
510: return(0);
511: }
513: static struct _MatOps MatOps_Values = {0,
514: 0,
515: 0,
516: 0,
517: /* 4*/ MatMultAdd_Shell,
518: 0,
519: MatMultTransposeAdd_Shell,
520: 0,
521: 0,
522: 0,
523: /*10*/ 0,
524: 0,
525: 0,
526: 0,
527: 0,
528: /*15*/ 0,
529: 0,
530: 0,
531: MatDiagonalScale_Shell,
532: 0,
533: /*20*/ 0,
534: MatAssemblyEnd_Shell,
535: 0,
536: 0,
537: /*24*/ 0,
538: 0,
539: 0,
540: 0,
541: 0,
542: /*29*/ 0,
543: 0,
544: 0,
545: 0,
546: 0,
547: /*34*/ MatDuplicate_Shell,
548: 0,
549: 0,
550: 0,
551: 0,
552: /*39*/ MatAXPY_Shell,
553: 0,
554: 0,
555: 0,
556: MatCopy_Shell,
557: /*44*/ 0,
558: MatScale_Shell,
559: MatShift_Shell,
560: MatDiagonalSet_Shell,
561: 0,
562: /*49*/ 0,
563: 0,
564: 0,
565: 0,
566: 0,
567: /*54*/ 0,
568: 0,
569: 0,
570: 0,
571: 0,
572: /*59*/ 0,
573: MatDestroy_Shell,
574: 0,
575: 0,
576: 0,
577: /*64*/ 0,
578: 0,
579: 0,
580: 0,
581: 0,
582: /*69*/ 0,
583: 0,
584: MatConvert_Shell,
585: 0,
586: 0,
587: /*74*/ 0,
588: 0,
589: 0,
590: 0,
591: 0,
592: /*79*/ 0,
593: 0,
594: 0,
595: 0,
596: 0,
597: /*84*/ 0,
598: 0,
599: 0,
600: 0,
601: 0,
602: /*89*/ 0,
603: 0,
604: 0,
605: 0,
606: 0,
607: /*94*/ 0,
608: 0,
609: 0,
610: 0,
611: 0,
612: /*99*/ 0,
613: 0,
614: 0,
615: 0,
616: 0,
617: /*104*/ 0,
618: 0,
619: 0,
620: 0,
621: 0,
622: /*109*/ 0,
623: 0,
624: 0,
625: 0,
626: MatMissingDiagonal_Shell,
627: /*114*/ 0,
628: 0,
629: 0,
630: 0,
631: 0,
632: /*119*/ 0,
633: 0,
634: 0,
635: 0,
636: 0,
637: /*124*/ 0,
638: 0,
639: 0,
640: 0,
641: 0,
642: /*129*/ 0,
643: 0,
644: 0,
645: 0,
646: 0,
647: /*134*/ 0,
648: 0,
649: 0,
650: 0,
651: 0,
652: /*139*/ 0,
653: 0,
654: 0
655: };
657: /*MC
658: MATSHELL - MATSHELL = "shell" - A matrix type to be used to define your own matrix type -- perhaps matrix free.
660: Level: advanced
662: .seealso: MatCreateShell()
663: M*/
665: PETSC_EXTERN PetscErrorCode MatCreate_Shell(Mat A)
666: {
667: Mat_Shell *b;
671: PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));
673: PetscNewLog(A,&b);
674: A->data = (void*)b;
676: PetscLayoutSetUp(A->rmap);
677: PetscLayoutSetUp(A->cmap);
679: b->ctx = 0;
680: b->vshift = 0.0;
681: b->vscale = 1.0;
682: b->managescalingshifts = PETSC_TRUE;
683: A->assembled = PETSC_TRUE;
684: A->preallocated = PETSC_FALSE;
686: PetscObjectChangeTypeName((PetscObject)A,MATSHELL);
687: return(0);
688: }
690: /*@C
691: MatCreateShell - Creates a new matrix class for use with a user-defined
692: private data storage format.
694: Collective on MPI_Comm
696: Input Parameters:
697: + comm - MPI communicator
698: . m - number of local rows (must be given)
699: . n - number of local columns (must be given)
700: . M - number of global rows (may be PETSC_DETERMINE)
701: . N - number of global columns (may be PETSC_DETERMINE)
702: - ctx - pointer to data needed by the shell matrix routines
704: Output Parameter:
705: . A - the matrix
707: Level: advanced
709: Usage:
710: $ extern int mult(Mat,Vec,Vec);
711: $ MatCreateShell(comm,m,n,M,N,ctx,&mat);
712: $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
713: $ [ Use matrix for operations that have been set ]
714: $ MatDestroy(mat);
716: Notes:
717: The shell matrix type is intended to provide a simple class to use
718: with KSP (such as, for use with matrix-free methods). You should not
719: use the shell type if you plan to define a complete matrix class.
721: Fortran Notes:
722: To use this from Fortran with a ctx you must write an interface definition for this
723: function and for MatShellGetContext() that tells Fortran the Fortran derived data type you are passing
724: in as the ctx argument.
726: PETSc requires that matrices and vectors being used for certain
727: operations are partitioned accordingly. For example, when
728: creating a shell matrix, A, that supports parallel matrix-vector
729: products using MatMult(A,x,y) the user should set the number
730: of local matrix rows to be the number of local elements of the
731: corresponding result vector, y. Note that this is information is
732: required for use of the matrix interface routines, even though
733: the shell matrix may not actually be physically partitioned.
734: For example,
736: $
737: $ Vec x, y
738: $ extern int mult(Mat,Vec,Vec);
739: $ Mat A
740: $
741: $ VecCreateMPI(comm,PETSC_DECIDE,M,&y);
742: $ VecCreateMPI(comm,PETSC_DECIDE,N,&x);
743: $ VecGetLocalSize(y,&m);
744: $ VecGetLocalSize(x,&n);
745: $ MatCreateShell(comm,m,n,M,N,ctx,&A);
746: $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
747: $ MatMult(A,x,y);
748: $ MatDestroy(A);
749: $ VecDestroy(y); VecDestroy(x);
750: $
753: MATSHELL handles MatShift(), MatDiagonalSet(), MatDiagonalScale(), MatAXPY(), and MatScale() internally so these
754: operations cannot be overwritten unless MatShellSetManageScalingShifts() is called.
757: For rectangular matrices do all the scalings and shifts make sense?
759: Developers Notes:
760: Regarding shifting and scaling. The general form is
762: diag(left)(vscale*A + diag(dshift) + vshift I)diag(right)
764: The order you apply the operations is important. For example if you have a dshift then
765: apply a MatScale(s) you get s*vscale*A + s*diag(shift). But if you first scale and then shift
766: you get s*vscale*A + diag(shift)
768: A is the user provided function.
770: KSP/PC uses changes in the Mat's "state" to decide if preconditioners need to be rebuilt: PCSetUp() only calls the setup() for
771: for the PC implementation if the Mat state has increased from the previous call. Thus to get changes in a MATSHELL to trigger
772: an update in the preconditioner you must call MatAssemblyBegin()/MatAssemblyEnd() or PetscObjectStateIncrease((PetscObject)mat);
773: each time the MATSHELL matrix has changed.
775: Calling MatAssemblyBegin()/MatAssemblyEnd() on a MATSHELL removes any previously supplied shift and scales that were provided
776: with MatDiagonalSet(), MatShift(), MatScale(), or MatDiagonalScale().
778: .keywords: matrix, shell, create
780: .seealso: MatShellSetOperation(), MatHasOperation(), MatShellGetContext(), MatShellSetContext(), MATSHELL, MatShellSetManageScalingShifts()
781: @*/
782: PetscErrorCode MatCreateShell(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,void *ctx,Mat *A)
783: {
787: MatCreate(comm,A);
788: MatSetSizes(*A,m,n,M,N);
789: MatSetType(*A,MATSHELL);
790: MatShellSetContext(*A,ctx);
791: MatSetUp(*A);
792: return(0);
793: }
795: /*@
796: MatShellSetContext - sets the context for a shell matrix
798: Logically Collective on Mat
800: Input Parameters:
801: + mat - the shell matrix
802: - ctx - the context
804: Level: advanced
806: Fortran Notes:
807: To use this from Fortran you must write a Fortran interface definition for this
808: function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.
810: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation()
811: @*/
812: PetscErrorCode MatShellSetContext(Mat mat,void *ctx)
813: {
814: Mat_Shell *shell = (Mat_Shell*)mat->data;
816: PetscBool flg;
820: PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
821: if (flg) {
822: shell->ctx = ctx;
823: } else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot attach context to non-shell matrix");
824: return(0);
825: }
827: /*@
828: MatShellSetManageScalingShifts - Allows the user to control the scaling and shift operations of the MATSHELL. Must be called immediately
829: after MatCreateShell()
831: Logically Collective on Mat
833: Input Parameter:
834: . mat - the shell matrix
836: Level: advanced
838: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatShellSetOperation()
839: @*/
840: PetscErrorCode MatShellSetManageScalingShifts(Mat A)
841: {
843: Mat_Shell *shell;
844: PetscBool flg;
848: PetscObjectTypeCompare((PetscObject)A,MATSHELL,&flg);
849: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Can only use with MATSHELL matrices");
850: shell = (Mat_Shell*)A->data;
851: shell->managescalingshifts = PETSC_FALSE;
852: A->ops->diagonalset = NULL;
853: A->ops->diagonalscale = NULL;
854: A->ops->scale = NULL;
855: A->ops->shift = NULL;
856: A->ops->axpy = NULL;
857: return(0);
858: }
860: /*@C
861: MatShellTestMult - Compares the multiply routine provided to the MATSHELL with differencing on a given function.
863: Logically Collective on Mat
865: Input Parameters:
866: + mat - the shell matrix
867: . f - the function
868: . base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated
869: - ctx - an optional context for the function
871: Output Parameter:
872: . flg - PETSC_TRUE if the multiply is likely correct
874: Options Database:
875: . -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference
877: Level: advanced
879: Fortran Notes:
880: Not supported from Fortran
882: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellTestMultTranspose()
883: @*/
884: PetscErrorCode MatShellTestMult(Mat mat,PetscErrorCode (*f)(void*,Vec,Vec),Vec base,void *ctx,PetscBool *flg)
885: {
887: PetscInt m,n;
888: Mat mf,Dmf,Dmat,Ddiff;
889: PetscReal Diffnorm,Dmfnorm;
890: PetscBool v = PETSC_FALSE, flag = PETSC_TRUE;
894: PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_view",&v);
895: MatGetLocalSize(mat,&m,&n);
896: MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,PETSC_DECIDE,PETSC_DECIDE,&mf);
897: MatMFFDSetFunction(mf,f,ctx);
898: MatMFFDSetBase(mf,base,NULL);
900: MatComputeExplicitOperator(mf,&Dmf);
901: MatComputeExplicitOperator(mat,&Dmat);
903: MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff);
904: MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN);
905: MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm);
906: MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm);
907: if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) {
908: flag = PETSC_FALSE;
909: if (v) {
910: 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));
911: MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_view");
912: MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_view");
913: MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_view");
914: }
915: } else if (v) {
916: PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL and matrix free multiple appear to produce the same results\n");
917: }
918: if (flg) *flg = flag;
919: MatDestroy(&Ddiff);
920: MatDestroy(&mf);
921: MatDestroy(&Dmf);
922: MatDestroy(&Dmat);
923: return(0);
924: }
926: /*@C
927: MatShellTestMultTranpose - Compares the multiply transpose routine provided to the MATSHELL with differencing on a given function.
929: Logically Collective on Mat
931: Input Parameters:
932: + mat - the shell matrix
933: . f - the function
934: . base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated
935: - ctx - an optional context for the function
937: Output Parameter:
938: . flg - PETSC_TRUE if the multiply is likely correct
940: Options Database:
941: . -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference
943: Level: advanced
945: Fortran Notes:
946: Not supported from Fortran
948: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellTestMult()
949: @*/
950: PetscErrorCode MatShellTestMultTranspose(Mat mat,PetscErrorCode (*f)(void*,Vec,Vec),Vec base,void *ctx,PetscBool *flg)
951: {
953: Vec x,y,z;
954: PetscInt m,n,M,N;
955: Mat mf,Dmf,Dmat,Ddiff;
956: PetscReal Diffnorm,Dmfnorm;
957: PetscBool v = PETSC_FALSE, flag = PETSC_TRUE;
961: PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_transpose_view",&v);
962: MatCreateVecs(mat,&x,&y);
963: VecDuplicate(y,&z);
964: MatGetLocalSize(mat,&m,&n);
965: MatGetSize(mat,&M,&N);
966: MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,M,N,&mf);
967: MatMFFDSetFunction(mf,f,ctx);
968: MatMFFDSetBase(mf,base,NULL);
969: MatComputeExplicitOperator(mf,&Dmf);
970: MatTranspose(Dmf,MAT_INPLACE_MATRIX,&Dmf);
971: MatComputeExplicitOperatorTranspose(mat,&Dmat);
973: MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff);
974: MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN);
975: MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm);
976: MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm);
977: if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) {
978: flag = PETSC_FALSE;
979: if (v) {
980: 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));
981: MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");
982: MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");
983: MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");
984: }
985: } else if (v) {
986: PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL transpose and matrix free multiple appear to produce the same results\n");
987: }
988: if (flg) *flg = flag;
989: MatDestroy(&mf);
990: MatDestroy(&Dmat);
991: MatDestroy(&Ddiff);
992: MatDestroy(&Dmf);
993: VecDestroy(&x);
994: VecDestroy(&y);
995: VecDestroy(&z);
996: return(0);
997: }
999: /*@C
1000: MatShellSetOperation - Allows user to set a matrix operation for a shell matrix.
1002: Logically Collective on Mat
1004: Input Parameters:
1005: + mat - the shell matrix
1006: . op - the name of the operation
1007: - f - the function that provides the operation.
1009: Level: advanced
1011: Usage:
1012: $ extern PetscErrorCode usermult(Mat,Vec,Vec);
1013: $ MatCreateShell(comm,m,n,M,N,ctx,&A);
1014: $ MatShellSetOperation(A,MATOP_MULT,(void(*)(void))usermult);
1016: Notes:
1017: See the file include/petscmat.h for a complete list of matrix
1018: operations, which all have the form MATOP_<OPERATION>, where
1019: <OPERATION> is the name (in all capital letters) of the
1020: user interface routine (e.g., MatMult() -> MATOP_MULT).
1022: All user-provided functions (except for MATOP_DESTROY) should have the same calling
1023: sequence as the usual matrix interface routines, since they
1024: are intended to be accessed via the usual matrix interface
1025: routines, e.g.,
1026: $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec)
1028: In particular each function MUST return an error code of 0 on success and
1029: nonzero on failure.
1031: Within each user-defined routine, the user should call
1032: MatShellGetContext() to obtain the user-defined context that was
1033: set by MatCreateShell().
1035: Fortran Notes:
1036: For MatCreateVecs() the user code should check if the input left or right matrix is -1 and in that case not
1037: generate a matrix. See src/mat/examples/tests/ex120f.F
1039: Use MatSetOperation() to set an operation for any matrix type
1041: .keywords: matrix, shell, set, operation
1043: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatSetOperation(), MatShellSetManageScalingShifts()
1044: @*/
1045: PetscErrorCode MatShellSetOperation(Mat mat,MatOperation op,void (*f)(void))
1046: {
1047: PetscBool flg;
1048: Mat_Shell *shell;
1053: PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
1054: if (!flg) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Can only use with MATSHELL matrices");
1055: shell = (Mat_Shell*)mat->data;
1057: switch (op) {
1058: case MATOP_DESTROY:
1059: shell->ops->destroy = (PetscErrorCode (*)(Mat))f;
1060: break;
1061: case MATOP_VIEW:
1062: if (!mat->ops->viewnative) {
1063: mat->ops->viewnative = mat->ops->view;
1064: }
1065: mat->ops->view = (PetscErrorCode (*)(Mat,PetscViewer))f;
1066: break;
1067: case MATOP_COPY:
1068: shell->ops->copy = (PetscErrorCode (*)(Mat,Mat,MatStructure))f;
1069: break;
1070: case MATOP_DIAGONAL_SET:
1071: case MATOP_DIAGONAL_SCALE:
1072: case MATOP_SHIFT:
1073: case MATOP_SCALE:
1074: case MATOP_AXPY:
1075: if (shell->managescalingshifts) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"MATSHELL is managing scalings and shifts, see MatShellSetManageScalingShifts()");
1076: (((void(**)(void))mat->ops)[op]) = f;
1077: break;
1078: case MATOP_GET_DIAGONAL:
1079: if (shell->managescalingshifts) {
1080: shell->ops->getdiagonal = (PetscErrorCode (*)(Mat,Vec))f;
1081: mat->ops->getdiagonal = MatGetDiagonal_Shell;
1082: } else {
1083: shell->ops->getdiagonal = NULL;
1084: mat->ops->getdiagonal = (PetscErrorCode (*)(Mat,Vec))f;
1085: }
1086: break;
1087: case MATOP_MULT:
1088: if (shell->managescalingshifts) {
1089: shell->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1090: mat->ops->mult = MatMult_Shell;
1091: } else {
1092: shell->ops->mult = NULL;
1093: mat->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1094: }
1095: break;
1096: case MATOP_MULT_TRANSPOSE:
1097: if (shell->managescalingshifts) {
1098: shell->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1099: mat->ops->multtranspose = MatMultTranspose_Shell;
1100: } else {
1101: shell->ops->multtranspose = NULL;
1102: mat->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1103: }
1104: break;
1105: default:
1106: (((void(**)(void))mat->ops)[op]) = f;
1107: break;
1108: }
1109: return(0);
1110: }
1112: /*@C
1113: MatShellGetOperation - Gets a matrix function for a shell matrix.
1115: Not Collective
1117: Input Parameters:
1118: + mat - the shell matrix
1119: - op - the name of the operation
1121: Output Parameter:
1122: . f - the function that provides the operation.
1124: Level: advanced
1126: Notes:
1127: See the file include/petscmat.h for a complete list of matrix
1128: operations, which all have the form MATOP_<OPERATION>, where
1129: <OPERATION> is the name (in all capital letters) of the
1130: user interface routine (e.g., MatMult() -> MATOP_MULT).
1132: All user-provided functions have the same calling
1133: sequence as the usual matrix interface routines, since they
1134: are intended to be accessed via the usual matrix interface
1135: routines, e.g.,
1136: $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec)
1138: Within each user-defined routine, the user should call
1139: MatShellGetContext() to obtain the user-defined context that was
1140: set by MatCreateShell().
1142: .keywords: matrix, shell, set, operation
1144: .seealso: MatCreateShell(), MatShellGetContext(), MatShellSetOperation(), MatShellSetContext()
1145: @*/
1146: PetscErrorCode MatShellGetOperation(Mat mat,MatOperation op,void(**f)(void))
1147: {
1148: PetscBool flg;
1149: Mat_Shell *shell;
1154: PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
1155: if (!flg) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Can only use with MATSHELL matrices");
1156: shell = (Mat_Shell*)mat->data;
1158: switch (op) {
1159: case MATOP_DESTROY:
1160: *f = (void (*)(void))shell->ops->destroy;
1161: break;
1162: case MATOP_VIEW:
1163: *f = (void (*)(void))mat->ops->view;
1164: break;
1165: case MATOP_COPY:
1166: *f = (void (*)(void))shell->ops->copy;
1167: break;
1168: case MATOP_DIAGONAL_SET:
1169: case MATOP_DIAGONAL_SCALE:
1170: case MATOP_SHIFT:
1171: case MATOP_SCALE:
1172: case MATOP_AXPY:
1173: *f = (((void (**)(void))mat->ops)[op]);
1174: break;
1175: case MATOP_GET_DIAGONAL:
1176: if (shell->ops->getdiagonal)
1177: *f = (void (*)(void))shell->ops->getdiagonal;
1178: else
1179: *f = (((void (**)(void))mat->ops)[op]);
1180: break;
1181: case MATOP_MULT:
1182: if (shell->ops->mult)
1183: *f = (void (*)(void))shell->ops->mult;
1184: else
1185: *f = (((void (**)(void))mat->ops)[op]);
1186: break;
1187: case MATOP_MULT_TRANSPOSE:
1188: if (shell->ops->multtranspose)
1189: *f = (void (*)(void))shell->ops->multtranspose;
1190: else
1191: *f = (((void (**)(void))mat->ops)[op]);
1192: break;
1193: default:
1194: *f = (((void (**)(void))mat->ops)[op]);
1195: }
1196: return(0);
1197: }