Actual source code: mffd.c

petsc-master 2019-09-18
Report Typos and Errors

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

  5: PetscFunctionList MatMFFDList              = 0;
  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: {

 25:   PetscFunctionListDestroy(&MatMFFDList);
 26:   MatMFFDPackageInitialized = PETSC_FALSE;
 27:   MatMFFDRegisterAllCalled  = PETSC_FALSE;
 28:   return(0);
 29: }

 31: /*@C
 32:   MatMFFDInitializePackage - This function initializes everything in the MatMFFD package. It is called
 33:   from MatInitializePackage().

 35:   Level: developer

 37: .seealso: PetscInitialize()
 38: @*/
 39: PetscErrorCode  MatMFFDInitializePackage(void)
 40: {
 41:   char           logList[256];
 42:   PetscBool      opt,pkg;

 46:   if (MatMFFDPackageInitialized) return(0);
 47:   MatMFFDPackageInitialized = PETSC_TRUE;
 48:   /* Register Classes */
 49:   PetscClassIdRegister("MatMFFD",&MATMFFD_CLASSID);
 50:   /* Register Constructors */
 51:   MatMFFDRegisterAll();
 52:   /* Register Events */
 53:   PetscLogEventRegister("MatMult MF",MATMFFD_CLASSID,&MATMFFD_Mult);
 54:   /* Process info exclusions */
 55:   PetscOptionsGetString(NULL,NULL,"-info_exclude",logList,sizeof(logList),&opt);
 56:   if (opt) {
 57:     PetscStrInList("matmffd",logList,',',&pkg);
 58:     if (pkg) {PetscInfoDeactivateClass(MATMFFD_CLASSID);}
 59:   }
 60:   /* Process summary exclusions */
 61:   PetscOptionsGetString(NULL,NULL,"-log_exclude",logList,sizeof(logList),&opt);
 62:   if (opt) {
 63:     PetscStrInList("matmffd",logList,',',&pkg);
 64:     if (pkg) {PetscLogEventExcludeClass(MATMFFD_CLASSID);}
 65:   }
 66:   /* Register package finalizer */
 67:   PetscRegisterFinalize(MatMFFDFinalizePackage);
 68:   return(0);
 69: }

 71: static PetscErrorCode  MatMFFDSetType_MFFD(Mat mat,MatMFFDType ftype)
 72: {
 73:   PetscErrorCode ierr,(*r)(MatMFFD);
 74:   MatMFFD        ctx;
 75:   PetscBool      match;

 80:   MatShellGetContext(mat,&ctx);

 82:   /* already set, so just return */
 83:   PetscObjectTypeCompare((PetscObject)ctx,ftype,&match);
 84:   if (match) return(0);

 86:   /* destroy the old one if it exists */
 87:   if (ctx->ops->destroy) {
 88:     (*ctx->ops->destroy)(ctx);
 89:   }

 91:    PetscFunctionListFind(MatMFFDList,ftype,&r);
 92:   if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown MatMFFD type %s given",ftype);
 93:   (*r)(ctx);
 94:   PetscObjectChangeTypeName((PetscObject)ctx,ftype);
 95:   return(0);
 96: }

 98: /*@C
 99:     MatMFFDSetType - Sets the method that is used to compute the
100:     differencing parameter for finite differene matrix-free formulations.

102:     Input Parameters:
103: +   mat - the "matrix-free" matrix created via MatCreateSNESMF(), or MatCreateMFFD()
104:           or MatSetType(mat,MATMFFD);
105: -   ftype - the type requested, either MATMFFD_WP or MATMFFD_DS

107:     Level: advanced

109:     Notes:
110:     For example, such routines can compute h for use in
111:     Jacobian-vector products of the form

113:                         F(x+ha) - F(x)
114:           F'(u)a  ~=  ----------------
115:                               h

117: .seealso: MatCreateSNESMF(), MatMFFDRegister(), MatMFFDSetFunction(), MatCreateMFFD()
118: @*/
119: PetscErrorCode  MatMFFDSetType(Mat mat,MatMFFDType ftype)
120: {

126:   PetscTryMethod(mat,"MatMFFDSetType_C",(Mat,MatMFFDType),(mat,ftype));
127:   return(0);
128: }

130: static PetscErrorCode MatGetDiagonal_MFFD(Mat,Vec);

132: typedef PetscErrorCode (*FCN1)(void*,Vec); /* force argument to next function to not be extern C*/
133: static PetscErrorCode  MatMFFDSetFunctioniBase_MFFD(Mat mat,FCN1 func)
134: {
135:   MatMFFD        ctx;

139:   MatShellGetContext(mat,&ctx);
140:   ctx->funcisetbase = func;
141:   return(0);
142: }

144: typedef PetscErrorCode (*FCN2)(void*,PetscInt,Vec,PetscScalar*); /* force argument to next function to not be extern C*/
145: static PetscErrorCode  MatMFFDSetFunctioni_MFFD(Mat mat,FCN2 funci)
146: {
147:   MatMFFD        ctx;

151:   MatShellGetContext(mat,&ctx);
152:   ctx->funci = funci;
153:   MatShellSetOperation(mat,MATOP_GET_DIAGONAL,(void (*)(void))MatGetDiagonal_MFFD);
154:   return(0);
155: }

157: static PetscErrorCode MatMFFDGetH_MFFD(Mat mat,PetscScalar *h)
158: {
159:   MatMFFD        ctx;

163:   MatShellGetContext(mat,&ctx);
164:   *h = ctx->currenth;
165:   return(0);
166: }

168: static PetscErrorCode  MatMFFDResetHHistory_MFFD(Mat J)
169: {
170:   MatMFFD        ctx;

174:   MatShellGetContext(J,&ctx);
175:   ctx->ncurrenth = 0;
176:   return(0);
177: }

179: /*@C
180:    MatMFFDRegister - Adds a method to the MatMFFD registry.

182:    Not Collective

184:    Input Parameters:
185: +  name_solver - name of a new user-defined compute-h module
186: -  routine_create - routine to create method context

188:    Level: developer

190:    Notes:
191:    MatMFFDRegister() may be called multiple times to add several user-defined solvers.

193:    Sample usage:
194: .vb
195:    MatMFFDRegister("my_h",MyHCreate);
196: .ve

198:    Then, your solver can be chosen with the procedural interface via
199: $     MatMFFDSetType(mfctx,"my_h")
200:    or at runtime via the option
201: $     -mat_mffd_type my_h

203: .seealso: MatMFFDRegisterAll(), MatMFFDRegisterDestroy()
204:  @*/
205: PetscErrorCode  MatMFFDRegister(const char sname[],PetscErrorCode (*function)(MatMFFD))
206: {

210:   MatInitializePackage();
211:   PetscFunctionListAdd(&MatMFFDList,sname,function);
212:   return(0);
213: }

215: /* ----------------------------------------------------------------------------------------*/
216: static PetscErrorCode MatDestroy_MFFD(Mat mat)
217: {
219:   MatMFFD        ctx;

222:   MatShellGetContext(mat,&ctx);
223:   VecDestroy(&ctx->w);
224:   VecDestroy(&ctx->current_u);
225:   if (ctx->current_f_allocated) {
226:     VecDestroy(&ctx->current_f);
227:   }
228:   if (ctx->ops->destroy) {(*ctx->ops->destroy)(ctx);}
229:   PetscHeaderDestroy(&ctx);

231:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetBase_C",NULL);
232:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctioniBase_C",NULL);
233:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctioni_C",NULL);
234:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunction_C",NULL);
235:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctionError_C",NULL);
236:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetCheckh_C",NULL);
237:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetPeriod_C",NULL);
238:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDResetHHistory_C",NULL);
239:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetHHistory_C",NULL);
240:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetType_C",NULL);
241:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDGetH_C",NULL);
242:   return(0);
243: }

245: /*
246:    MatMFFDView_MFFD - Views matrix-free parameters.

248: */
249: static PetscErrorCode MatView_MFFD(Mat J,PetscViewer viewer)
250: {
252:   MatMFFD        ctx;
253:   PetscBool      iascii, viewbase, viewfunction;
254:   const char     *prefix;

257:   MatShellGetContext(J,&ctx);
258:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
259:   if (iascii) {
260:     PetscViewerASCIIPrintf(viewer,"Matrix-free approximation:\n");
261:     PetscViewerASCIIPushTab(viewer);
262:     PetscViewerASCIIPrintf(viewer,"err=%g (relative error in function evaluation)\n",(double)ctx->error_rel);
263:     if (!((PetscObject)ctx)->type_name) {
264:       PetscViewerASCIIPrintf(viewer,"The compute h routine has not yet been set\n");
265:     } else {
266:       PetscViewerASCIIPrintf(viewer,"Using %s compute h routine\n",((PetscObject)ctx)->type_name);
267:     }
268: #if defined(PETSC_USE_COMPLEX)
269:     if (ctx->usecomplex) {
270:       PetscViewerASCIIPrintf(viewer,"Using Lyness complex number trick to compute the matrix-vector product\n");
271:     }
272: #endif
273:     if (ctx->ops->view) {
274:       (*ctx->ops->view)(ctx,viewer);
275:     }
276:     PetscObjectGetOptionsPrefix((PetscObject)J, &prefix);

278:     PetscOptionsHasName(((PetscObject)J)->options,prefix, "-mat_mffd_view_base", &viewbase);
279:     if (viewbase) {
280:       PetscViewerASCIIPrintf(viewer, "Base:\n");
281:       VecView(ctx->current_u, viewer);
282:     }
283:     PetscOptionsHasName(((PetscObject)J)->options,prefix, "-mat_mffd_view_function", &viewfunction);
284:     if (viewfunction) {
285:       PetscViewerASCIIPrintf(viewer, "Function:\n");
286:       VecView(ctx->current_f, viewer);
287:     }
288:     PetscViewerASCIIPopTab(viewer);
289:   }
290:   return(0);
291: }

293: /*
294:    MatAssemblyEnd_MFFD - Resets the ctx->ncurrenth to zero. This
295:    allows the user to indicate the beginning of a new linear solve by calling
296:    MatAssemblyXXX() on the matrix free matrix. This then allows the
297:    MatCreateMFFD_WP() to properly compute ||U|| only the first time
298:    in the linear solver rather than every time.

300:    This function is referenced directly from MatAssemblyEnd_SNESMF(), which may be in a different shared library hence
301:    it must be labeled as PETSC_EXTERN
302: */
303: PETSC_EXTERN PetscErrorCode MatAssemblyEnd_MFFD(Mat J,MatAssemblyType mt)
304: {
306:   MatMFFD        j;

309:   MatShellGetContext(J,&j);
310:   MatMFFDResetHHistory(J);
311:   return(0);
312: }

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

317:         y ~= (F(u + ha) - F(u))/h,
318:   where F = nonlinear function, as set by SNESSetFunction()
319:         u = current iterate
320:         h = difference interval
321: */
322: static PetscErrorCode MatMult_MFFD(Mat mat,Vec a,Vec y)
323: {
324:   MatMFFD        ctx;
325:   PetscScalar    h;
326:   Vec            w,U,F;
328:   PetscBool      zeroa;

331:   MatShellGetContext(mat,&ctx);
332:   if (!ctx->current_u) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"MatMFFDSetBase() has not been called, this is often caused by forgetting to call \n\t\tMatAssemblyBegin/End on the first Mat in the SNES compute function");
333:   /* We log matrix-free matrix-vector products separately, so that we can
334:      separate the performance monitoring from the cases that use conventional
335:      storage.  We may eventually modify event logging to associate events
336:      with particular objects, hence alleviating the more general problem. */
337:   PetscLogEventBegin(MATMFFD_Mult,a,y,0,0);

339:   w = ctx->w;
340:   U = ctx->current_u;
341:   F = ctx->current_f;
342:   /*
343:       Compute differencing parameter
344:   */
345:   if (!((PetscObject)ctx)->type_name) {
346:     MatMFFDSetType(mat,MATMFFD_WP);
347:     MatSetFromOptions(mat);
348:   }
349:   (*ctx->ops->compute)(ctx,U,a,&h,&zeroa);
350:   if (zeroa) {
351:     VecSet(y,0.0);
352:     return(0);
353:   }

355:   if (mat->erroriffailure && PetscIsInfOrNanScalar(h)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Computed Nan differencing parameter h");
356:   if (ctx->checkh) {
357:     (*ctx->checkh)(ctx->checkhctx,U,a,&h);
358:   }

360:   /* keep a record of the current differencing parameter h */
361:   ctx->currenth = h;
362: #if defined(PETSC_USE_COMPLEX)
363:   PetscInfo2(mat,"Current differencing parameter: %g + %g i\n",(double)PetscRealPart(h),(double)PetscImaginaryPart(h));
364: #else
365:   PetscInfo1(mat,"Current differencing parameter: %15.12e\n",h);
366: #endif
367:   if (ctx->historyh && ctx->ncurrenth < ctx->maxcurrenth) {
368:     ctx->historyh[ctx->ncurrenth] = h;
369:   }
370:   ctx->ncurrenth++;

372: #if defined(PETSC_USE_COMPLEX)
373:   if (ctx->usecomplex) h = PETSC_i*h;
374: #endif

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

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

385: #if defined(PETSC_USE_COMPLEX)  
386:   if (ctx->usecomplex) {
387:     VecImaginaryPart(y);
388:     h    = PetscImaginaryPart(h);
389:   } else {
390:     VecAXPY(y,-1.0,F);
391:   }
392: #else
393:   VecAXPY(y,-1.0,F);
394: #endif
395:   VecScale(y,1.0/h);
396:   if (mat->nullsp) {MatNullSpaceRemove(mat->nullsp,y);}

398:   PetscLogEventEnd(MATMFFD_Mult,a,y,0,0);
399:   return(0);
400: }

402: /*
403:   MatGetDiagonal_MFFD - Gets the diagonal for a matrix free matrix

405:         y ~= (F(u + ha) - F(u))/h,
406:   where F = nonlinear function, as set by SNESSetFunction()
407:         u = current iterate
408:         h = difference interval
409: */
410: PetscErrorCode MatGetDiagonal_MFFD(Mat mat,Vec a)
411: {
412:   MatMFFD        ctx;
413:   PetscScalar    h,*aa,*ww,v;
414:   PetscReal      epsilon = PETSC_SQRT_MACHINE_EPSILON,umin = 100.0*PETSC_SQRT_MACHINE_EPSILON;
415:   Vec            w,U;
417:   PetscInt       i,rstart,rend;

420:   MatShellGetContext(mat,&ctx);
421:   if (!ctx->func) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Requires calling MatMFFDSetFunction() first");
422:   if (!ctx->funci) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Requires calling MatMFFDSetFunctioni() first");
423:   w    = ctx->w;
424:   U    = ctx->current_u;
425:   (*ctx->func)(ctx->funcctx,U,a);
426:   if (ctx->funcisetbase) {
427:     (*ctx->funcisetbase)(ctx->funcctx,U);
428:   }
429:   VecCopy(U,w);

431:   VecGetOwnershipRange(a,&rstart,&rend);
432:   VecGetArray(a,&aa);
433:   for (i=rstart; i<rend; i++) {
434:     VecGetArray(w,&ww);
435:     h    = ww[i-rstart];
436:     if (h == 0.0) h = 1.0;
437:     if (PetscAbsScalar(h) < umin && PetscRealPart(h) >= 0.0)     h = umin;
438:     else if (PetscRealPart(h) < 0.0 && PetscAbsScalar(h) < umin) h = -umin;
439:     h *= epsilon;

441:     ww[i-rstart] += h;
442:     VecRestoreArray(w,&ww);
443:     (*ctx->funci)(ctx->funcctx,i,w,&v);
444:     aa[i-rstart]  = (v - aa[i-rstart])/h;

446:     VecGetArray(w,&ww);
447:     ww[i-rstart] -= h;
448:     VecRestoreArray(w,&ww);
449:   }
450:   VecRestoreArray(a,&aa);
451:   return(0);
452: }

454: PETSC_EXTERN PetscErrorCode MatMFFDSetBase_MFFD(Mat J,Vec U,Vec F)
455: {
457:   MatMFFD        ctx;

460:   MatShellGetContext(J,&ctx);
461:   MatMFFDResetHHistory(J);
462:   if (!ctx->current_u) {
463:     VecDuplicate(U,&ctx->current_u);
464:     VecLockReadPush(ctx->current_u);
465:   }
466:   VecLockReadPop(ctx->current_u);
467:   VecCopy(U,ctx->current_u);
468:   VecLockReadPush(ctx->current_u);
469:   if (F) {
470:     if (ctx->current_f_allocated) {VecDestroy(&ctx->current_f);}
471:     ctx->current_f           = F;
472:     ctx->current_f_allocated = PETSC_FALSE;
473:   } else if (!ctx->current_f_allocated) {
474:     MatCreateVecs(J,NULL,&ctx->current_f);

476:     ctx->current_f_allocated = PETSC_TRUE;
477:   }
478:   if (!ctx->w) {
479:     VecDuplicate(ctx->current_u,&ctx->w);
480:   }
481:   J->assembled = PETSC_TRUE;
482:   return(0);
483: }

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

487: static PetscErrorCode  MatMFFDSetCheckh_MFFD(Mat J,FCN3 fun,void *ectx)
488: {
489:   MatMFFD        ctx;

493:   MatShellGetContext(J,&ctx);
494:   ctx->checkh    = fun;
495:   ctx->checkhctx = ectx;
496:   return(0);
497: }

499: /*@C
500:    MatMFFDSetOptionsPrefix - Sets the prefix used for searching for all
501:    MatMFFD options in the database.

503:    Collective on Mat

505:    Input Parameter:
506: +  A - the Mat context
507: -  prefix - the prefix to prepend to all option names

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

513:    Level: advanced

515: .seealso: MatSetFromOptions(), MatCreateSNESMF(), MatCreateMFFD()
516: @*/
517: PetscErrorCode  MatMFFDSetOptionsPrefix(Mat mat,const char prefix[])

519: {
520:   MatMFFD        mfctx;

525:   MatShellGetContext(mat,&mfctx);
527:   PetscObjectSetOptionsPrefix((PetscObject)mfctx,prefix);
528:   return(0);
529: }

531: static PetscErrorCode  MatSetFromOptions_MFFD(PetscOptionItems *PetscOptionsObject,Mat mat)
532: {
533:   MatMFFD        mfctx;
535:   PetscBool      flg;
536:   char           ftype[256];

540:   MatShellGetContext(mat,&mfctx);
542:   PetscObjectOptionsBegin((PetscObject)mfctx);
543:   PetscOptionsFList("-mat_mffd_type","Matrix free type","MatMFFDSetType",MatMFFDList,((PetscObject)mfctx)->type_name,ftype,256,&flg);
544:   if (flg) {
545:     MatMFFDSetType(mat,ftype);
546:   }

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

551:   flg  = PETSC_FALSE;
552:   PetscOptionsBool("-mat_mffd_check_positivity","Insure that U + h*a is nonnegative","MatMFFDSetCheckh",flg,&flg,NULL);
553:   if (flg) {
554:     MatMFFDSetCheckh(mat,MatMFFDCheckPositivity,0);
555:   }
556: #if defined(PETSC_USE_COMPLEX)
557:   PetscOptionsBool("-mat_mffd_complex","Use Lyness complex number trick to compute the matrix-vector product","None",mfctx->usecomplex,&mfctx->usecomplex,NULL);
558: #endif
559:   if (mfctx->ops->setfromoptions) {
560:     (*mfctx->ops->setfromoptions)(PetscOptionsObject,mfctx);
561:   }
562:   PetscOptionsEnd();
563:   return(0);
564: }

566: static PetscErrorCode  MatMFFDSetPeriod_MFFD(Mat mat,PetscInt period)
567: {
568:   MatMFFD        ctx;

572:   MatShellGetContext(mat,&ctx);
573:   ctx->recomputeperiod = period;
574:   return(0);
575: }

577: static PetscErrorCode  MatMFFDSetFunction_MFFD(Mat mat,PetscErrorCode (*func)(void*,Vec,Vec),void *funcctx)
578: {
579:   MatMFFD        ctx;

583:   MatShellGetContext(mat,&ctx);
584:   ctx->func    = func;
585:   ctx->funcctx = funcctx;
586:   return(0);
587: }

589: static PetscErrorCode  MatMFFDSetFunctionError_MFFD(Mat mat,PetscReal error)
590: {
591:   MatMFFD        ctx;

595:   MatShellGetContext(mat,&ctx);
596:   if (error != PETSC_DEFAULT) ctx->error_rel = error;
597:   return(0);
598: }

600: PetscErrorCode  MatMFFDSetHHistory_MFFD(Mat J,PetscScalar history[],PetscInt nhistory)
601: {
602:   MatMFFD        ctx;

606:   MatShellGetContext(J,&ctx);
607:   ctx->historyh    = history;
608:   ctx->maxcurrenth = nhistory;
609:   ctx->currenth    = 0.;
610:   return(0);
611: }

613: /*MC
614:   MATMFFD - MATMFFD = "mffd" - A matrix free matrix type.

616:   Level: advanced

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

620: .seealso: MatCreateMFFD(), MatCreateSNESMF(), MatMFFDSetFunction(), MatMFFDSetType(),  
621:           MatMFFDSetFunctionError(), MatMFFDDSSetUmin(), MatMFFDSetFunction()
622:           MatMFFDSetHHistory(), MatMFFDResetHHistory(), MatCreateSNESMF(),
623:           MatMFFDGetH(),
624: M*/
625: PETSC_EXTERN PetscErrorCode MatCreate_MFFD(Mat A)
626: {
627:   MatMFFD        mfctx;

631:   MatMFFDInitializePackage();

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

635:   mfctx->error_rel                = PETSC_SQRT_MACHINE_EPSILON;
636:   mfctx->recomputeperiod          = 1;
637:   mfctx->count                    = 0;
638:   mfctx->currenth                 = 0.0;
639:   mfctx->historyh                 = NULL;
640:   mfctx->ncurrenth                = 0;
641:   mfctx->maxcurrenth              = 0;
642:   ((PetscObject)mfctx)->type_name = 0;

644:   /*
645:      Create the empty data structure to contain compute-h routines.
646:      These will be filled in below from the command line options or
647:      a later call with MatMFFDSetType() or if that is not called
648:      then it will default in the first use of MatMult_MFFD()
649:   */
650:   mfctx->ops->compute        = 0;
651:   mfctx->ops->destroy        = 0;
652:   mfctx->ops->view           = 0;
653:   mfctx->ops->setfromoptions = 0;
654:   mfctx->hctx                = 0;

656:   mfctx->func    = 0;
657:   mfctx->funcctx = 0;
658:   mfctx->w       = NULL;
659:   mfctx->mat     = A;

661:   MatSetType(A,MATSHELL);
662:   MatShellSetContext(A,mfctx);
663:   MatShellSetOperation(A,MATOP_MULT,(void (*)(void))MatMult_MFFD);
664:   MatShellSetOperation(A,MATOP_DESTROY,(void (*)(void))MatDestroy_MFFD);
665:   MatShellSetOperation(A,MATOP_VIEW,(void (*)(void))MatView_MFFD);
666:   MatShellSetOperation(A,MATOP_ASSEMBLY_END,(void (*)(void))MatAssemblyEnd_MFFD);
667:   MatShellSetOperation(A,MATOP_SET_FROM_OPTIONS,(void (*)(void))MatSetFromOptions_MFFD);

669:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetBase_C",MatMFFDSetBase_MFFD);
670:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctioniBase_C",MatMFFDSetFunctioniBase_MFFD);
671:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctioni_C",MatMFFDSetFunctioni_MFFD);
672:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunction_C",MatMFFDSetFunction_MFFD);
673:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetCheckh_C",MatMFFDSetCheckh_MFFD);
674:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetPeriod_C",MatMFFDSetPeriod_MFFD);
675:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctionError_C",MatMFFDSetFunctionError_MFFD);
676:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDResetHHistory_C",MatMFFDResetHHistory_MFFD);
677:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetHHistory_C",MatMFFDSetHHistory_MFFD);
678:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetType_C",MatMFFDSetType_MFFD);
679:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDGetH_C",MatMFFDGetH_MFFD);
680:   PetscObjectChangeTypeName((PetscObject)A,MATMFFD);
681:   return(0);
682: }

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

687:    Collective on Vec

689:    Input Parameters:
690: +  comm - MPI communicator
691: .  m - number of local rows (or PETSC_DECIDE to have calculated if M is given)
692:            This value should be the same as the local size used in creating the
693:            y vector for the matrix-vector product y = Ax.
694: .  n - This value should be the same as the local size used in creating the
695:        x vector for the matrix-vector product y = Ax. (or PETSC_DECIDE to have
696:        calculated if N is given) For square matrices n is almost always m.
697: .  M - number of global rows (or PETSC_DETERMINE to have calculated if m is given)
698: -  N - number of global columns (or PETSC_DETERMINE to have calculated if n is given)


701:    Output Parameter:
702: .  J - the matrix-free matrix

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


713:    Level: advanced

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

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

722: .vb
723:      F'(u)*a = [F(u+h*a) - F(u)]/h where
724:      h = error_rel*u'a/||a||^2                        if  |u'a| > umin*||a||_{1}
725:        = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2   otherwise
726:  where
727:      error_rel = square root of relative error in function evaluation
728:      umin = minimum iterate parameter
729: .ve

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

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

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

740:    Options Database Keys:
741: +  -mat_mffd_err <error_rel> - Sets error_rel
742: .  -mat_mffd_unim <umin> - Sets umin (for default PETSc routine that computes h only)
743: -  -mat_mffd_check_positivity

745: .seealso: MatDestroy(), MatMFFDSetFunctionError(), MatMFFDDSSetUmin(), MatMFFDSetFunction()
746:           MatMFFDSetHHistory(), MatMFFDResetHHistory(), MatCreateSNESMF(),
747:           MatMFFDGetH(), MatMFFDRegister(), MatMFFDComputeJacobian()

749: @*/
750: PetscErrorCode  MatCreateMFFD(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,Mat *J)
751: {

755:   MatCreate(comm,J);
756:   MatSetSizes(*J,m,n,M,N);
757:   MatSetType(*J,MATMFFD);
758:   MatSetUp(*J);
759:   return(0);
760: }

762: /*@
763:    MatMFFDGetH - Gets the last value that was used as the differencing
764:    parameter.

766:    Not Collective

768:    Input Parameters:
769: .  mat - the matrix obtained with MatCreateSNESMF()

771:    Output Paramter:
772: .  h - the differencing step size

774:    Level: advanced

776: .seealso: MatCreateSNESMF(),MatMFFDSetHHistory(), MatCreateMFFD(), MATMFFD, MatMFFDResetHHistory()
777: @*/
778: PetscErrorCode  MatMFFDGetH(Mat mat,PetscScalar *h)
779: {

785:   PetscUseMethod(mat,"MatMFFDGetH_C",(Mat,PetscScalar*),(mat,h));
786:   return(0);
787: }

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

792:    Logically Collective on Mat

794:    Input Parameters:
795: +  mat - the matrix free matrix created via MatCreateSNESMF() or MatCreateMFFD()
796: .  func - the function to use
797: -  funcctx - optional function context passed to function

799:    Calling Sequence of func:
800: $     func (void *funcctx, Vec x, Vec f)

802: +  funcctx - user provided context
803: .  x - input vector
804: -  f - computed output function

806:    Level: advanced

808:    Notes:
809:     If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
810:     matrix inside your compute Jacobian routine

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

814: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD,
815:           MatMFFDSetHHistory(), MatMFFDResetHHistory(), SNESetFunction()
816: @*/
817: PetscErrorCode  MatMFFDSetFunction(Mat mat,PetscErrorCode (*func)(void*,Vec,Vec),void *funcctx)
818: {

823:   PetscTryMethod(mat,"MatMFFDSetFunction_C",(Mat,PetscErrorCode (*)(void*,Vec,Vec),void*),(mat,func,funcctx));
824:   return(0);
825: }

827: /*@C
828:    MatMFFDSetFunctioni - Sets the function for a single component

830:    Logically Collective on Mat

832:    Input Parameters:
833: +  mat - the matrix free matrix created via MatCreateSNESMF()
834: -  funci - the function to use

836:    Level: advanced

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

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

846: @*/
847: PetscErrorCode  MatMFFDSetFunctioni(Mat mat,PetscErrorCode (*funci)(void*,PetscInt,Vec,PetscScalar*))
848: {

853:   PetscTryMethod(mat,"MatMFFDSetFunctioni_C",(Mat,PetscErrorCode (*)(void*,PetscInt,Vec,PetscScalar*)),(mat,funci));
854:   return(0);
855: }

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

860:    Logically Collective on Mat

862:    Input Parameters:
863: +  mat - the matrix free matrix created via MatCreateSNESMF()
864: -  func - the function to use

866:    Level: advanced

868:    Notes:
869:     If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
870:     matrix inside your compute Jacobian routine.
871:     This function is necessary to compute the diagonal of the matrix.


874: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD
875:           MatMFFDSetHHistory(), MatMFFDResetHHistory(), SNESetFunction(), MatGetDiagonal()
876: @*/
877: PetscErrorCode  MatMFFDSetFunctioniBase(Mat mat,PetscErrorCode (*func)(void*,Vec))
878: {

883:   PetscTryMethod(mat,"MatMFFDSetFunctioniBase_C",(Mat,PetscErrorCode (*)(void*,Vec)),(mat,func));
884:   return(0);
885: }

887: /*@
888:    MatMFFDSetPeriod - Sets how often h is recomputed, by default it is everytime

890:    Logically Collective on Mat

892:    Input Parameters:
893: +  mat - the matrix free matrix created via MatCreateSNESMF()
894: -  period - 1 for everytime, 2 for every second etc

896:    Options Database Keys:
897: .  -mat_mffd_period <period>

899:    Level: advanced


902: .seealso: MatCreateSNESMF(),MatMFFDGetH(),
903:           MatMFFDSetHHistory(), MatMFFDResetHHistory()
904: @*/
905: PetscErrorCode  MatMFFDSetPeriod(Mat mat,PetscInt period)
906: {

912:   PetscTryMethod(mat,"MatMFFDSetPeriod_C",(Mat,PetscInt),(mat,period));
913:   return(0);
914: }

916: /*@
917:    MatMFFDSetFunctionError - Sets the error_rel for the approximation of
918:    matrix-vector products using finite differences.

920:    Logically Collective on Mat

922:    Input Parameters:
923: +  mat - the matrix free matrix created via MatCreateMFFD() or MatCreateSNESMF()
924: -  error_rel - relative error (should be set to the square root of
925:                the relative error in the function evaluations)

927:    Options Database Keys:
928: .  -mat_mffd_err <error_rel> - Sets error_rel

930:    Level: advanced

932:    Notes:
933:    The default matrix-free matrix-vector product routine computes
934: .vb
935:      F'(u)*a = [F(u+h*a) - F(u)]/h where
936:      h = error_rel*u'a/||a||^2                        if  |u'a| > umin*||a||_{1}
937:        = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2   else
938: .ve

940: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD
941:           MatMFFDSetHHistory(), MatMFFDResetHHistory()
942: @*/
943: PetscErrorCode  MatMFFDSetFunctionError(Mat mat,PetscReal error)
944: {

950:   PetscTryMethod(mat,"MatMFFDSetFunctionError_C",(Mat,PetscReal),(mat,error));
951:   return(0);
952: }

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

958:    Logically Collective on Mat

960:    Input Parameters:
961: +  J - the matrix-free matrix context
962: .  histroy - space to hold the history
963: -  nhistory - number of entries in history, if more entries are generated than
964:               nhistory, then the later ones are discarded

966:    Level: advanced

968:    Notes:
969:    Use MatMFFDResetHHistory() to reset the history counter and collect
970:    a new batch of differencing parameters, h.

972: .seealso: MatMFFDGetH(), MatCreateSNESMF(),
973:           MatMFFDResetHHistory(), MatMFFDSetFunctionError()

975: @*/
976: PetscErrorCode  MatMFFDSetHHistory(Mat J,PetscScalar history[],PetscInt nhistory)
977: {

984:   PetscUseMethod(J,"MatMFFDSetHHistory_C",(Mat,PetscScalar[],PetscInt),(J,history,nhistory));
985:   return(0);
986: }

988: /*@
989:    MatMFFDResetHHistory - Resets the counter to zero to begin
990:    collecting a new set of differencing histories.

992:    Logically Collective on Mat

994:    Input Parameters:
995: .  J - the matrix-free matrix context

997:    Level: advanced

999:    Notes:
1000:    Use MatMFFDSetHHistory() to create the original history counter.

1002: .seealso: MatMFFDGetH(), MatCreateSNESMF(),
1003:           MatMFFDSetHHistory(), MatMFFDSetFunctionError()

1005: @*/
1006: PetscErrorCode  MatMFFDResetHHistory(Mat J)
1007: {

1012:   PetscTryMethod(J,"MatMFFDResetHHistory_C",(Mat),(J));
1013:   return(0);
1014: }

1016: /*@
1017:     MatMFFDSetBase - Sets the vector U at which matrix vector products of the
1018:         Jacobian are computed

1020:     Logically Collective on Mat

1022:     Input Parameters:
1023: +   J - the MatMFFD matrix
1024: .   U - the vector
1025: -   F - (optional) vector that contains F(u) if it has been already computed

1027:     Notes:
1028:     This is rarely used directly

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

1033:     Level: advanced

1035: @*/
1036: PetscErrorCode  MatMFFDSetBase(Mat J,Vec U,Vec F)
1037: {

1044:   PetscTryMethod(J,"MatMFFDSetBase_C",(Mat,Vec,Vec),(J,U,F));
1045:   return(0);
1046: }

1048: /*@C
1049:     MatMFFDSetCheckh - Sets a function that checks the computed h and adjusts
1050:         it to satisfy some criteria

1052:     Logically Collective on Mat

1054:     Input Parameters:
1055: +   J - the MatMFFD matrix
1056: .   fun - the function that checks h
1057: -   ctx - any context needed by the function

1059:     Options Database Keys:
1060: .   -mat_mffd_check_positivity

1062:     Level: advanced

1064:     Notes:
1065:     For example, MatMFFDCheckPositivity() insures that all entries
1066:        of U + h*a are non-negative

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

1071: .seealso:  MatMFFDCheckPositivity()
1072: @*/
1073: PetscErrorCode  MatMFFDSetCheckh(Mat J,PetscErrorCode (*fun)(void*,Vec,Vec,PetscScalar*),void *ctx)
1074: {

1079:   PetscTryMethod(J,"MatMFFDSetCheckh_C",(Mat,PetscErrorCode (*)(void*,Vec,Vec,PetscScalar*),void*),(J,fun,ctx));
1080:   return(0);
1081: }

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

1087:     Logically Collective on Vec

1089:     Input Parameters:
1090: +   U - base vector that is added to
1091: .   a - vector that is added
1092: .   h - scaling factor on a
1093: -   dummy - context variable (unused)

1095:     Options Database Keys:
1096: .   -mat_mffd_check_positivity

1098:     Level: advanced

1100:     Notes:
1101:     This is rarely used directly, rather it is passed as an argument to
1102:            MatMFFDSetCheckh()

1104: .seealso:  MatMFFDSetCheckh()
1105: @*/
1106: PetscErrorCode  MatMFFDCheckPositivity(void *dummy,Vec U,Vec a,PetscScalar *h)
1107: {
1108:   PetscReal      val, minval;
1109:   PetscScalar    *u_vec, *a_vec;
1111:   PetscInt       i,n;
1112:   MPI_Comm       comm;

1118:   PetscObjectGetComm((PetscObject)U,&comm);
1119:   VecGetArray(U,&u_vec);
1120:   VecGetArray(a,&a_vec);
1121:   VecGetLocalSize(U,&n);
1122:   minval = PetscAbsScalar(*h)*PetscRealConstant(1.01);
1123:   for (i=0; i<n; i++) {
1124:     if (PetscRealPart(u_vec[i] + *h*a_vec[i]) <= 0.0) {
1125:       val = PetscAbsScalar(u_vec[i]/a_vec[i]);
1126:       if (val < minval) minval = val;
1127:     }
1128:   }
1129:   VecRestoreArray(U,&u_vec);
1130:   VecRestoreArray(a,&a_vec);
1131:   MPIU_Allreduce(&minval,&val,1,MPIU_REAL,MPIU_MIN,comm);
1132:   if (val <= PetscAbsScalar(*h)) {
1133:     PetscInfo2(U,"Scaling back h from %g to %g\n",(double)PetscRealPart(*h),(double)(.99*val));
1134:     if (PetscRealPart(*h) > 0.0) *h =  0.99*val;
1135:     else                         *h = -0.99*val;
1136:   }
1137:   return(0);
1138: }