Actual source code: mffd.c

petsc-master 2017-02-25
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: .keywords: Petsc, destroy, package
 19: .seealso: PetscFinalize(), MatCreateMFFD(), MatCreateSNESMF()
 20: @*/
 21: PetscErrorCode  MatMFFDFinalizePackage(void)
 22: {

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

 32: /*@C
 33:   MatMFFDInitializePackage - This function initializes everything in the MatMFFD package. It is called
 34:   from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to MatCreate_MFFD()
 35:   when using static libraries.

 37:   Level: developer

 39: .keywords: Vec, initialize, package
 40: .seealso: PetscInitialize()
 41: @*/
 42: PetscErrorCode  MatMFFDInitializePackage(void)
 43: {
 44:   char           logList[256];
 45:   char           *className;
 46:   PetscBool      opt;

 50:   if (MatMFFDPackageInitialized) return(0);
 51:   MatMFFDPackageInitialized = PETSC_TRUE;
 52:   /* Register Classes */
 53:   PetscClassIdRegister("MatMFFD",&MATMFFD_CLASSID);
 54:   /* Register Constructors */
 55:   MatMFFDRegisterAll();
 56:   /* Register Events */
 57:   PetscLogEventRegister("MatMult MF",          MATMFFD_CLASSID,&MATMFFD_Mult);

 59:   /* Process info exclusions */
 60:   PetscOptionsGetString(NULL,NULL, "-info_exclude", logList, 256, &opt);
 61:   if (opt) {
 62:     PetscStrstr(logList, "matmffd", &className);
 63:     if (className) {
 64:       PetscInfoDeactivateClass(MATMFFD_CLASSID);
 65:     }
 66:   }
 67:   /* Process summary exclusions */
 68:   PetscOptionsGetString(NULL,NULL, "-log_exclude", logList, 256, &opt);
 69:   if (opt) {
 70:     PetscStrstr(logList, "matmffd", &className);
 71:     if (className) {
 72:       PetscLogEventDeactivateClass(MATMFFD_CLASSID);
 73:     }
 74:   }
 75:   PetscRegisterFinalize(MatMFFDFinalizePackage);
 76:   return(0);
 77: }

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

 83:     Input Parameters:
 84: +   mat - the "matrix-free" matrix created via MatCreateSNESMF(), or MatCreateMFFD()
 85:           or MatSetType(mat,MATMFFD);
 86: -   ftype - the type requested, either MATMFFD_WP or MATMFFD_DS

 88:     Level: advanced

 90:     Notes:
 91:     For example, such routines can compute h for use in
 92:     Jacobian-vector products of the form

 94:                         F(x+ha) - F(x)
 95:           F'(u)a  ~=  ----------------
 96:                               h

 98: .seealso: MatCreateSNESMF(), MatMFFDRegister(), MatMFFDSetFunction(), MatCreateMFFD()
 99: @*/
100: PetscErrorCode  MatMFFDSetType(Mat mat,MatMFFDType ftype)
101: {
102:   PetscErrorCode ierr,(*r)(MatMFFD);
103:   MatMFFD        ctx = (MatMFFD)mat->data;
104:   PetscBool      match;


110:   PetscObjectTypeCompare((PetscObject)mat,MATMFFD,&match);
111:   if (!match) return(0);

113:   /* already set, so just return */
114:   PetscObjectTypeCompare((PetscObject)ctx,ftype,&match);
115:   if (match) return(0);

117:   /* destroy the old one if it exists */
118:   if (ctx->ops->destroy) {
119:     (*ctx->ops->destroy)(ctx);
120:   }

122:    PetscFunctionListFind(MatMFFDList,ftype,&r);
123:   if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown MatMFFD type %s given",ftype);
124:   (*r)(ctx);
125:   PetscObjectChangeTypeName((PetscObject)ctx,ftype);
126:   return(0);
127: }

129: typedef PetscErrorCode (*FCN1)(void*,Vec); /* force argument to next function to not be extern C*/
130: static PetscErrorCode  MatMFFDSetFunctioniBase_MFFD(Mat mat,FCN1 func)
131: {
132:   MatMFFD ctx = (MatMFFD)mat->data;

135:   ctx->funcisetbase = func;
136:   return(0);
137: }

139: typedef PetscErrorCode (*FCN2)(void*,PetscInt,Vec,PetscScalar*); /* force argument to next function to not be extern C*/
140: static PetscErrorCode  MatMFFDSetFunctioni_MFFD(Mat mat,FCN2 funci)
141: {
142:   MatMFFD ctx = (MatMFFD)mat->data;

145:   ctx->funci = funci;
146:   return(0);
147: }

149: static PetscErrorCode  MatMFFDResetHHistory_MFFD(Mat J)
150: {
151:   MatMFFD ctx = (MatMFFD)J->data;

154:   ctx->ncurrenth = 0;
155:   return(0);
156: }

158: /*@C
159:    MatMFFDRegister - Adds a method to the MatMFFD registry.

161:    Not Collective

163:    Input Parameters:
164: +  name_solver - name of a new user-defined compute-h module
165: -  routine_create - routine to create method context

167:    Level: developer

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

172:    Sample usage:
173: .vb
174:    MatMFFDRegister("my_h",MyHCreate);
175: .ve

177:    Then, your solver can be chosen with the procedural interface via
178: $     MatMFFDSetType(mfctx,"my_h")
179:    or at runtime via the option
180: $     -mat_mffd_type my_h

182: .keywords: MatMFFD, register

184: .seealso: MatMFFDRegisterAll(), MatMFFDRegisterDestroy()
185:  @*/
186: PetscErrorCode  MatMFFDRegister(const char sname[],PetscErrorCode (*function)(MatMFFD))
187: {

191:   PetscFunctionListAdd(&MatMFFDList,sname,function);
192:   return(0);
193: }

195: /* ----------------------------------------------------------------------------------------*/
196: static PetscErrorCode MatDestroy_MFFD(Mat mat)
197: {
199:   MatMFFD        ctx = (MatMFFD)mat->data;

202:   VecDestroy(&ctx->w);
203:   VecDestroy(&ctx->drscale);
204:   VecDestroy(&ctx->dlscale);
205:   VecDestroy(&ctx->dshift);
206:   if (ctx->current_f_allocated) {
207:     VecDestroy(&ctx->current_f);
208:   }
209:   if (ctx->ops->destroy) {(*ctx->ops->destroy)(ctx);}
210:   PetscHeaderDestroy(&ctx);
211:   mat->data = 0;

213:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetBase_C",NULL);
214:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctioniBase_C",NULL);
215:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctioni_C",NULL);
216:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunction_C",NULL);
217:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctionError_C",NULL);
218:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetCheckh_C",NULL);
219:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetPeriod_C",NULL);
220:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDResetHHistory_C",NULL);
221:   return(0);
222: }

224: /*
225:    MatMFFDView_MFFD - Views matrix-free parameters.

227: */
228: static PetscErrorCode MatView_MFFD(Mat J,PetscViewer viewer)
229: {
231:   MatMFFD        ctx = (MatMFFD)J->data;
232:   PetscBool      iascii, viewbase, viewfunction;
233:   const char     *prefix;

236:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
237:   if (iascii) {
238:     PetscViewerASCIIPrintf(viewer,"Matrix-free approximation:\n");
239:     PetscViewerASCIIPushTab(viewer);
240:     PetscViewerASCIIPrintf(viewer,"err=%g (relative error in function evaluation)\n",(double)ctx->error_rel);
241:     if (!((PetscObject)ctx)->type_name) {
242:       PetscViewerASCIIPrintf(viewer,"The compute h routine has not yet been set\n");
243:     } else {
244:       PetscViewerASCIIPrintf(viewer,"Using %s compute h routine\n",((PetscObject)ctx)->type_name);
245:     }
246:     if (ctx->ops->view) {
247:       (*ctx->ops->view)(ctx,viewer);
248:     }
249:     PetscObjectGetOptionsPrefix((PetscObject)J, &prefix);

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

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

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

282:   MatMFFDResetHHistory(J);
283:   j->vshift = 0.0;
284:   j->vscale = 1.0;
285:   return(0);
286: }

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

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

305:   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");
306:   /* We log matrix-free matrix-vector products separately, so that we can
307:      separate the performance monitoring from the cases that use conventional
308:      storage.  We may eventually modify event logging to associate events
309:      with particular objects, hence alleviating the more general problem. */
310:   PetscLogEventBegin(MATMFFD_Mult,a,y,0,0);

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

328:   if (mat->erroriffailure && PetscIsInfOrNanScalar(h)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Computed Nan differencing parameter h");
329:   if (ctx->checkh) {
330:     (*ctx->checkh)(ctx->checkhctx,U,a,&h);
331:   }

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

345:   /* w = u + ha */
346:   if (ctx->drscale) {
347:     VecPointwiseMult(ctx->drscale,a,U);
348:     VecAYPX(U,h,w);
349:   } else {
350:     VecWAXPY(w,h,a,U);
351:   }

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

359:   VecAXPY(y,-1.0,F);
360:   VecScale(y,1.0/h);

362:   if ((ctx->vshift != 0.0) || (ctx->vscale != 1.0)) {
363:     VecAXPBY(y,ctx->vshift,ctx->vscale,a);
364:   }
365:   if (ctx->dlscale) {
366:     VecPointwiseMult(y,ctx->dlscale,y);
367:   }
368:   if (ctx->dshift) {
369:     VecPointwiseMult(ctx->dshift,a,U);
370:     VecAXPY(y,1.0,U);
371:   }

373:   if (mat->nullsp) {MatNullSpaceRemove(mat->nullsp,y);}

375:   PetscLogEventEnd(MATMFFD_Mult,a,y,0,0);
376:   return(0);
377: }

379: /*
380:   MatGetDiagonal_MFFD - Gets the diagonal for a matrix free matrix

382:         y ~= (F(u + ha) - F(u))/h,
383:   where F = nonlinear function, as set by SNESSetFunction()
384:         u = current iterate
385:         h = difference interval
386: */
387: static PetscErrorCode MatGetDiagonal_MFFD(Mat mat,Vec a)
388: {
389:   MatMFFD        ctx = (MatMFFD)mat->data;
390:   PetscScalar    h,*aa,*ww,v;
391:   PetscReal      epsilon = PETSC_SQRT_MACHINE_EPSILON,umin = 100.0*PETSC_SQRT_MACHINE_EPSILON;
392:   Vec            w,U;
394:   PetscInt       i,rstart,rend;

397:   if (!ctx->funci) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Requires calling MatMFFDSetFunctioni() first");

399:   w    = ctx->w;
400:   U    = ctx->current_u;
401:   (*ctx->func)(ctx->funcctx,U,a);
402:   (*ctx->funcisetbase)(ctx->funcctx,U);
403:   VecCopy(U,w);

405:   VecGetOwnershipRange(a,&rstart,&rend);
406:   VecGetArray(a,&aa);
407:   for (i=rstart; i<rend; i++) {
408:     VecGetArray(w,&ww);
409:     h    = ww[i-rstart];
410:     if (h == 0.0) h = 1.0;
411:     if (PetscAbsScalar(h) < umin && PetscRealPart(h) >= 0.0)     h = umin;
412:     else if (PetscRealPart(h) < 0.0 && PetscAbsScalar(h) < umin) h = -umin;
413:     h *= epsilon;

415:     ww[i-rstart] += h;
416:     VecRestoreArray(w,&ww);
417:     (*ctx->funci)(ctx->funcctx,i,w,&v);
418:     aa[i-rstart]  = (v - aa[i-rstart])/h;

420:     /* possibly shift and scale result */
421:     if ((ctx->vshift != 0.0) || (ctx->vscale != 1.0)) {
422:       aa[i - rstart] = ctx->vshift + ctx->vscale*aa[i-rstart];
423:     }

425:     VecGetArray(w,&ww);
426:     ww[i-rstart] -= h;
427:     VecRestoreArray(w,&ww);
428:   }
429:   VecRestoreArray(a,&aa);
430:   return(0);
431: }

433: static PetscErrorCode MatDiagonalScale_MFFD(Mat mat,Vec ll,Vec rr)
434: {
435:   MatMFFD        aij = (MatMFFD)mat->data;

439:   if (ll && !aij->dlscale) {
440:     VecDuplicate(ll,&aij->dlscale);
441:   }
442:   if (rr && !aij->drscale) {
443:     VecDuplicate(rr,&aij->drscale);
444:   }
445:   if (ll) {
446:     VecCopy(ll,aij->dlscale);
447:   }
448:   if (rr) {
449:     VecCopy(rr,aij->drscale);
450:   }
451:   return(0);
452: }

454: static PetscErrorCode MatDiagonalSet_MFFD(Mat mat,Vec ll,InsertMode mode)
455: {
456:   MatMFFD        aij = (MatMFFD)mat->data;

460:   if (mode == INSERT_VALUES) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"No diagonal set with INSERT_VALUES");
461:   if (!aij->dshift) {
462:     VecDuplicate(ll,&aij->dshift);
463:   }
464:   VecAXPY(aij->dshift,1.0,ll);
465:   return(0);
466: }

468: static PetscErrorCode MatShift_MFFD(Mat Y,PetscScalar a)
469: {
470:   MatMFFD shell = (MatMFFD)Y->data;

473:   shell->vshift += a;
474:   return(0);
475: }

477: static PetscErrorCode MatScale_MFFD(Mat Y,PetscScalar a)
478: {
479:   MatMFFD shell = (MatMFFD)Y->data;

482:   shell->vscale *= a;
483:   return(0);
484: }

486: PETSC_EXTERN PetscErrorCode MatMFFDSetBase_MFFD(Mat J,Vec U,Vec F)
487: {
489:   MatMFFD        ctx = (MatMFFD)J->data;

492:   MatMFFDResetHHistory(J);

494:   ctx->current_u = U;
495:   if (F) {
496:     if (ctx->current_f_allocated) {VecDestroy(&ctx->current_f);}
497:     ctx->current_f           = F;
498:     ctx->current_f_allocated = PETSC_FALSE;
499:   } else if (!ctx->current_f_allocated) {
500:     MatCreateVecs(J,NULL,&ctx->current_f);

502:     ctx->current_f_allocated = PETSC_TRUE;
503:   }
504:   if (!ctx->w) {
505:     VecDuplicate(ctx->current_u, &ctx->w);
506:   }
507:   J->assembled = PETSC_TRUE;
508:   return(0);
509: }

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

513: static PetscErrorCode  MatMFFDSetCheckh_MFFD(Mat J,FCN3 fun,void *ectx)
514: {
515:   MatMFFD ctx = (MatMFFD)J->data;

518:   ctx->checkh    = fun;
519:   ctx->checkhctx = ectx;
520:   return(0);
521: }

523: /*@C
524:    MatMFFDSetOptionsPrefix - Sets the prefix used for searching for all
525:    MatMFFD options in the database.

527:    Collective on Mat

529:    Input Parameter:
530: +  A - the Mat context
531: -  prefix - the prefix to prepend to all option names

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

537:    Level: advanced

539: .keywords: SNES, matrix-free, parameters

541: .seealso: MatSetFromOptions(), MatCreateSNESMF(), MatCreateMFFD()
542: @*/
543: PetscErrorCode  MatMFFDSetOptionsPrefix(Mat mat,const char prefix[])

545: {
546:   MatMFFD        mfctx = mat ? (MatMFFD)mat->data : (MatMFFD)NULL;

552:   PetscObjectSetOptionsPrefix((PetscObject)mfctx,prefix);
553:   return(0);
554: }

556: static PetscErrorCode  MatSetFromOptions_MFFD(PetscOptionItems *PetscOptionsObject,Mat mat)
557: {
558:   MatMFFD        mfctx = (MatMFFD)mat->data;
560:   PetscBool      flg;
561:   char           ftype[256];

566:   PetscObjectOptionsBegin((PetscObject)mfctx);
567:   PetscOptionsFList("-mat_mffd_type","Matrix free type","MatMFFDSetType",MatMFFDList,((PetscObject)mfctx)->type_name,ftype,256,&flg);
568:   if (flg) {
569:     MatMFFDSetType(mat,ftype);
570:   }

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

575:   flg  = PETSC_FALSE;
576:   PetscOptionsBool("-mat_mffd_check_positivity","Insure that U + h*a is nonnegative","MatMFFDSetCheckh",flg,&flg,NULL);
577:   if (flg) {
578:     MatMFFDSetCheckh(mat,MatMFFDCheckPositivity,0);
579:   }
580:   if (mfctx->ops->setfromoptions) {
581:     (*mfctx->ops->setfromoptions)(PetscOptionsObject,mfctx);
582:   }
583:   PetscOptionsEnd();
584:   return(0);
585: }

587: static PetscErrorCode  MatMFFDSetPeriod_MFFD(Mat mat,PetscInt period)
588: {
589:   MatMFFD ctx = (MatMFFD)mat->data;

593:   ctx->recomputeperiod = period;
594:   return(0);
595: }

597: static PetscErrorCode  MatMFFDSetFunction_MFFD(Mat mat,PetscErrorCode (*func)(void*,Vec,Vec),void *funcctx)
598: {
599:   MatMFFD ctx = (MatMFFD)mat->data;

602:   ctx->func    = func;
603:   ctx->funcctx = funcctx;
604:   return(0);
605: }

607: static PetscErrorCode  MatMFFDSetFunctionError_MFFD(Mat mat,PetscReal error)
608: {
609:   MatMFFD ctx = (MatMFFD)mat->data;

613:   if (error != PETSC_DEFAULT) ctx->error_rel = error;
614:   return(0);
615: }

617: static PetscErrorCode MatMissingDiagonal_MFFD(Mat A,PetscBool  *missing,PetscInt *d)
618: {
620:   *missing = PETSC_FALSE;
621:   return(0);
622: }

624: /*MC
625:   MATMFFD - MATMFFD = "mffd" - A matrix free matrix type.

627:   Level: advanced

629: .seealso: MatCreateMFFD(), MatCreateSNESMF(), MatMFFDSetFunction(), MatMFFDSetType(),  
630:           MatMFFDSetFunctionError(), MatMFFDDSSetUmin(), MatMFFDSetFunction()
631:           MatMFFDSetHHistory(), MatMFFDResetHHistory(), MatCreateSNESMF(),
632:           MatMFFDGetH(),
633: M*/
634: PETSC_EXTERN PetscErrorCode MatCreate_MFFD(Mat A)
635: {
636:   MatMFFD        mfctx;

640:   MatMFFDInitializePackage();

642:   PetscHeaderCreate(mfctx,MATMFFD_CLASSID,"MatMFFD","Matrix-free Finite Differencing","Mat",PetscObjectComm((PetscObject)A),MatDestroy_MFFD,MatView_MFFD);

644:   mfctx->error_rel                = PETSC_SQRT_MACHINE_EPSILON;
645:   mfctx->recomputeperiod          = 1;
646:   mfctx->count                    = 0;
647:   mfctx->currenth                 = 0.0;
648:   mfctx->historyh                 = NULL;
649:   mfctx->ncurrenth                = 0;
650:   mfctx->maxcurrenth              = 0;
651:   ((PetscObject)mfctx)->type_name = 0;

653:   mfctx->vshift = 0.0;
654:   mfctx->vscale = 1.0;

656:   /*
657:      Create the empty data structure to contain compute-h routines.
658:      These will be filled in below from the command line options or
659:      a later call with MatMFFDSetType() or if that is not called
660:      then it will default in the first use of MatMult_MFFD()
661:   */
662:   mfctx->ops->compute        = 0;
663:   mfctx->ops->destroy        = 0;
664:   mfctx->ops->view           = 0;
665:   mfctx->ops->setfromoptions = 0;
666:   mfctx->hctx                = 0;

668:   mfctx->func    = 0;
669:   mfctx->funcctx = 0;
670:   mfctx->w       = NULL;

672:   A->data = mfctx;

674:   A->ops->mult            = MatMult_MFFD;
675:   A->ops->destroy         = MatDestroy_MFFD;
676:   A->ops->view            = MatView_MFFD;
677:   A->ops->assemblyend     = MatAssemblyEnd_MFFD;
678:   A->ops->getdiagonal     = MatGetDiagonal_MFFD;
679:   A->ops->scale           = MatScale_MFFD;
680:   A->ops->shift           = MatShift_MFFD;
681:   A->ops->diagonalscale   = MatDiagonalScale_MFFD;
682:   A->ops->diagonalset     = MatDiagonalSet_MFFD;
683:   A->ops->setfromoptions  = MatSetFromOptions_MFFD;
684:   A->ops->missingdiagonal = MatMissingDiagonal_MFFD;
685:   A->assembled           = PETSC_TRUE;

687:   PetscLayoutSetUp(A->rmap);
688:   PetscLayoutSetUp(A->cmap);

690:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetBase_C",MatMFFDSetBase_MFFD);
691:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctioniBase_C",MatMFFDSetFunctioniBase_MFFD);
692:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctioni_C",MatMFFDSetFunctioni_MFFD);
693:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunction_C",MatMFFDSetFunction_MFFD);
694:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetCheckh_C",MatMFFDSetCheckh_MFFD);
695:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetPeriod_C",MatMFFDSetPeriod_MFFD);
696:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctionError_C",MatMFFDSetFunctionError_MFFD);
697:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDResetHHistory_C",MatMFFDResetHHistory_MFFD);

699:   mfctx->mat = A;

701:   PetscObjectChangeTypeName((PetscObject)A,MATMFFD);
702:   return(0);
703: }

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

708:    Collective on Vec

710:    Input Parameters:
711: +  comm - MPI communicator
712: .  m - number of local rows (or PETSC_DECIDE to have calculated if M is given)
713:            This value should be the same as the local size used in creating the
714:            y vector for the matrix-vector product y = Ax.
715: .  n - This value should be the same as the local size used in creating the
716:        x vector for the matrix-vector product y = Ax. (or PETSC_DECIDE to have
717:        calculated if N is given) For square matrices n is almost always m.
718: .  M - number of global rows (or PETSC_DETERMINE to have calculated if m is given)
719: -  N - number of global columns (or PETSC_DETERMINE to have calculated if n is given)


722:    Output Parameter:
723: .  J - the matrix-free matrix

725:    Options Database Keys: call MatSetFromOptions() to trigger these
726: +  -mat_mffd_type - wp or ds (see MATMFFD_WP or MATMFFD_DS)
727: -  -mat_mffd_err - square root of estimated relative error in function evaluation
728: -  -mat_mffd_period - how often h is recomputed, defaults to 1, everytime


731:    Level: advanced

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

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

740: .vb
741:      F'(u)*a = [F(u+h*a) - F(u)]/h where
742:      h = error_rel*u'a/||a||^2                        if  |u'a| > umin*||a||_{1}
743:        = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2   otherwise
744:  where
745:      error_rel = square root of relative error in function evaluation
746:      umin = minimum iterate parameter
747: .ve

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

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

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

758:    Options Database Keys:
759: +  -mat_mffd_err <error_rel> - Sets error_rel
760: .  -mat_mffd_unim <umin> - Sets umin (for default PETSc routine that computes h only)
761: -  -mat_mffd_check_positivity

763: .keywords: default, matrix-free, create, matrix

765: .seealso: MatDestroy(), MatMFFDSetFunctionError(), MatMFFDDSSetUmin(), MatMFFDSetFunction()
766:           MatMFFDSetHHistory(), MatMFFDResetHHistory(), MatCreateSNESMF(),
767:           MatMFFDGetH(), MatMFFDRegister(), MatMFFDComputeJacobian()

769: @*/
770: PetscErrorCode  MatCreateMFFD(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,Mat *J)
771: {

775:   MatCreate(comm,J);
776:   MatSetSizes(*J,m,n,M,N);
777:   MatSetType(*J,MATMFFD);
778:   MatSetUp(*J);
779:   return(0);
780: }


783: /*@
784:    MatMFFDGetH - Gets the last value that was used as the differencing
785:    parameter.

787:    Not Collective

789:    Input Parameters:
790: .  mat - the matrix obtained with MatCreateSNESMF()

792:    Output Paramter:
793: .  h - the differencing step size

795:    Level: advanced

797: .keywords: SNES, matrix-free, parameters

799: .seealso: MatCreateSNESMF(),MatMFFDSetHHistory(), MatCreateMFFD(), MATMFFD, MatMFFDResetHHistory()
800: @*/
801: PetscErrorCode  MatMFFDGetH(Mat mat,PetscScalar *h)
802: {
803:   MatMFFD        ctx = (MatMFFD)mat->data;
805:   PetscBool      match;

808:   PetscObjectTypeCompare((PetscObject)mat,MATMFFD,&match);
809:   if (!match) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONG,"Not a MFFD matrix");

811:   *h = ctx->currenth;
812:   return(0);
813: }

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

818:    Logically Collective on Mat

820:    Input Parameters:
821: +  mat - the matrix free matrix created via MatCreateSNESMF() or MatCreateMFFD()
822: .  func - the function to use
823: -  funcctx - optional function context passed to function

825:    Calling Sequence of func:
826: $     func (void *funcctx, Vec x, Vec f)

828: +  funcctx - user provided context
829: .  x - input vector
830: -  f - computed output function

832:    Level: advanced

834:    Notes:
835:     If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
836:     matrix inside your compute Jacobian routine

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

840: .keywords: SNES, matrix-free, function

842: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD,
843:           MatMFFDSetHHistory(), MatMFFDResetHHistory(), SNESetFunction()
844: @*/
845: PetscErrorCode  MatMFFDSetFunction(Mat mat,PetscErrorCode (*func)(void*,Vec,Vec),void *funcctx)
846: {

850:   PetscTryMethod(mat,"MatMFFDSetFunction_C",(Mat,PetscErrorCode (*)(void*,Vec,Vec),void*),(mat,func,funcctx));
851:   return(0);
852: }

854: /*@C
855:    MatMFFDSetFunctioni - Sets the function for a single component

857:    Logically Collective on Mat

859:    Input Parameters:
860: +  mat - the matrix free matrix created via MatCreateSNESMF()
861: -  funci - the function to use

863:    Level: advanced

865:    Notes:
866:     If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
867:     matrix inside your compute Jacobian routine


870: .keywords: SNES, matrix-free, function

872: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatMFFDSetHHistory(), MatMFFDResetHHistory(), SNESetFunction()

874: @*/
875: PetscErrorCode  MatMFFDSetFunctioni(Mat mat,PetscErrorCode (*funci)(void*,PetscInt,Vec,PetscScalar*))
876: {

881:   PetscTryMethod(mat,"MatMFFDSetFunctioni_C",(Mat,PetscErrorCode (*)(void*,PetscInt,Vec,PetscScalar*)),(mat,funci));
882:   return(0);
883: }


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

889:    Logically Collective on Mat

891:    Input Parameters:
892: +  mat - the matrix free matrix created via MatCreateSNESMF()
893: -  func - the function to use

895:    Level: advanced

897:    Notes:
898:     If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
899:     matrix inside your compute Jacobian routine


902: .keywords: SNES, matrix-free, function

904: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD
905:           MatMFFDSetHHistory(), MatMFFDResetHHistory(), SNESetFunction()
906: @*/
907: PetscErrorCode  MatMFFDSetFunctioniBase(Mat mat,PetscErrorCode (*func)(void*,Vec))
908: {

913:   PetscTryMethod(mat,"MatMFFDSetFunctioniBase_C",(Mat,PetscErrorCode (*)(void*,Vec)),(mat,func));
914:   return(0);
915: }

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

920:    Logically Collective on Mat

922:    Input Parameters:
923: +  mat - the matrix free matrix created via MatCreateSNESMF()
924: -  period - 1 for everytime, 2 for every second etc

926:    Options Database Keys:
927: +  -mat_mffd_period <period>

929:    Level: advanced


932: .keywords: SNES, matrix-free, parameters

934: .seealso: MatCreateSNESMF(),MatMFFDGetH(),
935:           MatMFFDSetHHistory(), MatMFFDResetHHistory()
936: @*/
937: PetscErrorCode  MatMFFDSetPeriod(Mat mat,PetscInt period)
938: {

942:   PetscTryMethod(mat,"MatMFFDSetPeriod_C",(Mat,PetscInt),(mat,period));
943:   return(0);
944: }

946: /*@
947:    MatMFFDSetFunctionError - Sets the error_rel for the approximation of
948:    matrix-vector products using finite differences.

950:    Logically Collective on Mat

952:    Input Parameters:
953: +  mat - the matrix free matrix created via MatCreateMFFD() or MatCreateSNESMF()
954: -  error_rel - relative error (should be set to the square root of
955:                the relative error in the function evaluations)

957:    Options Database Keys:
958: +  -mat_mffd_err <error_rel> - Sets error_rel

960:    Level: advanced

962:    Notes:
963:    The default matrix-free matrix-vector product routine computes
964: .vb
965:      F'(u)*a = [F(u+h*a) - F(u)]/h where
966:      h = error_rel*u'a/||a||^2                        if  |u'a| > umin*||a||_{1}
967:        = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2   else
968: .ve

970: .keywords: SNES, matrix-free, parameters

972: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD
973:           MatMFFDSetHHistory(), MatMFFDResetHHistory()
974: @*/
975: PetscErrorCode  MatMFFDSetFunctionError(Mat mat,PetscReal error)
976: {

980:   PetscTryMethod(mat,"MatMFFDSetFunctionError_C",(Mat,PetscReal),(mat,error));
981:   return(0);
982: }

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

988:    Logically Collective on Mat

990:    Input Parameters:
991: +  J - the matrix-free matrix context
992: .  histroy - space to hold the history
993: -  nhistory - number of entries in history, if more entries are generated than
994:               nhistory, then the later ones are discarded

996:    Level: advanced

998:    Notes:
999:    Use MatMFFDResetHHistory() to reset the history counter and collect
1000:    a new batch of differencing parameters, h.

1002: .keywords: SNES, matrix-free, h history, differencing history

1004: .seealso: MatMFFDGetH(), MatCreateSNESMF(),
1005:           MatMFFDResetHHistory(), MatMFFDSetFunctionError()

1007: @*/
1008: PetscErrorCode  MatMFFDSetHHistory(Mat J,PetscScalar history[],PetscInt nhistory)
1009: {
1010:   MatMFFD        ctx = (MatMFFD)J->data;
1012:   PetscBool      match;

1015:   PetscObjectTypeCompare((PetscObject)J,MATMFFD,&match);
1016:   if (!match) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_WRONG,"Not a MFFD matrix");
1017:   ctx->historyh    = history;
1018:   ctx->maxcurrenth = nhistory;
1019:   ctx->currenth    = 0.;
1020:   return(0);
1021: }


1024: /*@
1025:    MatMFFDResetHHistory - Resets the counter to zero to begin
1026:    collecting a new set of differencing histories.

1028:    Logically Collective on Mat

1030:    Input Parameters:
1031: .  J - the matrix-free matrix context

1033:    Level: advanced

1035:    Notes:
1036:    Use MatMFFDSetHHistory() to create the original history counter.

1038: .keywords: SNES, matrix-free, h history, differencing history

1040: .seealso: MatMFFDGetH(), MatCreateSNESMF(),
1041:           MatMFFDSetHHistory(), MatMFFDSetFunctionError()

1043: @*/
1044: PetscErrorCode  MatMFFDResetHHistory(Mat J)
1045: {

1049:   PetscTryMethod(J,"MatMFFDResetHHistory_C",(Mat),(J));
1050:   return(0);
1051: }


1054: /*@
1055:     MatMFFDSetBase - Sets the vector U at which matrix vector products of the
1056:         Jacobian are computed

1058:     Logically Collective on Mat

1060:     Input Parameters:
1061: +   J - the MatMFFD matrix
1062: .   U - the vector
1063: -   F - (optional) vector that contains F(u) if it has been already computed

1065:     Notes: This is rarely used directly

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

1070:     Level: advanced

1072: @*/
1073: PetscErrorCode  MatMFFDSetBase(Mat J,Vec U,Vec F)
1074: {

1081:   PetscTryMethod(J,"MatMFFDSetBase_C",(Mat,Vec,Vec),(J,U,F));
1082:   return(0);
1083: }

1085: /*@C
1086:     MatMFFDSetCheckh - Sets a function that checks the computed h and adjusts
1087:         it to satisfy some criteria

1089:     Logically Collective on Mat

1091:     Input Parameters:
1092: +   J - the MatMFFD matrix
1093: .   fun - the function that checks h
1094: -   ctx - any context needed by the function

1096:     Options Database Keys:
1097: .   -mat_mffd_check_positivity

1099:     Level: advanced

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

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

1107: .seealso:  MatMFFDCheckPositivity()
1108: @*/
1109: PetscErrorCode  MatMFFDSetCheckh(Mat J,PetscErrorCode (*fun)(void*,Vec,Vec,PetscScalar*),void *ctx)
1110: {

1115:   PetscTryMethod(J,"MatMFFDSetCheckh_C",(Mat,PetscErrorCode (*)(void*,Vec,Vec,PetscScalar*),void*),(J,fun,ctx));
1116:   return(0);
1117: }

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

1123:     Logically Collective on Vec

1125:     Input Parameters:
1126: +   U - base vector that is added to
1127: .   a - vector that is added
1128: .   h - scaling factor on a
1129: -   dummy - context variable (unused)

1131:     Options Database Keys:
1132: .   -mat_mffd_check_positivity

1134:     Level: advanced

1136:     Notes: This is rarely used directly, rather it is passed as an argument to
1137:            MatMFFDSetCheckh()

1139: .seealso:  MatMFFDSetCheckh()
1140: @*/
1141: PetscErrorCode  MatMFFDCheckPositivity(void *dummy,Vec U,Vec a,PetscScalar *h)
1142: {
1143:   PetscReal      val, minval;
1144:   PetscScalar    *u_vec, *a_vec;
1146:   PetscInt       i,n;
1147:   MPI_Comm       comm;

1150:   PetscObjectGetComm((PetscObject)U,&comm);
1151:   VecGetArray(U,&u_vec);
1152:   VecGetArray(a,&a_vec);
1153:   VecGetLocalSize(U,&n);
1154:   minval = PetscAbsScalar(*h*1.01);
1155:   for (i=0; i<n; i++) {
1156:     if (PetscRealPart(u_vec[i] + *h*a_vec[i]) <= 0.0) {
1157:       val = PetscAbsScalar(u_vec[i]/a_vec[i]);
1158:       if (val < minval) minval = val;
1159:     }
1160:   }
1161:   VecRestoreArray(U,&u_vec);
1162:   VecRestoreArray(a,&a_vec);
1163:   MPIU_Allreduce(&minval,&val,1,MPIU_REAL,MPIU_MIN,comm);
1164:   if (val <= PetscAbsScalar(*h)) {
1165:     PetscInfo2(U,"Scaling back h from %g to %g\n",(double)PetscRealPart(*h),(double)(.99*val));
1166:     if (PetscRealPart(*h) > 0.0) *h =  0.99*val;
1167:     else                         *h = -0.99*val;
1168:   }
1169:   return(0);
1170: }