Actual source code: mffd.c

petsc-3.9.2 2018-05-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 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:   PetscBool      opt,pkg;

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

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

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

 83:     Level: advanced

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

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

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


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

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

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

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

124: static PetscErrorCode MatGetDiagonal_MFFD(Mat,Vec);

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

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

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

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

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

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

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

168:    Not Collective

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

174:    Level: developer

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

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

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

189: .keywords: MatMFFD, register

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

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 (ctx->ops->view) {
256:       (*ctx->ops->view)(ctx,viewer);
257:     }
258:     PetscObjectGetOptionsPrefix((PetscObject)J, &prefix);

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

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

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

291:   MatMFFDResetHHistory(J);
292:   j->vshift = 0.0;
293:   j->vscale = 1.0;
294:   return(0);
295: }

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

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

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

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

337:   if (mat->erroriffailure && PetscIsInfOrNanScalar(h)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Computed Nan differencing parameter h");
338:   if (ctx->checkh) {
339:     (*ctx->checkh)(ctx->checkhctx,U,a,&h);
340:   }

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

354:   /* w = u + ha */
355:   if (ctx->drscale) {
356:     VecPointwiseMult(ctx->drscale,a,U);
357:     VecAYPX(U,h,w);
358:   } else {
359:     VecWAXPY(w,h,a,U);
360:   }

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

368:   VecAXPY(y,-1.0,F);
369:   VecScale(y,1.0/h);

371:   if ((ctx->vshift != 0.0) || (ctx->vscale != 1.0)) {
372:     VecAXPBY(y,ctx->vshift,ctx->vscale,a);
373:   }
374:   if (ctx->dlscale) {
375:     VecPointwiseMult(y,ctx->dlscale,y);
376:   }
377:   if (ctx->dshift) {
378:     if (!ctx->dshiftw) {
379:       VecDuplicate(y,&ctx->dshiftw);
380:     }
381:     VecPointwiseMult(ctx->dshift,a,ctx->dshiftw);
382:     VecAXPY(y,1.0,ctx->dshiftw);
383:   }

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

387:   PetscLogEventEnd(MATMFFD_Mult,a,y,0,0);
388:   return(0);
389: }

391: /*
392:   MatGetDiagonal_MFFD - Gets the diagonal for a matrix free matrix

394:         y ~= (F(u + ha) - F(u))/h,
395:   where F = nonlinear function, as set by SNESSetFunction()
396:         u = current iterate
397:         h = difference interval
398: */
399: PetscErrorCode MatGetDiagonal_MFFD(Mat mat,Vec a)
400: {
401:   MatMFFD        ctx = (MatMFFD)mat->data;
402:   PetscScalar    h,*aa,*ww,v;
403:   PetscReal      epsilon = PETSC_SQRT_MACHINE_EPSILON,umin = 100.0*PETSC_SQRT_MACHINE_EPSILON;
404:   Vec            w,U;
406:   PetscInt       i,rstart,rend;

409:   if (!ctx->func) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Requires calling MatMFFDSetFunction() first");
410:   if (!ctx->funci) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Requires calling MatMFFDSetFunctioni() first");
411:   if (!ctx->funcisetbase) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Requires calling MatMFFDSetFunctioniBase() first");
412:   w    = ctx->w;
413:   U    = ctx->current_u;
414:   (*ctx->func)(ctx->funcctx,U,a);
415:   (*ctx->funcisetbase)(ctx->funcctx,U);
416:   VecCopy(U,w);

418:   VecGetOwnershipRange(a,&rstart,&rend);
419:   VecGetArray(a,&aa);
420:   for (i=rstart; i<rend; i++) {
421:     VecGetArray(w,&ww);
422:     h    = ww[i-rstart];
423:     if (h == 0.0) h = 1.0;
424:     if (PetscAbsScalar(h) < umin && PetscRealPart(h) >= 0.0)     h = umin;
425:     else if (PetscRealPart(h) < 0.0 && PetscAbsScalar(h) < umin) h = -umin;
426:     h *= epsilon;

428:     ww[i-rstart] += h;
429:     VecRestoreArray(w,&ww);
430:     (*ctx->funci)(ctx->funcctx,i,w,&v);
431:     aa[i-rstart]  = (v - aa[i-rstart])/h;

433:     /* possibly shift and scale result */
434:     if ((ctx->vshift != 0.0) || (ctx->vscale != 1.0)) {
435:       aa[i - rstart] = ctx->vshift + ctx->vscale*aa[i-rstart];
436:     }

438:     VecGetArray(w,&ww);
439:     ww[i-rstart] -= h;
440:     VecRestoreArray(w,&ww);
441:   }
442:   VecRestoreArray(a,&aa);
443:   return(0);
444: }

446: static PetscErrorCode MatDiagonalScale_MFFD(Mat mat,Vec ll,Vec rr)
447: {
448:   MatMFFD        aij = (MatMFFD)mat->data;

452:   if (ll && !aij->dlscale) {
453:     VecDuplicate(ll,&aij->dlscale);
454:   }
455:   if (rr && !aij->drscale) {
456:     VecDuplicate(rr,&aij->drscale);
457:   }
458:   if (ll) {
459:     VecCopy(ll,aij->dlscale);
460:   }
461:   if (rr) {
462:     VecCopy(rr,aij->drscale);
463:   }
464:   return(0);
465: }

467: static PetscErrorCode MatDiagonalSet_MFFD(Mat mat,Vec ll,InsertMode mode)
468: {
469:   MatMFFD        aij = (MatMFFD)mat->data;

473:   if (mode == INSERT_VALUES) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"No diagonal set with INSERT_VALUES");
474:   if (!aij->dshift) {
475:     VecDuplicate(ll,&aij->dshift);
476:   }
477:   VecAXPY(aij->dshift,1.0,ll);
478:   return(0);
479: }

481: static PetscErrorCode MatShift_MFFD(Mat Y,PetscScalar a)
482: {
483:   MatMFFD shell = (MatMFFD)Y->data;

486:   shell->vshift += a;
487:   return(0);
488: }

490: static PetscErrorCode MatScale_MFFD(Mat Y,PetscScalar a)
491: {
492:   MatMFFD shell = (MatMFFD)Y->data;

495:   shell->vscale *= a;
496:   return(0);
497: }

499: PETSC_EXTERN PetscErrorCode MatMFFDSetBase_MFFD(Mat J,Vec U,Vec F)
500: {
502:   MatMFFD        ctx = (MatMFFD)J->data;

505:   MatMFFDResetHHistory(J);
506:   if (!ctx->current_u) {
507:     VecDuplicate(U,&ctx->current_u);
508:     VecLockPush(ctx->current_u);
509:   }
510:   VecLockPop(ctx->current_u);
511:   VecCopy(U,ctx->current_u);
512:   VecLockPush(ctx->current_u);
513:   if (F) {
514:     if (ctx->current_f_allocated) {VecDestroy(&ctx->current_f);}
515:     ctx->current_f           = F;
516:     ctx->current_f_allocated = PETSC_FALSE;
517:   } else if (!ctx->current_f_allocated) {
518:     MatCreateVecs(J,NULL,&ctx->current_f);

520:     ctx->current_f_allocated = PETSC_TRUE;
521:   }
522:   if (!ctx->w) {
523:     VecDuplicate(ctx->current_u,&ctx->w);
524:   }
525:   J->assembled = PETSC_TRUE;
526:   return(0);
527: }

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

531: static PetscErrorCode  MatMFFDSetCheckh_MFFD(Mat J,FCN3 fun,void *ectx)
532: {
533:   MatMFFD ctx = (MatMFFD)J->data;

536:   ctx->checkh    = fun;
537:   ctx->checkhctx = ectx;
538:   return(0);
539: }

541: /*@C
542:    MatMFFDSetOptionsPrefix - Sets the prefix used for searching for all
543:    MatMFFD options in the database.

545:    Collective on Mat

547:    Input Parameter:
548: +  A - the Mat context
549: -  prefix - the prefix to prepend to all option names

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

555:    Level: advanced

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

559: .seealso: MatSetFromOptions(), MatCreateSNESMF(), MatCreateMFFD()
560: @*/
561: PetscErrorCode  MatMFFDSetOptionsPrefix(Mat mat,const char prefix[])

563: {
564:   MatMFFD        mfctx = mat ? (MatMFFD)mat->data : (MatMFFD)NULL;

570:   PetscObjectSetOptionsPrefix((PetscObject)mfctx,prefix);
571:   return(0);
572: }

574: static PetscErrorCode  MatSetFromOptions_MFFD(PetscOptionItems *PetscOptionsObject,Mat mat)
575: {
576:   MatMFFD        mfctx = (MatMFFD)mat->data;
578:   PetscBool      flg;
579:   char           ftype[256];

584:   PetscObjectOptionsBegin((PetscObject)mfctx);
585:   PetscOptionsFList("-mat_mffd_type","Matrix free type","MatMFFDSetType",MatMFFDList,((PetscObject)mfctx)->type_name,ftype,256,&flg);
586:   if (flg) {
587:     MatMFFDSetType(mat,ftype);
588:   }

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

593:   flg  = PETSC_FALSE;
594:   PetscOptionsBool("-mat_mffd_check_positivity","Insure that U + h*a is nonnegative","MatMFFDSetCheckh",flg,&flg,NULL);
595:   if (flg) {
596:     MatMFFDSetCheckh(mat,MatMFFDCheckPositivity,0);
597:   }
598:   if (mfctx->ops->setfromoptions) {
599:     (*mfctx->ops->setfromoptions)(PetscOptionsObject,mfctx);
600:   }
601:   PetscOptionsEnd();
602:   return(0);
603: }

605: static PetscErrorCode  MatMFFDSetPeriod_MFFD(Mat mat,PetscInt period)
606: {
607:   MatMFFD ctx = (MatMFFD)mat->data;

610:   ctx->recomputeperiod = period;
611:   return(0);
612: }

614: static PetscErrorCode  MatMFFDSetFunction_MFFD(Mat mat,PetscErrorCode (*func)(void*,Vec,Vec),void *funcctx)
615: {
616:   MatMFFD ctx = (MatMFFD)mat->data;

619:   ctx->func    = func;
620:   ctx->funcctx = funcctx;
621:   return(0);
622: }

624: static PetscErrorCode  MatMFFDSetFunctionError_MFFD(Mat mat,PetscReal error)
625: {
626:   MatMFFD ctx = (MatMFFD)mat->data;

629:   if (error != PETSC_DEFAULT) ctx->error_rel = error;
630:   return(0);
631: }

633: static PetscErrorCode MatMissingDiagonal_MFFD(Mat A,PetscBool  *missing,PetscInt *d)
634: {
636:   *missing = PETSC_FALSE;
637:   return(0);
638: }

640: /*MC
641:   MATMFFD - MATMFFD = "mffd" - A matrix free matrix type.

643:   Level: advanced

645: .seealso: MatCreateMFFD(), MatCreateSNESMF(), MatMFFDSetFunction(), MatMFFDSetType(),  
646:           MatMFFDSetFunctionError(), MatMFFDDSSetUmin(), MatMFFDSetFunction()
647:           MatMFFDSetHHistory(), MatMFFDResetHHistory(), MatCreateSNESMF(),
648:           MatMFFDGetH(),
649: M*/
650: PETSC_EXTERN PetscErrorCode MatCreate_MFFD(Mat A)
651: {
652:   MatMFFD        mfctx;

656:   MatMFFDInitializePackage();

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

660:   mfctx->error_rel                = PETSC_SQRT_MACHINE_EPSILON;
661:   mfctx->recomputeperiod          = 1;
662:   mfctx->count                    = 0;
663:   mfctx->currenth                 = 0.0;
664:   mfctx->historyh                 = NULL;
665:   mfctx->ncurrenth                = 0;
666:   mfctx->maxcurrenth              = 0;
667:   ((PetscObject)mfctx)->type_name = 0;

669:   mfctx->vshift = 0.0;
670:   mfctx->vscale = 1.0;

672:   /*
673:      Create the empty data structure to contain compute-h routines.
674:      These will be filled in below from the command line options or
675:      a later call with MatMFFDSetType() or if that is not called
676:      then it will default in the first use of MatMult_MFFD()
677:   */
678:   mfctx->ops->compute        = 0;
679:   mfctx->ops->destroy        = 0;
680:   mfctx->ops->view           = 0;
681:   mfctx->ops->setfromoptions = 0;
682:   mfctx->hctx                = 0;

684:   mfctx->func    = 0;
685:   mfctx->funcctx = 0;
686:   mfctx->w       = NULL;

688:   A->data = mfctx;

690:   A->ops->mult            = MatMult_MFFD;
691:   A->ops->destroy         = MatDestroy_MFFD;
692:   A->ops->view            = MatView_MFFD;
693:   A->ops->assemblyend     = MatAssemblyEnd_MFFD;
694:   A->ops->scale           = MatScale_MFFD;
695:   A->ops->shift           = MatShift_MFFD;
696:   A->ops->diagonalscale   = MatDiagonalScale_MFFD;
697:   A->ops->diagonalset     = MatDiagonalSet_MFFD;
698:   A->ops->setfromoptions  = MatSetFromOptions_MFFD;
699:   A->ops->missingdiagonal = MatMissingDiagonal_MFFD;
700:   A->assembled            = PETSC_TRUE;

702:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetBase_C",MatMFFDSetBase_MFFD);
703:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctioniBase_C",MatMFFDSetFunctioniBase_MFFD);
704:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctioni_C",MatMFFDSetFunctioni_MFFD);
705:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunction_C",MatMFFDSetFunction_MFFD);
706:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetCheckh_C",MatMFFDSetCheckh_MFFD);
707:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetPeriod_C",MatMFFDSetPeriod_MFFD);
708:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctionError_C",MatMFFDSetFunctionError_MFFD);
709:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDResetHHistory_C",MatMFFDResetHHistory_MFFD);

711:   mfctx->mat = A;

713:   PetscObjectChangeTypeName((PetscObject)A,MATMFFD);
714:   return(0);
715: }

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

720:    Collective on Vec

722:    Input Parameters:
723: +  comm - MPI communicator
724: .  m - number of local rows (or PETSC_DECIDE to have calculated if M is given)
725:            This value should be the same as the local size used in creating the
726:            y vector for the matrix-vector product y = Ax.
727: .  n - This value should be the same as the local size used in creating the
728:        x vector for the matrix-vector product y = Ax. (or PETSC_DECIDE to have
729:        calculated if N is given) For square matrices n is almost always m.
730: .  M - number of global rows (or PETSC_DETERMINE to have calculated if m is given)
731: -  N - number of global columns (or PETSC_DETERMINE to have calculated if n is given)


734:    Output Parameter:
735: .  J - the matrix-free matrix

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


743:    Level: advanced

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

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

752: .vb
753:      F'(u)*a = [F(u+h*a) - F(u)]/h where
754:      h = error_rel*u'a/||a||^2                        if  |u'a| > umin*||a||_{1}
755:        = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2   otherwise
756:  where
757:      error_rel = square root of relative error in function evaluation
758:      umin = minimum iterate parameter
759: .ve

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

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

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

770:    Options Database Keys:
771: +  -mat_mffd_err <error_rel> - Sets error_rel
772: .  -mat_mffd_unim <umin> - Sets umin (for default PETSc routine that computes h only)
773: -  -mat_mffd_check_positivity

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

777: .seealso: MatDestroy(), MatMFFDSetFunctionError(), MatMFFDDSSetUmin(), MatMFFDSetFunction()
778:           MatMFFDSetHHistory(), MatMFFDResetHHistory(), MatCreateSNESMF(),
779:           MatMFFDGetH(), MatMFFDRegister(), MatMFFDComputeJacobian()

781: @*/
782: PetscErrorCode  MatCreateMFFD(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,Mat *J)
783: {

787:   MatCreate(comm,J);
788:   MatSetSizes(*J,m,n,M,N);
789:   MatSetType(*J,MATMFFD);
790:   MatSetUp(*J);
791:   return(0);
792: }

794: /*@
795:    MatMFFDGetH - Gets the last value that was used as the differencing
796:    parameter.

798:    Not Collective

800:    Input Parameters:
801: .  mat - the matrix obtained with MatCreateSNESMF()

803:    Output Paramter:
804: .  h - the differencing step size

806:    Level: advanced

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

810: .seealso: MatCreateSNESMF(),MatMFFDSetHHistory(), MatCreateMFFD(), MATMFFD, MatMFFDResetHHistory()
811: @*/
812: PetscErrorCode  MatMFFDGetH(Mat mat,PetscScalar *h)
813: {
814:   MatMFFD        ctx = (MatMFFD)mat->data;
816:   PetscBool      match;

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

824:   *h = ctx->currenth;
825:   return(0);
826: }

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

831:    Logically Collective on Mat

833:    Input Parameters:
834: +  mat - the matrix free matrix created via MatCreateSNESMF() or MatCreateMFFD()
835: .  func - the function to use
836: -  funcctx - optional function context passed to function

838:    Calling Sequence of func:
839: $     func (void *funcctx, Vec x, Vec f)

841: +  funcctx - user provided context
842: .  x - input vector
843: -  f - computed output function

845:    Level: advanced

847:    Notes:
848:     If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
849:     matrix inside your compute Jacobian routine

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

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

855: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD,
856:           MatMFFDSetHHistory(), MatMFFDResetHHistory(), SNESetFunction()
857: @*/
858: PetscErrorCode  MatMFFDSetFunction(Mat mat,PetscErrorCode (*func)(void*,Vec,Vec),void *funcctx)
859: {

864:   PetscTryMethod(mat,"MatMFFDSetFunction_C",(Mat,PetscErrorCode (*)(void*,Vec,Vec),void*),(mat,func,funcctx));
865:   return(0);
866: }

868: /*@C
869:    MatMFFDSetFunctioni - Sets the function for a single component

871:    Logically Collective on Mat

873:    Input Parameters:
874: +  mat - the matrix free matrix created via MatCreateSNESMF()
875: -  funci - the function to use

877:    Level: advanced

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

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

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

889: @*/
890: PetscErrorCode  MatMFFDSetFunctioni(Mat mat,PetscErrorCode (*funci)(void*,PetscInt,Vec,PetscScalar*))
891: {

896:   PetscTryMethod(mat,"MatMFFDSetFunctioni_C",(Mat,PetscErrorCode (*)(void*,PetscInt,Vec,PetscScalar*)),(mat,funci));
897:   return(0);
898: }

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

903:    Logically Collective on Mat

905:    Input Parameters:
906: +  mat - the matrix free matrix created via MatCreateSNESMF()
907: -  func - the function to use

909:    Level: advanced

911:    Notes:
912:     If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
913:     matrix inside your compute Jacobian routine.
914:     This function is necessary to compute the diagonal of the matrix.


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

919: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD
920:           MatMFFDSetHHistory(), MatMFFDResetHHistory(), SNESetFunction(), MatGetDiagonal()
921: @*/
922: PetscErrorCode  MatMFFDSetFunctioniBase(Mat mat,PetscErrorCode (*func)(void*,Vec))
923: {

928:   PetscTryMethod(mat,"MatMFFDSetFunctioniBase_C",(Mat,PetscErrorCode (*)(void*,Vec)),(mat,func));
929:   return(0);
930: }

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

935:    Logically Collective on Mat

937:    Input Parameters:
938: +  mat - the matrix free matrix created via MatCreateSNESMF()
939: -  period - 1 for everytime, 2 for every second etc

941:    Options Database Keys:
942: +  -mat_mffd_period <period>

944:    Level: advanced


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

949: .seealso: MatCreateSNESMF(),MatMFFDGetH(),
950:           MatMFFDSetHHistory(), MatMFFDResetHHistory()
951: @*/
952: PetscErrorCode  MatMFFDSetPeriod(Mat mat,PetscInt period)
953: {

959:   PetscTryMethod(mat,"MatMFFDSetPeriod_C",(Mat,PetscInt),(mat,period));
960:   return(0);
961: }

963: /*@
964:    MatMFFDSetFunctionError - Sets the error_rel for the approximation of
965:    matrix-vector products using finite differences.

967:    Logically Collective on Mat

969:    Input Parameters:
970: +  mat - the matrix free matrix created via MatCreateMFFD() or MatCreateSNESMF()
971: -  error_rel - relative error (should be set to the square root of
972:                the relative error in the function evaluations)

974:    Options Database Keys:
975: +  -mat_mffd_err <error_rel> - Sets error_rel

977:    Level: advanced

979:    Notes:
980:    The default matrix-free matrix-vector product routine computes
981: .vb
982:      F'(u)*a = [F(u+h*a) - F(u)]/h where
983:      h = error_rel*u'a/||a||^2                        if  |u'a| > umin*||a||_{1}
984:        = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2   else
985: .ve

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

989: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD
990:           MatMFFDSetHHistory(), MatMFFDResetHHistory()
991: @*/
992: PetscErrorCode  MatMFFDSetFunctionError(Mat mat,PetscReal error)
993: {

999:   PetscTryMethod(mat,"MatMFFDSetFunctionError_C",(Mat,PetscReal),(mat,error));
1000:   return(0);
1001: }

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

1007:    Logically Collective on Mat

1009:    Input Parameters:
1010: +  J - the matrix-free matrix context
1011: .  histroy - space to hold the history
1012: -  nhistory - number of entries in history, if more entries are generated than
1013:               nhistory, then the later ones are discarded

1015:    Level: advanced

1017:    Notes:
1018:    Use MatMFFDResetHHistory() to reset the history counter and collect
1019:    a new batch of differencing parameters, h.

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

1023: .seealso: MatMFFDGetH(), MatCreateSNESMF(),
1024:           MatMFFDResetHHistory(), MatMFFDSetFunctionError()

1026: @*/
1027: PetscErrorCode  MatMFFDSetHHistory(Mat J,PetscScalar history[],PetscInt nhistory)
1028: {
1029:   MatMFFD        ctx = (MatMFFD)J->data;
1031:   PetscBool      match;

1037:   PetscObjectTypeCompare((PetscObject)J,MATMFFD,&match);
1038:   if (!match) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_WRONG,"Not a MFFD matrix");
1039:   ctx->historyh    = history;
1040:   ctx->maxcurrenth = nhistory;
1041:   ctx->currenth    = 0.;
1042:   return(0);
1043: }

1045: /*@
1046:    MatMFFDResetHHistory - Resets the counter to zero to begin
1047:    collecting a new set of differencing histories.

1049:    Logically Collective on Mat

1051:    Input Parameters:
1052: .  J - the matrix-free matrix context

1054:    Level: advanced

1056:    Notes:
1057:    Use MatMFFDSetHHistory() to create the original history counter.

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

1061: .seealso: MatMFFDGetH(), MatCreateSNESMF(),
1062:           MatMFFDSetHHistory(), MatMFFDSetFunctionError()

1064: @*/
1065: PetscErrorCode  MatMFFDResetHHistory(Mat J)
1066: {

1071:   PetscTryMethod(J,"MatMFFDResetHHistory_C",(Mat),(J));
1072:   return(0);
1073: }

1075: /*@
1076:     MatMFFDSetBase - Sets the vector U at which matrix vector products of the
1077:         Jacobian are computed

1079:     Logically Collective on Mat

1081:     Input Parameters:
1082: +   J - the MatMFFD matrix
1083: .   U - the vector
1084: -   F - (optional) vector that contains F(u) if it has been already computed

1086:     Notes: This is rarely used directly

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

1091:     Level: advanced

1093: @*/
1094: PetscErrorCode  MatMFFDSetBase(Mat J,Vec U,Vec F)
1095: {

1102:   PetscTryMethod(J,"MatMFFDSetBase_C",(Mat,Vec,Vec),(J,U,F));
1103:   return(0);
1104: }

1106: /*@C
1107:     MatMFFDSetCheckh - Sets a function that checks the computed h and adjusts
1108:         it to satisfy some criteria

1110:     Logically Collective on Mat

1112:     Input Parameters:
1113: +   J - the MatMFFD matrix
1114: .   fun - the function that checks h
1115: -   ctx - any context needed by the function

1117:     Options Database Keys:
1118: .   -mat_mffd_check_positivity

1120:     Level: advanced

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

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

1128: .seealso:  MatMFFDCheckPositivity()
1129: @*/
1130: PetscErrorCode  MatMFFDSetCheckh(Mat J,PetscErrorCode (*fun)(void*,Vec,Vec,PetscScalar*),void *ctx)
1131: {

1136:   PetscTryMethod(J,"MatMFFDSetCheckh_C",(Mat,PetscErrorCode (*)(void*,Vec,Vec,PetscScalar*),void*),(J,fun,ctx));
1137:   return(0);
1138: }

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

1144:     Logically Collective on Vec

1146:     Input Parameters:
1147: +   U - base vector that is added to
1148: .   a - vector that is added
1149: .   h - scaling factor on a
1150: -   dummy - context variable (unused)

1152:     Options Database Keys:
1153: .   -mat_mffd_check_positivity

1155:     Level: advanced

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

1160: .seealso:  MatMFFDSetCheckh()
1161: @*/
1162: PetscErrorCode  MatMFFDCheckPositivity(void *dummy,Vec U,Vec a,PetscScalar *h)
1163: {
1164:   PetscReal      val, minval;
1165:   PetscScalar    *u_vec, *a_vec;
1167:   PetscInt       i,n;
1168:   MPI_Comm       comm;

1174:   PetscObjectGetComm((PetscObject)U,&comm);
1175:   VecGetArray(U,&u_vec);
1176:   VecGetArray(a,&a_vec);
1177:   VecGetLocalSize(U,&n);
1178:   minval = PetscAbsScalar(*h)*PetscRealConstant(1.01);
1179:   for (i=0; i<n; i++) {
1180:     if (PetscRealPart(u_vec[i] + *h*a_vec[i]) <= 0.0) {
1181:       val = PetscAbsScalar(u_vec[i]/a_vec[i]);
1182:       if (val < minval) minval = val;
1183:     }
1184:   }
1185:   VecRestoreArray(U,&u_vec);
1186:   VecRestoreArray(a,&a_vec);
1187:   MPIU_Allreduce(&minval,&val,1,MPIU_REAL,MPIU_MIN,comm);
1188:   if (val <= PetscAbsScalar(*h)) {
1189:     PetscInfo2(U,"Scaling back h from %g to %g\n",(double)PetscRealPart(*h),(double)(.99*val));
1190:     if (PetscRealPart(*h) > 0.0) *h =  0.99*val;
1191:     else                         *h = -0.99*val;
1192:   }
1193:   return(0);
1194: }