Actual source code: shell.c

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

8: #include <petsc-private/matimpl.h>        /*I "petscmat.h" I*/
9: #include <petsc-private/vecimpl.h>

11: typedef struct {
12:   PetscErrorCode (*destroy)(Mat);
13:   PetscErrorCode (*mult)(Mat,Vec,Vec);
14:   PetscErrorCode (*multtranspose)(Mat,Vec,Vec);
15:   PetscErrorCode (*getdiagonal)(Mat,Vec);
16:   PetscScalar    vscale,vshift;
17:   Vec            dshift;
18:   Vec            left,right;
19:   Vec            dshift_owned,left_owned,right_owned;
20:   Vec            left_work,right_work;
21:   PetscBool      usingscaled;
22:   void           *ctx;
23: } Mat_Shell;
24: /*
25:  The most general expression for the matrix is

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

29:  where
30:  A is the matrix defined by the user's function
31:  a is a scalar multiple
32:  L is left scaling
33:  R is right scaling
34:  B is a diagonal shift defined by
35:    diag(dshift) if the vector dshift is non-NULL
36:    vscale*identity otherwise

38:  The following identities apply:

40:  Scale by c:
41:   c [L (a A + B) R] = L [(a c) A + c B] R

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

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

48:  In the data structure:

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

55: static PetscErrorCode MatMult_Shell(Mat,Vec,Vec);
56: static PetscErrorCode MatMultTranspose_Shell(Mat,Vec,Vec);
57: static PetscErrorCode MatGetDiagonal_Shell(Mat,Vec);

61: static PetscErrorCode MatShellUseScaledMethods(Mat Y)
62: {
63:   Mat_Shell *shell = (Mat_Shell*)Y->data;

66:   if (shell->usingscaled) return(0);
67:   shell->mult  = Y->ops->mult;
68:   Y->ops->mult = MatMult_Shell;
69:   if (Y->ops->multtranspose) {
70:     shell->multtranspose  = Y->ops->multtranspose;
71:     Y->ops->multtranspose = MatMultTranspose_Shell;
72:   }
73:   if (Y->ops->getdiagonal) {
74:     shell->getdiagonal  = Y->ops->getdiagonal;
75:     Y->ops->getdiagonal = MatGetDiagonal_Shell;
76:   }
77:   shell->usingscaled = PETSC_TRUE;
78:   return(0);
79: }

83: static PetscErrorCode MatShellPreScaleLeft(Mat A,Vec x,Vec *xx)
84: {
85:   Mat_Shell *shell = (Mat_Shell*)A->data;

89:   *xx = PETSC_NULL;
90:   if (!shell->left) {
91:     *xx = x;
92:   } else {
93:     if (!shell->left_work) {VecDuplicate(shell->left,&shell->left_work);}
94:     VecPointwiseMult(shell->left_work,x,shell->left);
95:     *xx = shell->left_work;
96:   }
97:   return(0);
98: }

102: static PetscErrorCode MatShellPreScaleRight(Mat A,Vec x,Vec *xx)
103: {
104:   Mat_Shell *shell = (Mat_Shell*)A->data;

108:   *xx = PETSC_NULL;
109:   if (!shell->right) {
110:     *xx = x;
111:   } else {
112:     if (!shell->right_work) {VecDuplicate(shell->right,&shell->right_work);}
113:     VecPointwiseMult(shell->right_work,x,shell->right);
114:     *xx = shell->right_work;
115:   }
116:   return(0);
117: }

121: static PetscErrorCode MatShellPostScaleLeft(Mat A,Vec x)
122: {
123:   Mat_Shell *shell = (Mat_Shell*)A->data;

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

133: static PetscErrorCode MatShellPostScaleRight(Mat A,Vec x)
134: {
135:   Mat_Shell *shell = (Mat_Shell*)A->data;

139:   if (shell->right) {VecPointwiseMult(x,x,shell->right);}
140:   return(0);
141: }

145: static PetscErrorCode MatShellShiftAndScale(Mat A,Vec X,Vec Y)
146: {
147:   Mat_Shell *shell = (Mat_Shell*)A->data;

151:   if (shell->dshift) {          /* get arrays because there is no VecPointwiseMultAdd() */
152:     PetscInt          i,m;
153:     const PetscScalar *x,*d;
154:     PetscScalar       *y;
155:     VecGetLocalSize(X,&m);
158:     VecGetArray(Y,&y);
159:     for (i=0; i<m; i++) y[i] = shell->vscale*y[i] + d[i]*x[i];
162:     VecRestoreArray(Y,&y);
163:   } else if (PetscAbsScalar(shell->vshift) != 0) {
164:     VecAXPBY(Y,shell->vshift,shell->vscale,X);
165:   } else if (shell->vscale != 1.0) {
166:     VecScale(Y,shell->vscale);
167:   }
168:   return(0);
169: }

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

176:     Not Collective

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

181:     Output Parameter:
182: .   ctx - the user provided context

186:     Notes:
187:     This routine is intended for use within various shell matrix routines,
188:     as set with MatShellSetOperation().
189:
190: .keywords: matrix, shell, get, context

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

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

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

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

230: PetscErrorCode MatMult_Shell(Mat A,Vec x,Vec y)
231: {
232:   Mat_Shell      *shell = (Mat_Shell*)A->data;
234:   Vec            xx;

237:   MatShellPreScaleRight(A,x,&xx);
238:   (*shell->mult)(A,xx,y);
239:   MatShellShiftAndScale(A,xx,y);
240:   MatShellPostScaleLeft(A,y);
241:   return(0);
242: }

246: PetscErrorCode MatMultTranspose_Shell(Mat A,Vec x,Vec y)
247: {
248:   Mat_Shell      *shell = (Mat_Shell*)A->data;
250:   Vec            xx;

253:   MatShellPreScaleLeft(A,x,&xx);
254:   (*shell->multtranspose)(A,xx,y);
255:   MatShellShiftAndScale(A,xx,y);
256:   MatShellPostScaleRight(A,y);
257:   return(0);
258: }

262: PetscErrorCode MatGetDiagonal_Shell(Mat A,Vec v)
263: {
264:   Mat_Shell      *shell = (Mat_Shell*)A->data;

268:   (*shell->getdiagonal)(A,v);
269:   VecScale(v,shell->vscale);
270:   if (shell->dshift) {
271:     VecPointwiseMult(v,v,shell->dshift);
272:   } else {
273:     VecShift(v,shell->vshift);
274:   }
275:   if (shell->left)  {VecPointwiseMult(v,v,shell->left);}
276:   if (shell->right) {VecPointwiseMult(v,v,shell->right);}
277:   return(0);
278: }

282: PetscErrorCode MatShift_Shell(Mat Y,PetscScalar a)
283: {
284:   Mat_Shell *shell = (Mat_Shell*)Y->data;

288:   if (shell->left || shell->right || shell->dshift) {
289:     if (!shell->dshift) {
290:       if (!shell->dshift_owned) {VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift_owned);}
291:       shell->dshift = shell->dshift_owned;
292:       VecSet(shell->dshift,shell->vshift+a);
293:     } else {VecScale(shell->dshift,a);}
294:     if (shell->left)  {VecPointwiseDivide(shell->dshift,shell->dshift,shell->left);}
295:     if (shell->right) {VecPointwiseDivide(shell->dshift,shell->dshift,shell->right);}
296:   } else shell->vshift += a;
297:   MatShellUseScaledMethods(Y);
298:   return(0);
299: }

303: PetscErrorCode MatScale_Shell(Mat Y,PetscScalar a)
304: {
305:   Mat_Shell      *shell = (Mat_Shell*)Y->data;

309:   shell->vscale *= a;
310:   if (shell->dshift) {VecScale(shell->dshift,a);}
311:   else shell->vshift *= a;
312:   MatShellUseScaledMethods(Y);
313:   return(0);
314: }

318: static PetscErrorCode MatDiagonalScale_Shell(Mat Y,Vec left,Vec right)
319: {
320:   Mat_Shell *shell = (Mat_Shell*)Y->data;

324:   if (left) {
325:     if (!shell->left_owned) {VecDuplicate(left,&shell->left_owned);}
326:     if (shell->left) {VecPointwiseMult(shell->left,shell->left,left);}
327:     else {
328:       shell->left = shell->left_owned;
329:       VecCopy(left,shell->left);
330:     }
331:   }
332:   if (right) {
333:     if (!shell->right_owned) {VecDuplicate(right,&shell->right_owned);}
334:     if (shell->right) {VecPointwiseMult(shell->right,shell->right,right);}
335:     else {
336:       shell->right = shell->right_owned;
337:       VecCopy(right,shell->right);
338:     }
339:   }
340:   MatShellUseScaledMethods(Y);
341:   return(0);
342: }

346: PetscErrorCode MatAssemblyEnd_Shell(Mat Y,MatAssemblyType t)
347: {
348:   Mat_Shell *shell = (Mat_Shell*)Y->data;

351:   if (t == MAT_FINAL_ASSEMBLY) {
352:     shell->vshift      = 0.0;
353:     shell->vscale      = 1.0;
354:     shell->dshift      = PETSC_NULL;
355:     shell->left        = PETSC_NULL;
356:     shell->right       = PETSC_NULL;
357:     if (shell->mult) {
358:       Y->ops->mult = shell->mult;
359:       shell->mult  = PETSC_NULL;
360:     }
361:     if (shell->multtranspose) {
362:       Y->ops->multtranspose = shell->multtranspose;
363:       shell->multtranspose  = PETSC_NULL;
364:     }
365:     if (shell->getdiagonal) {
366:       Y->ops->getdiagonal = shell->getdiagonal;
367:       shell->getdiagonal  = PETSC_NULL;
368:     }
369:     shell->usingscaled = PETSC_FALSE;
370:   }
371:   return(0);
372: }

374: extern PetscErrorCode MatConvert_Shell(Mat, const MatType,MatReuse,Mat*);

376: static struct _MatOps MatOps_Values = {0,
377:        0,
378:        0,
379:        0,
380: /* 4*/ 0,
381:        0,
382:        0,
383:        0,
384:        0,
385:        0,
386: /*10*/ 0,
387:        0,
388:        0,
389:        0,
390:        0,
391: /*15*/ 0,
392:        0,
393:        0,
394:        MatDiagonalScale_Shell,
395:        0,
396: /*20*/ 0,
397:        MatAssemblyEnd_Shell,
398:        0,
399:        0,
400: /*24*/ 0,
401:        0,
402:        0,
403:        0,
404:        0,
405: /*29*/ 0,
406:        0,
407:        0,
408:        0,
409:        0,
410: /*34*/ 0,
411:        0,
412:        0,
413:        0,
414:        0,
415: /*39*/ 0,
416:        0,
417:        0,
418:        0,
419:        0,
420: /*44*/ 0,
421:        MatScale_Shell,
422:        MatShift_Shell,
423:        0,
424:        0,
425: /*49*/ 0,
426:        0,
427:        0,
428:        0,
429:        0,
430: /*54*/ 0,
431:        0,
432:        0,
433:        0,
434:        0,
435: /*59*/ 0,
436:        MatDestroy_Shell,
437:        0,
438:        0,
439:        0,
440: /*64*/ 0,
441:        0,
442:        0,
443:        0,
444:        0,
445: /*69*/ 0,
446:        0,
447:        MatConvert_Shell,
448:        0,
449:        0,
450: /*74*/ 0,
451:        0,
452:        0,
453:        0,
454:        0,
455: /*79*/ 0,
456:        0,
457:        0,
458:        0,
459:        0,
460: /*84*/ 0,
461:        0,
462:        0,
463:        0,
464:        0,
465: /*89*/ 0,
466:        0,
467:        0,
468:        0,
469:        0,
470: /*94*/ 0,
471:        0,
472:        0,
473:        0};

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

480: .seealso: MatCreateShell
481: M*/

483: EXTERN_C_BEGIN
486: PetscErrorCode  MatCreate_Shell(Mat A)
487: {
488:   Mat_Shell      *b;

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

494:   PetscNewLog(A,Mat_Shell,&b);
495:   A->data = (void*)b;

497:   PetscLayoutSetUp(A->rmap);
498:   PetscLayoutSetUp(A->cmap);

500:   b->ctx           = 0;
501:   b->vshift        = 0.0;
502:   b->vscale        = 1.0;
503:   b->mult          = 0;
504:   b->multtranspose = 0;
505:   b->getdiagonal   = 0;
506:   A->assembled     = PETSC_TRUE;
507:   A->preallocated  = PETSC_FALSE;
508:   PetscObjectChangeTypeName((PetscObject)A,MATSHELL);
509:   return(0);
510: }
511: EXTERN_C_END

515: /*@C
516:    MatCreateShell - Creates a new matrix class for use with a user-defined
517:    private data storage format.

519:   Collective on MPI_Comm

521:    Input Parameters:
522: +  comm - MPI communicator
523: .  m - number of local rows (must be given)
524: .  n - number of local columns (must be given)
525: .  M - number of global rows (may be PETSC_DETERMINE)
526: .  N - number of global columns (may be PETSC_DETERMINE)
527: -  ctx - pointer to data needed by the shell matrix routines

529:    Output Parameter:
530: .  A - the matrix

534:   Usage:
535: \$    extern int mult(Mat,Vec,Vec);
536: \$    MatCreateShell(comm,m,n,M,N,ctx,&mat);
537: \$    MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
538: \$    [ Use matrix for operations that have been set ]
539: \$    MatDestroy(mat);

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

546:    Fortran Notes: The context can only be an integer or a PetscObject
547:       unfortunately it cannot be a Fortran array or derived type.

549:    PETSc requires that matrices and vectors being used for certain
550:    operations are partitioned accordingly.  For example, when
551:    creating a shell matrix, A, that supports parallel matrix-vector
552:    products using MatMult(A,x,y) the user should set the number
553:    of local matrix rows to be the number of local elements of the
554:    corresponding result vector, y. Note that this is information is
555:    required for use of the matrix interface routines, even though
556:    the shell matrix may not actually be physically partitioned.
557:    For example,

559: \$
560: \$     Vec x, y
561: \$     extern int mult(Mat,Vec,Vec);
562: \$     Mat A
563: \$
564: \$     VecCreateMPI(comm,PETSC_DECIDE,M,&y);
565: \$     VecCreateMPI(comm,PETSC_DECIDE,N,&x);
566: \$     VecGetLocalSize(y,&m);
567: \$     VecGetLocalSize(x,&n);
568: \$     MatCreateShell(comm,m,n,M,N,ctx,&A);
569: \$     MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
570: \$     MatMult(A,x,y);
571: \$     MatDestroy(A);
572: \$     VecDestroy(y); VecDestroy(x);
573: \$

575: .keywords: matrix, shell, create

577: .seealso: MatShellSetOperation(), MatHasOperation(), MatShellGetContext(), MatShellSetContext()
578: @*/
579: PetscErrorCode  MatCreateShell(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,void *ctx,Mat *A)
580: {

584:   MatCreate(comm,A);
585:   MatSetSizes(*A,m,n,M,N);
586:   MatSetType(*A,MATSHELL);
587:   MatShellSetContext(*A,ctx);
588:   MatSetUp(*A);
589:   return(0);
590: }

594: /*@
595:     MatShellSetContext - sets the context for a shell matrix

597:    Logically Collective on Mat

599:     Input Parameters:
600: +   mat - the shell matrix
601: -   ctx - the context

605:    Fortran Notes: The context can only be an integer or a PetscObject
606:       unfortunately it cannot be a Fortran array or derived type.

608: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation()
609: @*/
610: PetscErrorCode  MatShellSetContext(Mat mat,void *ctx)
611: {
612:   Mat_Shell      *shell = (Mat_Shell*)mat->data;
614:   PetscBool      flg;

618:   PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
619:   if (flg) {
620:     shell->ctx = ctx;
621:   } else SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Cannot attach context to non-shell matrix");
622:   return(0);
623: }

627: /*@C
628:     MatShellSetOperation - Allows user to set a matrix operation for
629:                            a shell matrix.

631:    Logically Collective on Mat

633:     Input Parameters:
634: +   mat - the shell matrix
635: .   op - the name of the operation
636: -   f - the function that provides the operation.

640:     Usage:
641: \$      extern PetscErrorCode usermult(Mat,Vec,Vec);
642: \$      MatCreateShell(comm,m,n,M,N,ctx,&A);
643: \$      MatShellSetOperation(A,MATOP_MULT,(void(*)(void))usermult);

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

651:     All user-provided functions (execept for MATOP_DESTROY) should have the same calling
652:     sequence as the usual matrix interface routines, since they
653:     are intended to be accessed via the usual matrix interface
654:     routines, e.g.,
655: \$       MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec)

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

660:     Within each user-defined routine, the user should call
661:     MatShellGetContext() to obtain the user-defined context that was
662:     set by MatCreateShell().

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

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

669: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext()
670: @*/
671: PetscErrorCode  MatShellSetOperation(Mat mat,MatOperation op,void (*f)(void))
672: {
674:   PetscBool      flg;

678:   if (op == MATOP_DESTROY) {
679:     PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
680:     if (flg) {
681:        Mat_Shell *shell = (Mat_Shell*)mat->data;
682:        shell->destroy                 = (PetscErrorCode (*)(Mat)) f;
683:     } else mat->ops->destroy          = (PetscErrorCode (*)(Mat)) f;
684:   }
685:   else if (op == MATOP_VIEW) mat->ops->view  = (PetscErrorCode (*)(Mat,PetscViewer)) f;
686:   else                       (((void(**)(void))mat->ops)[op]) = f;

688:   return(0);
689: }

693: /*@C
694:     MatShellGetOperation - Gets a matrix function for a shell matrix.

696:     Not Collective

698:     Input Parameters:
699: +   mat - the shell matrix
700: -   op - the name of the operation

702:     Output Parameter:
703: .   f - the function that provides the operation.

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

713:     All user-provided functions have the same calling
714:     sequence as the usual matrix interface routines, since they
715:     are intended to be accessed via the usual matrix interface
716:     routines, e.g.,
717: \$       MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec)

719:     Within each user-defined routine, the user should call
720:     MatShellGetContext() to obtain the user-defined context that was
721:     set by MatCreateShell().

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

725: .seealso: MatCreateShell(), MatShellGetContext(), MatShellSetOperation(), MatShellSetContext()
726: @*/
727: PetscErrorCode  MatShellGetOperation(Mat mat,MatOperation op,void(**f)(void))
728: {
730:   PetscBool      flg;

734:   if (op == MATOP_DESTROY) {
735:     PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);
736:     if (flg) {
737:       Mat_Shell *shell = (Mat_Shell*)mat->data;
738:       *f = (void(*)(void))shell->destroy;
739:     } else {
740:       *f = (void(*)(void))mat->ops->destroy;
741:     }
742:   } else if (op == MATOP_VIEW) {
743:     *f = (void(*)(void))mat->ops->view;
744:   } else {
745:     *f = (((void(**)(void))mat->ops)[op]);
746:   }

748:   return(0);
749: }

```