Actual source code: shell.c

petsc-master 2018-02-21
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 (*axpy)(Mat,PetscScalar,Mat,MatStructure);
 24:   PetscErrorCode (*copy)(Mat,Mat,MatStructure);
 25:   /*  44 */
 26:   PetscErrorCode (*diagonalset)(Mat,Vec,InsertMode);
 27:   /*  49 */
 28:   /*  54 */
 29:   /*  59 */
 30:   PetscErrorCode (*destroy)(Mat);
 31:   /*  64 */
 32:   /*  69 */
 33:   /*  74 */
 34:   /*  79 */
 35:   /*  84 */
 36:   /*  89 */
 37:   /*  94 */
 38:   /*  99 */
 39:   /* 104 */
 40:   /* 109 */
 41:   /* 114 */
 42:   /* 119 */
 43:   /* 124 */
 44:   /* 129 */
 45:   /* 134 */
 46:   /* 139 */
 47:   /* 144 */
 48: };

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

 53:   PetscScalar vscale,vshift;
 54:   Vec         dshift;
 55:   Vec         left,right;
 56:   Vec         left_work,right_work;
 57:   Vec         left_add_work,right_add_work;
 58:   Mat         axpy;
 59:   PetscScalar axpy_vscale;
 60:   PetscBool   managescalingshifts;                   /* The user will manage the scaling and shifts for the MATSHELL, not the default */
 61:   void        *ctx;
 62: } Mat_Shell;

 64: /*
 65:       xx = diag(left)*x
 66: */
 67: static PetscErrorCode MatShellPreScaleLeft(Mat A,Vec x,Vec *xx)
 68: {
 69:   Mat_Shell      *shell = (Mat_Shell*)A->data;

 73:   *xx = NULL;
 74:   if (!shell->left) {
 75:     *xx = x;
 76:   } else {
 77:     if (!shell->left_work) {VecDuplicate(shell->left,&shell->left_work);}
 78:     VecPointwiseMult(shell->left_work,x,shell->left);
 79:     *xx  = shell->left_work;
 80:   }
 81:   return(0);
 82: }

 84: /*
 85:      xx = diag(right)*x
 86: */
 87: static PetscErrorCode MatShellPreScaleRight(Mat A,Vec x,Vec *xx)
 88: {
 89:   Mat_Shell      *shell = (Mat_Shell*)A->data;

 93:   *xx = NULL;
 94:   if (!shell->right) {
 95:     *xx = x;
 96:   } else {
 97:     if (!shell->right_work) {VecDuplicate(shell->right,&shell->right_work);}
 98:     VecPointwiseMult(shell->right_work,x,shell->right);
 99:     *xx  = shell->right_work;
100:   }
101:   return(0);
102: }

104: /*
105:     x = diag(left)*x
106: */
107: static PetscErrorCode MatShellPostScaleLeft(Mat A,Vec x)
108: {
109:   Mat_Shell      *shell = (Mat_Shell*)A->data;

113:   if (shell->left) {VecPointwiseMult(x,x,shell->left);}
114:   return(0);
115: }

117: /*
118:     x = diag(right)*x
119: */
120: static PetscErrorCode MatShellPostScaleRight(Mat A,Vec x)
121: {
122:   Mat_Shell      *shell = (Mat_Shell*)A->data;

126:   if (shell->right) {VecPointwiseMult(x,x,shell->right);}
127:   return(0);
128: }

130: /*
131:          Y = vscale*Y + diag(dshift)*X + vshift*X

133:          On input Y already contains A*x
134: */
135: static PetscErrorCode MatShellShiftAndScale(Mat A,Vec X,Vec Y)
136: {
137:   Mat_Shell      *shell = (Mat_Shell*)A->data;

141:   if (shell->dshift) {          /* get arrays because there is no VecPointwiseMultAdd() */
142:     PetscInt          i,m;
143:     const PetscScalar *x,*d;
144:     PetscScalar       *y;
145:     VecGetLocalSize(X,&m);
146:     VecGetArrayRead(shell->dshift,&d);
147:     VecGetArrayRead(X,&x);
148:     VecGetArray(Y,&y);
149:     for (i=0; i<m; i++) y[i] = shell->vscale*y[i] + d[i]*x[i];
150:     VecRestoreArrayRead(shell->dshift,&d);
151:     VecRestoreArrayRead(X,&x);
152:     VecRestoreArray(Y,&y);
153:   } else {
154:     VecScale(Y,shell->vscale);
155:   }
156:   if (shell->vshift != 0.0) {VecAXPY(Y,shell->vshift,X);} /* if test is for non-square matrices */
157:   return(0);
158: }

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

163:     Not Collective

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

168:     Output Parameter:
169: .   ctx - the user provided context

171:     Level: advanced

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

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

178: .seealso: MatCreateShell(), MatShellSetOperation(), MatShellSetContext()
179: @*/
180: PetscErrorCode  MatShellGetContext(Mat mat,void *ctx)
181: {
183:   PetscBool      flg;

188:   PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
189:   if (flg) *(void**)ctx = ((Mat_Shell*)(mat->data))->ctx;
190:   else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot get context from non-shell matrix");
191:   return(0);
192: }

194: PetscErrorCode MatDestroy_Shell(Mat mat)
195: {
197:   Mat_Shell      *shell = (Mat_Shell*)mat->data;

200:   if (shell->ops->destroy) {
201:     (*shell->ops->destroy)(mat);
202:   }
203:   VecDestroy(&shell->left);
204:   VecDestroy(&shell->right);
205:   VecDestroy(&shell->dshift);
206:   VecDestroy(&shell->left_work);
207:   VecDestroy(&shell->right_work);
208:   VecDestroy(&shell->left_add_work);
209:   VecDestroy(&shell->right_add_work);
210:   MatDestroy(&shell->axpy);
211:   PetscFree(mat->data);
212:   return(0);
213: }

215: PetscErrorCode MatCopy_Shell(Mat A,Mat B,MatStructure str)
216: {
217:   Mat_Shell       *shellA = (Mat_Shell*)A->data,*shellB = (Mat_Shell*)B->data;
218:   PetscErrorCode  ierr;
219:   PetscBool       matflg;

222:   PetscObjectTypeCompare((PetscObject)B,MATSHELL,&matflg);
223:   if (!matflg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_NOTSAMETYPE,"Output matrix must be a MATSHELL");

225:   PetscMemcpy(B->ops,A->ops,sizeof(struct _MatOps));
226:   PetscMemcpy(shellB->ops,shellA->ops,sizeof(struct _MatShellOps));

228:   if (shellA->ops->copy) {
229:     (*shellA->ops->copy)(A,B,str);
230:   }
231:   shellB->vscale = shellA->vscale;
232:   shellB->vshift = shellA->vshift;
233:   if (shellA->dshift) {
234:     if (!shellB->dshift) {
235:       VecDuplicate(shellA->dshift,&shellB->dshift);
236:     }
237:     VecCopy(shellA->dshift,shellB->dshift);
238:   } else {
239:     VecDestroy(&shellB->dshift);
240:   }
241:   if (shellA->left) {
242:     if (!shellB->left) {
243:       VecDuplicate(shellA->left,&shellB->left);
244:     }
245:     VecCopy(shellA->left,shellB->left);
246:   } else {
247:     VecDestroy(&shellB->left);
248:   }
249:   if (shellA->right) {
250:     if (!shellB->right) {
251:       VecDuplicate(shellA->right,&shellB->right);
252:     }
253:     VecCopy(shellA->right,shellB->right);
254:   } else {
255:     VecDestroy(&shellB->right);
256:   }
257:   MatDestroy(&shellB->axpy);
258:   if (shellA->axpy) {
259:     PetscObjectReference((PetscObject)shellA->axpy);
260:     shellB->axpy        = shellA->axpy;
261:     shellB->axpy_vscale = shellA->axpy_vscale;
262:   }
263:   return(0);
264: }

266: PetscErrorCode MatDuplicate_Shell(Mat mat,MatDuplicateOption op,Mat *M)
267: {
269:   void           *ctx;

272:   MatShellGetContext(mat,&ctx);
273:   MatCreateShell(PetscObjectComm((PetscObject)mat),mat->rmap->n,mat->cmap->n,mat->rmap->N,mat->cmap->N,ctx,M);
274:   MatCopy(mat,*M,SAME_NONZERO_PATTERN);
275:   return(0);
276: }

278: PetscErrorCode MatMult_Shell(Mat A,Vec x,Vec y)
279: {
280:   Mat_Shell        *shell = (Mat_Shell*)A->data;
281:   PetscErrorCode   ierr;
282:   Vec              xx;
283:   PetscObjectState instate,outstate;

286:   MatShellPreScaleRight(A,x,&xx);
287:   PetscObjectStateGet((PetscObject)y, &instate);
288:   (*shell->ops->mult)(A,xx,y);
289:   PetscObjectStateGet((PetscObject)y, &outstate);
290:   if (instate == outstate) {
291:     /* increase the state of the output vector since the user did not update its state themself as should have been done */
292:     PetscObjectStateIncrease((PetscObject)y);
293:   }
294:   MatShellShiftAndScale(A,xx,y);
295:   MatShellPostScaleLeft(A,y);
296: 
297:   if (shell->axpy) {
298:     if (!shell->left_work) {MatCreateVecs(A,&shell->left_work,NULL);}
299:     MatMult(shell->axpy,x,shell->left_work);
300:     VecAXPY(y,shell->axpy_vscale,shell->left_work);
301:   }
302:   return(0);
303: }

305: PetscErrorCode MatMultAdd_Shell(Mat A,Vec x,Vec y,Vec z)
306: {
307:   Mat_Shell      *shell = (Mat_Shell*)A->data;

311:   if (y == z) {
312:     if (!shell->right_add_work) {VecDuplicate(z,&shell->right_add_work);}
313:     MatMult(A,x,shell->right_add_work);
314:     VecAXPY(z,1.0,shell->right_add_work);
315:   } else {
316:     MatMult(A,x,z);
317:     VecAXPY(z,1.0,y);
318:   }
319:   return(0);
320: }

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

330:   MatShellPreScaleLeft(A,x,&xx);
331:   PetscObjectStateGet((PetscObject)y, &instate);
332:   if (!shell->ops->multtranspose) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Have not provided a MatMultTranspose() for this MATSHELL");
333:   (*shell->ops->multtranspose)(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:   MatShellPostScaleRight(A,y);
341:   return(0);
342: }

344: PetscErrorCode MatMultTransposeAdd_Shell(Mat A,Vec x,Vec y,Vec z)
345: {
346:   Mat_Shell      *shell = (Mat_Shell*)A->data;

350:   if (y == z) {
351:     if (!shell->left_add_work) {VecDuplicate(z,&shell->left_add_work);}
352:     MatMultTranspose(A,x,shell->left_add_work);
353:     VecWAXPY(z,1.0,shell->left_add_work,y);
354:   } else {
355:     MatMultTranspose(A,x,z);
356:     VecAXPY(z,1.0,y);
357:   }
358:   return(0);
359: }

361: /*
362:           diag(left)(vscale*A + diag(dshift) + vshift I)diag(right)
363: */
364: PetscErrorCode MatGetDiagonal_Shell(Mat A,Vec v)
365: {
366:   Mat_Shell      *shell = (Mat_Shell*)A->data;

370:   if (shell->ops->getdiagonal) {
371:     (*shell->ops->getdiagonal)(A,v);
372:   } else {
373:     VecSet(v,0.0);
374:   }
375:   VecScale(v,shell->vscale);
376:   if (shell->dshift) {
377:     VecAXPY(v,1.0,shell->dshift);
378:   }
379:   VecShift(v,shell->vshift);
380:   if (shell->left)  {VecPointwiseMult(v,v,shell->left);}
381:   if (shell->right) {VecPointwiseMult(v,v,shell->right);}
382:   if (shell->axpy) {
383:     if (!shell->left_work) {VecDuplicate(v,&shell->left_work);}
384:     MatGetDiagonal(shell->axpy,shell->left_work);
385:     VecAXPY(v,shell->axpy_vscale,shell->left_work);
386:   }
387:   return(0);
388: }

390: PetscErrorCode MatShift_Shell(Mat Y,PetscScalar a)
391: {
392:   Mat_Shell      *shell = (Mat_Shell*)Y->data;

396:   if (shell->left || shell->right) {
397:     if (!shell->dshift) {
398:       VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift);
399:       VecSet(shell->dshift,a);
400:     } else {
401:       if (shell->left)  {VecPointwiseMult(shell->dshift,shell->dshift,shell->left);}
402:       if (shell->right) {VecPointwiseMult(shell->dshift,shell->dshift,shell->right);}
403:       VecShift(shell->dshift,a);
404:     }
405:     if (shell->left)  {VecPointwiseDivide(shell->dshift,shell->dshift,shell->left);}
406:     if (shell->right) {VecPointwiseDivide(shell->dshift,shell->dshift,shell->right);}
407:   } else shell->vshift += a;
408:   return(0);
409: }

411: PetscErrorCode MatDiagonalSet_Shell(Mat A,Vec D,InsertMode ins)
412: {
413:   Mat_Shell      *shell = (Mat_Shell*)A->data;

417:   if (ins == INSERT_VALUES) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE, "Operation not supported with INSERT_VALUES");
418:   if (!shell->dshift) {VecDuplicate(D,&shell->dshift);}
419:   if (shell->left || shell->right) {
420:     if (!shell->right_work) {VecDuplicate(shell->left ? shell->left : shell->right, &shell->right_work);}
421:     if (shell->left && shell->right)  {
422:       VecPointwiseDivide(shell->right_work,D,shell->left);
423:       VecPointwiseDivide(shell->right_work,shell->right_work,shell->right);
424:     } else if (shell->left) {
425:       VecPointwiseDivide(shell->right_work,D,shell->left);
426:     } else {
427:       VecPointwiseDivide(shell->right_work,D,shell->right);
428:     }
429:     if (!shell->dshift) {
430:       VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift);
431:       VecCopy(shell->dshift,shell->right_work);
432:     } else {
433:       VecAXPY(shell->dshift,1.0,shell->right_work);
434:     }
435:   } else {
436:     VecAXPY(shell->dshift,1.0,D);
437:   }
438:   return(0);
439: }

441: PetscErrorCode MatScale_Shell(Mat Y,PetscScalar a)
442: {
443:   Mat_Shell      *shell = (Mat_Shell*)Y->data;
445: 
447:   shell->vscale *= a;
448:   shell->vshift *= a;
449:   if (shell->dshift) {
450:     VecScale(shell->dshift,a);
451:   }
452:   shell->axpy_vscale *= a;
453:   return(0);
454: }

456: static PetscErrorCode MatDiagonalScale_Shell(Mat Y,Vec left,Vec right)
457: {
458:   Mat_Shell      *shell = (Mat_Shell*)Y->data;

462:   if (shell->axpy) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_SUP,"Cannot diagonal scale MATSHELL after MatAXPY operation");
463:   if (left) {
464:     if (!shell->left) {
465:       VecDuplicate(left,&shell->left);
466:       VecCopy(left,shell->left);
467:     } else {
468:       VecPointwiseMult(shell->left,shell->left,left);
469:     }
470:   }
471:   if (right) {
472:     if (!shell->right) {
473:       VecDuplicate(right,&shell->right);
474:       VecCopy(right,shell->right);
475:     } else {
476:       VecPointwiseMult(shell->right,shell->right,right);
477:     }
478:   }
479:   return(0);
480: }

482: PetscErrorCode MatAssemblyEnd_Shell(Mat Y,MatAssemblyType t)
483: {
484:   Mat_Shell      *shell = (Mat_Shell*)Y->data;

488:   if (t == MAT_FINAL_ASSEMBLY) {
489:     shell->vshift = 0.0;
490:     shell->vscale = 1.0;
491:     VecDestroy(&shell->dshift);
492:     VecDestroy(&shell->left);
493:     VecDestroy(&shell->right);
494:     MatDestroy(&shell->axpy);
495:   }
496:   return(0);
497: }

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

501: static PetscErrorCode MatMissingDiagonal_Shell(Mat A,PetscBool  *missing,PetscInt *d)
502: {
504:   *missing = PETSC_FALSE;
505:   return(0);
506: }

508: PetscErrorCode MatAXPY_Shell(Mat Y,PetscScalar a,Mat X,MatStructure str)
509: {
510:   Mat_Shell      *shell = (Mat_Shell*)Y->data;

514:   PetscObjectReference((PetscObject)X);
515:   MatDestroy(&shell->axpy);
516:   shell->axpy        = X;
517:   shell->axpy_vscale = a;
518:   return(0);
519: }

521: static struct _MatOps MatOps_Values = {0,
522:                                        0,
523:                                        0,
524:                                        0,
525:                                 /* 4*/ MatMultAdd_Shell,
526:                                        0,
527:                                        MatMultTransposeAdd_Shell,
528:                                        0,
529:                                        0,
530:                                        0,
531:                                 /*10*/ 0,
532:                                        0,
533:                                        0,
534:                                        0,
535:                                        0,
536:                                 /*15*/ 0,
537:                                        0,
538:                                        MatGetDiagonal_Shell,
539:                                        MatDiagonalScale_Shell,
540:                                        0,
541:                                 /*20*/ 0,
542:                                        MatAssemblyEnd_Shell,
543:                                        0,
544:                                        0,
545:                                 /*24*/ 0,
546:                                        0,
547:                                        0,
548:                                        0,
549:                                        0,
550:                                 /*29*/ 0,
551:                                        0,
552:                                        0,
553:                                        0,
554:                                        0,
555:                                 /*34*/ MatDuplicate_Shell,
556:                                        0,
557:                                        0,
558:                                        0,
559:                                        0,
560:                                 /*39*/ MatAXPY_Shell,
561:                                        0,
562:                                        0,
563:                                        0,
564:                                        MatCopy_Shell,
565:                                 /*44*/ 0,
566:                                        MatScale_Shell,
567:                                        MatShift_Shell,
568:                                        MatDiagonalSet_Shell,
569:                                        0,
570:                                 /*49*/ 0,
571:                                        0,
572:                                        0,
573:                                        0,
574:                                        0,
575:                                 /*54*/ 0,
576:                                        0,
577:                                        0,
578:                                        0,
579:                                        0,
580:                                 /*59*/ 0,
581:                                        MatDestroy_Shell,
582:                                        0,
583:                                        0,
584:                                        0,
585:                                 /*64*/ 0,
586:                                        0,
587:                                        0,
588:                                        0,
589:                                        0,
590:                                 /*69*/ 0,
591:                                        0,
592:                                        MatConvert_Shell,
593:                                        0,
594:                                        0,
595:                                 /*74*/ 0,
596:                                        0,
597:                                        0,
598:                                        0,
599:                                        0,
600:                                 /*79*/ 0,
601:                                        0,
602:                                        0,
603:                                        0,
604:                                        0,
605:                                 /*84*/ 0,
606:                                        0,
607:                                        0,
608:                                        0,
609:                                        0,
610:                                 /*89*/ 0,
611:                                        0,
612:                                        0,
613:                                        0,
614:                                        0,
615:                                 /*94*/ 0,
616:                                        0,
617:                                        0,
618:                                        0,
619:                                        0,
620:                                 /*99*/ 0,
621:                                        0,
622:                                        0,
623:                                        0,
624:                                        0,
625:                                /*104*/ 0,
626:                                        0,
627:                                        0,
628:                                        0,
629:                                        0,
630:                                /*109*/ 0,
631:                                        0,
632:                                        0,
633:                                        0,
634:                                        MatMissingDiagonal_Shell,
635:                                /*114*/ 0,
636:                                        0,
637:                                        0,
638:                                        0,
639:                                        0,
640:                                /*119*/ 0,
641:                                        0,
642:                                        0,
643:                                        0,
644:                                        0,
645:                                /*124*/ 0,
646:                                        0,
647:                                        0,
648:                                        0,
649:                                        0,
650:                                /*129*/ 0,
651:                                        0,
652:                                        0,
653:                                        0,
654:                                        0,
655:                                /*134*/ 0,
656:                                        0,
657:                                        0,
658:                                        0,
659:                                        0,
660:                                /*139*/ 0,
661:                                        0,
662:                                        0
663: };

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

668:   Level: advanced

670: .seealso: MatCreateShell()
671: M*/

673: PETSC_EXTERN PetscErrorCode MatCreate_Shell(Mat A)
674: {
675:   Mat_Shell      *b;

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

681:   PetscNewLog(A,&b);
682:   A->data = (void*)b;

684:   PetscLayoutSetUp(A->rmap);
685:   PetscLayoutSetUp(A->cmap);

687:   b->ctx                 = 0;
688:   b->vshift              = 0.0;
689:   b->vscale              = 1.0;
690:   b->managescalingshifts = PETSC_TRUE;
691:   A->assembled           = PETSC_TRUE;
692:   A->preallocated        = PETSC_FALSE;

694:   PetscObjectChangeTypeName((PetscObject)A,MATSHELL);
695:   return(0);
696: }

698: /*@C
699:    MatCreateShell - Creates a new matrix class for use with a user-defined
700:    private data storage format.

702:   Collective on MPI_Comm

704:    Input Parameters:
705: +  comm - MPI communicator
706: .  m - number of local rows (must be given)
707: .  n - number of local columns (must be given)
708: .  M - number of global rows (may be PETSC_DETERMINE)
709: .  N - number of global columns (may be PETSC_DETERMINE)
710: -  ctx - pointer to data needed by the shell matrix routines

712:    Output Parameter:
713: .  A - the matrix

715:    Level: advanced

717:   Usage:
718: $    extern int mult(Mat,Vec,Vec);
719: $    MatCreateShell(comm,m,n,M,N,ctx,&mat);
720: $    MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
721: $    [ Use matrix for operations that have been set ]
722: $    MatDestroy(mat);

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

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

733:    PETSc requires that matrices and vectors being used for certain
734:    operations are partitioned accordingly.  For example, when
735:    creating a shell matrix, A, that supports parallel matrix-vector
736:    products using MatMult(A,x,y) the user should set the number
737:    of local matrix rows to be the number of local elements of the
738:    corresponding result vector, y. Note that this is information is
739:    required for use of the matrix interface routines, even though
740:    the shell matrix may not actually be physically partitioned.
741:    For example,

743: $
744: $     Vec x, y
745: $     extern int mult(Mat,Vec,Vec);
746: $     Mat A
747: $
748: $     VecCreateMPI(comm,PETSC_DECIDE,M,&y);
749: $     VecCreateMPI(comm,PETSC_DECIDE,N,&x);
750: $     VecGetLocalSize(y,&m);
751: $     VecGetLocalSize(x,&n);
752: $     MatCreateShell(comm,m,n,M,N,ctx,&A);
753: $     MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
754: $     MatMult(A,x,y);
755: $     MatDestroy(A);
756: $     VecDestroy(y); VecDestroy(x);
757: $


760:    MATSHELL handles MatShift(), MatDiagonalSet(), MatDiagonalScale(), MatAXPY(), and MatScale() internally so these
761:    operations cannot be overwritten unless MatShellSetManageScalingShifts() is called.


764:     For rectangular matrices do all the scalings and shifts make sense?

766:     Developers Notes: Regarding shifting and scaling. The general form is

768:           diag(left)(vscale*A + diag(dshift) + vshift I)diag(right)

770:       The order you apply the operations is important. For example if you have a dshift then
771:       apply a MatScale(s) you get s*vscale*A + s*diag(shift). But if you first scale and then shift
772:       you get s*vscale*A + diag(shift)

774:           A is the user provided function. 

776: .keywords: matrix, shell, create

778: .seealso: MatShellSetOperation(), MatHasOperation(), MatShellGetContext(), MatShellSetContext(), MATSHELL, MatShellSetManageScalingShifts()
779: @*/
780: PetscErrorCode  MatCreateShell(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,void *ctx,Mat *A)
781: {

785:   MatCreate(comm,A);
786:   MatSetSizes(*A,m,n,M,N);
787:   MatSetType(*A,MATSHELL);
788:   MatShellSetContext(*A,ctx);
789:   MatSetUp(*A);
790:   return(0);
791: }

793: /*@
794:     MatShellSetContext - sets the context for a shell matrix

796:    Logically Collective on Mat

798:     Input Parameters:
799: +   mat - the shell matrix
800: -   ctx - the context

802:    Level: advanced

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

807: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation()
808: @*/
809: PetscErrorCode  MatShellSetContext(Mat mat,void *ctx)
810: {
811:   Mat_Shell      *shell = (Mat_Shell*)mat->data;
813:   PetscBool      flg;

817:   PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
818:   if (flg) {
819:     shell->ctx = ctx;
820:   } else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot attach context to non-shell matrix");
821:   return(0);
822: }

824: /*@
825:     MatShellSetManageScalingShifts - Allows the user to control the scaling and shift operations of the MATSHELL. Must be called immediately
826:           after MatCreateShell()

828:    Logically Collective on Mat

830:     Input Parameter:
831: .   mat - the shell matrix

833:   Level: advanced

835: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatShellSetOperation()
836: @*/
837: PetscErrorCode MatShellSetManageScalingShifts(Mat A)
838: {
840:   Mat_Shell      *shell = (Mat_Shell*)A->data;
841:   PetscBool      flg;

845:   PetscObjectTypeCompare((PetscObject)A,MATSHELL,&flg);
846:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Can only use with MATSHELL matrices");
847:   shell->managescalingshifts = PETSC_FALSE;
848:   A->ops->diagonalscale = 0;
849:   A->ops->scale         = 0;
850:   A->ops->shift         = 0;
851:   A->ops->diagonalset   = 0;
852:   return(0);
853: }

855: /*@C
856:     MatShellTestMult - Compares the multiply routine provided to the MATSHELL with differencing on a given function.

858:    Logically Collective on Mat

860:     Input Parameters:
861: +   mat - the shell matrix
862: .   f - the function
863: .   base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated
864: -   ctx - an optional context for the function

866:    Output Parameter:
867: .   flg - PETSC_TRUE if the multiply is likely correct

869:    Options Database:
870: .   -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference

872:    Level: advanced

874:    Fortran Notes: Not supported from Fortran

876: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellTestMultTranspose()
877: @*/
878: PetscErrorCode  MatShellTestMult(Mat mat,PetscErrorCode (*f)(void*,Vec,Vec),Vec base,void *ctx,PetscBool *flg)
879: {
881:   PetscInt       m,n;
882:   Mat            mf,Dmf,Dmat,Ddiff;
883:   PetscReal      Diffnorm,Dmfnorm;
884:   PetscBool      v = PETSC_FALSE, flag = PETSC_TRUE;

888:   PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_view",&v);
889:   MatGetLocalSize(mat,&m,&n);
890:   MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,PETSC_DECIDE,PETSC_DECIDE,&mf);
891:   MatMFFDSetFunction(mf,f,ctx);
892:   MatMFFDSetBase(mf,base,NULL);

894:   MatComputeExplicitOperator(mf,&Dmf);
895:   MatComputeExplicitOperator(mat,&Dmat);

897:   MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff);
898:   MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN);
899:   MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm);
900:   MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm);
901:   if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) {
902:     flag = PETSC_FALSE;
903:     if (v) {
904:       PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL and matrix free multiple appear to produce different results.\n Norm Ratio %g Difference results followed by finite difference one\n",(double)(Diffnorm/Dmfnorm));
905:       MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_view");
906:       MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_view");
907:       MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_view");
908:     }
909:   } else if (v) {
910:     PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL and matrix free multiple appear to produce the same results\n");
911:   }
912:   if (flg) *flg = flag;
913:   MatDestroy(&Ddiff);
914:   MatDestroy(&mf);
915:   MatDestroy(&Dmf);
916:   MatDestroy(&Dmat);
917:   return(0);
918: }

920: /*@C
921:     MatShellTestMultTranpose - Compares the multiply transpose routine provided to the MATSHELL with differencing on a given function.

923:    Logically Collective on Mat

925:     Input Parameters:
926: +   mat - the shell matrix
927: .   f - the function
928: .   base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated
929: -   ctx - an optional context for the function

931:    Output Parameter:
932: .   flg - PETSC_TRUE if the multiply is likely correct

934:    Options Database:
935: .   -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference

937:    Level: advanced

939:    Fortran Notes: Not supported from Fortran

941: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellTestMult()
942: @*/
943: PetscErrorCode  MatShellTestMultTranspose(Mat mat,PetscErrorCode (*f)(void*,Vec,Vec),Vec base,void *ctx,PetscBool *flg)
944: {
946:   Vec            x,y,z;
947:   PetscInt       m,n,M,N;
948:   Mat            mf,Dmf,Dmat,Ddiff;
949:   PetscReal      Diffnorm,Dmfnorm;
950:   PetscBool      v = PETSC_FALSE, flag = PETSC_TRUE;

954:   PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_transpose_view",&v);
955:   MatCreateVecs(mat,&x,&y);
956:   VecDuplicate(y,&z);
957:   MatGetLocalSize(mat,&m,&n);
958:   MatGetSize(mat,&M,&N);
959:   MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,M,N,&mf);
960:   MatMFFDSetFunction(mf,f,ctx);
961:   MatMFFDSetBase(mf,base,NULL);
962:   MatComputeExplicitOperator(mf,&Dmf);
963:   MatTranspose(Dmf,MAT_INPLACE_MATRIX,&Dmf);
964:   MatComputeExplicitOperatorTranspose(mat,&Dmat);

966:   MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff);
967:   MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN);
968:   MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm);
969:   MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm);
970:   if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) {
971:     flag = PETSC_FALSE;
972:     if (v) {
973:       PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL and matrix free multiple appear to produce different results.\n Norm Ratio %g Difference results followed by finite difference one\n",(double)(Diffnorm/Dmfnorm));
974:       MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");
975:       MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");
976:       MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");
977:     }
978:   } else if (v) {
979:     PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL transpose and matrix free multiple appear to produce the same results\n");
980:   }
981:   if (flg) *flg = flag;
982:   MatDestroy(&mf);
983:   MatDestroy(&Dmat);
984:   MatDestroy(&Ddiff);
985:   MatDestroy(&Dmf);
986:   VecDestroy(&x);
987:   VecDestroy(&y);
988:   VecDestroy(&z);
989:   return(0);
990: }

992: /*@C
993:     MatShellSetOperation - Allows user to set a matrix operation for a shell matrix.

995:    Logically Collective on Mat

997:     Input Parameters:
998: +   mat - the shell matrix
999: .   op - the name of the operation
1000: -   f - the function that provides the operation.

1002:    Level: advanced

1004:     Usage:
1005: $      extern PetscErrorCode usermult(Mat,Vec,Vec);
1006: $      MatCreateShell(comm,m,n,M,N,ctx,&A);
1007: $      MatShellSetOperation(A,MATOP_MULT,(void(*)(void))usermult);

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

1015:     All user-provided functions (except for MATOP_DESTROY) should have the same calling
1016:     sequence as the usual matrix interface routines, since they
1017:     are intended to be accessed via the usual matrix interface
1018:     routines, e.g.,
1019: $       MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec)

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

1024:     Within each user-defined routine, the user should call
1025:     MatShellGetContext() to obtain the user-defined context that was
1026:     set by MatCreateShell().

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

1031:     Use MatSetOperation() to set an operation for any matrix type

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

1035: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatSetOperation(), MatShellSetManageScalingShifts()
1036: @*/
1037: PetscErrorCode  MatShellSetOperation(Mat mat,MatOperation op,void (*f)(void))
1038: {
1040:   PetscBool      flg;
1041:   Mat_Shell      *shell = (Mat_Shell*)mat->data;

1045:   PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
1046:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Can only use with MATSHELL matrices");
1047:   switch (op) {
1048:   case MATOP_DESTROY:
1049:     shell->ops->destroy = (PetscErrorCode (*)(Mat))f;
1050:     break;
1051:   case MATOP_DIAGONAL_SET:
1052:   case MATOP_DIAGONAL_SCALE:
1053:   case MATOP_SHIFT:
1054:   case MATOP_SCALE:
1055:   case MATOP_AXPY:
1056:     if (shell->managescalingshifts) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"MATSHELL is managing scalings and shifts, see MatShellSetManageScalingShifts()");
1057:     (((void(**)(void))mat->ops)[op]) = f;
1058:     break;
1059:   case MATOP_GET_DIAGONAL:
1060:     if (shell->managescalingshifts) shell->ops->getdiagonal = (PetscErrorCode (*)(Mat,Vec))f;
1061:     else mat->ops->getdiagonal = (PetscErrorCode (*)(Mat,Vec))f;
1062:   case MATOP_VIEW:
1063:     if (!mat->ops->viewnative) {
1064:       mat->ops->viewnative = mat->ops->view;
1065:     }
1066:     mat->ops->view = (PetscErrorCode (*)(Mat,PetscViewer))f;
1067:     break;
1068:   case MATOP_MULT:
1069:     if (shell->managescalingshifts){
1070:       shell->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1071:       mat->ops->mult   = MatMult_Shell;
1072:     } else mat->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1073:     break;
1074:   case MATOP_MULT_TRANSPOSE:
1075:     if (shell->managescalingshifts) {
1076:       shell->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1077:       mat->ops->multtranspose   = MatMultTranspose_Shell;
1078:     } else mat->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1079:     break;
1080:   default:
1081:     (((void(**)(void))mat->ops)[op]) = f;
1082:   }
1083:   return(0);
1084: }

1086: /*@C
1087:     MatShellGetOperation - Gets a matrix function for a shell matrix.

1089:     Not Collective

1091:     Input Parameters:
1092: +   mat - the shell matrix
1093: -   op - the name of the operation

1095:     Output Parameter:
1096: .   f - the function that provides the operation.

1098:     Level: advanced

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

1106:     All user-provided functions have the same calling
1107:     sequence as the usual matrix interface routines, since they
1108:     are intended to be accessed via the usual matrix interface
1109:     routines, e.g.,
1110: $       MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec)

1112:     Within each user-defined routine, the user should call
1113:     MatShellGetContext() to obtain the user-defined context that was
1114:     set by MatCreateShell().

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

1118: .seealso: MatCreateShell(), MatShellGetContext(), MatShellSetOperation(), MatShellSetContext()
1119: @*/
1120: PetscErrorCode  MatShellGetOperation(Mat mat,MatOperation op,void(**f)(void))
1121: {
1123:   PetscBool      flg;

1127:   switch (op) {
1128:   case MATOP_DESTROY:
1129:     PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
1130:     if (flg) {
1131:       Mat_Shell *shell = (Mat_Shell*)mat->data;
1132:       *f = (void (*)(void))shell->ops->destroy;
1133:     } else {
1134:       *f = (void (*)(void))mat->ops->destroy;
1135:     }
1136:     break;
1137:   case MATOP_DIAGONAL_SET:
1138:     PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
1139:     if (flg) {
1140:       Mat_Shell *shell = (Mat_Shell*)mat->data;
1141:       *f = (void (*)(void))shell->ops->diagonalset;
1142:     } else {
1143:       *f = (void (*)(void))mat->ops->diagonalset;
1144:     }
1145:     break;
1146:   case MATOP_VIEW:
1147:     *f = (void (*)(void))mat->ops->view;
1148:     break;
1149:   default:
1150:     *f = (((void (**)(void))mat->ops)[op]);
1151:   }
1152:   return(0);
1153: }

1155: /*@C
1156:     MatSetOperation - Allows user to set a matrix operation for any matrix type

1158:    Logically Collective on Mat

1160:     Input Parameters:
1161: +   mat - the shell matrix
1162: .   op - the name of the operation
1163: -   f - the function that provides the operation.

1165:    Level: developer

1167:     Usage:
1168: $      extern PetscErrorCode usermult(Mat,Vec,Vec);
1169: $      MatCreateXXX(comm,...&A);
1170: $      MatShellSetOperation(A,MATOP_MULT,(void(*)(void))usermult);

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

1178:     All user-provided functions (except for MATOP_DESTROY) should have the same calling
1179:     sequence as the usual matrix interface routines, since they
1180:     are intended to be accessed via the usual matrix interface
1181:     routines, e.g.,
1182: $       MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec)

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

1187:     This routine is distinct from MatShellSetOperation() in that it can be called on any matrix type.

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

1191: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatShellSetOperation()
1192: @*/
1193: PetscErrorCode  MatSetOperation(Mat mat,MatOperation op,void (*f)(void))
1194: {
1197:   (((void(**)(void))mat->ops)[op]) = f;
1198:   return(0);
1199: }