Actual source code: shell.c

      1: #define PETSCMAT_DLL

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

     9:  #include src/mat/matimpl.h
     10:  #include vecimpl.h

     12: typedef struct {
     13:   PetscErrorCode (*destroy)(Mat);
     14:   PetscErrorCode (*mult)(Mat,Vec,Vec);
     15:   PetscErrorCode (*multtranspose)(Mat,Vec,Vec);
     16:   PetscErrorCode (*getdiagonal)(Mat,Vec);
     17:   PetscTruth     scale,shift;
     18:   PetscScalar    vscale,vshift;
     19:   void           *ctx;
     20: } Mat_Shell;

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

     27:     Not Collective

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

     32:     Output Parameter:
     33: .   ctx - the user provided context

     35:     Level: advanced

     37:     Notes:
     38:     This routine is intended for use within various shell matrix routines,
     39:     as set with MatShellSetOperation().
     40:
     41: .keywords: matrix, shell, get, context

     43: .seealso: MatCreateShell(), MatShellSetOperation(), MatShellSetContext()
     44: @*/
     45: PetscErrorCode PETSCMAT_DLLEXPORT MatShellGetContext(Mat mat,void **ctx)
     46: {
     48:   PetscTruth     flg;

     53:   PetscTypeCompare((PetscObject)mat,MATSHELL,&flg);
     54:   if (!flg) *ctx = 0;
     55:   else      *ctx = ((Mat_Shell*)(mat->data))->ctx;
     56:   return(0);
     57: }

     61: PetscErrorCode MatDestroy_Shell(Mat mat)
     62: {
     64:   Mat_Shell      *shell;

     67:   shell = (Mat_Shell*)mat->data;
     68:   if (shell->destroy) {(*shell->destroy)(mat);}
     69:   PetscFree(shell);
     70:   return(0);
     71: }

     75: PetscErrorCode MatMult_Shell(Mat A,Vec x,Vec y)
     76: {
     77:   Mat_Shell      *shell = (Mat_Shell*)A->data;

     81:   (*shell->mult)(A,x,y);
     82:   if (shell->shift && shell->scale) {
     83:     VecAXPBY(y,shell->vshift,shell->vscale,x);
     84:   } else if (shell->scale) {
     85:     VecScale(y,shell->vscale);
     86:   } else {
     87:     VecAXPY(y,shell->vshift,x);
     88:   }
     89:   return(0);
     90: }

     94: PetscErrorCode MatMultTranspose_Shell(Mat A,Vec x,Vec y)
     95: {
     96:   Mat_Shell      *shell = (Mat_Shell*)A->data;

    100:   (*shell->multtranspose)(A,x,y);
    101:   if (shell->shift && shell->scale) {
    102:     VecAXPBY(y,shell->vshift,shell->vscale,x);
    103:   } else if (shell->scale) {
    104:     VecScale(y,shell->vscale);
    105:   } else {
    106:     VecAXPY(y,shell->vshift,x);
    107:   }
    108:   return(0);
    109: }

    113: PetscErrorCode MatGetDiagonal_Shell(Mat A,Vec v)
    114: {
    115:   Mat_Shell      *shell = (Mat_Shell*)A->data;

    119:   (*shell->getdiagonal)(A,v);
    120:   if (shell->scale) {
    121:     VecScale(v,shell->vscale);
    122:   }
    123:   if (shell->shift) {
    124:     VecShift(v,shell->vshift);
    125:   }
    126:   return(0);
    127: }

    131: PetscErrorCode MatShift_Shell(Mat Y,PetscScalar a)
    132: {
    133:   Mat_Shell *shell = (Mat_Shell*)Y->data;

    136:   if (shell->scale || shell->shift) {
    137:     shell->vshift += a;
    138:   } else {
    139:     shell->mult  = Y->ops->mult;
    140:     Y->ops->mult = MatMult_Shell;
    141:     if (Y->ops->multtranspose) {
    142:       shell->multtranspose  = Y->ops->multtranspose;
    143:       Y->ops->multtranspose = MatMultTranspose_Shell;
    144:     }
    145:     if (Y->ops->getdiagonal) {
    146:       shell->getdiagonal  = Y->ops->getdiagonal;
    147:       Y->ops->getdiagonal = MatGetDiagonal_Shell;
    148:     }
    149:     shell->vshift = a;
    150:   }
    151:   shell->shift  =  PETSC_TRUE;
    152:   return(0);
    153: }

    157: PetscErrorCode MatScale_Shell(Mat Y,PetscScalar a)
    158: {
    159:   Mat_Shell *shell = (Mat_Shell*)Y->data;

    162:   if (shell->scale || shell->shift) {
    163:     shell->vscale *= a;
    164:   } else {
    165:     shell->mult  = Y->ops->mult;
    166:     Y->ops->mult = MatMult_Shell;
    167:     if (Y->ops->multtranspose) {
    168:       shell->multtranspose  = Y->ops->multtranspose;
    169:       Y->ops->multtranspose = MatMultTranspose_Shell;
    170:     }
    171:     if (Y->ops->getdiagonal) {
    172:       shell->getdiagonal  = Y->ops->getdiagonal;
    173:       Y->ops->getdiagonal = MatGetDiagonal_Shell;
    174:     }
    175:     shell->vscale = a;
    176:   }
    177:   shell->scale  =  PETSC_TRUE;
    178:   return(0);
    179: }

    183: PetscErrorCode MatAssemblyEnd_Shell(Mat Y,MatAssemblyType t)
    184: {
    185:   Mat_Shell *shell = (Mat_Shell*)Y->data;

    188:   if ((shell->shift || shell->scale) && t == MAT_FINAL_ASSEMBLY) {
    189:     shell->scale  = PETSC_FALSE;
    190:     shell->shift  = PETSC_FALSE;
    191:     shell->vshift = 0.0;
    192:     shell->vscale = 1.0;
    193:     Y->ops->mult          = shell->mult;
    194:     Y->ops->multtranspose = shell->multtranspose;
    195:     Y->ops->getdiagonal   = shell->getdiagonal;
    196:   }
    197:   return(0);
    198: }

    200: EXTERN PetscErrorCode MatConvert_Shell(Mat, MatType,MatReuse,Mat*);

    204: PetscErrorCode MatSetBlockSize_Shell(Mat A,PetscInt bs)
    205: {
    207:   return(0);
    208: }

    210: static struct _MatOps MatOps_Values = {0,
    211:        0,
    212:        0,
    213:        0,
    214: /* 4*/ 0,
    215:        0,
    216:        0,
    217:        0,
    218:        0,
    219:        0,
    220: /*10*/ 0,
    221:        0,
    222:        0,
    223:        0,
    224:        0,
    225: /*15*/ 0,
    226:        0,
    227:        0,
    228:        0,
    229:        0,
    230: /*20*/ 0,
    231:        MatAssemblyEnd_Shell,
    232:        0,
    233:        0,
    234:        0,
    235: /*25*/ 0,
    236:        0,
    237:        0,
    238:        0,
    239:        0,
    240: /*30*/ 0,
    241:        0,
    242:        0,
    243:        0,
    244:        0,
    245: /*35*/ 0,
    246:        0,
    247:        0,
    248:        0,
    249:        0,
    250: /*40*/ 0,
    251:        0,
    252:        0,
    253:        0,
    254:        0,
    255: /*45*/ 0,
    256:        MatScale_Shell,
    257:        MatShift_Shell,
    258:        0,
    259:        0,
    260: /*50*/ MatSetBlockSize_Shell,
    261:        0,
    262:        0,
    263:        0,
    264:        0,
    265: /*55*/ 0,
    266:        0,
    267:        0,
    268:        0,
    269:        0,
    270: /*60*/ 0,
    271:        MatDestroy_Shell,
    272:        0,
    273:        MatGetPetscMaps_Petsc,
    274:        0,
    275: /*65*/ 0,
    276:        0,
    277:        0,
    278:        0,
    279:        0,
    280: /*70*/ 0,
    281:        MatConvert_Shell,
    282:        0,
    283:        0,
    284:        0,
    285: /*75*/ 0,
    286:        0,
    287:        0,
    288:        0,
    289:        0,
    290: /*80*/ 0,
    291:        0,
    292:        0,
    293:        0,
    294:        0,
    295: /*85*/ 0,
    296:        0,
    297:        0,
    298:        0,
    299:        0,
    300: /*90*/ 0,
    301:        0,
    302:        0,
    303:        0,
    304:        0,
    305: /*95*/ 0,
    306:        0,
    307:        0,
    308:        0};

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

    313:   Level: advanced

    315: .seealso: MatCreateShell
    316: M*/

    321: PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_Shell(Mat A)
    322: {
    323:   Mat_Shell      *b;

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

    329:   PetscNew(Mat_Shell,&b);
    330:   PetscLogObjectMemory(A,sizeof(struct _p_Mat)+sizeof(Mat_Shell));
    331:   A->data = (void*)b;

    333:   if (A->m == PETSC_DECIDE || A->n == PETSC_DECIDE) {
    334:     SETERRQ(PETSC_ERR_ARG_WRONG,"Must give local row and column count for matrix");
    335:   }

    337:   PetscSplitOwnership(A->comm,&A->m,&A->M);
    338:   PetscSplitOwnership(A->comm,&A->n,&A->N);

    340:   PetscMapCreateMPI(A->comm,A->m,A->M,&A->rmap);
    341:   PetscMapCreateMPI(A->comm,A->n,A->N,&A->cmap);

    343:   b->ctx           = 0;
    344:   b->scale         = PETSC_FALSE;
    345:   b->shift         = PETSC_FALSE;
    346:   b->vshift        = 0.0;
    347:   b->vscale        = 1.0;
    348:   b->mult          = 0;
    349:   b->multtranspose = 0;
    350:   b->getdiagonal   = 0;
    351:   A->assembled     = PETSC_TRUE;
    352:   A->preallocated  = PETSC_TRUE;
    353:   return(0);
    354: }

    359: /*@C
    360:    MatCreateShell - Creates a new matrix class for use with a user-defined
    361:    private data storage format.

    363:   Collective on MPI_Comm

    365:    Input Parameters:
    366: +  comm - MPI communicator
    367: .  m - number of local rows (must be given)
    368: .  n - number of local columns (must be given)
    369: .  M - number of global rows (may be PETSC_DETERMINE)
    370: .  N - number of global columns (may be PETSC_DETERMINE)
    371: -  ctx - pointer to data needed by the shell matrix routines

    373:    Output Parameter:
    374: .  A - the matrix

    376:    Level: advanced

    378:   Usage:
    380: $    MatCreateShell(comm,m,n,M,N,ctx,&mat);
    381: $    MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
    382: $    [ Use matrix for operations that have been set ]
    383: $    MatDestroy(mat);

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

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

    393:    PETSc requires that matrices and vectors being used for certain
    394:    operations are partitioned accordingly.  For example, when
    395:    creating a shell matrix, A, that supports parallel matrix-vector
    396:    products using MatMult(A,x,y) the user should set the number
    397:    of local matrix rows to be the number of local elements of the
    398:    corresponding result vector, y. Note that this is information is
    399:    required for use of the matrix interface routines, even though
    400:    the shell matrix may not actually be physically partitioned.
    401:    For example,

    403: $
    404: $     Vec x, y
    406: $     Mat A
    407: $
    408: $     VecCreateMPI(comm,PETSC_DECIDE,M,&y);
    409: $     VecCreateMPI(comm,PETSC_DECIDE,N,&x);
    410: $     VecGetLocalSize(y,&m);
    411: $     VecGetLocalSize(x,&n);
    412: $     MatCreateShell(comm,m,n,M,N,ctx,&A);
    413: $     MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
    414: $     MatMult(A,x,y);
    415: $     MatDestroy(A);
    416: $     VecDestroy(y); VecDestroy(x);
    417: $

    419: .keywords: matrix, shell, create

    421: .seealso: MatShellSetOperation(), MatHasOperation(), MatShellGetContext(), MatShellSetContext()
    422: @*/
    423: PetscErrorCode PETSCMAT_DLLEXPORT MatCreateShell(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,void *ctx,Mat *A)
    424: {

    428:   MatCreate(comm,A);
    429:   MatSetSizes(*A,m,n,M,N);
    430:   MatSetType(*A,MATSHELL);
    431:   MatShellSetContext(*A,ctx);
    432:   return(0);
    433: }

    437: /*@C
    438:     MatShellSetContext - sets the context for a shell matrix

    440:    Collective on Mat

    442:     Input Parameters:
    443: +   mat - the shell matrix
    444: -   ctx - the context

    446:    Level: advanced

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

    451: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation()
    452: @*/
    453: PetscErrorCode PETSCMAT_DLLEXPORT MatShellSetContext(Mat mat,void *ctx)
    454: {
    455:   Mat_Shell      *shell = (Mat_Shell*)mat->data;
    457:   PetscTruth     flg;

    461:   PetscTypeCompare((PetscObject)mat,MATSHELL,&flg);
    462:   if (flg) {
    463:     shell->ctx = ctx;
    464:   }
    465:   return(0);
    466: }

    470: /*@C
    471:     MatShellSetOperation - Allows user to set a matrix operation for
    472:                            a shell matrix.

    474:    Collective on Mat

    476:     Input Parameters:
    477: +   mat - the shell matrix
    478: .   op - the name of the operation
    479: -   f - the function that provides the operation.

    481:    Level: advanced

    483:     Usage:
    485: $      MatCreateShell(comm,m,n,M,N,ctx,&A);
    486: $      MatShellSetOperation(A,MATOP_MULT,(void(*)(void))usermult);

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

    494:     All user-provided functions should have the same calling
    495:     sequence as the usual matrix interface routines, since they
    496:     are intended to be accessed via the usual matrix interface
    497:     routines, e.g.,
    498: $       MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec)

    500:     Within each user-defined routine, the user should call
    501:     MatShellGetContext() to obtain the user-defined context that was
    502:     set by MatCreateShell().

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

    506: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext()
    507: @*/
    508: PetscErrorCode PETSCMAT_DLLEXPORT MatShellSetOperation(Mat mat,MatOperation op,void (*f)(void))
    509: {
    511:   PetscTruth     flg;

    515:   if (op == MATOP_DESTROY) {
    516:     PetscTypeCompare((PetscObject)mat,MATSHELL,&flg);
    517:     if (flg) {
    518:        Mat_Shell *shell = (Mat_Shell*)mat->data;
    519:        shell->destroy                 = (PetscErrorCode (*)(Mat)) f;
    520:     } else mat->ops->destroy          = (PetscErrorCode (*)(Mat)) f;
    521:   }
    522:   else if (op == MATOP_VIEW) mat->ops->view  = (PetscErrorCode (*)(Mat,PetscViewer)) f;
    523:   else                       (((void(**)(void))mat->ops)[op]) = f;

    525:   return(0);
    526: }

    530: /*@C
    531:     MatShellGetOperation - Gets a matrix function for a shell matrix.

    533:     Not Collective

    535:     Input Parameters:
    536: +   mat - the shell matrix
    537: -   op - the name of the operation

    539:     Output Parameter:
    540: .   f - the function that provides the operation.

    542:     Level: advanced

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

    550:     All user-provided functions have the same calling
    551:     sequence as the usual matrix interface routines, since they
    552:     are intended to be accessed via the usual matrix interface
    553:     routines, e.g.,
    554: $       MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec)

    556:     Within each user-defined routine, the user should call
    557:     MatShellGetContext() to obtain the user-defined context that was
    558:     set by MatCreateShell().

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

    562: .seealso: MatCreateShell(), MatShellGetContext(), MatShellSetOperation(), MatShellSetContext()
    563: @*/
    564: PetscErrorCode PETSCMAT_DLLEXPORT MatShellGetOperation(Mat mat,MatOperation op,void(**f)(void))
    565: {
    567:   PetscTruth     flg;

    571:   if (op == MATOP_DESTROY) {
    572:     PetscTypeCompare((PetscObject)mat,MATSHELL,&flg);
    573:     if (flg) {
    574:       Mat_Shell *shell = (Mat_Shell*)mat->data;
    575:       *f = (void(*)(void))shell->destroy;
    576:     } else {
    577:       *f = (void(*)(void))mat->ops->destroy;
    578:     }
    579:   } else if (op == MATOP_VIEW) {
    580:     *f = (void(*)(void))mat->ops->view;
    581:   } else {
    582:     *f = (((void(**)(void))mat->ops)[op]);
    583:   }

    585:   return(0);
    586: }