Actual source code: shell.c

petsc-master 2019-11-13
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:   Mat         axpy;
 28:   PetscScalar axpy_vscale;
 29:   PetscBool   managescalingshifts;                   /* The user will manage the scaling and shifts for the MATSHELL, not the default */
 30:   /* support for ZeroRows/Columns operations */
 31:   IS          zrows;
 32:   IS          zcols;
 33:   Vec         zvals;
 34:   Vec         zvals_w;
 35:   VecScatter  zvals_sct_r;
 36:   VecScatter  zvals_sct_c;
 37:   void        *ctx;
 38: } Mat_Shell;


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

195: /*
196:          Y = vscale*Y + diag(dshift)*X + vshift*X

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

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

225: PetscErrorCode  MatShellGetContext_Shell(Mat mat,void *ctx)
226: {
228:   *(void**)ctx = ((Mat_Shell*)(mat->data))->ctx;
229:   return(0);
230: }

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

235:     Not Collective

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

240:     Output Parameter:
241: .   ctx - the user provided context

243:     Level: advanced

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

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

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

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

273:   if (shell->axpy) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot diagonal scale MATSHELL after MatAXPY operation");
274:   if (!shell->zvals) {
275:     MatCreateVecs(mat,NULL,&shell->zvals);
276:   }
277:   if (!shell->zvals_w) {
278:     VecDuplicate(shell->zvals,&shell->zvals_w);
279:   }
280:   MatGetOwnershipRange(mat,&rst,NULL);
281:   MatGetOwnershipRangeColumn(mat,&cst,NULL);

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

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

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

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

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

335: static PetscErrorCode MatZeroRows_Shell(Mat mat,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
336: {
337:   PetscInt       nr, *lrows;

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

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

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

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

379: static PetscErrorCode MatZeroRowsColumns_Shell(Mat mat,PetscInt n,const PetscInt rowscols[],PetscScalar diag,Vec x,Vec b)
380: {
381:   PetscInt       *lrows, *lcols;
382:   PetscInt       nr, nc;
383:   PetscBool      congruent;

387:   if (x && b) {
388:     Vec          xt, bt;
389:     PetscScalar *vals;
390:     PetscInt    *grows,*gcols,i,st,nl;

392:     PetscMalloc2(n,&grows,n,&gcols);
393:     for (i = 0, nr = 0; i < n; i++) if (rowscols[i] < mat->rmap->N) grows[nr++] = rowscols[i];
394:     for (i = 0, nc = 0; i < n; i++) if (rowscols[i] < mat->cmap->N) gcols[nc++] = rowscols[i];
395:     PetscCalloc1(n,&vals);

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

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

437:     PetscMalloc1(n,&t);
438:     for (i = 0, nt = 0; i < n; i++) if (rowscols[i] < mat->cmap->N) t[nt++] = rowscols[i];
439:     PetscLayoutMapLocal(mat->cmap,nt,t,&nc,&lcols,NULL);
440:     PetscFree(t);
441:   }
442:   MatZeroRowsColumns_Local_Shell(mat,nr,lrows,nc,lcols,diag,PETSC_TRUE);
443:   if (!congruent) {
444:     PetscFree(lcols);
445:   }
446:   PetscFree(lrows);
447:   return(0);
448: }

450: PetscErrorCode MatDestroy_Shell(Mat mat)
451: {
453:   Mat_Shell      *shell = (Mat_Shell*)mat->data;

456:   if (shell->ops->destroy) {
457:     (*shell->ops->destroy)(mat);
458:   }
459:   VecDestroy(&shell->left);
460:   VecDestroy(&shell->right);
461:   VecDestroy(&shell->dshift);
462:   VecDestroy(&shell->left_work);
463:   VecDestroy(&shell->right_work);
464:   VecDestroy(&shell->left_add_work);
465:   VecDestroy(&shell->right_add_work);
466:   MatDestroy(&shell->axpy);
467:   VecDestroy(&shell->zvals_w);
468:   VecDestroy(&shell->zvals);
469:   VecScatterDestroy(&shell->zvals_sct_c);
470:   VecScatterDestroy(&shell->zvals_sct_r);
471:   ISDestroy(&shell->zrows);
472:   ISDestroy(&shell->zcols);
473:   PetscFree(mat->data);
474:   PetscObjectComposeFunction((PetscObject)mat,"MatShellGetContext_C",NULL);
475:   PetscObjectComposeFunction((PetscObject)mat,"MatShellSetContext_C",NULL);
476:   PetscObjectComposeFunction((PetscObject)mat,"MatShellSetVecType_C",NULL);
477:   PetscObjectComposeFunction((PetscObject)mat,"MatShellSetManageScalingShifts_C",NULL);
478:   PetscObjectComposeFunction((PetscObject)mat,"MatShellSetOperation_C",NULL);
479:   PetscObjectComposeFunction((PetscObject)mat,"MatShellGetOperation_C",NULL);
480:   return(0);
481: }

483: PetscErrorCode MatCopy_Shell(Mat A,Mat B,MatStructure str)
484: {
485:   Mat_Shell       *shellA = (Mat_Shell*)A->data,*shellB = (Mat_Shell*)B->data;
486:   PetscErrorCode  ierr;
487:   PetscBool       matflg;

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

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

496:   if (shellA->ops->copy) {
497:     (*shellA->ops->copy)(A,B,str);
498:   }
499:   shellB->vscale = shellA->vscale;
500:   shellB->vshift = shellA->vshift;
501:   if (shellA->dshift) {
502:     if (!shellB->dshift) {
503:       VecDuplicate(shellA->dshift,&shellB->dshift);
504:     }
505:     VecCopy(shellA->dshift,shellB->dshift);
506:   } else {
507:     VecDestroy(&shellB->dshift);
508:   }
509:   if (shellA->left) {
510:     if (!shellB->left) {
511:       VecDuplicate(shellA->left,&shellB->left);
512:     }
513:     VecCopy(shellA->left,shellB->left);
514:   } else {
515:     VecDestroy(&shellB->left);
516:   }
517:   if (shellA->right) {
518:     if (!shellB->right) {
519:       VecDuplicate(shellA->right,&shellB->right);
520:     }
521:     VecCopy(shellA->right,shellB->right);
522:   } else {
523:     VecDestroy(&shellB->right);
524:   }
525:   MatDestroy(&shellB->axpy);
526:   if (shellA->axpy) {
527:     PetscObjectReference((PetscObject)shellA->axpy);
528:     shellB->axpy        = shellA->axpy;
529:     shellB->axpy_vscale = shellA->axpy_vscale;
530:   }
531:   if (shellA->zrows) {
532:     ISDuplicate(shellA->zrows,&shellB->zrows);
533:     if (shellA->zcols) {
534:       ISDuplicate(shellA->zcols,&shellB->zcols);
535:     }
536:     VecDuplicate(shellA->zvals,&shellB->zvals);
537:     VecCopy(shellA->zvals,shellB->zvals);
538:     VecDuplicate(shellA->zvals_w,&shellB->zvals_w);
539:     PetscObjectReference((PetscObject)shellA->zvals_sct_r);
540:     PetscObjectReference((PetscObject)shellA->zvals_sct_c);
541:     shellB->zvals_sct_r = shellA->zvals_sct_r;
542:     shellB->zvals_sct_c = shellA->zvals_sct_c;
543:   }
544:   return(0);
545: }

547: PetscErrorCode MatDuplicate_Shell(Mat mat,MatDuplicateOption op,Mat *M)
548: {
550:   void           *ctx;

553:   MatShellGetContext(mat,&ctx);
554:   MatCreateShell(PetscObjectComm((PetscObject)mat),mat->rmap->n,mat->cmap->n,mat->rmap->N,mat->cmap->N,ctx,M);
555:   if (op != MAT_DO_NOT_COPY_VALUES) {
556:     MatCopy(mat,*M,SAME_NONZERO_PATTERN);
557:   }
558:   return(0);
559: }

561: PetscErrorCode MatMult_Shell(Mat A,Vec x,Vec y)
562: {
563:   Mat_Shell        *shell = (Mat_Shell*)A->data;
564:   PetscErrorCode   ierr;
565:   Vec              xx;
566:   PetscObjectState instate,outstate;

569:   if (!shell->ops->mult) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Have not provided a MatMult() for this MATSHELL");
570:   MatShellPreZeroRight(A,x,&xx);
571:   MatShellPreScaleRight(A,xx,&xx);
572:   PetscObjectStateGet((PetscObject)y, &instate);
573:   (*shell->ops->mult)(A,xx,y);
574:   PetscObjectStateGet((PetscObject)y, &outstate);
575:   if (instate == outstate) {
576:     /* increase the state of the output vector since the user did not update its state themself as should have been done */
577:     PetscObjectStateIncrease((PetscObject)y);
578:   }
579:   MatShellShiftAndScale(A,xx,y);
580:   MatShellPostScaleLeft(A,y);
581:   MatShellPostZeroLeft(A,y);

583:   if (shell->axpy) {
584:     if (!shell->left_work) {MatCreateVecs(A,&shell->left_work,NULL);}
585:     MatMult(shell->axpy,x,shell->left_work);
586:     VecAXPY(y,shell->axpy_vscale,shell->left_work);
587:   }
588:   return(0);
589: }

591: PetscErrorCode MatMultAdd_Shell(Mat A,Vec x,Vec y,Vec z)
592: {
593:   Mat_Shell      *shell = (Mat_Shell*)A->data;

597:   if (y == z) {
598:     if (!shell->right_add_work) {VecDuplicate(z,&shell->right_add_work);}
599:     MatMult(A,x,shell->right_add_work);
600:     VecAXPY(z,1.0,shell->right_add_work);
601:   } else {
602:     MatMult(A,x,z);
603:     VecAXPY(z,1.0,y);
604:   }
605:   return(0);
606: }

608: PetscErrorCode MatMultTranspose_Shell(Mat A,Vec x,Vec y)
609: {
610:   Mat_Shell        *shell = (Mat_Shell*)A->data;
611:   PetscErrorCode   ierr;
612:   Vec              xx;
613:   PetscObjectState instate,outstate;

616:   MatShellPreZeroLeft(A,x,&xx);
617:   MatShellPreScaleLeft(A,xx,&xx);
618:   PetscObjectStateGet((PetscObject)y, &instate);
619:   if (!shell->ops->multtranspose) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Have not provided a MatMultTranspose() for this MATSHELL");
620:   (*shell->ops->multtranspose)(A,xx,y);
621:   PetscObjectStateGet((PetscObject)y, &outstate);
622:   if (instate == outstate) {
623:     /* increase the state of the output vector since the user did not update its state themself as should have been done */
624:     PetscObjectStateIncrease((PetscObject)y);
625:   }
626:   MatShellShiftAndScale(A,xx,y);
627:   MatShellPostScaleRight(A,y);
628:   MatShellPostZeroRight(A,y);

630:   if (shell->axpy) {
631:     if (!shell->right_work) {MatCreateVecs(A,NULL,&shell->right_work);}
632:     MatMultTranspose(shell->axpy,x,shell->right_work);
633:     VecAXPY(y,shell->axpy_vscale,shell->right_work);
634:   }
635:   return(0);
636: }

638: PetscErrorCode MatMultTransposeAdd_Shell(Mat A,Vec x,Vec y,Vec z)
639: {
640:   Mat_Shell      *shell = (Mat_Shell*)A->data;

644:   if (y == z) {
645:     if (!shell->left_add_work) {VecDuplicate(z,&shell->left_add_work);}
646:     MatMultTranspose(A,x,shell->left_add_work);
647:     VecAXPY(z,1.0,shell->left_add_work);
648:   } else {
649:     MatMultTranspose(A,x,z);
650:     VecAXPY(z,1.0,y);
651:   }
652:   return(0);
653: }

655: /*
656:           diag(left)(vscale*A + diag(dshift) + vshift I)diag(right)
657: */
658: PetscErrorCode MatGetDiagonal_Shell(Mat A,Vec v)
659: {
660:   Mat_Shell      *shell = (Mat_Shell*)A->data;

664:   if (shell->ops->getdiagonal) {
665:     (*shell->ops->getdiagonal)(A,v);
666:   } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Must provide shell matrix with routine to return diagonal using\nMatShellSetOperation(S,MATOP_GET_DIAGONAL,...)");
667:   VecScale(v,shell->vscale);
668:   if (shell->dshift) {
669:     VecAXPY(v,1.0,shell->dshift);
670:   }
671:   VecShift(v,shell->vshift);
672:   if (shell->left)  {VecPointwiseMult(v,v,shell->left);}
673:   if (shell->right) {VecPointwiseMult(v,v,shell->right);}
674:   if (shell->zrows) {
675:     VecScatterBegin(shell->zvals_sct_r,shell->zvals,v,INSERT_VALUES,SCATTER_REVERSE);
676:     VecScatterEnd(shell->zvals_sct_r,shell->zvals,v,INSERT_VALUES,SCATTER_REVERSE);
677:   }
678:   if (shell->axpy) {
679:     if (!shell->left_work) {VecDuplicate(v,&shell->left_work);}
680:     MatGetDiagonal(shell->axpy,shell->left_work);
681:     VecAXPY(v,shell->axpy_vscale,shell->left_work);
682:   }
683:   return(0);
684: }

686: PetscErrorCode MatShift_Shell(Mat Y,PetscScalar a)
687: {
688:   Mat_Shell      *shell = (Mat_Shell*)Y->data;
690:   PetscBool      flg;

693:   MatHasCongruentLayouts(Y,&flg);
694:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_SUP,"Cannot shift shell matrix if it is not congruent");
695:   if (shell->left || shell->right) {
696:     if (!shell->dshift) {
697:       VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift);
698:       VecSet(shell->dshift,a);
699:     } else {
700:       if (shell->left)  {VecPointwiseMult(shell->dshift,shell->dshift,shell->left);}
701:       if (shell->right) {VecPointwiseMult(shell->dshift,shell->dshift,shell->right);}
702:       VecShift(shell->dshift,a);
703:     }
704:     if (shell->left)  {VecPointwiseDivide(shell->dshift,shell->dshift,shell->left);}
705:     if (shell->right) {VecPointwiseDivide(shell->dshift,shell->dshift,shell->right);}
706:   } else shell->vshift += a;
707:   if (shell->zrows) {
708:     VecShift(shell->zvals,a);
709:   }
710:   return(0);
711: }

713: PetscErrorCode MatDiagonalSet_Shell_Private(Mat A,Vec D,PetscScalar s)
714: {
715:   Mat_Shell      *shell = (Mat_Shell*)A->data;

719:   if (!shell->dshift) {VecDuplicate(D,&shell->dshift);}
720:   if (shell->left || shell->right) {
721:     if (!shell->right_work) {VecDuplicate(shell->left ? shell->left : shell->right, &shell->right_work);}
722:     if (shell->left && shell->right)  {
723:       VecPointwiseDivide(shell->right_work,D,shell->left);
724:       VecPointwiseDivide(shell->right_work,shell->right_work,shell->right);
725:     } else if (shell->left) {
726:       VecPointwiseDivide(shell->right_work,D,shell->left);
727:     } else {
728:       VecPointwiseDivide(shell->right_work,D,shell->right);
729:     }
730:     VecAXPY(shell->dshift,s,shell->right_work);
731:   } else {
732:     VecAXPY(shell->dshift,s,D);
733:   }
734:   return(0);
735: }

737: PetscErrorCode MatDiagonalSet_Shell(Mat A,Vec D,InsertMode ins)
738: {
739:   Mat_Shell      *shell = (Mat_Shell*)A->data;
740:   Vec            d;
742:   PetscBool      flg;

745:   MatHasCongruentLayouts(A,&flg);
746:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Cannot diagonal set or shift shell matrix if it is not congruent");
747:   if (ins == INSERT_VALUES) {
748:     if (!A->ops->getdiagonal) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Operation MATOP_GETDIAGONAL must be set first");
749:     VecDuplicate(D,&d);
750:     MatGetDiagonal(A,d);
751:     MatDiagonalSet_Shell_Private(A,d,-1.);
752:     MatDiagonalSet_Shell_Private(A,D,1.);
753:     VecDestroy(&d);
754:     if (shell->zrows) {
755:       VecCopy(D,shell->zvals);
756:     }
757:   } else {
758:     MatDiagonalSet_Shell_Private(A,D,1.);
759:     if (shell->zrows) {
760:       VecAXPY(shell->zvals,1.0,D);
761:     }
762:   }
763:   return(0);
764: }

766: PetscErrorCode MatScale_Shell(Mat Y,PetscScalar a)
767: {
768:   Mat_Shell      *shell = (Mat_Shell*)Y->data;

772:   shell->vscale *= a;
773:   shell->vshift *= a;
774:   if (shell->dshift) {
775:     VecScale(shell->dshift,a);
776:   }
777:   shell->axpy_vscale *= a;
778:   if (shell->zrows) {
779:     VecScale(shell->zvals,a);
780:   }
781:   return(0);
782: }

784: static PetscErrorCode MatDiagonalScale_Shell(Mat Y,Vec left,Vec right)
785: {
786:   Mat_Shell      *shell = (Mat_Shell*)Y->data;

790:   if (shell->axpy) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_SUP,"Cannot diagonal scale MATSHELL after MatAXPY operation");
791:   if (left) {
792:     if (!shell->left) {
793:       VecDuplicate(left,&shell->left);
794:       VecCopy(left,shell->left);
795:     } else {
796:       VecPointwiseMult(shell->left,shell->left,left);
797:     }
798:     if (shell->zrows) {
799:       VecPointwiseMult(shell->zvals,shell->zvals,left);
800:     }
801:   }
802:   if (right) {
803:     if (!shell->right) {
804:       VecDuplicate(right,&shell->right);
805:       VecCopy(right,shell->right);
806:     } else {
807:       VecPointwiseMult(shell->right,shell->right,right);
808:     }
809:     if (shell->zrows) {
810:       if (!shell->left_work) {
811:         MatCreateVecs(Y,NULL,&shell->left_work);
812:       }
813:       VecSet(shell->zvals_w,1.0);
814:       VecScatterBegin(shell->zvals_sct_c,right,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);
815:       VecScatterEnd(shell->zvals_sct_c,right,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);
816:       VecPointwiseMult(shell->zvals,shell->zvals,shell->zvals_w);
817:     }
818:   }
819:   return(0);
820: }

822: PetscErrorCode MatAssemblyEnd_Shell(Mat Y,MatAssemblyType t)
823: {
824:   Mat_Shell      *shell = (Mat_Shell*)Y->data;

828:   if (t == MAT_FINAL_ASSEMBLY) {
829:     shell->vshift = 0.0;
830:     shell->vscale = 1.0;
831:     VecDestroy(&shell->dshift);
832:     VecDestroy(&shell->left);
833:     VecDestroy(&shell->right);
834:     MatDestroy(&shell->axpy);
835:     VecScatterDestroy(&shell->zvals_sct_c);
836:     VecScatterDestroy(&shell->zvals_sct_r);
837:     ISDestroy(&shell->zrows);
838:     ISDestroy(&shell->zcols);
839:   }
840:   return(0);
841: }

843: static PetscErrorCode MatMissingDiagonal_Shell(Mat A,PetscBool  *missing,PetscInt *d)
844: {
846:   *missing = PETSC_FALSE;
847:   return(0);
848: }

850: PetscErrorCode MatAXPY_Shell(Mat Y,PetscScalar a,Mat X,MatStructure str)
851: {
852:   Mat_Shell      *shell = (Mat_Shell*)Y->data;

856:   PetscObjectReference((PetscObject)X);
857:   MatDestroy(&shell->axpy);
858:   shell->axpy        = X;
859:   shell->axpy_vscale = a;
860:   return(0);
861: }

863: static struct _MatOps MatOps_Values = {0,
864:                                        0,
865:                                        0,
866:                                        0,
867:                                 /* 4*/ MatMultAdd_Shell,
868:                                        0,
869:                                        MatMultTransposeAdd_Shell,
870:                                        0,
871:                                        0,
872:                                        0,
873:                                 /*10*/ 0,
874:                                        0,
875:                                        0,
876:                                        0,
877:                                        0,
878:                                 /*15*/ 0,
879:                                        0,
880:                                        0,
881:                                        MatDiagonalScale_Shell,
882:                                        0,
883:                                 /*20*/ 0,
884:                                        MatAssemblyEnd_Shell,
885:                                        0,
886:                                        0,
887:                                 /*24*/ MatZeroRows_Shell,
888:                                        0,
889:                                        0,
890:                                        0,
891:                                        0,
892:                                 /*29*/ 0,
893:                                        0,
894:                                        0,
895:                                        0,
896:                                        0,
897:                                 /*34*/ MatDuplicate_Shell,
898:                                        0,
899:                                        0,
900:                                        0,
901:                                        0,
902:                                 /*39*/ MatAXPY_Shell,
903:                                        0,
904:                                        0,
905:                                        0,
906:                                        MatCopy_Shell,
907:                                 /*44*/ 0,
908:                                        MatScale_Shell,
909:                                        MatShift_Shell,
910:                                        MatDiagonalSet_Shell,
911:                                        MatZeroRowsColumns_Shell,
912:                                 /*49*/ 0,
913:                                        0,
914:                                        0,
915:                                        0,
916:                                        0,
917:                                 /*54*/ 0,
918:                                        0,
919:                                        0,
920:                                        0,
921:                                        0,
922:                                 /*59*/ 0,
923:                                        MatDestroy_Shell,
924:                                        0,
925:                                        MatConvertFrom_Shell,
926:                                        0,
927:                                 /*64*/ 0,
928:                                        0,
929:                                        0,
930:                                        0,
931:                                        0,
932:                                 /*69*/ 0,
933:                                        0,
934:                                        MatConvert_Shell,
935:                                        0,
936:                                        0,
937:                                 /*74*/ 0,
938:                                        0,
939:                                        0,
940:                                        0,
941:                                        0,
942:                                 /*79*/ 0,
943:                                        0,
944:                                        0,
945:                                        0,
946:                                        0,
947:                                 /*84*/ 0,
948:                                        0,
949:                                        0,
950:                                        0,
951:                                        0,
952:                                 /*89*/ 0,
953:                                        0,
954:                                        0,
955:                                        0,
956:                                        0,
957:                                 /*94*/ 0,
958:                                        0,
959:                                        0,
960:                                        0,
961:                                        0,
962:                                 /*99*/ 0,
963:                                        0,
964:                                        0,
965:                                        0,
966:                                        0,
967:                                /*104*/ 0,
968:                                        0,
969:                                        0,
970:                                        0,
971:                                        0,
972:                                /*109*/ 0,
973:                                        0,
974:                                        0,
975:                                        0,
976:                                        MatMissingDiagonal_Shell,
977:                                /*114*/ 0,
978:                                        0,
979:                                        0,
980:                                        0,
981:                                        0,
982:                                /*119*/ 0,
983:                                        0,
984:                                        0,
985:                                        0,
986:                                        0,
987:                                /*124*/ 0,
988:                                        0,
989:                                        0,
990:                                        0,
991:                                        0,
992:                                /*129*/ 0,
993:                                        0,
994:                                        0,
995:                                        0,
996:                                        0,
997:                                /*134*/ 0,
998:                                        0,
999:                                        0,
1000:                                        0,
1001:                                        0,
1002:                                /*139*/ 0,
1003:                                        0,
1004:                                        0
1005: };

1007: PetscErrorCode  MatShellSetContext_Shell(Mat mat,void *ctx)
1008: {
1009:   Mat_Shell      *shell = (Mat_Shell*)mat->data;

1012:   shell->ctx = ctx;
1013:   return(0);
1014: }

1016: static PetscErrorCode MatShellSetVecType_Shell(Mat mat,VecType vtype)
1017: {

1021:   PetscFree(mat->defaultvectype);
1022:   PetscStrallocpy(vtype,(char**)&mat->defaultvectype);
1023:   return(0);
1024: }

1026: PetscErrorCode MatShellSetManageScalingShifts_Shell(Mat A)
1027: {
1028:   Mat_Shell      *shell = (Mat_Shell*)A->data;

1031:   shell->managescalingshifts = PETSC_FALSE;
1032:   A->ops->diagonalset   = NULL;
1033:   A->ops->diagonalscale = NULL;
1034:   A->ops->scale         = NULL;
1035:   A->ops->shift         = NULL;
1036:   A->ops->axpy          = NULL;
1037:   return(0);
1038: }

1040: PetscErrorCode MatShellSetOperation_Shell(Mat mat,MatOperation op,void (*f)(void))
1041: {
1042:   Mat_Shell      *shell = (Mat_Shell*)mat->data;

1045:   switch (op) {
1046:   case MATOP_DESTROY:
1047:     shell->ops->destroy = (PetscErrorCode (*)(Mat))f;
1048:     break;
1049:   case MATOP_VIEW:
1050:     if (!mat->ops->viewnative) {
1051:       mat->ops->viewnative = mat->ops->view;
1052:     }
1053:     mat->ops->view = (PetscErrorCode (*)(Mat,PetscViewer))f;
1054:     break;
1055:   case MATOP_COPY:
1056:     shell->ops->copy = (PetscErrorCode (*)(Mat,Mat,MatStructure))f;
1057:     break;
1058:   case MATOP_DIAGONAL_SET:
1059:   case MATOP_DIAGONAL_SCALE:
1060:   case MATOP_SHIFT:
1061:   case MATOP_SCALE:
1062:   case MATOP_AXPY:
1063:   case MATOP_ZERO_ROWS:
1064:   case MATOP_ZERO_ROWS_COLUMNS:
1065:     if (shell->managescalingshifts) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"MATSHELL is managing scalings and shifts, see MatShellSetManageScalingShifts()");
1066:     (((void(**)(void))mat->ops)[op]) = f;
1067:     break;
1068:   case MATOP_GET_DIAGONAL:
1069:     if (shell->managescalingshifts) {
1070:       shell->ops->getdiagonal = (PetscErrorCode (*)(Mat,Vec))f;
1071:       mat->ops->getdiagonal   = MatGetDiagonal_Shell;
1072:     } else {
1073:       shell->ops->getdiagonal = NULL;
1074:       mat->ops->getdiagonal   = (PetscErrorCode (*)(Mat,Vec))f;
1075:     }
1076:     break;
1077:   case MATOP_MULT:
1078:     if (shell->managescalingshifts) {
1079:       shell->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1080:       mat->ops->mult   = MatMult_Shell;
1081:     } else {
1082:       shell->ops->mult = NULL;
1083:       mat->ops->mult   = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1084:     }
1085:     break;
1086:   case MATOP_MULT_TRANSPOSE:
1087:     if (shell->managescalingshifts) {
1088:       shell->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1089:       mat->ops->multtranspose   = MatMultTranspose_Shell;
1090:     } else {
1091:       shell->ops->multtranspose = NULL;
1092:       mat->ops->multtranspose   = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1093:     }
1094:     break;
1095:   default:
1096:     (((void(**)(void))mat->ops)[op]) = f;
1097:     break;
1098:   }
1099:   return(0);
1100: }

1102: PetscErrorCode MatShellGetOperation_Shell(Mat mat,MatOperation op,void(**f)(void))
1103: {
1104:   Mat_Shell      *shell = (Mat_Shell*)mat->data;

1107:   switch (op) {
1108:   case MATOP_DESTROY:
1109:     *f = (void (*)(void))shell->ops->destroy;
1110:     break;
1111:   case MATOP_VIEW:
1112:     *f = (void (*)(void))mat->ops->view;
1113:     break;
1114:   case MATOP_COPY:
1115:     *f = (void (*)(void))shell->ops->copy;
1116:     break;
1117:   case MATOP_DIAGONAL_SET:
1118:   case MATOP_DIAGONAL_SCALE:
1119:   case MATOP_SHIFT:
1120:   case MATOP_SCALE:
1121:   case MATOP_AXPY:
1122:   case MATOP_ZERO_ROWS:
1123:   case MATOP_ZERO_ROWS_COLUMNS:
1124:     *f = (((void (**)(void))mat->ops)[op]);
1125:     break;
1126:   case MATOP_GET_DIAGONAL:
1127:     if (shell->ops->getdiagonal)
1128:       *f = (void (*)(void))shell->ops->getdiagonal;
1129:     else
1130:       *f = (((void (**)(void))mat->ops)[op]);
1131:     break;
1132:   case MATOP_MULT:
1133:     if (shell->ops->mult)
1134:       *f = (void (*)(void))shell->ops->mult;
1135:     else
1136:       *f = (((void (**)(void))mat->ops)[op]);
1137:     break;
1138:   case MATOP_MULT_TRANSPOSE:
1139:     if (shell->ops->multtranspose)
1140:       *f = (void (*)(void))shell->ops->multtranspose;
1141:     else
1142:       *f = (((void (**)(void))mat->ops)[op]);
1143:     break;
1144:   default:
1145:     *f = (((void (**)(void))mat->ops)[op]);
1146:   }
1147:   return(0);
1148: }

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

1153:   Level: advanced

1155: .seealso: MatCreateShell()
1156: M*/

1158: PETSC_EXTERN PetscErrorCode MatCreate_Shell(Mat A)
1159: {
1160:   Mat_Shell      *b;

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

1166:   PetscNewLog(A,&b);
1167:   A->data = (void*)b;

1169:   PetscLayoutSetUp(A->rmap);
1170:   PetscLayoutSetUp(A->cmap);

1172:   b->ctx                 = 0;
1173:   b->vshift              = 0.0;
1174:   b->vscale              = 1.0;
1175:   b->managescalingshifts = PETSC_TRUE;
1176:   A->assembled           = PETSC_TRUE;
1177:   A->preallocated        = PETSC_FALSE;

1179:   PetscObjectComposeFunction((PetscObject)A,"MatShellGetContext_C",MatShellGetContext_Shell);
1180:   PetscObjectComposeFunction((PetscObject)A,"MatShellSetContext_C",MatShellSetContext_Shell);
1181:   PetscObjectComposeFunction((PetscObject)A,"MatShellSetVecType_C",MatShellSetVecType_Shell);
1182:   PetscObjectComposeFunction((PetscObject)A,"MatShellSetManageScalingShifts_C",MatShellSetManageScalingShifts_Shell);
1183:   PetscObjectComposeFunction((PetscObject)A,"MatShellSetOperation_C",MatShellSetOperation_Shell);
1184:   PetscObjectComposeFunction((PetscObject)A,"MatShellGetOperation_C",MatShellGetOperation_Shell);
1185:   PetscObjectChangeTypeName((PetscObject)A,MATSHELL);
1186:   return(0);
1187: }

1189: /*@C
1190:    MatCreateShell - Creates a new matrix class for use with a user-defined
1191:    private data storage format.

1193:   Collective

1195:    Input Parameters:
1196: +  comm - MPI communicator
1197: .  m - number of local rows (must be given)
1198: .  n - number of local columns (must be given)
1199: .  M - number of global rows (may be PETSC_DETERMINE)
1200: .  N - number of global columns (may be PETSC_DETERMINE)
1201: -  ctx - pointer to data needed by the shell matrix routines

1203:    Output Parameter:
1204: .  A - the matrix

1206:    Level: advanced

1208:   Usage:
1209: $    extern PetscErrorCode mult(Mat,Vec,Vec);
1210: $    MatCreateShell(comm,m,n,M,N,ctx,&mat);
1211: $    MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
1212: $    [ Use matrix for operations that have been set ]
1213: $    MatDestroy(mat);

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

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

1225:    PETSc requires that matrices and vectors being used for certain
1226:    operations are partitioned accordingly.  For example, when
1227:    creating a shell matrix, A, that supports parallel matrix-vector
1228:    products using MatMult(A,x,y) the user should set the number
1229:    of local matrix rows to be the number of local elements of the
1230:    corresponding result vector, y. Note that this is information is
1231:    required for use of the matrix interface routines, even though
1232:    the shell matrix may not actually be physically partitioned.
1233:    For example,

1235: $
1236: $     Vec x, y
1237: $     extern PetscErrorCode mult(Mat,Vec,Vec);
1238: $     Mat A
1239: $
1240: $     VecCreateMPI(comm,PETSC_DECIDE,M,&y);
1241: $     VecCreateMPI(comm,PETSC_DECIDE,N,&x);
1242: $     VecGetLocalSize(y,&m);
1243: $     VecGetLocalSize(x,&n);
1244: $     MatCreateShell(comm,m,n,M,N,ctx,&A);
1245: $     MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
1246: $     MatMult(A,x,y);
1247: $     MatDestroy(&A);
1248: $     VecDestroy(&y);
1249: $     VecDestroy(&x);
1250: $


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


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

1259:     Developers Notes:
1260:     Regarding shifting and scaling. The general form is

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

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

1268:           A is the user provided function.

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

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

1278: .seealso: MatShellSetOperation(), MatHasOperation(), MatShellGetContext(), MatShellSetContext(), MATSHELL, MatShellSetManageScalingShifts()
1279: @*/
1280: PetscErrorCode  MatCreateShell(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,void *ctx,Mat *A)
1281: {

1285:   MatCreate(comm,A);
1286:   MatSetSizes(*A,m,n,M,N);
1287:   MatSetType(*A,MATSHELL);
1288:   MatShellSetContext(*A,ctx);
1289:   MatSetUp(*A);
1290:   return(0);
1291: }


1294: /*@
1295:     MatShellSetContext - sets the context for a shell matrix

1297:    Logically Collective on Mat

1299:     Input Parameters:
1300: +   mat - the shell matrix
1301: -   ctx - the context

1303:    Level: advanced

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

1309: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation()
1310: @*/
1311: PetscErrorCode  MatShellSetContext(Mat mat,void *ctx)
1312: {

1317:   PetscUseMethod(mat,"MatShellSetContext_C",(Mat,void*),(mat,ctx));
1318:   return(0);
1319: }

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

1324:  Logically collective

1326:     Input Parameters:
1327: +   mat   - the shell matrix
1328: -   vtype - type to use for creating vectors

1330:  Notes:

1332:  Level: advanced

1334: .seealso: MatCreateVecs()
1335: @*/
1336: PetscErrorCode  MatShellSetVecType(Mat mat,VecType vtype)
1337: {

1341:   PetscTryMethod(mat,"MatShellSetVecType_C",(Mat,VecType),(mat,vtype));
1342:   return(0);
1343: }

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

1349:    Logically Collective on Mat

1351:     Input Parameter:
1352: .   mat - the shell matrix

1354:   Level: advanced

1356: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatShellSetOperation()
1357: @*/
1358: PetscErrorCode MatShellSetManageScalingShifts(Mat A)
1359: {

1364:   PetscUseMethod(A,"MatShellSetManageScalingShifts_C",(Mat),(A));
1365:   return(0);
1366: }

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

1371:    Logically Collective on Mat

1373:     Input Parameters:
1374: +   mat - the shell matrix
1375: .   f - the function
1376: .   base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated
1377: -   ctx - an optional context for the function

1379:    Output Parameter:
1380: .   flg - PETSC_TRUE if the multiply is likely correct

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

1385:    Level: advanced

1387:    Fortran Notes:
1388:     Not supported from Fortran

1390: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellTestMultTranspose()
1391: @*/
1392: PetscErrorCode  MatShellTestMult(Mat mat,PetscErrorCode (*f)(void*,Vec,Vec),Vec base,void *ctx,PetscBool *flg)
1393: {
1395:   PetscInt       m,n;
1396:   Mat            mf,Dmf,Dmat,Ddiff;
1397:   PetscReal      Diffnorm,Dmfnorm;
1398:   PetscBool      v = PETSC_FALSE, flag = PETSC_TRUE;

1402:   PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_view",&v);
1403:   MatGetLocalSize(mat,&m,&n);
1404:   MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,PETSC_DECIDE,PETSC_DECIDE,&mf);
1405:   MatMFFDSetFunction(mf,f,ctx);
1406:   MatMFFDSetBase(mf,base,NULL);

1408:   MatComputeOperator(mf,MATAIJ,&Dmf);
1409:   MatComputeOperator(mat,MATAIJ,&Dmat);

1411:   MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff);
1412:   MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN);
1413:   MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm);
1414:   MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm);
1415:   if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) {
1416:     flag = PETSC_FALSE;
1417:     if (v) {
1418:       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));
1419:       MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_view");
1420:       MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_view");
1421:       MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_view");
1422:     }
1423:   } else if (v) {
1424:     PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL and matrix free multiple appear to produce the same results\n");
1425:   }
1426:   if (flg) *flg = flag;
1427:   MatDestroy(&Ddiff);
1428:   MatDestroy(&mf);
1429:   MatDestroy(&Dmf);
1430:   MatDestroy(&Dmat);
1431:   return(0);
1432: }

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

1437:    Logically Collective on Mat

1439:     Input Parameters:
1440: +   mat - the shell matrix
1441: .   f - the function
1442: .   base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated
1443: -   ctx - an optional context for the function

1445:    Output Parameter:
1446: .   flg - PETSC_TRUE if the multiply is likely correct

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

1451:    Level: advanced

1453:    Fortran Notes:
1454:     Not supported from Fortran

1456: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellTestMult()
1457: @*/
1458: PetscErrorCode  MatShellTestMultTranspose(Mat mat,PetscErrorCode (*f)(void*,Vec,Vec),Vec base,void *ctx,PetscBool *flg)
1459: {
1461:   Vec            x,y,z;
1462:   PetscInt       m,n,M,N;
1463:   Mat            mf,Dmf,Dmat,Ddiff;
1464:   PetscReal      Diffnorm,Dmfnorm;
1465:   PetscBool      v = PETSC_FALSE, flag = PETSC_TRUE;

1469:   PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_transpose_view",&v);
1470:   MatCreateVecs(mat,&x,&y);
1471:   VecDuplicate(y,&z);
1472:   MatGetLocalSize(mat,&m,&n);
1473:   MatGetSize(mat,&M,&N);
1474:   MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,M,N,&mf);
1475:   MatMFFDSetFunction(mf,f,ctx);
1476:   MatMFFDSetBase(mf,base,NULL);
1477:   MatComputeOperator(mf,MATAIJ,&Dmf);
1478:   MatTranspose(Dmf,MAT_INPLACE_MATRIX,&Dmf);
1479:   MatComputeOperatorTranspose(mat,MATAIJ,&Dmat);

1481:   MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff);
1482:   MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN);
1483:   MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm);
1484:   MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm);
1485:   if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) {
1486:     flag = PETSC_FALSE;
1487:     if (v) {
1488:       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));
1489:       MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");
1490:       MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");
1491:       MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");
1492:     }
1493:   } else if (v) {
1494:     PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL transpose and matrix free multiple appear to produce the same results\n");
1495:   }
1496:   if (flg) *flg = flag;
1497:   MatDestroy(&mf);
1498:   MatDestroy(&Dmat);
1499:   MatDestroy(&Ddiff);
1500:   MatDestroy(&Dmf);
1501:   VecDestroy(&x);
1502:   VecDestroy(&y);
1503:   VecDestroy(&z);
1504:   return(0);
1505: }

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

1510:    Logically Collective on Mat

1512:     Input Parameters:
1513: +   mat - the shell matrix
1514: .   op - the name of the operation
1515: -   g - the function that provides the operation.

1517:    Level: advanced

1519:     Usage:
1520: $      extern PetscErrorCode usermult(Mat,Vec,Vec);
1521: $      MatCreateShell(comm,m,n,M,N,ctx,&A);
1522: $      MatShellSetOperation(A,MATOP_MULT,(void(*)(void))usermult);

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

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

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

1539:     Within each user-defined routine, the user should call
1540:     MatShellGetContext() to obtain the user-defined context that was
1541:     set by MatCreateShell().

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

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

1549: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatSetOperation(), MatShellSetManageScalingShifts()
1550: @*/
1551: PetscErrorCode MatShellSetOperation(Mat mat,MatOperation op,void (*g)(void))
1552: {

1557:   PetscUseMethod(mat,"MatShellSetOperation_C",(Mat,MatOperation,void (*)(void)),(mat,op,g));
1558:   return(0);
1559: }

1561: /*@C
1562:     MatShellGetOperation - Gets a matrix function for a shell matrix.

1564:     Not Collective

1566:     Input Parameters:
1567: +   mat - the shell matrix
1568: -   op - the name of the operation

1570:     Output Parameter:
1571: .   g - the function that provides the operation.

1573:     Level: advanced

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

1581:     All user-provided functions have the same calling
1582:     sequence as the usual matrix interface routines, since they
1583:     are intended to be accessed via the usual matrix interface
1584:     routines, e.g.,
1585: $       MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec)

1587:     Within each user-defined routine, the user should call
1588:     MatShellGetContext() to obtain the user-defined context that was
1589:     set by MatCreateShell().

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

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