Actual source code: shell.c

petsc-3.11.4 2019-09-28
Report Typos and Errors

  2: /*
  3:    This provides a simple shell for Fortran (and C programmers) to
  4:   create a very simple matrix class for use with KSP without coding
  5:   much of anything.
  6: */

  8:  #include <petsc/private/matimpl.h>

 10: struct _MatShellOps {
 11:   /*  3 */ PetscErrorCode (*mult)(Mat,Vec,Vec);
 12:   /*  5 */ PetscErrorCode (*multtranspose)(Mat,Vec,Vec);
 13:   /* 17 */ PetscErrorCode (*getdiagonal)(Mat,Vec);
 14:   /* 43 */ PetscErrorCode (*copy)(Mat,Mat,MatStructure);
 15:   /* 60 */ PetscErrorCode (*destroy)(Mat);
 16: };

 18: typedef struct {
 19:   struct _MatShellOps ops[1];

 21:   PetscScalar vscale,vshift;
 22:   Vec         dshift;
 23:   Vec         left,right;
 24:   Vec         left_work,right_work;
 25:   Vec         left_add_work,right_add_work;
 26:   Mat         axpy;
 27:   PetscScalar axpy_vscale;
 28:   PetscBool   managescalingshifts;                   /* The user will manage the scaling and shifts for the MATSHELL, not the default */
 29:   void        *ctx;
 30: } Mat_Shell;

 32: /*
 33:       xx = diag(left)*x
 34: */
 35: static PetscErrorCode MatShellPreScaleLeft(Mat A,Vec x,Vec *xx)
 36: {
 37:   Mat_Shell      *shell = (Mat_Shell*)A->data;

 41:   *xx = NULL;
 42:   if (!shell->left) {
 43:     *xx = x;
 44:   } else {
 45:     if (!shell->left_work) {VecDuplicate(shell->left,&shell->left_work);}
 46:     VecPointwiseMult(shell->left_work,x,shell->left);
 47:     *xx  = shell->left_work;
 48:   }
 49:   return(0);
 50: }

 52: /*
 53:      xx = diag(right)*x
 54: */
 55: static PetscErrorCode MatShellPreScaleRight(Mat A,Vec x,Vec *xx)
 56: {
 57:   Mat_Shell      *shell = (Mat_Shell*)A->data;

 61:   *xx = NULL;
 62:   if (!shell->right) {
 63:     *xx = x;
 64:   } else {
 65:     if (!shell->right_work) {VecDuplicate(shell->right,&shell->right_work);}
 66:     VecPointwiseMult(shell->right_work,x,shell->right);
 67:     *xx  = shell->right_work;
 68:   }
 69:   return(0);
 70: }

 72: /*
 73:     x = diag(left)*x
 74: */
 75: static PetscErrorCode MatShellPostScaleLeft(Mat A,Vec x)
 76: {
 77:   Mat_Shell      *shell = (Mat_Shell*)A->data;

 81:   if (shell->left) {VecPointwiseMult(x,x,shell->left);}
 82:   return(0);
 83: }

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

 94:   if (shell->right) {VecPointwiseMult(x,x,shell->right);}
 95:   return(0);
 96: }

 98: /*
 99:          Y = vscale*Y + diag(dshift)*X + vshift*X

101:          On input Y already contains A*x
102: */
103: static PetscErrorCode MatShellShiftAndScale(Mat A,Vec X,Vec Y)
104: {
105:   Mat_Shell      *shell = (Mat_Shell*)A->data;

109:   if (shell->dshift) {          /* get arrays because there is no VecPointwiseMultAdd() */
110:     PetscInt          i,m;
111:     const PetscScalar *x,*d;
112:     PetscScalar       *y;
113:     VecGetLocalSize(X,&m);
114:     VecGetArrayRead(shell->dshift,&d);
115:     VecGetArrayRead(X,&x);
116:     VecGetArray(Y,&y);
117:     for (i=0; i<m; i++) y[i] = shell->vscale*y[i] + d[i]*x[i];
118:     VecRestoreArrayRead(shell->dshift,&d);
119:     VecRestoreArrayRead(X,&x);
120:     VecRestoreArray(Y,&y);
121:   } else {
122:     VecScale(Y,shell->vscale);
123:   }
124:   if (shell->vshift != 0.0) {VecAXPY(Y,shell->vshift,X);} /* if test is for non-square matrices */
125:   return(0);
126: }

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

131:     Not Collective

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

136:     Output Parameter:
137: .   ctx - the user provided context

139:     Level: advanced

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

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

147: .seealso: MatCreateShell(), MatShellSetOperation(), MatShellSetContext()
148: @*/
149: PetscErrorCode  MatShellGetContext(Mat mat,void *ctx)
150: {
152:   PetscBool      flg;

157:   PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
158:   if (flg) *(void**)ctx = ((Mat_Shell*)(mat->data))->ctx;
159:   else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot get context from non-shell matrix");
160:   return(0);
161: }

163: PetscErrorCode MatDestroy_Shell(Mat mat)
164: {
166:   Mat_Shell      *shell = (Mat_Shell*)mat->data;

169:   if (shell->ops->destroy) {
170:     (*shell->ops->destroy)(mat);
171:   }
172:   VecDestroy(&shell->left);
173:   VecDestroy(&shell->right);
174:   VecDestroy(&shell->dshift);
175:   VecDestroy(&shell->left_work);
176:   VecDestroy(&shell->right_work);
177:   VecDestroy(&shell->left_add_work);
178:   VecDestroy(&shell->right_add_work);
179:   MatDestroy(&shell->axpy);
180:   PetscFree(mat->data);
181:   return(0);
182: }

184: PetscErrorCode MatCopy_Shell(Mat A,Mat B,MatStructure str)
185: {
186:   Mat_Shell       *shellA = (Mat_Shell*)A->data,*shellB = (Mat_Shell*)B->data;
187:   PetscErrorCode  ierr;
188:   PetscBool       matflg;

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

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

197:   if (shellA->ops->copy) {
198:     (*shellA->ops->copy)(A,B,str);
199:   }
200:   shellB->vscale = shellA->vscale;
201:   shellB->vshift = shellA->vshift;
202:   if (shellA->dshift) {
203:     if (!shellB->dshift) {
204:       VecDuplicate(shellA->dshift,&shellB->dshift);
205:     }
206:     VecCopy(shellA->dshift,shellB->dshift);
207:   } else {
208:     VecDestroy(&shellB->dshift);
209:   }
210:   if (shellA->left) {
211:     if (!shellB->left) {
212:       VecDuplicate(shellA->left,&shellB->left);
213:     }
214:     VecCopy(shellA->left,shellB->left);
215:   } else {
216:     VecDestroy(&shellB->left);
217:   }
218:   if (shellA->right) {
219:     if (!shellB->right) {
220:       VecDuplicate(shellA->right,&shellB->right);
221:     }
222:     VecCopy(shellA->right,shellB->right);
223:   } else {
224:     VecDestroy(&shellB->right);
225:   }
226:   MatDestroy(&shellB->axpy);
227:   if (shellA->axpy) {
228:     PetscObjectReference((PetscObject)shellA->axpy);
229:     shellB->axpy        = shellA->axpy;
230:     shellB->axpy_vscale = shellA->axpy_vscale;
231:   }
232:   return(0);
233: }

235: PetscErrorCode MatDuplicate_Shell(Mat mat,MatDuplicateOption op,Mat *M)
236: {
238:   void           *ctx;

241:   MatShellGetContext(mat,&ctx);
242:   MatCreateShell(PetscObjectComm((PetscObject)mat),mat->rmap->n,mat->cmap->n,mat->rmap->N,mat->cmap->N,ctx,M);
243:   if (op != MAT_DO_NOT_COPY_VALUES) {
244:     MatCopy(mat,*M,SAME_NONZERO_PATTERN);
245:   }
246:   return(0);
247: }

249: PetscErrorCode MatMult_Shell(Mat A,Vec x,Vec y)
250: {
251:   Mat_Shell        *shell = (Mat_Shell*)A->data;
252:   PetscErrorCode   ierr;
253:   Vec              xx;
254:   PetscObjectState instate,outstate;

257:   MatShellPreScaleRight(A,x,&xx);
258:   PetscObjectStateGet((PetscObject)y, &instate);
259:   if (!shell->ops->mult) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Have not provided a MatMult() for this MATSHELL");
260:   (*shell->ops->mult)(A,xx,y);
261:   PetscObjectStateGet((PetscObject)y, &outstate);
262:   if (instate == outstate) {
263:     /* increase the state of the output vector since the user did not update its state themself as should have been done */
264:     PetscObjectStateIncrease((PetscObject)y);
265:   }
266:   MatShellShiftAndScale(A,xx,y);
267:   MatShellPostScaleLeft(A,y);

269:   if (shell->axpy) {
270:     if (!shell->left_work) {MatCreateVecs(A,&shell->left_work,NULL);}
271:     MatMult(shell->axpy,x,shell->left_work);
272:     VecAXPY(y,shell->axpy_vscale,shell->left_work);
273:   }
274:   return(0);
275: }

277: PetscErrorCode MatMultAdd_Shell(Mat A,Vec x,Vec y,Vec z)
278: {
279:   Mat_Shell      *shell = (Mat_Shell*)A->data;

283:   if (y == z) {
284:     if (!shell->right_add_work) {VecDuplicate(z,&shell->right_add_work);}
285:     MatMult(A,x,shell->right_add_work);
286:     VecAXPY(z,1.0,shell->right_add_work);
287:   } else {
288:     MatMult(A,x,z);
289:     VecAXPY(z,1.0,y);
290:   }
291:   return(0);
292: }

294: PetscErrorCode MatMultTranspose_Shell(Mat A,Vec x,Vec y)
295: {
296:   Mat_Shell        *shell = (Mat_Shell*)A->data;
297:   PetscErrorCode   ierr;
298:   Vec              xx;
299:   PetscObjectState instate,outstate;

302:   MatShellPreScaleLeft(A,x,&xx);
303:   PetscObjectStateGet((PetscObject)y, &instate);
304:   if (!shell->ops->multtranspose) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Have not provided a MatMultTranspose() for this MATSHELL");
305:   (*shell->ops->multtranspose)(A,xx,y);
306:   PetscObjectStateGet((PetscObject)y, &outstate);
307:   if (instate == outstate) {
308:     /* increase the state of the output vector since the user did not update its state themself as should have been done */
309:     PetscObjectStateIncrease((PetscObject)y);
310:   }
311:   MatShellShiftAndScale(A,xx,y);
312:   MatShellPostScaleRight(A,y);

314:   if (shell->axpy) {
315:     if (!shell->right_work) {MatCreateVecs(A,NULL,&shell->right_work);}
316:     MatMultTranspose(shell->axpy,x,shell->right_work);
317:     VecAXPY(y,shell->axpy_vscale,shell->right_work);
318:   }
319:   return(0);
320: }

322: PetscErrorCode MatMultTransposeAdd_Shell(Mat A,Vec x,Vec y,Vec z)
323: {
324:   Mat_Shell      *shell = (Mat_Shell*)A->data;

328:   if (y == z) {
329:     if (!shell->left_add_work) {VecDuplicate(z,&shell->left_add_work);}
330:     MatMultTranspose(A,x,shell->left_add_work);
331:     VecAXPY(z,1.0,shell->left_add_work);
332:   } else {
333:     MatMultTranspose(A,x,z);
334:     VecAXPY(z,1.0,y);
335:   }
336:   return(0);
337: }

339: /*
340:           diag(left)(vscale*A + diag(dshift) + vshift I)diag(right)
341: */
342: PetscErrorCode MatGetDiagonal_Shell(Mat A,Vec v)
343: {
344:   Mat_Shell      *shell = (Mat_Shell*)A->data;

348:   if (shell->ops->getdiagonal) {
349:     (*shell->ops->getdiagonal)(A,v);
350:   } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Must provide shell matrix with routine to return diagonal using\nMatShellSetOperation(S,MATOP_GET_DIAGONAL,...)");
351:   VecScale(v,shell->vscale);
352:   if (shell->dshift) {
353:     VecAXPY(v,1.0,shell->dshift);
354:   }
355:   VecShift(v,shell->vshift);
356:   if (shell->left)  {VecPointwiseMult(v,v,shell->left);}
357:   if (shell->right) {VecPointwiseMult(v,v,shell->right);}
358:   if (shell->axpy) {
359:     if (!shell->left_work) {VecDuplicate(v,&shell->left_work);}
360:     MatGetDiagonal(shell->axpy,shell->left_work);
361:     VecAXPY(v,shell->axpy_vscale,shell->left_work);
362:   }
363:   return(0);
364: }

366: PetscErrorCode MatShift_Shell(Mat Y,PetscScalar a)
367: {
368:   Mat_Shell      *shell = (Mat_Shell*)Y->data;

372:   if (shell->left || shell->right) {
373:     if (!shell->dshift) {
374:       VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift);
375:       VecSet(shell->dshift,a);
376:     } else {
377:       if (shell->left)  {VecPointwiseMult(shell->dshift,shell->dshift,shell->left);}
378:       if (shell->right) {VecPointwiseMult(shell->dshift,shell->dshift,shell->right);}
379:       VecShift(shell->dshift,a);
380:     }
381:     if (shell->left)  {VecPointwiseDivide(shell->dshift,shell->dshift,shell->left);}
382:     if (shell->right) {VecPointwiseDivide(shell->dshift,shell->dshift,shell->right);}
383:   } else shell->vshift += a;
384:   return(0);
385: }

387: PetscErrorCode MatDiagonalSet_Shell_Private(Mat A,Vec D,PetscScalar s)
388: {
389:   Mat_Shell      *shell = (Mat_Shell*)A->data;

393:   if (!shell->dshift) {VecDuplicate(D,&shell->dshift);}
394:   if (shell->left || shell->right) {
395:     if (!shell->right_work) {VecDuplicate(shell->left ? shell->left : shell->right, &shell->right_work);}
396:     if (shell->left && shell->right)  {
397:       VecPointwiseDivide(shell->right_work,D,shell->left);
398:       VecPointwiseDivide(shell->right_work,shell->right_work,shell->right);
399:     } else if (shell->left) {
400:       VecPointwiseDivide(shell->right_work,D,shell->left);
401:     } else {
402:       VecPointwiseDivide(shell->right_work,D,shell->right);
403:     }
404:     if (!shell->dshift) {
405:       VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift);
406:       VecCopy(shell->dshift,shell->right_work);
407:     } else {
408:       VecAXPY(shell->dshift,s,shell->right_work);
409:     }
410:   } else {
411:     VecAXPY(shell->dshift,s,D);
412:   }
413:   return(0);
414: }

416: PetscErrorCode MatDiagonalSet_Shell(Mat A,Vec D,InsertMode ins)
417: {
418:   Vec            d;

422:   if (ins == INSERT_VALUES) {
423:     if (!A->ops->getdiagonal) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Operation MATOP_GETDIAGONAL must be set first");
424:     VecDuplicate(D,&d);
425:     MatGetDiagonal(A,d);
426:     MatDiagonalSet_Shell_Private(A,d,-1.);
427:     MatDiagonalSet_Shell_Private(A,D,1.);
428:     VecDestroy(&d);
429:   } else {
430:     MatDiagonalSet_Shell_Private(A,D,1.);
431:   }
432:   return(0);
433: }

435: PetscErrorCode MatScale_Shell(Mat Y,PetscScalar a)
436: {
437:   Mat_Shell      *shell = (Mat_Shell*)Y->data;

441:   shell->vscale *= a;
442:   shell->vshift *= a;
443:   if (shell->dshift) {
444:     VecScale(shell->dshift,a);
445:   }
446:   shell->axpy_vscale *= a;
447:   return(0);
448: }

450: static PetscErrorCode MatDiagonalScale_Shell(Mat Y,Vec left,Vec right)
451: {
452:   Mat_Shell      *shell = (Mat_Shell*)Y->data;

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

476: PetscErrorCode MatAssemblyEnd_Shell(Mat Y,MatAssemblyType t)
477: {
478:   Mat_Shell      *shell = (Mat_Shell*)Y->data;

482:   if (t == MAT_FINAL_ASSEMBLY) {
483:     shell->vshift = 0.0;
484:     shell->vscale = 1.0;
485:     VecDestroy(&shell->dshift);
486:     VecDestroy(&shell->left);
487:     VecDestroy(&shell->right);
488:     MatDestroy(&shell->axpy);
489:   }
490:   return(0);
491: }

493: static PetscErrorCode MatMissingDiagonal_Shell(Mat A,PetscBool  *missing,PetscInt *d)
494: {
496:   *missing = PETSC_FALSE;
497:   return(0);
498: }

500: PetscErrorCode MatAXPY_Shell(Mat Y,PetscScalar a,Mat X,MatStructure str)
501: {
502:   Mat_Shell      *shell = (Mat_Shell*)Y->data;

506:   PetscObjectReference((PetscObject)X);
507:   MatDestroy(&shell->axpy);
508:   shell->axpy        = X;
509:   shell->axpy_vscale = a;
510:   return(0);
511: }

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

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

660:   Level: advanced

662: .seealso: MatCreateShell()
663: M*/

665: PETSC_EXTERN PetscErrorCode MatCreate_Shell(Mat A)
666: {
667:   Mat_Shell      *b;

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

673:   PetscNewLog(A,&b);
674:   A->data = (void*)b;

676:   PetscLayoutSetUp(A->rmap);
677:   PetscLayoutSetUp(A->cmap);

679:   b->ctx                 = 0;
680:   b->vshift              = 0.0;
681:   b->vscale              = 1.0;
682:   b->managescalingshifts = PETSC_TRUE;
683:   A->assembled           = PETSC_TRUE;
684:   A->preallocated        = PETSC_FALSE;

686:   PetscObjectChangeTypeName((PetscObject)A,MATSHELL);
687:   return(0);
688: }

690: /*@C
691:    MatCreateShell - Creates a new matrix class for use with a user-defined
692:    private data storage format.

694:   Collective on MPI_Comm

696:    Input Parameters:
697: +  comm - MPI communicator
698: .  m - number of local rows (must be given)
699: .  n - number of local columns (must be given)
700: .  M - number of global rows (may be PETSC_DETERMINE)
701: .  N - number of global columns (may be PETSC_DETERMINE)
702: -  ctx - pointer to data needed by the shell matrix routines

704:    Output Parameter:
705: .  A - the matrix

707:    Level: advanced

709:   Usage:
710: $    extern int mult(Mat,Vec,Vec);
711: $    MatCreateShell(comm,m,n,M,N,ctx,&mat);
712: $    MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
713: $    [ Use matrix for operations that have been set ]
714: $    MatDestroy(mat);

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

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

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

736: $
737: $     Vec x, y
738: $     extern int mult(Mat,Vec,Vec);
739: $     Mat A
740: $
741: $     VecCreateMPI(comm,PETSC_DECIDE,M,&y);
742: $     VecCreateMPI(comm,PETSC_DECIDE,N,&x);
743: $     VecGetLocalSize(y,&m);
744: $     VecGetLocalSize(x,&n);
745: $     MatCreateShell(comm,m,n,M,N,ctx,&A);
746: $     MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
747: $     MatMult(A,x,y);
748: $     MatDestroy(A);
749: $     VecDestroy(y); VecDestroy(x);
750: $


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


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

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

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

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

768:           A is the user provided function.

770:    KSP/PC uses changes in the Mat's "state" to decide if preconditioners need to be rebuilt: PCSetUp() only calls the setup() for
771:    for the PC implementation if the Mat state has increased from the previous call. Thus to get changes in a MATSHELL to trigger
772:    an update in the preconditioner you must call MatAssemblyBegin()/MatAssemblyEnd() or PetscObjectStateIncrease((PetscObject)mat);
773:    each time the MATSHELL matrix has changed.

775:    Calling MatAssemblyBegin()/MatAssemblyEnd() on a MATSHELL removes any previously supplied shift and scales that were provided
776:    with MatDiagonalSet(), MatShift(), MatScale(), or MatDiagonalScale().

778: .keywords: matrix, shell, create

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

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

795: /*@
796:     MatShellSetContext - sets the context for a shell matrix

798:    Logically Collective on Mat

800:     Input Parameters:
801: +   mat - the shell matrix
802: -   ctx - the context

804:    Level: advanced

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

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

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

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

831:    Logically Collective on Mat

833:     Input Parameter:
834: .   mat - the shell matrix

836:   Level: advanced

838: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatShellSetOperation()
839: @*/
840: PetscErrorCode MatShellSetManageScalingShifts(Mat A)
841: {
843:   Mat_Shell      *shell;
844:   PetscBool      flg;

848:   PetscObjectTypeCompare((PetscObject)A,MATSHELL,&flg);
849:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Can only use with MATSHELL matrices");
850:   shell = (Mat_Shell*)A->data;
851:   shell->managescalingshifts = PETSC_FALSE;
852:   A->ops->diagonalset   = NULL;
853:   A->ops->diagonalscale = NULL;
854:   A->ops->scale         = NULL;
855:   A->ops->shift         = NULL;
856:   A->ops->axpy          = NULL;
857:   return(0);
858: }

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

863:    Logically Collective on Mat

865:     Input Parameters:
866: +   mat - the shell matrix
867: .   f - the function
868: .   base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated
869: -   ctx - an optional context for the function

871:    Output Parameter:
872: .   flg - PETSC_TRUE if the multiply is likely correct

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

877:    Level: advanced

879:    Fortran Notes:
880:     Not supported from Fortran

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

894:   PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_view",&v);
895:   MatGetLocalSize(mat,&m,&n);
896:   MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,PETSC_DECIDE,PETSC_DECIDE,&mf);
897:   MatMFFDSetFunction(mf,f,ctx);
898:   MatMFFDSetBase(mf,base,NULL);

900:   MatComputeExplicitOperator(mf,&Dmf);
901:   MatComputeExplicitOperator(mat,&Dmat);

903:   MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff);
904:   MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN);
905:   MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm);
906:   MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm);
907:   if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) {
908:     flag = PETSC_FALSE;
909:     if (v) {
910:       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));
911:       MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_view");
912:       MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_view");
913:       MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_view");
914:     }
915:   } else if (v) {
916:     PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL and matrix free multiple appear to produce the same results\n");
917:   }
918:   if (flg) *flg = flag;
919:   MatDestroy(&Ddiff);
920:   MatDestroy(&mf);
921:   MatDestroy(&Dmf);
922:   MatDestroy(&Dmat);
923:   return(0);
924: }

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

929:    Logically Collective on Mat

931:     Input Parameters:
932: +   mat - the shell matrix
933: .   f - the function
934: .   base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated
935: -   ctx - an optional context for the function

937:    Output Parameter:
938: .   flg - PETSC_TRUE if the multiply is likely correct

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

943:    Level: advanced

945:    Fortran Notes:
946:     Not supported from Fortran

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

961:   PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_transpose_view",&v);
962:   MatCreateVecs(mat,&x,&y);
963:   VecDuplicate(y,&z);
964:   MatGetLocalSize(mat,&m,&n);
965:   MatGetSize(mat,&M,&N);
966:   MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,M,N,&mf);
967:   MatMFFDSetFunction(mf,f,ctx);
968:   MatMFFDSetBase(mf,base,NULL);
969:   MatComputeExplicitOperator(mf,&Dmf);
970:   MatTranspose(Dmf,MAT_INPLACE_MATRIX,&Dmf);
971:   MatComputeExplicitOperatorTranspose(mat,&Dmat);

973:   MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff);
974:   MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN);
975:   MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm);
976:   MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm);
977:   if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) {
978:     flag = PETSC_FALSE;
979:     if (v) {
980:       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));
981:       MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");
982:       MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");
983:       MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");
984:     }
985:   } else if (v) {
986:     PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL transpose and matrix free multiple appear to produce the same results\n");
987:   }
988:   if (flg) *flg = flag;
989:   MatDestroy(&mf);
990:   MatDestroy(&Dmat);
991:   MatDestroy(&Ddiff);
992:   MatDestroy(&Dmf);
993:   VecDestroy(&x);
994:   VecDestroy(&y);
995:   VecDestroy(&z);
996:   return(0);
997: }

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

1002:    Logically Collective on Mat

1004:     Input Parameters:
1005: +   mat - the shell matrix
1006: .   op - the name of the operation
1007: -   f - the function that provides the operation.

1009:    Level: advanced

1011:     Usage:
1012: $      extern PetscErrorCode usermult(Mat,Vec,Vec);
1013: $      MatCreateShell(comm,m,n,M,N,ctx,&A);
1014: $      MatShellSetOperation(A,MATOP_MULT,(void(*)(void))usermult);

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

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

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

1031:     Within each user-defined routine, the user should call
1032:     MatShellGetContext() to obtain the user-defined context that was
1033:     set by MatCreateShell().

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

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

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

1043: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatSetOperation(), MatShellSetManageScalingShifts()
1044: @*/
1045: PetscErrorCode MatShellSetOperation(Mat mat,MatOperation op,void (*f)(void))
1046: {
1047:   PetscBool      flg;
1048:   Mat_Shell      *shell;

1053:   PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
1054:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Can only use with MATSHELL matrices");
1055:   shell = (Mat_Shell*)mat->data;

1057:   switch (op) {
1058:   case MATOP_DESTROY:
1059:     shell->ops->destroy = (PetscErrorCode (*)(Mat))f;
1060:     break;
1061:   case MATOP_VIEW:
1062:     if (!mat->ops->viewnative) {
1063:       mat->ops->viewnative = mat->ops->view;
1064:     }
1065:     mat->ops->view = (PetscErrorCode (*)(Mat,PetscViewer))f;
1066:     break;
1067:   case MATOP_COPY:
1068:     shell->ops->copy = (PetscErrorCode (*)(Mat,Mat,MatStructure))f;
1069:     break;
1070:   case MATOP_DIAGONAL_SET:
1071:   case MATOP_DIAGONAL_SCALE:
1072:   case MATOP_SHIFT:
1073:   case MATOP_SCALE:
1074:   case MATOP_AXPY:
1075:     if (shell->managescalingshifts) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"MATSHELL is managing scalings and shifts, see MatShellSetManageScalingShifts()");
1076:     (((void(**)(void))mat->ops)[op]) = f;
1077:     break;
1078:   case MATOP_GET_DIAGONAL:
1079:     if (shell->managescalingshifts) {
1080:       shell->ops->getdiagonal = (PetscErrorCode (*)(Mat,Vec))f;
1081:       mat->ops->getdiagonal   = MatGetDiagonal_Shell;
1082:     } else {
1083:       shell->ops->getdiagonal = NULL;
1084:       mat->ops->getdiagonal   = (PetscErrorCode (*)(Mat,Vec))f;
1085:     }
1086:     break;
1087:   case MATOP_MULT:
1088:     if (shell->managescalingshifts) {
1089:       shell->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1090:       mat->ops->mult   = MatMult_Shell;
1091:     } else {
1092:       shell->ops->mult = NULL;
1093:       mat->ops->mult   = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1094:     }
1095:     break;
1096:   case MATOP_MULT_TRANSPOSE:
1097:     if (shell->managescalingshifts) {
1098:       shell->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1099:       mat->ops->multtranspose   = MatMultTranspose_Shell;
1100:     } else {
1101:       shell->ops->multtranspose = NULL;
1102:       mat->ops->multtranspose   = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1103:     }
1104:     break;
1105:   default:
1106:     (((void(**)(void))mat->ops)[op]) = f;
1107:     break;
1108:   }
1109:   return(0);
1110: }

1112: /*@C
1113:     MatShellGetOperation - Gets a matrix function for a shell matrix.

1115:     Not Collective

1117:     Input Parameters:
1118: +   mat - the shell matrix
1119: -   op - the name of the operation

1121:     Output Parameter:
1122: .   f - the function that provides the operation.

1124:     Level: advanced

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

1132:     All user-provided functions have the same calling
1133:     sequence as the usual matrix interface routines, since they
1134:     are intended to be accessed via the usual matrix interface
1135:     routines, e.g.,
1136: $       MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec)

1138:     Within each user-defined routine, the user should call
1139:     MatShellGetContext() to obtain the user-defined context that was
1140:     set by MatCreateShell().

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

1144: .seealso: MatCreateShell(), MatShellGetContext(), MatShellSetOperation(), MatShellSetContext()
1145: @*/
1146: PetscErrorCode MatShellGetOperation(Mat mat,MatOperation op,void(**f)(void))
1147: {
1148:   PetscBool      flg;
1149:   Mat_Shell      *shell;

1154:   PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
1155:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Can only use with MATSHELL matrices");
1156:   shell = (Mat_Shell*)mat->data;

1158:   switch (op) {
1159:   case MATOP_DESTROY:
1160:     *f = (void (*)(void))shell->ops->destroy;
1161:     break;
1162:   case MATOP_VIEW:
1163:     *f = (void (*)(void))mat->ops->view;
1164:     break;
1165:   case MATOP_COPY:
1166:     *f = (void (*)(void))shell->ops->copy;
1167:     break;
1168:   case MATOP_DIAGONAL_SET:
1169:   case MATOP_DIAGONAL_SCALE:
1170:   case MATOP_SHIFT:
1171:   case MATOP_SCALE:
1172:   case MATOP_AXPY:
1173:     *f = (((void (**)(void))mat->ops)[op]);
1174:     break;
1175:   case MATOP_GET_DIAGONAL:
1176:     if (shell->ops->getdiagonal)
1177:       *f = (void (*)(void))shell->ops->getdiagonal;
1178:     else
1179:       *f = (((void (**)(void))mat->ops)[op]);
1180:     break;
1181:   case MATOP_MULT:
1182:     if (shell->ops->mult)
1183:       *f = (void (*)(void))shell->ops->mult;
1184:     else
1185:       *f = (((void (**)(void))mat->ops)[op]);
1186:     break;
1187:   case MATOP_MULT_TRANSPOSE:
1188:     if (shell->ops->multtranspose)
1189:       *f = (void (*)(void))shell->ops->multtranspose;
1190:     else
1191:       *f = (((void (**)(void))mat->ops)[op]);
1192:     break;
1193:   default:
1194:     *f = (((void (**)(void))mat->ops)[op]);
1195:   }
1196:   return(0);
1197: }