Actual source code: shell.c
petsc-3.3-p7 2013-05-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> /*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);
16: PetscScalar vscale,vshift;
17: Vec dshift;
18: Vec left,right;
19: Vec dshift_owned,left_owned,right_owned;
20: Vec left_work,right_work;
21: PetscBool usingscaled;
22: void *ctx;
23: } Mat_Shell;
24: /*
25: The most general expression for the matrix is
27: S = L (a A + B) R
29: where
30: A is the matrix defined by the user's function
31: a is a scalar multiple
32: L is left scaling
33: R is right scaling
34: B is a diagonal shift defined by
35: diag(dshift) if the vector dshift is non-NULL
36: vscale*identity otherwise
38: The following identities apply:
40: Scale by c:
41: c [L (a A + B) R] = L [(a c) A + c B] R
43: Shift by c:
44: [L (a A + B) R] + c = L [a A + (B + c Linv Rinv)] R
46: Diagonal scaling is achieved by simply multiplying with existing L and R vectors
48: In the data structure:
50: vscale=1.0 means no special scaling will be applied
51: dshift=NULL means a constant diagonal shift (fall back to vshift)
52: vshift=0.0 means no constant diagonal shift, note that vshift is only used if dshift is NULL
53: */
55: static PetscErrorCode MatMult_Shell(Mat,Vec,Vec);
56: static PetscErrorCode MatMultTranspose_Shell(Mat,Vec,Vec);
57: static PetscErrorCode MatGetDiagonal_Shell(Mat,Vec);
61: static PetscErrorCode MatShellUseScaledMethods(Mat Y)
62: {
63: Mat_Shell *shell = (Mat_Shell*)Y->data;
66: if (shell->usingscaled) return(0);
67: shell->mult = Y->ops->mult;
68: Y->ops->mult = MatMult_Shell;
69: if (Y->ops->multtranspose) {
70: shell->multtranspose = Y->ops->multtranspose;
71: Y->ops->multtranspose = MatMultTranspose_Shell;
72: }
73: if (Y->ops->getdiagonal) {
74: shell->getdiagonal = Y->ops->getdiagonal;
75: Y->ops->getdiagonal = MatGetDiagonal_Shell;
76: }
77: shell->usingscaled = PETSC_TRUE;
78: return(0);
79: }
83: static PetscErrorCode MatShellPreScaleLeft(Mat A,Vec x,Vec *xx)
84: {
85: Mat_Shell *shell = (Mat_Shell*)A->data;
89: *xx = PETSC_NULL;
90: if (!shell->left) {
91: *xx = x;
92: } else {
93: if (!shell->left_work) {VecDuplicate(shell->left,&shell->left_work);}
94: VecPointwiseMult(shell->left_work,x,shell->left);
95: *xx = shell->left_work;
96: }
97: return(0);
98: }
102: static PetscErrorCode MatShellPreScaleRight(Mat A,Vec x,Vec *xx)
103: {
104: Mat_Shell *shell = (Mat_Shell*)A->data;
108: *xx = PETSC_NULL;
109: if (!shell->right) {
110: *xx = x;
111: } else {
112: if (!shell->right_work) {VecDuplicate(shell->right,&shell->right_work);}
113: VecPointwiseMult(shell->right_work,x,shell->right);
114: *xx = shell->right_work;
115: }
116: return(0);
117: }
121: static PetscErrorCode MatShellPostScaleLeft(Mat A,Vec x)
122: {
123: Mat_Shell *shell = (Mat_Shell*)A->data;
127: if (shell->left) {VecPointwiseMult(x,x,shell->left);}
128: return(0);
129: }
133: static PetscErrorCode MatShellPostScaleRight(Mat A,Vec x)
134: {
135: Mat_Shell *shell = (Mat_Shell*)A->data;
139: if (shell->right) {VecPointwiseMult(x,x,shell->right);}
140: return(0);
141: }
145: static PetscErrorCode MatShellShiftAndScale(Mat A,Vec X,Vec Y)
146: {
147: Mat_Shell *shell = (Mat_Shell*)A->data;
151: if (shell->dshift) { /* get arrays because there is no VecPointwiseMultAdd() */
152: PetscInt i,m;
153: const PetscScalar *x,*d;
154: PetscScalar *y;
155: VecGetLocalSize(X,&m);
156: VecGetArrayRead(shell->dshift,&d);
157: VecGetArrayRead(X,&x);
158: VecGetArray(Y,&y);
159: for (i=0; i<m; i++) y[i] = shell->vscale*y[i] + d[i]*x[i];
160: VecRestoreArrayRead(shell->dshift,&d);
161: VecRestoreArrayRead(X,&x);
162: VecRestoreArray(Y,&y);
163: } else if (PetscAbsScalar(shell->vshift) != 0) {
164: VecAXPBY(Y,shell->vshift,shell->vscale,X);
165: } else if (shell->vscale != 1.0) {
166: VecScale(Y,shell->vscale);
167: }
168: return(0);
169: }
173: /*@
174: MatShellGetContext - Returns the user-provided context associated with a shell matrix.
176: Not Collective
178: Input Parameter:
179: . mat - the matrix, should have been created with MatCreateShell()
181: Output Parameter:
182: . ctx - the user provided context
184: Level: advanced
186: Notes:
187: This routine is intended for use within various shell matrix routines,
188: as set with MatShellSetOperation().
189:
190: .keywords: matrix, shell, get, context
192: .seealso: MatCreateShell(), MatShellSetOperation(), MatShellSetContext()
193: @*/
194: PetscErrorCode MatShellGetContext(Mat mat,void *ctx)
195: {
197: PetscBool flg;
202: PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
203: if (flg) *(void**)ctx = ((Mat_Shell*)(mat->data))->ctx;
204: else SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Cannot get context from non-shell matrix");
205: return(0);
206: }
210: PetscErrorCode MatDestroy_Shell(Mat mat)
211: {
213: Mat_Shell *shell = (Mat_Shell*)mat->data;
216: if (shell->destroy) {
217: (*shell->destroy)(mat);
218: }
219: VecDestroy(&shell->left_owned);
220: VecDestroy(&shell->right_owned);
221: VecDestroy(&shell->dshift_owned);
222: VecDestroy(&shell->left_work);
223: VecDestroy(&shell->right_work);
224: PetscFree(mat->data);
225: return(0);
226: }
230: PetscErrorCode MatMult_Shell(Mat A,Vec x,Vec y)
231: {
232: Mat_Shell *shell = (Mat_Shell*)A->data;
234: Vec xx;
237: MatShellPreScaleRight(A,x,&xx);
238: (*shell->mult)(A,xx,y);
239: MatShellShiftAndScale(A,xx,y);
240: MatShellPostScaleLeft(A,y);
241: return(0);
242: }
246: PetscErrorCode MatMultTranspose_Shell(Mat A,Vec x,Vec y)
247: {
248: Mat_Shell *shell = (Mat_Shell*)A->data;
250: Vec xx;
253: MatShellPreScaleLeft(A,x,&xx);
254: (*shell->multtranspose)(A,xx,y);
255: MatShellShiftAndScale(A,xx,y);
256: MatShellPostScaleRight(A,y);
257: return(0);
258: }
262: PetscErrorCode MatGetDiagonal_Shell(Mat A,Vec v)
263: {
264: Mat_Shell *shell = (Mat_Shell*)A->data;
268: (*shell->getdiagonal)(A,v);
269: VecScale(v,shell->vscale);
270: if (shell->dshift) {
271: VecPointwiseMult(v,v,shell->dshift);
272: } else {
273: VecShift(v,shell->vshift);
274: }
275: if (shell->left) {VecPointwiseMult(v,v,shell->left);}
276: if (shell->right) {VecPointwiseMult(v,v,shell->right);}
277: return(0);
278: }
282: PetscErrorCode MatShift_Shell(Mat Y,PetscScalar a)
283: {
284: Mat_Shell *shell = (Mat_Shell*)Y->data;
288: if (shell->left || shell->right || shell->dshift) {
289: if (!shell->dshift) {
290: if (!shell->dshift_owned) {VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift_owned);}
291: shell->dshift = shell->dshift_owned;
292: VecSet(shell->dshift,shell->vshift+a);
293: } else {VecScale(shell->dshift,a);}
294: if (shell->left) {VecPointwiseDivide(shell->dshift,shell->dshift,shell->left);}
295: if (shell->right) {VecPointwiseDivide(shell->dshift,shell->dshift,shell->right);}
296: } else shell->vshift += a;
297: MatShellUseScaledMethods(Y);
298: return(0);
299: }
303: PetscErrorCode MatScale_Shell(Mat Y,PetscScalar a)
304: {
305: Mat_Shell *shell = (Mat_Shell*)Y->data;
309: shell->vscale *= a;
310: if (shell->dshift) {VecScale(shell->dshift,a);}
311: else shell->vshift *= a;
312: MatShellUseScaledMethods(Y);
313: return(0);
314: }
318: static PetscErrorCode MatDiagonalScale_Shell(Mat Y,Vec left,Vec right)
319: {
320: Mat_Shell *shell = (Mat_Shell*)Y->data;
324: if (left) {
325: if (!shell->left_owned) {VecDuplicate(left,&shell->left_owned);}
326: if (shell->left) {VecPointwiseMult(shell->left,shell->left,left);}
327: else {
328: shell->left = shell->left_owned;
329: VecCopy(left,shell->left);
330: }
331: }
332: if (right) {
333: if (!shell->right_owned) {VecDuplicate(right,&shell->right_owned);}
334: if (shell->right) {VecPointwiseMult(shell->right,shell->right,right);}
335: else {
336: shell->right = shell->right_owned;
337: VecCopy(right,shell->right);
338: }
339: }
340: MatShellUseScaledMethods(Y);
341: return(0);
342: }
346: PetscErrorCode MatAssemblyEnd_Shell(Mat Y,MatAssemblyType t)
347: {
348: Mat_Shell *shell = (Mat_Shell*)Y->data;
351: if (t == MAT_FINAL_ASSEMBLY) {
352: shell->vshift = 0.0;
353: shell->vscale = 1.0;
354: shell->dshift = PETSC_NULL;
355: shell->left = PETSC_NULL;
356: shell->right = PETSC_NULL;
357: if (shell->mult) {
358: Y->ops->mult = shell->mult;
359: shell->mult = PETSC_NULL;
360: }
361: if (shell->multtranspose) {
362: Y->ops->multtranspose = shell->multtranspose;
363: shell->multtranspose = PETSC_NULL;
364: }
365: if (shell->getdiagonal) {
366: Y->ops->getdiagonal = shell->getdiagonal;
367: shell->getdiagonal = PETSC_NULL;
368: }
369: shell->usingscaled = PETSC_FALSE;
370: }
371: return(0);
372: }
374: extern PetscErrorCode MatConvert_Shell(Mat, const MatType,MatReuse,Mat*);
376: static struct _MatOps MatOps_Values = {0,
377: 0,
378: 0,
379: 0,
380: /* 4*/ 0,
381: 0,
382: 0,
383: 0,
384: 0,
385: 0,
386: /*10*/ 0,
387: 0,
388: 0,
389: 0,
390: 0,
391: /*15*/ 0,
392: 0,
393: 0,
394: MatDiagonalScale_Shell,
395: 0,
396: /*20*/ 0,
397: MatAssemblyEnd_Shell,
398: 0,
399: 0,
400: /*24*/ 0,
401: 0,
402: 0,
403: 0,
404: 0,
405: /*29*/ 0,
406: 0,
407: 0,
408: 0,
409: 0,
410: /*34*/ 0,
411: 0,
412: 0,
413: 0,
414: 0,
415: /*39*/ 0,
416: 0,
417: 0,
418: 0,
419: 0,
420: /*44*/ 0,
421: MatScale_Shell,
422: MatShift_Shell,
423: 0,
424: 0,
425: /*49*/ 0,
426: 0,
427: 0,
428: 0,
429: 0,
430: /*54*/ 0,
431: 0,
432: 0,
433: 0,
434: 0,
435: /*59*/ 0,
436: MatDestroy_Shell,
437: 0,
438: 0,
439: 0,
440: /*64*/ 0,
441: 0,
442: 0,
443: 0,
444: 0,
445: /*69*/ 0,
446: 0,
447: MatConvert_Shell,
448: 0,
449: 0,
450: /*74*/ 0,
451: 0,
452: 0,
453: 0,
454: 0,
455: /*79*/ 0,
456: 0,
457: 0,
458: 0,
459: 0,
460: /*84*/ 0,
461: 0,
462: 0,
463: 0,
464: 0,
465: /*89*/ 0,
466: 0,
467: 0,
468: 0,
469: 0,
470: /*94*/ 0,
471: 0,
472: 0,
473: 0};
475: /*MC
476: MATSHELL - MATSHELL = "shell" - A matrix type to be used to define your own matrix type -- perhaps matrix free.
478: Level: advanced
480: .seealso: MatCreateShell
481: M*/
483: EXTERN_C_BEGIN
486: PetscErrorCode MatCreate_Shell(Mat A)
487: {
488: Mat_Shell *b;
492: PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));
494: PetscNewLog(A,Mat_Shell,&b);
495: A->data = (void*)b;
497: PetscLayoutSetUp(A->rmap);
498: PetscLayoutSetUp(A->cmap);
500: b->ctx = 0;
501: b->vshift = 0.0;
502: b->vscale = 1.0;
503: b->mult = 0;
504: b->multtranspose = 0;
505: b->getdiagonal = 0;
506: A->assembled = PETSC_TRUE;
507: A->preallocated = PETSC_FALSE;
508: PetscObjectChangeTypeName((PetscObject)A,MATSHELL);
509: return(0);
510: }
511: EXTERN_C_END
515: /*@C
516: MatCreateShell - Creates a new matrix class for use with a user-defined
517: private data storage format.
519: Collective on MPI_Comm
521: Input Parameters:
522: + comm - MPI communicator
523: . m - number of local rows (must be given)
524: . n - number of local columns (must be given)
525: . M - number of global rows (may be PETSC_DETERMINE)
526: . N - number of global columns (may be PETSC_DETERMINE)
527: - ctx - pointer to data needed by the shell matrix routines
529: Output Parameter:
530: . A - the matrix
532: Level: advanced
534: Usage:
535: $ extern int mult(Mat,Vec,Vec);
536: $ MatCreateShell(comm,m,n,M,N,ctx,&mat);
537: $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
538: $ [ Use matrix for operations that have been set ]
539: $ MatDestroy(mat);
541: Notes:
542: The shell matrix type is intended to provide a simple class to use
543: with KSP (such as, for use with matrix-free methods). You should not
544: use the shell type if you plan to define a complete matrix class.
546: Fortran Notes: The context can only be an integer or a PetscObject
547: unfortunately it cannot be a Fortran array or derived type.
549: PETSc requires that matrices and vectors being used for certain
550: operations are partitioned accordingly. For example, when
551: creating a shell matrix, A, that supports parallel matrix-vector
552: products using MatMult(A,x,y) the user should set the number
553: of local matrix rows to be the number of local elements of the
554: corresponding result vector, y. Note that this is information is
555: required for use of the matrix interface routines, even though
556: the shell matrix may not actually be physically partitioned.
557: For example,
559: $
560: $ Vec x, y
561: $ extern int mult(Mat,Vec,Vec);
562: $ Mat A
563: $
564: $ VecCreateMPI(comm,PETSC_DECIDE,M,&y);
565: $ VecCreateMPI(comm,PETSC_DECIDE,N,&x);
566: $ VecGetLocalSize(y,&m);
567: $ VecGetLocalSize(x,&n);
568: $ MatCreateShell(comm,m,n,M,N,ctx,&A);
569: $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
570: $ MatMult(A,x,y);
571: $ MatDestroy(A);
572: $ VecDestroy(y); VecDestroy(x);
573: $
575: .keywords: matrix, shell, create
577: .seealso: MatShellSetOperation(), MatHasOperation(), MatShellGetContext(), MatShellSetContext()
578: @*/
579: PetscErrorCode MatCreateShell(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,void *ctx,Mat *A)
580: {
584: MatCreate(comm,A);
585: MatSetSizes(*A,m,n,M,N);
586: MatSetType(*A,MATSHELL);
587: MatShellSetContext(*A,ctx);
588: MatSetUp(*A);
589: return(0);
590: }
594: /*@
595: MatShellSetContext - sets the context for a shell matrix
597: Logically Collective on Mat
599: Input Parameters:
600: + mat - the shell matrix
601: - ctx - the context
603: Level: advanced
605: Fortran Notes: The context can only be an integer or a PetscObject
606: unfortunately it cannot be a Fortran array or derived type.
608: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation()
609: @*/
610: PetscErrorCode MatShellSetContext(Mat mat,void *ctx)
611: {
612: Mat_Shell *shell = (Mat_Shell*)mat->data;
614: PetscBool flg;
618: PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
619: if (flg) {
620: shell->ctx = ctx;
621: } else SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Cannot attach context to non-shell matrix");
622: return(0);
623: }
627: /*@C
628: MatShellSetOperation - Allows user to set a matrix operation for
629: a shell matrix.
631: Logically Collective on Mat
633: Input Parameters:
634: + mat - the shell matrix
635: . op - the name of the operation
636: - f - the function that provides the operation.
638: Level: advanced
640: Usage:
641: $ extern PetscErrorCode usermult(Mat,Vec,Vec);
642: $ MatCreateShell(comm,m,n,M,N,ctx,&A);
643: $ MatShellSetOperation(A,MATOP_MULT,(void(*)(void))usermult);
645: Notes:
646: See the file include/petscmat.h for a complete list of matrix
647: operations, which all have the form MATOP_<OPERATION>, where
648: <OPERATION> is the name (in all capital letters) of the
649: user interface routine (e.g., MatMult() -> MATOP_MULT).
651: All user-provided functions (execept for MATOP_DESTROY) should have the same calling
652: sequence as the usual matrix interface routines, since they
653: are intended to be accessed via the usual matrix interface
654: routines, e.g.,
655: $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec)
657: In particular each function MUST return an error code of 0 on success and
658: nonzero on failure.
660: Within each user-defined routine, the user should call
661: MatShellGetContext() to obtain the user-defined context that was
662: set by MatCreateShell().
664: Fortran Notes: For MatGetVecs() the user code should check if the input left or right matrix is -1 and in that case not
665: generate a matrix. See src/mat/examples/tests/ex120f.F
667: .keywords: matrix, shell, set, operation
669: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext()
670: @*/
671: PetscErrorCode MatShellSetOperation(Mat mat,MatOperation op,void (*f)(void))
672: {
674: PetscBool flg;
678: if (op == MATOP_DESTROY) {
679: PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
680: if (flg) {
681: Mat_Shell *shell = (Mat_Shell*)mat->data;
682: shell->destroy = (PetscErrorCode (*)(Mat)) f;
683: } else mat->ops->destroy = (PetscErrorCode (*)(Mat)) f;
684: }
685: else if (op == MATOP_VIEW) mat->ops->view = (PetscErrorCode (*)(Mat,PetscViewer)) f;
686: else (((void(**)(void))mat->ops)[op]) = f;
688: return(0);
689: }
693: /*@C
694: MatShellGetOperation - Gets a matrix function for a shell matrix.
696: Not Collective
698: Input Parameters:
699: + mat - the shell matrix
700: - op - the name of the operation
702: Output Parameter:
703: . f - the function that provides the operation.
705: Level: advanced
707: Notes:
708: See the file include/petscmat.h for a complete list of matrix
709: operations, which all have the form MATOP_<OPERATION>, where
710: <OPERATION> is the name (in all capital letters) of the
711: user interface routine (e.g., MatMult() -> MATOP_MULT).
713: All user-provided functions have the same calling
714: sequence as the usual matrix interface routines, since they
715: are intended to be accessed via the usual matrix interface
716: routines, e.g.,
717: $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec)
719: Within each user-defined routine, the user should call
720: MatShellGetContext() to obtain the user-defined context that was
721: set by MatCreateShell().
723: .keywords: matrix, shell, set, operation
725: .seealso: MatCreateShell(), MatShellGetContext(), MatShellSetOperation(), MatShellSetContext()
726: @*/
727: PetscErrorCode MatShellGetOperation(Mat mat,MatOperation op,void(**f)(void))
728: {
730: PetscBool flg;
734: if (op == MATOP_DESTROY) {
735: PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
736: if (flg) {
737: Mat_Shell *shell = (Mat_Shell*)mat->data;
738: *f = (void(*)(void))shell->destroy;
739: } else {
740: *f = (void(*)(void))mat->ops->destroy;
741: }
742: } else if (op == MATOP_VIEW) {
743: *f = (void(*)(void))mat->ops->view;
744: } else {
745: *f = (((void(**)(void))mat->ops)[op]);
746: }
748: return(0);
749: }