Actual source code: shell.c

petsc-3.13.4 2020-08-01
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: typedef struct {
 20:   struct _MatShellOps ops[1];

 22:   PetscScalar vscale,vshift;
 23:   Vec         dshift;
 24:   Vec         left,right;
 25:   Vec         left_work,right_work;
 26:   Vec         left_add_work,right_add_work;
 27:   /* support for MatAXPY */
 28:   Mat              axpy;
 29:   PetscScalar      axpy_vscale;
 30:   PetscObjectState axpy_state;
 31:   PetscBool        managescalingshifts; /* The user will manage the scaling and shifts for the MATSHELL, not the default */
 32:   /* support for ZeroRows/Columns operations */
 33:   IS          zrows;
 34:   IS          zcols;
 35:   Vec         zvals;
 36:   Vec         zvals_w;
 37:   VecScatter  zvals_sct_r;
 38:   VecScatter  zvals_sct_c;
 39:   void        *ctx;
 40: } Mat_Shell;


 43: /*
 44:      Store and scale values on zeroed rows
 45:      xx = [x_1, 0], 0 on zeroed columns
 46: */
 47: static PetscErrorCode MatShellPreZeroRight(Mat A,Vec x,Vec *xx)
 48: {
 49:   Mat_Shell      *shell = (Mat_Shell*)A->data;

 53:   *xx = x;
 54:   if (shell->zrows) {
 55:     VecSet(shell->zvals_w,0.0);
 56:     VecScatterBegin(shell->zvals_sct_c,x,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);
 57:     VecScatterEnd(shell->zvals_sct_c,x,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);
 58:     VecPointwiseMult(shell->zvals_w,shell->zvals_w,shell->zvals);
 59:   }
 60:   if (shell->zcols) {
 61:     if (!shell->right_work) {
 62:       MatCreateVecs(A,&shell->right_work,NULL);
 63:     }
 64:     VecCopy(x,shell->right_work);
 65:     VecISSet(shell->right_work,shell->zcols,0.0);
 66:     *xx  = shell->right_work;
 67:   }
 68:   return(0);
 69: }

 71: /* Insert properly diagonally scaled values stored in MatShellPreZeroRight */
 72: static PetscErrorCode MatShellPostZeroLeft(Mat A,Vec x)
 73: {
 74:   Mat_Shell      *shell = (Mat_Shell*)A->data;

 78:   if (shell->zrows) {
 79:     VecScatterBegin(shell->zvals_sct_r,shell->zvals_w,x,INSERT_VALUES,SCATTER_REVERSE);
 80:     VecScatterEnd(shell->zvals_sct_r,shell->zvals_w,x,INSERT_VALUES,SCATTER_REVERSE);
 81:   }
 82:   return(0);
 83: }

 85: /*
 86:      Store and scale values on zeroed rows
 87:      xx = [x_1, 0], 0 on zeroed rows
 88: */
 89: static PetscErrorCode MatShellPreZeroLeft(Mat A,Vec x,Vec *xx)
 90: {
 91:   Mat_Shell      *shell = (Mat_Shell*)A->data;

 95:   *xx = NULL;
 96:   if (!shell->zrows) {
 97:     *xx = x;
 98:   } else {
 99:     if (!shell->left_work) {
100:       MatCreateVecs(A,NULL,&shell->left_work);
101:     }
102:     VecCopy(x,shell->left_work);
103:     VecSet(shell->zvals_w,0.0);
104:     VecScatterBegin(shell->zvals_sct_r,shell->zvals_w,shell->left_work,INSERT_VALUES,SCATTER_REVERSE);
105:     VecScatterEnd(shell->zvals_sct_r,shell->zvals_w,shell->left_work,INSERT_VALUES,SCATTER_REVERSE);
106:     VecScatterBegin(shell->zvals_sct_r,x,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);
107:     VecScatterEnd(shell->zvals_sct_r,x,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);
108:     VecPointwiseMult(shell->zvals_w,shell->zvals_w,shell->zvals);
109:     *xx  = shell->left_work;
110:   }
111:   return(0);
112: }

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

121:   if (shell->zcols) {
122:     VecISSet(x,shell->zcols,0.0);
123:   }
124:   if (shell->zrows) {
125:     VecScatterBegin(shell->zvals_sct_c,shell->zvals_w,x,ADD_VALUES,SCATTER_REVERSE);
126:     VecScatterEnd(shell->zvals_sct_c,shell->zvals_w,x,ADD_VALUES,SCATTER_REVERSE);
127:   }
128:   return(0);
129: }

131: /*
132:       xx = diag(left)*x
133: */
134: static PetscErrorCode MatShellPreScaleLeft(Mat A,Vec x,Vec *xx)
135: {
136:   Mat_Shell      *shell = (Mat_Shell*)A->data;

140:   *xx = NULL;
141:   if (!shell->left) {
142:     *xx = x;
143:   } else {
144:     if (!shell->left_work) {VecDuplicate(shell->left,&shell->left_work);}
145:     VecPointwiseMult(shell->left_work,x,shell->left);
146:     *xx  = shell->left_work;
147:   }
148:   return(0);
149: }

151: /*
152:      xx = diag(right)*x
153: */
154: static PetscErrorCode MatShellPreScaleRight(Mat A,Vec x,Vec *xx)
155: {
156:   Mat_Shell      *shell = (Mat_Shell*)A->data;

160:   *xx = NULL;
161:   if (!shell->right) {
162:     *xx = x;
163:   } else {
164:     if (!shell->right_work) {VecDuplicate(shell->right,&shell->right_work);}
165:     VecPointwiseMult(shell->right_work,x,shell->right);
166:     *xx  = shell->right_work;
167:   }
168:   return(0);
169: }

171: /*
172:     x = diag(left)*x
173: */
174: static PetscErrorCode MatShellPostScaleLeft(Mat A,Vec x)
175: {
176:   Mat_Shell      *shell = (Mat_Shell*)A->data;

180:   if (shell->left) {VecPointwiseMult(x,x,shell->left);}
181:   return(0);
182: }

184: /*
185:     x = diag(right)*x
186: */
187: static PetscErrorCode MatShellPostScaleRight(Mat A,Vec x)
188: {
189:   Mat_Shell      *shell = (Mat_Shell*)A->data;

193:   if (shell->right) {VecPointwiseMult(x,x,shell->right);}
194:   return(0);
195: }

197: /*
198:          Y = vscale*Y + diag(dshift)*X + vshift*X

200:          On input Y already contains A*x
201: */
202: static PetscErrorCode MatShellShiftAndScale(Mat A,Vec X,Vec Y)
203: {
204:   Mat_Shell      *shell = (Mat_Shell*)A->data;

208:   if (shell->dshift) {          /* get arrays because there is no VecPointwiseMultAdd() */
209:     PetscInt          i,m;
210:     const PetscScalar *x,*d;
211:     PetscScalar       *y;
212:     VecGetLocalSize(X,&m);
213:     VecGetArrayRead(shell->dshift,&d);
214:     VecGetArrayRead(X,&x);
215:     VecGetArray(Y,&y);
216:     for (i=0; i<m; i++) y[i] = shell->vscale*y[i] + d[i]*x[i];
217:     VecRestoreArrayRead(shell->dshift,&d);
218:     VecRestoreArrayRead(X,&x);
219:     VecRestoreArray(Y,&y);
220:   } else {
221:     VecScale(Y,shell->vscale);
222:   }
223:   if (shell->vshift != 0.0) {VecAXPY(Y,shell->vshift,X);} /* if test is for non-square matrices */
224:   return(0);
225: }

227: PetscErrorCode MatShellGetContext_Shell(Mat mat,void *ctx)
228: {
230:   *(void**)ctx = ((Mat_Shell*)(mat->data))->ctx;
231:   return(0);
232: }

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

237:     Not Collective

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

242:     Output Parameter:
243: .   ctx - the user provided context

245:     Level: advanced

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

251: .seealso: MatCreateShell(), MatShellSetOperation(), MatShellSetContext()
252: @*/
253: PetscErrorCode  MatShellGetContext(Mat mat,void *ctx)
254: {

260:   PetscUseMethod(mat,"MatShellGetContext_C",(Mat,void*),(mat,ctx));
261:   return(0);
262: }

264: static PetscErrorCode MatZeroRowsColumns_Local_Shell(Mat mat,PetscInt nr,PetscInt rows[],PetscInt nc,PetscInt cols[],PetscScalar diag,PetscBool rc)
265: {
267:   Mat_Shell      *shell = (Mat_Shell*)mat->data;
268:   Vec            x = NULL,b = NULL;
269:   IS             is1, is2;
270:   const PetscInt *ridxs;
271:   PetscInt       *idxs,*gidxs;
272:   PetscInt       cum,rst,cst,i;

275:   if (!shell->zvals) {
276:     MatCreateVecs(mat,NULL,&shell->zvals);
277:   }
278:   if (!shell->zvals_w) {
279:     VecDuplicate(shell->zvals,&shell->zvals_w);
280:   }
281:   MatGetOwnershipRange(mat,&rst,NULL);
282:   MatGetOwnershipRangeColumn(mat,&cst,NULL);

284:   /* Expand/create index set of zeroed rows */
285:   PetscMalloc1(nr,&idxs);
286:   for (i = 0; i < nr; i++) idxs[i] = rows[i] + rst;
287:   ISCreateGeneral(PETSC_COMM_SELF,nr,idxs,PETSC_OWN_POINTER,&is1);
288:   ISSort(is1);
289:   VecISSet(shell->zvals,is1,diag);
290:   if (shell->zrows) {
291:     ISSum(shell->zrows,is1,&is2);
292:     ISDestroy(&shell->zrows);
293:     ISDestroy(&is1);
294:     shell->zrows = is2;
295:   } else shell->zrows = is1;

297:   /* Create scatters for diagonal values communications */
298:   VecScatterDestroy(&shell->zvals_sct_c);
299:   VecScatterDestroy(&shell->zvals_sct_r);

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

305:   /* col scatter: from right vector to left vector */
306:   ISGetIndices(shell->zrows,&ridxs);
307:   ISGetLocalSize(shell->zrows,&nr);
308:   PetscMalloc1(nr,&gidxs);
309:   for (i = 0, cum  = 0; i < nr; i++) {
310:     if (ridxs[i] >= mat->cmap->N) continue;
311:     gidxs[cum] = ridxs[i];
312:     cum++;
313:   }
314:   ISCreateGeneral(PetscObjectComm((PetscObject)mat),cum,gidxs,PETSC_OWN_POINTER,&is1);
315:   VecScatterCreate(x,is1,shell->zvals_w,is1,&shell->zvals_sct_c);
316:   ISDestroy(&is1);
317:   VecDestroy(&x);
318:   VecDestroy(&b);

320:   /* Expand/create index set of zeroed columns */
321:   if (rc) {
322:     PetscMalloc1(nc,&idxs);
323:     for (i = 0; i < nc; i++) idxs[i] = cols[i] + cst;
324:     ISCreateGeneral(PETSC_COMM_SELF,nc,idxs,PETSC_OWN_POINTER,&is1);
325:     ISSort(is1);
326:     if (shell->zcols) {
327:       ISSum(shell->zcols,is1,&is2);
328:       ISDestroy(&shell->zcols);
329:       ISDestroy(&is1);
330:       shell->zcols = is2;
331:     } else shell->zcols = is1;
332:   }
333:   return(0);
334: }

336: static PetscErrorCode MatZeroRows_Shell(Mat mat,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
337: {
338:   Mat_Shell      *shell = (Mat_Shell*)mat->data;
339:   PetscInt       nr, *lrows;

343:   if (x && b) {
344:     Vec          xt;
345:     PetscScalar *vals;
346:     PetscInt    *gcols,i,st,nl,nc;

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

351:     MatCreateVecs(mat,&xt,NULL);
352:     VecCopy(x,xt);
353:     PetscCalloc1(nc,&vals);
354:     VecSetValues(xt,nc,gcols,vals,INSERT_VALUES); /* xt = [x1, 0] */
355:     PetscFree(vals);
356:     VecAssemblyBegin(xt);
357:     VecAssemblyEnd(xt);
358:     VecAYPX(xt,-1.0,x);                           /* xt = [0, x2] */

360:     VecGetOwnershipRange(xt,&st,NULL);
361:     VecGetLocalSize(xt,&nl);
362:     VecGetArray(xt,&vals);
363:     for (i = 0; i < nl; i++) {
364:       PetscInt g = i + st;
365:       if (g > mat->rmap->N) continue;
366:       if (PetscAbsScalar(vals[i]) == 0.0) continue;
367:       VecSetValue(b,g,diag*vals[i],INSERT_VALUES);
368:     }
369:     VecRestoreArray(xt,&vals);
370:     VecAssemblyBegin(b);
371:     VecAssemblyEnd(b);                            /* b  = [b1, x2 * diag] */
372:     VecDestroy(&xt);
373:     PetscFree(gcols);
374:   }
375:   PetscLayoutMapLocal(mat->rmap,n,rows,&nr,&lrows,NULL);
376:   MatZeroRowsColumns_Local_Shell(mat,nr,lrows,0,NULL,diag,PETSC_FALSE);
377:   if (shell->axpy) {
378:     MatZeroRows(shell->axpy,n,rows,0.0,NULL,NULL);
379:   }
380:   PetscFree(lrows);
381:   return(0);
382: }

384: static PetscErrorCode MatZeroRowsColumns_Shell(Mat mat,PetscInt n,const PetscInt rowscols[],PetscScalar diag,Vec x,Vec b)
385: {
386:   Mat_Shell      *shell = (Mat_Shell*)mat->data;
387:   PetscInt       *lrows, *lcols;
388:   PetscInt       nr, nc;
389:   PetscBool      congruent;

393:   if (x && b) {
394:     Vec          xt, bt;
395:     PetscScalar *vals;
396:     PetscInt    *grows,*gcols,i,st,nl;

398:     PetscMalloc2(n,&grows,n,&gcols);
399:     for (i = 0, nr = 0; i < n; i++) if (rowscols[i] < mat->rmap->N) grows[nr++] = rowscols[i];
400:     for (i = 0, nc = 0; i < n; i++) if (rowscols[i] < mat->cmap->N) gcols[nc++] = rowscols[i];
401:     PetscCalloc1(n,&vals);

403:     MatCreateVecs(mat,&xt,&bt);
404:     VecCopy(x,xt);
405:     VecSetValues(xt,nc,gcols,vals,INSERT_VALUES); /* xt = [x1, 0] */
406:     VecAssemblyBegin(xt);
407:     VecAssemblyEnd(xt);
408:     VecAXPY(xt,-1.0,x);                           /* xt = [0, -x2] */
409:     MatMult(mat,xt,bt);                           /* bt = [-A12*x2,-A22*x2] */
410:     VecSetValues(bt,nr,grows,vals,INSERT_VALUES); /* bt = [-A12*x2,0] */
411:     VecAssemblyBegin(bt);
412:     VecAssemblyEnd(bt);
413:     VecAXPY(b,1.0,bt);                            /* b  = [b1 - A12*x2, b2] */
414:     VecSetValues(bt,nr,grows,vals,INSERT_VALUES); /* b  = [b1 - A12*x2, 0] */
415:     VecAssemblyBegin(bt);
416:     VecAssemblyEnd(bt);
417:     PetscFree(vals);

419:     VecGetOwnershipRange(xt,&st,NULL);
420:     VecGetLocalSize(xt,&nl);
421:     VecGetArray(xt,&vals);
422:     for (i = 0; i < nl; i++) {
423:       PetscInt g = i + st;
424:       if (g > mat->rmap->N) continue;
425:       if (PetscAbsScalar(vals[i]) == 0.0) continue;
426:       VecSetValue(b,g,-diag*vals[i],INSERT_VALUES);
427:     }
428:     VecRestoreArray(xt,&vals);
429:     VecAssemblyBegin(b);
430:     VecAssemblyEnd(b);                            /* b  = [b1 - A12*x2, x2 * diag] */
431:     VecDestroy(&xt);
432:     VecDestroy(&bt);
433:     PetscFree2(grows,gcols);
434:   }
435:   PetscLayoutMapLocal(mat->rmap,n,rowscols,&nr,&lrows,NULL);
436:   MatHasCongruentLayouts(mat,&congruent);
437:   if (congruent) {
438:     nc    = nr;
439:     lcols = lrows;
440:   } else { /* MatZeroRowsColumns implicitly assumes the rowscols indices are for a square matrix, here we handle a more general case */
441:     PetscInt i,nt,*t;

443:     PetscMalloc1(n,&t);
444:     for (i = 0, nt = 0; i < n; i++) if (rowscols[i] < mat->cmap->N) t[nt++] = rowscols[i];
445:     PetscLayoutMapLocal(mat->cmap,nt,t,&nc,&lcols,NULL);
446:     PetscFree(t);
447:   }
448:   MatZeroRowsColumns_Local_Shell(mat,nr,lrows,nc,lcols,diag,PETSC_TRUE);
449:   if (!congruent) {
450:     PetscFree(lcols);
451:   }
452:   PetscFree(lrows);
453:   if (shell->axpy) {
454:     MatZeroRowsColumns(shell->axpy,n,rowscols,0.0,NULL,NULL);
455:   }
456:   return(0);
457: }

459: PetscErrorCode MatDestroy_Shell(Mat mat)
460: {
462:   Mat_Shell      *shell = (Mat_Shell*)mat->data;

465:   if (shell->ops->destroy) {
466:     (*shell->ops->destroy)(mat);
467:   }
468:   PetscMemzero(shell->ops,sizeof(struct _MatShellOps));
469:   VecDestroy(&shell->left);
470:   VecDestroy(&shell->right);
471:   VecDestroy(&shell->dshift);
472:   VecDestroy(&shell->left_work);
473:   VecDestroy(&shell->right_work);
474:   VecDestroy(&shell->left_add_work);
475:   VecDestroy(&shell->right_add_work);
476:   MatDestroy(&shell->axpy);
477:   VecDestroy(&shell->zvals_w);
478:   VecDestroy(&shell->zvals);
479:   VecScatterDestroy(&shell->zvals_sct_c);
480:   VecScatterDestroy(&shell->zvals_sct_r);
481:   ISDestroy(&shell->zrows);
482:   ISDestroy(&shell->zcols);
483:   PetscFree(mat->data);
484:   PetscObjectComposeFunction((PetscObject)mat,"MatShellGetContext_C",NULL);
485:   PetscObjectComposeFunction((PetscObject)mat,"MatShellSetContext_C",NULL);
486:   PetscObjectComposeFunction((PetscObject)mat,"MatShellSetVecType_C",NULL);
487:   PetscObjectComposeFunction((PetscObject)mat,"MatShellSetManageScalingShifts_C",NULL);
488:   PetscObjectComposeFunction((PetscObject)mat,"MatShellSetOperation_C",NULL);
489:   PetscObjectComposeFunction((PetscObject)mat,"MatShellGetOperation_C",NULL);
490:   return(0);
491: }

493: PetscErrorCode MatCopy_Shell(Mat A,Mat B,MatStructure str)
494: {
495:   Mat_Shell       *shellA = (Mat_Shell*)A->data,*shellB = (Mat_Shell*)B->data;
496:   PetscErrorCode  ierr;
497:   PetscBool       matflg;

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

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

506:   if (shellA->ops->copy) {
507:     (*shellA->ops->copy)(A,B,str);
508:   }
509:   shellB->vscale = shellA->vscale;
510:   shellB->vshift = shellA->vshift;
511:   if (shellA->dshift) {
512:     if (!shellB->dshift) {
513:       VecDuplicate(shellA->dshift,&shellB->dshift);
514:     }
515:     VecCopy(shellA->dshift,shellB->dshift);
516:   } else {
517:     VecDestroy(&shellB->dshift);
518:   }
519:   if (shellA->left) {
520:     if (!shellB->left) {
521:       VecDuplicate(shellA->left,&shellB->left);
522:     }
523:     VecCopy(shellA->left,shellB->left);
524:   } else {
525:     VecDestroy(&shellB->left);
526:   }
527:   if (shellA->right) {
528:     if (!shellB->right) {
529:       VecDuplicate(shellA->right,&shellB->right);
530:     }
531:     VecCopy(shellA->right,shellB->right);
532:   } else {
533:     VecDestroy(&shellB->right);
534:   }
535:   MatDestroy(&shellB->axpy);
536:   shellB->axpy_vscale = 0.0;
537:   shellB->axpy_state  = 0;
538:   if (shellA->axpy) {
539:     PetscObjectReference((PetscObject)shellA->axpy);
540:     shellB->axpy        = shellA->axpy;
541:     shellB->axpy_vscale = shellA->axpy_vscale;
542:     shellB->axpy_state  = shellA->axpy_state;
543:   }
544:   if (shellA->zrows) {
545:     ISDuplicate(shellA->zrows,&shellB->zrows);
546:     if (shellA->zcols) {
547:       ISDuplicate(shellA->zcols,&shellB->zcols);
548:     }
549:     VecDuplicate(shellA->zvals,&shellB->zvals);
550:     VecCopy(shellA->zvals,shellB->zvals);
551:     VecDuplicate(shellA->zvals_w,&shellB->zvals_w);
552:     PetscObjectReference((PetscObject)shellA->zvals_sct_r);
553:     PetscObjectReference((PetscObject)shellA->zvals_sct_c);
554:     shellB->zvals_sct_r = shellA->zvals_sct_r;
555:     shellB->zvals_sct_c = shellA->zvals_sct_c;
556:   }
557:   return(0);
558: }

560: PetscErrorCode MatDuplicate_Shell(Mat mat,MatDuplicateOption op,Mat *M)
561: {
563:   void           *ctx;

566:   MatShellGetContext(mat,&ctx);
567:   MatCreateShell(PetscObjectComm((PetscObject)mat),mat->rmap->n,mat->cmap->n,mat->rmap->N,mat->cmap->N,ctx,M);
568:   if (op != MAT_DO_NOT_COPY_VALUES) {
569:     MatCopy(mat,*M,SAME_NONZERO_PATTERN);
570:   }
571:   return(0);
572: }

574: PetscErrorCode MatMult_Shell(Mat A,Vec x,Vec y)
575: {
576:   Mat_Shell        *shell = (Mat_Shell*)A->data;
577:   PetscErrorCode   ierr;
578:   Vec              xx;
579:   PetscObjectState instate,outstate;

582:   if (!shell->ops->mult) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Have not provided a MatMult() for this MATSHELL");
583:   MatShellPreZeroRight(A,x,&xx);
584:   MatShellPreScaleRight(A,xx,&xx);
585:   PetscObjectStateGet((PetscObject)y, &instate);
586:   (*shell->ops->mult)(A,xx,y);
587:   PetscObjectStateGet((PetscObject)y, &outstate);
588:   if (instate == outstate) {
589:     /* increase the state of the output vector since the user did not update its state themself as should have been done */
590:     PetscObjectStateIncrease((PetscObject)y);
591:   }
592:   MatShellShiftAndScale(A,xx,y);
593:   MatShellPostScaleLeft(A,y);
594:   MatShellPostZeroLeft(A,y);

596:   if (shell->axpy) {
597:     Mat              X;
598:     PetscObjectState axpy_state;

600:     MatShellGetContext(shell->axpy,(void *)&X);
601:     PetscObjectStateGet((PetscObject)X,&axpy_state);
602:     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,...)");
603:     if (!shell->left_work) {MatCreateVecs(A,&shell->left_work,NULL);}
604:     MatMult(shell->axpy,x,shell->left_work);
605:     VecAXPY(y,shell->axpy_vscale,shell->left_work);
606:   }
607:   return(0);
608: }

610: PetscErrorCode MatMultAdd_Shell(Mat A,Vec x,Vec y,Vec z)
611: {
612:   Mat_Shell      *shell = (Mat_Shell*)A->data;

616:   if (y == z) {
617:     if (!shell->right_add_work) {VecDuplicate(z,&shell->right_add_work);}
618:     MatMult(A,x,shell->right_add_work);
619:     VecAXPY(z,1.0,shell->right_add_work);
620:   } else {
621:     MatMult(A,x,z);
622:     VecAXPY(z,1.0,y);
623:   }
624:   return(0);
625: }

627: PetscErrorCode MatMultTranspose_Shell(Mat A,Vec x,Vec y)
628: {
629:   Mat_Shell        *shell = (Mat_Shell*)A->data;
630:   PetscErrorCode   ierr;
631:   Vec              xx;
632:   PetscObjectState instate,outstate;

635:   MatShellPreZeroLeft(A,x,&xx);
636:   MatShellPreScaleLeft(A,xx,&xx);
637:   PetscObjectStateGet((PetscObject)y, &instate);
638:   if (!shell->ops->multtranspose) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Have not provided a MatMultTranspose() for this MATSHELL");
639:   (*shell->ops->multtranspose)(A,xx,y);
640:   PetscObjectStateGet((PetscObject)y, &outstate);
641:   if (instate == outstate) {
642:     /* increase the state of the output vector since the user did not update its state themself as should have been done */
643:     PetscObjectStateIncrease((PetscObject)y);
644:   }
645:   MatShellShiftAndScale(A,xx,y);
646:   MatShellPostScaleRight(A,y);
647:   MatShellPostZeroRight(A,y);

649:   if (shell->axpy) {
650:     Mat              X;
651:     PetscObjectState axpy_state;

653:     MatShellGetContext(shell->axpy,(void *)&X);
654:     PetscObjectStateGet((PetscObject)X,&axpy_state);
655:     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,...)");
656:     if (!shell->right_work) {MatCreateVecs(A,NULL,&shell->right_work);}
657:     MatMultTranspose(shell->axpy,x,shell->right_work);
658:     VecAXPY(y,shell->axpy_vscale,shell->right_work);
659:   }
660:   return(0);
661: }

663: PetscErrorCode MatMultTransposeAdd_Shell(Mat A,Vec x,Vec y,Vec z)
664: {
665:   Mat_Shell      *shell = (Mat_Shell*)A->data;

669:   if (y == z) {
670:     if (!shell->left_add_work) {VecDuplicate(z,&shell->left_add_work);}
671:     MatMultTranspose(A,x,shell->left_add_work);
672:     VecAXPY(z,1.0,shell->left_add_work);
673:   } else {
674:     MatMultTranspose(A,x,z);
675:     VecAXPY(z,1.0,y);
676:   }
677:   return(0);
678: }

680: /*
681:           diag(left)(vscale*A + diag(dshift) + vshift I)diag(right)
682: */
683: PetscErrorCode MatGetDiagonal_Shell(Mat A,Vec v)
684: {
685:   Mat_Shell      *shell = (Mat_Shell*)A->data;

689:   if (shell->ops->getdiagonal) {
690:     (*shell->ops->getdiagonal)(A,v);
691:   } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Must provide shell matrix with routine to return diagonal using\nMatShellSetOperation(S,MATOP_GET_DIAGONAL,...)");
692:   VecScale(v,shell->vscale);
693:   if (shell->dshift) {
694:     VecAXPY(v,1.0,shell->dshift);
695:   }
696:   VecShift(v,shell->vshift);
697:   if (shell->left)  {VecPointwiseMult(v,v,shell->left);}
698:   if (shell->right) {VecPointwiseMult(v,v,shell->right);}
699:   if (shell->zrows) {
700:     VecScatterBegin(shell->zvals_sct_r,shell->zvals,v,INSERT_VALUES,SCATTER_REVERSE);
701:     VecScatterEnd(shell->zvals_sct_r,shell->zvals,v,INSERT_VALUES,SCATTER_REVERSE);
702:   }
703:   if (shell->axpy) {
704:     Mat              X;
705:     PetscObjectState axpy_state;

707:     MatShellGetContext(shell->axpy,(void *)&X);
708:     PetscObjectStateGet((PetscObject)X,&axpy_state);
709:     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,...)");
710:     if (!shell->left_work) {VecDuplicate(v,&shell->left_work);}
711:     MatGetDiagonal(shell->axpy,shell->left_work);
712:     VecAXPY(v,shell->axpy_vscale,shell->left_work);
713:   }
714:   return(0);
715: }

717: PetscErrorCode MatShift_Shell(Mat Y,PetscScalar a)
718: {
719:   Mat_Shell      *shell = (Mat_Shell*)Y->data;
721:   PetscBool      flg;

724:   MatHasCongruentLayouts(Y,&flg);
725:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_SUP,"Cannot shift shell matrix if it is not congruent");
726:   if (shell->left || shell->right) {
727:     if (!shell->dshift) {
728:       VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift);
729:       VecSet(shell->dshift,a);
730:     } else {
731:       if (shell->left)  {VecPointwiseMult(shell->dshift,shell->dshift,shell->left);}
732:       if (shell->right) {VecPointwiseMult(shell->dshift,shell->dshift,shell->right);}
733:       VecShift(shell->dshift,a);
734:     }
735:     if (shell->left)  {VecPointwiseDivide(shell->dshift,shell->dshift,shell->left);}
736:     if (shell->right) {VecPointwiseDivide(shell->dshift,shell->dshift,shell->right);}
737:   } else shell->vshift += a;
738:   if (shell->zrows) {
739:     VecShift(shell->zvals,a);
740:   }
741:   return(0);
742: }

744: PetscErrorCode MatDiagonalSet_Shell_Private(Mat A,Vec D,PetscScalar s)
745: {
746:   Mat_Shell      *shell = (Mat_Shell*)A->data;

750:   if (!shell->dshift) {VecDuplicate(D,&shell->dshift);}
751:   if (shell->left || shell->right) {
752:     if (!shell->right_work) {VecDuplicate(shell->left ? shell->left : shell->right, &shell->right_work);}
753:     if (shell->left && shell->right)  {
754:       VecPointwiseDivide(shell->right_work,D,shell->left);
755:       VecPointwiseDivide(shell->right_work,shell->right_work,shell->right);
756:     } else if (shell->left) {
757:       VecPointwiseDivide(shell->right_work,D,shell->left);
758:     } else {
759:       VecPointwiseDivide(shell->right_work,D,shell->right);
760:     }
761:     VecAXPY(shell->dshift,s,shell->right_work);
762:   } else {
763:     VecAXPY(shell->dshift,s,D);
764:   }
765:   return(0);
766: }

768: PetscErrorCode MatDiagonalSet_Shell(Mat A,Vec D,InsertMode ins)
769: {
770:   Mat_Shell      *shell = (Mat_Shell*)A->data;
771:   Vec            d;
773:   PetscBool      flg;

776:   MatHasCongruentLayouts(A,&flg);
777:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Cannot diagonal set or shift shell matrix if it is not congruent");
778:   if (ins == INSERT_VALUES) {
779:     if (!A->ops->getdiagonal) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Operation MATOP_GETDIAGONAL must be set first");
780:     VecDuplicate(D,&d);
781:     MatGetDiagonal(A,d);
782:     MatDiagonalSet_Shell_Private(A,d,-1.);
783:     MatDiagonalSet_Shell_Private(A,D,1.);
784:     VecDestroy(&d);
785:     if (shell->zrows) {
786:       VecCopy(D,shell->zvals);
787:     }
788:   } else {
789:     MatDiagonalSet_Shell_Private(A,D,1.);
790:     if (shell->zrows) {
791:       VecAXPY(shell->zvals,1.0,D);
792:     }
793:   }
794:   return(0);
795: }

797: PetscErrorCode MatScale_Shell(Mat Y,PetscScalar a)
798: {
799:   Mat_Shell      *shell = (Mat_Shell*)Y->data;

803:   shell->vscale *= a;
804:   shell->vshift *= a;
805:   if (shell->dshift) {
806:     VecScale(shell->dshift,a);
807:   }
808:   shell->axpy_vscale *= a;
809:   if (shell->zrows) {
810:     VecScale(shell->zvals,a);
811:   }
812:   return(0);
813: }

815: static PetscErrorCode MatDiagonalScale_Shell(Mat Y,Vec left,Vec right)
816: {
817:   Mat_Shell      *shell = (Mat_Shell*)Y->data;

821:   if (left) {
822:     if (!shell->left) {
823:       VecDuplicate(left,&shell->left);
824:       VecCopy(left,shell->left);
825:     } else {
826:       VecPointwiseMult(shell->left,shell->left,left);
827:     }
828:     if (shell->zrows) {
829:       VecPointwiseMult(shell->zvals,shell->zvals,left);
830:     }
831:   }
832:   if (right) {
833:     if (!shell->right) {
834:       VecDuplicate(right,&shell->right);
835:       VecCopy(right,shell->right);
836:     } else {
837:       VecPointwiseMult(shell->right,shell->right,right);
838:     }
839:     if (shell->zrows) {
840:       if (!shell->left_work) {
841:         MatCreateVecs(Y,NULL,&shell->left_work);
842:       }
843:       VecSet(shell->zvals_w,1.0);
844:       VecScatterBegin(shell->zvals_sct_c,right,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);
845:       VecScatterEnd(shell->zvals_sct_c,right,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);
846:       VecPointwiseMult(shell->zvals,shell->zvals,shell->zvals_w);
847:     }
848:   }
849:   if (shell->axpy) {
850:     MatDiagonalScale(shell->axpy,left,right);
851:   }
852:   return(0);
853: }

855: PetscErrorCode MatAssemblyEnd_Shell(Mat Y,MatAssemblyType t)
856: {
857:   Mat_Shell      *shell = (Mat_Shell*)Y->data;

861:   if (t == MAT_FINAL_ASSEMBLY) {
862:     shell->vshift = 0.0;
863:     shell->vscale = 1.0;
864:     shell->axpy_vscale = 0.0;
865:     shell->axpy_state  = 0;
866:     VecDestroy(&shell->dshift);
867:     VecDestroy(&shell->left);
868:     VecDestroy(&shell->right);
869:     MatDestroy(&shell->axpy);
870:     VecScatterDestroy(&shell->zvals_sct_c);
871:     VecScatterDestroy(&shell->zvals_sct_r);
872:     ISDestroy(&shell->zrows);
873:     ISDestroy(&shell->zcols);
874:   }
875:   return(0);
876: }

878: static PetscErrorCode MatMissingDiagonal_Shell(Mat A,PetscBool *missing,PetscInt *d)
879: {
881:   *missing = PETSC_FALSE;
882:   return(0);
883: }

885: PetscErrorCode MatAXPY_Shell(Mat Y,PetscScalar a,Mat X,MatStructure str)
886: {
887:   Mat_Shell      *shell = (Mat_Shell*)Y->data;

891:   if (X == Y) {
892:     MatScale(Y,1.0 + a);
893:     return(0);
894:   }
895:   if (!shell->axpy) {
896:     MatConvertFrom_Shell(X,MATSHELL,MAT_INITIAL_MATRIX,&shell->axpy);
897:     shell->axpy_vscale = a;
898:     PetscObjectStateGet((PetscObject)X,&shell->axpy_state);
899:   } else {
900:     MatAXPY(shell->axpy,a/shell->axpy_vscale,X,str);
901:   }
902:   return(0);
903: }

905: static struct _MatOps MatOps_Values = {0,
906:                                        0,
907:                                        0,
908:                                        0,
909:                                 /* 4*/ MatMultAdd_Shell,
910:                                        0,
911:                                        MatMultTransposeAdd_Shell,
912:                                        0,
913:                                        0,
914:                                        0,
915:                                 /*10*/ 0,
916:                                        0,
917:                                        0,
918:                                        0,
919:                                        0,
920:                                 /*15*/ 0,
921:                                        0,
922:                                        0,
923:                                        MatDiagonalScale_Shell,
924:                                        0,
925:                                 /*20*/ 0,
926:                                        MatAssemblyEnd_Shell,
927:                                        0,
928:                                        0,
929:                                 /*24*/ MatZeroRows_Shell,
930:                                        0,
931:                                        0,
932:                                        0,
933:                                        0,
934:                                 /*29*/ 0,
935:                                        0,
936:                                        0,
937:                                        0,
938:                                        0,
939:                                 /*34*/ MatDuplicate_Shell,
940:                                        0,
941:                                        0,
942:                                        0,
943:                                        0,
944:                                 /*39*/ MatAXPY_Shell,
945:                                        0,
946:                                        0,
947:                                        0,
948:                                        MatCopy_Shell,
949:                                 /*44*/ 0,
950:                                        MatScale_Shell,
951:                                        MatShift_Shell,
952:                                        MatDiagonalSet_Shell,
953:                                        MatZeroRowsColumns_Shell,
954:                                 /*49*/ 0,
955:                                        0,
956:                                        0,
957:                                        0,
958:                                        0,
959:                                 /*54*/ 0,
960:                                        0,
961:                                        0,
962:                                        0,
963:                                        0,
964:                                 /*59*/ 0,
965:                                        MatDestroy_Shell,
966:                                        0,
967:                                        MatConvertFrom_Shell,
968:                                        0,
969:                                 /*64*/ 0,
970:                                        0,
971:                                        0,
972:                                        0,
973:                                        0,
974:                                 /*69*/ 0,
975:                                        0,
976:                                        MatConvert_Shell,
977:                                        0,
978:                                        0,
979:                                 /*74*/ 0,
980:                                        0,
981:                                        0,
982:                                        0,
983:                                        0,
984:                                 /*79*/ 0,
985:                                        0,
986:                                        0,
987:                                        0,
988:                                        0,
989:                                 /*84*/ 0,
990:                                        0,
991:                                        0,
992:                                        0,
993:                                        0,
994:                                 /*89*/ 0,
995:                                        0,
996:                                        0,
997:                                        0,
998:                                        0,
999:                                 /*94*/ 0,
1000:                                        0,
1001:                                        0,
1002:                                        0,
1003:                                        0,
1004:                                 /*99*/ 0,
1005:                                        0,
1006:                                        0,
1007:                                        0,
1008:                                        0,
1009:                                /*104*/ 0,
1010:                                        0,
1011:                                        0,
1012:                                        0,
1013:                                        0,
1014:                                /*109*/ 0,
1015:                                        0,
1016:                                        0,
1017:                                        0,
1018:                                        MatMissingDiagonal_Shell,
1019:                                /*114*/ 0,
1020:                                        0,
1021:                                        0,
1022:                                        0,
1023:                                        0,
1024:                                /*119*/ 0,
1025:                                        0,
1026:                                        0,
1027:                                        0,
1028:                                        0,
1029:                                /*124*/ 0,
1030:                                        0,
1031:                                        0,
1032:                                        0,
1033:                                        0,
1034:                                /*129*/ 0,
1035:                                        0,
1036:                                        0,
1037:                                        0,
1038:                                        0,
1039:                                /*134*/ 0,
1040:                                        0,
1041:                                        0,
1042:                                        0,
1043:                                        0,
1044:                                /*139*/ 0,
1045:                                        0,
1046:                                        0
1047: };

1049: PetscErrorCode  MatShellSetContext_Shell(Mat mat,void *ctx)
1050: {
1051:   Mat_Shell *shell = (Mat_Shell*)mat->data;

1054:   shell->ctx = ctx;
1055:   return(0);
1056: }

1058: static PetscErrorCode MatShellSetVecType_Shell(Mat mat,VecType vtype)
1059: {

1063:   PetscFree(mat->defaultvectype);
1064:   PetscStrallocpy(vtype,(char**)&mat->defaultvectype);
1065:   return(0);
1066: }

1068: PetscErrorCode MatShellSetManageScalingShifts_Shell(Mat A)
1069: {
1070:   Mat_Shell      *shell = (Mat_Shell*)A->data;

1073:   shell->managescalingshifts = PETSC_FALSE;
1074:   A->ops->diagonalset   = NULL;
1075:   A->ops->diagonalscale = NULL;
1076:   A->ops->scale         = NULL;
1077:   A->ops->shift         = NULL;
1078:   A->ops->axpy          = NULL;
1079:   return(0);
1080: }

1082: PetscErrorCode MatShellSetOperation_Shell(Mat mat,MatOperation op,void (*f)(void))
1083: {
1084:   Mat_Shell      *shell = (Mat_Shell*)mat->data;

1087:   switch (op) {
1088:   case MATOP_DESTROY:
1089:     shell->ops->destroy = (PetscErrorCode (*)(Mat))f;
1090:     break;
1091:   case MATOP_VIEW:
1092:     if (!mat->ops->viewnative) {
1093:       mat->ops->viewnative = mat->ops->view;
1094:     }
1095:     mat->ops->view = (PetscErrorCode (*)(Mat,PetscViewer))f;
1096:     break;
1097:   case MATOP_COPY:
1098:     shell->ops->copy = (PetscErrorCode (*)(Mat,Mat,MatStructure))f;
1099:     break;
1100:   case MATOP_DIAGONAL_SET:
1101:   case MATOP_DIAGONAL_SCALE:
1102:   case MATOP_SHIFT:
1103:   case MATOP_SCALE:
1104:   case MATOP_AXPY:
1105:   case MATOP_ZERO_ROWS:
1106:   case MATOP_ZERO_ROWS_COLUMNS:
1107:     if (shell->managescalingshifts) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"MATSHELL is managing scalings and shifts, see MatShellSetManageScalingShifts()");
1108:     (((void(**)(void))mat->ops)[op]) = f;
1109:     break;
1110:   case MATOP_GET_DIAGONAL:
1111:     if (shell->managescalingshifts) {
1112:       shell->ops->getdiagonal = (PetscErrorCode (*)(Mat,Vec))f;
1113:       mat->ops->getdiagonal   = MatGetDiagonal_Shell;
1114:     } else {
1115:       shell->ops->getdiagonal = NULL;
1116:       mat->ops->getdiagonal   = (PetscErrorCode (*)(Mat,Vec))f;
1117:     }
1118:     break;
1119:   case MATOP_MULT:
1120:     if (shell->managescalingshifts) {
1121:       shell->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1122:       mat->ops->mult   = MatMult_Shell;
1123:     } else {
1124:       shell->ops->mult = NULL;
1125:       mat->ops->mult   = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1126:     }
1127:     break;
1128:   case MATOP_MULT_TRANSPOSE:
1129:     if (shell->managescalingshifts) {
1130:       shell->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1131:       mat->ops->multtranspose   = MatMultTranspose_Shell;
1132:     } else {
1133:       shell->ops->multtranspose = NULL;
1134:       mat->ops->multtranspose   = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1135:     }
1136:     break;
1137:   default:
1138:     (((void(**)(void))mat->ops)[op]) = f;
1139:     break;
1140:   }
1141:   return(0);
1142: }

1144: PetscErrorCode MatShellGetOperation_Shell(Mat mat,MatOperation op,void(**f)(void))
1145: {
1146:   Mat_Shell      *shell = (Mat_Shell*)mat->data;

1149:   switch (op) {
1150:   case MATOP_DESTROY:
1151:     *f = (void (*)(void))shell->ops->destroy;
1152:     break;
1153:   case MATOP_VIEW:
1154:     *f = (void (*)(void))mat->ops->view;
1155:     break;
1156:   case MATOP_COPY:
1157:     *f = (void (*)(void))shell->ops->copy;
1158:     break;
1159:   case MATOP_DIAGONAL_SET:
1160:   case MATOP_DIAGONAL_SCALE:
1161:   case MATOP_SHIFT:
1162:   case MATOP_SCALE:
1163:   case MATOP_AXPY:
1164:   case MATOP_ZERO_ROWS:
1165:   case MATOP_ZERO_ROWS_COLUMNS:
1166:     *f = (((void (**)(void))mat->ops)[op]);
1167:     break;
1168:   case MATOP_GET_DIAGONAL:
1169:     if (shell->ops->getdiagonal)
1170:       *f = (void (*)(void))shell->ops->getdiagonal;
1171:     else
1172:       *f = (((void (**)(void))mat->ops)[op]);
1173:     break;
1174:   case MATOP_MULT:
1175:     if (shell->ops->mult)
1176:       *f = (void (*)(void))shell->ops->mult;
1177:     else
1178:       *f = (((void (**)(void))mat->ops)[op]);
1179:     break;
1180:   case MATOP_MULT_TRANSPOSE:
1181:     if (shell->ops->multtranspose)
1182:       *f = (void (*)(void))shell->ops->multtranspose;
1183:     else
1184:       *f = (((void (**)(void))mat->ops)[op]);
1185:     break;
1186:   default:
1187:     *f = (((void (**)(void))mat->ops)[op]);
1188:   }
1189:   return(0);
1190: }

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

1195:   Level: advanced

1197: .seealso: MatCreateShell()
1198: M*/

1200: PETSC_EXTERN PetscErrorCode MatCreate_Shell(Mat A)
1201: {
1202:   Mat_Shell      *b;

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

1208:   PetscNewLog(A,&b);
1209:   A->data = (void*)b;

1211:   PetscLayoutSetUp(A->rmap);
1212:   PetscLayoutSetUp(A->cmap);

1214:   b->ctx                 = 0;
1215:   b->vshift              = 0.0;
1216:   b->vscale              = 1.0;
1217:   b->managescalingshifts = PETSC_TRUE;
1218:   A->assembled           = PETSC_TRUE;
1219:   A->preallocated        = PETSC_FALSE;

1221:   PetscObjectComposeFunction((PetscObject)A,"MatShellGetContext_C",MatShellGetContext_Shell);
1222:   PetscObjectComposeFunction((PetscObject)A,"MatShellSetContext_C",MatShellSetContext_Shell);
1223:   PetscObjectComposeFunction((PetscObject)A,"MatShellSetVecType_C",MatShellSetVecType_Shell);
1224:   PetscObjectComposeFunction((PetscObject)A,"MatShellSetManageScalingShifts_C",MatShellSetManageScalingShifts_Shell);
1225:   PetscObjectComposeFunction((PetscObject)A,"MatShellSetOperation_C",MatShellSetOperation_Shell);
1226:   PetscObjectComposeFunction((PetscObject)A,"MatShellGetOperation_C",MatShellGetOperation_Shell);
1227:   PetscObjectChangeTypeName((PetscObject)A,MATSHELL);
1228:   return(0);
1229: }

1231: /*@C
1232:    MatCreateShell - Creates a new matrix class for use with a user-defined
1233:    private data storage format.

1235:   Collective

1237:    Input Parameters:
1238: +  comm - MPI communicator
1239: .  m - number of local rows (must be given)
1240: .  n - number of local columns (must be given)
1241: .  M - number of global rows (may be PETSC_DETERMINE)
1242: .  N - number of global columns (may be PETSC_DETERMINE)
1243: -  ctx - pointer to data needed by the shell matrix routines

1245:    Output Parameter:
1246: .  A - the matrix

1248:    Level: advanced

1250:   Usage:
1251: $    extern PetscErrorCode mult(Mat,Vec,Vec);
1252: $    MatCreateShell(comm,m,n,M,N,ctx,&mat);
1253: $    MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
1254: $    [ Use matrix for operations that have been set ]
1255: $    MatDestroy(mat);

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

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

1267:    PETSc requires that matrices and vectors being used for certain
1268:    operations are partitioned accordingly.  For example, when
1269:    creating a shell matrix, A, that supports parallel matrix-vector
1270:    products using MatMult(A,x,y) the user should set the number
1271:    of local matrix rows to be the number of local elements of the
1272:    corresponding result vector, y. Note that this is information is
1273:    required for use of the matrix interface routines, even though
1274:    the shell matrix may not actually be physically partitioned.
1275:    For example,

1277: $
1278: $     Vec x, y
1279: $     extern PetscErrorCode mult(Mat,Vec,Vec);
1280: $     Mat A
1281: $
1282: $     VecCreateMPI(comm,PETSC_DECIDE,M,&y);
1283: $     VecCreateMPI(comm,PETSC_DECIDE,N,&x);
1284: $     VecGetLocalSize(y,&m);
1285: $     VecGetLocalSize(x,&n);
1286: $     MatCreateShell(comm,m,n,M,N,ctx,&A);
1287: $     MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
1288: $     MatMult(A,x,y);
1289: $     MatDestroy(&A);
1290: $     VecDestroy(&y);
1291: $     VecDestroy(&x);
1292: $


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


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

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

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

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

1310:           A is the user provided function.

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

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

1320: .seealso: MatShellSetOperation(), MatHasOperation(), MatShellGetContext(), MatShellSetContext(), MATSHELL, MatShellSetManageScalingShifts()
1321: @*/
1322: PetscErrorCode  MatCreateShell(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,void *ctx,Mat *A)
1323: {

1327:   MatCreate(comm,A);
1328:   MatSetSizes(*A,m,n,M,N);
1329:   MatSetType(*A,MATSHELL);
1330:   MatShellSetContext(*A,ctx);
1331:   MatSetUp(*A);
1332:   return(0);
1333: }


1336: /*@
1337:     MatShellSetContext - sets the context for a shell matrix

1339:    Logically Collective on Mat

1341:     Input Parameters:
1342: +   mat - the shell matrix
1343: -   ctx - the context

1345:    Level: advanced

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

1351: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation()
1352: @*/
1353: PetscErrorCode  MatShellSetContext(Mat mat,void *ctx)
1354: {

1359:   PetscTryMethod(mat,"MatShellSetContext_C",(Mat,void*),(mat,ctx));
1360:   return(0);
1361: }

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

1366:  Logically collective

1368:     Input Parameters:
1369: +   mat   - the shell matrix
1370: -   vtype - type to use for creating vectors

1372:  Notes:

1374:  Level: advanced

1376: .seealso: MatCreateVecs()
1377: @*/
1378: PetscErrorCode  MatShellSetVecType(Mat mat,VecType vtype)
1379: {

1383:   PetscTryMethod(mat,"MatShellSetVecType_C",(Mat,VecType),(mat,vtype));
1384:   return(0);
1385: }

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

1391:    Logically Collective on Mat

1393:     Input Parameter:
1394: .   mat - the shell matrix

1396:   Level: advanced

1398: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatShellSetOperation()
1399: @*/
1400: PetscErrorCode MatShellSetManageScalingShifts(Mat A)
1401: {

1406:   PetscTryMethod(A,"MatShellSetManageScalingShifts_C",(Mat),(A));
1407:   return(0);
1408: }

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

1413:    Logically Collective on Mat

1415:     Input Parameters:
1416: +   mat - the shell matrix
1417: .   f - the function
1418: .   base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated
1419: -   ctx - an optional context for the function

1421:    Output Parameter:
1422: .   flg - PETSC_TRUE if the multiply is likely correct

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

1427:    Level: advanced

1429:    Fortran Notes:
1430:     Not supported from Fortran

1432: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellTestMultTranspose()
1433: @*/
1434: PetscErrorCode  MatShellTestMult(Mat mat,PetscErrorCode (*f)(void*,Vec,Vec),Vec base,void *ctx,PetscBool *flg)
1435: {
1437:   PetscInt       m,n;
1438:   Mat            mf,Dmf,Dmat,Ddiff;
1439:   PetscReal      Diffnorm,Dmfnorm;
1440:   PetscBool      v = PETSC_FALSE, flag = PETSC_TRUE;

1444:   PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_view",&v);
1445:   MatGetLocalSize(mat,&m,&n);
1446:   MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,PETSC_DECIDE,PETSC_DECIDE,&mf);
1447:   MatMFFDSetFunction(mf,f,ctx);
1448:   MatMFFDSetBase(mf,base,NULL);

1450:   MatComputeOperator(mf,MATAIJ,&Dmf);
1451:   MatComputeOperator(mat,MATAIJ,&Dmat);

1453:   MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff);
1454:   MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN);
1455:   MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm);
1456:   MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm);
1457:   if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) {
1458:     flag = PETSC_FALSE;
1459:     if (v) {
1460:       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));
1461:       MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_view");
1462:       MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_view");
1463:       MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_view");
1464:     }
1465:   } else if (v) {
1466:     PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL and matrix free multiple appear to produce the same results\n");
1467:   }
1468:   if (flg) *flg = flag;
1469:   MatDestroy(&Ddiff);
1470:   MatDestroy(&mf);
1471:   MatDestroy(&Dmf);
1472:   MatDestroy(&Dmat);
1473:   return(0);
1474: }

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

1479:    Logically Collective on Mat

1481:     Input Parameters:
1482: +   mat - the shell matrix
1483: .   f - the function
1484: .   base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated
1485: -   ctx - an optional context for the function

1487:    Output Parameter:
1488: .   flg - PETSC_TRUE if the multiply is likely correct

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

1493:    Level: advanced

1495:    Fortran Notes:
1496:     Not supported from Fortran

1498: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellTestMult()
1499: @*/
1500: PetscErrorCode  MatShellTestMultTranspose(Mat mat,PetscErrorCode (*f)(void*,Vec,Vec),Vec base,void *ctx,PetscBool *flg)
1501: {
1503:   Vec            x,y,z;
1504:   PetscInt       m,n,M,N;
1505:   Mat            mf,Dmf,Dmat,Ddiff;
1506:   PetscReal      Diffnorm,Dmfnorm;
1507:   PetscBool      v = PETSC_FALSE, flag = PETSC_TRUE;

1511:   PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_transpose_view",&v);
1512:   MatCreateVecs(mat,&x,&y);
1513:   VecDuplicate(y,&z);
1514:   MatGetLocalSize(mat,&m,&n);
1515:   MatGetSize(mat,&M,&N);
1516:   MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,M,N,&mf);
1517:   MatMFFDSetFunction(mf,f,ctx);
1518:   MatMFFDSetBase(mf,base,NULL);
1519:   MatComputeOperator(mf,MATAIJ,&Dmf);
1520:   MatTranspose(Dmf,MAT_INPLACE_MATRIX,&Dmf);
1521:   MatComputeOperatorTranspose(mat,MATAIJ,&Dmat);

1523:   MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff);
1524:   MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN);
1525:   MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm);
1526:   MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm);
1527:   if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) {
1528:     flag = PETSC_FALSE;
1529:     if (v) {
1530:       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));
1531:       MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");
1532:       MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");
1533:       MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");
1534:     }
1535:   } else if (v) {
1536:     PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL transpose and matrix free multiple appear to produce the same results\n");
1537:   }
1538:   if (flg) *flg = flag;
1539:   MatDestroy(&mf);
1540:   MatDestroy(&Dmat);
1541:   MatDestroy(&Ddiff);
1542:   MatDestroy(&Dmf);
1543:   VecDestroy(&x);
1544:   VecDestroy(&y);
1545:   VecDestroy(&z);
1546:   return(0);
1547: }

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

1552:    Logically Collective on Mat

1554:     Input Parameters:
1555: +   mat - the shell matrix
1556: .   op - the name of the operation
1557: -   g - the function that provides the operation.

1559:    Level: advanced

1561:     Usage:
1562: $      extern PetscErrorCode usermult(Mat,Vec,Vec);
1563: $      MatCreateShell(comm,m,n,M,N,ctx,&A);
1564: $      MatShellSetOperation(A,MATOP_MULT,(void(*)(void))usermult);

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

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

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

1581:     Within each user-defined routine, the user should call
1582:     MatShellGetContext() to obtain the user-defined context that was
1583:     set by MatCreateShell().

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

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

1591: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatSetOperation(), MatShellSetManageScalingShifts()
1592: @*/
1593: PetscErrorCode MatShellSetOperation(Mat mat,MatOperation op,void (*g)(void))
1594: {

1599:   PetscTryMethod(mat,"MatShellSetOperation_C",(Mat,MatOperation,void (*)(void)),(mat,op,g));
1600:   return(0);
1601: }

1603: /*@C
1604:     MatShellGetOperation - Gets a matrix function for a shell matrix.

1606:     Not Collective

1608:     Input Parameters:
1609: +   mat - the shell matrix
1610: -   op - the name of the operation

1612:     Output Parameter:
1613: .   g - the function that provides the operation.

1615:     Level: advanced

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

1623:     All user-provided functions have the same calling
1624:     sequence as the usual matrix interface routines, since they
1625:     are intended to be accessed via the usual matrix interface
1626:     routines, e.g.,
1627: $       MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec)

1629:     Within each user-defined routine, the user should call
1630:     MatShellGetContext() to obtain the user-defined context that was
1631:     set by MatCreateShell().

1633: .seealso: MatCreateShell(), MatShellGetContext(), MatShellSetOperation(), MatShellSetContext()
1634: @*/
1635: PetscErrorCode MatShellGetOperation(Mat mat,MatOperation op,void(**g)(void))
1636: {

1641:   PetscUseMethod(mat,"MatShellGetOperation_C",(Mat,MatOperation,void (**)(void)),(mat,op,g));
1642:   return(0);
1643: }