Actual source code: shell.c

petsc-3.15.0 2021-04-05
Report Typos and Errors

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

  8: #include <petsc/private/matimpl.h>
  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: }