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: }