Actual source code: shell.c


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

  8: #include <petsc/private/matimpl.h>
  9: #include <petsc/private/vecimpl.h>

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

259:     Not Collective

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

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

267:     Level: advanced

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

703:         MatShellGetContext(shell->axpy,&X);
704:         PetscObjectStateGet((PetscObject)X,&axpy_state);
705:         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,...)");
706:         if (!mdata->axpy) {
707:           str  = DIFFERENT_NONZERO_PATTERN;
708:           MatProductCreate(shell->axpy,B,NULL,&mdata->axpy);
709:           MatProductSetType(mdata->axpy,product->type);
710:           MatProductSetFromOptions(mdata->axpy);
711:           MatProductSymbolic(mdata->axpy);
712:         } else { /* May be that shell->axpy has changed */
713:           PetscBool flg;

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

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

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

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

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

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

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

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

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

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

883:    Logically Collective on Mat

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

894:    Level: advanced

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1080:     MatShellGetContext(shell->axpy,&X);
1081:     PetscObjectStateGet((PetscObject)X,&axpy_state);
1082:     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,...)");

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

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

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

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

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

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

1135:     MatShellGetContext(shell->axpy,&X);
1136:     PetscObjectStateGet((PetscObject)X,&axpy_state);
1137:     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,...)");
1138:     MatCreateVecs(shell->axpy,shell->axpy_right ? NULL : &shell->axpy_right,shell->axpy_left ? NULL : &shell->axpy_left);
1139:     VecCopy(x,shell->axpy_left);
1140:     MatMultTranspose(shell->axpy,shell->axpy_left,shell->axpy_right);
1141:     VecAXPY(y,shell->axpy_vscale,shell->axpy_right);
1142:   }
1143:   return(0);
1144: }

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

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

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

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

1190:     MatShellGetContext(shell->axpy,&X);
1191:     PetscObjectStateGet((PetscObject)X,&axpy_state);
1192:     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,...)");
1193:     MatCreateVecs(shell->axpy,NULL,shell->axpy_left ? NULL : &shell->axpy_left);
1194:     MatGetDiagonal(shell->axpy,shell->axpy_left);
1195:     VecAXPY(v,shell->axpy_vscale,shell->axpy_left);
1196:   }
1197:   return(0);
1198: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1680:   Level: advanced

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

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

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

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

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

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

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

1718:   Collective

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

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

1731:    Level: advanced

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

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

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

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

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

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

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

1782:     Developers Notes:
1783:     Regarding shifting and scaling. The general form is

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

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

1791:           A is the user provided function.

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

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

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

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

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

1818: /*@
1819:     MatShellSetContext - sets the context for a shell matrix

1821:    Logically Collective on Mat

1823:     Input Parameters:
1824: +   mat - the shell matrix
1825: -   ctx - the context

1827:    Level: advanced

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

1833: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation()
1834: @*/
1835: PetscErrorCode  MatShellSetContext(Mat mat,void *ctx)
1836: {

1841:   PetscTryMethod(mat,"MatShellSetContext_C",(Mat,void*),(mat,ctx));
1842:   return(0);
1843: }

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

1848:  Logically collective

1850:     Input Parameters:
1851: +   mat   - the shell matrix
1852: -   vtype - type to use for creating vectors

1854:  Notes:

1856:  Level: advanced

1858: .seealso: MatCreateVecs()
1859: @*/
1860: PetscErrorCode  MatShellSetVecType(Mat mat,VecType vtype)
1861: {

1865:   PetscTryMethod(mat,"MatShellSetVecType_C",(Mat,VecType),(mat,vtype));
1866:   return(0);
1867: }

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

1873:    Logically Collective on Mat

1875:     Input Parameter:
1876: .   mat - the shell matrix

1878:   Level: advanced

1880: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatShellSetOperation()
1881: @*/
1882: PetscErrorCode MatShellSetManageScalingShifts(Mat A)
1883: {

1888:   PetscTryMethod(A,"MatShellSetManageScalingShifts_C",(Mat),(A));
1889:   return(0);
1890: }

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

1895:    Logically Collective on Mat

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

1903:    Output Parameter:
1904: .   flg - PETSC_TRUE if the multiply is likely correct

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

1909:    Level: advanced

1911:    Fortran Notes:
1912:     Not supported from Fortran

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

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

1932:   MatComputeOperator(mf,MATAIJ,&Dmf);
1933:   MatComputeOperator(mat,MATAIJ,&Dmat);

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

1958: /*@C
1959:     MatShellTestMultTranspose - Compares the multiply transpose routine provided to the MATSHELL with differencing on a given function.

1961:    Logically Collective on Mat

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

1969:    Output Parameter:
1970: .   flg - PETSC_TRUE if the multiply is likely correct

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

1975:    Level: advanced

1977:    Fortran Notes:
1978:     Not supported from Fortran

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

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

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

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

2034:    Logically Collective on Mat

2036:     Input Parameters:
2037: +   mat - the shell matrix
2038: .   op - the name of the operation
2039: -   g - the function that provides the operation.

2041:    Level: advanced

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

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

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

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

2063:     Within each user-defined routine, the user should call
2064:     MatShellGetContext() to obtain the user-defined context that was
2065:     set by MatCreateShell().

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

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

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

2081:   PetscTryMethod(mat,"MatShellSetOperation_C",(Mat,MatOperation,void (*)(void)),(mat,op,g));
2082:   return(0);
2083: }

2085: /*@C
2086:     MatShellGetOperation - Gets a matrix function for a shell matrix.

2088:     Not Collective

2090:     Input Parameters:
2091: +   mat - the shell matrix
2092: -   op - the name of the operation

2094:     Output Parameter:
2095: .   g - the function that provides the operation.

2097:     Level: advanced

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

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

2111:     Within each user-defined routine, the user should call
2112:     MatShellGetContext() to obtain the user-defined context that was
2113:     set by MatCreateShell().

2115: .seealso: MatCreateShell(), MatShellGetContext(), MatShellSetOperation(), MatShellSetContext()
2116: @*/
2117: PetscErrorCode MatShellGetOperation(Mat mat,MatOperation op,void(**g)(void))
2118: {

2123:   PetscUseMethod(mat,"MatShellGetOperation_C",(Mat,MatOperation,void (**)(void)),(mat,op,g));
2124:   return(0);
2125: }

2127: /*@
2128:     MatIsShell - Inquires if a matrix is derived from MATSHELL

2130:     Input Parameter:
2131: .   mat - the matrix

2133:     Output Parameter:
2134: .   flg - the boolean value

2136:     Level: developer

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

2140: .seealso: MatCreateShell()
2141: @*/
2142: PetscErrorCode MatIsShell(Mat mat, PetscBool *flg)
2143: {
2147:   *flg = (PetscBool)(mat->ops->destroy == MatDestroy_Shell);
2148:   return(0);
2149: }