Actual source code: shell.c

petsc-3.7.3 2016-08-01
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:    Fortran Notes: To use this from Fortran you must write a Fortran interface definition for this
189:     function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

590:   Level: advanced

592: .seealso: MatCreateShell
593: M*/

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

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

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

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

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

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

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

630:   Collective on MPI_Comm

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

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

643:    Level: advanced

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

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

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

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

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

687: .keywords: matrix, shell, create

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

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

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

709:    Logically Collective on Mat

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

715:    Level: advanced

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

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

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

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

743:    Logically Collective on Mat

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

750:    Level: advanced

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

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

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

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

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

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

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

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

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

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

820:     Not Collective

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

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

829:     Level: advanced

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

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

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

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

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

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