Actual source code: shell.c

petsc-master 2018-04-16
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: To use this from Fortran you must write a Fortran interface definition for this
142:     function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.

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

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

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

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

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

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

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

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

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

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

240:   MatShellGetContext(mat,&ctx);
241:   MatCreateShell(PetscObjectComm((PetscObject)mat),mat->rmap->n,mat->cmap->n,mat->rmap->N,mat->cmap->N,ctx,M);
242:   MatCopy(mat,*M,SAME_NONZERO_PATTERN);
243:   return(0);
244: }

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

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

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

274: PetscErrorCode MatMultAdd_Shell(Mat A,Vec x,Vec y,Vec z)
275: {
276:   Mat_Shell      *shell = (Mat_Shell*)A->data;

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

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

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

313: PetscErrorCode MatMultTransposeAdd_Shell(Mat A,Vec x,Vec y,Vec z)
314: {
315:   Mat_Shell      *shell = (Mat_Shell*)A->data;

319:   if (y == z) {
320:     if (!shell->left_add_work) {VecDuplicate(z,&shell->left_add_work);}
321:     MatMultTranspose(A,x,shell->left_add_work);
322:     VecWAXPY(z,1.0,shell->left_add_work,y);
323:   } else {
324:     MatMultTranspose(A,x,z);
325:     VecAXPY(z,1.0,y);
326:   }
327:   return(0);
328: }

330: /*
331:           diag(left)(vscale*A + diag(dshift) + vshift I)diag(right)
332: */
333: PetscErrorCode MatGetDiagonal_Shell(Mat A,Vec v)
334: {
335:   Mat_Shell      *shell = (Mat_Shell*)A->data;

339:   if (shell->ops->getdiagonal) {
340:     (*shell->ops->getdiagonal)(A,v);
341:   } else {
342:     VecSet(v,0.0);
343:   }
344:   VecScale(v,shell->vscale);
345:   if (shell->dshift) {
346:     VecAXPY(v,1.0,shell->dshift);
347:   }
348:   VecShift(v,shell->vshift);
349:   if (shell->left)  {VecPointwiseMult(v,v,shell->left);}
350:   if (shell->right) {VecPointwiseMult(v,v,shell->right);}
351:   if (shell->axpy) {
352:     if (!shell->left_work) {VecDuplicate(v,&shell->left_work);}
353:     MatGetDiagonal(shell->axpy,shell->left_work);
354:     VecAXPY(v,shell->axpy_vscale,shell->left_work);
355:   }
356:   return(0);
357: }

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

365:   if (shell->left || shell->right) {
366:     if (!shell->dshift) {
367:       VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift);
368:       VecSet(shell->dshift,a);
369:     } else {
370:       if (shell->left)  {VecPointwiseMult(shell->dshift,shell->dshift,shell->left);}
371:       if (shell->right) {VecPointwiseMult(shell->dshift,shell->dshift,shell->right);}
372:       VecShift(shell->dshift,a);
373:     }
374:     if (shell->left)  {VecPointwiseDivide(shell->dshift,shell->dshift,shell->left);}
375:     if (shell->right) {VecPointwiseDivide(shell->dshift,shell->dshift,shell->right);}
376:   } else shell->vshift += a;
377:   return(0);
378: }

380: PetscErrorCode MatDiagonalSet_Shell(Mat A,Vec D,InsertMode ins)
381: {
382:   Mat_Shell      *shell = (Mat_Shell*)A->data;

386:   if (ins == INSERT_VALUES) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE, "Operation not supported with INSERT_VALUES");
387:   if (!shell->dshift) {VecDuplicate(D,&shell->dshift);}
388:   if (shell->left || shell->right) {
389:     if (!shell->right_work) {VecDuplicate(shell->left ? shell->left : shell->right, &shell->right_work);}
390:     if (shell->left && shell->right)  {
391:       VecPointwiseDivide(shell->right_work,D,shell->left);
392:       VecPointwiseDivide(shell->right_work,shell->right_work,shell->right);
393:     } else if (shell->left) {
394:       VecPointwiseDivide(shell->right_work,D,shell->left);
395:     } else {
396:       VecPointwiseDivide(shell->right_work,D,shell->right);
397:     }
398:     if (!shell->dshift) {
399:       VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift);
400:       VecCopy(shell->dshift,shell->right_work);
401:     } else {
402:       VecAXPY(shell->dshift,1.0,shell->right_work);
403:     }
404:   } else {
405:     VecAXPY(shell->dshift,1.0,D);
406:   }
407:   return(0);
408: }

410: PetscErrorCode MatScale_Shell(Mat Y,PetscScalar a)
411: {
412:   Mat_Shell      *shell = (Mat_Shell*)Y->data;

416:   shell->vscale *= a;
417:   shell->vshift *= a;
418:   if (shell->dshift) {
419:     VecScale(shell->dshift,a);
420:   }
421:   shell->axpy_vscale *= a;
422:   return(0);
423: }

425: static PetscErrorCode MatDiagonalScale_Shell(Mat Y,Vec left,Vec right)
426: {
427:   Mat_Shell      *shell = (Mat_Shell*)Y->data;

431:   if (shell->axpy) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_SUP,"Cannot diagonal scale MATSHELL after MatAXPY operation");
432:   if (left) {
433:     if (!shell->left) {
434:       VecDuplicate(left,&shell->left);
435:       VecCopy(left,shell->left);
436:     } else {
437:       VecPointwiseMult(shell->left,shell->left,left);
438:     }
439:   }
440:   if (right) {
441:     if (!shell->right) {
442:       VecDuplicate(right,&shell->right);
443:       VecCopy(right,shell->right);
444:     } else {
445:       VecPointwiseMult(shell->right,shell->right,right);
446:     }
447:   }
448:   return(0);
449: }

451: PetscErrorCode MatAssemblyEnd_Shell(Mat Y,MatAssemblyType t)
452: {
453:   Mat_Shell      *shell = (Mat_Shell*)Y->data;

457:   if (t == MAT_FINAL_ASSEMBLY) {
458:     shell->vshift = 0.0;
459:     shell->vscale = 1.0;
460:     VecDestroy(&shell->dshift);
461:     VecDestroy(&shell->left);
462:     VecDestroy(&shell->right);
463:     MatDestroy(&shell->axpy);
464:   }
465:   return(0);
466: }

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

470: static PetscErrorCode MatMissingDiagonal_Shell(Mat A,PetscBool  *missing,PetscInt *d)
471: {
473:   *missing = PETSC_FALSE;
474:   return(0);
475: }

477: PetscErrorCode MatAXPY_Shell(Mat Y,PetscScalar a,Mat X,MatStructure str)
478: {
479:   Mat_Shell      *shell = (Mat_Shell*)Y->data;

483:   PetscObjectReference((PetscObject)X);
484:   MatDestroy(&shell->axpy);
485:   shell->axpy        = X;
486:   shell->axpy_vscale = a;
487:   return(0);
488: }

490: static struct _MatOps MatOps_Values = {0,
491:                                        0,
492:                                        0,
493:                                        0,
494:                                 /* 4*/ MatMultAdd_Shell,
495:                                        0,
496:                                        MatMultTransposeAdd_Shell,
497:                                        0,
498:                                        0,
499:                                        0,
500:                                 /*10*/ 0,
501:                                        0,
502:                                        0,
503:                                        0,
504:                                        0,
505:                                 /*15*/ 0,
506:                                        0,
507:                                        MatGetDiagonal_Shell,
508:                                        MatDiagonalScale_Shell,
509:                                        0,
510:                                 /*20*/ 0,
511:                                        MatAssemblyEnd_Shell,
512:                                        0,
513:                                        0,
514:                                 /*24*/ 0,
515:                                        0,
516:                                        0,
517:                                        0,
518:                                        0,
519:                                 /*29*/ 0,
520:                                        0,
521:                                        0,
522:                                        0,
523:                                        0,
524:                                 /*34*/ MatDuplicate_Shell,
525:                                        0,
526:                                        0,
527:                                        0,
528:                                        0,
529:                                 /*39*/ MatAXPY_Shell,
530:                                        0,
531:                                        0,
532:                                        0,
533:                                        MatCopy_Shell,
534:                                 /*44*/ 0,
535:                                        MatScale_Shell,
536:                                        MatShift_Shell,
537:                                        MatDiagonalSet_Shell,
538:                                        0,
539:                                 /*49*/ 0,
540:                                        0,
541:                                        0,
542:                                        0,
543:                                        0,
544:                                 /*54*/ 0,
545:                                        0,
546:                                        0,
547:                                        0,
548:                                        0,
549:                                 /*59*/ 0,
550:                                        MatDestroy_Shell,
551:                                        0,
552:                                        0,
553:                                        0,
554:                                 /*64*/ 0,
555:                                        0,
556:                                        0,
557:                                        0,
558:                                        0,
559:                                 /*69*/ 0,
560:                                        0,
561:                                        MatConvert_Shell,
562:                                        0,
563:                                        0,
564:                                 /*74*/ 0,
565:                                        0,
566:                                        0,
567:                                        0,
568:                                        0,
569:                                 /*79*/ 0,
570:                                        0,
571:                                        0,
572:                                        0,
573:                                        0,
574:                                 /*84*/ 0,
575:                                        0,
576:                                        0,
577:                                        0,
578:                                        0,
579:                                 /*89*/ 0,
580:                                        0,
581:                                        0,
582:                                        0,
583:                                        0,
584:                                 /*94*/ 0,
585:                                        0,
586:                                        0,
587:                                        0,
588:                                        0,
589:                                 /*99*/ 0,
590:                                        0,
591:                                        0,
592:                                        0,
593:                                        0,
594:                                /*104*/ 0,
595:                                        0,
596:                                        0,
597:                                        0,
598:                                        0,
599:                                /*109*/ 0,
600:                                        0,
601:                                        0,
602:                                        0,
603:                                        MatMissingDiagonal_Shell,
604:                                /*114*/ 0,
605:                                        0,
606:                                        0,
607:                                        0,
608:                                        0,
609:                                /*119*/ 0,
610:                                        0,
611:                                        0,
612:                                        0,
613:                                        0,
614:                                /*124*/ 0,
615:                                        0,
616:                                        0,
617:                                        0,
618:                                        0,
619:                                /*129*/ 0,
620:                                        0,
621:                                        0,
622:                                        0,
623:                                        0,
624:                                /*134*/ 0,
625:                                        0,
626:                                        0,
627:                                        0,
628:                                        0,
629:                                /*139*/ 0,
630:                                        0,
631:                                        0
632: };

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

637:   Level: advanced

639: .seealso: MatCreateShell()
640: M*/

642: PETSC_EXTERN PetscErrorCode MatCreate_Shell(Mat A)
643: {
644:   Mat_Shell      *b;

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

650:   PetscNewLog(A,&b);
651:   A->data = (void*)b;

653:   PetscLayoutSetUp(A->rmap);
654:   PetscLayoutSetUp(A->cmap);

656:   b->ctx                 = 0;
657:   b->vshift              = 0.0;
658:   b->vscale              = 1.0;
659:   b->managescalingshifts = PETSC_TRUE;
660:   A->assembled           = PETSC_TRUE;
661:   A->preallocated        = PETSC_FALSE;

663:   PetscObjectChangeTypeName((PetscObject)A,MATSHELL);
664:   return(0);
665: }

667: /*@C
668:    MatCreateShell - Creates a new matrix class for use with a user-defined
669:    private data storage format.

671:   Collective on MPI_Comm

673:    Input Parameters:
674: +  comm - MPI communicator
675: .  m - number of local rows (must be given)
676: .  n - number of local columns (must be given)
677: .  M - number of global rows (may be PETSC_DETERMINE)
678: .  N - number of global columns (may be PETSC_DETERMINE)
679: -  ctx - pointer to data needed by the shell matrix routines

681:    Output Parameter:
682: .  A - the matrix

684:    Level: advanced

686:   Usage:
687: $    extern int mult(Mat,Vec,Vec);
688: $    MatCreateShell(comm,m,n,M,N,ctx,&mat);
689: $    MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
690: $    [ Use matrix for operations that have been set ]
691: $    MatDestroy(mat);

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

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

702:    PETSc requires that matrices and vectors being used for certain
703:    operations are partitioned accordingly.  For example, when
704:    creating a shell matrix, A, that supports parallel matrix-vector
705:    products using MatMult(A,x,y) the user should set the number
706:    of local matrix rows to be the number of local elements of the
707:    corresponding result vector, y. Note that this is information is
708:    required for use of the matrix interface routines, even though
709:    the shell matrix may not actually be physically partitioned.
710:    For example,

712: $
713: $     Vec x, y
714: $     extern int mult(Mat,Vec,Vec);
715: $     Mat A
716: $
717: $     VecCreateMPI(comm,PETSC_DECIDE,M,&y);
718: $     VecCreateMPI(comm,PETSC_DECIDE,N,&x);
719: $     VecGetLocalSize(y,&m);
720: $     VecGetLocalSize(x,&n);
721: $     MatCreateShell(comm,m,n,M,N,ctx,&A);
722: $     MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
723: $     MatMult(A,x,y);
724: $     MatDestroy(A);
725: $     VecDestroy(y); VecDestroy(x);
726: $


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


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

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

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

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

743:           A is the user provided function.

745: .keywords: matrix, shell, create

747: .seealso: MatShellSetOperation(), MatHasOperation(), MatShellGetContext(), MatShellSetContext(), MATSHELL, MatShellSetManageScalingShifts()
748: @*/
749: PetscErrorCode  MatCreateShell(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,void *ctx,Mat *A)
750: {

754:   MatCreate(comm,A);
755:   MatSetSizes(*A,m,n,M,N);
756:   MatSetType(*A,MATSHELL);
757:   MatShellSetContext(*A,ctx);
758:   MatSetUp(*A);
759:   return(0);
760: }

762: /*@
763:     MatShellSetContext - sets the context for a shell matrix

765:    Logically Collective on Mat

767:     Input Parameters:
768: +   mat - the shell matrix
769: -   ctx - the context

771:    Level: advanced

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

776: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation()
777: @*/
778: PetscErrorCode  MatShellSetContext(Mat mat,void *ctx)
779: {
780:   Mat_Shell      *shell = (Mat_Shell*)mat->data;
782:   PetscBool      flg;

786:   PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
787:   if (flg) {
788:     shell->ctx = ctx;
789:   } else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot attach context to non-shell matrix");
790:   return(0);
791: }

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

797:    Logically Collective on Mat

799:     Input Parameter:
800: .   mat - the shell matrix

802:   Level: advanced

804: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatShellSetOperation()
805: @*/
806: PetscErrorCode MatShellSetManageScalingShifts(Mat A)
807: {
809:   Mat_Shell      *shell;
810:   PetscBool      flg;

814:   PetscObjectTypeCompare((PetscObject)A,MATSHELL,&flg);
815:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Can only use with MATSHELL matrices");
816:   shell = (Mat_Shell*)A->data;
817:   shell->managescalingshifts = PETSC_FALSE;
818:   A->ops->diagonalset   = NULL;
819:   A->ops->diagonalscale = NULL;
820:   A->ops->scale         = NULL;
821:   A->ops->shift         = NULL;
822:   A->ops->axpy          = NULL;
823:   return(0);
824: }

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

829:    Logically Collective on Mat

831:     Input Parameters:
832: +   mat - the shell matrix
833: .   f - the function
834: .   base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated
835: -   ctx - an optional context for the function

837:    Output Parameter:
838: .   flg - PETSC_TRUE if the multiply is likely correct

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

843:    Level: advanced

845:    Fortran Notes: Not supported from Fortran

847: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellTestMultTranspose()
848: @*/
849: PetscErrorCode  MatShellTestMult(Mat mat,PetscErrorCode (*f)(void*,Vec,Vec),Vec base,void *ctx,PetscBool *flg)
850: {
852:   PetscInt       m,n;
853:   Mat            mf,Dmf,Dmat,Ddiff;
854:   PetscReal      Diffnorm,Dmfnorm;
855:   PetscBool      v = PETSC_FALSE, flag = PETSC_TRUE;

859:   PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_view",&v);
860:   MatGetLocalSize(mat,&m,&n);
861:   MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,PETSC_DECIDE,PETSC_DECIDE,&mf);
862:   MatMFFDSetFunction(mf,f,ctx);
863:   MatMFFDSetBase(mf,base,NULL);

865:   MatComputeExplicitOperator(mf,&Dmf);
866:   MatComputeExplicitOperator(mat,&Dmat);

868:   MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff);
869:   MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN);
870:   MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm);
871:   MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm);
872:   if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) {
873:     flag = PETSC_FALSE;
874:     if (v) {
875:       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));
876:       MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_view");
877:       MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_view");
878:       MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_view");
879:     }
880:   } else if (v) {
881:     PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL and matrix free multiple appear to produce the same results\n");
882:   }
883:   if (flg) *flg = flag;
884:   MatDestroy(&Ddiff);
885:   MatDestroy(&mf);
886:   MatDestroy(&Dmf);
887:   MatDestroy(&Dmat);
888:   return(0);
889: }

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

894:    Logically Collective on Mat

896:     Input Parameters:
897: +   mat - the shell matrix
898: .   f - the function
899: .   base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated
900: -   ctx - an optional context for the function

902:    Output Parameter:
903: .   flg - PETSC_TRUE if the multiply is likely correct

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

908:    Level: advanced

910:    Fortran Notes: Not supported from Fortran

912: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellTestMult()
913: @*/
914: PetscErrorCode  MatShellTestMultTranspose(Mat mat,PetscErrorCode (*f)(void*,Vec,Vec),Vec base,void *ctx,PetscBool *flg)
915: {
917:   Vec            x,y,z;
918:   PetscInt       m,n,M,N;
919:   Mat            mf,Dmf,Dmat,Ddiff;
920:   PetscReal      Diffnorm,Dmfnorm;
921:   PetscBool      v = PETSC_FALSE, flag = PETSC_TRUE;

925:   PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_transpose_view",&v);
926:   MatCreateVecs(mat,&x,&y);
927:   VecDuplicate(y,&z);
928:   MatGetLocalSize(mat,&m,&n);
929:   MatGetSize(mat,&M,&N);
930:   MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,M,N,&mf);
931:   MatMFFDSetFunction(mf,f,ctx);
932:   MatMFFDSetBase(mf,base,NULL);
933:   MatComputeExplicitOperator(mf,&Dmf);
934:   MatTranspose(Dmf,MAT_INPLACE_MATRIX,&Dmf);
935:   MatComputeExplicitOperatorTranspose(mat,&Dmat);

937:   MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff);
938:   MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN);
939:   MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm);
940:   MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm);
941:   if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) {
942:     flag = PETSC_FALSE;
943:     if (v) {
944:       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));
945:       MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");
946:       MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");
947:       MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");
948:     }
949:   } else if (v) {
950:     PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL transpose and matrix free multiple appear to produce the same results\n");
951:   }
952:   if (flg) *flg = flag;
953:   MatDestroy(&mf);
954:   MatDestroy(&Dmat);
955:   MatDestroy(&Ddiff);
956:   MatDestroy(&Dmf);
957:   VecDestroy(&x);
958:   VecDestroy(&y);
959:   VecDestroy(&z);
960:   return(0);
961: }

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

966:    Logically Collective on Mat

968:     Input Parameters:
969: +   mat - the shell matrix
970: .   op - the name of the operation
971: -   f - the function that provides the operation.

973:    Level: advanced

975:     Usage:
976: $      extern PetscErrorCode usermult(Mat,Vec,Vec);
977: $      MatCreateShell(comm,m,n,M,N,ctx,&A);
978: $      MatShellSetOperation(A,MATOP_MULT,(void(*)(void))usermult);

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

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

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

995:     Within each user-defined routine, the user should call
996:     MatShellGetContext() to obtain the user-defined context that was
997:     set by MatCreateShell().

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

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

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

1006: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatSetOperation(), MatShellSetManageScalingShifts()
1007: @*/
1008: PetscErrorCode MatShellSetOperation(Mat mat,MatOperation op,void (*f)(void))
1009: {
1010:   PetscBool      flg;
1011:   Mat_Shell      *shell;

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

1020:   switch (op) {
1021:   case MATOP_DESTROY:
1022:     shell->ops->destroy = (PetscErrorCode (*)(Mat))f;
1023:     break;
1024:   case MATOP_VIEW:
1025:     if (!mat->ops->viewnative) {
1026:       mat->ops->viewnative = mat->ops->view;
1027:     }
1028:     mat->ops->view = (PetscErrorCode (*)(Mat,PetscViewer))f;
1029:     break;
1030:   case MATOP_COPY:
1031:     shell->ops->copy = (PetscErrorCode (*)(Mat,Mat,MatStructure))f;
1032:     break;
1033:   case MATOP_DIAGONAL_SET:
1034:   case MATOP_DIAGONAL_SCALE:
1035:   case MATOP_SHIFT:
1036:   case MATOP_SCALE:
1037:   case MATOP_AXPY:
1038:     if (shell->managescalingshifts) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"MATSHELL is managing scalings and shifts, see MatShellSetManageScalingShifts()");
1039:     (((void(**)(void))mat->ops)[op]) = f;
1040:     break;
1041:   case MATOP_GET_DIAGONAL:
1042:     if (shell->managescalingshifts) {
1043:       shell->ops->getdiagonal = (PetscErrorCode (*)(Mat,Vec))f;
1044:       mat->ops->getdiagonal   = MatGetDiagonal_Shell;
1045:     } else {
1046:       shell->ops->getdiagonal = NULL;
1047:       mat->ops->getdiagonal   = (PetscErrorCode (*)(Mat,Vec))f;
1048:     }
1049:     break;
1050:   case MATOP_MULT:
1051:     if (shell->managescalingshifts) {
1052:       shell->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1053:       mat->ops->mult   = MatMult_Shell;
1054:     } else {
1055:       shell->ops->mult = NULL;
1056:       mat->ops->mult   = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1057:     }
1058:     break;
1059:   case MATOP_MULT_TRANSPOSE:
1060:     if (shell->managescalingshifts) {
1061:       shell->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1062:       mat->ops->multtranspose   = MatMultTranspose_Shell;
1063:     } else {
1064:       shell->ops->multtranspose = NULL;
1065:       mat->ops->multtranspose   = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1066:     }
1067:     break;
1068:   default:
1069:     (((void(**)(void))mat->ops)[op]) = f;
1070:     break;
1071:   }
1072:   return(0);
1073: }

1075: /*@C
1076:     MatShellGetOperation - Gets a matrix function for a shell matrix.

1078:     Not Collective

1080:     Input Parameters:
1081: +   mat - the shell matrix
1082: -   op - the name of the operation

1084:     Output Parameter:
1085: .   f - the function that provides the operation.

1087:     Level: advanced

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

1095:     All user-provided functions have the same calling
1096:     sequence as the usual matrix interface routines, since they
1097:     are intended to be accessed via the usual matrix interface
1098:     routines, e.g.,
1099: $       MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec)

1101:     Within each user-defined routine, the user should call
1102:     MatShellGetContext() to obtain the user-defined context that was
1103:     set by MatCreateShell().

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

1107: .seealso: MatCreateShell(), MatShellGetContext(), MatShellSetOperation(), MatShellSetContext()
1108: @*/
1109: PetscErrorCode MatShellGetOperation(Mat mat,MatOperation op,void(**f)(void))
1110: {
1111:   PetscBool      flg;
1112:   Mat_Shell      *shell;

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

1121:   switch (op) {
1122:   case MATOP_DESTROY:
1123:     *f = (void (*)(void))shell->ops->destroy;
1124:     break;
1125:   case MATOP_VIEW:
1126:     *f = (void (*)(void))mat->ops->view;
1127:     break;
1128:   case MATOP_COPY:
1129:     *f = (void (*)(void))shell->ops->copy;
1130:     break;
1131:   case MATOP_DIAGONAL_SET:
1132:   case MATOP_DIAGONAL_SCALE:
1133:   case MATOP_SHIFT:
1134:   case MATOP_SCALE:
1135:   case MATOP_AXPY:
1136:     *f = (((void (**)(void))mat->ops)[op]);
1137:     break;
1138:   case MATOP_GET_DIAGONAL:
1139:     if (shell->ops->getdiagonal)
1140:       *f = (void (*)(void))shell->ops->getdiagonal;
1141:     else
1142:       *f = (((void (**)(void))mat->ops)[op]);
1143:     break;
1144:   case MATOP_MULT:
1145:     if (shell->ops->mult)
1146:       *f = (void (*)(void))shell->ops->mult;
1147:     else
1148:       *f = (((void (**)(void))mat->ops)[op]);
1149:     break;
1150:   case MATOP_MULT_TRANSPOSE:
1151:     if (shell->ops->multtranspose)
1152:       *f = (void (*)(void))shell->ops->multtranspose;
1153:     else
1154:       *f = (((void (**)(void))mat->ops)[op]);
1155:     break;
1156:   default:
1157:     *f = (((void (**)(void))mat->ops)[op]);
1158:   }
1159:   return(0);
1160: }