Actual source code: shell.c

petsc-3.14.6 2021-03-30
Report Typos and Errors

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

260:     Not Collective

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

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

268:     Level: advanced

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

704:         MatShellGetContext(shell->axpy,(void *)&X);
705:         PetscObjectStateGet((PetscObject)X,&axpy_state);
706:         if (shell->axpy_state != axpy_state) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Invalid AXPY state: cannot modify the X matrix passed to MatAXPY(Y,a,X,...)");
707:         if (!mdata->axpy) {
708:           str  = DIFFERENT_NONZERO_PATTERN;
709:           MatProductCreate(shell->axpy,B,NULL,&mdata->axpy);
710:           MatProductSetType(mdata->axpy,product->type);
711:           MatProductSetFromOptions(mdata->axpy);
712:           MatProductSymbolic(mdata->axpy);
713:         } else { /* May be that shell->axpy has changed */
714:           PetscBool flg;

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

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

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

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

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

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

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

848:   /* add product callback */
849:   shell = (Mat_Shell*)A->data;
850:   matmat = shell->matmat;
851:   if (!matmat) {
852:     PetscNew(&shell->matmat);
853:     matmat = shell->matmat;
854:   } else {
855:     MatShellMatFunctionList entry = matmat;
856:     while (entry) {
857:       PetscStrcmp(composedname,entry->composedname,&flg);
858:       flg  = (PetscBool)(flg && (entry->ptype == ptype));
859:       if (flg) break;
860:       matmat = entry;
861:       entry = entry->next;
862:     }
863:     if (!flg) {
864:       PetscNew(&matmat->next);
865:       matmat = matmat->next;
866:     } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"This should not happen");
867:   }

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

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

885:    Logically Collective on Mat

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

896:    Level: advanced

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

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

920: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatSetOperation(), MatProductType, MatType, MatSetUp()
921: @*/
922: 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)
923: {

929:   if (ptype == MATPRODUCT_ABC) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Not for product type %s",MatProductTypes[ptype]);
930:   if (!numeric) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Missing numeric routine, argument 4");
933:   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));
934:   return(0);
935: }

937: 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)
938: {
939:   PetscBool      flg;
941:   char           composedname[256];
942:   MatRootName    Bnames = MatRootNameList, Cnames = MatRootNameList;
943:   PetscMPIInt    size;

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

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

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

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

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

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

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

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

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

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

1078:   if (shell->axpy) {
1079:     Mat              X;
1080:     PetscObjectState axpy_state;

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

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

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

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

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

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

1133:   if (shell->axpy) {
1134:     Mat              X;
1135:     PetscObjectState axpy_state;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1365: static PetscErrorCode MatMissingDiagonal_Shell(Mat A,PetscBool *missing,PetscInt *d)
1366: {
1368:   *missing = PETSC_FALSE;
1369:   return(0);
1370: }

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

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

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

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

1541:   shell->ctx = ctx;
1542:   return(0);
1543: }

1545: static PetscErrorCode MatShellSetVecType_Shell(Mat mat,VecType vtype)
1546: {

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

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

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

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

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

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

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

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

1682:   Level: advanced

1684: .seealso: MatCreateShell()
1685: M*/

1687: PETSC_EXTERN PetscErrorCode MatCreate_Shell(Mat A)
1688: {
1689:   Mat_Shell      *b;

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

1695:   PetscNewLog(A,&b);
1696:   A->data = (void*)b;

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

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

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

1720:   Collective

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

1730:    Output Parameter:
1731: .  A - the matrix

1733:    Level: advanced

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

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

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

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

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


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


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

1786:     Developers Notes:
1787:     Regarding shifting and scaling. The general form is

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

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

1795:           A is the user provided function.

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

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

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

1807: .seealso: MatShellSetOperation(), MatHasOperation(), MatShellGetContext(), MatShellSetContext(), MATSHELL, MatShellSetManageScalingShifts(), MatShellSetMatProductOperation()
1808: @*/
1809: PetscErrorCode  MatCreateShell(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,void *ctx,Mat *A)
1810: {

1814:   MatCreate(comm,A);
1815:   MatSetSizes(*A,m,n,M,N);
1816:   MatSetType(*A,MATSHELL);
1817:   MatShellSetContext(*A,ctx);
1818:   MatSetUp(*A);
1819:   return(0);
1820: }

1822: /*@
1823:     MatShellSetContext - sets the context for a shell matrix

1825:    Logically Collective on Mat

1827:     Input Parameters:
1828: +   mat - the shell matrix
1829: -   ctx - the context

1831:    Level: advanced

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

1837: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation()
1838: @*/
1839: PetscErrorCode  MatShellSetContext(Mat mat,void *ctx)
1840: {

1845:   PetscTryMethod(mat,"MatShellSetContext_C",(Mat,void*),(mat,ctx));
1846:   return(0);
1847: }

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

1852:  Logically collective

1854:     Input Parameters:
1855: +   mat   - the shell matrix
1856: -   vtype - type to use for creating vectors

1858:  Notes:

1860:  Level: advanced

1862: .seealso: MatCreateVecs()
1863: @*/
1864: PetscErrorCode  MatShellSetVecType(Mat mat,VecType vtype)
1865: {

1869:   PetscTryMethod(mat,"MatShellSetVecType_C",(Mat,VecType),(mat,vtype));
1870:   return(0);
1871: }

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

1877:    Logically Collective on Mat

1879:     Input Parameter:
1880: .   mat - the shell matrix

1882:   Level: advanced

1884: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatShellSetOperation()
1885: @*/
1886: PetscErrorCode MatShellSetManageScalingShifts(Mat A)
1887: {

1892:   PetscTryMethod(A,"MatShellSetManageScalingShifts_C",(Mat),(A));
1893:   return(0);
1894: }

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

1899:    Logically Collective on Mat

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

1907:    Output Parameter:
1908: .   flg - PETSC_TRUE if the multiply is likely correct

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

1913:    Level: advanced

1915:    Fortran Notes:
1916:     Not supported from Fortran

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

1930:   PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_view",&v);
1931:   MatGetLocalSize(mat,&m,&n);
1932:   MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,PETSC_DECIDE,PETSC_DECIDE,&mf);
1933:   MatMFFDSetFunction(mf,f,ctx);
1934:   MatMFFDSetBase(mf,base,NULL);

1936:   MatComputeOperator(mf,MATAIJ,&Dmf);
1937:   MatComputeOperator(mat,MATAIJ,&Dmat);

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

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

1965:    Logically Collective on Mat

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

1973:    Output Parameter:
1974: .   flg - PETSC_TRUE if the multiply is likely correct

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

1979:    Level: advanced

1981:    Fortran Notes:
1982:     Not supported from Fortran

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

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

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

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

2038:    Logically Collective on Mat

2040:     Input Parameters:
2041: +   mat - the shell matrix
2042: .   op - the name of the operation
2043: -   g - the function that provides the operation.

2045:    Level: advanced

2047:     Usage:
2048: $      extern PetscErrorCode usermult(Mat,Vec,Vec);
2049: $      MatCreateShell(comm,m,n,M,N,ctx,&A);
2050: $      MatShellSetOperation(A,MATOP_MULT,(void(*)(void))usermult);

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

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

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

2067:     Within each user-defined routine, the user should call
2068:     MatShellGetContext() to obtain the user-defined context that was
2069:     set by MatCreateShell().

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

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

2077: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatSetOperation(), MatShellSetManageScalingShifts(), MatShellSetMatProductOperation()
2078: @*/
2079: PetscErrorCode MatShellSetOperation(Mat mat,MatOperation op,void (*g)(void))
2080: {

2085:   PetscTryMethod(mat,"MatShellSetOperation_C",(Mat,MatOperation,void (*)(void)),(mat,op,g));
2086:   return(0);
2087: }

2089: /*@C
2090:     MatShellGetOperation - Gets a matrix function for a shell matrix.

2092:     Not Collective

2094:     Input Parameters:
2095: +   mat - the shell matrix
2096: -   op - the name of the operation

2098:     Output Parameter:
2099: .   g - the function that provides the operation.

2101:     Level: advanced

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

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

2115:     Within each user-defined routine, the user should call
2116:     MatShellGetContext() to obtain the user-defined context that was
2117:     set by MatCreateShell().

2119: .seealso: MatCreateShell(), MatShellGetContext(), MatShellSetOperation(), MatShellSetContext()
2120: @*/
2121: PetscErrorCode MatShellGetOperation(Mat mat,MatOperation op,void(**g)(void))
2122: {

2127:   PetscUseMethod(mat,"MatShellGetOperation_C",(Mat,MatOperation,void (**)(void)),(mat,op,g));
2128:   return(0);
2129: }

2131: /*@
2132:     MatIsShell - Inquires if a matrix is derived from MATSHELL

2134:     Input Parameter:
2135: .   mat - the matrix

2137:     Output Parameter:
2138: .   flg - the boolean value

2140:     Level: developer

2142:     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)

2144: .seealso: MatCreateShell()
2145: @*/
2146: PetscErrorCode MatIsShell(Mat mat, PetscBool *flg)
2147: {
2151:   *flg = (PetscBool)(mat->ops->destroy == MatDestroy_Shell);
2152:   return(0);
2153: }