Actual source code: mffd.c

petsc-master 2019-03-20
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 MatInitializePackage().

 36:   Level: developer

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

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

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

 77:     Input Parameters:
 78: +   mat - the "matrix-free" matrix created via MatCreateSNESMF(), or MatCreateMFFD()
 79:           or MatSetType(mat,MATMFFD);
 80: -   ftype - the type requested, either MATMFFD_WP or MATMFFD_DS

 82:     Level: advanced

 84:     Notes:
 85:     For example, such routines can compute h for use in
 86:     Jacobian-vector products of the form

 88:                         F(x+ha) - F(x)
 89:           F'(u)a  ~=  ----------------
 90:                               h

 92: .seealso: MatCreateSNESMF(), MatMFFDRegister(), MatMFFDSetFunction(), MatCreateMFFD()
 93: @*/
 94: PetscErrorCode  MatMFFDSetType(Mat mat,MatMFFDType ftype)
 95: {
 96:   PetscErrorCode ierr,(*r)(MatMFFD);
 97:   MatMFFD        ctx = (MatMFFD)mat->data;
 98:   PetscBool      match;


104:   PetscObjectTypeCompare((PetscObject)mat,MATMFFD,&match);
105:   if (!match) return(0);

107:   /* already set, so just return */
108:   PetscObjectTypeCompare((PetscObject)ctx,ftype,&match);
109:   if (match) return(0);

111:   /* destroy the old one if it exists */
112:   if (ctx->ops->destroy) {
113:     (*ctx->ops->destroy)(ctx);
114:   }

116:    PetscFunctionListFind(MatMFFDList,ftype,&r);
117:   if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown MatMFFD type %s given",ftype);
118:   (*r)(ctx);
119:   PetscObjectChangeTypeName((PetscObject)ctx,ftype);
120:   return(0);
121: }

123: static PetscErrorCode MatGetDiagonal_MFFD(Mat,Vec);

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

131:   ctx->funcisetbase = func;
132:   /* allow users to compose their own getdiagonal and allow MatHasOperation
133:      to return false if the two functions pointers are not set */
134:   if (!mat->ops->getdiagonal && func) {
135:     mat->ops->getdiagonal = MatGetDiagonal_MFFD;
136:   }
137:   return(0);
138: }

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

146:   ctx->funci = funci;
147:   /* allow users to compose their own getdiagonal and allow MatHasOperation
148:      to return false if the two functions pointers are not set */
149:   if (!mat->ops->getdiagonal && funci) {
150:     mat->ops->getdiagonal = MatGetDiagonal_MFFD;
151:   }
152:   return(0);
153: }

155: static PetscErrorCode  MatMFFDResetHHistory_MFFD(Mat J)
156: {
157:   MatMFFD ctx = (MatMFFD)J->data;

160:   ctx->ncurrenth = 0;
161:   return(0);
162: }

164: /*@C
165:    MatMFFDRegister - Adds a method to the MatMFFD registry.

167:    Not Collective

169:    Input Parameters:
170: +  name_solver - name of a new user-defined compute-h module
171: -  routine_create - routine to create method context

173:    Level: developer

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

178:    Sample usage:
179: .vb
180:    MatMFFDRegister("my_h",MyHCreate);
181: .ve

183:    Then, your solver can be chosen with the procedural interface via
184: $     MatMFFDSetType(mfctx,"my_h")
185:    or at runtime via the option
186: $     -mat_mffd_type my_h

188: .keywords: MatMFFD, register

190: .seealso: MatMFFDRegisterAll(), MatMFFDRegisterDestroy()
191:  @*/
192: PetscErrorCode  MatMFFDRegister(const char sname[],PetscErrorCode (*function)(MatMFFD))
193: {

197:   MatInitializePackage();
198:   PetscFunctionListAdd(&MatMFFDList,sname,function);
199:   return(0);
200: }

202: /* ----------------------------------------------------------------------------------------*/
203: static PetscErrorCode MatDestroy_MFFD(Mat mat)
204: {
206:   MatMFFD        ctx = (MatMFFD)mat->data;

209:   VecDestroy(&ctx->w);
210:   VecDestroy(&ctx->drscale);
211:   VecDestroy(&ctx->dlscale);
212:   VecDestroy(&ctx->dshift);
213:   VecDestroy(&ctx->dshiftw);
214:   VecDestroy(&ctx->current_u);
215:   if (ctx->current_f_allocated) {
216:     VecDestroy(&ctx->current_f);
217:   }
218:   if (ctx->ops->destroy) {(*ctx->ops->destroy)(ctx);}
219:   PetscHeaderDestroy(&ctx);
220:   mat->data = 0;

222:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetBase_C",NULL);
223:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctioniBase_C",NULL);
224:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctioni_C",NULL);
225:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunction_C",NULL);
226:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctionError_C",NULL);
227:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetCheckh_C",NULL);
228:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetPeriod_C",NULL);
229:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDResetHHistory_C",NULL);
230:   return(0);
231: }

233: /*
234:    MatMFFDView_MFFD - Views matrix-free parameters.

236: */
237: static PetscErrorCode MatView_MFFD(Mat J,PetscViewer viewer)
238: {
240:   MatMFFD        ctx = (MatMFFD)J->data;
241:   PetscBool      iascii, viewbase, viewfunction;
242:   const char     *prefix;

245:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
246:   if (iascii) {
247:     PetscViewerASCIIPrintf(viewer,"Matrix-free approximation:\n");
248:     PetscViewerASCIIPushTab(viewer);
249:     PetscViewerASCIIPrintf(viewer,"err=%g (relative error in function evaluation)\n",(double)ctx->error_rel);
250:     if (!((PetscObject)ctx)->type_name) {
251:       PetscViewerASCIIPrintf(viewer,"The compute h routine has not yet been set\n");
252:     } else {
253:       PetscViewerASCIIPrintf(viewer,"Using %s compute h routine\n",((PetscObject)ctx)->type_name);
254:     }
255: #if defined(PETSC_USE_COMPLEX)
256:     if (ctx->usecomplex) {
257:       PetscViewerASCIIPrintf(viewer,"Using Lyness complex number trick to compute the matrix-vector product\n");
258:     }
259: #endif
260:     if (ctx->ops->view) {
261:       (*ctx->ops->view)(ctx,viewer);
262:     }
263:     PetscObjectGetOptionsPrefix((PetscObject)J, &prefix);

265:     PetscOptionsHasName(((PetscObject)J)->options,prefix, "-mat_mffd_view_base", &viewbase);
266:     if (viewbase) {
267:       PetscViewerASCIIPrintf(viewer, "Base:\n");
268:       VecView(ctx->current_u, viewer);
269:     }
270:     PetscOptionsHasName(((PetscObject)J)->options,prefix, "-mat_mffd_view_function", &viewfunction);
271:     if (viewfunction) {
272:       PetscViewerASCIIPrintf(viewer, "Function:\n");
273:       VecView(ctx->current_f, viewer);
274:     }
275:     PetscViewerASCIIPopTab(viewer);
276:   }
277:   return(0);
278: }

280: /*
281:    MatAssemblyEnd_MFFD - Resets the ctx->ncurrenth to zero. This
282:    allows the user to indicate the beginning of a new linear solve by calling
283:    MatAssemblyXXX() on the matrix free matrix. This then allows the
284:    MatCreateMFFD_WP() to properly compute ||U|| only the first time
285:    in the linear solver rather than every time.

287:    This function is referenced directly from MatAssemblyEnd_SNESMF(), which may be in a different shared library hence
288:    it must be labeled as PETSC_EXTERN
289: */
290: PETSC_EXTERN PetscErrorCode MatAssemblyEnd_MFFD(Mat J,MatAssemblyType mt)
291: {
293:   MatMFFD        j = (MatMFFD)J->data;

296:   MatMFFDResetHHistory(J);
297:   j->vshift = 0.0;
298:   j->vscale = 1.0;
299:   return(0);
300: }

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

305:         y ~= (F(u + ha) - F(u))/h,
306:   where F = nonlinear function, as set by SNESSetFunction()
307:         u = current iterate
308:         h = difference interval
309: */
310: static PetscErrorCode MatMult_MFFD(Mat mat,Vec a,Vec y)
311: {
312:   MatMFFD        ctx = (MatMFFD)mat->data;
313:   PetscScalar    h;
314:   Vec            w,U,F;
316:   PetscBool      zeroa;

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

326:   w = ctx->w;
327:   U = ctx->current_u;
328:   F = ctx->current_f;
329:   /*
330:       Compute differencing parameter
331:   */
332:   if (!((PetscObject)ctx)->type_name) {
333:     MatMFFDSetType(mat,MATMFFD_WP);
334:     MatSetFromOptions(mat);
335:   }
336:   (*ctx->ops->compute)(ctx,U,a,&h,&zeroa);
337:   if (zeroa) {
338:     VecSet(y,0.0);
339:     return(0);
340:   }

342:   if (mat->erroriffailure && PetscIsInfOrNanScalar(h)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Computed Nan differencing parameter h");
343:   if (ctx->checkh) {
344:     (*ctx->checkh)(ctx->checkhctx,U,a,&h);
345:   }

347:   /* keep a record of the current differencing parameter h */
348:   ctx->currenth = h;
349: #if defined(PETSC_USE_COMPLEX)
350:   PetscInfo2(mat,"Current differencing parameter: %g + %g i\n",(double)PetscRealPart(h),(double)PetscImaginaryPart(h));
351: #else
352:   PetscInfo1(mat,"Current differencing parameter: %15.12e\n",h);
353: #endif
354:   if (ctx->historyh && ctx->ncurrenth < ctx->maxcurrenth) {
355:     ctx->historyh[ctx->ncurrenth] = h;
356:   }
357:   ctx->ncurrenth++;

359: #if defined(PETSC_USE_COMPLEX)
360:   if (ctx->usecomplex) h = PETSC_i*h;
361: #endif
362: 
363:   /* w = u + ha */
364:   if (ctx->drscale) {
365:     VecPointwiseMult(ctx->drscale,a,U);
366:     VecAYPX(U,h,w);
367:   } else {
368:     VecWAXPY(w,h,a,U);
369:   }

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

377: #if defined(PETSC_USE_COMPLEX)  
378:   if (ctx->usecomplex) {
379:     VecImaginaryPart(y);
380:     h    = PetscImaginaryPart(h);
381:   } else {
382:     VecAXPY(y,-1.0,F);
383:   }
384: #else
385:   VecAXPY(y,-1.0,F);
386: #endif
387:   VecScale(y,1.0/h);

389:   if ((ctx->vshift != 0.0) || (ctx->vscale != 1.0)) {
390:     VecAXPBY(y,ctx->vshift,ctx->vscale,a);
391:   }
392:   if (ctx->dlscale) {
393:     VecPointwiseMult(y,ctx->dlscale,y);
394:   }
395:   if (ctx->dshift) {
396:     if (!ctx->dshiftw) {
397:       VecDuplicate(y,&ctx->dshiftw);
398:     }
399:     VecPointwiseMult(ctx->dshift,a,ctx->dshiftw);
400:     VecAXPY(y,1.0,ctx->dshiftw);
401:   }

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

405:   PetscLogEventEnd(MATMFFD_Mult,a,y,0,0);
406:   return(0);
407: }

409: /*
410:   MatGetDiagonal_MFFD - Gets the diagonal for a matrix free matrix

412:         y ~= (F(u + ha) - F(u))/h,
413:   where F = nonlinear function, as set by SNESSetFunction()
414:         u = current iterate
415:         h = difference interval
416: */
417: PetscErrorCode MatGetDiagonal_MFFD(Mat mat,Vec a)
418: {
419:   MatMFFD        ctx = (MatMFFD)mat->data;
420:   PetscScalar    h,*aa,*ww,v;
421:   PetscReal      epsilon = PETSC_SQRT_MACHINE_EPSILON,umin = 100.0*PETSC_SQRT_MACHINE_EPSILON;
422:   Vec            w,U;
424:   PetscInt       i,rstart,rend;

427:   if (!ctx->func) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Requires calling MatMFFDSetFunction() first");
428:   if (!ctx->funci) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Requires calling MatMFFDSetFunctioni() first");
429:   if (!ctx->funcisetbase) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Requires calling MatMFFDSetFunctioniBase() first");
430:   w    = ctx->w;
431:   U    = ctx->current_u;
432:   (*ctx->func)(ctx->funcctx,U,a);
433:   (*ctx->funcisetbase)(ctx->funcctx,U);
434:   VecCopy(U,w);

436:   VecGetOwnershipRange(a,&rstart,&rend);
437:   VecGetArray(a,&aa);
438:   for (i=rstart; i<rend; i++) {
439:     VecGetArray(w,&ww);
440:     h    = ww[i-rstart];
441:     if (h == 0.0) h = 1.0;
442:     if (PetscAbsScalar(h) < umin && PetscRealPart(h) >= 0.0)     h = umin;
443:     else if (PetscRealPart(h) < 0.0 && PetscAbsScalar(h) < umin) h = -umin;
444:     h *= epsilon;

446:     ww[i-rstart] += h;
447:     VecRestoreArray(w,&ww);
448:     (*ctx->funci)(ctx->funcctx,i,w,&v);
449:     aa[i-rstart]  = (v - aa[i-rstart])/h;

451:     /* possibly shift and scale result */
452:     if ((ctx->vshift != 0.0) || (ctx->vscale != 1.0)) {
453:       aa[i - rstart] = ctx->vshift + ctx->vscale*aa[i-rstart];
454:     }

456:     VecGetArray(w,&ww);
457:     ww[i-rstart] -= h;
458:     VecRestoreArray(w,&ww);
459:   }
460:   VecRestoreArray(a,&aa);
461:   return(0);
462: }

464: static PetscErrorCode MatDiagonalScale_MFFD(Mat mat,Vec ll,Vec rr)
465: {
466:   MatMFFD        aij = (MatMFFD)mat->data;

470:   if (ll && !aij->dlscale) {
471:     VecDuplicate(ll,&aij->dlscale);
472:   }
473:   if (rr && !aij->drscale) {
474:     VecDuplicate(rr,&aij->drscale);
475:   }
476:   if (ll) {
477:     VecCopy(ll,aij->dlscale);
478:   }
479:   if (rr) {
480:     VecCopy(rr,aij->drscale);
481:   }
482:   return(0);
483: }

485: static PetscErrorCode MatDiagonalSet_MFFD(Mat mat,Vec ll,InsertMode mode)
486: {
487:   MatMFFD        aij = (MatMFFD)mat->data;

491:   if (mode == INSERT_VALUES) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"No diagonal set with INSERT_VALUES");
492:   if (!aij->dshift) {
493:     VecDuplicate(ll,&aij->dshift);
494:   }
495:   VecAXPY(aij->dshift,1.0,ll);
496:   return(0);
497: }

499: static PetscErrorCode MatShift_MFFD(Mat Y,PetscScalar a)
500: {
501:   MatMFFD shell = (MatMFFD)Y->data;

504:   shell->vshift += a;
505:   return(0);
506: }

508: static PetscErrorCode MatScale_MFFD(Mat Y,PetscScalar a)
509: {
510:   MatMFFD shell = (MatMFFD)Y->data;

513:   shell->vscale *= a;
514:   return(0);
515: }

517: PETSC_EXTERN PetscErrorCode MatMFFDSetBase_MFFD(Mat J,Vec U,Vec F)
518: {
520:   MatMFFD        ctx = (MatMFFD)J->data;

523:   MatMFFDResetHHistory(J);
524:   if (!ctx->current_u) {
525:     VecDuplicate(U,&ctx->current_u);
526:     VecLockReadPush(ctx->current_u);
527:   }
528:   VecLockReadPop(ctx->current_u);
529:   VecCopy(U,ctx->current_u);
530:   VecLockReadPush(ctx->current_u);
531:   if (F) {
532:     if (ctx->current_f_allocated) {VecDestroy(&ctx->current_f);}
533:     ctx->current_f           = F;
534:     ctx->current_f_allocated = PETSC_FALSE;
535:   } else if (!ctx->current_f_allocated) {
536:     MatCreateVecs(J,NULL,&ctx->current_f);

538:     ctx->current_f_allocated = PETSC_TRUE;
539:   }
540:   if (!ctx->w) {
541:     VecDuplicate(ctx->current_u,&ctx->w);
542:   }
543:   J->assembled = PETSC_TRUE;
544:   return(0);
545: }

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

549: static PetscErrorCode  MatMFFDSetCheckh_MFFD(Mat J,FCN3 fun,void *ectx)
550: {
551:   MatMFFD ctx = (MatMFFD)J->data;

554:   ctx->checkh    = fun;
555:   ctx->checkhctx = ectx;
556:   return(0);
557: }

559: /*@C
560:    MatMFFDSetOptionsPrefix - Sets the prefix used for searching for all
561:    MatMFFD options in the database.

563:    Collective on Mat

565:    Input Parameter:
566: +  A - the Mat context
567: -  prefix - the prefix to prepend to all option names

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

573:    Level: advanced

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

577: .seealso: MatSetFromOptions(), MatCreateSNESMF(), MatCreateMFFD()
578: @*/
579: PetscErrorCode  MatMFFDSetOptionsPrefix(Mat mat,const char prefix[])

581: {
582:   MatMFFD        mfctx = mat ? (MatMFFD)mat->data : (MatMFFD)NULL;

588:   PetscObjectSetOptionsPrefix((PetscObject)mfctx,prefix);
589:   return(0);
590: }

592: static PetscErrorCode  MatSetFromOptions_MFFD(PetscOptionItems *PetscOptionsObject,Mat mat)
593: {
594:   MatMFFD        mfctx = (MatMFFD)mat->data;
596:   PetscBool      flg;
597:   char           ftype[256];

602:   PetscObjectOptionsBegin((PetscObject)mfctx);
603:   PetscOptionsFList("-mat_mffd_type","Matrix free type","MatMFFDSetType",MatMFFDList,((PetscObject)mfctx)->type_name,ftype,256,&flg);
604:   if (flg) {
605:     MatMFFDSetType(mat,ftype);
606:   }

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

611:   flg  = PETSC_FALSE;
612:   PetscOptionsBool("-mat_mffd_check_positivity","Insure that U + h*a is nonnegative","MatMFFDSetCheckh",flg,&flg,NULL);
613:   if (flg) {
614:     MatMFFDSetCheckh(mat,MatMFFDCheckPositivity,0);
615:   }
616: #if defined(PETSC_USE_COMPLEX)
617:   PetscOptionsBool("-mat_mffd_complex","Use Lyness complex number trick to compute the matrix-vector product","None",mfctx->usecomplex,&mfctx->usecomplex,NULL);
618: #endif
619:   if (mfctx->ops->setfromoptions) {
620:     (*mfctx->ops->setfromoptions)(PetscOptionsObject,mfctx);
621:   }
622:   PetscOptionsEnd();
623:   return(0);
624: }

626: static PetscErrorCode  MatMFFDSetPeriod_MFFD(Mat mat,PetscInt period)
627: {
628:   MatMFFD ctx = (MatMFFD)mat->data;

631:   ctx->recomputeperiod = period;
632:   return(0);
633: }

635: static PetscErrorCode  MatMFFDSetFunction_MFFD(Mat mat,PetscErrorCode (*func)(void*,Vec,Vec),void *funcctx)
636: {
637:   MatMFFD ctx = (MatMFFD)mat->data;

640:   ctx->func    = func;
641:   ctx->funcctx = funcctx;
642:   return(0);
643: }

645: static PetscErrorCode  MatMFFDSetFunctionError_MFFD(Mat mat,PetscReal error)
646: {
647:   MatMFFD ctx = (MatMFFD)mat->data;

650:   if (error != PETSC_DEFAULT) ctx->error_rel = error;
651:   return(0);
652: }

654: static PetscErrorCode MatMissingDiagonal_MFFD(Mat A,PetscBool  *missing,PetscInt *d)
655: {
657:   *missing = PETSC_FALSE;
658:   return(0);
659: }

661: /*MC
662:   MATMFFD - MATMFFD = "mffd" - A matrix free matrix type.

664:   Level: advanced

666: .seealso: MatCreateMFFD(), MatCreateSNESMF(), MatMFFDSetFunction(), MatMFFDSetType(),  
667:           MatMFFDSetFunctionError(), MatMFFDDSSetUmin(), MatMFFDSetFunction()
668:           MatMFFDSetHHistory(), MatMFFDResetHHistory(), MatCreateSNESMF(),
669:           MatMFFDGetH(),
670: M*/
671: PETSC_EXTERN PetscErrorCode MatCreate_MFFD(Mat A)
672: {
673:   MatMFFD        mfctx;

677:   MatMFFDInitializePackage();

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

681:   mfctx->error_rel                = PETSC_SQRT_MACHINE_EPSILON;
682:   mfctx->recomputeperiod          = 1;
683:   mfctx->count                    = 0;
684:   mfctx->currenth                 = 0.0;
685:   mfctx->historyh                 = NULL;
686:   mfctx->ncurrenth                = 0;
687:   mfctx->maxcurrenth              = 0;
688:   ((PetscObject)mfctx)->type_name = 0;

690:   mfctx->vshift = 0.0;
691:   mfctx->vscale = 1.0;

693:   /*
694:      Create the empty data structure to contain compute-h routines.
695:      These will be filled in below from the command line options or
696:      a later call with MatMFFDSetType() or if that is not called
697:      then it will default in the first use of MatMult_MFFD()
698:   */
699:   mfctx->ops->compute        = 0;
700:   mfctx->ops->destroy        = 0;
701:   mfctx->ops->view           = 0;
702:   mfctx->ops->setfromoptions = 0;
703:   mfctx->hctx                = 0;

705:   mfctx->func    = 0;
706:   mfctx->funcctx = 0;
707:   mfctx->w       = NULL;

709:   A->data = mfctx;

711:   A->ops->mult            = MatMult_MFFD;
712:   A->ops->destroy         = MatDestroy_MFFD;
713:   A->ops->view            = MatView_MFFD;
714:   A->ops->assemblyend     = MatAssemblyEnd_MFFD;
715:   A->ops->scale           = MatScale_MFFD;
716:   A->ops->shift           = MatShift_MFFD;
717:   A->ops->diagonalscale   = MatDiagonalScale_MFFD;
718:   A->ops->diagonalset     = MatDiagonalSet_MFFD;
719:   A->ops->setfromoptions  = MatSetFromOptions_MFFD;
720:   A->ops->missingdiagonal = MatMissingDiagonal_MFFD;
721:   A->assembled            = PETSC_TRUE;

723:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetBase_C",MatMFFDSetBase_MFFD);
724:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctioniBase_C",MatMFFDSetFunctioniBase_MFFD);
725:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctioni_C",MatMFFDSetFunctioni_MFFD);
726:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunction_C",MatMFFDSetFunction_MFFD);
727:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetCheckh_C",MatMFFDSetCheckh_MFFD);
728:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetPeriod_C",MatMFFDSetPeriod_MFFD);
729:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctionError_C",MatMFFDSetFunctionError_MFFD);
730:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDResetHHistory_C",MatMFFDResetHHistory_MFFD);

732:   mfctx->mat = A;

734:   PetscObjectChangeTypeName((PetscObject)A,MATMFFD);
735:   return(0);
736: }

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

741:    Collective on Vec

743:    Input Parameters:
744: +  comm - MPI communicator
745: .  m - number of local rows (or PETSC_DECIDE to have calculated if M is given)
746:            This value should be the same as the local size used in creating the
747:            y vector for the matrix-vector product y = Ax.
748: .  n - This value should be the same as the local size used in creating the
749:        x vector for the matrix-vector product y = Ax. (or PETSC_DECIDE to have
750:        calculated if N is given) For square matrices n is almost always m.
751: .  M - number of global rows (or PETSC_DETERMINE to have calculated if m is given)
752: -  N - number of global columns (or PETSC_DETERMINE to have calculated if n is given)


755:    Output Parameter:
756: .  J - the matrix-free matrix

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


767:    Level: advanced

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

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

776: .vb
777:      F'(u)*a = [F(u+h*a) - F(u)]/h where
778:      h = error_rel*u'a/||a||^2                        if  |u'a| > umin*||a||_{1}
779:        = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2   otherwise
780:  where
781:      error_rel = square root of relative error in function evaluation
782:      umin = minimum iterate parameter
783: .ve

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

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

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

794:    Options Database Keys:
795: +  -mat_mffd_err <error_rel> - Sets error_rel
796: .  -mat_mffd_unim <umin> - Sets umin (for default PETSc routine that computes h only)
797: -  -mat_mffd_check_positivity

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

801: .seealso: MatDestroy(), MatMFFDSetFunctionError(), MatMFFDDSSetUmin(), MatMFFDSetFunction()
802:           MatMFFDSetHHistory(), MatMFFDResetHHistory(), MatCreateSNESMF(),
803:           MatMFFDGetH(), MatMFFDRegister(), MatMFFDComputeJacobian()

805: @*/
806: PetscErrorCode  MatCreateMFFD(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,Mat *J)
807: {

811:   MatCreate(comm,J);
812:   MatSetSizes(*J,m,n,M,N);
813:   MatSetType(*J,MATMFFD);
814:   MatSetUp(*J);
815:   return(0);
816: }

818: /*@
819:    MatMFFDGetH - Gets the last value that was used as the differencing
820:    parameter.

822:    Not Collective

824:    Input Parameters:
825: .  mat - the matrix obtained with MatCreateSNESMF()

827:    Output Paramter:
828: .  h - the differencing step size

830:    Level: advanced

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

834: .seealso: MatCreateSNESMF(),MatMFFDSetHHistory(), MatCreateMFFD(), MATMFFD, MatMFFDResetHHistory()
835: @*/
836: PetscErrorCode  MatMFFDGetH(Mat mat,PetscScalar *h)
837: {
838:   MatMFFD        ctx = (MatMFFD)mat->data;
840:   PetscBool      match;

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

848:   *h = ctx->currenth;
849:   return(0);
850: }

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

855:    Logically Collective on Mat

857:    Input Parameters:
858: +  mat - the matrix free matrix created via MatCreateSNESMF() or MatCreateMFFD()
859: .  func - the function to use
860: -  funcctx - optional function context passed to function

862:    Calling Sequence of func:
863: $     func (void *funcctx, Vec x, Vec f)

865: +  funcctx - user provided context
866: .  x - input vector
867: -  f - computed output function

869:    Level: advanced

871:    Notes:
872:     If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
873:     matrix inside your compute Jacobian routine

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

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

879: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD,
880:           MatMFFDSetHHistory(), MatMFFDResetHHistory(), SNESetFunction()
881: @*/
882: PetscErrorCode  MatMFFDSetFunction(Mat mat,PetscErrorCode (*func)(void*,Vec,Vec),void *funcctx)
883: {

888:   PetscTryMethod(mat,"MatMFFDSetFunction_C",(Mat,PetscErrorCode (*)(void*,Vec,Vec),void*),(mat,func,funcctx));
889:   return(0);
890: }

892: /*@C
893:    MatMFFDSetFunctioni - Sets the function for a single component

895:    Logically Collective on Mat

897:    Input Parameters:
898: +  mat - the matrix free matrix created via MatCreateSNESMF()
899: -  funci - the function to use

901:    Level: advanced

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

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

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

913: @*/
914: PetscErrorCode  MatMFFDSetFunctioni(Mat mat,PetscErrorCode (*funci)(void*,PetscInt,Vec,PetscScalar*))
915: {

920:   PetscTryMethod(mat,"MatMFFDSetFunctioni_C",(Mat,PetscErrorCode (*)(void*,PetscInt,Vec,PetscScalar*)),(mat,funci));
921:   return(0);
922: }

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

927:    Logically Collective on Mat

929:    Input Parameters:
930: +  mat - the matrix free matrix created via MatCreateSNESMF()
931: -  func - the function to use

933:    Level: advanced

935:    Notes:
936:     If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
937:     matrix inside your compute Jacobian routine.
938:     This function is necessary to compute the diagonal of the matrix.


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

943: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD
944:           MatMFFDSetHHistory(), MatMFFDResetHHistory(), SNESetFunction(), MatGetDiagonal()
945: @*/
946: PetscErrorCode  MatMFFDSetFunctioniBase(Mat mat,PetscErrorCode (*func)(void*,Vec))
947: {

952:   PetscTryMethod(mat,"MatMFFDSetFunctioniBase_C",(Mat,PetscErrorCode (*)(void*,Vec)),(mat,func));
953:   return(0);
954: }

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

959:    Logically Collective on Mat

961:    Input Parameters:
962: +  mat - the matrix free matrix created via MatCreateSNESMF()
963: -  period - 1 for everytime, 2 for every second etc

965:    Options Database Keys:
966: +  -mat_mffd_period <period>

968:    Level: advanced


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

973: .seealso: MatCreateSNESMF(),MatMFFDGetH(),
974:           MatMFFDSetHHistory(), MatMFFDResetHHistory()
975: @*/
976: PetscErrorCode  MatMFFDSetPeriod(Mat mat,PetscInt period)
977: {

983:   PetscTryMethod(mat,"MatMFFDSetPeriod_C",(Mat,PetscInt),(mat,period));
984:   return(0);
985: }

987: /*@
988:    MatMFFDSetFunctionError - Sets the error_rel for the approximation of
989:    matrix-vector products using finite differences.

991:    Logically Collective on Mat

993:    Input Parameters:
994: +  mat - the matrix free matrix created via MatCreateMFFD() or MatCreateSNESMF()
995: -  error_rel - relative error (should be set to the square root of
996:                the relative error in the function evaluations)

998:    Options Database Keys:
999: +  -mat_mffd_err <error_rel> - Sets error_rel

1001:    Level: advanced

1003:    Notes:
1004:    The default matrix-free matrix-vector product routine computes
1005: .vb
1006:      F'(u)*a = [F(u+h*a) - F(u)]/h where
1007:      h = error_rel*u'a/||a||^2                        if  |u'a| > umin*||a||_{1}
1008:        = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2   else
1009: .ve

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

1013: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD
1014:           MatMFFDSetHHistory(), MatMFFDResetHHistory()
1015: @*/
1016: PetscErrorCode  MatMFFDSetFunctionError(Mat mat,PetscReal error)
1017: {

1023:   PetscTryMethod(mat,"MatMFFDSetFunctionError_C",(Mat,PetscReal),(mat,error));
1024:   return(0);
1025: }

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

1031:    Logically Collective on Mat

1033:    Input Parameters:
1034: +  J - the matrix-free matrix context
1035: .  histroy - space to hold the history
1036: -  nhistory - number of entries in history, if more entries are generated than
1037:               nhistory, then the later ones are discarded

1039:    Level: advanced

1041:    Notes:
1042:    Use MatMFFDResetHHistory() to reset the history counter and collect
1043:    a new batch of differencing parameters, h.

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

1047: .seealso: MatMFFDGetH(), MatCreateSNESMF(),
1048:           MatMFFDResetHHistory(), MatMFFDSetFunctionError()

1050: @*/
1051: PetscErrorCode  MatMFFDSetHHistory(Mat J,PetscScalar history[],PetscInt nhistory)
1052: {
1053:   MatMFFD        ctx = (MatMFFD)J->data;
1055:   PetscBool      match;

1061:   PetscObjectTypeCompare((PetscObject)J,MATMFFD,&match);
1062:   if (!match) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_WRONG,"Not a MFFD matrix");
1063:   ctx->historyh    = history;
1064:   ctx->maxcurrenth = nhistory;
1065:   ctx->currenth    = 0.;
1066:   return(0);
1067: }

1069: /*@
1070:    MatMFFDResetHHistory - Resets the counter to zero to begin
1071:    collecting a new set of differencing histories.

1073:    Logically Collective on Mat

1075:    Input Parameters:
1076: .  J - the matrix-free matrix context

1078:    Level: advanced

1080:    Notes:
1081:    Use MatMFFDSetHHistory() to create the original history counter.

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

1085: .seealso: MatMFFDGetH(), MatCreateSNESMF(),
1086:           MatMFFDSetHHistory(), MatMFFDSetFunctionError()

1088: @*/
1089: PetscErrorCode  MatMFFDResetHHistory(Mat J)
1090: {

1095:   PetscTryMethod(J,"MatMFFDResetHHistory_C",(Mat),(J));
1096:   return(0);
1097: }

1099: /*@
1100:     MatMFFDSetBase - Sets the vector U at which matrix vector products of the
1101:         Jacobian are computed

1103:     Logically Collective on Mat

1105:     Input Parameters:
1106: +   J - the MatMFFD matrix
1107: .   U - the vector
1108: -   F - (optional) vector that contains F(u) if it has been already computed

1110:     Notes:
1111:     This is rarely used directly

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

1116:     Level: advanced

1118: @*/
1119: PetscErrorCode  MatMFFDSetBase(Mat J,Vec U,Vec F)
1120: {

1127:   PetscTryMethod(J,"MatMFFDSetBase_C",(Mat,Vec,Vec),(J,U,F));
1128:   return(0);
1129: }

1131: /*@C
1132:     MatMFFDSetCheckh - Sets a function that checks the computed h and adjusts
1133:         it to satisfy some criteria

1135:     Logically Collective on Mat

1137:     Input Parameters:
1138: +   J - the MatMFFD matrix
1139: .   fun - the function that checks h
1140: -   ctx - any context needed by the function

1142:     Options Database Keys:
1143: .   -mat_mffd_check_positivity

1145:     Level: advanced

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

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

1154: .seealso:  MatMFFDCheckPositivity()
1155: @*/
1156: PetscErrorCode  MatMFFDSetCheckh(Mat J,PetscErrorCode (*fun)(void*,Vec,Vec,PetscScalar*),void *ctx)
1157: {

1162:   PetscTryMethod(J,"MatMFFDSetCheckh_C",(Mat,PetscErrorCode (*)(void*,Vec,Vec,PetscScalar*),void*),(J,fun,ctx));
1163:   return(0);
1164: }

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

1170:     Logically Collective on Vec

1172:     Input Parameters:
1173: +   U - base vector that is added to
1174: .   a - vector that is added
1175: .   h - scaling factor on a
1176: -   dummy - context variable (unused)

1178:     Options Database Keys:
1179: .   -mat_mffd_check_positivity

1181:     Level: advanced

1183:     Notes:
1184:     This is rarely used directly, rather it is passed as an argument to
1185:            MatMFFDSetCheckh()

1187: .seealso:  MatMFFDSetCheckh()
1188: @*/
1189: PetscErrorCode  MatMFFDCheckPositivity(void *dummy,Vec U,Vec a,PetscScalar *h)
1190: {
1191:   PetscReal      val, minval;
1192:   PetscScalar    *u_vec, *a_vec;
1194:   PetscInt       i,n;
1195:   MPI_Comm       comm;

1201:   PetscObjectGetComm((PetscObject)U,&comm);
1202:   VecGetArray(U,&u_vec);
1203:   VecGetArray(a,&a_vec);
1204:   VecGetLocalSize(U,&n);
1205:   minval = PetscAbsScalar(*h)*PetscRealConstant(1.01);
1206:   for (i=0; i<n; i++) {
1207:     if (PetscRealPart(u_vec[i] + *h*a_vec[i]) <= 0.0) {
1208:       val = PetscAbsScalar(u_vec[i]/a_vec[i]);
1209:       if (val < minval) minval = val;
1210:     }
1211:   }
1212:   VecRestoreArray(U,&u_vec);
1213:   VecRestoreArray(a,&a_vec);
1214:   MPIU_Allreduce(&minval,&val,1,MPIU_REAL,MPIU_MIN,comm);
1215:   if (val <= PetscAbsScalar(*h)) {
1216:     PetscInfo2(U,"Scaling back h from %g to %g\n",(double)PetscRealPart(*h),(double)(.99*val));
1217:     if (PetscRealPart(*h) > 0.0) *h =  0.99*val;
1218:     else                         *h = -0.99*val;
1219:   }
1220:   return(0);
1221: }