Actual source code: mffd.c

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

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

  5: PetscFunctionList MatMFFDList              = 0;
  6: PetscBool         MatMFFDRegisterAllCalled = PETSC_FALSE;

  8: PetscClassId  MATMFFD_CLASSID;
  9: PetscLogEvent MATMFFD_Mult;

 11: static PetscBool MatMFFDPackageInitialized = PETSC_FALSE;
 12: /*@C
 13:   MatMFFDFinalizePackage - This function destroys everything in the MatMFFD package. It is
 14:   called from PetscFinalize().

 16:   Level: developer

 18: .seealso: PetscFinalize(), MatCreateMFFD(), MatCreateSNESMF()
 19: @*/
 20: PetscErrorCode  MatMFFDFinalizePackage(void)
 21: {

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

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

 35:   Level: developer

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

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

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

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

 80:     Level: advanced

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

 86:                         F(x+ha) - F(x)
 87:           F'(u)a  ~=  ----------------
 88:                               h

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


102:   PetscObjectTypeCompare((PetscObject)mat,MATMFFD,&match);
103:   if (!match) return(0);

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

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

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

121: static PetscErrorCode MatGetDiagonal_MFFD(Mat,Vec);

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

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

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

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

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

158:   ctx->ncurrenth = 0;
159:   return(0);
160: }

162: /*@C
163:    MatMFFDRegister - Adds a method to the MatMFFD registry.

165:    Not Collective

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

171:    Level: developer

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

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

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

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

193:   MatInitializePackage();
194:   PetscFunctionListAdd(&MatMFFDList,sname,function);
195:   return(0);
196: }

198: /* ----------------------------------------------------------------------------------------*/
199: static PetscErrorCode MatDestroy_MFFD(Mat mat)
200: {
202:   MatMFFD        ctx = (MatMFFD)mat->data;

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

218:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetBase_C",NULL);
219:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctioniBase_C",NULL);
220:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctioni_C",NULL);
221:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunction_C",NULL);
222:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctionError_C",NULL);
223:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetCheckh_C",NULL);
224:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetPeriod_C",NULL);
225:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDResetHHistory_C",NULL);
226:   return(0);
227: }

229: /*
230:    MatMFFDView_MFFD - Views matrix-free parameters.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

385:   /* This "if" prevents PETSc from erroring when the mat is rectangular */
386:   if ((ctx->vshift != 0.0) || (ctx->vscale != 1.0)) {
387:     VecAXPBY(y,ctx->vshift,ctx->vscale,a);
388:   }

390:   if (ctx->dlscale) {
391:     VecPointwiseMult(y,ctx->dlscale,y);
392:   }
393:   if (ctx->dshift) {
394:     if (!ctx->dshiftw) {
395:       VecDuplicate(y,&ctx->dshiftw);
396:     }
397:     VecPointwiseMult(ctx->dshift,a,ctx->dshiftw);
398:     VecAXPY(y,1.0,ctx->dshiftw);
399:   }

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

403:   PetscLogEventEnd(MATMFFD_Mult,a,y,0,0);
404:   return(0);
405: }

407: /*
408:   MatGetDiagonal_MFFD - Gets the diagonal for a matrix free matrix

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

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

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

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

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

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

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

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

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

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

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

502:   shell->vshift += a;
503:   return(0);
504: }

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

511:   shell->vscale *= a;
512:   return(0);
513: }

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

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

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

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

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

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

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

561:    Collective on Mat

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

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

571:    Level: advanced

573: .seealso: MatSetFromOptions(), MatCreateSNESMF(), MatCreateMFFD()
574: @*/
575: PetscErrorCode  MatMFFDSetOptionsPrefix(Mat mat,const char prefix[])

577: {
578:   MatMFFD        mfctx = mat ? (MatMFFD)mat->data : (MatMFFD)NULL;

584:   PetscObjectSetOptionsPrefix((PetscObject)mfctx,prefix);
585:   return(0);
586: }

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

598:   PetscObjectOptionsBegin((PetscObject)mfctx);
599:   PetscOptionsFList("-mat_mffd_type","Matrix free type","MatMFFDSetType",MatMFFDList,((PetscObject)mfctx)->type_name,ftype,256,&flg);
600:   if (flg) {
601:     MatMFFDSetType(mat,ftype);
602:   }

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

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

622: static PetscErrorCode  MatMFFDSetPeriod_MFFD(Mat mat,PetscInt period)
623: {
624:   MatMFFD ctx = (MatMFFD)mat->data;

627:   ctx->recomputeperiod = period;
628:   return(0);
629: }

631: static PetscErrorCode  MatMFFDSetFunction_MFFD(Mat mat,PetscErrorCode (*func)(void*,Vec,Vec),void *funcctx)
632: {
633:   MatMFFD ctx = (MatMFFD)mat->data;

636:   ctx->func    = func;
637:   ctx->funcctx = funcctx;
638:   return(0);
639: }

641: static PetscErrorCode  MatMFFDSetFunctionError_MFFD(Mat mat,PetscReal error)
642: {
643:   MatMFFD ctx = (MatMFFD)mat->data;

646:   if (error != PETSC_DEFAULT) ctx->error_rel = error;
647:   return(0);
648: }

650: static PetscErrorCode MatMissingDiagonal_MFFD(Mat A,PetscBool  *missing,PetscInt *d)
651: {
653:   *missing = PETSC_FALSE;
654:   return(0);
655: }

657: /*MC
658:   MATMFFD - MATMFFD = "mffd" - A matrix free matrix type.

660:   Level: advanced

662: .seealso: MatCreateMFFD(), MatCreateSNESMF(), MatMFFDSetFunction(), MatMFFDSetType(),  
663:           MatMFFDSetFunctionError(), MatMFFDDSSetUmin(), MatMFFDSetFunction()
664:           MatMFFDSetHHistory(), MatMFFDResetHHistory(), MatCreateSNESMF(),
665:           MatMFFDGetH(),
666: M*/
667: PETSC_EXTERN PetscErrorCode MatCreate_MFFD(Mat A)
668: {
669:   MatMFFD        mfctx;

673:   MatMFFDInitializePackage();

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

677:   mfctx->error_rel                = PETSC_SQRT_MACHINE_EPSILON;
678:   mfctx->recomputeperiod          = 1;
679:   mfctx->count                    = 0;
680:   mfctx->currenth                 = 0.0;
681:   mfctx->historyh                 = NULL;
682:   mfctx->ncurrenth                = 0;
683:   mfctx->maxcurrenth              = 0;
684:   ((PetscObject)mfctx)->type_name = 0;

686:   mfctx->vshift = 0.0;
687:   mfctx->vscale = 1.0;

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

701:   mfctx->func    = 0;
702:   mfctx->funcctx = 0;
703:   mfctx->w       = NULL;

705:   A->data = mfctx;

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

719:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetBase_C",MatMFFDSetBase_MFFD);
720:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctioniBase_C",MatMFFDSetFunctioniBase_MFFD);
721:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctioni_C",MatMFFDSetFunctioni_MFFD);
722:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunction_C",MatMFFDSetFunction_MFFD);
723:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetCheckh_C",MatMFFDSetCheckh_MFFD);
724:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetPeriod_C",MatMFFDSetPeriod_MFFD);
725:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctionError_C",MatMFFDSetFunctionError_MFFD);
726:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDResetHHistory_C",MatMFFDResetHHistory_MFFD);

728:   mfctx->mat = A;

730:   PetscObjectChangeTypeName((PetscObject)A,MATMFFD);
731:   return(0);
732: }

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

737:    Collective on Vec

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


751:    Output Parameter:
752: .  J - the matrix-free matrix

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


763:    Level: advanced

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

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

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

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

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

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

790:    Options Database Keys:
791: +  -mat_mffd_err <error_rel> - Sets error_rel
792: .  -mat_mffd_unim <umin> - Sets umin (for default PETSc routine that computes h only)
793: -  -mat_mffd_check_positivity

795: .seealso: MatDestroy(), MatMFFDSetFunctionError(), MatMFFDDSSetUmin(), MatMFFDSetFunction()
796:           MatMFFDSetHHistory(), MatMFFDResetHHistory(), MatCreateSNESMF(),
797:           MatMFFDGetH(), MatMFFDRegister(), MatMFFDComputeJacobian()

799: @*/
800: PetscErrorCode  MatCreateMFFD(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,Mat *J)
801: {

805:   MatCreate(comm,J);
806:   MatSetSizes(*J,m,n,M,N);
807:   MatSetType(*J,MATMFFD);
808:   MatSetUp(*J);
809:   return(0);
810: }

812: /*@
813:    MatMFFDGetH - Gets the last value that was used as the differencing
814:    parameter.

816:    Not Collective

818:    Input Parameters:
819: .  mat - the matrix obtained with MatCreateSNESMF()

821:    Output Paramter:
822: .  h - the differencing step size

824:    Level: advanced

826: .seealso: MatCreateSNESMF(),MatMFFDSetHHistory(), MatCreateMFFD(), MATMFFD, MatMFFDResetHHistory()
827: @*/
828: PetscErrorCode  MatMFFDGetH(Mat mat,PetscScalar *h)
829: {
830:   MatMFFD        ctx = (MatMFFD)mat->data;
832:   PetscBool      match;

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

840:   *h = ctx->currenth;
841:   return(0);
842: }

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

847:    Logically Collective on Mat

849:    Input Parameters:
850: +  mat - the matrix free matrix created via MatCreateSNESMF() or MatCreateMFFD()
851: .  func - the function to use
852: -  funcctx - optional function context passed to function

854:    Calling Sequence of func:
855: $     func (void *funcctx, Vec x, Vec f)

857: +  funcctx - user provided context
858: .  x - input vector
859: -  f - computed output function

861:    Level: advanced

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

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

869: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD,
870:           MatMFFDSetHHistory(), MatMFFDResetHHistory(), SNESetFunction()
871: @*/
872: PetscErrorCode  MatMFFDSetFunction(Mat mat,PetscErrorCode (*func)(void*,Vec,Vec),void *funcctx)
873: {

878:   PetscTryMethod(mat,"MatMFFDSetFunction_C",(Mat,PetscErrorCode (*)(void*,Vec,Vec),void*),(mat,func,funcctx));
879:   return(0);
880: }

882: /*@C
883:    MatMFFDSetFunctioni - Sets the function for a single component

885:    Logically Collective on Mat

887:    Input Parameters:
888: +  mat - the matrix free matrix created via MatCreateSNESMF()
889: -  funci - the function to use

891:    Level: advanced

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

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

901: @*/
902: PetscErrorCode  MatMFFDSetFunctioni(Mat mat,PetscErrorCode (*funci)(void*,PetscInt,Vec,PetscScalar*))
903: {

908:   PetscTryMethod(mat,"MatMFFDSetFunctioni_C",(Mat,PetscErrorCode (*)(void*,PetscInt,Vec,PetscScalar*)),(mat,funci));
909:   return(0);
910: }

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

915:    Logically Collective on Mat

917:    Input Parameters:
918: +  mat - the matrix free matrix created via MatCreateSNESMF()
919: -  func - the function to use

921:    Level: advanced

923:    Notes:
924:     If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
925:     matrix inside your compute Jacobian routine.
926:     This function is necessary to compute the diagonal of the matrix.


929: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD
930:           MatMFFDSetHHistory(), MatMFFDResetHHistory(), SNESetFunction(), MatGetDiagonal()
931: @*/
932: PetscErrorCode  MatMFFDSetFunctioniBase(Mat mat,PetscErrorCode (*func)(void*,Vec))
933: {

938:   PetscTryMethod(mat,"MatMFFDSetFunctioniBase_C",(Mat,PetscErrorCode (*)(void*,Vec)),(mat,func));
939:   return(0);
940: }

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

945:    Logically Collective on Mat

947:    Input Parameters:
948: +  mat - the matrix free matrix created via MatCreateSNESMF()
949: -  period - 1 for everytime, 2 for every second etc

951:    Options Database Keys:
952: .  -mat_mffd_period <period>

954:    Level: advanced


957: .seealso: MatCreateSNESMF(),MatMFFDGetH(),
958:           MatMFFDSetHHistory(), MatMFFDResetHHistory()
959: @*/
960: PetscErrorCode  MatMFFDSetPeriod(Mat mat,PetscInt period)
961: {

967:   PetscTryMethod(mat,"MatMFFDSetPeriod_C",(Mat,PetscInt),(mat,period));
968:   return(0);
969: }

971: /*@
972:    MatMFFDSetFunctionError - Sets the error_rel for the approximation of
973:    matrix-vector products using finite differences.

975:    Logically Collective on Mat

977:    Input Parameters:
978: +  mat - the matrix free matrix created via MatCreateMFFD() or MatCreateSNESMF()
979: -  error_rel - relative error (should be set to the square root of
980:                the relative error in the function evaluations)

982:    Options Database Keys:
983: .  -mat_mffd_err <error_rel> - Sets error_rel

985:    Level: advanced

987:    Notes:
988:    The default matrix-free matrix-vector product routine computes
989: .vb
990:      F'(u)*a = [F(u+h*a) - F(u)]/h where
991:      h = error_rel*u'a/||a||^2                        if  |u'a| > umin*||a||_{1}
992:        = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2   else
993: .ve

995: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD
996:           MatMFFDSetHHistory(), MatMFFDResetHHistory()
997: @*/
998: PetscErrorCode  MatMFFDSetFunctionError(Mat mat,PetscReal error)
999: {

1005:   PetscTryMethod(mat,"MatMFFDSetFunctionError_C",(Mat,PetscReal),(mat,error));
1006:   return(0);
1007: }

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

1013:    Logically Collective on Mat

1015:    Input Parameters:
1016: +  J - the matrix-free matrix context
1017: .  histroy - space to hold the history
1018: -  nhistory - number of entries in history, if more entries are generated than
1019:               nhistory, then the later ones are discarded

1021:    Level: advanced

1023:    Notes:
1024:    Use MatMFFDResetHHistory() to reset the history counter and collect
1025:    a new batch of differencing parameters, h.

1027: .seealso: MatMFFDGetH(), MatCreateSNESMF(),
1028:           MatMFFDResetHHistory(), MatMFFDSetFunctionError()

1030: @*/
1031: PetscErrorCode  MatMFFDSetHHistory(Mat J,PetscScalar history[],PetscInt nhistory)
1032: {
1033:   MatMFFD        ctx = (MatMFFD)J->data;
1035:   PetscBool      match;

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

1049: /*@
1050:    MatMFFDResetHHistory - Resets the counter to zero to begin
1051:    collecting a new set of differencing histories.

1053:    Logically Collective on Mat

1055:    Input Parameters:
1056: .  J - the matrix-free matrix context

1058:    Level: advanced

1060:    Notes:
1061:    Use MatMFFDSetHHistory() to create the original history counter.

1063: .seealso: MatMFFDGetH(), MatCreateSNESMF(),
1064:           MatMFFDSetHHistory(), MatMFFDSetFunctionError()

1066: @*/
1067: PetscErrorCode  MatMFFDResetHHistory(Mat J)
1068: {

1073:   PetscTryMethod(J,"MatMFFDResetHHistory_C",(Mat),(J));
1074:   return(0);
1075: }

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

1081:     Logically Collective on Mat

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

1088:     Notes:
1089:     This is rarely used directly

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

1094:     Level: advanced

1096: @*/
1097: PetscErrorCode  MatMFFDSetBase(Mat J,Vec U,Vec F)
1098: {

1105:   PetscTryMethod(J,"MatMFFDSetBase_C",(Mat,Vec,Vec),(J,U,F));
1106:   return(0);
1107: }

1109: /*@C
1110:     MatMFFDSetCheckh - Sets a function that checks the computed h and adjusts
1111:         it to satisfy some criteria

1113:     Logically Collective on Mat

1115:     Input Parameters:
1116: +   J - the MatMFFD matrix
1117: .   fun - the function that checks h
1118: -   ctx - any context needed by the function

1120:     Options Database Keys:
1121: .   -mat_mffd_check_positivity

1123:     Level: advanced

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

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

1132: .seealso:  MatMFFDCheckPositivity()
1133: @*/
1134: PetscErrorCode  MatMFFDSetCheckh(Mat J,PetscErrorCode (*fun)(void*,Vec,Vec,PetscScalar*),void *ctx)
1135: {

1140:   PetscTryMethod(J,"MatMFFDSetCheckh_C",(Mat,PetscErrorCode (*)(void*,Vec,Vec,PetscScalar*),void*),(J,fun,ctx));
1141:   return(0);
1142: }

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

1148:     Logically Collective on Vec

1150:     Input Parameters:
1151: +   U - base vector that is added to
1152: .   a - vector that is added
1153: .   h - scaling factor on a
1154: -   dummy - context variable (unused)

1156:     Options Database Keys:
1157: .   -mat_mffd_check_positivity

1159:     Level: advanced

1161:     Notes:
1162:     This is rarely used directly, rather it is passed as an argument to
1163:            MatMFFDSetCheckh()

1165: .seealso:  MatMFFDSetCheckh()
1166: @*/
1167: PetscErrorCode  MatMFFDCheckPositivity(void *dummy,Vec U,Vec a,PetscScalar *h)
1168: {
1169:   PetscReal      val, minval;
1170:   PetscScalar    *u_vec, *a_vec;
1172:   PetscInt       i,n;
1173:   MPI_Comm       comm;

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