Actual source code: shell.c

petsc-dev 2014-08-28
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*/
  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);

 17:   PetscScalar vscale,vshift;
 18:   Vec         dshift;
 19:   Vec         left,right;
 20:   Vec         dshift_owned,left_owned,right_owned;
 21:   Vec         left_work,right_work;
 22:   Vec         left_add_work,right_add_work;
 23:   PetscBool   usingscaled;
 24:   void        *ctx;
 25: } Mat_Shell;
 26: /*
 27:  The most general expression for the matrix is

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

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

 40:  The following identities apply:

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

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

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

 50:  In the data structure:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

178:     Not Collective

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

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

186:     Level: advanced

188:     Notes:
189:     This routine is intended for use within various shell matrix routines,
190:     as set with MatShellSetOperation().

192: .keywords: matrix, shell, get, context

194: .seealso: MatCreateShell(), MatShellSetOperation(), MatShellSetContext()
195: @*/
196: PetscErrorCode  MatShellGetContext(Mat mat,void *ctx)
197: {
199:   PetscBool      flg;

204:   PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
205:   if (flg) *(void**)ctx = ((Mat_Shell*)(mat->data))->ctx;
206:   else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot get context from non-shell matrix");
207:   return(0);
208: }

212: PetscErrorCode MatDestroy_Shell(Mat mat)
213: {
215:   Mat_Shell      *shell = (Mat_Shell*)mat->data;

218:   if (shell->destroy) {
219:     (*shell->destroy)(mat);
220:   }
221:   VecDestroy(&shell->left_owned);
222:   VecDestroy(&shell->right_owned);
223:   VecDestroy(&shell->dshift_owned);
224:   VecDestroy(&shell->left_work);
225:   VecDestroy(&shell->right_work);
226:   VecDestroy(&shell->left_add_work);
227:   VecDestroy(&shell->right_add_work);
228:   PetscFree(mat->data);
229:   return(0);
230: }

234: PetscErrorCode MatMult_Shell(Mat A,Vec x,Vec y)
235: {
236:   Mat_Shell      *shell = (Mat_Shell*)A->data;
238:   Vec            xx;

241:   MatShellPreScaleRight(A,x,&xx);
242:   (*shell->mult)(A,xx,y);
243:   MatShellShiftAndScale(A,xx,y);
244:   MatShellPostScaleLeft(A,y);
245:   return(0);
246: }

250: PetscErrorCode MatMultAdd_Shell(Mat A,Vec x,Vec y,Vec z)
251: {
252:   Mat_Shell      *shell = (Mat_Shell*)A->data;

256:   if (y == z) {
257:     if (!shell->right_add_work) {VecDuplicate(z,&shell->right_add_work);}
258:     MatMult(A,x,shell->right_add_work);
259:     VecWAXPY(z,1.0,shell->right_add_work,y);
260:   } else {
261:     MatMult(A,x,z);
262:     VecAXPY(z,1.0,y);
263:   }
264:   return(0);
265: }

269: PetscErrorCode MatMultTranspose_Shell(Mat A,Vec x,Vec y)
270: {
271:   Mat_Shell      *shell = (Mat_Shell*)A->data;
273:   Vec            xx;

276:   MatShellPreScaleLeft(A,x,&xx);
277:   (*shell->multtranspose)(A,xx,y);
278:   MatShellShiftAndScale(A,xx,y);
279:   MatShellPostScaleRight(A,y);
280:   return(0);
281: }

285: PetscErrorCode MatMultTransposeAdd_Shell(Mat A,Vec x,Vec y,Vec z)
286: {
287:   Mat_Shell      *shell = (Mat_Shell*)A->data;

291:   if (y == z) {
292:     if (!shell->left_add_work) {VecDuplicate(z,&shell->left_add_work);}
293:     MatMultTranspose(A,x,shell->left_add_work);
294:     VecWAXPY(z,1.0,shell->left_add_work,y);
295:   } else {
296:     MatMultTranspose(A,x,z);
297:     VecAXPY(z,1.0,y);
298:   }
299:   return(0);
300: }

304: PetscErrorCode MatGetDiagonal_Shell(Mat A,Vec v)
305: {
306:   Mat_Shell      *shell = (Mat_Shell*)A->data;

310:   (*shell->getdiagonal)(A,v);
311:   VecScale(v,shell->vscale);
312:   if (shell->dshift) {
313:     VecPointwiseMult(v,v,shell->dshift);
314:   } else {
315:     VecShift(v,shell->vshift);
316:   }
317:   if (shell->left)  {VecPointwiseMult(v,v,shell->left);}
318:   if (shell->right) {VecPointwiseMult(v,v,shell->right);}
319:   return(0);
320: }

324: PetscErrorCode MatShift_Shell(Mat Y,PetscScalar a)
325: {
326:   Mat_Shell      *shell = (Mat_Shell*)Y->data;

330:   if (shell->left || shell->right || shell->dshift) {
331:     if (!shell->dshift) {
332:       if (!shell->dshift_owned) {VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift_owned);}
333:       shell->dshift = shell->dshift_owned;
334:       VecSet(shell->dshift,shell->vshift+a);
335:     } else {VecScale(shell->dshift,a);}
336:     if (shell->left)  {VecPointwiseDivide(shell->dshift,shell->dshift,shell->left);}
337:     if (shell->right) {VecPointwiseDivide(shell->dshift,shell->dshift,shell->right);}
338:   } else shell->vshift += a;
339:   MatShellUseScaledMethods(Y);
340:   return(0);
341: }

345: PetscErrorCode MatScale_Shell(Mat Y,PetscScalar a)
346: {
347:   Mat_Shell      *shell = (Mat_Shell*)Y->data;

351:   shell->vscale *= a;
352:   if (shell->dshift) {
353:     VecScale(shell->dshift,a);
354:   } else shell->vshift *= a;
355:   MatShellUseScaledMethods(Y);
356:   return(0);
357: }

361: static PetscErrorCode MatDiagonalScale_Shell(Mat Y,Vec left,Vec right)
362: {
363:   Mat_Shell      *shell = (Mat_Shell*)Y->data;

367:   if (left) {
368:     if (!shell->left_owned) {VecDuplicate(left,&shell->left_owned);}
369:     if (shell->left) {
370:       VecPointwiseMult(shell->left,shell->left,left);
371:     } else {
372:       shell->left = shell->left_owned;
373:       VecCopy(left,shell->left);
374:     }
375:   }
376:   if (right) {
377:     if (!shell->right_owned) {VecDuplicate(right,&shell->right_owned);}
378:     if (shell->right) {
379:       VecPointwiseMult(shell->right,shell->right,right);
380:     } else {
381:       shell->right = shell->right_owned;
382:       VecCopy(right,shell->right);
383:     }
384:   }
385:   MatShellUseScaledMethods(Y);
386:   return(0);
387: }

391: PetscErrorCode MatAssemblyEnd_Shell(Mat Y,MatAssemblyType t)
392: {
393:   Mat_Shell *shell = (Mat_Shell*)Y->data;

396:   if (t == MAT_FINAL_ASSEMBLY) {
397:     shell->vshift = 0.0;
398:     shell->vscale = 1.0;
399:     shell->dshift = NULL;
400:     shell->left   = NULL;
401:     shell->right  = NULL;
402:     if (shell->mult) {
403:       Y->ops->mult = shell->mult;
404:       shell->mult  = NULL;
405:     }
406:     if (shell->multtranspose) {
407:       Y->ops->multtranspose = shell->multtranspose;
408:       shell->multtranspose  = NULL;
409:     }
410:     if (shell->getdiagonal) {
411:       Y->ops->getdiagonal = shell->getdiagonal;
412:       shell->getdiagonal  = NULL;
413:     }
414:     shell->usingscaled = PETSC_FALSE;
415:   }
416:   return(0);
417: }

419: extern PetscErrorCode MatConvert_Shell(Mat, MatType,MatReuse,Mat*);

421: static struct _MatOps MatOps_Values = {0,
422:                                        0,
423:                                        0,
424:                                        0,
425:                                 /* 4*/ 0,
426:                                        0,
427:                                        0,
428:                                        0,
429:                                        0,
430:                                        0,
431:                                 /*10*/ 0,
432:                                        0,
433:                                        0,
434:                                        0,
435:                                        0,
436:                                 /*15*/ 0,
437:                                        0,
438:                                        0,
439:                                        MatDiagonalScale_Shell,
440:                                        0,
441:                                 /*20*/ 0,
442:                                        MatAssemblyEnd_Shell,
443:                                        0,
444:                                        0,
445:                                 /*24*/ 0,
446:                                        0,
447:                                        0,
448:                                        0,
449:                                        0,
450:                                 /*29*/ 0,
451:                                        0,
452:                                        0,
453:                                        0,
454:                                        0,
455:                                 /*34*/ 0,
456:                                        0,
457:                                        0,
458:                                        0,
459:                                        0,
460:                                 /*39*/ 0,
461:                                        0,
462:                                        0,
463:                                        0,
464:                                        0,
465:                                 /*44*/ 0,
466:                                        MatScale_Shell,
467:                                        MatShift_Shell,
468:                                        0,
469:                                        0,
470:                                 /*49*/ 0,
471:                                        0,
472:                                        0,
473:                                        0,
474:                                        0,
475:                                 /*54*/ 0,
476:                                        0,
477:                                        0,
478:                                        0,
479:                                        0,
480:                                 /*59*/ 0,
481:                                        MatDestroy_Shell,
482:                                        0,
483:                                        0,
484:                                        0,
485:                                 /*64*/ 0,
486:                                        0,
487:                                        0,
488:                                        0,
489:                                        0,
490:                                 /*69*/ 0,
491:                                        0,
492:                                        MatConvert_Shell,
493:                                        0,
494:                                        0,
495:                                 /*74*/ 0,
496:                                        0,
497:                                        0,
498:                                        0,
499:                                        0,
500:                                 /*79*/ 0,
501:                                        0,
502:                                        0,
503:                                        0,
504:                                        0,
505:                                 /*84*/ 0,
506:                                        0,
507:                                        0,
508:                                        0,
509:                                        0,
510:                                 /*89*/ 0,
511:                                        0,
512:                                        0,
513:                                        0,
514:                                        0,
515:                                 /*94*/ 0,
516:                                        0,
517:                                        0,
518:                                        0,
519:                                        0,
520:                                 /*99*/ 0,
521:                                        0,
522:                                        0,
523:                                        0,
524:                                        0,
525:                                /*104*/ 0,
526:                                        0,
527:                                        0,
528:                                        0,
529:                                        0,
530:                                /*109*/ 0,
531:                                        0,
532:                                        0,
533:                                        0,
534:                                        0,
535:                                /*114*/ 0,
536:                                        0,
537:                                        0,
538:                                        0,
539:                                        0,
540:                                /*119*/ 0,
541:                                        0,
542:                                        0,
543:                                        0,
544:                                        0,
545:                                /*124*/ 0,
546:                                        0,
547:                                        0,
548:                                        0,
549:                                        0,
550:                                /*129*/ 0,
551:                                        0,
552:                                        0,
553:                                        0,
554:                                        0,
555:                                /*134*/ 0,
556:                                        0,
557:                                        0,
558:                                        0,
559:                                        0,
560:                                /*139*/ 0,
561:                                        0,
562:                                        0
563: };

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

568:   Level: advanced

570: .seealso: MatCreateShell
571: M*/

575: PETSC_EXTERN PetscErrorCode MatCreate_Shell(Mat A)
576: {
577:   Mat_Shell      *b;

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

583:   PetscNewLog(A,&b);
584:   A->data = (void*)b;

586:   PetscLayoutSetUp(A->rmap);
587:   PetscLayoutSetUp(A->cmap);

589:   b->ctx           = 0;
590:   b->vshift        = 0.0;
591:   b->vscale        = 1.0;
592:   b->mult          = 0;
593:   b->multtranspose = 0;
594:   b->getdiagonal   = 0;
595:   A->assembled     = PETSC_TRUE;
596:   A->preallocated  = PETSC_FALSE;

598:   PetscObjectChangeTypeName((PetscObject)A,MATSHELL);
599:   return(0);
600: }

604: /*@C
605:    MatCreateShell - Creates a new matrix class for use with a user-defined
606:    private data storage format.

608:   Collective on MPI_Comm

610:    Input Parameters:
611: +  comm - MPI communicator
612: .  m - number of local rows (must be given)
613: .  n - number of local columns (must be given)
614: .  M - number of global rows (may be PETSC_DETERMINE)
615: .  N - number of global columns (may be PETSC_DETERMINE)
616: -  ctx - pointer to data needed by the shell matrix routines

618:    Output Parameter:
619: .  A - the matrix

621:    Level: advanced

623:   Usage:
624: $    extern int mult(Mat,Vec,Vec);
625: $    MatCreateShell(comm,m,n,M,N,ctx,&mat);
626: $    MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
627: $    [ Use matrix for operations that have been set ]
628: $    MatDestroy(mat);

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

635:    Fortran Notes: The context can only be an integer or a PetscObject
636:       unfortunately it cannot be a Fortran array or derived type.

638:    PETSc requires that matrices and vectors being used for certain
639:    operations are partitioned accordingly.  For example, when
640:    creating a shell matrix, A, that supports parallel matrix-vector
641:    products using MatMult(A,x,y) the user should set the number
642:    of local matrix rows to be the number of local elements of the
643:    corresponding result vector, y. Note that this is information is
644:    required for use of the matrix interface routines, even though
645:    the shell matrix may not actually be physically partitioned.
646:    For example,

648: $
649: $     Vec x, y
650: $     extern int mult(Mat,Vec,Vec);
651: $     Mat A
652: $
653: $     VecCreateMPI(comm,PETSC_DECIDE,M,&y);
654: $     VecCreateMPI(comm,PETSC_DECIDE,N,&x);
655: $     VecGetLocalSize(y,&m);
656: $     VecGetLocalSize(x,&n);
657: $     MatCreateShell(comm,m,n,M,N,ctx,&A);
658: $     MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
659: $     MatMult(A,x,y);
660: $     MatDestroy(A);
661: $     VecDestroy(y); VecDestroy(x);
662: $

664: .keywords: matrix, shell, create

666: .seealso: MatShellSetOperation(), MatHasOperation(), MatShellGetContext(), MatShellSetContext()
667: @*/
668: PetscErrorCode  MatCreateShell(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,void *ctx,Mat *A)
669: {

673:   MatCreate(comm,A);
674:   MatSetSizes(*A,m,n,M,N);
675:   MatSetType(*A,MATSHELL);
676:   MatShellSetContext(*A,ctx);
677:   MatSetUp(*A);
678:   return(0);
679: }

683: /*@
684:     MatShellSetContext - sets the context for a shell matrix

686:    Logically Collective on Mat

688:     Input Parameters:
689: +   mat - the shell matrix
690: -   ctx - the context

692:    Level: advanced

694:    Fortran Notes: The context can only be an integer or a PetscObject
695:       unfortunately it cannot be a Fortran array or derived type.

697: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation()
698: @*/
699: PetscErrorCode  MatShellSetContext(Mat mat,void *ctx)
700: {
701:   Mat_Shell      *shell = (Mat_Shell*)mat->data;
703:   PetscBool      flg;

707:   PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
708:   if (flg) {
709:     shell->ctx = ctx;
710:   } else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot attach context to non-shell matrix");
711:   return(0);
712: }

716: /*@C
717:     MatShellSetOperation - Allows user to set a matrix operation for
718:                            a shell matrix.

720:    Logically Collective on Mat

722:     Input Parameters:
723: +   mat - the shell matrix
724: .   op - the name of the operation
725: -   f - the function that provides the operation.

727:    Level: advanced

729:     Usage:
730: $      extern PetscErrorCode usermult(Mat,Vec,Vec);
731: $      MatCreateShell(comm,m,n,M,N,ctx,&A);
732: $      MatShellSetOperation(A,MATOP_MULT,(void(*)(void))usermult);

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

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

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

749:     Within each user-defined routine, the user should call
750:     MatShellGetContext() to obtain the user-defined context that was
751:     set by MatCreateShell().

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

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

758: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext()
759: @*/
760: PetscErrorCode  MatShellSetOperation(Mat mat,MatOperation op,void (*f)(void))
761: {
763:   PetscBool      flg;

767:   switch (op) {
768:   case MATOP_DESTROY:
769:     PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
770:     if (flg) {
771:       Mat_Shell *shell = (Mat_Shell*)mat->data;
772:       shell->destroy = (PetscErrorCode (*)(Mat))f;
773:     } else mat->ops->destroy = (PetscErrorCode (*)(Mat))f;
774:     break;
775:   case MATOP_VIEW:
776:     mat->ops->view = (PetscErrorCode (*)(Mat,PetscViewer))f;
777:     break;
778:   case MATOP_MULT:
779:     mat->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f;
780:     if (!mat->ops->multadd) mat->ops->multadd = MatMultAdd_Shell;
781:     break;
782:   case MATOP_MULT_TRANSPOSE:
783:     mat->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f;
784:     if (!mat->ops->multtransposeadd) mat->ops->multtransposeadd = MatMultTransposeAdd_Shell;
785:     break;
786:   default:
787:     (((void(**)(void))mat->ops)[op]) = f;
788:   }
789:   return(0);
790: }

794: /*@C
795:     MatShellGetOperation - Gets a matrix function for a shell matrix.

797:     Not Collective

799:     Input Parameters:
800: +   mat - the shell matrix
801: -   op - the name of the operation

803:     Output Parameter:
804: .   f - the function that provides the operation.

806:     Level: advanced

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

814:     All user-provided functions have the same calling
815:     sequence as the usual matrix interface routines, since they
816:     are intended to be accessed via the usual matrix interface
817:     routines, e.g.,
818: $       MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec)

820:     Within each user-defined routine, the user should call
821:     MatShellGetContext() to obtain the user-defined context that was
822:     set by MatCreateShell().

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

826: .seealso: MatCreateShell(), MatShellGetContext(), MatShellSetOperation(), MatShellSetContext()
827: @*/
828: PetscErrorCode  MatShellGetOperation(Mat mat,MatOperation op,void(**f)(void))
829: {
831:   PetscBool      flg;

835:   if (op == MATOP_DESTROY) {
836:     PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
837:     if (flg) {
838:       Mat_Shell *shell = (Mat_Shell*)mat->data;
839:       *f = (void (*)(void))shell->destroy;
840:     } else {
841:       *f = (void (*)(void))mat->ops->destroy;
842:     }
843:   } else if (op == MATOP_VIEW) {
844:     *f = (void (*)(void))mat->ops->view;
845:   } else {
846:     *f = (((void (**)(void))mat->ops)[op]);
847:   }
848:   return(0);
849: }