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:   PetscContainer ctxcontainer;
 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;

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

 89: /* Insert properly diagonally scaled values stored in MatShellPreZeroRight */
 90: static PetscErrorCode MatShellPostZeroLeft(Mat A, Vec x)
 91: {
 92:   Mat_Shell *shell = (Mat_Shell *)A->data;

 94:   if (shell->zrows) {
 95:     VecScatterBegin(shell->zvals_sct_r, shell->zvals_w, x, INSERT_VALUES, SCATTER_REVERSE);
 96:     VecScatterEnd(shell->zvals_sct_r, shell->zvals_w, x, INSERT_VALUES, SCATTER_REVERSE);
 97:   }
 98:   return 0;
 99: }

101: /*
102:      Store and scale values on zeroed rows
103:      xx = [x_1, 0], 0 on zeroed rows
104: */
105: static PetscErrorCode MatShellPreZeroLeft(Mat A, Vec x, Vec *xx)
106: {
107:   Mat_Shell *shell = (Mat_Shell *)A->data;

109:   *xx = NULL;
110:   if (!shell->zrows) {
111:     *xx = x;
112:   } else {
113:     if (!shell->left_work) MatCreateVecs(A, NULL, &shell->left_work);
114:     VecCopy(x, shell->left_work);
115:     VecSet(shell->zvals_w, 0.0);
116:     VecScatterBegin(shell->zvals_sct_r, shell->zvals_w, shell->left_work, INSERT_VALUES, SCATTER_REVERSE);
117:     VecScatterEnd(shell->zvals_sct_r, shell->zvals_w, shell->left_work, INSERT_VALUES, SCATTER_REVERSE);
118:     VecScatterBegin(shell->zvals_sct_r, x, shell->zvals_w, INSERT_VALUES, SCATTER_FORWARD);
119:     VecScatterEnd(shell->zvals_sct_r, x, shell->zvals_w, INSERT_VALUES, SCATTER_FORWARD);
120:     VecPointwiseMult(shell->zvals_w, shell->zvals_w, shell->zvals);
121:     *xx = shell->left_work;
122:   }
123:   return 0;
124: }

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

131:   if (shell->zcols) VecISSet(x, shell->zcols, 0.0);
132:   if (shell->zrows) {
133:     VecScatterBegin(shell->zvals_sct_c, shell->zvals_w, x, ADD_VALUES, SCATTER_REVERSE);
134:     VecScatterEnd(shell->zvals_sct_c, shell->zvals_w, x, ADD_VALUES, SCATTER_REVERSE);
135:   }
136:   return 0;
137: }

139: /*
140:       xx = diag(left)*x
141: */
142: static PetscErrorCode MatShellPreScaleLeft(Mat A, Vec x, Vec *xx)
143: {
144:   Mat_Shell *shell = (Mat_Shell *)A->data;

146:   *xx = NULL;
147:   if (!shell->left) {
148:     *xx = x;
149:   } else {
150:     if (!shell->left_work) VecDuplicate(shell->left, &shell->left_work);
151:     VecPointwiseMult(shell->left_work, x, shell->left);
152:     *xx = shell->left_work;
153:   }
154:   return 0;
155: }

157: /*
158:      xx = diag(right)*x
159: */
160: static PetscErrorCode MatShellPreScaleRight(Mat A, Vec x, Vec *xx)
161: {
162:   Mat_Shell *shell = (Mat_Shell *)A->data;

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

175: /*
176:     x = diag(left)*x
177: */
178: static PetscErrorCode MatShellPostScaleLeft(Mat A, Vec x)
179: {
180:   Mat_Shell *shell = (Mat_Shell *)A->data;

182:   if (shell->left) VecPointwiseMult(x, x, shell->left);
183:   return 0;
184: }

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

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

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

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

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: static PetscErrorCode MatShellGetContext_Shell(Mat mat, void *ctx)
226: {
227:   Mat_Shell *shell = (Mat_Shell *)mat->data;

229:   if (shell->ctxcontainer) PetscContainerGetPointer(shell->ctxcontainer, (void **)ctx);
230:   else *(void **)ctx = NULL;
231:   return 0;
232: }

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

237:     Not Collective

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

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

245:     Level: advanced

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

251: .seealso: `MATSHELL`, `MatCreateShell()`, `MatShellSetOperation()`, `MatShellSetContext()`
252: @*/
253: PetscErrorCode MatShellGetContext(Mat mat, void *ctx)
254: {
257:   PetscUseMethod(mat, "MatShellGetContext_C", (Mat, void *), (mat, ctx));
258:   return 0;
259: }

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

270:   if (!shell->zvals) MatCreateVecs(mat, NULL, &shell->zvals);
271:   if (!shell->zvals_w) VecDuplicate(shell->zvals, &shell->zvals_w);
272:   MatGetOwnershipRange(mat, &rst, NULL);
273:   MatGetOwnershipRangeColumn(mat, &cst, NULL);

275:   /* Expand/create index set of zeroed rows */
276:   PetscMalloc1(nr, &idxs);
277:   for (i = 0; i < nr; i++) idxs[i] = rows[i] + rst;
278:   ISCreateGeneral(PETSC_COMM_SELF, nr, idxs, PETSC_OWN_POINTER, &is1);
279:   ISSort(is1);
280:   VecISSet(shell->zvals, is1, diag);
281:   if (shell->zrows) {
282:     ISSum(shell->zrows, is1, &is2);
283:     ISDestroy(&shell->zrows);
284:     ISDestroy(&is1);
285:     shell->zrows = is2;
286:   } else shell->zrows = is1;

288:   /* Create scatters for diagonal values communications */
289:   VecScatterDestroy(&shell->zvals_sct_c);
290:   VecScatterDestroy(&shell->zvals_sct_r);

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

296:   /* col scatter: from right vector to left vector */
297:   ISGetIndices(shell->zrows, &ridxs);
298:   ISGetLocalSize(shell->zrows, &nr);
299:   PetscMalloc1(nr, &gidxs);
300:   for (i = 0, cum = 0; i < nr; i++) {
301:     if (ridxs[i] >= mat->cmap->N) continue;
302:     gidxs[cum] = ridxs[i];
303:     cum++;
304:   }
305:   ISCreateGeneral(PetscObjectComm((PetscObject)mat), cum, gidxs, PETSC_OWN_POINTER, &is1);
306:   VecScatterCreate(x, is1, shell->zvals_w, is1, &shell->zvals_sct_c);
307:   ISDestroy(&is1);
308:   VecDestroy(&x);
309:   VecDestroy(&b);

311:   /* Expand/create index set of zeroed columns */
312:   if (rc) {
313:     PetscMalloc1(nc, &idxs);
314:     for (i = 0; i < nc; i++) idxs[i] = cols[i] + cst;
315:     ISCreateGeneral(PETSC_COMM_SELF, nc, idxs, PETSC_OWN_POINTER, &is1);
316:     ISSort(is1);
317:     if (shell->zcols) {
318:       ISSum(shell->zcols, is1, &is2);
319:       ISDestroy(&shell->zcols);
320:       ISDestroy(&is1);
321:       shell->zcols = is2;
322:     } else shell->zcols = is1;
323:   }
324:   return 0;
325: }

327: static PetscErrorCode MatZeroRows_Shell(Mat mat, PetscInt n, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
328: {
329:   Mat_Shell *shell = (Mat_Shell *)mat->data;
330:   PetscInt   nr, *lrows;

332:   if (x && b) {
333:     Vec          xt;
334:     PetscScalar *vals;
335:     PetscInt    *gcols, i, st, nl, nc;

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

341:     MatCreateVecs(mat, &xt, NULL);
342:     VecCopy(x, xt);
343:     PetscCalloc1(nc, &vals);
344:     VecSetValues(xt, nc, gcols, vals, INSERT_VALUES); /* xt = [x1, 0] */
345:     PetscFree(vals);
346:     VecAssemblyBegin(xt);
347:     VecAssemblyEnd(xt);
348:     VecAYPX(xt, -1.0, x); /* xt = [0, x2] */

350:     VecGetOwnershipRange(xt, &st, NULL);
351:     VecGetLocalSize(xt, &nl);
352:     VecGetArray(xt, &vals);
353:     for (i = 0; i < nl; i++) {
354:       PetscInt g = i + st;
355:       if (g > mat->rmap->N) continue;
356:       if (PetscAbsScalar(vals[i]) == 0.0) continue;
357:       VecSetValue(b, g, diag * vals[i], INSERT_VALUES);
358:     }
359:     VecRestoreArray(xt, &vals);
360:     VecAssemblyBegin(b);
361:     VecAssemblyEnd(b); /* b  = [b1, x2 * diag] */
362:     VecDestroy(&xt);
363:     PetscFree(gcols);
364:   }
365:   PetscLayoutMapLocal(mat->rmap, n, rows, &nr, &lrows, NULL);
366:   MatZeroRowsColumns_Local_Shell(mat, nr, lrows, 0, NULL, diag, PETSC_FALSE);
367:   if (shell->axpy) MatZeroRows(shell->axpy, n, rows, 0.0, NULL, NULL);
368:   PetscFree(lrows);
369:   return 0;
370: }

372: static PetscErrorCode MatZeroRowsColumns_Shell(Mat mat, PetscInt n, const PetscInt rowscols[], PetscScalar diag, Vec x, Vec b)
373: {
374:   Mat_Shell *shell = (Mat_Shell *)mat->data;
375:   PetscInt  *lrows, *lcols;
376:   PetscInt   nr, nc;
377:   PetscBool  congruent;

379:   if (x && b) {
380:     Vec          xt, bt;
381:     PetscScalar *vals;
382:     PetscInt    *grows, *gcols, i, st, nl;

384:     PetscMalloc2(n, &grows, n, &gcols);
385:     for (i = 0, nr = 0; i < n; i++)
386:       if (rowscols[i] < mat->rmap->N) grows[nr++] = rowscols[i];
387:     for (i = 0, nc = 0; i < n; i++)
388:       if (rowscols[i] < mat->cmap->N) gcols[nc++] = rowscols[i];
389:     PetscCalloc1(n, &vals);

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

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

431:     PetscMalloc1(n, &t);
432:     for (i = 0, nt = 0; i < n; i++)
433:       if (rowscols[i] < mat->cmap->N) t[nt++] = rowscols[i];
434:     PetscLayoutMapLocal(mat->cmap, nt, t, &nc, &lcols, NULL);
435:     PetscFree(t);
436:   }
437:   MatZeroRowsColumns_Local_Shell(mat, nr, lrows, nc, lcols, diag, PETSC_TRUE);
438:   if (!congruent) PetscFree(lcols);
439:   PetscFree(lrows);
440:   if (shell->axpy) MatZeroRowsColumns(shell->axpy, n, rowscols, 0.0, NULL, NULL);
441:   return 0;
442: }

444: static PetscErrorCode MatDestroy_Shell(Mat mat)
445: {
446:   Mat_Shell              *shell = (Mat_Shell *)mat->data;
447:   MatShellMatFunctionList matmat;

449:   if (shell->ops->destroy) (*shell->ops->destroy)(mat);
450:   PetscMemzero(shell->ops, sizeof(struct _MatShellOps));
451:   VecDestroy(&shell->left);
452:   VecDestroy(&shell->right);
453:   VecDestroy(&shell->dshift);
454:   VecDestroy(&shell->left_work);
455:   VecDestroy(&shell->right_work);
456:   VecDestroy(&shell->left_add_work);
457:   VecDestroy(&shell->right_add_work);
458:   VecDestroy(&shell->axpy_left);
459:   VecDestroy(&shell->axpy_right);
460:   MatDestroy(&shell->axpy);
461:   VecDestroy(&shell->zvals_w);
462:   VecDestroy(&shell->zvals);
463:   VecScatterDestroy(&shell->zvals_sct_c);
464:   VecScatterDestroy(&shell->zvals_sct_r);
465:   ISDestroy(&shell->zrows);
466:   ISDestroy(&shell->zcols);

468:   matmat = shell->matmat;
469:   while (matmat) {
470:     MatShellMatFunctionList next = matmat->next;

472:     PetscObjectComposeFunction((PetscObject)mat, matmat->composedname, NULL);
473:     PetscFree(matmat->composedname);
474:     PetscFree(matmat->resultname);
475:     PetscFree(matmat);
476:     matmat = next;
477:   }
478:   MatShellSetContext(mat, NULL);
479:   PetscObjectComposeFunction((PetscObject)mat, "MatShellGetContext_C", NULL);
480:   PetscObjectComposeFunction((PetscObject)mat, "MatShellSetContext_C", NULL);
481:   PetscObjectComposeFunction((PetscObject)mat, "MatShellSetContextDestroy_C", NULL);
482:   PetscObjectComposeFunction((PetscObject)mat, "MatShellSetVecType_C", NULL);
483:   PetscObjectComposeFunction((PetscObject)mat, "MatShellSetManageScalingShifts_C", NULL);
484:   PetscObjectComposeFunction((PetscObject)mat, "MatShellSetOperation_C", NULL);
485:   PetscObjectComposeFunction((PetscObject)mat, "MatShellGetOperation_C", NULL);
486:   PetscObjectComposeFunction((PetscObject)mat, "MatShellSetMatProductOperation_C", NULL);
487:   PetscFree(mat->data);
488:   return 0;
489: }

491: typedef struct {
492:   PetscErrorCode (*numeric)(Mat, Mat, Mat, void *);
493:   PetscErrorCode (*destroy)(void *);
494:   void *userdata;
495:   Mat   B;
496:   Mat   Bt;
497:   Mat   axpy;
498: } MatMatDataShell;

500: static PetscErrorCode DestroyMatMatDataShell(void *data)
501: {
502:   MatMatDataShell *mmdata = (MatMatDataShell *)data;

504:   if (mmdata->destroy) (*mmdata->destroy)(mmdata->userdata);
505:   MatDestroy(&mmdata->B);
506:   MatDestroy(&mmdata->Bt);
507:   MatDestroy(&mmdata->axpy);
508:   PetscFree(mmdata);
509:   return 0;
510: }

512: static PetscErrorCode MatProductNumeric_Shell_X(Mat D)
513: {
514:   Mat_Product     *product;
515:   Mat              A, B;
516:   MatMatDataShell *mdata;
517:   PetscScalar      zero = 0.0;

519:   MatCheckProduct(D, 1);
520:   product = D->product;
522:   A     = product->A;
523:   B     = product->B;
524:   mdata = (MatMatDataShell *)product->data;
525:   if (mdata->numeric) {
526:     Mat_Shell *shell                = (Mat_Shell *)A->data;
527:     PetscErrorCode (*stashsym)(Mat) = D->ops->productsymbolic;
528:     PetscErrorCode (*stashnum)(Mat) = D->ops->productnumeric;
529:     PetscBool useBmdata = PETSC_FALSE, newB = PETSC_TRUE;

531:     if (shell->managescalingshifts) {
533:       if (shell->right || shell->left) {
534:         useBmdata = PETSC_TRUE;
535:         if (!mdata->B) {
536:           MatDuplicate(B, MAT_SHARE_NONZERO_PATTERN, &mdata->B);
537:         } else {
538:           newB = PETSC_FALSE;
539:         }
540:         MatCopy(B, mdata->B, SAME_NONZERO_PATTERN);
541:       }
542:       switch (product->type) {
543:       case MATPRODUCT_AB: /* s L A R B + v L R B + L D R B */
544:         if (shell->right) MatDiagonalScale(mdata->B, shell->right, NULL);
545:         break;
546:       case MATPRODUCT_AtB: /* s R A^t L B + v R L B + R D L B */
547:         if (shell->left) MatDiagonalScale(mdata->B, shell->left, NULL);
548:         break;
549:       case MATPRODUCT_ABt: /* s L A R B^t + v L R B^t + L D R B^t */
550:         if (shell->right) MatDiagonalScale(mdata->B, NULL, shell->right);
551:         break;
552:       case MATPRODUCT_RARt: /* s B L A R B^t + v B L R B^t + B L D R B^t */
553:         if (shell->right && shell->left) {
554:           PetscBool flg;

556:           VecEqual(shell->right, shell->left, &flg);
558:                      ((PetscObject)B)->type_name);
559:         }
560:         if (shell->right) MatDiagonalScale(mdata->B, NULL, shell->right);
561:         break;
562:       case MATPRODUCT_PtAP: /* s B^t L A R B + v B^t L R B + B^t L D R B */
563:         if (shell->right && shell->left) {
564:           PetscBool flg;

566:           VecEqual(shell->right, shell->left, &flg);
568:                      ((PetscObject)B)->type_name);
569:         }
570:         if (shell->right) MatDiagonalScale(mdata->B, shell->right, NULL);
571:         break;
572:       default:
573:         SETERRQ(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);
574:       }
575:     }
576:     /* allow the user to call MatMat operations on D */
577:     D->product              = NULL;
578:     D->ops->productsymbolic = NULL;
579:     D->ops->productnumeric  = NULL;

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

583:     /* clear any leftover user data and restore D pointers */
584:     MatProductClear(D);
585:     D->ops->productsymbolic = stashsym;
586:     D->ops->productnumeric  = stashnum;
587:     D->product              = product;

589:     if (shell->managescalingshifts) {
590:       MatScale(D, shell->vscale);
591:       switch (product->type) {
592:       case MATPRODUCT_AB:  /* s L A R B + v L R B + L D R B */
593:       case MATPRODUCT_ABt: /* s L A R B^t + v L R B^t + L D R B^t */
594:         if (shell->left) {
595:           MatDiagonalScale(D, shell->left, NULL);
596:           if (shell->dshift || shell->vshift != zero) {
597:             if (!shell->left_work) MatCreateVecs(A, NULL, &shell->left_work);
598:             if (shell->dshift) {
599:               VecCopy(shell->dshift, shell->left_work);
600:               VecShift(shell->left_work, shell->vshift);
601:               VecPointwiseMult(shell->left_work, shell->left_work, shell->left);
602:             } else {
603:               VecSet(shell->left_work, shell->vshift);
604:             }
605:             if (product->type == MATPRODUCT_ABt) {
606:               MatReuse     reuse = mdata->Bt ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX;
607:               MatStructure str   = mdata->Bt ? SUBSET_NONZERO_PATTERN : DIFFERENT_NONZERO_PATTERN;

609:               MatTranspose(mdata->B, reuse, &mdata->Bt);
610:               MatDiagonalScale(mdata->Bt, shell->left_work, NULL);
611:               MatAXPY(D, 1.0, mdata->Bt, str);
612:             } else {
613:               MatStructure str = newB ? DIFFERENT_NONZERO_PATTERN : SUBSET_NONZERO_PATTERN;

615:               MatDiagonalScale(mdata->B, shell->left_work, NULL);
616:               MatAXPY(D, 1.0, mdata->B, str);
617:             }
618:           }
619:         }
620:         break;
621:       case MATPRODUCT_AtB: /* s R A^t L B + v R L B + R D L B */
622:         if (shell->right) {
623:           MatDiagonalScale(D, shell->right, NULL);
624:           if (shell->dshift || shell->vshift != zero) {
625:             MatStructure str = newB ? DIFFERENT_NONZERO_PATTERN : SUBSET_NONZERO_PATTERN;

627:             if (!shell->right_work) MatCreateVecs(A, &shell->right_work, NULL);
628:             if (shell->dshift) {
629:               VecCopy(shell->dshift, shell->right_work);
630:               VecShift(shell->right_work, shell->vshift);
631:               VecPointwiseMult(shell->right_work, shell->right_work, shell->right);
632:             } else {
633:               VecSet(shell->right_work, shell->vshift);
634:             }
635:             MatDiagonalScale(mdata->B, shell->right_work, NULL);
636:             MatAXPY(D, 1.0, mdata->B, str);
637:           }
638:         }
639:         break;
640:       case MATPRODUCT_PtAP: /* s B^t L A R B + v B^t L R B + B^t L D R B */
641:       case MATPRODUCT_RARt: /* s B L A R B^t + v B L R B^t + B L D R B^t */
643:                    ((PetscObject)B)->type_name);
644:         break;
645:       default:
646:         SETERRQ(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);
647:       }
648:       if (shell->axpy && shell->axpy_vscale != zero) {
649:         Mat              X;
650:         PetscObjectState axpy_state;
651:         MatStructure     str = DIFFERENT_NONZERO_PATTERN; /* not sure it is safe to ever use SUBSET_NONZERO_PATTERN */

653:         MatShellGetContext(shell->axpy, &X);
654:         PetscObjectStateGet((PetscObject)X, &axpy_state);
656:         if (!mdata->axpy) {
657:           str = DIFFERENT_NONZERO_PATTERN;
658:           MatProductCreate(shell->axpy, B, NULL, &mdata->axpy);
659:           MatProductSetType(mdata->axpy, product->type);
660:           MatProductSetFromOptions(mdata->axpy);
661:           MatProductSymbolic(mdata->axpy);
662:         } else { /* May be that shell->axpy has changed */
663:           PetscBool flg;

665:           MatProductReplaceMats(shell->axpy, B, NULL, mdata->axpy);
666:           MatHasOperation(mdata->axpy, MATOP_PRODUCTSYMBOLIC, &flg);
667:           if (!flg) {
668:             str = DIFFERENT_NONZERO_PATTERN;
669:             MatProductSetFromOptions(mdata->axpy);
670:             MatProductSymbolic(mdata->axpy);
671:           }
672:         }
673:         MatProductNumeric(mdata->axpy);
674:         MatAXPY(D, shell->axpy_vscale, mdata->axpy, str);
675:       }
676:     }
677:   } else SETERRQ(PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "Missing numeric operation");
678:   return 0;
679: }

681: static PetscErrorCode MatProductSymbolic_Shell_X(Mat D)
682: {
683:   Mat_Product            *product;
684:   Mat                     A, B;
685:   MatShellMatFunctionList matmat;
686:   Mat_Shell              *shell;
687:   PetscBool               flg;
688:   char                    composedname[256];
689:   MatMatDataShell        *mdata;

691:   MatCheckProduct(D, 1);
692:   product = D->product;
694:   A      = product->A;
695:   B      = product->B;
696:   shell  = (Mat_Shell *)A->data;
697:   matmat = shell->matmat;
698:   PetscSNPrintf(composedname, sizeof(composedname), "MatProductSetFromOptions_%s_%s_C", ((PetscObject)A)->type_name, ((PetscObject)B)->type_name);
699:   while (matmat) {
700:     PetscStrcmp(composedname, matmat->composedname, &flg);
701:     flg = (PetscBool)(flg && (matmat->ptype == product->type));
702:     if (flg) break;
703:     matmat = matmat->next;
704:   }
706:   switch (product->type) {
707:   case MATPRODUCT_AB:
708:     MatSetSizes(D, A->rmap->n, B->cmap->n, A->rmap->N, B->cmap->N);
709:     break;
710:   case MATPRODUCT_AtB:
711:     MatSetSizes(D, A->cmap->n, B->cmap->n, A->cmap->N, B->cmap->N);
712:     break;
713:   case MATPRODUCT_ABt:
714:     MatSetSizes(D, A->rmap->n, B->rmap->n, A->rmap->N, B->rmap->N);
715:     break;
716:   case MATPRODUCT_RARt:
717:     MatSetSizes(D, B->rmap->n, B->rmap->n, B->rmap->N, B->rmap->N);
718:     break;
719:   case MATPRODUCT_PtAP:
720:     MatSetSizes(D, B->cmap->n, B->cmap->n, B->cmap->N, B->cmap->N);
721:     break;
722:   default:
723:     SETERRQ(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);
724:   }
725:   /* respect users who passed in a matrix for which resultname is the base type */
726:   if (matmat->resultname) {
727:     PetscObjectBaseTypeCompare((PetscObject)D, matmat->resultname, &flg);
728:     if (!flg) MatSetType(D, matmat->resultname);
729:   }
730:   /* If matrix type was not set or different, we need to reset this pointers */
731:   D->ops->productsymbolic = MatProductSymbolic_Shell_X;
732:   D->ops->productnumeric  = MatProductNumeric_Shell_X;
733:   /* attach product data */
734:   PetscNew(&mdata);
735:   mdata->numeric = matmat->numeric;
736:   mdata->destroy = matmat->destroy;
737:   if (matmat->symbolic) {
738:     (*matmat->symbolic)(A, B, D, &mdata->userdata);
739:   } else { /* call general setup if symbolic operation not provided */
740:     MatSetUp(D);
741:   }
744:   D->product->data    = mdata;
745:   D->product->destroy = DestroyMatMatDataShell;
746:   /* Be sure to reset these pointers if the user did something unexpected */
747:   D->ops->productsymbolic = MatProductSymbolic_Shell_X;
748:   D->ops->productnumeric  = MatProductNumeric_Shell_X;
749:   return 0;
750: }

752: static PetscErrorCode MatProductSetFromOptions_Shell_X(Mat D)
753: {
754:   Mat_Product            *product;
755:   Mat                     A, B;
756:   MatShellMatFunctionList matmat;
757:   Mat_Shell              *shell;
758:   PetscBool               flg;
759:   char                    composedname[256];

761:   MatCheckProduct(D, 1);
762:   product = D->product;
763:   A       = product->A;
764:   B       = product->B;
765:   MatIsShell(A, &flg);
766:   if (!flg) return 0;
767:   shell  = (Mat_Shell *)A->data;
768:   matmat = shell->matmat;
769:   PetscSNPrintf(composedname, sizeof(composedname), "MatProductSetFromOptions_%s_%s_C", ((PetscObject)A)->type_name, ((PetscObject)B)->type_name);
770:   while (matmat) {
771:     PetscStrcmp(composedname, matmat->composedname, &flg);
772:     flg = (PetscBool)(flg && (matmat->ptype == product->type));
773:     if (flg) break;
774:     matmat = matmat->next;
775:   }
776:   if (flg) {
777:     D->ops->productsymbolic = MatProductSymbolic_Shell_X;
778:   } else PetscInfo(D, "  symbolic product %s not registered for product type %s\n", composedname, MatProductTypes[product->type]);
779:   return 0;
780: }

782: 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)
783: {
784:   PetscBool               flg;
785:   Mat_Shell              *shell;
786:   MatShellMatFunctionList matmat;


791:   /* add product callback */
792:   shell  = (Mat_Shell *)A->data;
793:   matmat = shell->matmat;
794:   if (!matmat) {
795:     PetscNew(&shell->matmat);
796:     matmat = shell->matmat;
797:   } else {
798:     MatShellMatFunctionList entry = matmat;
799:     while (entry) {
800:       PetscStrcmp(composedname, entry->composedname, &flg);
801:       flg    = (PetscBool)(flg && (entry->ptype == ptype));
802:       matmat = entry;
803:       if (flg) goto set;
804:       entry = entry->next;
805:     }
806:     PetscNew(&matmat->next);
807:     matmat = matmat->next;
808:   }

810: set:
811:   matmat->symbolic = symbolic;
812:   matmat->numeric  = numeric;
813:   matmat->destroy  = destroy;
814:   matmat->ptype    = ptype;
815:   PetscFree(matmat->composedname);
816:   PetscFree(matmat->resultname);
817:   PetscStrallocpy(composedname, &matmat->composedname);
818:   PetscStrallocpy(resultname, &matmat->resultname);
819:   PetscInfo(A, "Composing %s for product type %s with result %s\n", matmat->composedname, MatProductTypes[matmat->ptype], matmat->resultname ? matmat->resultname : "not specified");
820:   PetscObjectComposeFunction((PetscObject)A, matmat->composedname, MatProductSetFromOptions_Shell_X);
821:   return 0;
822: }

824: /*@C
825:     MatShellSetMatProductOperation - Allows user to set a matrix matrix operation for a `MATSHELL` shell matrix.

827:    Logically Collective on A

829:     Input Parameters:
830: +   A - the `MATSHELL` shell matrix
831: .   ptype - the product type
832: .   symbolic - the function for the symbolic phase (can be NULL)
833: .   numeric - the function for the numerical phase
834: .   destroy - the function for the destruction of the needed data generated during the symbolic phase (can be NULL)
835: .   Btype - the matrix type for the matrix to be multiplied against
836: -   Ctype - the matrix type for the result (can be NULL)

838:    Level: advanced

840:     Usage:
841: $      extern PetscErrorCode usersymbolic(Mat,Mat,Mat,void**);
842: $      extern PetscErrorCode usernumeric(Mat,Mat,Mat,void*);
843: $      extern PetscErrorCode userdestroy(void*);
844: $      MatCreateShell(comm,m,n,M,N,ctx,&A);
845: $      MatShellSetMatProductOperation(A,MATPRODUCT_AB,usersymbolic,usernumeric,userdestroy,MATSEQAIJ,MATDENSE);
846: $      [ create B of type SEQAIJ etc..]
847: $      MatProductCreate(A,B,NULL,&C);
848: $      MatProductSetType(C,MATPRODUCT_AB);
849: $      MatProductSetFromOptions(C);
850: $      MatProductSymbolic(C); -> actually runs the user defined symbolic operation
851: $      MatProductNumeric(C); -> actually runs the user defined numeric operation
852: $      [ use C = A*B ]

854:     Notes:
855:     `MATPRODUCT_ABC` is not supported yet. Not supported in Fortran.

857:     If the symbolic phase is not specified, `MatSetUp()` is called on the result matrix that must have its type set if Ctype is NULL.

859:     Any additional data needed by the matrix product needs to be returned during the symbolic phase and destroyed with the destroy callback.
860:     PETSc will take care of calling the user-defined callbacks.
861:     It is allowed to specify the same callbacks for different Btype matrix types.
862:     The couple (Btype,ptype) uniquely identifies the operation: the last specified callbacks takes precedence.

864: .seealso: `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellSetContext()`, `MatSetOperation()`, `MatProductType`, `MatType`, `MatSetUp()`
865: @*/
866: 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)
867: {
874:   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));
875:   return 0;
876: }

878: static 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)
879: {
880:   PetscBool   flg;
881:   char        composedname[256];
882:   MatRootName Bnames = MatRootNameList, Cnames = MatRootNameList;
883:   PetscMPIInt size;

886:   while (Bnames) { /* user passed in the root name */
887:     PetscStrcmp(Btype, Bnames->rname, &flg);
888:     if (flg) break;
889:     Bnames = Bnames->next;
890:   }
891:   while (Cnames) { /* user passed in the root name */
892:     PetscStrcmp(Ctype, Cnames->rname, &flg);
893:     if (flg) break;
894:     Cnames = Cnames->next;
895:   }
896:   MPI_Comm_size(PetscObjectComm((PetscObject)A), &size);
897:   Btype = Bnames ? (size > 1 ? Bnames->mname : Bnames->sname) : Btype;
898:   Ctype = Cnames ? (size > 1 ? Cnames->mname : Cnames->sname) : Ctype;
899:   PetscSNPrintf(composedname, sizeof(composedname), "MatProductSetFromOptions_%s_%s_C", ((PetscObject)A)->type_name, Btype);
900:   MatShellSetMatProductOperation_Private(A, ptype, symbolic, numeric, destroy, composedname, Ctype);
901:   return 0;
902: }

904: static PetscErrorCode MatCopy_Shell(Mat A, Mat B, MatStructure str)
905: {
906:   Mat_Shell              *shellA = (Mat_Shell *)A->data, *shellB = (Mat_Shell *)B->data;
907:   PetscBool               matflg;
908:   MatShellMatFunctionList matmatA;

910:   MatIsShell(B, &matflg);

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

916:   if (shellA->ops->copy) (*shellA->ops->copy)(A, B, str);
917:   shellB->vscale = shellA->vscale;
918:   shellB->vshift = shellA->vshift;
919:   if (shellA->dshift) {
920:     if (!shellB->dshift) VecDuplicate(shellA->dshift, &shellB->dshift);
921:     VecCopy(shellA->dshift, shellB->dshift);
922:   } else {
923:     VecDestroy(&shellB->dshift);
924:   }
925:   if (shellA->left) {
926:     if (!shellB->left) VecDuplicate(shellA->left, &shellB->left);
927:     VecCopy(shellA->left, shellB->left);
928:   } else {
929:     VecDestroy(&shellB->left);
930:   }
931:   if (shellA->right) {
932:     if (!shellB->right) VecDuplicate(shellA->right, &shellB->right);
933:     VecCopy(shellA->right, shellB->right);
934:   } else {
935:     VecDestroy(&shellB->right);
936:   }
937:   MatDestroy(&shellB->axpy);
938:   shellB->axpy_vscale = 0.0;
939:   shellB->axpy_state  = 0;
940:   if (shellA->axpy) {
941:     PetscObjectReference((PetscObject)shellA->axpy);
942:     shellB->axpy        = shellA->axpy;
943:     shellB->axpy_vscale = shellA->axpy_vscale;
944:     shellB->axpy_state  = shellA->axpy_state;
945:   }
946:   if (shellA->zrows) {
947:     ISDuplicate(shellA->zrows, &shellB->zrows);
948:     if (shellA->zcols) ISDuplicate(shellA->zcols, &shellB->zcols);
949:     VecDuplicate(shellA->zvals, &shellB->zvals);
950:     VecCopy(shellA->zvals, shellB->zvals);
951:     VecDuplicate(shellA->zvals_w, &shellB->zvals_w);
952:     PetscObjectReference((PetscObject)shellA->zvals_sct_r);
953:     PetscObjectReference((PetscObject)shellA->zvals_sct_c);
954:     shellB->zvals_sct_r = shellA->zvals_sct_r;
955:     shellB->zvals_sct_c = shellA->zvals_sct_c;
956:   }

958:   matmatA = shellA->matmat;
959:   if (matmatA) {
960:     while (matmatA->next) {
961:       MatShellSetMatProductOperation_Private(B, matmatA->ptype, matmatA->symbolic, matmatA->numeric, matmatA->destroy, matmatA->composedname, matmatA->resultname);
962:       matmatA = matmatA->next;
963:     }
964:   }
965:   return 0;
966: }

968: static PetscErrorCode MatDuplicate_Shell(Mat mat, MatDuplicateOption op, Mat *M)
969: {
970:   MatCreateShell(PetscObjectComm((PetscObject)mat), mat->rmap->n, mat->cmap->n, mat->rmap->N, mat->cmap->N, NULL, M);
971:   ((Mat_Shell *)(*M)->data)->ctxcontainer = ((Mat_Shell *)mat->data)->ctxcontainer;
972:   PetscObjectCompose((PetscObject)(*M), "MatShell ctx", (PetscObject)((Mat_Shell *)(*M)->data)->ctxcontainer);
973:   PetscObjectChangeTypeName((PetscObject)(*M), ((PetscObject)mat)->type_name);
974:   if (op != MAT_DO_NOT_COPY_VALUES) MatCopy(mat, *M, SAME_NONZERO_PATTERN);
975:   return 0;
976: }

978: PetscErrorCode MatMult_Shell(Mat A, Vec x, Vec y)
979: {
980:   Mat_Shell       *shell = (Mat_Shell *)A->data;
981:   Vec              xx;
982:   PetscObjectState instate, outstate;

984:   MatShellPreZeroRight(A, x, &xx);
985:   MatShellPreScaleRight(A, xx, &xx);
986:   PetscObjectStateGet((PetscObject)y, &instate);
987:   (*shell->ops->mult)(A, xx, y);
988:   PetscObjectStateGet((PetscObject)y, &outstate);
989:   if (instate == outstate) {
990:     /* increase the state of the output vector since the user did not update its state themself as should have been done */
991:     PetscObjectStateIncrease((PetscObject)y);
992:   }
993:   MatShellShiftAndScale(A, xx, y);
994:   MatShellPostScaleLeft(A, y);
995:   MatShellPostZeroLeft(A, y);

997:   if (shell->axpy) {
998:     Mat              X;
999:     PetscObjectState axpy_state;

1001:     MatShellGetContext(shell->axpy, &X);
1002:     PetscObjectStateGet((PetscObject)X, &axpy_state);

1005:     MatCreateVecs(shell->axpy, shell->axpy_right ? NULL : &shell->axpy_right, shell->axpy_left ? NULL : &shell->axpy_left);
1006:     VecCopy(x, shell->axpy_right);
1007:     MatMult(shell->axpy, shell->axpy_right, shell->axpy_left);
1008:     VecAXPY(y, shell->axpy_vscale, shell->axpy_left);
1009:   }
1010:   return 0;
1011: }

1013: static PetscErrorCode MatMultAdd_Shell(Mat A, Vec x, Vec y, Vec z)
1014: {
1015:   Mat_Shell *shell = (Mat_Shell *)A->data;

1017:   if (y == z) {
1018:     if (!shell->right_add_work) VecDuplicate(z, &shell->right_add_work);
1019:     MatMult(A, x, shell->right_add_work);
1020:     VecAXPY(z, 1.0, shell->right_add_work);
1021:   } else {
1022:     MatMult(A, x, z);
1023:     VecAXPY(z, 1.0, y);
1024:   }
1025:   return 0;
1026: }

1028: static PetscErrorCode MatMultTranspose_Shell(Mat A, Vec x, Vec y)
1029: {
1030:   Mat_Shell       *shell = (Mat_Shell *)A->data;
1031:   Vec              xx;
1032:   PetscObjectState instate, outstate;

1034:   MatShellPreZeroLeft(A, x, &xx);
1035:   MatShellPreScaleLeft(A, xx, &xx);
1036:   PetscObjectStateGet((PetscObject)y, &instate);
1037:   (*shell->ops->multtranspose)(A, xx, y);
1038:   PetscObjectStateGet((PetscObject)y, &outstate);
1039:   if (instate == outstate) {
1040:     /* increase the state of the output vector since the user did not update its state themself as should have been done */
1041:     PetscObjectStateIncrease((PetscObject)y);
1042:   }
1043:   MatShellShiftAndScale(A, xx, y);
1044:   MatShellPostScaleRight(A, y);
1045:   MatShellPostZeroRight(A, y);

1047:   if (shell->axpy) {
1048:     Mat              X;
1049:     PetscObjectState axpy_state;

1051:     MatShellGetContext(shell->axpy, &X);
1052:     PetscObjectStateGet((PetscObject)X, &axpy_state);
1054:     MatCreateVecs(shell->axpy, shell->axpy_right ? NULL : &shell->axpy_right, shell->axpy_left ? NULL : &shell->axpy_left);
1055:     VecCopy(x, shell->axpy_left);
1056:     MatMultTranspose(shell->axpy, shell->axpy_left, shell->axpy_right);
1057:     VecAXPY(y, shell->axpy_vscale, shell->axpy_right);
1058:   }
1059:   return 0;
1060: }

1062: static PetscErrorCode MatMultTransposeAdd_Shell(Mat A, Vec x, Vec y, Vec z)
1063: {
1064:   Mat_Shell *shell = (Mat_Shell *)A->data;

1066:   if (y == z) {
1067:     if (!shell->left_add_work) VecDuplicate(z, &shell->left_add_work);
1068:     MatMultTranspose(A, x, shell->left_add_work);
1069:     VecAXPY(z, 1.0, shell->left_add_work);
1070:   } else {
1071:     MatMultTranspose(A, x, z);
1072:     VecAXPY(z, 1.0, y);
1073:   }
1074:   return 0;
1075: }

1077: /*
1078:           diag(left)(vscale*A + diag(dshift) + vshift I)diag(right)
1079: */
1080: static PetscErrorCode MatGetDiagonal_Shell(Mat A, Vec v)
1081: {
1082:   Mat_Shell *shell = (Mat_Shell *)A->data;

1084:   if (shell->ops->getdiagonal) {
1085:     (*shell->ops->getdiagonal)(A, v);
1086:   } else SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Must provide shell matrix with routine to return diagonal using\nMatShellSetOperation(S,MATOP_GET_DIAGONAL,...)");
1087:   VecScale(v, shell->vscale);
1088:   if (shell->dshift) VecAXPY(v, 1.0, shell->dshift);
1089:   VecShift(v, shell->vshift);
1090:   if (shell->left) VecPointwiseMult(v, v, shell->left);
1091:   if (shell->right) VecPointwiseMult(v, v, shell->right);
1092:   if (shell->zrows) {
1093:     VecScatterBegin(shell->zvals_sct_r, shell->zvals, v, INSERT_VALUES, SCATTER_REVERSE);
1094:     VecScatterEnd(shell->zvals_sct_r, shell->zvals, v, INSERT_VALUES, SCATTER_REVERSE);
1095:   }
1096:   if (shell->axpy) {
1097:     Mat              X;
1098:     PetscObjectState axpy_state;

1100:     MatShellGetContext(shell->axpy, &X);
1101:     PetscObjectStateGet((PetscObject)X, &axpy_state);
1103:     MatCreateVecs(shell->axpy, NULL, shell->axpy_left ? NULL : &shell->axpy_left);
1104:     MatGetDiagonal(shell->axpy, shell->axpy_left);
1105:     VecAXPY(v, shell->axpy_vscale, shell->axpy_left);
1106:   }
1107:   return 0;
1108: }

1110: static PetscErrorCode MatShift_Shell(Mat Y, PetscScalar a)
1111: {
1112:   Mat_Shell *shell = (Mat_Shell *)Y->data;
1113:   PetscBool  flg;

1115:   MatHasCongruentLayouts(Y, &flg);
1117:   if (shell->left || shell->right) {
1118:     if (!shell->dshift) {
1119:       VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift);
1120:       VecSet(shell->dshift, a);
1121:     } else {
1122:       if (shell->left) VecPointwiseMult(shell->dshift, shell->dshift, shell->left);
1123:       if (shell->right) VecPointwiseMult(shell->dshift, shell->dshift, shell->right);
1124:       VecShift(shell->dshift, a);
1125:     }
1126:     if (shell->left) VecPointwiseDivide(shell->dshift, shell->dshift, shell->left);
1127:     if (shell->right) VecPointwiseDivide(shell->dshift, shell->dshift, shell->right);
1128:   } else shell->vshift += a;
1129:   if (shell->zrows) VecShift(shell->zvals, a);
1130:   return 0;
1131: }

1133: static PetscErrorCode MatDiagonalSet_Shell_Private(Mat A, Vec D, PetscScalar s)
1134: {
1135:   Mat_Shell *shell = (Mat_Shell *)A->data;

1137:   if (!shell->dshift) VecDuplicate(D, &shell->dshift);
1138:   if (shell->left || shell->right) {
1139:     if (!shell->right_work) VecDuplicate(shell->left ? shell->left : shell->right, &shell->right_work);
1140:     if (shell->left && shell->right) {
1141:       VecPointwiseDivide(shell->right_work, D, shell->left);
1142:       VecPointwiseDivide(shell->right_work, shell->right_work, shell->right);
1143:     } else if (shell->left) {
1144:       VecPointwiseDivide(shell->right_work, D, shell->left);
1145:     } else {
1146:       VecPointwiseDivide(shell->right_work, D, shell->right);
1147:     }
1148:     VecAXPY(shell->dshift, s, shell->right_work);
1149:   } else {
1150:     VecAXPY(shell->dshift, s, D);
1151:   }
1152:   return 0;
1153: }

1155: PetscErrorCode MatDiagonalSet_Shell(Mat A, Vec D, InsertMode ins)
1156: {
1157:   Mat_Shell *shell = (Mat_Shell *)A->data;
1158:   Vec        d;
1159:   PetscBool  flg;

1161:   MatHasCongruentLayouts(A, &flg);
1163:   if (ins == INSERT_VALUES) {
1164:     VecDuplicate(D, &d);
1165:     MatGetDiagonal(A, d);
1166:     MatDiagonalSet_Shell_Private(A, d, -1.);
1167:     MatDiagonalSet_Shell_Private(A, D, 1.);
1168:     VecDestroy(&d);
1169:     if (shell->zrows) VecCopy(D, shell->zvals);
1170:   } else {
1171:     MatDiagonalSet_Shell_Private(A, D, 1.);
1172:     if (shell->zrows) VecAXPY(shell->zvals, 1.0, D);
1173:   }
1174:   return 0;
1175: }

1177: static PetscErrorCode MatScale_Shell(Mat Y, PetscScalar a)
1178: {
1179:   Mat_Shell *shell = (Mat_Shell *)Y->data;

1181:   shell->vscale *= a;
1182:   shell->vshift *= a;
1183:   if (shell->dshift) VecScale(shell->dshift, a);
1184:   shell->axpy_vscale *= a;
1185:   if (shell->zrows) VecScale(shell->zvals, a);
1186:   return 0;
1187: }

1189: static PetscErrorCode MatDiagonalScale_Shell(Mat Y, Vec left, Vec right)
1190: {
1191:   Mat_Shell *shell = (Mat_Shell *)Y->data;

1193:   if (left) {
1194:     if (!shell->left) {
1195:       VecDuplicate(left, &shell->left);
1196:       VecCopy(left, shell->left);
1197:     } else {
1198:       VecPointwiseMult(shell->left, shell->left, left);
1199:     }
1200:     if (shell->zrows) VecPointwiseMult(shell->zvals, shell->zvals, left);
1201:   }
1202:   if (right) {
1203:     if (!shell->right) {
1204:       VecDuplicate(right, &shell->right);
1205:       VecCopy(right, shell->right);
1206:     } else {
1207:       VecPointwiseMult(shell->right, shell->right, right);
1208:     }
1209:     if (shell->zrows) {
1210:       if (!shell->left_work) MatCreateVecs(Y, NULL, &shell->left_work);
1211:       VecSet(shell->zvals_w, 1.0);
1212:       VecScatterBegin(shell->zvals_sct_c, right, shell->zvals_w, INSERT_VALUES, SCATTER_FORWARD);
1213:       VecScatterEnd(shell->zvals_sct_c, right, shell->zvals_w, INSERT_VALUES, SCATTER_FORWARD);
1214:       VecPointwiseMult(shell->zvals, shell->zvals, shell->zvals_w);
1215:     }
1216:   }
1217:   if (shell->axpy) MatDiagonalScale(shell->axpy, left, right);
1218:   return 0;
1219: }

1221: PetscErrorCode MatAssemblyEnd_Shell(Mat Y, MatAssemblyType t)
1222: {
1223:   Mat_Shell *shell = (Mat_Shell *)Y->data;

1225:   if (t == MAT_FINAL_ASSEMBLY) {
1226:     shell->vshift      = 0.0;
1227:     shell->vscale      = 1.0;
1228:     shell->axpy_vscale = 0.0;
1229:     shell->axpy_state  = 0;
1230:     VecDestroy(&shell->dshift);
1231:     VecDestroy(&shell->left);
1232:     VecDestroy(&shell->right);
1233:     MatDestroy(&shell->axpy);
1234:     VecDestroy(&shell->axpy_left);
1235:     VecDestroy(&shell->axpy_right);
1236:     VecScatterDestroy(&shell->zvals_sct_c);
1237:     VecScatterDestroy(&shell->zvals_sct_r);
1238:     ISDestroy(&shell->zrows);
1239:     ISDestroy(&shell->zcols);
1240:   }
1241:   return 0;
1242: }

1244: static PetscErrorCode MatMissingDiagonal_Shell(Mat A, PetscBool *missing, PetscInt *d)
1245: {
1246:   *missing = PETSC_FALSE;
1247:   return 0;
1248: }

1250: PetscErrorCode MatAXPY_Shell(Mat Y, PetscScalar a, Mat X, MatStructure str)
1251: {
1252:   Mat_Shell *shell = (Mat_Shell *)Y->data;

1254:   if (X == Y) {
1255:     MatScale(Y, 1.0 + a);
1256:     return 0;
1257:   }
1258:   if (!shell->axpy) {
1259:     MatConvertFrom_Shell(X, MATSHELL, MAT_INITIAL_MATRIX, &shell->axpy);
1260:     shell->axpy_vscale = a;
1261:     PetscObjectStateGet((PetscObject)X, &shell->axpy_state);
1262:   } else {
1263:     MatAXPY(shell->axpy, a / shell->axpy_vscale, X, str);
1264:   }
1265:   return 0;
1266: }

1268: static struct _MatOps MatOps_Values = {NULL,
1269:                                        NULL,
1270:                                        NULL,
1271:                                        NULL,
1272:                                        /* 4*/ MatMultAdd_Shell,
1273:                                        NULL,
1274:                                        MatMultTransposeAdd_Shell,
1275:                                        NULL,
1276:                                        NULL,
1277:                                        NULL,
1278:                                        /*10*/ NULL,
1279:                                        NULL,
1280:                                        NULL,
1281:                                        NULL,
1282:                                        NULL,
1283:                                        /*15*/ NULL,
1284:                                        NULL,
1285:                                        NULL,
1286:                                        MatDiagonalScale_Shell,
1287:                                        NULL,
1288:                                        /*20*/ NULL,
1289:                                        MatAssemblyEnd_Shell,
1290:                                        NULL,
1291:                                        NULL,
1292:                                        /*24*/ MatZeroRows_Shell,
1293:                                        NULL,
1294:                                        NULL,
1295:                                        NULL,
1296:                                        NULL,
1297:                                        /*29*/ NULL,
1298:                                        NULL,
1299:                                        NULL,
1300:                                        NULL,
1301:                                        NULL,
1302:                                        /*34*/ MatDuplicate_Shell,
1303:                                        NULL,
1304:                                        NULL,
1305:                                        NULL,
1306:                                        NULL,
1307:                                        /*39*/ MatAXPY_Shell,
1308:                                        NULL,
1309:                                        NULL,
1310:                                        NULL,
1311:                                        MatCopy_Shell,
1312:                                        /*44*/ NULL,
1313:                                        MatScale_Shell,
1314:                                        MatShift_Shell,
1315:                                        MatDiagonalSet_Shell,
1316:                                        MatZeroRowsColumns_Shell,
1317:                                        /*49*/ NULL,
1318:                                        NULL,
1319:                                        NULL,
1320:                                        NULL,
1321:                                        NULL,
1322:                                        /*54*/ NULL,
1323:                                        NULL,
1324:                                        NULL,
1325:                                        NULL,
1326:                                        NULL,
1327:                                        /*59*/ NULL,
1328:                                        MatDestroy_Shell,
1329:                                        NULL,
1330:                                        MatConvertFrom_Shell,
1331:                                        NULL,
1332:                                        /*64*/ NULL,
1333:                                        NULL,
1334:                                        NULL,
1335:                                        NULL,
1336:                                        NULL,
1337:                                        /*69*/ NULL,
1338:                                        NULL,
1339:                                        MatConvert_Shell,
1340:                                        NULL,
1341:                                        NULL,
1342:                                        /*74*/ NULL,
1343:                                        NULL,
1344:                                        NULL,
1345:                                        NULL,
1346:                                        NULL,
1347:                                        /*79*/ NULL,
1348:                                        NULL,
1349:                                        NULL,
1350:                                        NULL,
1351:                                        NULL,
1352:                                        /*84*/ NULL,
1353:                                        NULL,
1354:                                        NULL,
1355:                                        NULL,
1356:                                        NULL,
1357:                                        /*89*/ NULL,
1358:                                        NULL,
1359:                                        NULL,
1360:                                        NULL,
1361:                                        NULL,
1362:                                        /*94*/ NULL,
1363:                                        NULL,
1364:                                        NULL,
1365:                                        NULL,
1366:                                        NULL,
1367:                                        /*99*/ NULL,
1368:                                        NULL,
1369:                                        NULL,
1370:                                        NULL,
1371:                                        NULL,
1372:                                        /*104*/ NULL,
1373:                                        NULL,
1374:                                        NULL,
1375:                                        NULL,
1376:                                        NULL,
1377:                                        /*109*/ NULL,
1378:                                        NULL,
1379:                                        NULL,
1380:                                        NULL,
1381:                                        MatMissingDiagonal_Shell,
1382:                                        /*114*/ NULL,
1383:                                        NULL,
1384:                                        NULL,
1385:                                        NULL,
1386:                                        NULL,
1387:                                        /*119*/ NULL,
1388:                                        NULL,
1389:                                        NULL,
1390:                                        NULL,
1391:                                        NULL,
1392:                                        /*124*/ NULL,
1393:                                        NULL,
1394:                                        NULL,
1395:                                        NULL,
1396:                                        NULL,
1397:                                        /*129*/ NULL,
1398:                                        NULL,
1399:                                        NULL,
1400:                                        NULL,
1401:                                        NULL,
1402:                                        /*134*/ NULL,
1403:                                        NULL,
1404:                                        NULL,
1405:                                        NULL,
1406:                                        NULL,
1407:                                        /*139*/ NULL,
1408:                                        NULL,
1409:                                        NULL,
1410:                                        NULL,
1411:                                        NULL,
1412:                                        /*144*/ NULL,
1413:                                        NULL,
1414:                                        NULL,
1415:                                        NULL,
1416:                                        NULL,
1417:                                        NULL,
1418:                                        /*150*/ NULL};

1420: static PetscErrorCode MatShellSetContext_Shell(Mat mat, void *ctx)
1421: {
1422:   Mat_Shell *shell = (Mat_Shell *)mat->data;

1424:   if (ctx) {
1425:     PetscContainer ctxcontainer;
1426:     PetscContainerCreate(PetscObjectComm((PetscObject)mat), &ctxcontainer);
1427:     PetscContainerSetPointer(ctxcontainer, ctx);
1428:     PetscObjectCompose((PetscObject)mat, "MatShell ctx", (PetscObject)ctxcontainer);
1429:     shell->ctxcontainer = ctxcontainer;
1430:     PetscContainerDestroy(&ctxcontainer);
1431:   } else {
1432:     PetscObjectCompose((PetscObject)mat, "MatShell ctx", NULL);
1433:     shell->ctxcontainer = NULL;
1434:   }

1436:   return 0;
1437: }

1439: PetscErrorCode MatShellSetContextDestroy_Shell(Mat mat, PetscErrorCode (*f)(void *))
1440: {
1441:   Mat_Shell *shell = (Mat_Shell *)mat->data;

1443:   if (shell->ctxcontainer) PetscContainerSetUserDestroy(shell->ctxcontainer, f);
1444:   return 0;
1445: }

1447: static PetscErrorCode MatShellSetVecType_Shell(Mat mat, VecType vtype)
1448: {
1449:   PetscFree(mat->defaultvectype);
1450:   PetscStrallocpy(vtype, (char **)&mat->defaultvectype);
1451:   return 0;
1452: }

1454: PetscErrorCode MatShellSetManageScalingShifts_Shell(Mat A)
1455: {
1456:   Mat_Shell *shell = (Mat_Shell *)A->data;

1458:   shell->managescalingshifts = PETSC_FALSE;
1459:   A->ops->diagonalset        = NULL;
1460:   A->ops->diagonalscale      = NULL;
1461:   A->ops->scale              = NULL;
1462:   A->ops->shift              = NULL;
1463:   A->ops->axpy               = NULL;
1464:   return 0;
1465: }

1467: static PetscErrorCode MatShellSetOperation_Shell(Mat mat, MatOperation op, void (*f)(void))
1468: {
1469:   Mat_Shell *shell = (Mat_Shell *)mat->data;

1471:   switch (op) {
1472:   case MATOP_DESTROY:
1473:     shell->ops->destroy = (PetscErrorCode(*)(Mat))f;
1474:     break;
1475:   case MATOP_VIEW:
1476:     if (!mat->ops->viewnative) mat->ops->viewnative = mat->ops->view;
1477:     mat->ops->view = (PetscErrorCode(*)(Mat, PetscViewer))f;
1478:     break;
1479:   case MATOP_COPY:
1480:     shell->ops->copy = (PetscErrorCode(*)(Mat, Mat, MatStructure))f;
1481:     break;
1482:   case MATOP_DIAGONAL_SET:
1483:   case MATOP_DIAGONAL_SCALE:
1484:   case MATOP_SHIFT:
1485:   case MATOP_SCALE:
1486:   case MATOP_AXPY:
1487:   case MATOP_ZERO_ROWS:
1488:   case MATOP_ZERO_ROWS_COLUMNS:
1490:     (((void (**)(void))mat->ops)[op]) = f;
1491:     break;
1492:   case MATOP_GET_DIAGONAL:
1493:     if (shell->managescalingshifts) {
1494:       shell->ops->getdiagonal = (PetscErrorCode(*)(Mat, Vec))f;
1495:       mat->ops->getdiagonal   = MatGetDiagonal_Shell;
1496:     } else {
1497:       shell->ops->getdiagonal = NULL;
1498:       mat->ops->getdiagonal   = (PetscErrorCode(*)(Mat, Vec))f;
1499:     }
1500:     break;
1501:   case MATOP_MULT:
1502:     if (shell->managescalingshifts) {
1503:       shell->ops->mult = (PetscErrorCode(*)(Mat, Vec, Vec))f;
1504:       mat->ops->mult   = MatMult_Shell;
1505:     } else {
1506:       shell->ops->mult = NULL;
1507:       mat->ops->mult   = (PetscErrorCode(*)(Mat, Vec, Vec))f;
1508:     }
1509:     break;
1510:   case MATOP_MULT_TRANSPOSE:
1511:     if (shell->managescalingshifts) {
1512:       shell->ops->multtranspose = (PetscErrorCode(*)(Mat, Vec, Vec))f;
1513:       mat->ops->multtranspose   = MatMultTranspose_Shell;
1514:     } else {
1515:       shell->ops->multtranspose = NULL;
1516:       mat->ops->multtranspose   = (PetscErrorCode(*)(Mat, Vec, Vec))f;
1517:     }
1518:     break;
1519:   default:
1520:     (((void (**)(void))mat->ops)[op]) = f;
1521:     break;
1522:   }
1523:   return 0;
1524: }

1526: static PetscErrorCode MatShellGetOperation_Shell(Mat mat, MatOperation op, void (**f)(void))
1527: {
1528:   Mat_Shell *shell = (Mat_Shell *)mat->data;

1530:   switch (op) {
1531:   case MATOP_DESTROY:
1532:     *f = (void (*)(void))shell->ops->destroy;
1533:     break;
1534:   case MATOP_VIEW:
1535:     *f = (void (*)(void))mat->ops->view;
1536:     break;
1537:   case MATOP_COPY:
1538:     *f = (void (*)(void))shell->ops->copy;
1539:     break;
1540:   case MATOP_DIAGONAL_SET:
1541:   case MATOP_DIAGONAL_SCALE:
1542:   case MATOP_SHIFT:
1543:   case MATOP_SCALE:
1544:   case MATOP_AXPY:
1545:   case MATOP_ZERO_ROWS:
1546:   case MATOP_ZERO_ROWS_COLUMNS:
1547:     *f = (((void (**)(void))mat->ops)[op]);
1548:     break;
1549:   case MATOP_GET_DIAGONAL:
1550:     if (shell->ops->getdiagonal) *f = (void (*)(void))shell->ops->getdiagonal;
1551:     else *f = (((void (**)(void))mat->ops)[op]);
1552:     break;
1553:   case MATOP_MULT:
1554:     if (shell->ops->mult) *f = (void (*)(void))shell->ops->mult;
1555:     else *f = (((void (**)(void))mat->ops)[op]);
1556:     break;
1557:   case MATOP_MULT_TRANSPOSE:
1558:     if (shell->ops->multtranspose) *f = (void (*)(void))shell->ops->multtranspose;
1559:     else *f = (((void (**)(void))mat->ops)[op]);
1560:     break;
1561:   default:
1562:     *f = (((void (**)(void))mat->ops)[op]);
1563:   }
1564:   return 0;
1565: }

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

1570:   Level: advanced

1572: .seealso: `Mat`, `MatCreateShell()`
1573: M*/

1575: PETSC_EXTERN PetscErrorCode MatCreate_Shell(Mat A)
1576: {
1577:   Mat_Shell *b;

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

1581:   PetscNew(&b);
1582:   A->data = (void *)b;

1584:   b->ctxcontainer        = NULL;
1585:   b->vshift              = 0.0;
1586:   b->vscale              = 1.0;
1587:   b->managescalingshifts = PETSC_TRUE;
1588:   A->assembled           = PETSC_TRUE;
1589:   A->preallocated        = PETSC_FALSE;

1591:   PetscObjectComposeFunction((PetscObject)A, "MatShellGetContext_C", MatShellGetContext_Shell);
1592:   PetscObjectComposeFunction((PetscObject)A, "MatShellSetContext_C", MatShellSetContext_Shell);
1593:   PetscObjectComposeFunction((PetscObject)A, "MatShellSetContextDestroy_C", MatShellSetContextDestroy_Shell);
1594:   PetscObjectComposeFunction((PetscObject)A, "MatShellSetVecType_C", MatShellSetVecType_Shell);
1595:   PetscObjectComposeFunction((PetscObject)A, "MatShellSetManageScalingShifts_C", MatShellSetManageScalingShifts_Shell);
1596:   PetscObjectComposeFunction((PetscObject)A, "MatShellSetOperation_C", MatShellSetOperation_Shell);
1597:   PetscObjectComposeFunction((PetscObject)A, "MatShellGetOperation_C", MatShellGetOperation_Shell);
1598:   PetscObjectComposeFunction((PetscObject)A, "MatShellSetMatProductOperation_C", MatShellSetMatProductOperation_Shell);
1599:   PetscObjectChangeTypeName((PetscObject)A, MATSHELL);
1600:   return 0;
1601: }

1603: /*@C
1604:    MatCreateShell - Creates a new matrix of `MatType` `MATSHELL` for use with a user-defined
1605:    private data storage format.

1607:   Collective

1609:    Input Parameters:
1610: +  comm - MPI communicator
1611: .  m - number of local rows (must be given)
1612: .  n - number of local columns (must be given)
1613: .  M - number of global rows (may be PETSC_DETERMINE)
1614: .  N - number of global columns (may be PETSC_DETERMINE)
1615: -  ctx - pointer to data needed by the shell matrix routines

1617:    Output Parameter:
1618: .  A - the matrix

1620:    Level: advanced

1622:   Usage:
1623: $    extern PetscErrorCode mult(Mat,Vec,Vec);
1624: $    MatCreateShell(comm,m,n,M,N,ctx,&mat);
1625: $    MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
1626: $    [ Use matrix for operations that have been set ]
1627: $    MatDestroy(mat);

1629:    Notes:
1630:    The shell matrix type is intended to provide a simple class to use
1631:    with `KSP` (such as, for use with matrix-free methods). You should not
1632:    use the shell type if you plan to define a complete matrix class.

1634:    PETSc requires that matrices and vectors being used for certain
1635:    operations are partitioned accordingly.  For example, when
1636:    creating a shell matrix, A, that supports parallel matrix-vector
1637:    products using `MatMult`(A,x,y) the user should set the number
1638:    of local matrix rows to be the number of local elements of the
1639:    corresponding result vector, y. Note that this is information is
1640:    required for use of the matrix interface routines, even though
1641:    the shell matrix may not actually be physically partitioned.
1642:    For example,

1644: $
1645: $     Vec x, y
1646: $     extern PetscErrorCode mult(Mat,Vec,Vec);
1647: $     Mat A
1648: $
1649: $     VecCreateMPI(comm,PETSC_DECIDE,M,&y);
1650: $     VecCreateMPI(comm,PETSC_DECIDE,N,&x);
1651: $     VecGetLocalSize(y,&m);
1652: $     VecGetLocalSize(x,&n);
1653: $     MatCreateShell(comm,m,n,M,N,ctx,&A);
1654: $     MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
1655: $     MatMult(A,x,y);
1656: $     MatDestroy(&A);
1657: $     VecDestroy(&y);
1658: $     VecDestroy(&x);
1659: $

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

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

1666:     Developers Notes:
1667:     Regarding shifting and scaling. The general form is

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

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

1675:           A is the user provided function.

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

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

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

1687:    Fortran Note:
1688:     To use this from Fortran with a ctx you must write an interface definition for this
1689:     function and for `MatShellGetContext()` that tells Fortran the Fortran derived data type you are passing
1690:     in as the ctx argument.

1692: .seealso: `MATSHELL`, `MatShellSetOperation()`, `MatHasOperation()`, `MatShellGetContext()`, `MatShellSetContext()`, `MATSHELL`, `MatShellSetManageScalingShifts()`, `MatShellSetMatProductOperation()`
1693: @*/
1694: PetscErrorCode MatCreateShell(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, PetscInt N, void *ctx, Mat *A)
1695: {
1696:   MatCreate(comm, A);
1697:   MatSetSizes(*A, m, n, M, N);
1698:   MatSetType(*A, MATSHELL);
1699:   MatShellSetContext(*A, ctx);
1700:   MatSetUp(*A);
1701:   return 0;
1702: }

1704: /*@
1705:     MatShellSetContext - sets the context for a `MATSHELL` shell matrix

1707:    Logically Collective on mat

1709:     Input Parameters:
1710: +   mat - the `MATSHELL` shell matrix
1711: -   ctx - the context

1713:    Level: advanced

1715:    Fortran Note:
1716:     To use this from Fortran you must write a Fortran interface definition for this
1717:     function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.

1719: .seealso: `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`
1720: @*/
1721: PetscErrorCode MatShellSetContext(Mat mat, void *ctx)
1722: {
1724:   PetscTryMethod(mat, "MatShellSetContext_C", (Mat, void *), (mat, ctx));
1725:   return 0;
1726: }

1728: /*@
1729:     MatShellSetContextDestroy - sets the destroy function for a `MATSHELL` shell matrix context

1731:    Logically Collective on mat

1733:     Input Parameters:
1734: +   mat - the shell matrix
1735: -   f - the context destroy function

1737:     Note:
1738:     If the `MatShell` is never duplicated, the behavior of this function is equivalent
1739:     to `MatShellSetOperation`(`Mat`,`MATOP_DESTROY`,f). However, `MatShellSetContextDestroy()`
1740:     ensures proper reference counting for the user provided context data in the case that
1741:     the `MATSHELL` is duplicated.

1743:    Level: advanced

1745: .seealso: `MATSHELL`, `MatCreateShell()`, `MatShellSetContext()`
1746: @*/
1747: PetscErrorCode MatShellSetContextDestroy(Mat mat, PetscErrorCode (*f)(void *))
1748: {
1750:   PetscTryMethod(mat, "MatShellSetContextDestroy_C", (Mat, PetscErrorCode(*)(void *)), (mat, f));
1751:   return 0;
1752: }

1754: /*@C
1755:  MatShellSetVecType - Sets the type of `Vec` returned by `MatCreateVecs()`

1757:  Logically collective

1759:     Input Parameters:
1760: +   mat   - the `MATSHELL` shell matrix
1761: -   vtype - type to use for creating vectors

1763:  Level: advanced

1765: .seealso: `MATSHELL`, `MatCreateVecs()`
1766: @*/
1767: PetscErrorCode MatShellSetVecType(Mat mat, VecType vtype)
1768: {
1769:   PetscTryMethod(mat, "MatShellSetVecType_C", (Mat, VecType), (mat, vtype));
1770:   return 0;
1771: }

1773: /*@
1774:     MatShellSetManageScalingShifts - Allows the user to control the scaling and shift operations of the `MATSHELL`. Must be called immediately
1775:           after `MatCreateShell()`

1777:    Logically Collective on A

1779:     Input Parameter:
1780: .   mat - the `MATSHELL` shell matrix

1782:   Level: advanced

1784: .seealso: `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellSetContext()`, `MatShellSetOperation()`
1785: @*/
1786: PetscErrorCode MatShellSetManageScalingShifts(Mat A)
1787: {
1789:   PetscTryMethod(A, "MatShellSetManageScalingShifts_C", (Mat), (A));
1790:   return 0;
1791: }

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

1796:    Logically Collective on mat

1798:     Input Parameters:
1799: +   mat - the `MATSHELL` shell matrix
1800: .   f - the function
1801: .   base - differences are computed around this vector, see `MatMFFDSetBase()`, for Jacobians this is the point at which the Jacobian is being evaluated
1802: -   ctx - an optional context for the function

1804:    Output Parameter:
1805: .   flg - `PETSC_TRUE` if the multiply is likely correct

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

1810:    Level: advanced

1812:    Fortran Note:
1813:     Not supported from Fortran

1815: .seealso: `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMultTranspose()`
1816: @*/
1817: PetscErrorCode MatShellTestMult(Mat mat, PetscErrorCode (*f)(void *, Vec, Vec), Vec base, void *ctx, PetscBool *flg)
1818: {
1819:   PetscInt  m, n;
1820:   Mat       mf, Dmf, Dmat, Ddiff;
1821:   PetscReal Diffnorm, Dmfnorm;
1822:   PetscBool v = PETSC_FALSE, flag = PETSC_TRUE;

1825:   PetscOptionsHasName(NULL, ((PetscObject)mat)->prefix, "-mat_shell_test_mult_view", &v);
1826:   MatGetLocalSize(mat, &m, &n);
1827:   MatCreateMFFD(PetscObjectComm((PetscObject)mat), m, n, PETSC_DECIDE, PETSC_DECIDE, &mf);
1828:   MatMFFDSetFunction(mf, f, ctx);
1829:   MatMFFDSetBase(mf, base, NULL);

1831:   MatComputeOperator(mf, MATAIJ, &Dmf);
1832:   MatComputeOperator(mat, MATAIJ, &Dmat);

1834:   MatDuplicate(Dmat, MAT_COPY_VALUES, &Ddiff);
1835:   MatAXPY(Ddiff, -1.0, Dmf, DIFFERENT_NONZERO_PATTERN);
1836:   MatNorm(Ddiff, NORM_FROBENIUS, &Diffnorm);
1837:   MatNorm(Dmf, NORM_FROBENIUS, &Dmfnorm);
1838:   if (Diffnorm / Dmfnorm > 10 * PETSC_SQRT_MACHINE_EPSILON) {
1839:     flag = PETSC_FALSE;
1840:     if (v) {
1841:       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));
1842:       MatViewFromOptions(Ddiff, (PetscObject)mat, "-mat_shell_test_mult_view");
1843:       MatViewFromOptions(Dmf, (PetscObject)mat, "-mat_shell_test_mult_view");
1844:       MatViewFromOptions(Dmat, (PetscObject)mat, "-mat_shell_test_mult_view");
1845:     }
1846:   } else if (v) {
1847:     PetscPrintf(PetscObjectComm((PetscObject)mat), "MATSHELL and matrix free multiple appear to produce the same results\n");
1848:   }
1849:   if (flg) *flg = flag;
1850:   MatDestroy(&Ddiff);
1851:   MatDestroy(&mf);
1852:   MatDestroy(&Dmf);
1853:   MatDestroy(&Dmat);
1854:   return 0;
1855: }

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

1860:    Logically Collective on mat

1862:     Input Parameters:
1863: +   mat - the `MATSHELL` shell matrix
1864: .   f - the function
1865: .   base - differences are computed around this vector, see `MatMFFDSetBase()`, for Jacobians this is the point at which the Jacobian is being evaluated
1866: -   ctx - an optional context for the function

1868:    Output Parameter:
1869: .   flg - `PETSC_TRUE` if the multiply is likely correct

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

1874:    Level: advanced

1876:    Fortran Note:
1877:     Not supported from Fortran

1879: .seealso: `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMult()`
1880: @*/
1881: PetscErrorCode MatShellTestMultTranspose(Mat mat, PetscErrorCode (*f)(void *, Vec, Vec), Vec base, void *ctx, PetscBool *flg)
1882: {
1883:   Vec       x, y, z;
1884:   PetscInt  m, n, M, N;
1885:   Mat       mf, Dmf, Dmat, Ddiff;
1886:   PetscReal Diffnorm, Dmfnorm;
1887:   PetscBool v = PETSC_FALSE, flag = PETSC_TRUE;

1890:   PetscOptionsHasName(NULL, ((PetscObject)mat)->prefix, "-mat_shell_test_mult_transpose_view", &v);
1891:   MatCreateVecs(mat, &x, &y);
1892:   VecDuplicate(y, &z);
1893:   MatGetLocalSize(mat, &m, &n);
1894:   MatGetSize(mat, &M, &N);
1895:   MatCreateMFFD(PetscObjectComm((PetscObject)mat), m, n, M, N, &mf);
1896:   MatMFFDSetFunction(mf, f, ctx);
1897:   MatMFFDSetBase(mf, base, NULL);
1898:   MatComputeOperator(mf, MATAIJ, &Dmf);
1899:   MatTranspose(Dmf, MAT_INPLACE_MATRIX, &Dmf);
1900:   MatComputeOperatorTranspose(mat, MATAIJ, &Dmat);

1902:   MatDuplicate(Dmat, MAT_COPY_VALUES, &Ddiff);
1903:   MatAXPY(Ddiff, -1.0, Dmf, DIFFERENT_NONZERO_PATTERN);
1904:   MatNorm(Ddiff, NORM_FROBENIUS, &Diffnorm);
1905:   MatNorm(Dmf, NORM_FROBENIUS, &Dmfnorm);
1906:   if (Diffnorm / Dmfnorm > 10 * PETSC_SQRT_MACHINE_EPSILON) {
1907:     flag = PETSC_FALSE;
1908:     if (v) {
1909:       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));
1910:       MatViewFromOptions(Ddiff, (PetscObject)mat, "-mat_shell_test_mult_transpose_view");
1911:       MatViewFromOptions(Dmf, (PetscObject)mat, "-mat_shell_test_mult_transpose_view");
1912:       MatViewFromOptions(Dmat, (PetscObject)mat, "-mat_shell_test_mult_transpose_view");
1913:     }
1914:   } else if (v) {
1915:     PetscPrintf(PetscObjectComm((PetscObject)mat), "MATSHELL transpose and matrix free multiple appear to produce the same results\n");
1916:   }
1917:   if (flg) *flg = flag;
1918:   MatDestroy(&mf);
1919:   MatDestroy(&Dmat);
1920:   MatDestroy(&Ddiff);
1921:   MatDestroy(&Dmf);
1922:   VecDestroy(&x);
1923:   VecDestroy(&y);
1924:   VecDestroy(&z);
1925:   return 0;
1926: }

1928: /*@C
1929:     MatShellSetOperation - Allows user to set a matrix operation for a `MATSHELL` shell matrix.

1931:    Logically Collective on mat

1933:     Input Parameters:
1934: +   mat - the `MATSHELL` shell matrix
1935: .   op - the name of the operation
1936: -   g - the function that provides the operation.

1938:    Level: advanced

1940:     Usage:
1941: $      extern PetscErrorCode usermult(Mat,Vec,Vec);
1942: $      MatCreateShell(comm,m,n,M,N,ctx,&A);
1943: $      MatShellSetOperation(A,MATOP_MULT,(void(*)(void))usermult);

1945:     Notes:
1946:     See the file include/petscmat.h for a complete list of matrix
1947:     operations, which all have the form MATOP_<OPERATION>, where
1948:     <OPERATION> is the name (in all capital letters) of the
1949:     user interface routine (e.g., `MatMult()` -> `MATOP_MULT`).

1951:     All user-provided functions (except for `MATOP_DESTROY`) should have the same calling
1952:     sequence as the usual matrix interface routines, since they
1953:     are intended to be accessed via the usual matrix interface
1954:     routines, e.g.,
1955: $       MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec)

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

1960:     Within each user-defined routine, the user should call
1961:     `MatShellGetContext()` to obtain the user-defined context that was
1962:     set by `MatCreateShell()`.

1964:     Use `MatSetOperation()` to set an operation for any matrix type. For matrix product operations (i.e. `MatMatXXX()`, `MatTransposeMatXXX()` etc) use `MatShellSetMatProductOperation()`

1966:     Fortran Note:
1967:     For `MatCreateVecs()` the user code should check if the input left or right matrix is -1 and in that case not
1968:        generate a matrix. See src/mat/tests/ex120f.F

1970: .seealso: `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellSetContext()`, `MatSetOperation()`, `MatShellSetManageScalingShifts()`, `MatShellSetMatProductOperation()`
1971: @*/
1972: PetscErrorCode MatShellSetOperation(Mat mat, MatOperation op, void (*g)(void))
1973: {
1975:   PetscTryMethod(mat, "MatShellSetOperation_C", (Mat, MatOperation, void (*)(void)), (mat, op, g));
1976:   return 0;
1977: }

1979: /*@C
1980:     MatShellGetOperation - Gets a matrix function for a `MATSHELL` shell matrix.

1982:     Not Collective

1984:     Input Parameters:
1985: +   mat - the `MATSHELL` shell matrix
1986: -   op - the name of the operation

1988:     Output Parameter:
1989: .   g - the function that provides the operation.

1991:     Level: advanced

1993:     Note:
1994:     See the file include/petscmat.h for a complete list of matrix
1995:     operations, which all have the form MATOP_<OPERATION>, where
1996:     <OPERATION> is the name (in all capital letters) of the
1997:     user interface routine (e.g., `MatMult()` -> `MATOP_MULT`).

1999:     All user-provided functions have the same calling
2000:     sequence as the usual matrix interface routines, since they
2001:     are intended to be accessed via the usual matrix interface
2002:     routines, e.g.,
2003: $       MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec)

2005:     Within each user-defined routine, the user should call
2006:     `MatShellGetContext()` to obtain the user-defined context that was
2007:     set by `MatCreateShell()`.

2009: .seealso: `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellSetOperation()`, `MatShellSetContext()`
2010: @*/
2011: PetscErrorCode MatShellGetOperation(Mat mat, MatOperation op, void (**g)(void))
2012: {
2014:   PetscUseMethod(mat, "MatShellGetOperation_C", (Mat, MatOperation, void (**)(void)), (mat, op, g));
2015:   return 0;
2016: }

2018: /*@
2019:     MatIsShell - Inquires if a matrix is derived from `MATSHELL`

2021:     Input Parameter:
2022: .   mat - the matrix

2024:     Output Parameter:
2025: .   flg - the boolean value

2027:     Level: developer

2029:     Developer Note:
2030:     In the future, we should allow the object type name to be changed still using the `MATSHELL` data structure for other matrices (i.e. `MATTRANSPOSEVIRTUAL`, `MATSCHURCOMPLEMENT` etc)

2032: .seealso: `MATSHELL`, `MATMFFD`, `MatCreateShell()`, `MATTRANSPOSEVIRTUAL`, `MATSCHURCOMPLEMENT`
2033: @*/
2034: PetscErrorCode MatIsShell(Mat mat, PetscBool *flg)
2035: {
2038:   *flg = (PetscBool)(mat->ops->destroy == MatDestroy_Shell);
2039:   return 0;
2040: }