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: 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: 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:   if (ctx->ops->destroy) (*ctx->ops->destroy)(ctx);

 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 differene 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:     Notes:
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: 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:    Notes:
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: 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) {
203:     VecDestroy(&ctx->current_f);
204:   }
205:   if (ctx->ops->destroy) (*ctx->ops->destroy)(ctx);
206:   PetscHeaderDestroy(&ctx);

208:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetBase_C",NULL);
209:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctioniBase_C",NULL);
210:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctioni_C",NULL);
211:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunction_C",NULL);
212:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctionError_C",NULL);
213:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetCheckh_C",NULL);
214:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetPeriod_C",NULL);
215:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDResetHHistory_C",NULL);
216:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetHHistory_C",NULL);
217:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetType_C",NULL);
218:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDGetH_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) {
245:       PetscViewerASCIIPrintf(viewer,"Using Lyness complex number trick to compute the matrix-vector product\n");
246:     }
247: #endif
248:     if (ctx->ops->view) {
249:       (*ctx->ops->view)(ctx,viewer);
250:     }
251:     PetscObjectGetOptionsPrefix((PetscObject)J, &prefix);

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

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

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

282:   MatShellGetContext(J,&j);
283:   MatMFFDResetHHistory(J);
284:   return 0;
285: }

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

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

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

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

327:   if (ctx->checkh) {
328:     (*ctx->checkh)(ctx->checkhctx,U,a,&h);
329:   }

331:   /* keep a record of the current differencing parameter h */
332:   ctx->currenth = h;
333: #if defined(PETSC_USE_COMPLEX)
334:   PetscInfo(mat,"Current differencing parameter: %g + %g i\n",(double)PetscRealPart(h),(double)PetscImaginaryPart(h));
335: #else
336:   PetscInfo(mat,"Current differencing parameter: %15.12e\n",(double)PetscRealPart(h));
337: #endif
338:   if (ctx->historyh && ctx->ncurrenth < ctx->maxcurrenth) {
339:     ctx->historyh[ctx->ncurrenth] = h;
340:   }
341:   ctx->ncurrenth++;

343: #if defined(PETSC_USE_COMPLEX)
344:   if (ctx->usecomplex) h = PETSC_i*h;
345: #endif

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

350:   /* compute func(U) as base for differencing; only needed first time in and not when provided by user */
351:   if (ctx->ncurrenth == 1 && ctx->current_f_allocated) {
352:     (*ctx->func)(ctx->funcctx,U,F);
353:   }
354:   (*ctx->func)(ctx->funcctx,w,y);

356: #if defined(PETSC_USE_COMPLEX)
357:   if (ctx->usecomplex) {
358:     VecImaginaryPart(y);
359:     h    = PetscImaginaryPart(h);
360:   } else {
361:     VecAXPY(y,-1.0,F);
362:   }
363: #else
364:   VecAXPY(y,-1.0,F);
365: #endif
366:   VecScale(y,1.0/h);
367:   if (mat->nullsp) MatNullSpaceRemove(mat->nullsp,y);

369:   PetscLogEventEnd(MATMFFD_Mult,a,y,0,0);
370:   return 0;
371: }

373: /*
374:   MatGetDiagonal_MFFD - Gets the diagonal for a matrix free matrix

376:         y ~= (F(u + ha) - F(u))/h,
377:   where F = nonlinear function, as set by SNESSetFunction()
378:         u = current iterate
379:         h = difference interval
380: */
381: PetscErrorCode MatGetDiagonal_MFFD(Mat mat,Vec a)
382: {
383:   MatMFFD        ctx;
384:   PetscScalar    h,*aa,*ww,v;
385:   PetscReal      epsilon = PETSC_SQRT_MACHINE_EPSILON,umin = 100.0*PETSC_SQRT_MACHINE_EPSILON;
386:   Vec            w,U;
387:   PetscInt       i,rstart,rend;

389:   MatShellGetContext(mat,&ctx);
392:   w    = ctx->w;
393:   U    = ctx->current_u;
394:   (*ctx->func)(ctx->funcctx,U,a);
395:   if (ctx->funcisetbase) {
396:     (*ctx->funcisetbase)(ctx->funcctx,U);
397:   }
398:   VecCopy(U,w);

400:   VecGetOwnershipRange(a,&rstart,&rend);
401:   VecGetArray(a,&aa);
402:   for (i=rstart; i<rend; i++) {
403:     VecGetArray(w,&ww);
404:     h    = ww[i-rstart];
405:     if (h == 0.0) h = 1.0;
406:     if (PetscAbsScalar(h) < umin && PetscRealPart(h) >= 0.0)     h = umin;
407:     else if (PetscRealPart(h) < 0.0 && PetscAbsScalar(h) < umin) h = -umin;
408:     h *= epsilon;

410:     ww[i-rstart] += h;
411:     VecRestoreArray(w,&ww);
412:     (*ctx->funci)(ctx->funcctx,i,w,&v);
413:     aa[i-rstart]  = (v - aa[i-rstart])/h;

415:     VecGetArray(w,&ww);
416:     ww[i-rstart] -= h;
417:     VecRestoreArray(w,&ww);
418:   }
419:   VecRestoreArray(a,&aa);
420:   return 0;
421: }

423: PETSC_EXTERN PetscErrorCode MatMFFDSetBase_MFFD(Mat J,Vec U,Vec F)
424: {
425:   MatMFFD        ctx;

427:   MatShellGetContext(J,&ctx);
428:   MatMFFDResetHHistory(J);
429:   if (!ctx->current_u) {
430:     VecDuplicate(U,&ctx->current_u);
431:     VecLockReadPush(ctx->current_u);
432:   }
433:   VecLockReadPop(ctx->current_u);
434:   VecCopy(U,ctx->current_u);
435:   VecLockReadPush(ctx->current_u);
436:   if (F) {
437:     if (ctx->current_f_allocated) VecDestroy(&ctx->current_f);
438:     ctx->current_f           = F;
439:     ctx->current_f_allocated = PETSC_FALSE;
440:   } else if (!ctx->current_f_allocated) {
441:     MatCreateVecs(J,NULL,&ctx->current_f);
442:     ctx->current_f_allocated = PETSC_TRUE;
443:   }
444:   if (!ctx->w) {
445:     VecDuplicate(ctx->current_u,&ctx->w);
446:   }
447:   J->assembled = PETSC_TRUE;
448:   return 0;
449: }

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

453: static PetscErrorCode  MatMFFDSetCheckh_MFFD(Mat J,FCN3 fun,void *ectx)
454: {
455:   MatMFFD        ctx;

457:   MatShellGetContext(J,&ctx);
458:   ctx->checkh    = fun;
459:   ctx->checkhctx = ectx;
460:   return 0;
461: }

463: /*@C
464:    MatMFFDSetOptionsPrefix - Sets the prefix used for searching for all
465:    MatMFFD options in the database.

467:    Collective on Mat

469:    Input Parameters:
470: +  A - the Mat context
471: -  prefix - the prefix to prepend to all option names

473:    Notes:
474:    A hyphen (-) must NOT be given at the beginning of the prefix name.
475:    The first character of all runtime options is AUTOMATICALLY the hyphen.

477:    Level: advanced

479: .seealso: MatSetFromOptions(), MatCreateSNESMF(), MatCreateMFFD()
480: @*/
481: PetscErrorCode  MatMFFDSetOptionsPrefix(Mat mat,const char prefix[])
482: {
483:   MatMFFD mfctx;

486:   MatShellGetContext(mat,&mfctx);
488:   PetscObjectSetOptionsPrefix((PetscObject)mfctx,prefix);
489:   return 0;
490: }

492: static PetscErrorCode  MatSetFromOptions_MFFD(PetscOptionItems *PetscOptionsObject,Mat mat)
493: {
494:   MatMFFD        mfctx;
496:   PetscBool      flg;
497:   char           ftype[256];

500:   MatShellGetContext(mat,&mfctx);
502:   PetscObjectOptionsBegin((PetscObject)mfctx);
503:   PetscOptionsFList("-mat_mffd_type","Matrix free type","MatMFFDSetType",MatMFFDList,((PetscObject)mfctx)->type_name,ftype,256,&flg);
504:   if (flg) {
505:     MatMFFDSetType(mat,ftype);
506:   }

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

511:   flg  = PETSC_FALSE;
512:   PetscOptionsBool("-mat_mffd_check_positivity","Insure that U + h*a is nonnegative","MatMFFDSetCheckh",flg,&flg,NULL);
513:   if (flg) {
514:     MatMFFDSetCheckh(mat,MatMFFDCheckPositivity,NULL);
515:   }
516: #if defined(PETSC_USE_COMPLEX)
517:   PetscOptionsBool("-mat_mffd_complex","Use Lyness complex number trick to compute the matrix-vector product","None",mfctx->usecomplex,&mfctx->usecomplex,NULL);
518: #endif
519:   if (mfctx->ops->setfromoptions) {
520:     (*mfctx->ops->setfromoptions)(PetscOptionsObject,mfctx);
521:   }
522:   PetscOptionsEnd();
523:   return 0;
524: }

526: static PetscErrorCode  MatMFFDSetPeriod_MFFD(Mat mat,PetscInt period)
527: {
528:   MatMFFD        ctx;

530:   MatShellGetContext(mat,&ctx);
531:   ctx->recomputeperiod = period;
532:   return 0;
533: }

535: static PetscErrorCode  MatMFFDSetFunction_MFFD(Mat mat,PetscErrorCode (*func)(void*,Vec,Vec),void *funcctx)
536: {
537:   MatMFFD        ctx;

539:   MatShellGetContext(mat,&ctx);
540:   ctx->func    = func;
541:   ctx->funcctx = funcctx;
542:   return 0;
543: }

545: static PetscErrorCode  MatMFFDSetFunctionError_MFFD(Mat mat,PetscReal error)
546: {
547:   MatMFFD        ctx;

549:   MatShellGetContext(mat,&ctx);
550:   if (error != PETSC_DEFAULT) ctx->error_rel = error;
551:   return 0;
552: }

554: PetscErrorCode  MatMFFDSetHHistory_MFFD(Mat J,PetscScalar history[],PetscInt nhistory)
555: {
556:   MatMFFD        ctx;

558:   MatShellGetContext(J,&ctx);
559:   ctx->historyh    = history;
560:   ctx->maxcurrenth = nhistory;
561:   ctx->currenth    = 0.;
562:   return 0;
563: }

565: /*MC
566:   MATMFFD - MATMFFD = "mffd" - A matrix free matrix type.

568:   Level: advanced

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

572: .seealso: MatCreateMFFD(), MatCreateSNESMF(), MatMFFDSetFunction(), MatMFFDSetType(),
573:           MatMFFDSetFunctionError(), MatMFFDDSSetUmin(), MatMFFDSetFunction()
574:           MatMFFDSetHHistory(), MatMFFDResetHHistory(), MatCreateSNESMF(),
575:           MatMFFDGetH(),
576: M*/
577: PETSC_EXTERN PetscErrorCode MatCreate_MFFD(Mat A)
578: {
579:   MatMFFD        mfctx;

581:   MatMFFDInitializePackage();

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

585:   mfctx->error_rel                = PETSC_SQRT_MACHINE_EPSILON;
586:   mfctx->recomputeperiod          = 1;
587:   mfctx->count                    = 0;
588:   mfctx->currenth                 = 0.0;
589:   mfctx->historyh                 = NULL;
590:   mfctx->ncurrenth                = 0;
591:   mfctx->maxcurrenth              = 0;
592:   ((PetscObject)mfctx)->type_name = NULL;

594:   /*
595:      Create the empty data structure to contain compute-h routines.
596:      These will be filled in below from the command line options or
597:      a later call with MatMFFDSetType() or if that is not called
598:      then it will default in the first use of MatMult_MFFD()
599:   */
600:   mfctx->ops->compute        = NULL;
601:   mfctx->ops->destroy        = NULL;
602:   mfctx->ops->view           = NULL;
603:   mfctx->ops->setfromoptions = NULL;
604:   mfctx->hctx                = NULL;

606:   mfctx->func    = NULL;
607:   mfctx->funcctx = NULL;
608:   mfctx->w       = NULL;
609:   mfctx->mat     = A;

611:   MatSetType(A,MATSHELL);
612:   MatShellSetContext(A,mfctx);
613:   MatShellSetOperation(A,MATOP_MULT,(void (*)(void))MatMult_MFFD);
614:   MatShellSetOperation(A,MATOP_DESTROY,(void (*)(void))MatDestroy_MFFD);
615:   MatShellSetOperation(A,MATOP_VIEW,(void (*)(void))MatView_MFFD);
616:   MatShellSetOperation(A,MATOP_ASSEMBLY_END,(void (*)(void))MatAssemblyEnd_MFFD);
617:   MatShellSetOperation(A,MATOP_SET_FROM_OPTIONS,(void (*)(void))MatSetFromOptions_MFFD);

619:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetBase_C",MatMFFDSetBase_MFFD);
620:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctioniBase_C",MatMFFDSetFunctioniBase_MFFD);
621:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctioni_C",MatMFFDSetFunctioni_MFFD);
622:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunction_C",MatMFFDSetFunction_MFFD);
623:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetCheckh_C",MatMFFDSetCheckh_MFFD);
624:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetPeriod_C",MatMFFDSetPeriod_MFFD);
625:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctionError_C",MatMFFDSetFunctionError_MFFD);
626:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDResetHHistory_C",MatMFFDResetHHistory_MFFD);
627:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetHHistory_C",MatMFFDSetHHistory_MFFD);
628:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetType_C",MatMFFDSetType_MFFD);
629:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDGetH_C",MatMFFDGetH_MFFD);
630:   PetscObjectChangeTypeName((PetscObject)A,MATMFFD);
631:   return 0;
632: }

634: /*@
635:    MatCreateMFFD - Creates a matrix-free matrix. See also MatCreateSNESMF()

637:    Collective on Vec

639:    Input Parameters:
640: +  comm - MPI communicator
641: .  m - number of local rows (or PETSC_DECIDE to have calculated if M is given)
642:            This value should be the same as the local size used in creating the
643:            y vector for the matrix-vector product y = Ax.
644: .  n - This value should be the same as the local size used in creating the
645:        x vector for the matrix-vector product y = Ax. (or PETSC_DECIDE to have
646:        calculated if N is given) For square matrices n is almost always m.
647: .  M - number of global rows (or PETSC_DETERMINE to have calculated if m is given)
648: -  N - number of global columns (or PETSC_DETERMINE to have calculated if n is given)

650:    Output Parameter:
651: .  J - the matrix-free matrix

653:    Options Database Keys: call MatSetFromOptions() to trigger these
654: +  -mat_mffd_type - wp or ds (see MATMFFD_WP or MATMFFD_DS)
655: .  -mat_mffd_err - square root of estimated relative error in function evaluation
656: .  -mat_mffd_period - how often h is recomputed, defaults to 1, every time
657: .  -mat_mffd_check_positivity - possibly decrease h until U + h*a has only positive values
658: -  -mat_mffd_complex - use the Lyness trick with complex numbers to compute the matrix-vector product instead of differencing
659:                        (requires real valued functions but that PETSc be configured for complex numbers)

661:    Level: advanced

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

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

670: .vb
671:      F'(u)*a = [F(u+h*a) - F(u)]/h where
672:      h = error_rel*u'a/||a||^2                        if  |u'a| > umin*||a||_{1}
673:        = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2   otherwise
674:  where
675:      error_rel = square root of relative error in function evaluation
676:      umin = minimum iterate parameter
677: .ve

679:    You can call SNESSetJacobian() with MatMFFDComputeJacobian() if you are using matrix and not a different
680:    preconditioner matrix

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

685:    The user should call MatDestroy() when finished with the matrix-free
686:    matrix context.

688:    Options Database Keys:
689: +  -mat_mffd_err <error_rel> - Sets error_rel
690: .  -mat_mffd_umin <umin> - Sets umin (for default PETSc routine that computes h only)
691: -  -mat_mffd_check_positivity - check positivity

693: .seealso: MatDestroy(), MatMFFDSetFunctionError(), MatMFFDDSSetUmin(), MatMFFDSetFunction()
694:           MatMFFDSetHHistory(), MatMFFDResetHHistory(), MatCreateSNESMF(),
695:           MatMFFDGetH(), MatMFFDRegister(), MatMFFDComputeJacobian()

697: @*/
698: PetscErrorCode  MatCreateMFFD(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,Mat *J)
699: {
700:   MatCreate(comm,J);
701:   MatSetSizes(*J,m,n,M,N);
702:   MatSetType(*J,MATMFFD);
703:   MatSetUp(*J);
704:   return 0;
705: }

707: /*@
708:    MatMFFDGetH - Gets the last value that was used as the differencing
709:    parameter.

711:    Not Collective

713:    Input Parameters:
714: .  mat - the matrix obtained with MatCreateSNESMF()

716:    Output Parameter:
717: .  h - the differencing step size

719:    Level: advanced

721: .seealso: MatCreateSNESMF(),MatMFFDSetHHistory(), MatCreateMFFD(), MATMFFD, MatMFFDResetHHistory()
722: @*/
723: PetscErrorCode  MatMFFDGetH(Mat mat,PetscScalar *h)
724: {
727:   PetscUseMethod(mat,"MatMFFDGetH_C",(Mat,PetscScalar*),(mat,h));
728:   return 0;
729: }

731: /*@C
732:    MatMFFDSetFunction - Sets the function used in applying the matrix free.

734:    Logically Collective on Mat

736:    Input Parameters:
737: +  mat - the matrix free matrix created via MatCreateSNESMF() or MatCreateMFFD()
738: .  func - the function to use
739: -  funcctx - optional function context passed to function

741:    Calling Sequence of func:
742: $     func (void *funcctx, Vec x, Vec f)

744: +  funcctx - user provided context
745: .  x - input vector
746: -  f - computed output function

748:    Level: advanced

750:    Notes:
751:     If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
752:     matrix inside your compute Jacobian routine

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

756: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD,
757:           MatMFFDSetHHistory(), MatMFFDResetHHistory(), SNESetFunction()
758: @*/
759: PetscErrorCode  MatMFFDSetFunction(Mat mat,PetscErrorCode (*func)(void*,Vec,Vec),void *funcctx)
760: {
762:   PetscTryMethod(mat,"MatMFFDSetFunction_C",(Mat,PetscErrorCode (*)(void*,Vec,Vec),void*),(mat,func,funcctx));
763:   return 0;
764: }

766: /*@C
767:    MatMFFDSetFunctioni - Sets the function for a single component

769:    Logically Collective on Mat

771:    Input Parameters:
772: +  mat - the matrix free matrix created via MatCreateSNESMF()
773: -  funci - the function to use

775:    Level: advanced

777:    Notes:
778:     If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
779:     matrix inside your compute Jacobian routine.
780:     This function is necessary to compute the diagonal of the matrix.
781:     funci must not contain any MPI call as it is called inside a loop on the local portion of the vector.

783: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatMFFDSetHHistory(), MatMFFDResetHHistory(), SNESetFunction(), MatGetDiagonal()

785: @*/
786: PetscErrorCode  MatMFFDSetFunctioni(Mat mat,PetscErrorCode (*funci)(void*,PetscInt,Vec,PetscScalar*))
787: {
789:   PetscTryMethod(mat,"MatMFFDSetFunctioni_C",(Mat,PetscErrorCode (*)(void*,PetscInt,Vec,PetscScalar*)),(mat,funci));
790:   return 0;
791: }

793: /*@C
794:    MatMFFDSetFunctioniBase - Sets the base vector for a single component function evaluation

796:    Logically Collective on Mat

798:    Input Parameters:
799: +  mat - the matrix free matrix created via MatCreateSNESMF()
800: -  func - the function to use

802:    Level: advanced

804:    Notes:
805:     If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
806:     matrix inside your compute Jacobian routine.
807:     This function is necessary to compute the diagonal of the matrix.

809: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD
810:           MatMFFDSetHHistory(), MatMFFDResetHHistory(), SNESetFunction(), MatGetDiagonal()
811: @*/
812: PetscErrorCode  MatMFFDSetFunctioniBase(Mat mat,PetscErrorCode (*func)(void*,Vec))
813: {
815:   PetscTryMethod(mat,"MatMFFDSetFunctioniBase_C",(Mat,PetscErrorCode (*)(void*,Vec)),(mat,func));
816:   return 0;
817: }

819: /*@
820:    MatMFFDSetPeriod - Sets how often h is recomputed, by default it is every time

822:    Logically Collective on Mat

824:    Input Parameters:
825: +  mat - the matrix free matrix created via MatCreateSNESMF()
826: -  period - 1 for every time, 2 for every second etc

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

831:    Level: advanced

833: .seealso: MatCreateSNESMF(),MatMFFDGetH(),
834:           MatMFFDSetHHistory(), MatMFFDResetHHistory()
835: @*/
836: PetscErrorCode  MatMFFDSetPeriod(Mat mat,PetscInt period)
837: {
840:   PetscTryMethod(mat,"MatMFFDSetPeriod_C",(Mat,PetscInt),(mat,period));
841:   return 0;
842: }

844: /*@
845:    MatMFFDSetFunctionError - Sets the error_rel for the approximation of
846:    matrix-vector products using finite differences.

848:    Logically Collective on Mat

850:    Input Parameters:
851: +  mat - the matrix free matrix created via MatCreateMFFD() or MatCreateSNESMF()
852: -  error_rel - relative error (should be set to the square root of
853:                the relative error in the function evaluations)

855:    Options Database Keys:
856: .  -mat_mffd_err <error_rel> - Sets error_rel

858:    Level: advanced

860:    Notes:
861:    The default matrix-free matrix-vector product routine computes
862: .vb
863:      F'(u)*a = [F(u+h*a) - F(u)]/h where
864:      h = error_rel*u'a/||a||^2                        if  |u'a| > umin*||a||_{1}
865:        = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2   else
866: .ve

868: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD
869:           MatMFFDSetHHistory(), MatMFFDResetHHistory()
870: @*/
871: PetscErrorCode  MatMFFDSetFunctionError(Mat mat,PetscReal error)
872: {
875:   PetscTryMethod(mat,"MatMFFDSetFunctionError_C",(Mat,PetscReal),(mat,error));
876:   return 0;
877: }

879: /*@
880:    MatMFFDSetHHistory - Sets an array to collect a history of the
881:    differencing values (h) computed for the matrix-free product.

883:    Logically Collective on Mat

885:    Input Parameters:
886: +  J - the matrix-free matrix context
887: .  history - space to hold the history
888: -  nhistory - number of entries in history, if more entries are generated than
889:               nhistory, then the later ones are discarded

891:    Level: advanced

893:    Notes:
894:    Use MatMFFDResetHHistory() to reset the history counter and collect
895:    a new batch of differencing parameters, h.

897: .seealso: MatMFFDGetH(), MatCreateSNESMF(),
898:           MatMFFDResetHHistory(), MatMFFDSetFunctionError()

900: @*/
901: PetscErrorCode  MatMFFDSetHHistory(Mat J,PetscScalar history[],PetscInt nhistory)
902: {
906:   PetscUseMethod(J,"MatMFFDSetHHistory_C",(Mat,PetscScalar[],PetscInt),(J,history,nhistory));
907:   return 0;
908: }

910: /*@
911:    MatMFFDResetHHistory - Resets the counter to zero to begin
912:    collecting a new set of differencing histories.

914:    Logically Collective on Mat

916:    Input Parameters:
917: .  J - the matrix-free matrix context

919:    Level: advanced

921:    Notes:
922:    Use MatMFFDSetHHistory() to create the original history counter.

924: .seealso: MatMFFDGetH(), MatCreateSNESMF(),
925:           MatMFFDSetHHistory(), MatMFFDSetFunctionError()

927: @*/
928: PetscErrorCode  MatMFFDResetHHistory(Mat J)
929: {
931:   PetscTryMethod(J,"MatMFFDResetHHistory_C",(Mat),(J));
932:   return 0;
933: }

935: /*@
936:     MatMFFDSetBase - Sets the vector U at which matrix vector products of the
937:         Jacobian are computed

939:     Logically Collective on Mat

941:     Input Parameters:
942: +   J - the MatMFFD matrix
943: .   U - the vector
944: -   F - (optional) vector that contains F(u) if it has been already computed

946:     Notes:
947:     This is rarely used directly

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

952:     Level: advanced

954: @*/
955: PetscErrorCode  MatMFFDSetBase(Mat J,Vec U,Vec F)
956: {
960:   PetscTryMethod(J,"MatMFFDSetBase_C",(Mat,Vec,Vec),(J,U,F));
961:   return 0;
962: }

964: /*@C
965:     MatMFFDSetCheckh - Sets a function that checks the computed h and adjusts
966:         it to satisfy some criteria

968:     Logically Collective on Mat

970:     Input Parameters:
971: +   J - the MatMFFD matrix
972: .   fun - the function that checks h
973: -   ctx - any context needed by the function

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

978:     Level: advanced

980:     Notes:
981:     For example, MatMFFDCheckPositivity() insures that all entries
982:        of U + h*a are non-negative

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

987: .seealso:  MatMFFDCheckPositivity()
988: @*/
989: PetscErrorCode  MatMFFDSetCheckh(Mat J,PetscErrorCode (*fun)(void*,Vec,Vec,PetscScalar*),void *ctx)
990: {
992:   PetscTryMethod(J,"MatMFFDSetCheckh_C",(Mat,PetscErrorCode (*)(void*,Vec,Vec,PetscScalar*),void*),(J,fun,ctx));
993:   return 0;
994: }

996: /*@
997:     MatMFFDCheckPositivity - Checks that all entries in U + h*a are positive or
998:         zero, decreases h until this is satisfied.

1000:     Logically Collective on Vec

1002:     Input Parameters:
1003: +   U - base vector that is added to
1004: .   a - vector that is added
1005: .   h - scaling factor on a
1006: -   dummy - context variable (unused)

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

1011:     Level: advanced

1013:     Notes:
1014:     This is rarely used directly, rather it is passed as an argument to
1015:            MatMFFDSetCheckh()

1017: .seealso:  MatMFFDSetCheckh()
1018: @*/
1019: PetscErrorCode  MatMFFDCheckPositivity(void *dummy,Vec U,Vec a,PetscScalar *h)
1020: {
1021:   PetscReal      val, minval;
1022:   PetscScalar    *u_vec, *a_vec;
1023:   PetscInt       i,n;
1024:   MPI_Comm       comm;

1029:   PetscObjectGetComm((PetscObject)U,&comm);
1030:   VecGetArray(U,&u_vec);
1031:   VecGetArray(a,&a_vec);
1032:   VecGetLocalSize(U,&n);
1033:   minval = PetscAbsScalar(*h)*PetscRealConstant(1.01);
1034:   for (i=0; i<n; i++) {
1035:     if (PetscRealPart(u_vec[i] + *h*a_vec[i]) <= 0.0) {
1036:       val = PetscAbsScalar(u_vec[i]/a_vec[i]);
1037:       if (val < minval) minval = val;
1038:     }
1039:   }
1040:   VecRestoreArray(U,&u_vec);
1041:   VecRestoreArray(a,&a_vec);
1042:   MPIU_Allreduce(&minval,&val,1,MPIU_REAL,MPIU_MIN,comm);
1043:   if (val <= PetscAbsScalar(*h)) {
1044:     PetscInfo(U,"Scaling back h from %g to %g\n",(double)PetscRealPart(*h),(double)(.99*val));
1045:     if (PetscRealPart(*h) > 0.0) *h =  0.99*val;
1046:     else                         *h = -0.99*val;
1047:   }
1048:   return 0;
1049: }