Actual source code: shell.c

  1: /*$Id: shell.c,v 1.88 2001/09/07 20:09:41 bsmith Exp $*/

  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 src/mat/matimpl.h
 10:  #include src/vec/vecimpl.h

 12: typedef struct {
 13:   int         (*destroy)(Mat);
 14:   int         (*mult)(Mat,Vec,Vec);
 15:   PetscTruth  scale,shift;
 16:   PetscScalar vscale,vshift;
 17:   void        *ctx;
 18: } Mat_Shell;

 20: #undef __FUNCT__  
 22: /*@
 23:     MatShellGetContext - Returns the user-provided context associated with a shell matrix.

 25:     Not Collective

 27:     Input Parameter:
 28: .   mat - the matrix, should have been created with MatCreateShell()

 30:     Output Parameter:
 31: .   ctx - the user provided context

 33:     Level: advanced

 35:     Notes:
 36:     This routine is intended for use within various shell matrix routines,
 37:     as set with MatShellSetOperation().
 38:     
 39: .keywords: matrix, shell, get, context

 41: .seealso: MatCreateShell(), MatShellSetOperation(), MatShellSetContext()
 42: @*/
 43: int MatShellGetContext(Mat mat,void **ctx)
 44: {
 45:   int        ierr;
 46:   PetscTruth flg;

 50:   PetscTypeCompare((PetscObject)mat,MATSHELL,&flg);
 51:   if (!flg) *ctx = 0;
 52:   else      *ctx = ((Mat_Shell*)(mat->data))->ctx;
 53:   return(0);
 54: }

 56: #undef __FUNCT__  
 58: int MatDestroy_Shell(Mat mat)
 59: {
 60:   int       ierr;
 61:   Mat_Shell *shell;

 64:   shell = (Mat_Shell*)mat->data;
 65:   if (shell->destroy) {(*shell->destroy)(mat);}
 66:   PetscFree(shell);
 67:   return(0);
 68: }

 70: #undef __FUNCT__  
 72: int MatMult_Shell(Mat A,Vec x,Vec y)
 73: {
 74:   Mat_Shell   *shell = (Mat_Shell*)A->data;
 75:   int         ierr;

 78:   (*shell->mult)(A,x,y);
 79:   if (shell->shift && shell->scale) {
 80:     VecAXPBY(&shell->vshift,&shell->vscale,x,y);
 81:   } else if (shell->scale) {
 82:     VecScale(&shell->vscale,y);
 83:   } else {
 84:     VecAXPY(&shell->vshift,x,y);
 85:   }
 86:   return(0);
 87: }

 89: #undef __FUNCT__  
 91: int MatShift_Shell(PetscScalar *a,Mat Y)
 92: {
 93:   Mat_Shell *shell = (Mat_Shell*)Y->data;
 95:   if (shell->scale || shell->shift) {
 96:     shell->vshift += *a;
 97:   } else {
 98:     shell->mult   = Y->ops->mult;
 99:     Y->ops->mult  = MatMult_Shell;
100:     shell->vshift = *a;
101:   }
102:   shell->shift  =  PETSC_TRUE;
103:   return(0);
104: }

106: #undef __FUNCT__  
108: int MatScale_Shell(PetscScalar *a,Mat Y)
109: {
110:   Mat_Shell *shell = (Mat_Shell*)Y->data;
112:   if (shell->scale || shell->shift) {
113:     shell->vscale *= *a;
114:   } else {
115:     shell->mult   = Y->ops->mult;
116:     Y->ops->mult  = MatMult_Shell;
117:     shell->vscale = *a;
118:   }
119:   shell->scale  =  PETSC_TRUE;
120:   return(0);
121: }

123: #undef __FUNCT__  
125: int MatAssemblyEnd_Shell(Mat Y,MatAssemblyType t)
126: {
127:   Mat_Shell *shell = (Mat_Shell*)Y->data;

130:   if ((shell->shift || shell->scale) && t == MAT_FINAL_ASSEMBLY) {
131:     shell->scale  = PETSC_FALSE;
132:     shell->shift  = PETSC_FALSE;
133:     shell->vshift = 0.0;
134:     shell->vscale = 1.0;
135:     Y->ops->mult  = shell->mult;
136:   }
137:   return(0);
138: }

140: extern int MatConvert_Shell(Mat,MatType,Mat*);

142: static struct _MatOps MatOps_Values = {0,
143:        0,
144:        0,
145:        0,
146:        0,
147:        0,
148:        0,
149:        0,
150:        0,
151:        0,
152:        0,
153:        0,
154:        0,
155:        0,
156:        0,
157:        0,
158:        0,
159:        0,
160:        0,
161:        0,
162:        0,
163:        MatAssemblyEnd_Shell,
164:        0,
165:        0,
166:        0,
167:        0,
168:        0,
169:        0,
170:        0,
171:        0,
172:        0,
173:        0,
174:        0,
175:        0,
176:        0,
177:        0,
178:        0,
179:        0,
180:        0,
181:        0,
182:        0,
183:        0,
184:        0,
185:        0,
186:        0,
187:        0,
188:        MatScale_Shell,
189:        MatShift_Shell,
190:        0,
191:        0,
192:        0,
193:        0,
194:        0,
195:        0,
196:        0,
197:        0,
198:        0,
199:        0,
200:        0,
201:        0,
202:        0,
203:        MatDestroy_Shell,
204:        0,
205:        MatGetPetscMaps_Petsc,
206:        0,
207:        0,
208:        0,
209:        0,
210:        0,
211:        0,
212:        0,
213:        MatConvert_Shell};

215: EXTERN_C_BEGIN
216: #undef __FUNCT__  
218: int MatCreate_Shell(Mat A)
219: {
220:   Mat_Shell *b;
221:   int       ierr;

224:   ierr            = PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));

226:   PetscNew(Mat_Shell,&b);
227:   PetscLogObjectMemory(A,sizeof(struct _p_Mat)+sizeof(Mat_Shell));
228:   ierr    = PetscMemzero(b,sizeof(Mat_Shell));
229:   A->data = (void*)b;

231:   if (A->m == PETSC_DECIDE || A->n == PETSC_DECIDE) {
232:     SETERRQ(1,"Must give local row and column count for matrix");
233:   }

235:   PetscSplitOwnership(A->comm,&A->m,&A->M);
236:   PetscSplitOwnership(A->comm,&A->n,&A->N);

238:   PetscMapCreateMPI(A->comm,A->m,A->M,&A->rmap);
239:   PetscMapCreateMPI(A->comm,A->n,A->N,&A->cmap);

241:   b->ctx          = 0;
242:   b->scale        = PETSC_FALSE;
243:   b->shift        = PETSC_FALSE;
244:   b->vshift       = 0.0;
245:   b->vscale       = 1.0;
246:   b->mult         = 0;
247:   A->assembled    = PETSC_TRUE;
248:   A->preallocated = PETSC_TRUE;
249:   return(0);
250: }
251: EXTERN_C_END

253: #undef __FUNCT__  
255: /*@C
256:    MatCreateShell - Creates a new matrix class for use with a user-defined
257:    private data storage format. 

259:   Collective on MPI_Comm

261:    Input Parameters:
262: +  comm - MPI communicator
263: .  m - number of local rows (must be given)
264: .  n - number of local columns (must be given)
265: .  M - number of global rows (may be PETSC_DETERMINE)
266: .  N - number of global columns (may be PETSC_DETERMINE)
267: -  ctx - pointer to data needed by the shell matrix routines

269:    Output Parameter:
270: .  A - the matrix

272:    Level: advanced

274:   Usage:
275: $    extern int mult(Mat,Vec,Vec);
276: $    MatCreateShell(comm,m,n,M,N,ctx,&mat);
277: $    MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
278: $    [ Use matrix for operations that have been set ]
279: $    MatDestroy(mat);

281:    Notes:
282:    The shell matrix type is intended to provide a simple class to use
283:    with KSP (such as, for use with matrix-free methods). You should not
284:    use the shell type if you plan to define a complete matrix class.

286:    PETSc requires that matrices and vectors being used for certain
287:    operations are partitioned accordingly.  For example, when
288:    creating a shell matrix, A, that supports parallel matrix-vector
289:    products using MatMult(A,x,y) the user should set the number
290:    of local matrix rows to be the number of local elements of the
291:    corresponding result vector, y. Note that this is information is
292:    required for use of the matrix interface routines, even though
293:    the shell matrix may not actually be physically partitioned.
294:    For example,

296: $
297: $     Vec x, y
298: $     extern int mult(Mat,Vec,Vec);
299: $     Mat A
300: $
301: $     VecCreateMPI(comm,PETSC_DECIDE,M,&y);
302: $     VecCreateMPI(comm,PETSC_DECIDE,N,&x);
303: $     VecGetLocalSize(y,&m);
304: $     VecGetLocalSize(x,&n);
305: $     MatCreateShell(comm,m,n,M,N,ctx,&A);
306: $     MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
307: $     MatMult(A,x,y);
308: $     MatDestroy(A);
309: $     VecDestroy(y); VecDestroy(x);
310: $

312: .keywords: matrix, shell, create

314: .seealso: MatShellSetOperation(), MatHasOperation(), MatShellGetContext(), MatShellSetContext()
315: @*/
316: int MatCreateShell(MPI_Comm comm,int m,int n,int M,int N,void *ctx,Mat *A)
317: {
318:   int       ierr;

321:   MatCreate(comm,m,n,M,N,A);
322:   MatSetType(*A,MATSHELL);
323:   MatShellSetContext(*A,ctx);
324:   return(0);
325: }

327: #undef __FUNCT__  
329: /*@C
330:     MatShellSetContext - sets the context for a shell matrix

332:    Collective on Mat

334:     Input Parameters:
335: +   mat - the shell matrix
336: -   ctx - the context

338:    Level: advanced


341: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation()
342: @*/
343: int MatShellSetContext(Mat mat,void *ctx)
344: {
345:   Mat_Shell  *shell = (Mat_Shell*)mat->data;
346:   int        ierr;
347:   PetscTruth flg;

351:   PetscTypeCompare((PetscObject)mat,MATSHELL,&flg);
352:   if (flg) {
353:     shell->ctx = ctx;
354:   }
355:   return(0);
356: }

358: #undef __FUNCT__  
360: /*@C
361:     MatShellSetOperation - Allows user to set a matrix operation for
362:                            a shell matrix.

364:    Collective on Mat

366:     Input Parameters:
367: +   mat - the shell matrix
368: .   op - the name of the operation
369: -   f - the function that provides the operation.

371:    Level: advanced

373:     Usage:
374: $      extern int usermult(Mat,Vec,Vec);
375: $      MatCreateShell(comm,m,n,M,N,ctx,&A);
376: $      MatShellSetOperation(A,MATOP_MULT,(void(*)(void))usermult);

378:     Notes:
379:     See the file include/petscmat.h for a complete list of matrix
380:     operations, which all have the form MATOP_<OPERATION>, where
381:     <OPERATION> is the name (in all capital letters) of the
382:     user interface routine (e.g., MatMult() -> MATOP_MULT).

384:     All user-provided functions should have the same calling
385:     sequence as the usual matrix interface routines, since they
386:     are intended to be accessed via the usual matrix interface
387:     routines, e.g., 
388: $       MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec)

390:     Within each user-defined routine, the user should call
391:     MatShellGetContext() to obtain the user-defined context that was
392:     set by MatCreateShell().

394: .keywords: matrix, shell, set, operation

396: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext()
397: @*/
398: int MatShellSetOperation(Mat mat,MatOperation op,void (*f)(void))
399: {
400:   int        ierr;
401:   PetscTruth flg;

405:   if (op == MATOP_DESTROY) {
406:     PetscTypeCompare((PetscObject)mat,MATSHELL,&flg);
407:     if (flg) {
408:        Mat_Shell *shell = (Mat_Shell*)mat->data;
409:        shell->destroy                 = (int (*)(Mat)) f;
410:     } else mat->ops->destroy            = (int (*)(Mat)) f;
411:   }
412:   else if (op == MATOP_VIEW) mat->ops->view  = (int (*)(Mat,PetscViewer)) f;
413:   else                       (((void(**)(void))mat->ops)[op]) = f;

415:   return(0);
416: }

418: #undef __FUNCT__  
420: /*@C
421:     MatShellGetOperation - Gets a matrix function for a shell matrix.

423:     Not Collective

425:     Input Parameters:
426: +   mat - the shell matrix
427: -   op - the name of the operation

429:     Output Parameter:
430: .   f - the function that provides the operation.

432:     Level: advanced

434:     Notes:
435:     See the file include/petscmat.h for a complete list of matrix
436:     operations, which all have the form MATOP_<OPERATION>, where
437:     <OPERATION> is the name (in all capital letters) of the
438:     user interface routine (e.g., MatMult() -> MATOP_MULT).

440:     All user-provided functions have the same calling
441:     sequence as the usual matrix interface routines, since they
442:     are intended to be accessed via the usual matrix interface
443:     routines, e.g., 
444: $       MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec)

446:     Within each user-defined routine, the user should call
447:     MatShellGetContext() to obtain the user-defined context that was
448:     set by MatCreateShell().

450: .keywords: matrix, shell, set, operation

452: .seealso: MatCreateShell(), MatShellGetContext(), MatShellSetOperation(), MatShellSetContext()
453: @*/
454: int MatShellGetOperation(Mat mat,MatOperation op,void(**f)(void))
455: {
456:   int        ierr;
457:   PetscTruth flg;

461:   if (op == MATOP_DESTROY) {
462:     PetscTypeCompare((PetscObject)mat,MATSHELL,&flg);
463:     if (flg) {
464:       Mat_Shell *shell = (Mat_Shell*)mat->data;
465:       *f = (void(*)(void))shell->destroy;
466:     } else {
467:       *f = (void(*)(void))mat->ops->destroy;
468:     }
469:   } else if (op == MATOP_VIEW) {
470:     *f = (void(*)(void))mat->ops->view;
471:   } else {
472:     *f = (((void(**)(void))mat->ops)[op]);
473:   }

475:   return(0);
476: }