Actual source code: shell.c

petsc-master 2016-09-25
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>

 10: struct _MatShellOps {
 11:   /*   0 */
 12:   PetscErrorCode (*mult)(Mat,Vec,Vec);
 13:   /*   5 */
 14:   PetscErrorCode (*multtranspose)(Mat,Vec,Vec);
 15:   /*  10 */
 16:   /*  15 */
 17:   PetscErrorCode (*getdiagonal)(Mat,Vec);
 18:   /*  20 */
 19:   /*  24 */
 20:   /*  29 */
 21:   /*  34 */
 22:   /*  39 */
 23:   PetscErrorCode (*copy)(Mat,Mat,MatStructure);
 24:   /*  44 */
 25:   PetscErrorCode (*diagonalset)(Mat,Vec,InsertMode);
 26:   /*  49 */
 27:   /*  54 */
 28:   /*  59 */
 29:   PetscErrorCode (*destroy)(Mat);
 30:   /*  64 */
 31:   /*  69 */
 32:   /*  74 */
 33:   /*  79 */
 34:   /*  84 */
 35:   /*  89 */
 36:   /*  94 */
 37:   /*  99 */
 38:   /* 104 */
 39:   /* 109 */
 40:   /* 114 */
 41:   /* 119 */
 42:   /* 124 */
 43:   /* 129 */
 44:   /* 134 */
 45:   /* 139 */
 46:   /* 144 */
 47: };

 49: typedef struct {
 50:   struct _MatShellOps ops[1];

 52:   PetscScalar vscale,vshift;
 53:   Vec         dshift;
 54:   Vec         left,right;
 55:   Vec         dshift_owned,left_owned,right_owned;
 56:   Vec         left_work,right_work;
 57:   Vec         left_add_work,right_add_work;
 58:   PetscBool   usingscaled;
 59:   void        *ctx;
 60: } Mat_Shell;
 61: /*
 62:  The most general expression for the matrix is

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

 66:  where
 67:  A is the matrix defined by the user's function
 68:  a is a scalar multiple
 69:  L is left scaling
 70:  R is right scaling
 71:  B is a diagonal shift defined by
 72:    diag(dshift) if the vector dshift is non-NULL
 73:    vscale*identity otherwise

 75:  The following identities apply:

 77:  Scale by c:
 78:   c [L (a A + B) R] = L [(a c) A + c B] R

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

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

 85:  In the data structure:

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

 92: static PetscErrorCode MatMult_Shell(Mat,Vec,Vec);
 93: static PetscErrorCode MatMultTranspose_Shell(Mat,Vec,Vec);
 94: static PetscErrorCode MatGetDiagonal_Shell(Mat,Vec);
 95: static PetscErrorCode MatCopy_Shell(Mat,Mat,MatStructure);

 99: static PetscErrorCode MatShellUseScaledMethods(Mat Y)
100: {
101:   Mat_Shell *shell = (Mat_Shell*)Y->data;

104:   if (shell->usingscaled) return(0);
105:   shell->ops->mult  = Y->ops->mult;
106:   Y->ops->mult = MatMult_Shell;
107:   if (Y->ops->multtranspose) {
108:     shell->ops->multtranspose  = Y->ops->multtranspose;
109:     Y->ops->multtranspose = MatMultTranspose_Shell;
110:   }
111:   if (Y->ops->getdiagonal) {
112:     shell->ops->getdiagonal  = Y->ops->getdiagonal;
113:     Y->ops->getdiagonal = MatGetDiagonal_Shell;
114:   }
115:   if (Y->ops->copy) {
116:     shell->ops->copy  = Y->ops->copy;
117:     Y->ops->copy = MatCopy_Shell;
118:   }
119:   shell->usingscaled = PETSC_TRUE;
120:   return(0);
121: }

125: static PetscErrorCode MatShellPreScaleLeft(Mat A,Vec x,Vec *xx)
126: {
127:   Mat_Shell      *shell = (Mat_Shell*)A->data;

131:   *xx = NULL;
132:   if (!shell->left) {
133:     *xx = x;
134:   } else {
135:     if (!shell->left_work) {VecDuplicate(shell->left,&shell->left_work);}
136:     VecPointwiseMult(shell->left_work,x,shell->left);
137:     *xx  = shell->left_work;
138:   }
139:   return(0);
140: }

144: static PetscErrorCode MatShellPreScaleRight(Mat A,Vec x,Vec *xx)
145: {
146:   Mat_Shell      *shell = (Mat_Shell*)A->data;

150:   *xx = NULL;
151:   if (!shell->right) {
152:     *xx = x;
153:   } else {
154:     if (!shell->right_work) {VecDuplicate(shell->right,&shell->right_work);}
155:     VecPointwiseMult(shell->right_work,x,shell->right);
156:     *xx  = shell->right_work;
157:   }
158:   return(0);
159: }

163: static PetscErrorCode MatShellPostScaleLeft(Mat A,Vec x)
164: {
165:   Mat_Shell      *shell = (Mat_Shell*)A->data;

169:   if (shell->left) {VecPointwiseMult(x,x,shell->left);}
170:   return(0);
171: }

175: static PetscErrorCode MatShellPostScaleRight(Mat A,Vec x)
176: {
177:   Mat_Shell      *shell = (Mat_Shell*)A->data;

181:   if (shell->right) {VecPointwiseMult(x,x,shell->right);}
182:   return(0);
183: }

187: static PetscErrorCode MatShellShiftAndScale(Mat A,Vec X,Vec Y)
188: {
189:   Mat_Shell      *shell = (Mat_Shell*)A->data;

193:   if (shell->dshift) {          /* get arrays because there is no VecPointwiseMultAdd() */
194:     PetscInt          i,m;
195:     const PetscScalar *x,*d;
196:     PetscScalar       *y;
197:     VecGetLocalSize(X,&m);
198:     VecGetArrayRead(shell->dshift,&d);
199:     VecGetArrayRead(X,&x);
200:     VecGetArray(Y,&y);
201:     for (i=0; i<m; i++) y[i] = shell->vscale*y[i] + d[i]*x[i];
202:     VecRestoreArrayRead(shell->dshift,&d);
203:     VecRestoreArrayRead(X,&x);
204:     VecRestoreArray(Y,&y);
205:   } else if (PetscAbsScalar(shell->vshift) != 0) {
206:     VecAXPBY(Y,shell->vshift,shell->vscale,X);
207:   } else if (shell->vscale != 1.0) {
208:     VecScale(Y,shell->vscale);
209:   }
210:   return(0);
211: }

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

218:     Not Collective

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

223:     Output Parameter:
224: .   ctx - the user provided context

226:     Level: advanced

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

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

233: .seealso: MatCreateShell(), MatShellSetOperation(), MatShellSetContext()
234: @*/
235: PetscErrorCode  MatShellGetContext(Mat mat,void *ctx)
236: {
238:   PetscBool      flg;

243:   PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
244:   if (flg) *(void**)ctx = ((Mat_Shell*)(mat->data))->ctx;
245:   else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot get context from non-shell matrix");
246:   return(0);
247: }

251: PetscErrorCode MatDestroy_Shell(Mat mat)
252: {
254:   Mat_Shell      *shell = (Mat_Shell*)mat->data;

257:   if (shell->ops->destroy) {
258:     (*shell->ops->destroy)(mat);
259:   }
260:   VecDestroy(&shell->left_owned);
261:   VecDestroy(&shell->right_owned);
262:   VecDestroy(&shell->dshift_owned);
263:   VecDestroy(&shell->left_work);
264:   VecDestroy(&shell->right_work);
265:   VecDestroy(&shell->left_add_work);
266:   VecDestroy(&shell->right_add_work);
267:   PetscFree(mat->data);
268:   return(0);
269: }

273: PetscErrorCode MatCopy_Shell(Mat A,Mat B,MatStructure str)
274: {
275:   Mat_Shell       *shellA = (Mat_Shell*)A->data,*shellB = (Mat_Shell*)B->data;
276:   PetscBool       matflg,shellflg;
277:   PetscErrorCode  ierr;

280:   PetscObjectTypeCompare((PetscObject)B,MATSHELL,&matflg);
281:   if(!matflg) { SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_NOTSAMETYPE,"Output matrix must be a MATSHELL"); }

283:   MatShellUseScaledMethods(B);
284:   PetscMemcmp(A->ops,B->ops,sizeof(struct _MatOps),&matflg);
285:   PetscMemcmp(shellA->ops,shellB->ops,sizeof(struct _MatShellOps),&shellflg);
286:   if (!matflg || !shellflg) { SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Matrices cannot be copied because they have different operations"); }

288:   (*shellA->ops->copy)(A,B,str);
289:   shellB->vscale = shellA->vscale;
290:   shellB->vshift = shellA->vshift;
291:   if (shellA->dshift_owned) {
292:     if (!shellB->dshift_owned) {
293:       VecDuplicate(shellA->dshift_owned,&shellB->dshift_owned);
294:     }
295:     VecCopy(shellA->dshift_owned,shellB->dshift_owned);
296:     shellB->dshift = shellB->dshift_owned;
297:   } else {
298:     shellB->dshift = NULL;
299:   }
300:   if (shellA->left_owned) {
301:     if (!shellB->left_owned) {
302:       VecDuplicate(shellA->left_owned,&shellB->left_owned);
303:     }
304:     VecCopy(shellA->left_owned,shellB->left_owned);
305:     shellB->left = shellB->left_owned;
306:   } else {
307:     shellB->left = NULL;
308:   }
309:   if (shellA->right_owned) {
310:     if (!shellB->right_owned) {
311:       VecDuplicate(shellA->right_owned,&shellB->right_owned);
312:     }
313:     VecCopy(shellA->right_owned,shellB->right_owned);
314:     shellB->right = shellB->right_owned;
315:   } else {
316:     shellB->right = NULL;
317:   }
318:   return(0);
319: }

323: PetscErrorCode MatMult_Shell(Mat A,Vec x,Vec y)
324: {
325:   Mat_Shell        *shell = (Mat_Shell*)A->data;
326:   PetscErrorCode   ierr;
327:   Vec              xx;
328:   PetscObjectState instate,outstate;

331:   MatShellPreScaleRight(A,x,&xx);
332:   PetscObjectStateGet((PetscObject)y, &instate);
333:   (*shell->ops->mult)(A,xx,y);
334:   PetscObjectStateGet((PetscObject)y, &outstate);
335:   if (instate == outstate) {
336:     /* increase the state of the output vector since the user did not update its state themself as should have been done */
337:     PetscObjectStateIncrease((PetscObject)y);
338:   }
339:   MatShellShiftAndScale(A,xx,y);
340:   MatShellPostScaleLeft(A,y);
341:   return(0);
342: }

346: PetscErrorCode MatMultAdd_Shell(Mat A,Vec x,Vec y,Vec z)
347: {
348:   Mat_Shell      *shell = (Mat_Shell*)A->data;

352:   if (y == z) {
353:     if (!shell->right_add_work) {VecDuplicate(z,&shell->right_add_work);}
354:     MatMult(A,x,shell->right_add_work);
355:     VecAXPY(z,1.0,shell->right_add_work);
356:   } else {
357:     MatMult(A,x,z);
358:     VecAXPY(z,1.0,y);
359:   }
360:   return(0);
361: }

365: PetscErrorCode MatMultTranspose_Shell(Mat A,Vec x,Vec y)
366: {
367:   Mat_Shell        *shell = (Mat_Shell*)A->data;
368:   PetscErrorCode   ierr;
369:   Vec              xx;
370:   PetscObjectState instate,outstate;

373:   MatShellPreScaleLeft(A,x,&xx);
374:   PetscObjectStateGet((PetscObject)y, &instate);
375:   (*shell->ops->multtranspose)(A,xx,y);
376:   PetscObjectStateGet((PetscObject)y, &outstate);
377:   if (instate == outstate) {
378:     /* increase the state of the output vector since the user did not update its state themself as should have been done */
379:     PetscObjectStateIncrease((PetscObject)y);
380:   }
381:   MatShellShiftAndScale(A,xx,y);
382:   MatShellPostScaleRight(A,y);
383:   return(0);
384: }

388: PetscErrorCode MatMultTransposeAdd_Shell(Mat A,Vec x,Vec y,Vec z)
389: {
390:   Mat_Shell      *shell = (Mat_Shell*)A->data;

394:   if (y == z) {
395:     if (!shell->left_add_work) {VecDuplicate(z,&shell->left_add_work);}
396:     MatMultTranspose(A,x,shell->left_add_work);
397:     VecWAXPY(z,1.0,shell->left_add_work,y);
398:   } else {
399:     MatMultTranspose(A,x,z);
400:     VecAXPY(z,1.0,y);
401:   }
402:   return(0);
403: }

407: PetscErrorCode MatGetDiagonal_Shell(Mat A,Vec v)
408: {
409:   Mat_Shell      *shell = (Mat_Shell*)A->data;

413:   (*shell->ops->getdiagonal)(A,v);
414:   VecScale(v,shell->vscale);
415:   if (shell->dshift) {
416:     VecPointwiseMult(v,v,shell->dshift);
417:   } else {
418:     VecShift(v,shell->vshift);
419:   }
420:   if (shell->left)  {VecPointwiseMult(v,v,shell->left);}
421:   if (shell->right) {VecPointwiseMult(v,v,shell->right);}
422:   return(0);
423: }

427: PetscErrorCode MatShift_Shell(Mat Y,PetscScalar a)
428: {
429:   Mat_Shell      *shell = (Mat_Shell*)Y->data;

433:   if (shell->left || shell->right || shell->dshift) {
434:     if (!shell->dshift) {
435:       if (!shell->dshift_owned) {VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift_owned);}
436:       shell->dshift = shell->dshift_owned;
437:       VecSet(shell->dshift,shell->vshift+a);
438:     } else {VecScale(shell->dshift,a);}
439:     if (shell->left)  {VecPointwiseDivide(shell->dshift,shell->dshift,shell->left);}
440:     if (shell->right) {VecPointwiseDivide(shell->dshift,shell->dshift,shell->right);}
441:   } else shell->vshift += a;
442:   MatShellUseScaledMethods(Y);
443:   return(0);
444: }

448: PetscErrorCode MatDiagonalSet_Shell(Mat A,Vec D,InsertMode ins)
449: {
450:   Mat_Shell      *shell = (Mat_Shell*)A->data;

454:   if (shell->vscale != 1.0 || shell->left || shell->right) {
455:     SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE, "Operation not supported on scaled matrix");
456:   }

458:   if (shell->ops->diagonalset) {
459:     (*shell->ops->diagonalset)(A,D,ins);
460:   }
461:   shell->vshift = 0.0;
462:   if (shell->dshift) { VecZeroEntries(shell->dshift); }
463:   return(0);
464: }

468: PetscErrorCode MatScale_Shell(Mat Y,PetscScalar a)
469: {
470:   Mat_Shell      *shell = (Mat_Shell*)Y->data;

474:   shell->vscale *= a;
475:   if (shell->dshift) {
476:     VecScale(shell->dshift,a);
477:   } else shell->vshift *= a;
478:   MatShellUseScaledMethods(Y);
479:   return(0);
480: }

484: static PetscErrorCode MatDiagonalScale_Shell(Mat Y,Vec left,Vec right)
485: {
486:   Mat_Shell      *shell = (Mat_Shell*)Y->data;

490:   if (left) {
491:     if (!shell->left_owned) {VecDuplicate(left,&shell->left_owned);}
492:     if (shell->left) {
493:       VecPointwiseMult(shell->left,shell->left,left);
494:     } else {
495:       shell->left = shell->left_owned;
496:       VecCopy(left,shell->left);
497:     }
498:   }
499:   if (right) {
500:     if (!shell->right_owned) {VecDuplicate(right,&shell->right_owned);}
501:     if (shell->right) {
502:       VecPointwiseMult(shell->right,shell->right,right);
503:     } else {
504:       shell->right = shell->right_owned;
505:       VecCopy(right,shell->right);
506:     }
507:   }
508:   MatShellUseScaledMethods(Y);
509:   return(0);
510: }

514: PetscErrorCode MatAssemblyEnd_Shell(Mat Y,MatAssemblyType t)
515: {
516:   Mat_Shell *shell = (Mat_Shell*)Y->data;

519:   if (t == MAT_FINAL_ASSEMBLY) {
520:     shell->vshift = 0.0;
521:     shell->vscale = 1.0;
522:     shell->dshift = NULL;
523:     shell->left   = NULL;
524:     shell->right  = NULL;
525:     if (shell->ops->mult) {
526:       Y->ops->mult = shell->ops->mult;
527:       shell->ops->mult  = NULL;
528:     }
529:     if (shell->ops->multtranspose) {
530:       Y->ops->multtranspose = shell->ops->multtranspose;
531:       shell->ops->multtranspose = NULL;
532:     }
533:     if (shell->ops->getdiagonal) {
534:       Y->ops->getdiagonal = shell->ops->getdiagonal;
535:       shell->ops->getdiagonal = NULL;
536:     }
537:     if (shell->ops->copy) {
538:       Y->ops->copy = shell->ops->copy;
539:       shell->ops->copy = NULL;
540:     }
541:     shell->usingscaled = PETSC_FALSE;
542:   }
543:   return(0);
544: }

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

550: static PetscErrorCode MatMissingDiagonal_Shell(Mat A,PetscBool  *missing,PetscInt *d)
551: {
553:   *missing = PETSC_FALSE;
554:   return(0);
555: }

557: static struct _MatOps MatOps_Values = {0,
558:                                        0,
559:                                        0,
560:                                        0,
561:                                 /* 4*/ 0,
562:                                        0,
563:                                        0,
564:                                        0,
565:                                        0,
566:                                        0,
567:                                 /*10*/ 0,
568:                                        0,
569:                                        0,
570:                                        0,
571:                                        0,
572:                                 /*15*/ 0,
573:                                        0,
574:                                        0,
575:                                        MatDiagonalScale_Shell,
576:                                        0,
577:                                 /*20*/ 0,
578:                                        MatAssemblyEnd_Shell,
579:                                        0,
580:                                        0,
581:                                 /*24*/ 0,
582:                                        0,
583:                                        0,
584:                                        0,
585:                                        0,
586:                                 /*29*/ 0,
587:                                        0,
588:                                        0,
589:                                        0,
590:                                        0,
591:                                 /*34*/ 0,
592:                                        0,
593:                                        0,
594:                                        0,
595:                                        0,
596:                                 /*39*/ 0,
597:                                        0,
598:                                        0,
599:                                        0,
600:                                        0,
601:                                 /*44*/ 0,
602:                                        MatScale_Shell,
603:                                        MatShift_Shell,
604:                                        MatDiagonalSet_Shell,
605:                                        0,
606:                                 /*49*/ 0,
607:                                        0,
608:                                        0,
609:                                        0,
610:                                        0,
611:                                 /*54*/ 0,
612:                                        0,
613:                                        0,
614:                                        0,
615:                                        0,
616:                                 /*59*/ 0,
617:                                        MatDestroy_Shell,
618:                                        0,
619:                                        0,
620:                                        0,
621:                                 /*64*/ 0,
622:                                        0,
623:                                        0,
624:                                        0,
625:                                        0,
626:                                 /*69*/ 0,
627:                                        0,
628:                                        MatConvert_Shell,
629:                                        0,
630:                                        0,
631:                                 /*74*/ 0,
632:                                        0,
633:                                        0,
634:                                        0,
635:                                        0,
636:                                 /*79*/ 0,
637:                                        0,
638:                                        0,
639:                                        0,
640:                                        0,
641:                                 /*84*/ 0,
642:                                        0,
643:                                        0,
644:                                        0,
645:                                        0,
646:                                 /*89*/ 0,
647:                                        0,
648:                                        0,
649:                                        0,
650:                                        0,
651:                                 /*94*/ 0,
652:                                        0,
653:                                        0,
654:                                        0,
655:                                        0,
656:                                 /*99*/ 0,
657:                                        0,
658:                                        0,
659:                                        0,
660:                                        0,
661:                                /*104*/ 0,
662:                                        0,
663:                                        0,
664:                                        0,
665:                                        0,
666:                                /*109*/ 0,
667:                                        0,
668:                                        0,
669:                                        0,
670:                                        MatMissingDiagonal_Shell,
671:                                /*114*/ 0,
672:                                        0,
673:                                        0,
674:                                        0,
675:                                        0,
676:                                /*119*/ 0,
677:                                        0,
678:                                        0,
679:                                        0,
680:                                        0,
681:                                /*124*/ 0,
682:                                        0,
683:                                        0,
684:                                        0,
685:                                        0,
686:                                /*129*/ 0,
687:                                        0,
688:                                        0,
689:                                        0,
690:                                        0,
691:                                /*134*/ 0,
692:                                        0,
693:                                        0,
694:                                        0,
695:                                        0,
696:                                /*139*/ 0,
697:                                        0,
698:                                        0
699: };

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

704:   Level: advanced

706: .seealso: MatCreateShell
707: M*/

711: PETSC_EXTERN PetscErrorCode MatCreate_Shell(Mat A)
712: {
713:   Mat_Shell      *b;

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

719:   PetscNewLog(A,&b);
720:   A->data = (void*)b;

722:   PetscLayoutSetUp(A->rmap);
723:   PetscLayoutSetUp(A->cmap);

725:   b->ctx           = 0;
726:   b->vshift        = 0.0;
727:   b->vscale        = 1.0;
728:   A->assembled     = PETSC_TRUE;
729:   A->preallocated  = PETSC_FALSE;

731:   PetscObjectChangeTypeName((PetscObject)A,MATSHELL);
732:   return(0);
733: }

737: /*@C
738:    MatCreateShell - Creates a new matrix class for use with a user-defined
739:    private data storage format.

741:   Collective on MPI_Comm

743:    Input Parameters:
744: +  comm - MPI communicator
745: .  m - number of local rows (must be given)
746: .  n - number of local columns (must be given)
747: .  M - number of global rows (may be PETSC_DETERMINE)
748: .  N - number of global columns (may be PETSC_DETERMINE)
749: -  ctx - pointer to data needed by the shell matrix routines

751:    Output Parameter:
752: .  A - the matrix

754:    Level: advanced

756:   Usage:
757: $    extern int mult(Mat,Vec,Vec);
758: $    MatCreateShell(comm,m,n,M,N,ctx,&mat);
759: $    MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
760: $    [ Use matrix for operations that have been set ]
761: $    MatDestroy(mat);

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

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

772:    PETSc requires that matrices and vectors being used for certain
773:    operations are partitioned accordingly.  For example, when
774:    creating a shell matrix, A, that supports parallel matrix-vector
775:    products using MatMult(A,x,y) the user should set the number
776:    of local matrix rows to be the number of local elements of the
777:    corresponding result vector, y. Note that this is information is
778:    required for use of the matrix interface routines, even though
779:    the shell matrix may not actually be physically partitioned.
780:    For example,

782: $
783: $     Vec x, y
784: $     extern int mult(Mat,Vec,Vec);
785: $     Mat A
786: $
787: $     VecCreateMPI(comm,PETSC_DECIDE,M,&y);
788: $     VecCreateMPI(comm,PETSC_DECIDE,N,&x);
789: $     VecGetLocalSize(y,&m);
790: $     VecGetLocalSize(x,&n);
791: $     MatCreateShell(comm,m,n,M,N,ctx,&A);
792: $     MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
793: $     MatMult(A,x,y);
794: $     MatDestroy(A);
795: $     VecDestroy(y); VecDestroy(x);
796: $

798: .keywords: matrix, shell, create

800: .seealso: MatShellSetOperation(), MatHasOperation(), MatShellGetContext(), MatShellSetContext()
801: @*/
802: PetscErrorCode  MatCreateShell(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,void *ctx,Mat *A)
803: {

807:   MatCreate(comm,A);
808:   MatSetSizes(*A,m,n,M,N);
809:   MatSetType(*A,MATSHELL);
810:   MatShellSetContext(*A,ctx);
811:   MatSetUp(*A);
812:   return(0);
813: }

817: /*@
818:     MatShellSetContext - sets the context for a shell matrix

820:    Logically Collective on Mat

822:     Input Parameters:
823: +   mat - the shell matrix
824: -   ctx - the context

826:    Level: advanced

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

831: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation()
832: @*/
833: PetscErrorCode  MatShellSetContext(Mat mat,void *ctx)
834: {
835:   Mat_Shell      *shell = (Mat_Shell*)mat->data;
837:   PetscBool      flg;

841:   PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
842:   if (flg) {
843:     shell->ctx = ctx;
844:   } else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot attach context to non-shell matrix");
845:   return(0);
846: }

850: /*@C
851:     MatShellSetOperation - Allows user to set a matrix operation for
852:                            a shell matrix.

854:    Logically Collective on Mat

856:     Input Parameters:
857: +   mat - the shell matrix
858: .   op - the name of the operation
859: -   f - the function that provides the operation.

861:    Level: advanced

863:     Usage:
864: $      extern PetscErrorCode usermult(Mat,Vec,Vec);
865: $      MatCreateShell(comm,m,n,M,N,ctx,&A);
866: $      MatShellSetOperation(A,MATOP_MULT,(void(*)(void))usermult);

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

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

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

883:     Within each user-defined routine, the user should call
884:     MatShellGetContext() to obtain the user-defined context that was
885:     set by MatCreateShell().

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

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

892: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext()
893: @*/
894: PetscErrorCode  MatShellSetOperation(Mat mat,MatOperation op,void (*f)(void))
895: {
897:   PetscBool      flg;

901:   switch (op) {
902:   case MATOP_DESTROY:
903:     PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
904:     if (flg) {
905:       Mat_Shell *shell = (Mat_Shell*)mat->data;
906:       shell->ops->destroy = (PetscErrorCode (*)(Mat))f;
907:     } else mat->ops->destroy = (PetscErrorCode (*)(Mat))f;
908:     break;
909:   case MATOP_DIAGONAL_SET:
910:     PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
911:     if (flg) {
912:       Mat_Shell *shell = (Mat_Shell*)mat->data;
913:       shell->ops->diagonalset = (PetscErrorCode (*)(Mat,Vec,InsertMode))f;
914:     } else {
915:       mat->ops->diagonalset = (PetscErrorCode (*)(Mat,Vec,InsertMode))f;
916:     }
917:     break;
918:   case MATOP_VIEW:
919:     mat->ops->view = (PetscErrorCode (*)(Mat,PetscViewer))f;
920:     break;
921:   case MATOP_MULT:
922:     mat->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f;
923:     if (!mat->ops->multadd) mat->ops->multadd = MatMultAdd_Shell;
924:     break;
925:   case MATOP_MULT_TRANSPOSE:
926:     mat->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f;
927:     if (!mat->ops->multtransposeadd) mat->ops->multtransposeadd = MatMultTransposeAdd_Shell;
928:     break;
929:   default:
930:     (((void(**)(void))mat->ops)[op]) = f;
931:   }
932:   return(0);
933: }

937: /*@C
938:     MatShellGetOperation - Gets a matrix function for a shell matrix.

940:     Not Collective

942:     Input Parameters:
943: +   mat - the shell matrix
944: -   op - the name of the operation

946:     Output Parameter:
947: .   f - the function that provides the operation.

949:     Level: advanced

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

957:     All user-provided functions have the same calling
958:     sequence as the usual matrix interface routines, since they
959:     are intended to be accessed via the usual matrix interface
960:     routines, e.g.,
961: $       MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec)

963:     Within each user-defined routine, the user should call
964:     MatShellGetContext() to obtain the user-defined context that was
965:     set by MatCreateShell().

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

969: .seealso: MatCreateShell(), MatShellGetContext(), MatShellSetOperation(), MatShellSetContext()
970: @*/
971: PetscErrorCode  MatShellGetOperation(Mat mat,MatOperation op,void(**f)(void))
972: {
974:   PetscBool      flg;

978:   switch (op) {
979:   case MATOP_DESTROY:
980:     PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
981:     if (flg) {
982:       Mat_Shell *shell = (Mat_Shell*)mat->data;
983:       *f = (void (*)(void))shell->ops->destroy;
984:     } else {
985:       *f = (void (*)(void))mat->ops->destroy;
986:     }
987:     break;
988:   case MATOP_DIAGONAL_SET:
989:     PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
990:     if (flg) {
991:       Mat_Shell *shell = (Mat_Shell*)mat->data;
992:       *f = (void (*)(void))shell->ops->diagonalset;
993:     } else {
994:       *f = (void (*)(void))mat->ops->diagonalset;
995:     }
996:     break;
997:   case MATOP_VIEW:
998:     *f = (void (*)(void))mat->ops->view;
999:     break;
1000:   default:
1001:     *f = (((void (**)(void))mat->ops)[op]);
1002:   }
1003:   return(0);
1004: }