Actual source code: mffd.c


  2: #include <petsc/private/matimpl.h>
  3: #include <../src/mat/impls/mffd/mffdimpl.h>

  5: PetscFunctionList MatMFFDList              = NULL;
  6: PetscBool         MatMFFDRegisterAllCalled = PETSC_FALSE;

  8: PetscClassId  MATMFFD_CLASSID;
  9: PetscLogEvent MATMFFD_Mult;

 11: static PetscBool MatMFFDPackageInitialized = PETSC_FALSE;
 12: /*@C
 13:   MatMFFDFinalizePackage - This function destroys everything in the MATMFFD` package. It is
 14:   called from `PetscFinalize()`.

 16:   Level: developer

 18: .seealso: `MATMFFD`, `PetscFinalize()`, `MatCreateMFFD()`, `MatCreateSNESMF()`
 19: @*/
 20: PetscErrorCode MatMFFDFinalizePackage(void)
 21: {
 22:   PetscFunctionListDestroy(&MatMFFDList);
 23:   MatMFFDPackageInitialized = PETSC_FALSE;
 24:   MatMFFDRegisterAllCalled  = PETSC_FALSE;
 25:   return 0;
 26: }

 28: /*@C
 29:   MatMFFDInitializePackage - This function initializes everything in the MATMFFD` package. It is called
 30:   from `MatInitializePackage()`.

 32:   Level: developer

 34: .seealso: `MATMFFD`, `PetscInitialize()`
 35: @*/
 36: PetscErrorCode MatMFFDInitializePackage(void)
 37: {
 38:   char      logList[256];
 39:   PetscBool opt, pkg;

 41:   if (MatMFFDPackageInitialized) return 0;
 42:   MatMFFDPackageInitialized = PETSC_TRUE;
 43:   /* Register Classes */
 44:   PetscClassIdRegister("MatMFFD", &MATMFFD_CLASSID);
 45:   /* Register Constructors */
 46:   MatMFFDRegisterAll();
 47:   /* Register Events */
 48:   PetscLogEventRegister("MatMult MF", MATMFFD_CLASSID, &MATMFFD_Mult);
 49:   /* Process Info */
 50:   {
 51:     PetscClassId classids[1];

 53:     classids[0] = MATMFFD_CLASSID;
 54:     PetscInfoProcessClass("matmffd", 1, classids);
 55:   }
 56:   /* Process summary exclusions */
 57:   PetscOptionsGetString(NULL, NULL, "-log_exclude", logList, sizeof(logList), &opt);
 58:   if (opt) {
 59:     PetscStrInList("matmffd", logList, ',', &pkg);
 60:     if (pkg) PetscLogEventExcludeClass(MATMFFD_CLASSID);
 61:   }
 62:   /* Register package finalizer */
 63:   PetscRegisterFinalize(MatMFFDFinalizePackage);
 64:   return 0;
 65: }

 67: static PetscErrorCode MatMFFDSetType_MFFD(Mat mat, MatMFFDType ftype)
 68: {
 69:   MatMFFD   ctx;
 70:   PetscBool match;
 71:   PetscErrorCode (*r)(MatMFFD);

 75:   MatShellGetContext(mat, &ctx);

 77:   /* already set, so just return */
 78:   PetscObjectTypeCompare((PetscObject)ctx, ftype, &match);
 79:   if (match) return 0;

 81:   /* destroy the old one if it exists */
 82:   PetscTryTypeMethod(ctx, destroy);

 84:   PetscFunctionListFind(MatMFFDList, ftype, &r);
 86:   (*r)(ctx);
 87:   PetscObjectChangeTypeName((PetscObject)ctx, ftype);
 88:   return 0;
 89: }

 91: /*@C
 92:     MatMFFDSetType - Sets the method that is used to compute the
 93:     differencing parameter for finite difference matrix-free formulations.

 95:     Input Parameters:
 96: +   mat - the "matrix-free" matrix created via `MatCreateSNESMF()`, or `MatCreateMFFD()`
 97:           or `MatSetType`(mat,`MATMFFD`);
 98: -   ftype - the type requested, either `MATMFFD_WP` or `MATMFFD_DS`

100:     Level: advanced

102:     Note:
103:     For example, such routines can compute h for use in
104:     Jacobian-vector products of the form

106:                         F(x+ha) - F(x)
107:           F'(u)a  ~=  ----------------
108:                               h

110: .seealso: `MATMFFD`, `MATMFFD_WP`, `MATMFFD_DS`, `MatCreateSNESMF()`, `MatMFFDRegister()`, `MatMFFDSetFunction()`, `MatCreateMFFD()`
111: @*/
112: PetscErrorCode MatMFFDSetType(Mat mat, MatMFFDType ftype)
113: {
116:   PetscTryMethod(mat, "MatMFFDSetType_C", (Mat, MatMFFDType), (mat, ftype));
117:   return 0;
118: }

120: static PetscErrorCode MatGetDiagonal_MFFD(Mat, Vec);

122: typedef PetscErrorCode (*FCN1)(void *, Vec); /* force argument to next function to not be extern C*/
123: static PetscErrorCode MatMFFDSetFunctioniBase_MFFD(Mat mat, FCN1 func)
124: {
125:   MatMFFD ctx;

127:   MatShellGetContext(mat, &ctx);
128:   ctx->funcisetbase = func;
129:   return 0;
130: }

132: typedef PetscErrorCode (*FCN2)(void *, PetscInt, Vec, PetscScalar *); /* force argument to next function to not be extern C*/
133: static PetscErrorCode MatMFFDSetFunctioni_MFFD(Mat mat, FCN2 funci)
134: {
135:   MatMFFD ctx;

137:   MatShellGetContext(mat, &ctx);
138:   ctx->funci = funci;
139:   MatShellSetOperation(mat, MATOP_GET_DIAGONAL, (void (*)(void))MatGetDiagonal_MFFD);
140:   return 0;
141: }

143: static PetscErrorCode MatMFFDGetH_MFFD(Mat mat, PetscScalar *h)
144: {
145:   MatMFFD ctx;

147:   MatShellGetContext(mat, &ctx);
148:   *h = ctx->currenth;
149:   return 0;
150: }

152: static PetscErrorCode MatMFFDResetHHistory_MFFD(Mat J)
153: {
154:   MatMFFD ctx;

156:   MatShellGetContext(J, &ctx);
157:   ctx->ncurrenth = 0;
158:   return 0;
159: }

161: /*@C
162:    MatMFFDRegister - Adds a method to the MATMFFD` registry.

164:    Not Collective

166:    Input Parameters:
167: +  name_solver - name of a new user-defined compute-h module
168: -  routine_create - routine to create method context

170:    Level: developer

172:    Note:
173:    `MatMFFDRegister()` may be called multiple times to add several user-defined solvers.

175:    Sample usage:
176: .vb
177:    MatMFFDRegister("my_h",MyHCreate);
178: .ve

180:    Then, your solver can be chosen with the procedural interface via
181: $     `MatMFFDSetType`(mfctx,"my_h")
182:    or at runtime via the option
183: $     -mat_mffd_type my_h

185: .seealso: `MATMFFD`, `MatMFFDRegisterAll()`, `MatMFFDRegisterDestroy()`
186:  @*/
187: PetscErrorCode MatMFFDRegister(const char sname[], PetscErrorCode (*function)(MatMFFD))
188: {
189:   MatInitializePackage();
190:   PetscFunctionListAdd(&MatMFFDList, sname, function);
191:   return 0;
192: }

194: /* ----------------------------------------------------------------------------------------*/
195: static PetscErrorCode MatDestroy_MFFD(Mat mat)
196: {
197:   MatMFFD ctx;

199:   MatShellGetContext(mat, &ctx);
200:   VecDestroy(&ctx->w);
201:   VecDestroy(&ctx->current_u);
202:   if (ctx->current_f_allocated) VecDestroy(&ctx->current_f);
203:   PetscTryTypeMethod(ctx, destroy);
204:   PetscHeaderDestroy(&ctx);

206:   PetscObjectComposeFunction((PetscObject)mat, "MatMFFDSetBase_C", NULL);
207:   PetscObjectComposeFunction((PetscObject)mat, "MatMFFDSetFunctioniBase_C", NULL);
208:   PetscObjectComposeFunction((PetscObject)mat, "MatMFFDSetFunctioni_C", NULL);
209:   PetscObjectComposeFunction((PetscObject)mat, "MatMFFDSetFunction_C", NULL);
210:   PetscObjectComposeFunction((PetscObject)mat, "MatMFFDSetFunctionError_C", NULL);
211:   PetscObjectComposeFunction((PetscObject)mat, "MatMFFDSetCheckh_C", NULL);
212:   PetscObjectComposeFunction((PetscObject)mat, "MatMFFDSetPeriod_C", NULL);
213:   PetscObjectComposeFunction((PetscObject)mat, "MatMFFDResetHHistory_C", NULL);
214:   PetscObjectComposeFunction((PetscObject)mat, "MatMFFDSetHHistory_C", NULL);
215:   PetscObjectComposeFunction((PetscObject)mat, "MatMFFDSetType_C", NULL);
216:   PetscObjectComposeFunction((PetscObject)mat, "MatMFFDGetH_C", NULL);
217:   PetscObjectComposeFunction((PetscObject)mat, "MatSNESMFSetReuseBase_C", NULL);
218:   PetscObjectComposeFunction((PetscObject)mat, "MatSNESMFGetReuseBase_C", NULL);
219:   return 0;
220: }

222: /*
223:    MatMFFDView_MFFD - Views matrix-free parameters.

225: */
226: static PetscErrorCode MatView_MFFD(Mat J, PetscViewer viewer)
227: {
228:   MatMFFD     ctx;
229:   PetscBool   iascii, viewbase, viewfunction;
230:   const char *prefix;

232:   MatShellGetContext(J, &ctx);
233:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
234:   if (iascii) {
235:     PetscViewerASCIIPrintf(viewer, "Matrix-free approximation:\n");
236:     PetscViewerASCIIPushTab(viewer);
237:     PetscViewerASCIIPrintf(viewer, "err=%g (relative error in function evaluation)\n", (double)ctx->error_rel);
238:     if (!((PetscObject)ctx)->type_name) {
239:       PetscViewerASCIIPrintf(viewer, "The compute h routine has not yet been set\n");
240:     } else {
241:       PetscViewerASCIIPrintf(viewer, "Using %s compute h routine\n", ((PetscObject)ctx)->type_name);
242:     }
243: #if defined(PETSC_USE_COMPLEX)
244:     if (ctx->usecomplex) PetscViewerASCIIPrintf(viewer, "Using Lyness complex number trick to compute the matrix-vector product\n");
245: #endif
246:     PetscTryTypeMethod(ctx, view, viewer);
247:     PetscObjectGetOptionsPrefix((PetscObject)J, &prefix);

249:     PetscOptionsHasName(((PetscObject)J)->options, prefix, "-mat_mffd_view_base", &viewbase);
250:     if (viewbase) {
251:       PetscViewerASCIIPrintf(viewer, "Base:\n");
252:       VecView(ctx->current_u, viewer);
253:     }
254:     PetscOptionsHasName(((PetscObject)J)->options, prefix, "-mat_mffd_view_function", &viewfunction);
255:     if (viewfunction) {
256:       PetscViewerASCIIPrintf(viewer, "Function:\n");
257:       VecView(ctx->current_f, viewer);
258:     }
259:     PetscViewerASCIIPopTab(viewer);
260:   }
261:   return 0;
262: }

264: /*
265:    MatAssemblyEnd_MFFD - Resets the ctx->ncurrenth to zero. This
266:    allows the user to indicate the beginning of a new linear solve by calling
267:    MatAssemblyXXX() on the matrix free matrix. This then allows the
268:    MatCreateMFFD_WP() to properly compute ||U|| only the first time
269:    in the linear solver rather than every time.

271:    This function is referenced directly from MatAssemblyEnd_SNESMF(), which may be in a different shared library hence
272:    it must be labeled as PETSC_EXTERN
273: */
274: PETSC_EXTERN PetscErrorCode MatAssemblyEnd_MFFD(Mat J, MatAssemblyType mt)
275: {
276:   MatMFFD j;

278:   MatShellGetContext(J, &j);
279:   MatMFFDResetHHistory(J);
280:   return 0;
281: }

283: /*
284:   MatMult_MFFD - Default matrix-free form for Jacobian-vector product, y = F'(u)*a:

286:         y ~= (F(u + ha) - F(u))/h,
287:   where F = nonlinear function, as set by SNESSetFunction()
288:         u = current iterate
289:         h = difference interval
290: */
291: static PetscErrorCode MatMult_MFFD(Mat mat, Vec a, Vec y)
292: {
293:   MatMFFD     ctx;
294:   PetscScalar h;
295:   Vec         w, U, F;
296:   PetscBool   zeroa;

298:   MatShellGetContext(mat, &ctx);
300:   /* We log matrix-free matrix-vector products separately, so that we can
301:      separate the performance monitoring from the cases that use conventional
302:      storage.  We may eventually modify event logging to associate events
303:      with particular objects, hence alleviating the more general problem. */
304:   PetscLogEventBegin(MATMFFD_Mult, a, y, 0, 0);

306:   w = ctx->w;
307:   U = ctx->current_u;
308:   F = ctx->current_f;
309:   /*
310:       Compute differencing parameter
311:   */
312:   if (!((PetscObject)ctx)->type_name) {
313:     MatMFFDSetType(mat, MATMFFD_WP);
314:     MatSetFromOptions(mat);
315:   }
316:   PetscUseTypeMethod(ctx, compute, U, a, &h, &zeroa);
317:   if (zeroa) {
318:     VecSet(y, 0.0);
319:     PetscLogEventEnd(MATMFFD_Mult, a, y, 0, 0);
320:     return 0;
321:   }

324:   if (ctx->checkh) (*ctx->checkh)(ctx->checkhctx, U, a, &h);

326:   /* keep a record of the current differencing parameter h */
327:   ctx->currenth = h;
328: #if defined(PETSC_USE_COMPLEX)
329:   PetscInfo(mat, "Current differencing parameter: %g + %g i\n", (double)PetscRealPart(h), (double)PetscImaginaryPart(h));
330: #else
331:   PetscInfo(mat, "Current differencing parameter: %15.12e\n", (double)PetscRealPart(h));
332: #endif
333:   if (ctx->historyh && ctx->ncurrenth < ctx->maxcurrenth) ctx->historyh[ctx->ncurrenth] = h;
334:   ctx->ncurrenth++;

336: #if defined(PETSC_USE_COMPLEX)
337:   if (ctx->usecomplex) h = PETSC_i * h;
338: #endif

340:   /* w = u + ha */
341:   VecWAXPY(w, h, a, U);

343:   /* compute func(U) as base for differencing; only needed first time in and not when provided by user */
344:   if (ctx->ncurrenth == 1 && ctx->current_f_allocated) (*ctx->func)(ctx->funcctx, U, F);
345:   (*ctx->func)(ctx->funcctx, w, y);

347: #if defined(PETSC_USE_COMPLEX)
348:   if (ctx->usecomplex) {
349:     VecImaginaryPart(y);
350:     h = PetscImaginaryPart(h);
351:   } else {
352:     VecAXPY(y, -1.0, F);
353:   }
354: #else
355:   VecAXPY(y, -1.0, F);
356: #endif
357:   VecScale(y, 1.0 / h);
358:   if (mat->nullsp) MatNullSpaceRemove(mat->nullsp, y);

360:   PetscLogEventEnd(MATMFFD_Mult, a, y, 0, 0);
361:   return 0;
362: }

364: /*
365:   MatGetDiagonal_MFFD - Gets the diagonal for a matrix free matrix

367:         y ~= (F(u + ha) - F(u))/h,
368:   where F = nonlinear function, as set by SNESSetFunction()
369:         u = current iterate
370:         h = difference interval
371: */
372: PetscErrorCode MatGetDiagonal_MFFD(Mat mat, Vec a)
373: {
374:   MatMFFD     ctx;
375:   PetscScalar h, *aa, *ww, v;
376:   PetscReal   epsilon = PETSC_SQRT_MACHINE_EPSILON, umin = 100.0 * PETSC_SQRT_MACHINE_EPSILON;
377:   Vec         w, U;
378:   PetscInt    i, rstart, rend;

380:   MatShellGetContext(mat, &ctx);
383:   w = ctx->w;
384:   U = ctx->current_u;
385:   (*ctx->func)(ctx->funcctx, U, a);
386:   if (ctx->funcisetbase) (*ctx->funcisetbase)(ctx->funcctx, U);
387:   VecCopy(U, w);

389:   VecGetOwnershipRange(a, &rstart, &rend);
390:   VecGetArray(a, &aa);
391:   for (i = rstart; i < rend; i++) {
392:     VecGetArray(w, &ww);
393:     h = ww[i - rstart];
394:     if (h == 0.0) h = 1.0;
395:     if (PetscAbsScalar(h) < umin && PetscRealPart(h) >= 0.0) h = umin;
396:     else if (PetscRealPart(h) < 0.0 && PetscAbsScalar(h) < umin) h = -umin;
397:     h *= epsilon;

399:     ww[i - rstart] += h;
400:     VecRestoreArray(w, &ww);
401:     (*ctx->funci)(ctx->funcctx, i, w, &v);
402:     aa[i - rstart] = (v - aa[i - rstart]) / h;

404:     VecGetArray(w, &ww);
405:     ww[i - rstart] -= h;
406:     VecRestoreArray(w, &ww);
407:   }
408:   VecRestoreArray(a, &aa);
409:   return 0;
410: }

412: PETSC_EXTERN PetscErrorCode MatMFFDSetBase_MFFD(Mat J, Vec U, Vec F)
413: {
414:   MatMFFD ctx;

416:   MatShellGetContext(J, &ctx);
417:   MatMFFDResetHHistory(J);
418:   if (!ctx->current_u) {
419:     VecDuplicate(U, &ctx->current_u);
420:     VecLockReadPush(ctx->current_u);
421:   }
422:   VecLockReadPop(ctx->current_u);
423:   VecCopy(U, ctx->current_u);
424:   VecLockReadPush(ctx->current_u);
425:   if (F) {
426:     if (ctx->current_f_allocated) VecDestroy(&ctx->current_f);
427:     ctx->current_f           = F;
428:     ctx->current_f_allocated = PETSC_FALSE;
429:   } else if (!ctx->current_f_allocated) {
430:     MatCreateVecs(J, NULL, &ctx->current_f);
431:     ctx->current_f_allocated = PETSC_TRUE;
432:   }
433:   if (!ctx->w) VecDuplicate(ctx->current_u, &ctx->w);
434:   J->assembled = PETSC_TRUE;
435:   return 0;
436: }

438: typedef PetscErrorCode (*FCN3)(void *, Vec, Vec, PetscScalar *); /* force argument to next function to not be extern C*/

440: static PetscErrorCode MatMFFDSetCheckh_MFFD(Mat J, FCN3 fun, void *ectx)
441: {
442:   MatMFFD ctx;

444:   MatShellGetContext(J, &ctx);
445:   ctx->checkh    = fun;
446:   ctx->checkhctx = ectx;
447:   return 0;
448: }

450: /*@C
451:    MatMFFDSetOptionsPrefix - Sets the prefix used for searching for all
452:    MATMFFD` options in the database.

454:    Collective

456:    Input Parameters:
457: +  A - the `MATMFFD` context
458: -  prefix - the prefix to prepend to all option names

460:    Note:
461:    A hyphen (-) must NOT be given at the beginning of the prefix name.
462:    The first character of all runtime options is AUTOMATICALLY the hyphen.

464:    Level: advanced

466: .seealso: `MATMFFD`, `MatSetFromOptions()`, `MatCreateSNESMF()`, `MatCreateMFFD()`
467: @*/
468: PetscErrorCode MatMFFDSetOptionsPrefix(Mat mat, const char prefix[])
469: {
470:   MatMFFD mfctx;

473:   MatShellGetContext(mat, &mfctx);
475:   PetscObjectSetOptionsPrefix((PetscObject)mfctx, prefix);
476:   return 0;
477: }

479: static PetscErrorCode MatSetFromOptions_MFFD(Mat mat, PetscOptionItems *PetscOptionsObject)
480: {
481:   MatMFFD   mfctx;
482:   PetscBool flg;
483:   char      ftype[256];

485:   MatShellGetContext(mat, &mfctx);
487:   PetscObjectOptionsBegin((PetscObject)mfctx);
488:   PetscOptionsFList("-mat_mffd_type", "Matrix free type", "MatMFFDSetType", MatMFFDList, ((PetscObject)mfctx)->type_name, ftype, 256, &flg);
489:   if (flg) MatMFFDSetType(mat, ftype);

491:   PetscOptionsReal("-mat_mffd_err", "set sqrt relative error in function", "MatMFFDSetFunctionError", mfctx->error_rel, &mfctx->error_rel, NULL);
492:   PetscOptionsInt("-mat_mffd_period", "how often h is recomputed", "MatMFFDSetPeriod", mfctx->recomputeperiod, &mfctx->recomputeperiod, NULL);

494:   flg = PETSC_FALSE;
495:   PetscOptionsBool("-mat_mffd_check_positivity", "Insure that U + h*a is nonnegative", "MatMFFDSetCheckh", flg, &flg, NULL);
496:   if (flg) MatMFFDSetCheckh(mat, MatMFFDCheckPositivity, NULL);
497: #if defined(PETSC_USE_COMPLEX)
498:   PetscOptionsBool("-mat_mffd_complex", "Use Lyness complex number trick to compute the matrix-vector product", "None", mfctx->usecomplex, &mfctx->usecomplex, NULL);
499: #endif
500:   PetscTryTypeMethod(mfctx, setfromoptions, PetscOptionsObject);
501:   PetscOptionsEnd();
502:   return 0;
503: }

505: static PetscErrorCode MatMFFDSetPeriod_MFFD(Mat mat, PetscInt period)
506: {
507:   MatMFFD ctx;

509:   MatShellGetContext(mat, &ctx);
510:   ctx->recomputeperiod = period;
511:   return 0;
512: }

514: static PetscErrorCode MatMFFDSetFunction_MFFD(Mat mat, PetscErrorCode (*func)(void *, Vec, Vec), void *funcctx)
515: {
516:   MatMFFD ctx;

518:   MatShellGetContext(mat, &ctx);
519:   ctx->func    = func;
520:   ctx->funcctx = funcctx;
521:   return 0;
522: }

524: static PetscErrorCode MatMFFDSetFunctionError_MFFD(Mat mat, PetscReal error)
525: {
526:   MatMFFD ctx;

528:   MatShellGetContext(mat, &ctx);
529:   if (error != PETSC_DEFAULT) ctx->error_rel = error;
530:   return 0;
531: }

533: PetscErrorCode MatMFFDSetHHistory_MFFD(Mat J, PetscScalar history[], PetscInt nhistory)
534: {
535:   MatMFFD ctx;

537:   MatShellGetContext(J, &ctx);
538:   ctx->historyh    = history;
539:   ctx->maxcurrenth = nhistory;
540:   ctx->currenth    = 0.;
541:   return 0;
542: }

544: /*MC
545:   MATMFFD - MATMFFD = "mffd" - A matrix free matrix type.

547:   Level: advanced

549:   Developers Note: This is implemented on top of `MATSHELL` to get support for scaling and shifting without requiring duplicate code

551: .seealso: `MatCreateMFFD()`, `MatCreateSNESMF()`, `MatMFFDSetFunction()`, `MatMFFDSetType()`,
552:           `MatMFFDSetFunctionError()`, `MatMFFDDSSetUmin()`, `MatMFFDSetFunction()`
553:           `MatMFFDSetHHistory()`, `MatMFFDResetHHistory()`, `MatCreateSNESMF()`,
554:           `MatMFFDGetH()`,
555: M*/
556: PETSC_EXTERN PetscErrorCode MatCreate_MFFD(Mat A)
557: {
558:   MatMFFD mfctx;

560:   MatMFFDInitializePackage();

562:   PetscHeaderCreate(mfctx, MATMFFD_CLASSID, "MatMFFD", "Matrix-free Finite Differencing", "Mat", PetscObjectComm((PetscObject)A), NULL, NULL);

564:   mfctx->error_rel                = PETSC_SQRT_MACHINE_EPSILON;
565:   mfctx->recomputeperiod          = 1;
566:   mfctx->count                    = 0;
567:   mfctx->currenth                 = 0.0;
568:   mfctx->historyh                 = NULL;
569:   mfctx->ncurrenth                = 0;
570:   mfctx->maxcurrenth              = 0;
571:   ((PetscObject)mfctx)->type_name = NULL;

573:   /*
574:      Create the empty data structure to contain compute-h routines.
575:      These will be filled in below from the command line options or
576:      a later call with MatMFFDSetType() or if that is not called
577:      then it will default in the first use of MatMult_MFFD()
578:   */
579:   mfctx->ops->compute        = NULL;
580:   mfctx->ops->destroy        = NULL;
581:   mfctx->ops->view           = NULL;
582:   mfctx->ops->setfromoptions = NULL;
583:   mfctx->hctx                = NULL;

585:   mfctx->func    = NULL;
586:   mfctx->funcctx = NULL;
587:   mfctx->w       = NULL;
588:   mfctx->mat     = A;

590:   MatSetType(A, MATSHELL);
591:   MatShellSetContext(A, mfctx);
592:   MatShellSetOperation(A, MATOP_MULT, (void (*)(void))MatMult_MFFD);
593:   MatShellSetOperation(A, MATOP_DESTROY, (void (*)(void))MatDestroy_MFFD);
594:   MatShellSetOperation(A, MATOP_VIEW, (void (*)(void))MatView_MFFD);
595:   MatShellSetOperation(A, MATOP_ASSEMBLY_END, (void (*)(void))MatAssemblyEnd_MFFD);
596:   MatShellSetOperation(A, MATOP_SET_FROM_OPTIONS, (void (*)(void))MatSetFromOptions_MFFD);

598:   PetscObjectComposeFunction((PetscObject)A, "MatMFFDSetBase_C", MatMFFDSetBase_MFFD);
599:   PetscObjectComposeFunction((PetscObject)A, "MatMFFDSetFunctioniBase_C", MatMFFDSetFunctioniBase_MFFD);
600:   PetscObjectComposeFunction((PetscObject)A, "MatMFFDSetFunctioni_C", MatMFFDSetFunctioni_MFFD);
601:   PetscObjectComposeFunction((PetscObject)A, "MatMFFDSetFunction_C", MatMFFDSetFunction_MFFD);
602:   PetscObjectComposeFunction((PetscObject)A, "MatMFFDSetCheckh_C", MatMFFDSetCheckh_MFFD);
603:   PetscObjectComposeFunction((PetscObject)A, "MatMFFDSetPeriod_C", MatMFFDSetPeriod_MFFD);
604:   PetscObjectComposeFunction((PetscObject)A, "MatMFFDSetFunctionError_C", MatMFFDSetFunctionError_MFFD);
605:   PetscObjectComposeFunction((PetscObject)A, "MatMFFDResetHHistory_C", MatMFFDResetHHistory_MFFD);
606:   PetscObjectComposeFunction((PetscObject)A, "MatMFFDSetHHistory_C", MatMFFDSetHHistory_MFFD);
607:   PetscObjectComposeFunction((PetscObject)A, "MatMFFDSetType_C", MatMFFDSetType_MFFD);
608:   PetscObjectComposeFunction((PetscObject)A, "MatMFFDGetH_C", MatMFFDGetH_MFFD);
609:   PetscObjectChangeTypeName((PetscObject)A, MATMFFD);
610:   return 0;
611: }

613: /*@
614:    MatCreateMFFD - Creates a matrix-free matrix of type `MATMFFD`. See also `MatCreateSNESMF()`

616:    Collective

618:    Input Parameters:
619: +  comm - MPI communicator
620: .  m - number of local rows (or `PETSC_DECIDE` to have calculated if M is given)
621:            This value should be the same as the local size used in creating the
622:            y vector for the matrix-vector product y = Ax.
623: .  n - This value should be the same as the local size used in creating the
624:        x vector for the matrix-vector product y = Ax. (or `PETSC_DECIDE` to have
625:        calculated if N is given) For square matrices n is almost always m.
626: .  M - number of global rows (or `PETSC_DETERMINE` to have calculated if m is given)
627: -  N - number of global columns (or `PETSC_DETERMINE` to have calculated if n is given)

629:    Output Parameter:
630: .  J - the matrix-free matrix

632:    Options Database Keys:
633: +  -mat_mffd_type - wp or ds (see `MATMFFD_WP` or `MATMFFD_DS`)
634: .  -mat_mffd_err - square root of estimated relative error in function evaluation
635: .  -mat_mffd_period - how often h is recomputed, defaults to 1, every time
636: .  -mat_mffd_check_positivity - possibly decrease h until U + h*a has only positive values
637: .  -mat_mffd_umin <umin> - Sets umin (for default PETSc routine that computes h only)
638: -  -mat_mffd_complex - use the Lyness trick with complex numbers to compute the matrix-vector product instead of differencing
639:                        (requires real valued functions but that PETSc be configured for complex numbers)

641:    Level: advanced

643:    Notes:
644:    The matrix-free matrix context merely contains the function pointers
645:    and work space for performing finite difference approximations of
646:    Jacobian-vector products, F'(u)*a,

648:    The default code uses the following approach to compute h

650: .vb
651:      F'(u)*a = [F(u+h*a) - F(u)]/h where
652:      h = error_rel*u'a/||a||^2                        if  |u'a| > umin*||a||_{1}
653:        = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2   otherwise
654:  where
655:      error_rel = square root of relative error in function evaluation
656:      umin = minimum iterate parameter
657: .ve

659:    You can call `SNESSetJacobian()` with `MatMFFDComputeJacobian()` if you are using matrix and not a different
660:    preconditioner matrix

662:    The user can set the error_rel via `MatMFFDSetFunctionError()` and
663:    umin via `MatMFFDDSSetUmin()`; see Users-Manual: ch_snes for details.

665:    The user should call `MatDestroy()` when finished with the matrix-free
666:    matrix context.

668: .seealso: `MATMFFD`, `MatDestroy()`, `MatMFFDSetFunctionError()`, `MatMFFDDSSetUmin()`, `MatMFFDSetFunction()`
669:           `MatMFFDSetHHistory()`, `MatMFFDResetHHistory()`, `MatCreateSNESMF()`,
670:           `MatMFFDGetH()`, `MatMFFDRegister()`, `MatMFFDComputeJacobian()`
671: @*/
672: PetscErrorCode MatCreateMFFD(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, PetscInt N, Mat *J)
673: {
674:   MatCreate(comm, J);
675:   MatSetSizes(*J, m, n, M, N);
676:   MatSetType(*J, MATMFFD);
677:   MatSetUp(*J);
678:   return 0;
679: }

681: /*@
682:    MatMFFDGetH - Gets the last value that was used as the differencing for a `MATMFFD` matrix
683:    parameter.

685:    Not Collective

687:    Input Parameters:
688: .  mat - the `MATMFFD` matrix

690:    Output Parameter:
691: .  h - the differencing step size

693:    Level: advanced

695: .seealso: `MATMFFD`, `MatCreateSNESMF()`, `MatMFFDSetHHistory()`, `MatCreateMFFD()`, `MATMFFD`, `MatMFFDResetHHistory()`
696: @*/
697: PetscErrorCode MatMFFDGetH(Mat mat, PetscScalar *h)
698: {
701:   PetscUseMethod(mat, "MatMFFDGetH_C", (Mat, PetscScalar *), (mat, h));
702:   return 0;
703: }

705: /*@C
706:    MatMFFDSetFunction - Sets the function used in applying the matrix free `MATMFFD` matrix.

708:    Logically Collective

710:    Input Parameters:
711: +  mat - the matrix free matrix `MATMFFD` created via `MatCreateSNESMF()` or `MatCreateMFFD()`
712: .  func - the function to use
713: -  funcctx - optional function context passed to function

715:    Calling Sequence of func:
716: $     func (void *funcctx, Vec x, Vec f)

718: +  funcctx - user provided context
719: .  x - input vector
720: -  f - computed output function

722:    Level: advanced

724:    Notes:
725:     If you use this you MUST call `MatAssemblyBegin()` and `MatAssemblyEnd()` on the matrix free
726:     matrix inside your compute Jacobian routine

728:     If this is not set then it will use the function set with `SNESSetFunction()` if `MatCreateSNESMF()` was used.

730: .seealso: `MATMFFD`, `MatCreateSNESMF()`, `MatMFFDGetH()`, `MatCreateMFFD()`, `MATMFFD`,
731:           `MatMFFDSetHHistory()`, `MatMFFDResetHHistory()`, `SNESetFunction()`
732: @*/
733: PetscErrorCode MatMFFDSetFunction(Mat mat, PetscErrorCode (*func)(void *, Vec, Vec), void *funcctx)
734: {
736:   PetscTryMethod(mat, "MatMFFDSetFunction_C", (Mat, PetscErrorCode(*)(void *, Vec, Vec), void *), (mat, func, funcctx));
737:   return 0;
738: }

740: /*@C
741:    MatMFFDSetFunctioni - Sets the function for a single component for a `MATMFFD` matrix

743:    Logically Collective

745:    Input Parameters:
746: +  mat - the matrix free matrix `MATMFFD`
747: -  funci - the function to use

749:    Level: advanced

751:    Notes:
752:     If you use this you MUST call `MatAssemblyBegin()` and `MatAssemblyEnd()` on the matrix free
753:     matrix inside your compute Jacobian routine.

755:     This function is necessary to compute the diagonal of the matrix.
756:     funci must not contain any MPI call as it is called inside a loop on the local portion of the vector.

758: .seealso: `MATMFFD`, `MatCreateSNESMF()`, `MatMFFDGetH()`, `MatMFFDSetHHistory()`, `MatMFFDResetHHistory()`, `SNESetFunction()`, `MatGetDiagonal()`
759: @*/
760: PetscErrorCode MatMFFDSetFunctioni(Mat mat, PetscErrorCode (*funci)(void *, PetscInt, Vec, PetscScalar *))
761: {
763:   PetscTryMethod(mat, "MatMFFDSetFunctioni_C", (Mat, PetscErrorCode(*)(void *, PetscInt, Vec, PetscScalar *)), (mat, funci));
764:   return 0;
765: }

767: /*@C
768:    MatMFFDSetFunctioniBase - Sets the base vector for a single component function evaluation for a `MATMFFD` matrix

770:    Logically Collective

772:    Input Parameters:
773: +  mat - the `MATMFFD` matrix free matrix
774: -  func - the function to use

776:    Level: advanced

778:    Notes:
779:     If you use this you MUST call `MatAssemblyBegin()` and `MatAssemblyEnd()` on the matrix free
780:     matrix inside your compute Jacobian routine.

782:     This function is necessary to compute the diagonal of the matrix, used for example with `PCJACOBI`

784: .seealso: `MATMFFD`, `MatCreateSNESMF()`, `MatMFFDGetH()`, `MatCreateMFFD()`, `MATMFFD`
785:           `MatMFFDSetHHistory()`, `MatMFFDResetHHistory()`, `SNESetFunction()`, `MatGetDiagonal()`
786: @*/
787: PetscErrorCode MatMFFDSetFunctioniBase(Mat mat, PetscErrorCode (*func)(void *, Vec))
788: {
790:   PetscTryMethod(mat, "MatMFFDSetFunctioniBase_C", (Mat, PetscErrorCode(*)(void *, Vec)), (mat, func));
791:   return 0;
792: }

794: /*@
795:    MatMFFDSetPeriod - Sets how often h is recomputed for a `MATMFFD` matrix, by default it is every time

797:    Logically Collective

799:    Input Parameters:
800: +  mat - the `MATMFFD` matrix free matrix
801: -  period - 1 for every time, 2 for every second etc

803:    Options Database Keys:
804: .  -mat_mffd_period <period> - Sets how often h is recomputed

806:    Level: advanced

808: .seealso: `MATMFFD`, `MatCreateSNESMF()`, `MatMFFDGetH()`,
809:           `MatMFFDSetHHistory()`, `MatMFFDResetHHistory()`
810: @*/
811: PetscErrorCode MatMFFDSetPeriod(Mat mat, PetscInt period)
812: {
815:   PetscTryMethod(mat, "MatMFFDSetPeriod_C", (Mat, PetscInt), (mat, period));
816:   return 0;
817: }

819: /*@
820:    MatMFFDSetFunctionError - Sets the error_rel for the approximation of matrix-vector products using finite differences with the `MATMFFD` matrix

822:    Logically Collective

824:    Input Parameters:
825: +  mat - the `MATMFFD` matrix free matrix
826: -  error_rel - relative error (should be set to the square root of the relative error in the function evaluations)

828:    Options Database Keys:
829: .  -mat_mffd_err <error_rel> - Sets error_rel

831:    Level: advanced

833:    Note:
834:    The default matrix-free matrix-vector product routine computes
835: .vb
836:      F'(u)*a = [F(u+h*a) - F(u)]/h where
837:      h = error_rel*u'a/||a||^2                        if  |u'a| > umin*||a||_{1}
838:        = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2   else
839: .ve

841: .seealso: `MATMFFD`, `MatCreateSNESMF()`, `MatMFFDGetH()`, `MatCreateMFFD()`, `MATMFFD`
842:           `MatMFFDSetHHistory()`, `MatMFFDResetHHistory()`
843: @*/
844: PetscErrorCode MatMFFDSetFunctionError(Mat mat, PetscReal error)
845: {
848:   PetscTryMethod(mat, "MatMFFDSetFunctionError_C", (Mat, PetscReal), (mat, error));
849:   return 0;
850: }

852: /*@
853:    MatMFFDSetHHistory - Sets an array to collect a history of the
854:    differencing values (h) computed for the matrix-free product `MATMFFD` matrix

856:    Logically Collective

858:    Input Parameters:
859: +  J - the `MATMFFD` matrix-free matrix
860: .  history - space to hold the history
861: -  nhistory - number of entries in history, if more entries are generated than
862:               nhistory, then the later ones are discarded

864:    Level: advanced

866:    Note:
867:    Use `MatMFFDResetHHistory()` to reset the history counter and collect
868:    a new batch of differencing parameters, h.

870: .seealso: `MATMFFD`, `MatMFFDGetH()`, `MatCreateSNESMF()`,
871:           `MatMFFDResetHHistory()`, `MatMFFDSetFunctionError()`
872: @*/
873: PetscErrorCode MatMFFDSetHHistory(Mat J, PetscScalar history[], PetscInt nhistory)
874: {
878:   PetscUseMethod(J, "MatMFFDSetHHistory_C", (Mat, PetscScalar[], PetscInt), (J, history, nhistory));
879:   return 0;
880: }

882: /*@
883:    MatMFFDResetHHistory - Resets the counter to zero to begin
884:    collecting a new set of differencing histories for the `MATMFFD` matrix

886:    Logically Collective

888:    Input Parameters:
889: .  J - the matrix-free matrix context

891:    Level: advanced

893:    Note:
894:    Use `MatMFFDSetHHistory()` to create the original history counter.

896: .seealso: `MATMFFD`, `MatMFFDGetH()`, `MatCreateSNESMF()`,
897:           `MatMFFDSetHHistory()`, `MatMFFDSetFunctionError()`
898: @*/
899: PetscErrorCode MatMFFDResetHHistory(Mat J)
900: {
902:   PetscTryMethod(J, "MatMFFDResetHHistory_C", (Mat), (J));
903:   return 0;
904: }

906: /*@
907:     MatMFFDSetBase - Sets the vector U at which matrix vector products of the
908:         Jacobian are computed for the `MATMFFD` matrix

910:     Logically Collective

912:     Input Parameters:
913: +   J - the `MATMFFD` matrix
914: .   U - the vector
915: -   F - (optional) vector that contains F(u) if it has been already computed

917:     Notes:
918:     This is rarely used directly

920:     If F is provided then it is not recomputed. Otherwise the function is evaluated at the base
921:     point during the first `MatMult()` after each call to `MatMFFDSetBase()`.

923:     Level: advanced

925: .seealso: `MATMFFD`, `MatMult()`, `MatMFFDSetBase()`
926: @*/
927: PetscErrorCode MatMFFDSetBase(Mat J, Vec U, Vec F)
928: {
932:   PetscTryMethod(J, "MatMFFDSetBase_C", (Mat, Vec, Vec), (J, U, F));
933:   return 0;
934: }

936: /*@C
937:     MatMFFDSetCheckh - Sets a function that checks the computed h and adjusts
938:         it to satisfy some criteria for the `MATMFFD` matrix

940:     Logically Collective

942:     Input Parameters:
943: +   J - the `MATMFFD` matrix
944: .   fun - the function that checks h
945: -   ctx - any context needed by the function

947:     Options Database Keys:
948: .   -mat_mffd_check_positivity <bool> - Insure that U + h*a is non-negative

950:     Level: advanced

952:     Notes:
953:     For example, `MatMFFDCheckPositivity()` insures that all entries
954:        of U + h*a are non-negative

956:      The function you provide is called after the default h has been computed and allows you to
957:      modify it.

959: .seealso: `MATMFFD`, `MatMFFDCheckPositivity()`
960: @*/
961: PetscErrorCode MatMFFDSetCheckh(Mat J, PetscErrorCode (*fun)(void *, Vec, Vec, PetscScalar *), void *ctx)
962: {
964:   PetscTryMethod(J, "MatMFFDSetCheckh_C", (Mat, PetscErrorCode(*)(void *, Vec, Vec, PetscScalar *), void *), (J, fun, ctx));
965:   return 0;
966: }

968: /*@
969:     MatMFFDCheckPositivity - Checks that all entries in U + h*a are positive or
970:         zero, decreases h until this is satisfied for a `MATMFFD` matrix

972:     Logically Collective

974:     Input Parameters:
975: +   U - base vector that is added to
976: .   a - vector that is added
977: .   h - scaling factor on a
978: -   dummy - context variable (unused)

980:     Options Database Keys:
981: .   -mat_mffd_check_positivity <bool> - Insure that U + h*a is nonnegative

983:     Level: advanced

985:     Note:
986:     This is rarely used directly, rather it is passed as an argument to
987:            `MatMFFDSetCheckh()`

989: .seealso: `MATMFFD`, `MatMFFDSetCheckh()`
990: @*/
991: PetscErrorCode MatMFFDCheckPositivity(void *dummy, Vec U, Vec a, PetscScalar *h)
992: {
993:   PetscReal    val, minval;
994:   PetscScalar *u_vec, *a_vec;
995:   PetscInt     i, n;
996:   MPI_Comm     comm;

1001:   PetscObjectGetComm((PetscObject)U, &comm);
1002:   VecGetArray(U, &u_vec);
1003:   VecGetArray(a, &a_vec);
1004:   VecGetLocalSize(U, &n);
1005:   minval = PetscAbsScalar(*h) * PetscRealConstant(1.01);
1006:   for (i = 0; i < n; i++) {
1007:     if (PetscRealPart(u_vec[i] + *h * a_vec[i]) <= 0.0) {
1008:       val = PetscAbsScalar(u_vec[i] / a_vec[i]);
1009:       if (val < minval) minval = val;
1010:     }
1011:   }
1012:   VecRestoreArray(U, &u_vec);
1013:   VecRestoreArray(a, &a_vec);
1014:   MPIU_Allreduce(&minval, &val, 1, MPIU_REAL, MPIU_MIN, comm);
1015:   if (val <= PetscAbsScalar(*h)) {
1016:     PetscInfo(U, "Scaling back h from %g to %g\n", (double)PetscRealPart(*h), (double)(.99 * val));
1017:     if (PetscRealPart(*h) > 0.0) *h = 0.99 * val;
1018:     else *h = -0.99 * val;
1019:   }
1020:   return 0;
1021: }