Actual source code: shell.c

petsc-master 2016-08-23
Report Typos and Errors
  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*/

 10: typedef struct {
 11:   PetscErrorCode (*destroy)(Mat);
 12:   PetscErrorCode (*mult)(Mat,Vec,Vec);
 13:   PetscErrorCode (*multtranspose)(Mat,Vec,Vec);
 14:   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:   Vec         left_add_work,right_add_work;
 22:   PetscBool   usingscaled;
 23:   void        *ctx;
 24: } Mat_Shell;
 25: /*
 26:  The most general expression for the matrix is

 28:  S = L (a A + B) R

 30:  where
 31:  A is the matrix defined by the user's function
 32:  a is a scalar multiple
 33:  L is left scaling
 34:  R is right scaling
 35:  B is a diagonal shift defined by
 36:    diag(dshift) if the vector dshift is non-NULL
 37:    vscale*identity otherwise

 39:  The following identities apply:

 41:  Scale by c:
 42:   c [L (a A + B) R] = L [(a c) A + c B] R

 44:  Shift by c:
 45:   [L (a A + B) R] + c = L [a A + (B + c Linv Rinv)] R

 47:  Diagonal scaling is achieved by simply multiplying with existing L and R vectors

 49:  In the data structure:

 51:  vscale=1.0  means no special scaling will be applied
 52:  dshift=NULL means a constant diagonal shift (fall back to vshift)
 53:  vshift=0.0  means no constant diagonal shift, note that vshift is only used if dshift is NULL
 54: */

 56: static PetscErrorCode MatMult_Shell(Mat,Vec,Vec);
 57: static PetscErrorCode MatMultTranspose_Shell(Mat,Vec,Vec);
 58: static PetscErrorCode MatGetDiagonal_Shell(Mat,Vec);

 62: static PetscErrorCode MatShellUseScaledMethods(Mat Y)
 63: {
 64:   Mat_Shell *shell = (Mat_Shell*)Y->data;

 67:   if (shell->usingscaled) return(0);
 68:   shell->mult  = Y->ops->mult;
 69:   Y->ops->mult = MatMult_Shell;
 70:   if (Y->ops->multtranspose) {
 71:     shell->multtranspose  = Y->ops->multtranspose;
 72:     Y->ops->multtranspose = MatMultTranspose_Shell;
 73:   }
 74:   if (Y->ops->getdiagonal) {
 75:     shell->getdiagonal  = Y->ops->getdiagonal;
 76:     Y->ops->getdiagonal = MatGetDiagonal_Shell;
 77:   }
 78:   shell->usingscaled = PETSC_TRUE;
 79:   return(0);
 80: }

 84: static PetscErrorCode MatShellPreScaleLeft(Mat A,Vec x,Vec *xx)
 85: {
 86:   Mat_Shell      *shell = (Mat_Shell*)A->data;

 90:   *xx = NULL;
 91:   if (!shell->left) {
 92:     *xx = x;
 93:   } else {
 94:     if (!shell->left_work) {VecDuplicate(shell->left,&shell->left_work);}
 95:     VecPointwiseMult(shell->left_work,x,shell->left);
 96:     *xx  = shell->left_work;
 97:   }
 98:   return(0);
 99: }

103: static PetscErrorCode MatShellPreScaleRight(Mat A,Vec x,Vec *xx)
104: {
105:   Mat_Shell      *shell = (Mat_Shell*)A->data;

109:   *xx = NULL;
110:   if (!shell->right) {
111:     *xx = x;
112:   } else {
113:     if (!shell->right_work) {VecDuplicate(shell->right,&shell->right_work);}
114:     VecPointwiseMult(shell->right_work,x,shell->right);
115:     *xx  = shell->right_work;
116:   }
117:   return(0);
118: }

122: static PetscErrorCode MatShellPostScaleLeft(Mat A,Vec x)
123: {
124:   Mat_Shell      *shell = (Mat_Shell*)A->data;

128:   if (shell->left) {VecPointwiseMult(x,x,shell->left);}
129:   return(0);
130: }

134: static PetscErrorCode MatShellPostScaleRight(Mat A,Vec x)
135: {
136:   Mat_Shell      *shell = (Mat_Shell*)A->data;

140:   if (shell->right) {VecPointwiseMult(x,x,shell->right);}
141:   return(0);
142: }

146: static PetscErrorCode MatShellShiftAndScale(Mat A,Vec X,Vec Y)
147: {
148:   Mat_Shell      *shell = (Mat_Shell*)A->data;

152:   if (shell->dshift) {          /* get arrays because there is no VecPointwiseMultAdd() */
153:     PetscInt          i,m;
154:     const PetscScalar *x,*d;
155:     PetscScalar       *y;
156:     VecGetLocalSize(X,&m);
157:     VecGetArrayRead(shell->dshift,&d);
158:     VecGetArrayRead(X,&x);
159:     VecGetArray(Y,&y);
160:     for (i=0; i<m; i++) y[i] = shell->vscale*y[i] + d[i]*x[i];
161:     VecRestoreArrayRead(shell->dshift,&d);
162:     VecRestoreArrayRead(X,&x);
163:     VecRestoreArray(Y,&y);
164:   } else if (PetscAbsScalar(shell->vshift) != 0) {
165:     VecAXPBY(Y,shell->vshift,shell->vscale,X);
166:   } else if (shell->vscale != 1.0) {
167:     VecScale(Y,shell->vscale);
168:   }
169:   return(0);
170: }

174: /*@
175:     MatShellGetContext - Returns the user-provided context associated with a shell matrix.

177:     Not Collective

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

182:     Output Parameter:
183: .   ctx - the user provided context

185:     Level: advanced

187:    Fortran Notes: To use this from Fortran you must write a Fortran interface definition for this
188:     function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.

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(PetscObjectComm((PetscObject)mat),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:   VecDestroy(&shell->left_add_work);
225:   VecDestroy(&shell->right_add_work);
226:   PetscFree(mat->data);
227:   return(0);
228: }

232: PetscErrorCode MatMult_Shell(Mat A,Vec x,Vec y)
233: {
234:   Mat_Shell        *shell = (Mat_Shell*)A->data;
235:   PetscErrorCode   ierr;
236:   Vec              xx;
237:   PetscObjectState instate,outstate;

240:   MatShellPreScaleRight(A,x,&xx);
241:   PetscObjectStateGet((PetscObject)y, &instate);
242:   (*shell->mult)(A,xx,y);
243:   PetscObjectStateGet((PetscObject)y, &outstate);
244:   if (instate == outstate) {
245:     /* increase the state of the output vector since the user did not update its state themself as should have been done */
246:     PetscObjectStateIncrease((PetscObject)y);
247:   }
248:   MatShellShiftAndScale(A,xx,y);
249:   MatShellPostScaleLeft(A,y);
250:   return(0);
251: }

255: PetscErrorCode MatMultAdd_Shell(Mat A,Vec x,Vec y,Vec z)
256: {
257:   Mat_Shell      *shell = (Mat_Shell*)A->data;

261:   if (y == z) {
262:     if (!shell->right_add_work) {VecDuplicate(z,&shell->right_add_work);}
263:     MatMult(A,x,shell->right_add_work);
264:     VecAXPY(z,1.0,shell->right_add_work);
265:   } else {
266:     MatMult(A,x,z);
267:     VecAXPY(z,1.0,y);
268:   }
269:   return(0);
270: }

274: PetscErrorCode MatMultTranspose_Shell(Mat A,Vec x,Vec y)
275: {
276:   Mat_Shell        *shell = (Mat_Shell*)A->data;
277:   PetscErrorCode   ierr;
278:   Vec              xx;
279:   PetscObjectState instate,outstate;

282:   MatShellPreScaleLeft(A,x,&xx);
283:   PetscObjectStateGet((PetscObject)y, &instate);
284:   (*shell->multtranspose)(A,xx,y);
285:   PetscObjectStateGet((PetscObject)y, &outstate);
286:   if (instate == outstate) {
287:     /* increase the state of the output vector since the user did not update its state themself as should have been done */
288:     PetscObjectStateIncrease((PetscObject)y);
289:   }
290:   MatShellShiftAndScale(A,xx,y);
291:   MatShellPostScaleRight(A,y);
292:   return(0);
293: }

297: PetscErrorCode MatMultTransposeAdd_Shell(Mat A,Vec x,Vec y,Vec z)
298: {
299:   Mat_Shell      *shell = (Mat_Shell*)A->data;

303:   if (y == z) {
304:     if (!shell->left_add_work) {VecDuplicate(z,&shell->left_add_work);}
305:     MatMultTranspose(A,x,shell->left_add_work);
306:     VecWAXPY(z,1.0,shell->left_add_work,y);
307:   } else {
308:     MatMultTranspose(A,x,z);
309:     VecAXPY(z,1.0,y);
310:   }
311:   return(0);
312: }

316: PetscErrorCode MatGetDiagonal_Shell(Mat A,Vec v)
317: {
318:   Mat_Shell      *shell = (Mat_Shell*)A->data;

322:   (*shell->getdiagonal)(A,v);
323:   VecScale(v,shell->vscale);
324:   if (shell->dshift) {
325:     VecPointwiseMult(v,v,shell->dshift);
326:   } else {
327:     VecShift(v,shell->vshift);
328:   }
329:   if (shell->left)  {VecPointwiseMult(v,v,shell->left);}
330:   if (shell->right) {VecPointwiseMult(v,v,shell->right);}
331:   return(0);
332: }

336: PetscErrorCode MatShift_Shell(Mat Y,PetscScalar a)
337: {
338:   Mat_Shell      *shell = (Mat_Shell*)Y->data;

342:   if (shell->left || shell->right || shell->dshift) {
343:     if (!shell->dshift) {
344:       if (!shell->dshift_owned) {VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift_owned);}
345:       shell->dshift = shell->dshift_owned;
346:       VecSet(shell->dshift,shell->vshift+a);
347:     } else {VecScale(shell->dshift,a);}
348:     if (shell->left)  {VecPointwiseDivide(shell->dshift,shell->dshift,shell->left);}
349:     if (shell->right) {VecPointwiseDivide(shell->dshift,shell->dshift,shell->right);}
350:   } else shell->vshift += a;
351:   MatShellUseScaledMethods(Y);
352:   return(0);
353: }

357: PetscErrorCode MatScale_Shell(Mat Y,PetscScalar a)
358: {
359:   Mat_Shell      *shell = (Mat_Shell*)Y->data;

363:   shell->vscale *= a;
364:   if (shell->dshift) {
365:     VecScale(shell->dshift,a);
366:   } else shell->vshift *= a;
367:   MatShellUseScaledMethods(Y);
368:   return(0);
369: }

373: static PetscErrorCode MatDiagonalScale_Shell(Mat Y,Vec left,Vec right)
374: {
375:   Mat_Shell      *shell = (Mat_Shell*)Y->data;

379:   if (left) {
380:     if (!shell->left_owned) {VecDuplicate(left,&shell->left_owned);}
381:     if (shell->left) {
382:       VecPointwiseMult(shell->left,shell->left,left);
383:     } else {
384:       shell->left = shell->left_owned;
385:       VecCopy(left,shell->left);
386:     }
387:   }
388:   if (right) {
389:     if (!shell->right_owned) {VecDuplicate(right,&shell->right_owned);}
390:     if (shell->right) {
391:       VecPointwiseMult(shell->right,shell->right,right);
392:     } else {
393:       shell->right = shell->right_owned;
394:       VecCopy(right,shell->right);
395:     }
396:   }
397:   MatShellUseScaledMethods(Y);
398:   return(0);
399: }

403: PetscErrorCode MatAssemblyEnd_Shell(Mat Y,MatAssemblyType t)
404: {
405:   Mat_Shell *shell = (Mat_Shell*)Y->data;

408:   if (t == MAT_FINAL_ASSEMBLY) {
409:     shell->vshift = 0.0;
410:     shell->vscale = 1.0;
411:     shell->dshift = NULL;
412:     shell->left   = NULL;
413:     shell->right  = NULL;
414:     if (shell->mult) {
415:       Y->ops->mult = shell->mult;
416:       shell->mult  = NULL;
417:     }
418:     if (shell->multtranspose) {
419:       Y->ops->multtranspose = shell->multtranspose;
420:       shell->multtranspose  = NULL;
421:     }
422:     if (shell->getdiagonal) {
423:       Y->ops->getdiagonal = shell->getdiagonal;
424:       shell->getdiagonal  = NULL;
425:     }
426:     shell->usingscaled = PETSC_FALSE;
427:   }
428:   return(0);
429: }

431: PETSC_INTERN PetscErrorCode MatConvert_Shell(Mat, MatType,MatReuse,Mat*);

435: static PetscErrorCode MatMissingDiagonal_Shell(Mat A,PetscBool  *missing,PetscInt *d)
436: {
438:   *missing = PETSC_FALSE;
439:   return(0);
440: }

442: static struct _MatOps MatOps_Values = {0,
443:                                        0,
444:                                        0,
445:                                        0,
446:                                 /* 4*/ 0,
447:                                        0,
448:                                        0,
449:                                        0,
450:                                        0,
451:                                        0,
452:                                 /*10*/ 0,
453:                                        0,
454:                                        0,
455:                                        0,
456:                                        0,
457:                                 /*15*/ 0,
458:                                        0,
459:                                        0,
460:                                        MatDiagonalScale_Shell,
461:                                        0,
462:                                 /*20*/ 0,
463:                                        MatAssemblyEnd_Shell,
464:                                        0,
465:                                        0,
466:                                 /*24*/ 0,
467:                                        0,
468:                                        0,
469:                                        0,
470:                                        0,
471:                                 /*29*/ 0,
472:                                        0,
473:                                        0,
474:                                        0,
475:                                        0,
476:                                 /*34*/ 0,
477:                                        0,
478:                                        0,
479:                                        0,
480:                                        0,
481:                                 /*39*/ 0,
482:                                        0,
483:                                        0,
484:                                        0,
485:                                        0,
486:                                 /*44*/ 0,
487:                                        MatScale_Shell,
488:                                        MatShift_Shell,
489:                                        0,
490:                                        0,
491:                                 /*49*/ 0,
492:                                        0,
493:                                        0,
494:                                        0,
495:                                        0,
496:                                 /*54*/ 0,
497:                                        0,
498:                                        0,
499:                                        0,
500:                                        0,
501:                                 /*59*/ 0,
502:                                        MatDestroy_Shell,
503:                                        0,
504:                                        0,
505:                                        0,
506:                                 /*64*/ 0,
507:                                        0,
508:                                        0,
509:                                        0,
510:                                        0,
511:                                 /*69*/ 0,
512:                                        0,
513:                                        MatConvert_Shell,
514:                                        0,
515:                                        0,
516:                                 /*74*/ 0,
517:                                        0,
518:                                        0,
519:                                        0,
520:                                        0,
521:                                 /*79*/ 0,
522:                                        0,
523:                                        0,
524:                                        0,
525:                                        0,
526:                                 /*84*/ 0,
527:                                        0,
528:                                        0,
529:                                        0,
530:                                        0,
531:                                 /*89*/ 0,
532:                                        0,
533:                                        0,
534:                                        0,
535:                                        0,
536:                                 /*94*/ 0,
537:                                        0,
538:                                        0,
539:                                        0,
540:                                        0,
541:                                 /*99*/ 0,
542:                                        0,
543:                                        0,
544:                                        0,
545:                                        0,
546:                                /*104*/ 0,
547:                                        0,
548:                                        0,
549:                                        0,
550:                                        0,
551:                                /*109*/ 0,
552:                                        0,
553:                                        0,
554:                                        0,
555:                                        MatMissingDiagonal_Shell,
556:                                /*114*/ 0,
557:                                        0,
558:                                        0,
559:                                        0,
560:                                        0,
561:                                /*119*/ 0,
562:                                        0,
563:                                        0,
564:                                        0,
565:                                        0,
566:                                /*124*/ 0,
567:                                        0,
568:                                        0,
569:                                        0,
570:                                        0,
571:                                /*129*/ 0,
572:                                        0,
573:                                        0,
574:                                        0,
575:                                        0,
576:                                /*134*/ 0,
577:                                        0,
578:                                        0,
579:                                        0,
580:                                        0,
581:                                /*139*/ 0,
582:                                        0,
583:                                        0
584: };

586: /*MC
587:    MATSHELL - MATSHELL = "shell" - A matrix type to be used to define your own matrix type -- perhaps matrix free.

589:   Level: advanced

591: .seealso: MatCreateShell
592: M*/

596: PETSC_EXTERN PetscErrorCode MatCreate_Shell(Mat A)
597: {
598:   Mat_Shell      *b;

602:   PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));

604:   PetscNewLog(A,&b);
605:   A->data = (void*)b;

607:   PetscLayoutSetUp(A->rmap);
608:   PetscLayoutSetUp(A->cmap);

610:   b->ctx           = 0;
611:   b->vshift        = 0.0;
612:   b->vscale        = 1.0;
613:   b->mult          = 0;
614:   b->multtranspose = 0;
615:   b->getdiagonal   = 0;
616:   A->assembled     = PETSC_TRUE;
617:   A->preallocated  = PETSC_FALSE;

619:   PetscObjectChangeTypeName((PetscObject)A,MATSHELL);
620:   return(0);
621: }

625: /*@C
626:    MatCreateShell - Creates a new matrix class for use with a user-defined
627:    private data storage format.

629:   Collective on MPI_Comm

631:    Input Parameters:
632: +  comm - MPI communicator
633: .  m - number of local rows (must be given)
634: .  n - number of local columns (must be given)
635: .  M - number of global rows (may be PETSC_DETERMINE)
636: .  N - number of global columns (may be PETSC_DETERMINE)
637: -  ctx - pointer to data needed by the shell matrix routines

639:    Output Parameter:
640: .  A - the matrix

642:    Level: advanced

644:   Usage:
645: $    extern int mult(Mat,Vec,Vec);
646: $    MatCreateShell(comm,m,n,M,N,ctx,&mat);
647: $    MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
648: $    [ Use matrix for operations that have been set ]
649: $    MatDestroy(mat);

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

656:    Fortran Notes: To use this from Fortran with a ctx you must write an interface definition for this
657:     function and for MatShellGetContext() that tells Fortran the Fortran derived data type you are passing
658:     in as the ctx argument.

660:    PETSc requires that matrices and vectors being used for certain
661:    operations are partitioned accordingly.  For example, when
662:    creating a shell matrix, A, that supports parallel matrix-vector
663:    products using MatMult(A,x,y) the user should set the number
664:    of local matrix rows to be the number of local elements of the
665:    corresponding result vector, y. Note that this is information is
666:    required for use of the matrix interface routines, even though
667:    the shell matrix may not actually be physically partitioned.
668:    For example,

670: $
671: $     Vec x, y
672: $     extern int mult(Mat,Vec,Vec);
673: $     Mat A
674: $
675: $     VecCreateMPI(comm,PETSC_DECIDE,M,&y);
676: $     VecCreateMPI(comm,PETSC_DECIDE,N,&x);
677: $     VecGetLocalSize(y,&m);
678: $     VecGetLocalSize(x,&n);
679: $     MatCreateShell(comm,m,n,M,N,ctx,&A);
680: $     MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
681: $     MatMult(A,x,y);
682: $     MatDestroy(A);
683: $     VecDestroy(y); VecDestroy(x);
684: $

686: .keywords: matrix, shell, create

688: .seealso: MatShellSetOperation(), MatHasOperation(), MatShellGetContext(), MatShellSetContext()
689: @*/
690: PetscErrorCode  MatCreateShell(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,void *ctx,Mat *A)
691: {

695:   MatCreate(comm,A);
696:   MatSetSizes(*A,m,n,M,N);
697:   MatSetType(*A,MATSHELL);
698:   MatShellSetContext(*A,ctx);
699:   MatSetUp(*A);
700:   return(0);
701: }

705: /*@
706:     MatShellSetContext - sets the context for a shell matrix

708:    Logically Collective on Mat

710:     Input Parameters:
711: +   mat - the shell matrix
712: -   ctx - the context

714:    Level: advanced

716:    Fortran Notes: To use this from Fortran you must write a Fortran interface definition for this
717:     function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.

719: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation()
720: @*/
721: PetscErrorCode  MatShellSetContext(Mat mat,void *ctx)
722: {
723:   Mat_Shell      *shell = (Mat_Shell*)mat->data;
725:   PetscBool      flg;

729:   PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
730:   if (flg) {
731:     shell->ctx = ctx;
732:   } else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot attach context to non-shell matrix");
733:   return(0);
734: }

738: /*@C
739:     MatShellSetOperation - Allows user to set a matrix operation for
740:                            a shell matrix.

742:    Logically Collective on Mat

744:     Input Parameters:
745: +   mat - the shell matrix
746: .   op - the name of the operation
747: -   f - the function that provides the operation.

749:    Level: advanced

751:     Usage:
752: $      extern PetscErrorCode usermult(Mat,Vec,Vec);
753: $      MatCreateShell(comm,m,n,M,N,ctx,&A);
754: $      MatShellSetOperation(A,MATOP_MULT,(void(*)(void))usermult);

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

762:     All user-provided functions (execept for MATOP_DESTROY) should have the same calling
763:     sequence as the usual matrix interface routines, since they
764:     are intended to be accessed via the usual matrix interface
765:     routines, e.g.,
766: $       MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec)

768:     In particular each function MUST return an error code of 0 on success and
769:     nonzero on failure.

771:     Within each user-defined routine, the user should call
772:     MatShellGetContext() to obtain the user-defined context that was
773:     set by MatCreateShell().

775:     Fortran Notes: For MatCreateVecs() the user code should check if the input left or right matrix is -1 and in that case not
776:        generate a matrix. See src/mat/examples/tests/ex120f.F

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

780: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext()
781: @*/
782: PetscErrorCode  MatShellSetOperation(Mat mat,MatOperation op,void (*f)(void))
783: {
785:   PetscBool      flg;

789:   switch (op) {
790:   case MATOP_DESTROY:
791:     PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
792:     if (flg) {
793:       Mat_Shell *shell = (Mat_Shell*)mat->data;
794:       shell->destroy = (PetscErrorCode (*)(Mat))f;
795:     } else mat->ops->destroy = (PetscErrorCode (*)(Mat))f;
796:     break;
797:   case MATOP_VIEW:
798:     mat->ops->view = (PetscErrorCode (*)(Mat,PetscViewer))f;
799:     break;
800:   case MATOP_MULT:
801:     mat->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f;
802:     if (!mat->ops->multadd) mat->ops->multadd = MatMultAdd_Shell;
803:     break;
804:   case MATOP_MULT_TRANSPOSE:
805:     mat->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f;
806:     if (!mat->ops->multtransposeadd) mat->ops->multtransposeadd = MatMultTransposeAdd_Shell;
807:     break;
808:   default:
809:     (((void(**)(void))mat->ops)[op]) = f;
810:   }
811:   return(0);
812: }

816: /*@C
817:     MatShellGetOperation - Gets a matrix function for a shell matrix.

819:     Not Collective

821:     Input Parameters:
822: +   mat - the shell matrix
823: -   op - the name of the operation

825:     Output Parameter:
826: .   f - the function that provides the operation.

828:     Level: advanced

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

836:     All user-provided functions have the same calling
837:     sequence as the usual matrix interface routines, since they
838:     are intended to be accessed via the usual matrix interface
839:     routines, e.g.,
840: $       MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec)

842:     Within each user-defined routine, the user should call
843:     MatShellGetContext() to obtain the user-defined context that was
844:     set by MatCreateShell().

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

848: .seealso: MatCreateShell(), MatShellGetContext(), MatShellSetOperation(), MatShellSetContext()
849: @*/
850: PetscErrorCode  MatShellGetOperation(Mat mat,MatOperation op,void(**f)(void))
851: {
853:   PetscBool      flg;

857:   if (op == MATOP_DESTROY) {
858:     PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
859:     if (flg) {
860:       Mat_Shell *shell = (Mat_Shell*)mat->data;
861:       *f = (void (*)(void))shell->destroy;
862:     } else {
863:       *f = (void (*)(void))mat->ops->destroy;
864:     }
865:   } else if (op == MATOP_VIEW) {
866:     *f = (void (*)(void))mat->ops->view;
867:   } else {
868:     *f = (((void (**)(void))mat->ops)[op]);
869:   }
870:   return(0);
871: }