Actual source code: shell.c

      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>
      9: #include <petsc/private/vecimpl.h>

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

     19: struct _n_MatShellMatFunctionList {
     20:   PetscErrorCode  (*symbolic)(Mat,Mat,Mat,void**);
     21:   PetscErrorCode  (*numeric)(Mat,Mat,Mat,void*);
     22:   PetscErrorCode  (*destroy)(void*);
     23:   MatProductType  ptype;
     24:   char            *composedname;  /* string to identify routine with double dispatch */
     25:   char            *resultname; /* result matrix type */

     27:   struct _n_MatShellMatFunctionList *next;
     28: };
     29: typedef struct _n_MatShellMatFunctionList *MatShellMatFunctionList;

     31: typedef struct {
     32:   struct _MatShellOps ops[1];

     34:   /* The user will manage the scaling and shifts for the MATSHELL, not the default */
     35:   PetscBool managescalingshifts;

     37:   /* support for MatScale, MatShift and MatMultAdd */
     38:   PetscScalar vscale,vshift;
     39:   Vec         dshift;
     40:   Vec         left,right;
     41:   Vec         left_work,right_work;
     42:   Vec         left_add_work,right_add_work;

     44:   /* support for MatAXPY */
     45:   Mat              axpy;
     46:   PetscScalar      axpy_vscale;
     47:   Vec              axpy_left,axpy_right;
     48:   PetscObjectState axpy_state;

     50:   /* support for ZeroRows/Columns operations */
     51:   IS         zrows;
     52:   IS         zcols;
     53:   Vec        zvals;
     54:   Vec        zvals_w;
     55:   VecScatter zvals_sct_r;
     56:   VecScatter zvals_sct_c;

     58:   /* MatMat operations */
     59:   MatShellMatFunctionList matmat;

     61:   /* user defined context */
     62:   void *ctx;
     63: } Mat_Shell;

     66: /*
     67:      Store and scale values on zeroed rows
     68:      xx = [x_1, 0], 0 on zeroed columns
     69: */
     70: static PetscErrorCode MatShellPreZeroRight(Mat A,Vec x,Vec *xx)
     71: {
     72:   Mat_Shell      *shell = (Mat_Shell*)A->data;

     76:   *xx = x;
     77:   if (shell->zrows) {
     78:     VecSet(shell->zvals_w,0.0);
     79:     VecScatterBegin(shell->zvals_sct_c,x,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);
     80:     VecScatterEnd(shell->zvals_sct_c,x,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);
     81:     VecPointwiseMult(shell->zvals_w,shell->zvals_w,shell->zvals);
     82:   }
     83:   if (shell->zcols) {
     84:     if (!shell->right_work) {
     85:       MatCreateVecs(A,&shell->right_work,NULL);
     86:     }
     87:     VecCopy(x,shell->right_work);
     88:     VecISSet(shell->right_work,shell->zcols,0.0);
     89:     *xx  = shell->right_work;
     90:   }
     91:   return(0);
     92: }

     94: /* Insert properly diagonally scaled values stored in MatShellPreZeroRight */
     95: static PetscErrorCode MatShellPostZeroLeft(Mat A,Vec x)
     96: {
     97:   Mat_Shell      *shell = (Mat_Shell*)A->data;

    101:   if (shell->zrows) {
    102:     VecScatterBegin(shell->zvals_sct_r,shell->zvals_w,x,INSERT_VALUES,SCATTER_REVERSE);
    103:     VecScatterEnd(shell->zvals_sct_r,shell->zvals_w,x,INSERT_VALUES,SCATTER_REVERSE);
    104:   }
    105:   return(0);
    106: }

    108: /*
    109:      Store and scale values on zeroed rows
    110:      xx = [x_1, 0], 0 on zeroed rows
    111: */
    112: static PetscErrorCode MatShellPreZeroLeft(Mat A,Vec x,Vec *xx)
    113: {
    114:   Mat_Shell      *shell = (Mat_Shell*)A->data;

    118:   *xx = NULL;
    119:   if (!shell->zrows) {
    120:     *xx = x;
    121:   } else {
    122:     if (!shell->left_work) {
    123:       MatCreateVecs(A,NULL,&shell->left_work);
    124:     }
    125:     VecCopy(x,shell->left_work);
    126:     VecSet(shell->zvals_w,0.0);
    127:     VecScatterBegin(shell->zvals_sct_r,shell->zvals_w,shell->left_work,INSERT_VALUES,SCATTER_REVERSE);
    128:     VecScatterEnd(shell->zvals_sct_r,shell->zvals_w,shell->left_work,INSERT_VALUES,SCATTER_REVERSE);
    129:     VecScatterBegin(shell->zvals_sct_r,x,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);
    130:     VecScatterEnd(shell->zvals_sct_r,x,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);
    131:     VecPointwiseMult(shell->zvals_w,shell->zvals_w,shell->zvals);
    132:     *xx  = shell->left_work;
    133:   }
    134:   return(0);
    135: }

    137: /* Zero zero-columns contributions, sum contributions from properly scaled values stored in MatShellPreZeroLeft */
    138: static PetscErrorCode MatShellPostZeroRight(Mat A,Vec x)
    139: {
    140:   Mat_Shell      *shell = (Mat_Shell*)A->data;

    144:   if (shell->zcols) {
    145:     VecISSet(x,shell->zcols,0.0);
    146:   }
    147:   if (shell->zrows) {
    148:     VecScatterBegin(shell->zvals_sct_c,shell->zvals_w,x,ADD_VALUES,SCATTER_REVERSE);
    149:     VecScatterEnd(shell->zvals_sct_c,shell->zvals_w,x,ADD_VALUES,SCATTER_REVERSE);
    150:   }
    151:   return(0);
    152: }

    154: /*
    155:       xx = diag(left)*x
    156: */
    157: static PetscErrorCode MatShellPreScaleLeft(Mat A,Vec x,Vec *xx)
    158: {
    159:   Mat_Shell      *shell = (Mat_Shell*)A->data;

    163:   *xx = NULL;
    164:   if (!shell->left) {
    165:     *xx = x;
    166:   } else {
    167:     if (!shell->left_work) {VecDuplicate(shell->left,&shell->left_work);}
    168:     VecPointwiseMult(shell->left_work,x,shell->left);
    169:     *xx  = shell->left_work;
    170:   }
    171:   return(0);
    172: }

    174: /*
    175:      xx = diag(right)*x
    176: */
    177: static PetscErrorCode MatShellPreScaleRight(Mat A,Vec x,Vec *xx)
    178: {
    179:   Mat_Shell      *shell = (Mat_Shell*)A->data;

    183:   *xx = NULL;
    184:   if (!shell->right) {
    185:     *xx = x;
    186:   } else {
    187:     if (!shell->right_work) {VecDuplicate(shell->right,&shell->right_work);}
    188:     VecPointwiseMult(shell->right_work,x,shell->right);
    189:     *xx  = shell->right_work;
    190:   }
    191:   return(0);
    192: }

    194: /*
    195:     x = diag(left)*x
    196: */
    197: static PetscErrorCode MatShellPostScaleLeft(Mat A,Vec x)
    198: {
    199:   Mat_Shell      *shell = (Mat_Shell*)A->data;

    203:   if (shell->left) {VecPointwiseMult(x,x,shell->left);}
    204:   return(0);
    205: }

    207: /*
    208:     x = diag(right)*x
    209: */
    210: static PetscErrorCode MatShellPostScaleRight(Mat A,Vec x)
    211: {
    212:   Mat_Shell      *shell = (Mat_Shell*)A->data;

    216:   if (shell->right) {VecPointwiseMult(x,x,shell->right);}
    217:   return(0);
    218: }

    220: /*
    221:          Y = vscale*Y + diag(dshift)*X + vshift*X

    223:          On input Y already contains A*x
    224: */
    225: static PetscErrorCode MatShellShiftAndScale(Mat A,Vec X,Vec Y)
    226: {
    227:   Mat_Shell      *shell = (Mat_Shell*)A->data;

    231:   if (shell->dshift) {          /* get arrays because there is no VecPointwiseMultAdd() */
    232:     PetscInt          i,m;
    233:     const PetscScalar *x,*d;
    234:     PetscScalar       *y;
    235:     VecGetLocalSize(X,&m);
    236:     VecGetArrayRead(shell->dshift,&d);
    237:     VecGetArrayRead(X,&x);
    238:     VecGetArray(Y,&y);
    239:     for (i=0; i<m; i++) y[i] = shell->vscale*y[i] + d[i]*x[i];
    240:     VecRestoreArrayRead(shell->dshift,&d);
    241:     VecRestoreArrayRead(X,&x);
    242:     VecRestoreArray(Y,&y);
    243:   } else {
    244:     VecScale(Y,shell->vscale);
    245:   }
    246:   if (shell->vshift != 0.0) {VecAXPY(Y,shell->vshift,X);} /* if test is for non-square matrices */
    247:   return(0);
    248: }

    250: PetscErrorCode MatShellGetContext_Shell(Mat mat,void *ctx)
    251: {
    253:   *(void**)ctx = ((Mat_Shell*)(mat->data))->ctx;
    254:   return(0);
    255: }

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

    260:     Not Collective

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

    265:     Output Parameter:
    266: .   ctx - the user provided context

    268:     Level: advanced

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

    274: .seealso: MatCreateShell(), MatShellSetOperation(), MatShellSetContext()
    275: @*/
    276: PetscErrorCode  MatShellGetContext(Mat mat,void *ctx)
    277: {

    283:   PetscUseMethod(mat,"MatShellGetContext_C",(Mat,void*),(mat,ctx));
    284:   return(0);
    285: }

    287: static PetscErrorCode MatZeroRowsColumns_Local_Shell(Mat mat,PetscInt nr,PetscInt rows[],PetscInt nc,PetscInt cols[],PetscScalar diag,PetscBool rc)
    288: {
    290:   Mat_Shell      *shell = (Mat_Shell*)mat->data;
    291:   Vec            x = NULL,b = NULL;
    292:   IS             is1, is2;
    293:   const PetscInt *ridxs;
    294:   PetscInt       *idxs,*gidxs;
    295:   PetscInt       cum,rst,cst,i;

    298:   if (!shell->zvals) {
    299:     MatCreateVecs(mat,NULL,&shell->zvals);
    300:   }
    301:   if (!shell->zvals_w) {
    302:     VecDuplicate(shell->zvals,&shell->zvals_w);
    303:   }
    304:   MatGetOwnershipRange(mat,&rst,NULL);
    305:   MatGetOwnershipRangeColumn(mat,&cst,NULL);

    307:   /* Expand/create index set of zeroed rows */
    308:   PetscMalloc1(nr,&idxs);
    309:   for (i = 0; i < nr; i++) idxs[i] = rows[i] + rst;
    310:   ISCreateGeneral(PETSC_COMM_SELF,nr,idxs,PETSC_OWN_POINTER,&is1);
    311:   ISSort(is1);
    312:   VecISSet(shell->zvals,is1,diag);
    313:   if (shell->zrows) {
    314:     ISSum(shell->zrows,is1,&is2);
    315:     ISDestroy(&shell->zrows);
    316:     ISDestroy(&is1);
    317:     shell->zrows = is2;
    318:   } else shell->zrows = is1;

    320:   /* Create scatters for diagonal values communications */
    321:   VecScatterDestroy(&shell->zvals_sct_c);
    322:   VecScatterDestroy(&shell->zvals_sct_r);

    324:   /* row scatter: from/to left vector */
    325:   MatCreateVecs(mat,&x,&b);
    326:   VecScatterCreate(b,shell->zrows,shell->zvals_w,shell->zrows,&shell->zvals_sct_r);

    328:   /* col scatter: from right vector to left vector */
    329:   ISGetIndices(shell->zrows,&ridxs);
    330:   ISGetLocalSize(shell->zrows,&nr);
    331:   PetscMalloc1(nr,&gidxs);
    332:   for (i = 0, cum  = 0; i < nr; i++) {
    333:     if (ridxs[i] >= mat->cmap->N) continue;
    334:     gidxs[cum] = ridxs[i];
    335:     cum++;
    336:   }
    337:   ISCreateGeneral(PetscObjectComm((PetscObject)mat),cum,gidxs,PETSC_OWN_POINTER,&is1);
    338:   VecScatterCreate(x,is1,shell->zvals_w,is1,&shell->zvals_sct_c);
    339:   ISDestroy(&is1);
    340:   VecDestroy(&x);
    341:   VecDestroy(&b);

    343:   /* Expand/create index set of zeroed columns */
    344:   if (rc) {
    345:     PetscMalloc1(nc,&idxs);
    346:     for (i = 0; i < nc; i++) idxs[i] = cols[i] + cst;
    347:     ISCreateGeneral(PETSC_COMM_SELF,nc,idxs,PETSC_OWN_POINTER,&is1);
    348:     ISSort(is1);
    349:     if (shell->zcols) {
    350:       ISSum(shell->zcols,is1,&is2);
    351:       ISDestroy(&shell->zcols);
    352:       ISDestroy(&is1);
    353:       shell->zcols = is2;
    354:     } else shell->zcols = is1;
    355:   }
    356:   return(0);
    357: }

    359: static PetscErrorCode MatZeroRows_Shell(Mat mat,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
    360: {
    361:   Mat_Shell      *shell = (Mat_Shell*)mat->data;
    362:   PetscInt       nr, *lrows;

    366:   if (x && b) {
    367:     Vec          xt;
    368:     PetscScalar *vals;
    369:     PetscInt    *gcols,i,st,nl,nc;

    371:     PetscMalloc1(n,&gcols);
    372:     for (i = 0, nc = 0; i < n; i++) if (rows[i] < mat->cmap->N) gcols[nc++] = rows[i];

    374:     MatCreateVecs(mat,&xt,NULL);
    375:     VecCopy(x,xt);
    376:     PetscCalloc1(nc,&vals);
    377:     VecSetValues(xt,nc,gcols,vals,INSERT_VALUES); /* xt = [x1, 0] */
    378:     PetscFree(vals);
    379:     VecAssemblyBegin(xt);
    380:     VecAssemblyEnd(xt);
    381:     VecAYPX(xt,-1.0,x);                           /* xt = [0, x2] */

    383:     VecGetOwnershipRange(xt,&st,NULL);
    384:     VecGetLocalSize(xt,&nl);
    385:     VecGetArray(xt,&vals);
    386:     for (i = 0; i < nl; i++) {
    387:       PetscInt g = i + st;
    388:       if (g > mat->rmap->N) continue;
    389:       if (PetscAbsScalar(vals[i]) == 0.0) continue;
    390:       VecSetValue(b,g,diag*vals[i],INSERT_VALUES);
    391:     }
    392:     VecRestoreArray(xt,&vals);
    393:     VecAssemblyBegin(b);
    394:     VecAssemblyEnd(b);                            /* b  = [b1, x2 * diag] */
    395:     VecDestroy(&xt);
    396:     PetscFree(gcols);
    397:   }
    398:   PetscLayoutMapLocal(mat->rmap,n,rows,&nr,&lrows,NULL);
    399:   MatZeroRowsColumns_Local_Shell(mat,nr,lrows,0,NULL,diag,PETSC_FALSE);
    400:   if (shell->axpy) {
    401:     MatZeroRows(shell->axpy,n,rows,0.0,NULL,NULL);
    402:   }
    403:   PetscFree(lrows);
    404:   return(0);
    405: }

    407: static PetscErrorCode MatZeroRowsColumns_Shell(Mat mat,PetscInt n,const PetscInt rowscols[],PetscScalar diag,Vec x,Vec b)
    408: {
    409:   Mat_Shell      *shell = (Mat_Shell*)mat->data;
    410:   PetscInt       *lrows, *lcols;
    411:   PetscInt       nr, nc;
    412:   PetscBool      congruent;

    416:   if (x && b) {
    417:     Vec          xt, bt;
    418:     PetscScalar *vals;
    419:     PetscInt    *grows,*gcols,i,st,nl;

    421:     PetscMalloc2(n,&grows,n,&gcols);
    422:     for (i = 0, nr = 0; i < n; i++) if (rowscols[i] < mat->rmap->N) grows[nr++] = rowscols[i];
    423:     for (i = 0, nc = 0; i < n; i++) if (rowscols[i] < mat->cmap->N) gcols[nc++] = rowscols[i];
    424:     PetscCalloc1(n,&vals);

    426:     MatCreateVecs(mat,&xt,&bt);
    427:     VecCopy(x,xt);
    428:     VecSetValues(xt,nc,gcols,vals,INSERT_VALUES); /* xt = [x1, 0] */
    429:     VecAssemblyBegin(xt);
    430:     VecAssemblyEnd(xt);
    431:     VecAXPY(xt,-1.0,x);                           /* xt = [0, -x2] */
    432:     MatMult(mat,xt,bt);                           /* bt = [-A12*x2,-A22*x2] */
    433:     VecSetValues(bt,nr,grows,vals,INSERT_VALUES); /* bt = [-A12*x2,0] */
    434:     VecAssemblyBegin(bt);
    435:     VecAssemblyEnd(bt);
    436:     VecAXPY(b,1.0,bt);                            /* b  = [b1 - A12*x2, b2] */
    437:     VecSetValues(bt,nr,grows,vals,INSERT_VALUES); /* b  = [b1 - A12*x2, 0] */
    438:     VecAssemblyBegin(bt);
    439:     VecAssemblyEnd(bt);
    440:     PetscFree(vals);

    442:     VecGetOwnershipRange(xt,&st,NULL);
    443:     VecGetLocalSize(xt,&nl);
    444:     VecGetArray(xt,&vals);
    445:     for (i = 0; i < nl; i++) {
    446:       PetscInt g = i + st;
    447:       if (g > mat->rmap->N) continue;
    448:       if (PetscAbsScalar(vals[i]) == 0.0) continue;
    449:       VecSetValue(b,g,-diag*vals[i],INSERT_VALUES);
    450:     }
    451:     VecRestoreArray(xt,&vals);
    452:     VecAssemblyBegin(b);
    453:     VecAssemblyEnd(b);                            /* b  = [b1 - A12*x2, x2 * diag] */
    454:     VecDestroy(&xt);
    455:     VecDestroy(&bt);
    456:     PetscFree2(grows,gcols);
    457:   }
    458:   PetscLayoutMapLocal(mat->rmap,n,rowscols,&nr,&lrows,NULL);
    459:   MatHasCongruentLayouts(mat,&congruent);
    460:   if (congruent) {
    461:     nc    = nr;
    462:     lcols = lrows;
    463:   } else { /* MatZeroRowsColumns implicitly assumes the rowscols indices are for a square matrix, here we handle a more general case */
    464:     PetscInt i,nt,*t;

    466:     PetscMalloc1(n,&t);
    467:     for (i = 0, nt = 0; i < n; i++) if (rowscols[i] < mat->cmap->N) t[nt++] = rowscols[i];
    468:     PetscLayoutMapLocal(mat->cmap,nt,t,&nc,&lcols,NULL);
    469:     PetscFree(t);
    470:   }
    471:   MatZeroRowsColumns_Local_Shell(mat,nr,lrows,nc,lcols,diag,PETSC_TRUE);
    472:   if (!congruent) {
    473:     PetscFree(lcols);
    474:   }
    475:   PetscFree(lrows);
    476:   if (shell->axpy) {
    477:     MatZeroRowsColumns(shell->axpy,n,rowscols,0.0,NULL,NULL);
    478:   }
    479:   return(0);
    480: }

    482: PetscErrorCode MatDestroy_Shell(Mat mat)
    483: {
    484:   PetscErrorCode          ierr;
    485:   Mat_Shell               *shell = (Mat_Shell*)mat->data;
    486:   MatShellMatFunctionList matmat;

    489:   if (shell->ops->destroy) {
    490:     (*shell->ops->destroy)(mat);
    491:   }
    492:   PetscMemzero(shell->ops,sizeof(struct _MatShellOps));
    493:   VecDestroy(&shell->left);
    494:   VecDestroy(&shell->right);
    495:   VecDestroy(&shell->dshift);
    496:   VecDestroy(&shell->left_work);
    497:   VecDestroy(&shell->right_work);
    498:   VecDestroy(&shell->left_add_work);
    499:   VecDestroy(&shell->right_add_work);
    500:   VecDestroy(&shell->axpy_left);
    501:   VecDestroy(&shell->axpy_right);
    502:   MatDestroy(&shell->axpy);
    503:   VecDestroy(&shell->zvals_w);
    504:   VecDestroy(&shell->zvals);
    505:   VecScatterDestroy(&shell->zvals_sct_c);
    506:   VecScatterDestroy(&shell->zvals_sct_r);
    507:   ISDestroy(&shell->zrows);
    508:   ISDestroy(&shell->zcols);

    510:   matmat = shell->matmat;
    511:   while (matmat) {
    512:     MatShellMatFunctionList next = matmat->next;

    514:     PetscObjectComposeFunction((PetscObject)mat,matmat->composedname,NULL);
    515:     PetscFree(matmat->composedname);
    516:     PetscFree(matmat->resultname);
    517:     PetscFree(matmat);
    518:     matmat = next;
    519:   }
    520:   PetscObjectComposeFunction((PetscObject)mat,"MatShellGetContext_C",NULL);
    521:   PetscObjectComposeFunction((PetscObject)mat,"MatShellSetContext_C",NULL);
    522:   PetscObjectComposeFunction((PetscObject)mat,"MatShellSetVecType_C",NULL);
    523:   PetscObjectComposeFunction((PetscObject)mat,"MatShellSetManageScalingShifts_C",NULL);
    524:   PetscObjectComposeFunction((PetscObject)mat,"MatShellSetOperation_C",NULL);
    525:   PetscObjectComposeFunction((PetscObject)mat,"MatShellGetOperation_C",NULL);
    526:   PetscObjectComposeFunction((PetscObject)mat,"MatShellSetMatProductOperation_C",NULL);
    527:   PetscFree(mat->data);
    528:   return(0);
    529: }

    531: typedef struct {
    532:   PetscErrorCode (*numeric)(Mat,Mat,Mat,void*);
    533:   PetscErrorCode (*destroy)(void*);
    534:   void           *userdata;
    535:   Mat            B;
    536:   Mat            Bt;
    537:   Mat            axpy;
    538: } MatMatDataShell;

    540: static PetscErrorCode DestroyMatMatDataShell(void *data)
    541: {
    542:   MatMatDataShell *mmdata = (MatMatDataShell *)data;

    546:   if (mmdata->destroy) {
    547:     (*mmdata->destroy)(mmdata->userdata);
    548:   }
    549:   MatDestroy(&mmdata->B);
    550:   MatDestroy(&mmdata->Bt);
    551:   MatDestroy(&mmdata->axpy);
    552:   PetscFree(mmdata);
    553:   return(0);
    554: }

    556: static PetscErrorCode MatProductNumeric_Shell_X(Mat D)
    557: {
    558:   PetscErrorCode  ierr;
    559:   Mat_Product     *product;
    560:   Mat             A, B;
    561:   MatMatDataShell *mdata;
    562:   PetscScalar     zero = 0.0;

    565:   MatCheckProduct(D,1);
    566:   product = D->product;
    567:   if (!product->data) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_PLIB,"Product data empty");
    568:   A = product->A;
    569:   B = product->B;
    570:   mdata = (MatMatDataShell*)product->data;
    571:   if (mdata->numeric) {
    572:     Mat_Shell      *shell = (Mat_Shell*)A->data;
    573:     PetscErrorCode (*stashsym)(Mat) = D->ops->productsymbolic;
    574:     PetscErrorCode (*stashnum)(Mat) = D->ops->productnumeric;
    575:     PetscBool      useBmdata = PETSC_FALSE, newB = PETSC_TRUE;

    577:     if (shell->managescalingshifts) {
    578:       if (shell->zcols || shell->zrows) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_SUP,"MatProduct not supported with zeroed rows/columns");
    579:       if (shell->right || shell->left) {
    580:         useBmdata = PETSC_TRUE;
    581:         if (!mdata->B) {
    582:           MatDuplicate(B,MAT_SHARE_NONZERO_PATTERN,&mdata->B);
    583:         } else {
    584:           newB = PETSC_FALSE;
    585:         }
    586:         MatCopy(B,mdata->B,SAME_NONZERO_PATTERN);
    587:       }
    588:       switch (product->type) {
    589:       case MATPRODUCT_AB: /* s L A R B + v L R B + L D R B */
    590:         if (shell->right) {
    591:           MatDiagonalScale(mdata->B,shell->right,NULL);
    592:         }
    593:         break;
    594:       case MATPRODUCT_AtB: /* s R A^t L B + v R L B + R D L B */
    595:         if (shell->left) {
    596:           MatDiagonalScale(mdata->B,shell->left,NULL);
    597:         }
    598:         break;
    599:       case MATPRODUCT_ABt: /* s L A R B^t + v L R B^t + L D R B^t */
    600:         if (shell->right) {
    601:           MatDiagonalScale(mdata->B,NULL,shell->right);
    602:         }
    603:         break;
    604:       case MATPRODUCT_RARt: /* s B L A R B^t + v B L R B^t + B L D R B^t */
    605:         if (shell->right && shell->left) {
    606:           PetscBool flg;

    608:           VecEqual(shell->right,shell->left,&flg);
    609:           if (!flg) SETERRQ3(PetscObjectComm((PetscObject)D),PETSC_ERR_SUP,"MatProductSymbolic type %s not supported for %s and %s matrices because left scaling != from right scaling",MatProductTypes[product->type],((PetscObject)A)->type_name,((PetscObject)B)->type_name);
    610:         }
    611:         if (shell->right) {
    612:           MatDiagonalScale(mdata->B,NULL,shell->right);
    613:         }
    614:         break;
    615:       case MATPRODUCT_PtAP: /* s B^t L A R B + v B^t L R B + B^t L D R B */
    616:         if (shell->right && shell->left) {
    617:           PetscBool flg;

    619:           VecEqual(shell->right,shell->left,&flg);
    620:           if (!flg) SETERRQ3(PetscObjectComm((PetscObject)D),PETSC_ERR_SUP,"MatProductSymbolic type %s not supported for %s and %s matrices because left scaling != from right scaling",MatProductTypes[product->type],((PetscObject)A)->type_name,((PetscObject)B)->type_name);
    621:         }
    622:         if (shell->right) {
    623:           MatDiagonalScale(mdata->B,shell->right,NULL);
    624:         }
    625:         break;
    626:       default: SETERRQ3(PetscObjectComm((PetscObject)D),PETSC_ERR_SUP,"MatProductSymbolic type %s not supported for %s and %s matrices",MatProductTypes[product->type],((PetscObject)A)->type_name,((PetscObject)B)->type_name);
    627:       }
    628:     }
    629:     /* allow the user to call MatMat operations on D */
    630:     D->product = NULL;
    631:     D->ops->productsymbolic = NULL;
    632:     D->ops->productnumeric  = NULL;

    634:     (*mdata->numeric)(A,useBmdata ? mdata->B : B,D,mdata->userdata);

    636:     /* clear any leftover user data and restore D pointers */
    637:     MatProductClear(D);
    638:     D->ops->productsymbolic = stashsym;
    639:     D->ops->productnumeric  = stashnum;
    640:     D->product = product;

    642:     if (shell->managescalingshifts) {
    643:       MatScale(D,shell->vscale);
    644:       switch (product->type) {
    645:       case MATPRODUCT_AB: /* s L A R B + v L R B + L D R B */
    646:       case MATPRODUCT_ABt: /* s L A R B^t + v L R B^t + L D R B^t */
    647:         if (shell->left) {
    648:           MatDiagonalScale(D,shell->left,NULL);
    649:           if (shell->dshift || shell->vshift != zero) {
    650:             if (!shell->left_work) {MatCreateVecs(A,NULL,&shell->left_work);}
    651:             if (shell->dshift) {
    652:               VecCopy(shell->dshift,shell->left_work);
    653:               VecShift(shell->left_work,shell->vshift);
    654:               VecPointwiseMult(shell->left_work,shell->left_work,shell->left);
    655:             } else {
    656:               VecSet(shell->left_work,shell->vshift);
    657:             }
    658:             if (product->type == MATPRODUCT_ABt) {
    659:               MatReuse     reuse = mdata->Bt ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX;
    660:               MatStructure str = mdata->Bt ? SUBSET_NONZERO_PATTERN : DIFFERENT_NONZERO_PATTERN;

    662:               MatTranspose(mdata->B,reuse,&mdata->Bt);
    663:               MatDiagonalScale(mdata->Bt,shell->left_work,NULL);
    664:               MatAXPY(D,1.0,mdata->Bt,str);
    665:             } else {
    666:               MatStructure str = newB ? DIFFERENT_NONZERO_PATTERN : SUBSET_NONZERO_PATTERN;

    668:               MatDiagonalScale(mdata->B,shell->left_work,NULL);
    669:               MatAXPY(D,1.0,mdata->B,str);
    670:             }
    671:           }
    672:         }
    673:         break;
    674:       case MATPRODUCT_AtB: /* s R A^t L B + v R L B + R D L B */
    675:         if (shell->right) {
    676:           MatDiagonalScale(D,shell->right,NULL);
    677:           if (shell->dshift || shell->vshift != zero) {
    678:             MatStructure str = newB ? DIFFERENT_NONZERO_PATTERN : SUBSET_NONZERO_PATTERN;

    680:             if (!shell->right_work) {MatCreateVecs(A,&shell->right_work,NULL);}
    681:             if (shell->dshift) {
    682:               VecCopy(shell->dshift,shell->right_work);
    683:               VecShift(shell->right_work,shell->vshift);
    684:               VecPointwiseMult(shell->right_work,shell->right_work,shell->right);
    685:             } else {
    686:               VecSet(shell->right_work,shell->vshift);
    687:             }
    688:             MatDiagonalScale(mdata->B,shell->right_work,NULL);
    689:             MatAXPY(D,1.0,mdata->B,str);
    690:           }
    691:         }
    692:         break;
    693:       case MATPRODUCT_PtAP: /* s B^t L A R B + v B^t L R B + B^t L D R B */
    694:       case MATPRODUCT_RARt: /* s B L A R B^t + v B L R B^t + B L D R B^t */
    695:         if (shell->dshift || shell->vshift != zero) SETERRQ3(PetscObjectComm((PetscObject)D),PETSC_ERR_SUP,"MatProductSymbolic type %s not supported for %s and %s matrices with diagonal shift",MatProductTypes[product->type],((PetscObject)A)->type_name,((PetscObject)B)->type_name);
    696:         break;
    697:       default: SETERRQ3(PetscObjectComm((PetscObject)D),PETSC_ERR_SUP,"MatProductSymbolic type %s not supported for %s and %s matrices",MatProductTypes[product->type],((PetscObject)A)->type_name,((PetscObject)B)->type_name);
    698:       }
    699:       if (shell->axpy && shell->axpy_vscale != zero) {
    700:         Mat              X;
    701:         PetscObjectState axpy_state;
    702:         MatStructure     str = DIFFERENT_NONZERO_PATTERN; /* not sure it is safe to ever use SUBSET_NONZERO_PATTERN */

    704:         MatShellGetContext(shell->axpy,(void *)&X);
    705:         PetscObjectStateGet((PetscObject)X,&axpy_state);
    706:         if (shell->axpy_state != axpy_state) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Invalid AXPY state: cannot modify the X matrix passed to MatAXPY(Y,a,X,...)");
    707:         if (!mdata->axpy) {
    708:           str  = DIFFERENT_NONZERO_PATTERN;
    709:           MatProductCreate(shell->axpy,B,NULL,&mdata->axpy);
    710:           MatProductSetType(mdata->axpy,product->type);
    711:           MatProductSetFromOptions(mdata->axpy);
    712:           MatProductSymbolic(mdata->axpy);
    713:         } else { /* May be that shell->axpy has changed */
    714:           PetscBool flg;

    716:           MatProductReplaceMats(shell->axpy,B,NULL,mdata->axpy);
    717:           MatHasOperation(mdata->axpy,MATOP_PRODUCTSYMBOLIC,&flg);
    718:           if (!flg) {
    719:             str  = DIFFERENT_NONZERO_PATTERN;
    720:             MatProductSetFromOptions(mdata->axpy);
    721:             MatProductSymbolic(mdata->axpy);
    722:           }
    723:         }
    724:         MatProductNumeric(mdata->axpy);
    725:         MatAXPY(D,shell->axpy_vscale,mdata->axpy,str);
    726:       }
    727:     }
    728:   } else SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_PLIB,"Missing numeric operation");
    729:   return(0);
    730: }

    732: static PetscErrorCode MatProductSymbolic_Shell_X(Mat D)
    733: {
    734:   PetscErrorCode          ierr;
    735:   Mat_Product             *product;
    736:   Mat                     A,B;
    737:   MatShellMatFunctionList matmat;
    738:   Mat_Shell               *shell;
    739:   PetscBool               flg;
    740:   char                    composedname[256];
    741:   MatMatDataShell         *mdata;

    744:   MatCheckProduct(D,1);
    745:   product = D->product;
    746:   if (product->data) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_PLIB,"Product data not empty");
    747:   A = product->A;
    748:   B = product->B;
    749:   shell = (Mat_Shell*)A->data;
    750:   matmat = shell->matmat;
    751:   PetscSNPrintf(composedname,sizeof(composedname),"MatProductSetFromOptions_%s_%s_C",((PetscObject)A)->type_name,((PetscObject)B)->type_name);
    752:   while (matmat) {
    753:     PetscStrcmp(composedname,matmat->composedname,&flg);
    754:     flg  = (PetscBool)(flg && (matmat->ptype == product->type));
    755:     if (flg) break;
    756:     matmat = matmat->next;
    757:   }
    758:   if (!flg) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Composedname \"%s\" for product type %s not found",composedname,MatProductTypes[product->type]);
    759:   switch (product->type) {
    760:   case MATPRODUCT_AB:
    761:     MatSetSizes(D,A->rmap->n,B->cmap->n,A->rmap->N,B->cmap->N);
    762:     break;
    763:   case MATPRODUCT_AtB:
    764:     MatSetSizes(D,A->cmap->n,B->cmap->n,A->cmap->N,B->cmap->N);
    765:     break;
    766:   case MATPRODUCT_ABt:
    767:     MatSetSizes(D,A->rmap->n,B->rmap->n,A->rmap->N,B->rmap->N);
    768:     break;
    769:   case MATPRODUCT_RARt:
    770:     MatSetSizes(D,B->rmap->n,B->rmap->n,B->rmap->N,B->rmap->N);
    771:     break;
    772:   case MATPRODUCT_PtAP:
    773:     MatSetSizes(D,B->cmap->n,B->cmap->n,B->cmap->N,B->cmap->N);
    774:     break;
    775:   default: SETERRQ3(PetscObjectComm((PetscObject)D),PETSC_ERR_SUP,"MatProductSymbolic type %s not supported for %s and %s matrices",MatProductTypes[product->type],((PetscObject)A)->type_name,((PetscObject)B)->type_name);
    776:   }
    777:   /* respect users who passed in a matrix for which resultname is the base type */
    778:   if (matmat->resultname) {
    779:     PetscObjectBaseTypeCompare((PetscObject)D,matmat->resultname,&flg);
    780:     if (!flg) {
    781:       MatSetType(D,matmat->resultname);
    782:     }
    783:   }
    784:   /* If matrix type was not set or different, we need to reset this pointers */
    785:   D->ops->productsymbolic = MatProductSymbolic_Shell_X;
    786:   D->ops->productnumeric  = MatProductNumeric_Shell_X;
    787:   /* attach product data */
    788:   PetscNew(&mdata);
    789:   mdata->numeric = matmat->numeric;
    790:   mdata->destroy = matmat->destroy;
    791:   if (matmat->symbolic) {
    792:     (*matmat->symbolic)(A,B,D,&mdata->userdata);
    793:   } else { /* call general setup if symbolic operation not provided */
    794:     MatSetUp(D);
    795:   }
    796:   if (!D->product) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_COR,"Product disappeared after user symbolic phase");
    797:   if (D->product->data) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_COR,"Product data not empty after user symbolic phase");
    798:   D->product->data = mdata;
    799:   D->product->destroy = DestroyMatMatDataShell;
    800:   /* Be sure to reset these pointers if the user did something unexpected */
    801:   D->ops->productsymbolic = MatProductSymbolic_Shell_X;
    802:   D->ops->productnumeric  = MatProductNumeric_Shell_X;
    803:   return(0);
    804: }

    806: static PetscErrorCode MatProductSetFromOptions_Shell_X(Mat D)
    807: {
    808:   PetscErrorCode          ierr;
    809:   Mat_Product             *product;
    810:   Mat                     A,B;
    811:   MatShellMatFunctionList matmat;
    812:   Mat_Shell               *shell;
    813:   PetscBool               flg;
    814:   char                    composedname[256];

    817:   MatCheckProduct(D,1);
    818:   product = D->product;
    819:   A = product->A;
    820:   B = product->B;
    821:   MatIsShell(A,&flg);
    822:   if (!flg) return(0);
    823:   shell = (Mat_Shell*)A->data;
    824:   matmat = shell->matmat;
    825:   PetscSNPrintf(composedname,sizeof(composedname),"MatProductSetFromOptions_%s_%s_C",((PetscObject)A)->type_name,((PetscObject)B)->type_name);
    826:   while (matmat) {
    827:     PetscStrcmp(composedname,matmat->composedname,&flg);
    828:     flg  = (PetscBool)(flg && (matmat->ptype == product->type));
    829:     if (flg) break;
    830:     matmat = matmat->next;
    831:   }
    832:   if (flg) { D->ops->productsymbolic = MatProductSymbolic_Shell_X; }
    833:   else { PetscInfo2(D,"  symbolic product %s not registered for product type %s\n",composedname,MatProductTypes[product->type]); }
    834:   return(0);
    835: }

    837: static PetscErrorCode MatShellSetMatProductOperation_Private(Mat A,MatProductType ptype,PetscErrorCode (*symbolic)(Mat,Mat,Mat,void**),PetscErrorCode (*numeric)(Mat,Mat,Mat,void*),PetscErrorCode (*destroy)(void*),char *composedname,const char *resultname)
    838: {
    839:   PetscBool               flg;
    840:   PetscErrorCode          ierr;
    841:   Mat_Shell               *shell;
    842:   MatShellMatFunctionList matmat;

    845:   if (!numeric) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Missing numeric routine");
    846:   if (!composedname) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Missing composed name");

    848:   /* add product callback */
    849:   shell = (Mat_Shell*)A->data;
    850:   matmat = shell->matmat;
    851:   if (!matmat) {
    852:     PetscNew(&shell->matmat);
    853:     matmat = shell->matmat;
    854:   } else {
    855:     MatShellMatFunctionList entry = matmat;
    856:     while (entry) {
    857:       PetscStrcmp(composedname,entry->composedname,&flg);
    858:       flg  = (PetscBool)(flg && (entry->ptype == ptype));
    859:       if (flg) goto set;
    860:       matmat = entry;
    861:       entry = entry->next;
    862:     }
    863:     PetscNew(&matmat->next);
    864:     matmat = matmat->next;
    865:   }

    867: set:
    868:   matmat->symbolic = symbolic;
    869:   matmat->numeric  = numeric;
    870:   matmat->destroy  = destroy;
    871:   matmat->ptype    = ptype;
    872:   PetscFree(matmat->composedname);
    873:   PetscFree(matmat->resultname);
    874:   PetscStrallocpy(composedname,&matmat->composedname);
    875:   PetscStrallocpy(resultname,&matmat->resultname);
    876:   PetscInfo3(A,"Composing %s for product type %s with result %s\n",matmat->composedname,MatProductTypes[matmat->ptype],matmat->resultname ? matmat->resultname : "not specified");
    877:   PetscObjectComposeFunction((PetscObject)A,matmat->composedname,MatProductSetFromOptions_Shell_X);
    878:   return(0);
    879: }

    881: /*@C
    882:     MatShellSetMatProductOperation - Allows user to set a matrix matrix operation for a shell matrix.

    884:    Logically Collective on Mat

    886:     Input Parameters:
    887: +   A - the shell matrix
    888: .   ptype - the product type
    889: .   symbolic - the function for the symbolic phase (can be NULL)
    890: .   numeric - the function for the numerical phase
    891: .   destroy - the function for the destruction of the needed data generated during the symbolic phase (can be NULL)
    892: .   Btype - the matrix type for the matrix to be multiplied against
    893: -   Ctype - the matrix type for the result (can be NULL)

    895:    Level: advanced

    897:     Usage:
    898: $      extern PetscErrorCode usersymbolic(Mat,Mat,Mat,void**);
    899: $      extern PetscErrorCode usernumeric(Mat,Mat,Mat,void*);
    900: $      extern PetscErrorCode userdestroy(void*);
    901: $      MatCreateShell(comm,m,n,M,N,ctx,&A);
    902: $      MatShellSetMatProductOperation(A,MATPRODUCT_AB,usersymbolic,usernumeric,userdestroy,MATSEQAIJ,MATDENSE);
    903: $      [ create B of type SEQAIJ etc..]
    904: $      MatProductCreate(A,B,NULL,&C);
    905: $      MatProductSetType(C,MATPRODUCT_AB);
    906: $      MatProductSetFromOptions(C);
    907: $      MatProductSymbolic(C); -> actually runs the user defined symbolic operation
    908: $      MatProductNumeric(C); -> actually runs the user defined numeric operation
    909: $      [ use C = A*B ]

    911:     Notes:
    912:     MATPRODUCT_ABC is not supported yet. Not supported in Fortran.
    913:     If the symbolic phase is not specified, MatSetUp() is called on the result matrix that must have its type set if Ctype is NULL.
    914:     Any additional data needed by the matrix product needs to be returned during the symbolic phase and destroyed with the destroy callback.
    915:     PETSc will take care of calling the user-defined callbacks.
    916:     It is allowed to specify the same callbacks for different Btype matrix types.
    917:     The couple (Btype,ptype) uniquely identifies the operation: the last specified callbacks takes precedence.

    919: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatSetOperation(), MatProductType, MatType, MatSetUp()
    920: @*/
    921: PetscErrorCode MatShellSetMatProductOperation(Mat A,MatProductType ptype,PetscErrorCode (*symbolic)(Mat,Mat,Mat,void**),PetscErrorCode (*numeric)(Mat,Mat,Mat,void*),PetscErrorCode (*destroy)(void *),MatType Btype,MatType Ctype)
    922: {

    928:   if (ptype == MATPRODUCT_ABC) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Not for product type %s",MatProductTypes[ptype]);
    929:   if (!numeric) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Missing numeric routine, argument 4");
    932:   PetscTryMethod(A,"MatShellSetMatProductOperation_C",(Mat,MatProductType,PetscErrorCode(*)(Mat,Mat,Mat,void**),PetscErrorCode(*)(Mat,Mat,Mat,void*),PetscErrorCode(*)(void*),MatType,MatType),(A,ptype,symbolic,numeric,destroy,Btype,Ctype));
    933:   return(0);
    934: }

    936: PetscErrorCode MatShellSetMatProductOperation_Shell(Mat A,MatProductType ptype,PetscErrorCode (*symbolic)(Mat,Mat,Mat,void**),PetscErrorCode (*numeric)(Mat,Mat,Mat,void*),PetscErrorCode (*destroy)(void *),MatType Btype,MatType Ctype)
    937: {
    938:   PetscBool      flg;
    940:   char           composedname[256];
    941:   MatRootName    Bnames = MatRootNameList, Cnames = MatRootNameList;
    942:   PetscMPIInt    size;

    946:   while (Bnames) { /* user passed in the root name */
    947:     PetscStrcmp(Btype,Bnames->rname,&flg);
    948:     if (flg) break;
    949:     Bnames = Bnames->next;
    950:   }
    951:   while (Cnames) { /* user passed in the root name */
    952:     PetscStrcmp(Ctype,Cnames->rname,&flg);
    953:     if (flg) break;
    954:     Cnames = Cnames->next;
    955:   }
    956:   MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);
    957:   Btype = Bnames ? (size > 1 ? Bnames->mname : Bnames->sname) : Btype;
    958:   Ctype = Cnames ? (size > 1 ? Cnames->mname : Cnames->sname) : Ctype;
    959:   PetscSNPrintf(composedname,sizeof(composedname),"MatProductSetFromOptions_%s_%s_C",((PetscObject)A)->type_name,Btype);
    960:   MatShellSetMatProductOperation_Private(A,ptype,symbolic,numeric,destroy,composedname,Ctype);
    961:   return(0);
    962: }

    964: PetscErrorCode MatCopy_Shell(Mat A,Mat B,MatStructure str)
    965: {
    966:   Mat_Shell               *shellA = (Mat_Shell*)A->data,*shellB = (Mat_Shell*)B->data;
    967:   PetscErrorCode          ierr;
    968:   PetscBool               matflg;
    969:   MatShellMatFunctionList matmatA;

    972:   MatIsShell(B,&matflg);
    973:   if (!matflg) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix %s not derived from MATSHELL",((PetscObject)B)->type_name);

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

    978:   if (shellA->ops->copy) {
    979:     (*shellA->ops->copy)(A,B,str);
    980:   }
    981:   shellB->vscale = shellA->vscale;
    982:   shellB->vshift = shellA->vshift;
    983:   if (shellA->dshift) {
    984:     if (!shellB->dshift) {
    985:       VecDuplicate(shellA->dshift,&shellB->dshift);
    986:     }
    987:     VecCopy(shellA->dshift,shellB->dshift);
    988:   } else {
    989:     VecDestroy(&shellB->dshift);
    990:   }
    991:   if (shellA->left) {
    992:     if (!shellB->left) {
    993:       VecDuplicate(shellA->left,&shellB->left);
    994:     }
    995:     VecCopy(shellA->left,shellB->left);
    996:   } else {
    997:     VecDestroy(&shellB->left);
    998:   }
    999:   if (shellA->right) {
    1000:     if (!shellB->right) {
    1001:       VecDuplicate(shellA->right,&shellB->right);
    1002:     }
    1003:     VecCopy(shellA->right,shellB->right);
    1004:   } else {
    1005:     VecDestroy(&shellB->right);
    1006:   }
    1007:   MatDestroy(&shellB->axpy);
    1008:   shellB->axpy_vscale = 0.0;
    1009:   shellB->axpy_state  = 0;
    1010:   if (shellA->axpy) {
    1011:     PetscObjectReference((PetscObject)shellA->axpy);
    1012:     shellB->axpy        = shellA->axpy;
    1013:     shellB->axpy_vscale = shellA->axpy_vscale;
    1014:     shellB->axpy_state  = shellA->axpy_state;
    1015:   }
    1016:   if (shellA->zrows) {
    1017:     ISDuplicate(shellA->zrows,&shellB->zrows);
    1018:     if (shellA->zcols) {
    1019:       ISDuplicate(shellA->zcols,&shellB->zcols);
    1020:     }
    1021:     VecDuplicate(shellA->zvals,&shellB->zvals);
    1022:     VecCopy(shellA->zvals,shellB->zvals);
    1023:     VecDuplicate(shellA->zvals_w,&shellB->zvals_w);
    1024:     PetscObjectReference((PetscObject)shellA->zvals_sct_r);
    1025:     PetscObjectReference((PetscObject)shellA->zvals_sct_c);
    1026:     shellB->zvals_sct_r = shellA->zvals_sct_r;
    1027:     shellB->zvals_sct_c = shellA->zvals_sct_c;
    1028:   }

    1030:   matmatA = shellA->matmat;
    1031:   if (matmatA) {
    1032:     while (matmatA->next) {
    1033:       MatShellSetMatProductOperation_Private(B,matmatA->ptype,matmatA->symbolic,matmatA->numeric,matmatA->destroy,matmatA->composedname,matmatA->resultname);
    1034:       matmatA = matmatA->next;
    1035:     }
    1036:   }
    1037:   return(0);
    1038: }

    1040: PetscErrorCode MatDuplicate_Shell(Mat mat,MatDuplicateOption op,Mat *M)
    1041: {
    1043:   void           *ctx;

    1046:   MatShellGetContext(mat,&ctx);
    1047:   MatCreateShell(PetscObjectComm((PetscObject)mat),mat->rmap->n,mat->cmap->n,mat->rmap->N,mat->cmap->N,ctx,M);
    1048:   PetscObjectChangeTypeName((PetscObject)(*M),((PetscObject)mat)->type_name);
    1049:   if (op != MAT_DO_NOT_COPY_VALUES) {
    1050:     MatCopy(mat,*M,SAME_NONZERO_PATTERN);
    1051:   }
    1052:   return(0);
    1053: }

    1055: PetscErrorCode MatMult_Shell(Mat A,Vec x,Vec y)
    1056: {
    1057:   Mat_Shell        *shell = (Mat_Shell*)A->data;
    1058:   PetscErrorCode   ierr;
    1059:   Vec              xx;
    1060:   PetscObjectState instate,outstate;

    1063:   if (!shell->ops->mult) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Have not provided a MatMult() for this MATSHELL");
    1064:   MatShellPreZeroRight(A,x,&xx);
    1065:   MatShellPreScaleRight(A,xx,&xx);
    1066:   PetscObjectStateGet((PetscObject)y, &instate);
    1067:   (*shell->ops->mult)(A,xx,y);
    1068:   PetscObjectStateGet((PetscObject)y, &outstate);
    1069:   if (instate == outstate) {
    1070:     /* increase the state of the output vector since the user did not update its state themself as should have been done */
    1071:     PetscObjectStateIncrease((PetscObject)y);
    1072:   }
    1073:   MatShellShiftAndScale(A,xx,y);
    1074:   MatShellPostScaleLeft(A,y);
    1075:   MatShellPostZeroLeft(A,y);

    1077:   if (shell->axpy) {
    1078:     Mat              X;
    1079:     PetscObjectState axpy_state;

    1081:     MatShellGetContext(shell->axpy,(void *)&X);
    1082:     PetscObjectStateGet((PetscObject)X,&axpy_state);
    1083:     if (shell->axpy_state != axpy_state) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Invalid AXPY state: cannot modify the X matrix passed to MatAXPY(Y,a,X,...)");

    1085:     MatCreateVecs(shell->axpy,shell->axpy_right ? NULL : &shell->axpy_right,shell->axpy_left ? NULL : &shell->axpy_left);
    1086:     VecCopy(x,shell->axpy_right);
    1087:     MatMult(shell->axpy,shell->axpy_right,shell->axpy_left);
    1088:     VecAXPY(y,shell->axpy_vscale,shell->axpy_left);
    1089:   }
    1090:   return(0);
    1091: }

    1093: PetscErrorCode MatMultAdd_Shell(Mat A,Vec x,Vec y,Vec z)
    1094: {
    1095:   Mat_Shell      *shell = (Mat_Shell*)A->data;

    1099:   if (y == z) {
    1100:     if (!shell->right_add_work) {VecDuplicate(z,&shell->right_add_work);}
    1101:     MatMult(A,x,shell->right_add_work);
    1102:     VecAXPY(z,1.0,shell->right_add_work);
    1103:   } else {
    1104:     MatMult(A,x,z);
    1105:     VecAXPY(z,1.0,y);
    1106:   }
    1107:   return(0);
    1108: }

    1110: PetscErrorCode MatMultTranspose_Shell(Mat A,Vec x,Vec y)
    1111: {
    1112:   Mat_Shell        *shell = (Mat_Shell*)A->data;
    1113:   PetscErrorCode   ierr;
    1114:   Vec              xx;
    1115:   PetscObjectState instate,outstate;

    1118:   if (!shell->ops->multtranspose) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Have not provided a MatMultTranspose() for this MATSHELL");
    1119:   MatShellPreZeroLeft(A,x,&xx);
    1120:   MatShellPreScaleLeft(A,xx,&xx);
    1121:   PetscObjectStateGet((PetscObject)y, &instate);
    1122:   (*shell->ops->multtranspose)(A,xx,y);
    1123:   PetscObjectStateGet((PetscObject)y, &outstate);
    1124:   if (instate == outstate) {
    1125:     /* increase the state of the output vector since the user did not update its state themself as should have been done */
    1126:     PetscObjectStateIncrease((PetscObject)y);
    1127:   }
    1128:   MatShellShiftAndScale(A,xx,y);
    1129:   MatShellPostScaleRight(A,y);
    1130:   MatShellPostZeroRight(A,y);

    1132:   if (shell->axpy) {
    1133:     Mat              X;
    1134:     PetscObjectState axpy_state;

    1136:     MatShellGetContext(shell->axpy,(void *)&X);
    1137:     PetscObjectStateGet((PetscObject)X,&axpy_state);
    1138:     if (shell->axpy_state != axpy_state) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Invalid AXPY state: cannot modify the X matrix passed to MatAXPY(Y,a,X,...)");
    1139:     MatCreateVecs(shell->axpy,shell->axpy_right ? NULL : &shell->axpy_right,shell->axpy_left ? NULL : &shell->axpy_left);
    1140:     VecCopy(x,shell->axpy_left);
    1141:     MatMultTranspose(shell->axpy,shell->axpy_left,shell->axpy_right);
    1142:     VecAXPY(y,shell->axpy_vscale,shell->axpy_right);
    1143:   }
    1144:   return(0);
    1145: }

    1147: PetscErrorCode MatMultTransposeAdd_Shell(Mat A,Vec x,Vec y,Vec z)
    1148: {
    1149:   Mat_Shell      *shell = (Mat_Shell*)A->data;

    1153:   if (y == z) {
    1154:     if (!shell->left_add_work) {VecDuplicate(z,&shell->left_add_work);}
    1155:     MatMultTranspose(A,x,shell->left_add_work);
    1156:     VecAXPY(z,1.0,shell->left_add_work);
    1157:   } else {
    1158:     MatMultTranspose(A,x,z);
    1159:     VecAXPY(z,1.0,y);
    1160:   }
    1161:   return(0);
    1162: }

    1164: /*
    1165:           diag(left)(vscale*A + diag(dshift) + vshift I)diag(right)
    1166: */
    1167: PetscErrorCode MatGetDiagonal_Shell(Mat A,Vec v)
    1168: {
    1169:   Mat_Shell      *shell = (Mat_Shell*)A->data;

    1173:   if (shell->ops->getdiagonal) {
    1174:     (*shell->ops->getdiagonal)(A,v);
    1175:   } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Must provide shell matrix with routine to return diagonal using\nMatShellSetOperation(S,MATOP_GET_DIAGONAL,...)");
    1176:   VecScale(v,shell->vscale);
    1177:   if (shell->dshift) {
    1178:     VecAXPY(v,1.0,shell->dshift);
    1179:   }
    1180:   VecShift(v,shell->vshift);
    1181:   if (shell->left)  {VecPointwiseMult(v,v,shell->left);}
    1182:   if (shell->right) {VecPointwiseMult(v,v,shell->right);}
    1183:   if (shell->zrows) {
    1184:     VecScatterBegin(shell->zvals_sct_r,shell->zvals,v,INSERT_VALUES,SCATTER_REVERSE);
    1185:     VecScatterEnd(shell->zvals_sct_r,shell->zvals,v,INSERT_VALUES,SCATTER_REVERSE);
    1186:   }
    1187:   if (shell->axpy) {
    1188:     Mat              X;
    1189:     PetscObjectState axpy_state;

    1191:     MatShellGetContext(shell->axpy,(void *)&X);
    1192:     PetscObjectStateGet((PetscObject)X,&axpy_state);
    1193:     if (shell->axpy_state != axpy_state) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Invalid AXPY state: cannot modify the X matrix passed to MatAXPY(Y,a,X,...)");
    1194:     MatCreateVecs(shell->axpy,NULL,shell->axpy_left ? NULL : &shell->axpy_left);
    1195:     MatGetDiagonal(shell->axpy,shell->axpy_left);
    1196:     VecAXPY(v,shell->axpy_vscale,shell->axpy_left);
    1197:   }
    1198:   return(0);
    1199: }

    1201: PetscErrorCode MatShift_Shell(Mat Y,PetscScalar a)
    1202: {
    1203:   Mat_Shell      *shell = (Mat_Shell*)Y->data;
    1205:   PetscBool      flg;

    1208:   MatHasCongruentLayouts(Y,&flg);
    1209:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_SUP,"Cannot shift shell matrix if it is not congruent");
    1210:   if (shell->left || shell->right) {
    1211:     if (!shell->dshift) {
    1212:       VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift);
    1213:       VecSet(shell->dshift,a);
    1214:     } else {
    1215:       if (shell->left)  {VecPointwiseMult(shell->dshift,shell->dshift,shell->left);}
    1216:       if (shell->right) {VecPointwiseMult(shell->dshift,shell->dshift,shell->right);}
    1217:       VecShift(shell->dshift,a);
    1218:     }
    1219:     if (shell->left)  {VecPointwiseDivide(shell->dshift,shell->dshift,shell->left);}
    1220:     if (shell->right) {VecPointwiseDivide(shell->dshift,shell->dshift,shell->right);}
    1221:   } else shell->vshift += a;
    1222:   if (shell->zrows) {
    1223:     VecShift(shell->zvals,a);
    1224:   }
    1225:   return(0);
    1226: }

    1228: PetscErrorCode MatDiagonalSet_Shell_Private(Mat A,Vec D,PetscScalar s)
    1229: {
    1230:   Mat_Shell      *shell = (Mat_Shell*)A->data;

    1234:   if (!shell->dshift) {VecDuplicate(D,&shell->dshift);}
    1235:   if (shell->left || shell->right) {
    1236:     if (!shell->right_work) {VecDuplicate(shell->left ? shell->left : shell->right, &shell->right_work);}
    1237:     if (shell->left && shell->right)  {
    1238:       VecPointwiseDivide(shell->right_work,D,shell->left);
    1239:       VecPointwiseDivide(shell->right_work,shell->right_work,shell->right);
    1240:     } else if (shell->left) {
    1241:       VecPointwiseDivide(shell->right_work,D,shell->left);
    1242:     } else {
    1243:       VecPointwiseDivide(shell->right_work,D,shell->right);
    1244:     }
    1245:     VecAXPY(shell->dshift,s,shell->right_work);
    1246:   } else {
    1247:     VecAXPY(shell->dshift,s,D);
    1248:   }
    1249:   return(0);
    1250: }

    1252: PetscErrorCode MatDiagonalSet_Shell(Mat A,Vec D,InsertMode ins)
    1253: {
    1254:   Mat_Shell      *shell = (Mat_Shell*)A->data;
    1255:   Vec            d;
    1257:   PetscBool      flg;

    1260:   MatHasCongruentLayouts(A,&flg);
    1261:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Cannot diagonal set or shift shell matrix if it is not congruent");
    1262:   if (ins == INSERT_VALUES) {
    1263:     if (!A->ops->getdiagonal) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Operation MATOP_GETDIAGONAL must be set first");
    1264:     VecDuplicate(D,&d);
    1265:     MatGetDiagonal(A,d);
    1266:     MatDiagonalSet_Shell_Private(A,d,-1.);
    1267:     MatDiagonalSet_Shell_Private(A,D,1.);
    1268:     VecDestroy(&d);
    1269:     if (shell->zrows) {
    1270:       VecCopy(D,shell->zvals);
    1271:     }
    1272:   } else {
    1273:     MatDiagonalSet_Shell_Private(A,D,1.);
    1274:     if (shell->zrows) {
    1275:       VecAXPY(shell->zvals,1.0,D);
    1276:     }
    1277:   }
    1278:   return(0);
    1279: }

    1281: PetscErrorCode MatScale_Shell(Mat Y,PetscScalar a)
    1282: {
    1283:   Mat_Shell      *shell = (Mat_Shell*)Y->data;

    1287:   shell->vscale *= a;
    1288:   shell->vshift *= a;
    1289:   if (shell->dshift) {
    1290:     VecScale(shell->dshift,a);
    1291:   }
    1292:   shell->axpy_vscale *= a;
    1293:   if (shell->zrows) {
    1294:     VecScale(shell->zvals,a);
    1295:   }
    1296:   return(0);
    1297: }

    1299: static PetscErrorCode MatDiagonalScale_Shell(Mat Y,Vec left,Vec right)
    1300: {
    1301:   Mat_Shell      *shell = (Mat_Shell*)Y->data;

    1305:   if (left) {
    1306:     if (!shell->left) {
    1307:       VecDuplicate(left,&shell->left);
    1308:       VecCopy(left,shell->left);
    1309:     } else {
    1310:       VecPointwiseMult(shell->left,shell->left,left);
    1311:     }
    1312:     if (shell->zrows) {
    1313:       VecPointwiseMult(shell->zvals,shell->zvals,left);
    1314:     }
    1315:   }
    1316:   if (right) {
    1317:     if (!shell->right) {
    1318:       VecDuplicate(right,&shell->right);
    1319:       VecCopy(right,shell->right);
    1320:     } else {
    1321:       VecPointwiseMult(shell->right,shell->right,right);
    1322:     }
    1323:     if (shell->zrows) {
    1324:       if (!shell->left_work) {
    1325:         MatCreateVecs(Y,NULL,&shell->left_work);
    1326:       }
    1327:       VecSet(shell->zvals_w,1.0);
    1328:       VecScatterBegin(shell->zvals_sct_c,right,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);
    1329:       VecScatterEnd(shell->zvals_sct_c,right,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);
    1330:       VecPointwiseMult(shell->zvals,shell->zvals,shell->zvals_w);
    1331:     }
    1332:   }
    1333:   if (shell->axpy) {
    1334:     MatDiagonalScale(shell->axpy,left,right);
    1335:   }
    1336:   return(0);
    1337: }

    1339: PetscErrorCode MatAssemblyEnd_Shell(Mat Y,MatAssemblyType t)
    1340: {
    1341:   Mat_Shell      *shell = (Mat_Shell*)Y->data;

    1345:   if (t == MAT_FINAL_ASSEMBLY) {
    1346:     shell->vshift = 0.0;
    1347:     shell->vscale = 1.0;
    1348:     shell->axpy_vscale = 0.0;
    1349:     shell->axpy_state  = 0;
    1350:     VecDestroy(&shell->dshift);
    1351:     VecDestroy(&shell->left);
    1352:     VecDestroy(&shell->right);
    1353:     MatDestroy(&shell->axpy);
    1354:     VecDestroy(&shell->axpy_left);
    1355:     VecDestroy(&shell->axpy_right);
    1356:     VecScatterDestroy(&shell->zvals_sct_c);
    1357:     VecScatterDestroy(&shell->zvals_sct_r);
    1358:     ISDestroy(&shell->zrows);
    1359:     ISDestroy(&shell->zcols);
    1360:   }
    1361:   return(0);
    1362: }

    1364: static PetscErrorCode MatMissingDiagonal_Shell(Mat A,PetscBool *missing,PetscInt *d)
    1365: {
    1367:   *missing = PETSC_FALSE;
    1368:   return(0);
    1369: }

    1371: PetscErrorCode MatAXPY_Shell(Mat Y,PetscScalar a,Mat X,MatStructure str)
    1372: {
    1373:   Mat_Shell      *shell = (Mat_Shell*)Y->data;

    1377:   if (X == Y) {
    1378:     MatScale(Y,1.0 + a);
    1379:     return(0);
    1380:   }
    1381:   if (!shell->axpy) {
    1382:     MatConvertFrom_Shell(X,MATSHELL,MAT_INITIAL_MATRIX,&shell->axpy);
    1383:     shell->axpy_vscale = a;
    1384:     PetscObjectStateGet((PetscObject)X,&shell->axpy_state);
    1385:   } else {
    1386:     MatAXPY(shell->axpy,a/shell->axpy_vscale,X,str);
    1387:   }
    1388:   return(0);
    1389: }

    1391: static struct _MatOps MatOps_Values = {NULL,
    1392:                                        NULL,
    1393:                                        NULL,
    1394:                                        NULL,
    1395:                                 /* 4*/ MatMultAdd_Shell,
    1396:                                        NULL,
    1397:                                        MatMultTransposeAdd_Shell,
    1398:                                        NULL,
    1399:                                        NULL,
    1400:                                        NULL,
    1401:                                 /*10*/ NULL,
    1402:                                        NULL,
    1403:                                        NULL,
    1404:                                        NULL,
    1405:                                        NULL,
    1406:                                 /*15*/ NULL,
    1407:                                        NULL,
    1408:                                        NULL,
    1409:                                        MatDiagonalScale_Shell,
    1410:                                        NULL,
    1411:                                 /*20*/ NULL,
    1412:                                        MatAssemblyEnd_Shell,
    1413:                                        NULL,
    1414:                                        NULL,
    1415:                                 /*24*/ MatZeroRows_Shell,
    1416:                                        NULL,
    1417:                                        NULL,
    1418:                                        NULL,
    1419:                                        NULL,
    1420:                                 /*29*/ NULL,
    1421:                                        NULL,
    1422:                                        NULL,
    1423:                                        NULL,
    1424:                                        NULL,
    1425:                                 /*34*/ MatDuplicate_Shell,
    1426:                                        NULL,
    1427:                                        NULL,
    1428:                                        NULL,
    1429:                                        NULL,
    1430:                                 /*39*/ MatAXPY_Shell,
    1431:                                        NULL,
    1432:                                        NULL,
    1433:                                        NULL,
    1434:                                        MatCopy_Shell,
    1435:                                 /*44*/ NULL,
    1436:                                        MatScale_Shell,
    1437:                                        MatShift_Shell,
    1438:                                        MatDiagonalSet_Shell,
    1439:                                        MatZeroRowsColumns_Shell,
    1440:                                 /*49*/ NULL,
    1441:                                        NULL,
    1442:                                        NULL,
    1443:                                        NULL,
    1444:                                        NULL,
    1445:                                 /*54*/ NULL,
    1446:                                        NULL,
    1447:                                        NULL,
    1448:                                        NULL,
    1449:                                        NULL,
    1450:                                 /*59*/ NULL,
    1451:                                        MatDestroy_Shell,
    1452:                                        NULL,
    1453:                                        MatConvertFrom_Shell,
    1454:                                        NULL,
    1455:                                 /*64*/ NULL,
    1456:                                        NULL,
    1457:                                        NULL,
    1458:                                        NULL,
    1459:                                        NULL,
    1460:                                 /*69*/ NULL,
    1461:                                        NULL,
    1462:                                        MatConvert_Shell,
    1463:                                        NULL,
    1464:                                        NULL,
    1465:                                 /*74*/ NULL,
    1466:                                        NULL,
    1467:                                        NULL,
    1468:                                        NULL,
    1469:                                        NULL,
    1470:                                 /*79*/ NULL,
    1471:                                        NULL,
    1472:                                        NULL,
    1473:                                        NULL,
    1474:                                        NULL,
    1475:                                 /*84*/ NULL,
    1476:                                        NULL,
    1477:                                        NULL,
    1478:                                        NULL,
    1479:                                        NULL,
    1480:                                 /*89*/ NULL,
    1481:                                        NULL,
    1482:                                        NULL,
    1483:                                        NULL,
    1484:                                        NULL,
    1485:                                 /*94*/ NULL,
    1486:                                        NULL,
    1487:                                        NULL,
    1488:                                        NULL,
    1489:                                        NULL,
    1490:                                 /*99*/ NULL,
    1491:                                        NULL,
    1492:                                        NULL,
    1493:                                        NULL,
    1494:                                        NULL,
    1495:                                /*104*/ NULL,
    1496:                                        NULL,
    1497:                                        NULL,
    1498:                                        NULL,
    1499:                                        NULL,
    1500:                                /*109*/ NULL,
    1501:                                        NULL,
    1502:                                        NULL,
    1503:                                        NULL,
    1504:                                        MatMissingDiagonal_Shell,
    1505:                                /*114*/ NULL,
    1506:                                        NULL,
    1507:                                        NULL,
    1508:                                        NULL,
    1509:                                        NULL,
    1510:                                /*119*/ NULL,
    1511:                                        NULL,
    1512:                                        NULL,
    1513:                                        NULL,
    1514:                                        NULL,
    1515:                                /*124*/ NULL,
    1516:                                        NULL,
    1517:                                        NULL,
    1518:                                        NULL,
    1519:                                        NULL,
    1520:                                /*129*/ NULL,
    1521:                                        NULL,
    1522:                                        NULL,
    1523:                                        NULL,
    1524:                                        NULL,
    1525:                                /*134*/ NULL,
    1526:                                        NULL,
    1527:                                        NULL,
    1528:                                        NULL,
    1529:                                        NULL,
    1530:                                /*139*/ NULL,
    1531:                                        NULL,
    1532:                                        NULL
    1533: };

    1535: PetscErrorCode  MatShellSetContext_Shell(Mat mat,void *ctx)
    1536: {
    1537:   Mat_Shell *shell = (Mat_Shell*)mat->data;

    1540:   shell->ctx = ctx;
    1541:   return(0);
    1542: }

    1544: static PetscErrorCode MatShellSetVecType_Shell(Mat mat,VecType vtype)
    1545: {

    1549:   PetscFree(mat->defaultvectype);
    1550:   PetscStrallocpy(vtype,(char**)&mat->defaultvectype);
    1551:   return(0);
    1552: }

    1554: PetscErrorCode MatShellSetManageScalingShifts_Shell(Mat A)
    1555: {
    1556:   Mat_Shell      *shell = (Mat_Shell*)A->data;

    1559:   shell->managescalingshifts = PETSC_FALSE;
    1560:   A->ops->diagonalset   = NULL;
    1561:   A->ops->diagonalscale = NULL;
    1562:   A->ops->scale         = NULL;
    1563:   A->ops->shift         = NULL;
    1564:   A->ops->axpy          = NULL;
    1565:   return(0);
    1566: }

    1568: PetscErrorCode MatShellSetOperation_Shell(Mat mat,MatOperation op,void (*f)(void))
    1569: {
    1570:   Mat_Shell      *shell = (Mat_Shell*)mat->data;

    1573:   switch (op) {
    1574:   case MATOP_DESTROY:
    1575:     shell->ops->destroy = (PetscErrorCode (*)(Mat))f;
    1576:     break;
    1577:   case MATOP_VIEW:
    1578:     if (!mat->ops->viewnative) {
    1579:       mat->ops->viewnative = mat->ops->view;
    1580:     }
    1581:     mat->ops->view = (PetscErrorCode (*)(Mat,PetscViewer))f;
    1582:     break;
    1583:   case MATOP_COPY:
    1584:     shell->ops->copy = (PetscErrorCode (*)(Mat,Mat,MatStructure))f;
    1585:     break;
    1586:   case MATOP_DIAGONAL_SET:
    1587:   case MATOP_DIAGONAL_SCALE:
    1588:   case MATOP_SHIFT:
    1589:   case MATOP_SCALE:
    1590:   case MATOP_AXPY:
    1591:   case MATOP_ZERO_ROWS:
    1592:   case MATOP_ZERO_ROWS_COLUMNS:
    1593:     if (shell->managescalingshifts) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"MATSHELL is managing scalings and shifts, see MatShellSetManageScalingShifts()");
    1594:     (((void(**)(void))mat->ops)[op]) = f;
    1595:     break;
    1596:   case MATOP_GET_DIAGONAL:
    1597:     if (shell->managescalingshifts) {
    1598:       shell->ops->getdiagonal = (PetscErrorCode (*)(Mat,Vec))f;
    1599:       mat->ops->getdiagonal   = MatGetDiagonal_Shell;
    1600:     } else {
    1601:       shell->ops->getdiagonal = NULL;
    1602:       mat->ops->getdiagonal   = (PetscErrorCode (*)(Mat,Vec))f;
    1603:     }
    1604:     break;
    1605:   case MATOP_MULT:
    1606:     if (shell->managescalingshifts) {
    1607:       shell->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f;
    1608:       mat->ops->mult   = MatMult_Shell;
    1609:     } else {
    1610:       shell->ops->mult = NULL;
    1611:       mat->ops->mult   = (PetscErrorCode (*)(Mat,Vec,Vec))f;
    1612:     }
    1613:     break;
    1614:   case MATOP_MULT_TRANSPOSE:
    1615:     if (shell->managescalingshifts) {
    1616:       shell->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f;
    1617:       mat->ops->multtranspose   = MatMultTranspose_Shell;
    1618:     } else {
    1619:       shell->ops->multtranspose = NULL;
    1620:       mat->ops->multtranspose   = (PetscErrorCode (*)(Mat,Vec,Vec))f;
    1621:     }
    1622:     break;
    1623:   default:
    1624:     (((void(**)(void))mat->ops)[op]) = f;
    1625:     break;
    1626:   }
    1627:   return(0);
    1628: }

    1630: PetscErrorCode MatShellGetOperation_Shell(Mat mat,MatOperation op,void(**f)(void))
    1631: {
    1632:   Mat_Shell      *shell = (Mat_Shell*)mat->data;

    1635:   switch (op) {
    1636:   case MATOP_DESTROY:
    1637:     *f = (void (*)(void))shell->ops->destroy;
    1638:     break;
    1639:   case MATOP_VIEW:
    1640:     *f = (void (*)(void))mat->ops->view;
    1641:     break;
    1642:   case MATOP_COPY:
    1643:     *f = (void (*)(void))shell->ops->copy;
    1644:     break;
    1645:   case MATOP_DIAGONAL_SET:
    1646:   case MATOP_DIAGONAL_SCALE:
    1647:   case MATOP_SHIFT:
    1648:   case MATOP_SCALE:
    1649:   case MATOP_AXPY:
    1650:   case MATOP_ZERO_ROWS:
    1651:   case MATOP_ZERO_ROWS_COLUMNS:
    1652:     *f = (((void (**)(void))mat->ops)[op]);
    1653:     break;
    1654:   case MATOP_GET_DIAGONAL:
    1655:     if (shell->ops->getdiagonal)
    1656:       *f = (void (*)(void))shell->ops->getdiagonal;
    1657:     else
    1658:       *f = (((void (**)(void))mat->ops)[op]);
    1659:     break;
    1660:   case MATOP_MULT:
    1661:     if (shell->ops->mult)
    1662:       *f = (void (*)(void))shell->ops->mult;
    1663:     else
    1664:       *f = (((void (**)(void))mat->ops)[op]);
    1665:     break;
    1666:   case MATOP_MULT_TRANSPOSE:
    1667:     if (shell->ops->multtranspose)
    1668:       *f = (void (*)(void))shell->ops->multtranspose;
    1669:     else
    1670:       *f = (((void (**)(void))mat->ops)[op]);
    1671:     break;
    1672:   default:
    1673:     *f = (((void (**)(void))mat->ops)[op]);
    1674:   }
    1675:   return(0);
    1676: }

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

    1681:   Level: advanced

    1683: .seealso: MatCreateShell()
    1684: M*/

    1686: PETSC_EXTERN PetscErrorCode MatCreate_Shell(Mat A)
    1687: {
    1688:   Mat_Shell      *b;

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

    1694:   PetscNewLog(A,&b);
    1695:   A->data = (void*)b;

    1697:   b->ctx                 = NULL;
    1698:   b->vshift              = 0.0;
    1699:   b->vscale              = 1.0;
    1700:   b->managescalingshifts = PETSC_TRUE;
    1701:   A->assembled           = PETSC_TRUE;
    1702:   A->preallocated        = PETSC_FALSE;

    1704:   PetscObjectComposeFunction((PetscObject)A,"MatShellGetContext_C",MatShellGetContext_Shell);
    1705:   PetscObjectComposeFunction((PetscObject)A,"MatShellSetContext_C",MatShellSetContext_Shell);
    1706:   PetscObjectComposeFunction((PetscObject)A,"MatShellSetVecType_C",MatShellSetVecType_Shell);
    1707:   PetscObjectComposeFunction((PetscObject)A,"MatShellSetManageScalingShifts_C",MatShellSetManageScalingShifts_Shell);
    1708:   PetscObjectComposeFunction((PetscObject)A,"MatShellSetOperation_C",MatShellSetOperation_Shell);
    1709:   PetscObjectComposeFunction((PetscObject)A,"MatShellGetOperation_C",MatShellGetOperation_Shell);
    1710:   PetscObjectComposeFunction((PetscObject)A,"MatShellSetMatProductOperation_C",MatShellSetMatProductOperation_Shell);
    1711:   PetscObjectChangeTypeName((PetscObject)A,MATSHELL);
    1712:   return(0);
    1713: }

    1715: /*@C
    1716:    MatCreateShell - Creates a new matrix class for use with a user-defined
    1717:    private data storage format.

    1719:   Collective

    1721:    Input Parameters:
    1722: +  comm - MPI communicator
    1723: .  m - number of local rows (must be given)
    1724: .  n - number of local columns (must be given)
    1725: .  M - number of global rows (may be PETSC_DETERMINE)
    1726: .  N - number of global columns (may be PETSC_DETERMINE)
    1727: -  ctx - pointer to data needed by the shell matrix routines

    1729:    Output Parameter:
    1730: .  A - the matrix

    1732:    Level: advanced

    1734:   Usage:
    1735: $    extern PetscErrorCode mult(Mat,Vec,Vec);
    1736: $    MatCreateShell(comm,m,n,M,N,ctx,&mat);
    1737: $    MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
    1738: $    [ Use matrix for operations that have been set ]
    1739: $    MatDestroy(mat);

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

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

    1751:    PETSc requires that matrices and vectors being used for certain
    1752:    operations are partitioned accordingly.  For example, when
    1753:    creating a shell matrix, A, that supports parallel matrix-vector
    1754:    products using MatMult(A,x,y) the user should set the number
    1755:    of local matrix rows to be the number of local elements of the
    1756:    corresponding result vector, y. Note that this is information is
    1757:    required for use of the matrix interface routines, even though
    1758:    the shell matrix may not actually be physically partitioned.
    1759:    For example,

    1761: $
    1762: $     Vec x, y
    1763: $     extern PetscErrorCode mult(Mat,Vec,Vec);
    1764: $     Mat A
    1765: $
    1766: $     VecCreateMPI(comm,PETSC_DECIDE,M,&y);
    1767: $     VecCreateMPI(comm,PETSC_DECIDE,N,&x);
    1768: $     VecGetLocalSize(y,&m);
    1769: $     VecGetLocalSize(x,&n);
    1770: $     MatCreateShell(comm,m,n,M,N,ctx,&A);
    1771: $     MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
    1772: $     MatMult(A,x,y);
    1773: $     MatDestroy(&A);
    1774: $     VecDestroy(&y);
    1775: $     VecDestroy(&x);
    1776: $

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

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

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

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

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

    1794:           A is the user provided function.

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

    1801:    Matrix product operations (i.e. MatMat, MatTranposeMat etc) can be specified using MatShellSetMatProductOperation()

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

    1806: .seealso: MatShellSetOperation(), MatHasOperation(), MatShellGetContext(), MatShellSetContext(), MATSHELL, MatShellSetManageScalingShifts(), MatShellSetMatProductOperation()
    1807: @*/
    1808: PetscErrorCode  MatCreateShell(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,void *ctx,Mat *A)
    1809: {

    1813:   MatCreate(comm,A);
    1814:   MatSetSizes(*A,m,n,M,N);
    1815:   MatSetType(*A,MATSHELL);
    1816:   MatShellSetContext(*A,ctx);
    1817:   MatSetUp(*A);
    1818:   return(0);
    1819: }

    1821: /*@
    1822:     MatShellSetContext - sets the context for a shell matrix

    1824:    Logically Collective on Mat

    1826:     Input Parameters:
    1827: +   mat - the shell matrix
    1828: -   ctx - the context

    1830:    Level: advanced

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

    1836: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation()
    1837: @*/
    1838: PetscErrorCode  MatShellSetContext(Mat mat,void *ctx)
    1839: {

    1844:   PetscTryMethod(mat,"MatShellSetContext_C",(Mat,void*),(mat,ctx));
    1845:   return(0);
    1846: }

    1848: /*@C
    1849:  MatShellSetVecType - Sets the type of Vec returned by MatCreateVecs()

    1851:  Logically collective

    1853:     Input Parameters:
    1854: +   mat   - the shell matrix
    1855: -   vtype - type to use for creating vectors

    1857:  Notes:

    1859:  Level: advanced

    1861: .seealso: MatCreateVecs()
    1862: @*/
    1863: PetscErrorCode  MatShellSetVecType(Mat mat,VecType vtype)
    1864: {

    1868:   PetscTryMethod(mat,"MatShellSetVecType_C",(Mat,VecType),(mat,vtype));
    1869:   return(0);
    1870: }

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

    1876:    Logically Collective on Mat

    1878:     Input Parameter:
    1879: .   mat - the shell matrix

    1881:   Level: advanced

    1883: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatShellSetOperation()
    1884: @*/
    1885: PetscErrorCode MatShellSetManageScalingShifts(Mat A)
    1886: {

    1891:   PetscTryMethod(A,"MatShellSetManageScalingShifts_C",(Mat),(A));
    1892:   return(0);
    1893: }

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

    1898:    Logically Collective on Mat

    1900:     Input Parameters:
    1901: +   mat - the shell matrix
    1902: .   f - the function
    1903: .   base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated
    1904: -   ctx - an optional context for the function

    1906:    Output Parameter:
    1907: .   flg - PETSC_TRUE if the multiply is likely correct

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

    1912:    Level: advanced

    1914:    Fortran Notes:
    1915:     Not supported from Fortran

    1917: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellTestMultTranspose()
    1918: @*/
    1919: PetscErrorCode  MatShellTestMult(Mat mat,PetscErrorCode (*f)(void*,Vec,Vec),Vec base,void *ctx,PetscBool *flg)
    1920: {
    1922:   PetscInt       m,n;
    1923:   Mat            mf,Dmf,Dmat,Ddiff;
    1924:   PetscReal      Diffnorm,Dmfnorm;
    1925:   PetscBool      v = PETSC_FALSE, flag = PETSC_TRUE;

    1929:   PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_view",&v);
    1930:   MatGetLocalSize(mat,&m,&n);
    1931:   MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,PETSC_DECIDE,PETSC_DECIDE,&mf);
    1932:   MatMFFDSetFunction(mf,f,ctx);
    1933:   MatMFFDSetBase(mf,base,NULL);

    1935:   MatComputeOperator(mf,MATAIJ,&Dmf);
    1936:   MatComputeOperator(mat,MATAIJ,&Dmat);

    1938:   MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff);
    1939:   MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN);
    1940:   MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm);
    1941:   MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm);
    1942:   if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) {
    1943:     flag = PETSC_FALSE;
    1944:     if (v) {
    1945:       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));
    1946:       MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_view");
    1947:       MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_view");
    1948:       MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_view");
    1949:     }
    1950:   } else if (v) {
    1951:     PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL and matrix free multiple appear to produce the same results\n");
    1952:   }
    1953:   if (flg) *flg = flag;
    1954:   MatDestroy(&Ddiff);
    1955:   MatDestroy(&mf);
    1956:   MatDestroy(&Dmf);
    1957:   MatDestroy(&Dmat);
    1958:   return(0);
    1959: }

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

    1964:    Logically Collective on Mat

    1966:     Input Parameters:
    1967: +   mat - the shell matrix
    1968: .   f - the function
    1969: .   base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated
    1970: -   ctx - an optional context for the function

    1972:    Output Parameter:
    1973: .   flg - PETSC_TRUE if the multiply is likely correct

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

    1978:    Level: advanced

    1980:    Fortran Notes:
    1981:     Not supported from Fortran

    1983: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellTestMult()
    1984: @*/
    1985: PetscErrorCode  MatShellTestMultTranspose(Mat mat,PetscErrorCode (*f)(void*,Vec,Vec),Vec base,void *ctx,PetscBool *flg)
    1986: {
    1988:   Vec            x,y,z;
    1989:   PetscInt       m,n,M,N;
    1990:   Mat            mf,Dmf,Dmat,Ddiff;
    1991:   PetscReal      Diffnorm,Dmfnorm;
    1992:   PetscBool      v = PETSC_FALSE, flag = PETSC_TRUE;

    1996:   PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_transpose_view",&v);
    1997:   MatCreateVecs(mat,&x,&y);
    1998:   VecDuplicate(y,&z);
    1999:   MatGetLocalSize(mat,&m,&n);
    2000:   MatGetSize(mat,&M,&N);
    2001:   MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,M,N,&mf);
    2002:   MatMFFDSetFunction(mf,f,ctx);
    2003:   MatMFFDSetBase(mf,base,NULL);
    2004:   MatComputeOperator(mf,MATAIJ,&Dmf);
    2005:   MatTranspose(Dmf,MAT_INPLACE_MATRIX,&Dmf);
    2006:   MatComputeOperatorTranspose(mat,MATAIJ,&Dmat);

    2008:   MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff);
    2009:   MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN);
    2010:   MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm);
    2011:   MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm);
    2012:   if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) {
    2013:     flag = PETSC_FALSE;
    2014:     if (v) {
    2015:       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));
    2016:       MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");
    2017:       MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");
    2018:       MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");
    2019:     }
    2020:   } else if (v) {
    2021:     PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL transpose and matrix free multiple appear to produce the same results\n");
    2022:   }
    2023:   if (flg) *flg = flag;
    2024:   MatDestroy(&mf);
    2025:   MatDestroy(&Dmat);
    2026:   MatDestroy(&Ddiff);
    2027:   MatDestroy(&Dmf);
    2028:   VecDestroy(&x);
    2029:   VecDestroy(&y);
    2030:   VecDestroy(&z);
    2031:   return(0);
    2032: }

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

    2037:    Logically Collective on Mat

    2039:     Input Parameters:
    2040: +   mat - the shell matrix
    2041: .   op - the name of the operation
    2042: -   g - the function that provides the operation.

    2044:    Level: advanced

    2046:     Usage:
    2047: $      extern PetscErrorCode usermult(Mat,Vec,Vec);
    2048: $      MatCreateShell(comm,m,n,M,N,ctx,&A);
    2049: $      MatShellSetOperation(A,MATOP_MULT,(void(*)(void))usermult);

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

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

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

    2066:     Within each user-defined routine, the user should call
    2067:     MatShellGetContext() to obtain the user-defined context that was
    2068:     set by MatCreateShell().

    2070:     Use MatSetOperation() to set an operation for any matrix type. For matrix product operations (i.e. MatMat, MatTransposeMat etc) use MatShellSetMatProductOperation()

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

    2076: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatSetOperation(), MatShellSetManageScalingShifts(), MatShellSetMatProductOperation()
    2077: @*/
    2078: PetscErrorCode MatShellSetOperation(Mat mat,MatOperation op,void (*g)(void))
    2079: {

    2084:   PetscTryMethod(mat,"MatShellSetOperation_C",(Mat,MatOperation,void (*)(void)),(mat,op,g));
    2085:   return(0);
    2086: }

    2088: /*@C
    2089:     MatShellGetOperation - Gets a matrix function for a shell matrix.

    2091:     Not Collective

    2093:     Input Parameters:
    2094: +   mat - the shell matrix
    2095: -   op - the name of the operation

    2097:     Output Parameter:
    2098: .   g - the function that provides the operation.

    2100:     Level: advanced

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

    2108:     All user-provided functions have the same calling
    2109:     sequence as the usual matrix interface routines, since they
    2110:     are intended to be accessed via the usual matrix interface
    2111:     routines, e.g.,
    2112: $       MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec)

    2114:     Within each user-defined routine, the user should call
    2115:     MatShellGetContext() to obtain the user-defined context that was
    2116:     set by MatCreateShell().

    2118: .seealso: MatCreateShell(), MatShellGetContext(), MatShellSetOperation(), MatShellSetContext()
    2119: @*/
    2120: PetscErrorCode MatShellGetOperation(Mat mat,MatOperation op,void(**g)(void))
    2121: {

    2126:   PetscUseMethod(mat,"MatShellGetOperation_C",(Mat,MatOperation,void (**)(void)),(mat,op,g));
    2127:   return(0);
    2128: }

    2130: /*@
    2131:     MatIsShell - Inquires if a matrix is derived from MATSHELL

    2133:     Input Parameter:
    2134: .   mat - the matrix

    2136:     Output Parameter:
    2137: .   flg - the boolean value

    2139:     Level: developer

    2141:     Notes: in the future, we should allow the object type name to be changed still using the MatShell data structure for other matrices (i.e. MATTRANSPOSEMAT, MATSCHURCOMPLEMENT etc)

    2143: .seealso: MatCreateShell()
    2144: @*/
    2145: PetscErrorCode MatIsShell(Mat mat, PetscBool *flg)
    2146: {
    2150:   *flg = (PetscBool)(mat->ops->destroy == MatDestroy_Shell);
    2151:   return(0);
    2152: }