Actual source code: shell.c
1: #define PETSCMAT_DLL
3: /*
4: This provides a simple shell for Fortran (and C programmers) to
5: create a very simple matrix class for use with KSP without coding
6: much of anything.
7: */
9: #include include/private/matimpl.h
10: #include private/vecimpl.h
12: typedef struct {
13: PetscErrorCode (*destroy)(Mat);
14: PetscErrorCode (*mult)(Mat,Vec,Vec);
15: PetscErrorCode (*multtranspose)(Mat,Vec,Vec);
16: PetscErrorCode (*getdiagonal)(Mat,Vec);
17: PetscTruth scale,shift;
18: PetscScalar vscale,vshift;
19: void *ctx;
20: } Mat_Shell;
24: /*@
25: MatShellGetContext - Returns the user-provided context associated with a shell matrix.
27: Not Collective
29: Input Parameter:
30: . mat - the matrix, should have been created with MatCreateShell()
32: Output Parameter:
33: . ctx - the user provided context
35: Level: advanced
37: Notes:
38: This routine is intended for use within various shell matrix routines,
39: as set with MatShellSetOperation().
40:
41: .keywords: matrix, shell, get, context
43: .seealso: MatCreateShell(), MatShellSetOperation(), MatShellSetContext()
44: @*/
45: PetscErrorCode MatShellGetContext(Mat mat,void **ctx)
46: {
48: PetscTruth flg;
53: PetscTypeCompare((PetscObject)mat,MATSHELL,&flg);
54: if (!flg) *ctx = 0;
55: else *ctx = ((Mat_Shell*)(mat->data))->ctx;
56: return(0);
57: }
61: PetscErrorCode MatDestroy_Shell(Mat mat)
62: {
64: Mat_Shell *shell;
67: shell = (Mat_Shell*)mat->data;
68: if (shell->destroy) {(*shell->destroy)(mat);}
69: PetscFree(shell);
70: return(0);
71: }
75: PetscErrorCode MatMult_Shell(Mat A,Vec x,Vec y)
76: {
77: Mat_Shell *shell = (Mat_Shell*)A->data;
81: (*shell->mult)(A,x,y);
82: if (shell->shift && shell->scale) {
83: VecAXPBY(y,shell->vshift,shell->vscale,x);
84: } else if (shell->scale) {
85: VecScale(y,shell->vscale);
86: } else {
87: VecAXPY(y,shell->vshift,x);
88: }
89: return(0);
90: }
94: PetscErrorCode MatMultTranspose_Shell(Mat A,Vec x,Vec y)
95: {
96: Mat_Shell *shell = (Mat_Shell*)A->data;
100: (*shell->multtranspose)(A,x,y);
101: if (shell->shift && shell->scale) {
102: VecAXPBY(y,shell->vshift,shell->vscale,x);
103: } else if (shell->scale) {
104: VecScale(y,shell->vscale);
105: } else {
106: VecAXPY(y,shell->vshift,x);
107: }
108: return(0);
109: }
113: PetscErrorCode MatGetDiagonal_Shell(Mat A,Vec v)
114: {
115: Mat_Shell *shell = (Mat_Shell*)A->data;
119: (*shell->getdiagonal)(A,v);
120: if (shell->scale) {
121: VecScale(v,shell->vscale);
122: }
123: if (shell->shift) {
124: VecShift(v,shell->vshift);
125: }
126: return(0);
127: }
131: PetscErrorCode MatShift_Shell(Mat Y,PetscScalar a)
132: {
133: Mat_Shell *shell = (Mat_Shell*)Y->data;
136: if (shell->scale || shell->shift) {
137: shell->vshift += a;
138: } else {
139: shell->mult = Y->ops->mult;
140: Y->ops->mult = MatMult_Shell;
141: if (Y->ops->multtranspose) {
142: shell->multtranspose = Y->ops->multtranspose;
143: Y->ops->multtranspose = MatMultTranspose_Shell;
144: }
145: if (Y->ops->getdiagonal) {
146: shell->getdiagonal = Y->ops->getdiagonal;
147: Y->ops->getdiagonal = MatGetDiagonal_Shell;
148: }
149: shell->vshift = a;
150: }
151: shell->shift = PETSC_TRUE;
152: return(0);
153: }
157: PetscErrorCode MatScale_Shell(Mat Y,PetscScalar a)
158: {
159: Mat_Shell *shell = (Mat_Shell*)Y->data;
162: if (shell->scale || shell->shift) {
163: shell->vscale *= a;
164: } else {
165: shell->mult = Y->ops->mult;
166: Y->ops->mult = MatMult_Shell;
167: if (Y->ops->multtranspose) {
168: shell->multtranspose = Y->ops->multtranspose;
169: Y->ops->multtranspose = MatMultTranspose_Shell;
170: }
171: if (Y->ops->getdiagonal) {
172: shell->getdiagonal = Y->ops->getdiagonal;
173: Y->ops->getdiagonal = MatGetDiagonal_Shell;
174: }
175: shell->vscale = a;
176: }
177: shell->scale = PETSC_TRUE;
178: return(0);
179: }
183: PetscErrorCode MatAssemblyEnd_Shell(Mat Y,MatAssemblyType t)
184: {
185: Mat_Shell *shell = (Mat_Shell*)Y->data;
188: if ((shell->shift || shell->scale) && t == MAT_FINAL_ASSEMBLY) {
189: shell->scale = PETSC_FALSE;
190: shell->shift = PETSC_FALSE;
191: shell->vshift = 0.0;
192: shell->vscale = 1.0;
193: Y->ops->mult = shell->mult;
194: Y->ops->multtranspose = shell->multtranspose;
195: Y->ops->getdiagonal = shell->getdiagonal;
196: }
197: return(0);
198: }
200: EXTERN PetscErrorCode MatConvert_Shell(Mat, MatType,MatReuse,Mat*);
204: PetscErrorCode MatSetBlockSize_Shell(Mat A,PetscInt bs)
205: {
207: return(0);
208: }
210: static struct _MatOps MatOps_Values = {0,
211: 0,
212: 0,
213: 0,
214: /* 4*/ 0,
215: 0,
216: 0,
217: 0,
218: 0,
219: 0,
220: /*10*/ 0,
221: 0,
222: 0,
223: 0,
224: 0,
225: /*15*/ 0,
226: 0,
227: 0,
228: 0,
229: 0,
230: /*20*/ 0,
231: MatAssemblyEnd_Shell,
232: 0,
233: 0,
234: 0,
235: /*25*/ 0,
236: 0,
237: 0,
238: 0,
239: 0,
240: /*30*/ 0,
241: 0,
242: 0,
243: 0,
244: 0,
245: /*35*/ 0,
246: 0,
247: 0,
248: 0,
249: 0,
250: /*40*/ 0,
251: 0,
252: 0,
253: 0,
254: 0,
255: /*45*/ 0,
256: MatScale_Shell,
257: MatShift_Shell,
258: 0,
259: 0,
260: /*50*/ MatSetBlockSize_Shell,
261: 0,
262: 0,
263: 0,
264: 0,
265: /*55*/ 0,
266: 0,
267: 0,
268: 0,
269: 0,
270: /*60*/ 0,
271: MatDestroy_Shell,
272: 0,
273: 0,
274: 0,
275: /*65*/ 0,
276: 0,
277: 0,
278: 0,
279: 0,
280: /*70*/ 0,
281: MatConvert_Shell,
282: 0,
283: 0,
284: 0,
285: /*75*/ 0,
286: 0,
287: 0,
288: 0,
289: 0,
290: /*80*/ 0,
291: 0,
292: 0,
293: 0,
294: 0,
295: /*85*/ 0,
296: 0,
297: 0,
298: 0,
299: 0,
300: /*90*/ 0,
301: 0,
302: 0,
303: 0,
304: 0,
305: /*95*/ 0,
306: 0,
307: 0,
308: 0};
310: /*MC
311: MATSHELL - MATSHELL = "shell" - A matrix type to be used to define your own matrix type -- perhaps matrix free.
313: Level: advanced
315: .seealso: MatCreateShell
316: M*/
321: PetscErrorCode MatCreate_Shell(Mat A)
322: {
323: Mat_Shell *b;
327: PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));
329: PetscNew(Mat_Shell,&b);
330: PetscLogObjectMemory(A,sizeof(struct _p_Mat)+sizeof(Mat_Shell));
331: A->data = (void*)b;
333: if (A->rmap.n == PETSC_DECIDE || A->cmap.n == PETSC_DECIDE) {
334: SETERRQ(PETSC_ERR_ARG_WRONG,"Must give local row and column count for matrix");
335: }
337: PetscMapInitialize(A->comm,&A->rmap);
338: PetscMapInitialize(A->comm,&A->cmap);
340: b->ctx = 0;
341: b->scale = PETSC_FALSE;
342: b->shift = PETSC_FALSE;
343: b->vshift = 0.0;
344: b->vscale = 1.0;
345: b->mult = 0;
346: b->multtranspose = 0;
347: b->getdiagonal = 0;
348: A->assembled = PETSC_TRUE;
349: A->preallocated = PETSC_FALSE;
350: PetscObjectChangeTypeName((PetscObject)A,MATSHELL);
351: return(0);
352: }
357: /*@C
358: MatCreateShell - Creates a new matrix class for use with a user-defined
359: private data storage format.
361: Collective on MPI_Comm
363: Input Parameters:
364: + comm - MPI communicator
365: . m - number of local rows (must be given)
366: . n - number of local columns (must be given)
367: . M - number of global rows (may be PETSC_DETERMINE)
368: . N - number of global columns (may be PETSC_DETERMINE)
369: - ctx - pointer to data needed by the shell matrix routines
371: Output Parameter:
372: . A - the matrix
374: Level: advanced
376: Usage:
378: $ MatCreateShell(comm,m,n,M,N,ctx,&mat);
379: $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
380: $ [ Use matrix for operations that have been set ]
381: $ MatDestroy(mat);
383: Notes:
384: The shell matrix type is intended to provide a simple class to use
385: with KSP (such as, for use with matrix-free methods). You should not
386: use the shell type if you plan to define a complete matrix class.
388: Fortran Notes: The context can only be an integer or a PetscObject
389: unfortunately it cannot be a Fortran array or derived type.
391: PETSc requires that matrices and vectors being used for certain
392: operations are partitioned accordingly. For example, when
393: creating a shell matrix, A, that supports parallel matrix-vector
394: products using MatMult(A,x,y) the user should set the number
395: of local matrix rows to be the number of local elements of the
396: corresponding result vector, y. Note that this is information is
397: required for use of the matrix interface routines, even though
398: the shell matrix may not actually be physically partitioned.
399: For example,
401: $
402: $ Vec x, y
404: $ Mat A
405: $
406: $ VecCreateMPI(comm,PETSC_DECIDE,M,&y);
407: $ VecCreateMPI(comm,PETSC_DECIDE,N,&x);
408: $ VecGetLocalSize(y,&m);
409: $ VecGetLocalSize(x,&n);
410: $ MatCreateShell(comm,m,n,M,N,ctx,&A);
411: $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
412: $ MatMult(A,x,y);
413: $ MatDestroy(A);
414: $ VecDestroy(y); VecDestroy(x);
415: $
417: .keywords: matrix, shell, create
419: .seealso: MatShellSetOperation(), MatHasOperation(), MatShellGetContext(), MatShellSetContext()
420: @*/
421: PetscErrorCode MatCreateShell(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,void *ctx,Mat *A)
422: {
426: MatCreate(comm,A);
427: MatSetSizes(*A,m,n,M,N);
428:
429: MatSetType(*A,MATSHELL);
430: MatShellSetContext(*A,ctx);
431: return(0);
432: }
436: /*@
437: MatShellSetContext - sets the context for a shell matrix
439: Collective on Mat
441: Input Parameters:
442: + mat - the shell matrix
443: - ctx - the context
445: Level: advanced
447: Fortran Notes: The context can only be an integer or a PetscObject
448: unfortunately it cannot be a Fortran array or derived type.
450: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation()
451: @*/
452: PetscErrorCode MatShellSetContext(Mat mat,void *ctx)
453: {
454: Mat_Shell *shell = (Mat_Shell*)mat->data;
456: PetscTruth flg;
460: PetscTypeCompare((PetscObject)mat,MATSHELL,&flg);
461: if (flg) {
462: shell->ctx = ctx;
463: }
464: return(0);
465: }
469: /*@C
470: MatShellSetOperation - Allows user to set a matrix operation for
471: a shell matrix.
473: Collective on Mat
475: Input Parameters:
476: + mat - the shell matrix
477: . op - the name of the operation
478: - f - the function that provides the operation.
480: Level: advanced
482: Usage:
484: $ MatCreateShell(comm,m,n,M,N,ctx,&A);
485: $ MatShellSetOperation(A,MATOP_MULT,(void(*)(void))usermult);
487: Notes:
488: See the file include/petscmat.h for a complete list of matrix
489: operations, which all have the form MATOP_<OPERATION>, where
490: <OPERATION> is the name (in all capital letters) of the
491: user interface routine (e.g., MatMult() -> MATOP_MULT).
493: All user-provided functions should have the same calling
494: sequence as the usual matrix interface routines, since they
495: are intended to be accessed via the usual matrix interface
496: routines, e.g.,
497: $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec)
499: Within each user-defined routine, the user should call
500: MatShellGetContext() to obtain the user-defined context that was
501: set by MatCreateShell().
503: .keywords: matrix, shell, set, operation
505: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext()
506: @*/
507: PetscErrorCode MatShellSetOperation(Mat mat,MatOperation op,void (*f)(void))
508: {
510: PetscTruth flg;
514: if (op == MATOP_DESTROY) {
515: PetscTypeCompare((PetscObject)mat,MATSHELL,&flg);
516: if (flg) {
517: Mat_Shell *shell = (Mat_Shell*)mat->data;
518: shell->destroy = (PetscErrorCode (*)(Mat)) f;
519: } else mat->ops->destroy = (PetscErrorCode (*)(Mat)) f;
520: }
521: else if (op == MATOP_VIEW) mat->ops->view = (PetscErrorCode (*)(Mat,PetscViewer)) f;
522: else (((void(**)(void))mat->ops)[op]) = f;
524: return(0);
525: }
529: /*@C
530: MatShellGetOperation - Gets a matrix function for a shell matrix.
532: Not Collective
534: Input Parameters:
535: + mat - the shell matrix
536: - op - the name of the operation
538: Output Parameter:
539: . f - the function that provides the operation.
541: Level: advanced
543: Notes:
544: See the file include/petscmat.h for a complete list of matrix
545: operations, which all have the form MATOP_<OPERATION>, where
546: <OPERATION> is the name (in all capital letters) of the
547: user interface routine (e.g., MatMult() -> MATOP_MULT).
549: All user-provided functions have the same calling
550: sequence as the usual matrix interface routines, since they
551: are intended to be accessed via the usual matrix interface
552: routines, e.g.,
553: $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec)
555: Within each user-defined routine, the user should call
556: MatShellGetContext() to obtain the user-defined context that was
557: set by MatCreateShell().
559: .keywords: matrix, shell, set, operation
561: .seealso: MatCreateShell(), MatShellGetContext(), MatShellSetOperation(), MatShellSetContext()
562: @*/
563: PetscErrorCode MatShellGetOperation(Mat mat,MatOperation op,void(**f)(void))
564: {
566: PetscTruth flg;
570: if (op == MATOP_DESTROY) {
571: PetscTypeCompare((PetscObject)mat,MATSHELL,&flg);
572: if (flg) {
573: Mat_Shell *shell = (Mat_Shell*)mat->data;
574: *f = (void(*)(void))shell->destroy;
575: } else {
576: *f = (void(*)(void))mat->ops->destroy;
577: }
578: } else if (op == MATOP_VIEW) {
579: *f = (void(*)(void))mat->ops->view;
580: } else {
581: *f = (((void(**)(void))mat->ops)[op]);
582: }
584: return(0);
585: }