Actual source code: mffd.c

petsc-3.5.0 2014-06-30
Report Typos and Errors
  2: #include <petsc-private/matimpl.h>
  3: #include <../src/mat/impls/mffd/mffdimpl.h>   /*I  "petscmat.h"   I*/

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

  8: PetscClassId  MATMFFD_CLASSID;
  9: PetscLogEvent MATMFFD_Mult;

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

 18:   Level: developer

 20: .keywords: Petsc, destroy, package
 21: .seealso: PetscFinalize()
 22: @*/
 23: PetscErrorCode  MatMFFDFinalizePackage(void)
 24: {

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

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

 41:   Level: developer

 43: .keywords: Vec, initialize, package
 44: .seealso: PetscInitialize()
 45: @*/
 46: PetscErrorCode  MatMFFDInitializePackage(void)
 47: {
 48:   char           logList[256];
 49:   char           *className;
 50:   PetscBool      opt;

 54:   if (MatMFFDPackageInitialized) return(0);
 55:   MatMFFDPackageInitialized = PETSC_TRUE;
 56:   /* Register Classes */
 57:   PetscClassIdRegister("MatMFFD",&MATMFFD_CLASSID);
 58:   /* Register Constructors */
 59:   MatMFFDRegisterAll();
 60:   /* Register Events */
 61:   PetscLogEventRegister("MatMult MF",          MATMFFD_CLASSID,&MATMFFD_Mult);

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

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

 89:     Input Parameters:
 90: +   mat - the "matrix-free" matrix created via MatCreateSNESMF(), or MatCreateMFFD()
 91:           or MatSetType(mat,MATMFFD);
 92: -   ftype - the type requested, either MATMFFD_WP or MATMFFD_DS

 94:     Level: advanced

 96:     Notes:
 97:     For example, such routines can compute h for use in
 98:     Jacobian-vector products of the form

100:                         F(x+ha) - F(x)
101:           F'(u)a  ~=  ----------------
102:                               h

104: .seealso: MatCreateSNESMF(), MatMFFDRegister(), MatMFFDSetFunction()
105: @*/
106: PetscErrorCode  MatMFFDSetType(Mat mat,MatMFFDType ftype)
107: {
108:   PetscErrorCode ierr,(*r)(MatMFFD);
109:   MatMFFD        ctx = (MatMFFD)mat->data;
110:   PetscBool      match;


116:   PetscObjectTypeCompare((PetscObject)mat,MATMFFD,&match);
117:   if (!match) return(0);

119:   /* already set, so just return */
120:   PetscObjectTypeCompare((PetscObject)ctx,ftype,&match);
121:   if (match) return(0);

123:   /* destroy the old one if it exists */
124:   if (ctx->ops->destroy) {
125:     (*ctx->ops->destroy)(ctx);
126:   }

128:    PetscFunctionListFind(MatMFFDList,ftype,&r);
129:   if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown MatMFFD type %s given",ftype);
130:   (*r)(ctx);
131:   PetscObjectChangeTypeName((PetscObject)ctx,ftype);
132:   return(0);
133: }

135: typedef PetscErrorCode (*FCN1)(void*,Vec); /* force argument to next function to not be extern C*/
138: PetscErrorCode  MatMFFDSetFunctioniBase_MFFD(Mat mat,FCN1 func)
139: {
140:   MatMFFD ctx = (MatMFFD)mat->data;

143:   ctx->funcisetbase = func;
144:   return(0);
145: }

147: typedef PetscErrorCode (*FCN2)(void*,PetscInt,Vec,PetscScalar*); /* force argument to next function to not be extern C*/
150: PetscErrorCode  MatMFFDSetFunctioni_MFFD(Mat mat,FCN2 funci)
151: {
152:   MatMFFD ctx = (MatMFFD)mat->data;

155:   ctx->funci = funci;
156:   return(0);
157: }

161: PetscErrorCode  MatMFFDResetHHistory_MFFD(Mat J)
162: {
163:   MatMFFD ctx = (MatMFFD)J->data;

166:   ctx->ncurrenth = 0;
167:   return(0);
168: }

172: /*@C
173:    MatMFFDRegister - Adds a method to the MatMFFD registry.

175:    Not Collective

177:    Input Parameters:
178: +  name_solver - name of a new user-defined compute-h module
179: -  routine_create - routine to create method context

181:    Level: developer

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

186:    Sample usage:
187: .vb
188:    MatMFFDRegister("my_h",MyHCreate);
189: .ve

191:    Then, your solver can be chosen with the procedural interface via
192: $     MatMFFDSetType(mfctx,"my_h")
193:    or at runtime via the option
194: $     -mat_mffd_type my_h

196: .keywords: MatMFFD, register

198: .seealso: MatMFFDRegisterAll(), MatMFFDRegisterDestroy()
199:  @*/
200: PetscErrorCode  MatMFFDRegister(const char sname[],PetscErrorCode (*function)(MatMFFD))
201: {

205:   PetscFunctionListAdd(&MatMFFDList,sname,function);
206:   return(0);
207: }

211: PetscErrorCode  MatMFFDAddNullSpace_MFFD(Mat J,MatNullSpace nullsp)
212: {
214:   MatMFFD        ctx = (MatMFFD)J->data;

217:   PetscObjectReference((PetscObject)nullsp);
218:   if (ctx->sp) { MatNullSpaceDestroy(&ctx->sp); }
219:   ctx->sp = nullsp;
220:   return(0);
221: }

223: /* ----------------------------------------------------------------------------------------*/
226: PetscErrorCode MatDestroy_MFFD(Mat mat)
227: {
229:   MatMFFD        ctx = (MatMFFD)mat->data;

232:   VecDestroy(&ctx->w);
233:   VecDestroy(&ctx->drscale);
234:   VecDestroy(&ctx->dlscale);
235:   VecDestroy(&ctx->dshift);
236:   if (ctx->current_f_allocated) {
237:     VecDestroy(&ctx->current_f);
238:   }
239:   if (ctx->ops->destroy) {(*ctx->ops->destroy)(ctx);}
240:   MatNullSpaceDestroy(&ctx->sp);
241:   PetscHeaderDestroy(&ctx);
242:   mat->data = 0;

244:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetBase_C",NULL);
245:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctioniBase_C",NULL);
246:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctioni_C",NULL);
247:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunction_C",NULL);
248:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctionError_C",NULL);
249:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetCheckh_C",NULL);
250:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetPeriod_C",NULL);
251:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDResetHHistory_C",NULL);
252:   PetscObjectComposeFunction((PetscObject)mat,"MatMFFDAddNullSpace_C",NULL);
253:   return(0);
254: }

258: /*
259:    MatMFFDView_MFFD - Views matrix-free parameters.

261: */
262: PetscErrorCode MatView_MFFD(Mat J,PetscViewer viewer)
263: {
265:   MatMFFD        ctx = (MatMFFD)J->data;
266:   PetscBool      iascii, viewbase, viewfunction;
267:   const char     *prefix;

270:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
271:   if (iascii) {
272:     PetscViewerASCIIPrintf(viewer,"Matrix-free approximation:\n");
273:     PetscViewerASCIIPushTab(viewer);
274:     PetscViewerASCIIPrintf(viewer,"err=%g (relative error in function evaluation)\n",(double)ctx->error_rel);
275:     if (!((PetscObject)ctx)->type_name) {
276:       PetscViewerASCIIPrintf(viewer,"The compute h routine has not yet been set\n");
277:     } else {
278:       PetscViewerASCIIPrintf(viewer,"Using %s compute h routine\n",((PetscObject)ctx)->type_name);
279:     }
280:     if (ctx->ops->view) {
281:       (*ctx->ops->view)(ctx,viewer);
282:     }
283:     PetscObjectGetOptionsPrefix((PetscObject)J, &prefix);

285:     PetscOptionsHasName(prefix, "-mat_mffd_view_base", &viewbase);
286:     if (viewbase) {
287:       PetscViewerASCIIPrintf(viewer, "Base:\n");
288:       VecView(ctx->current_u, viewer);
289:     }
290:     PetscOptionsHasName(prefix, "-mat_mffd_view_function", &viewfunction);
291:     if (viewfunction) {
292:       PetscViewerASCIIPrintf(viewer, "Function:\n");
293:       VecView(ctx->current_f, viewer);
294:     }
295:     PetscViewerASCIIPopTab(viewer);
296:   }
297:   return(0);
298: }

302: /*
303:    MatAssemblyEnd_MFFD - Resets the ctx->ncurrenth to zero. This
304:    allows the user to indicate the beginning of a new linear solve by calling
305:    MatAssemblyXXX() on the matrix free matrix. This then allows the
306:    MatCreateMFFD_WP() to properly compute ||U|| only the first time
307:    in the linear solver rather than every time.

309:    This function is referenced directly from MatAssemblyEnd_SNESMF(), which may be in a different shared library.
310: */
311: PETSC_EXTERN PetscErrorCode MatAssemblyEnd_MFFD(Mat J,MatAssemblyType mt)
312: {
314:   MatMFFD        j = (MatMFFD)J->data;

317:   MatMFFDResetHHistory(J);
318:   j->vshift = 0.0;
319:   j->vscale = 1.0;
320:   return(0);
321: }

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

328:         y ~= (F(u + ha) - F(u))/h,
329:   where F = nonlinear function, as set by SNESSetFunction()
330:         u = current iterate
331:         h = difference interval
332: */
333: PetscErrorCode MatMult_MFFD(Mat mat,Vec a,Vec y)
334: {
335:   MatMFFD        ctx = (MatMFFD)mat->data;
336:   PetscScalar    h;
337:   Vec            w,U,F;
339:   PetscBool      zeroa;

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

349:   w = ctx->w;
350:   U = ctx->current_u;
351:   F = ctx->current_f;
352:   /*
353:       Compute differencing parameter
354:   */
355:   if (!ctx->ops->compute) {
356:     MatMFFDSetType(mat,MATMFFD_WP);
357:     MatSetFromOptions(mat);
358:   }
359:   (*ctx->ops->compute)(ctx,U,a,&h,&zeroa);
360:   if (zeroa) {
361:     VecSet(y,0.0);
362:     return(0);
363:   }

365:   if (PetscIsInfOrNanScalar(h)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Computed Nan differencing parameter h");
366:   if (ctx->checkh) {
367:     (*ctx->checkh)(ctx->checkhctx,U,a,&h);
368:   }

370:   /* keep a record of the current differencing parameter h */
371:   ctx->currenth = h;
372: #if defined(PETSC_USE_COMPLEX)
373:   PetscInfo2(mat,"Current differencing parameter: %g + %g i\n",(double)PetscRealPart(h),(double)PetscImaginaryPart(h));
374: #else
375:   PetscInfo1(mat,"Current differencing parameter: %15.12e\n",h);
376: #endif
377:   if (ctx->historyh && ctx->ncurrenth < ctx->maxcurrenth) {
378:     ctx->historyh[ctx->ncurrenth] = h;
379:   }
380:   ctx->ncurrenth++;

382:   /* w = u + ha */
383:   if (ctx->drscale) {
384:     VecPointwiseMult(ctx->drscale,a,U);
385:     VecAYPX(U,h,w);
386:   } else {
387:     VecWAXPY(w,h,a,U);
388:   }

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

396:   VecAXPY(y,-1.0,F);
397:   VecScale(y,1.0/h);

399:   if ((ctx->vshift != 0.0) || (ctx->vscale != 1.0)) {
400:     VecAXPBY(y,ctx->vshift,ctx->vscale,a);
401:   }
402:   if (ctx->dlscale) {
403:     VecPointwiseMult(y,ctx->dlscale,y);
404:   }
405:   if (ctx->dshift) {
406:     VecPointwiseMult(ctx->dshift,a,U);
407:     VecAXPY(y,1.0,U);
408:   }

410:   if (ctx->sp) {MatNullSpaceRemove(ctx->sp,y);}

412:   PetscLogEventEnd(MATMFFD_Mult,a,y,0,0);
413:   return(0);
414: }

418: /*
419:   MatGetDiagonal_MFFD - Gets the diagonal for a matrix free matrix

421:         y ~= (F(u + ha) - F(u))/h,
422:   where F = nonlinear function, as set by SNESSetFunction()
423:         u = current iterate
424:         h = difference interval
425: */
426: PetscErrorCode MatGetDiagonal_MFFD(Mat mat,Vec a)
427: {
428:   MatMFFD        ctx = (MatMFFD)mat->data;
429:   PetscScalar    h,*aa,*ww,v;
430:   PetscReal      epsilon = PETSC_SQRT_MACHINE_EPSILON,umin = 100.0*PETSC_SQRT_MACHINE_EPSILON;
431:   Vec            w,U;
433:   PetscInt       i,rstart,rend;

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

438:   w    = ctx->w;
439:   U    = ctx->current_u;
440:   (*ctx->func)(ctx->funcctx,U,a);
441:   (*ctx->funcisetbase)(ctx->funcctx,U);
442:   VecCopy(U,w);

444:   VecGetOwnershipRange(a,&rstart,&rend);
445:   VecGetArray(a,&aa);
446:   for (i=rstart; i<rend; i++) {
447:     VecGetArray(w,&ww);
448:     h    = ww[i-rstart];
449:     if (h == 0.0) h = 1.0;
450:     if (PetscAbsScalar(h) < umin && PetscRealPart(h) >= 0.0)     h = umin;
451:     else if (PetscRealPart(h) < 0.0 && PetscAbsScalar(h) < umin) h = -umin;
452:     h *= epsilon;

454:     ww[i-rstart] += h;
455:     VecRestoreArray(w,&ww);
456:     (*ctx->funci)(ctx->funcctx,i,w,&v);
457:     aa[i-rstart]  = (v - aa[i-rstart])/h;

459:     /* possibly shift and scale result */
460:     if ((ctx->vshift != 0.0) || (ctx->vscale != 1.0)) {
461:       aa[i - rstart] = ctx->vshift + ctx->vscale*aa[i-rstart];
462:     }

464:     VecGetArray(w,&ww);
465:     ww[i-rstart] -= h;
466:     VecRestoreArray(w,&ww);
467:   }
468:   VecRestoreArray(a,&aa);
469:   return(0);
470: }

474: PetscErrorCode MatDiagonalScale_MFFD(Mat mat,Vec ll,Vec rr)
475: {
476:   MatMFFD        aij = (MatMFFD)mat->data;

480:   if (ll && !aij->dlscale) {
481:     VecDuplicate(ll,&aij->dlscale);
482:   }
483:   if (rr && !aij->drscale) {
484:     VecDuplicate(rr,&aij->drscale);
485:   }
486:   if (ll) {
487:     VecCopy(ll,aij->dlscale);
488:   }
489:   if (rr) {
490:     VecCopy(rr,aij->drscale);
491:   }
492:   return(0);
493: }

497: PetscErrorCode MatDiagonalSet_MFFD(Mat mat,Vec ll,InsertMode mode)
498: {
499:   MatMFFD        aij = (MatMFFD)mat->data;

503:   if (mode == INSERT_VALUES) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"No diagonal set with INSERT_VALUES");
504:   if (!aij->dshift) {
505:     VecDuplicate(ll,&aij->dshift);
506:   }
507:   VecAXPY(aij->dshift,1.0,ll);
508:   return(0);
509: }

513: PetscErrorCode MatShift_MFFD(Mat Y,PetscScalar a)
514: {
515:   MatMFFD shell = (MatMFFD)Y->data;

518:   shell->vshift += a;
519:   return(0);
520: }

524: PetscErrorCode MatScale_MFFD(Mat Y,PetscScalar a)
525: {
526:   MatMFFD shell = (MatMFFD)Y->data;

529:   shell->vscale *= a;
530:   return(0);
531: }

535: PETSC_EXTERN PetscErrorCode MatMFFDSetBase_MFFD(Mat J,Vec U,Vec F)
536: {
538:   MatMFFD        ctx = (MatMFFD)J->data;

541:   MatMFFDResetHHistory(J);

543:   ctx->current_u = U;
544:   if (F) {
545:     if (ctx->current_f_allocated) {VecDestroy(&ctx->current_f);}
546:     ctx->current_f           = F;
547:     ctx->current_f_allocated = PETSC_FALSE;
548:   } else if (!ctx->current_f_allocated) {
549:     VecDuplicate(ctx->current_u, &ctx->current_f);

551:     ctx->current_f_allocated = PETSC_TRUE;
552:   }
553:   if (!ctx->w) {
554:     VecDuplicate(ctx->current_u, &ctx->w);
555:   }
556:   J->assembled = PETSC_TRUE;
557:   return(0);
558: }

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

564: PetscErrorCode  MatMFFDSetCheckh_MFFD(Mat J,FCN3 fun,void *ectx)
565: {
566:   MatMFFD ctx = (MatMFFD)J->data;

569:   ctx->checkh    = fun;
570:   ctx->checkhctx = ectx;
571:   return(0);
572: }

576: /*@C
577:    MatMFFDSetOptionsPrefix - Sets the prefix used for searching for all
578:    MatMFFD options in the database.

580:    Collective on Mat

582:    Input Parameter:
583: +  A - the Mat context
584: -  prefix - the prefix to prepend to all option names

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

590:    Level: advanced

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

594: .seealso: MatSetFromOptions(), MatCreateSNESMF()
595: @*/
596: PetscErrorCode  MatMFFDSetOptionsPrefix(Mat mat,const char prefix[])

598: {
599:   MatMFFD        mfctx = mat ? (MatMFFD)mat->data : (MatMFFD)NULL;

605:   PetscObjectSetOptionsPrefix((PetscObject)mfctx,prefix);
606:   return(0);
607: }

611: PetscErrorCode  MatSetFromOptions_MFFD(Mat mat)
612: {
613:   MatMFFD        mfctx = (MatMFFD)mat->data;
615:   PetscBool      flg;
616:   char           ftype[256];

621:   PetscObjectOptionsBegin((PetscObject)mfctx);
622:   PetscOptionsFList("-mat_mffd_type","Matrix free type","MatMFFDSetType",MatMFFDList,((PetscObject)mfctx)->type_name,ftype,256,&flg);
623:   if (flg) {
624:     MatMFFDSetType(mat,ftype);
625:   }

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

630:   flg  = PETSC_FALSE;
631:   PetscOptionsBool("-mat_mffd_check_positivity","Insure that U + h*a is nonnegative","MatMFFDSetCheckh",flg,&flg,NULL);
632:   if (flg) {
633:     MatMFFDSetCheckh(mat,MatMFFDCheckPositivity,0);
634:   }
635:   if (mfctx->ops->setfromoptions) {
636:     (*mfctx->ops->setfromoptions)(mfctx);
637:   }
638:   PetscOptionsEnd();
639:   return(0);
640: }

644: PetscErrorCode  MatMFFDSetPeriod_MFFD(Mat mat,PetscInt period)
645: {
646:   MatMFFD ctx = (MatMFFD)mat->data;

650:   ctx->recomputeperiod = period;
651:   return(0);
652: }

656: PetscErrorCode  MatMFFDSetFunction_MFFD(Mat mat,PetscErrorCode (*func)(void*,Vec,Vec),void *funcctx)
657: {
658:   MatMFFD ctx = (MatMFFD)mat->data;

661:   ctx->func    = func;
662:   ctx->funcctx = funcctx;
663:   return(0);
664: }

668: PetscErrorCode  MatMFFDSetFunctionError_MFFD(Mat mat,PetscReal error)
669: {
670:   MatMFFD ctx = (MatMFFD)mat->data;

674:   if (error != PETSC_DEFAULT) ctx->error_rel = error;
675:   return(0);
676: }

678: /*MC
679:   MATMFFD - MATMFFD = "mffd" - A matrix free matrix type.

681:   Level: advanced

683: .seealso: MatCreateMFFD(), MatCreateSNESMF(), MatMFFDSetFunction()
684: M*/
687: PETSC_EXTERN PetscErrorCode MatCreate_MFFD(Mat A)
688: {
689:   MatMFFD        mfctx;

693:   MatMFFDInitializePackage();

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

697:   mfctx->sp                       = 0;
698:   mfctx->error_rel                = PETSC_SQRT_MACHINE_EPSILON;
699:   mfctx->recomputeperiod          = 1;
700:   mfctx->count                    = 0;
701:   mfctx->currenth                 = 0.0;
702:   mfctx->historyh                 = NULL;
703:   mfctx->ncurrenth                = 0;
704:   mfctx->maxcurrenth              = 0;
705:   ((PetscObject)mfctx)->type_name = 0;

707:   mfctx->vshift = 0.0;
708:   mfctx->vscale = 1.0;

710:   /*
711:      Create the empty data structure to contain compute-h routines.
712:      These will be filled in below from the command line options or
713:      a later call with MatMFFDSetType() or if that is not called
714:      then it will default in the first use of MatMult_MFFD()
715:   */
716:   mfctx->ops->compute        = 0;
717:   mfctx->ops->destroy        = 0;
718:   mfctx->ops->view           = 0;
719:   mfctx->ops->setfromoptions = 0;
720:   mfctx->hctx                = 0;

722:   mfctx->func    = 0;
723:   mfctx->funcctx = 0;
724:   mfctx->w       = NULL;

726:   A->data = mfctx;

728:   A->ops->mult           = MatMult_MFFD;
729:   A->ops->destroy        = MatDestroy_MFFD;
730:   A->ops->view           = MatView_MFFD;
731:   A->ops->assemblyend    = MatAssemblyEnd_MFFD;
732:   A->ops->getdiagonal    = MatGetDiagonal_MFFD;
733:   A->ops->scale          = MatScale_MFFD;
734:   A->ops->shift          = MatShift_MFFD;
735:   A->ops->diagonalscale  = MatDiagonalScale_MFFD;
736:   A->ops->diagonalset    = MatDiagonalSet_MFFD;
737:   A->ops->setfromoptions = MatSetFromOptions_MFFD;
738:   A->assembled           = PETSC_TRUE;

740:   PetscLayoutSetUp(A->rmap);
741:   PetscLayoutSetUp(A->cmap);

743:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetBase_C",MatMFFDSetBase_MFFD);
744:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctioniBase_C",MatMFFDSetFunctioniBase_MFFD);
745:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctioni_C",MatMFFDSetFunctioni_MFFD);
746:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunction_C",MatMFFDSetFunction_MFFD);
747:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetCheckh_C",MatMFFDSetCheckh_MFFD);
748:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetPeriod_C",MatMFFDSetPeriod_MFFD);
749:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctionError_C",MatMFFDSetFunctionError_MFFD);
750:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDResetHHistory_C",MatMFFDResetHHistory_MFFD);
751:   PetscObjectComposeFunction((PetscObject)A,"MatMFFDAddNullSpace_C",MatMFFDAddNullSpace_MFFD);

753:   mfctx->mat = A;

755:   PetscObjectChangeTypeName((PetscObject)A,MATMFFD);
756:   return(0);
757: }

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

764:    Collective on Vec

766:    Input Parameters:
767: +  comm - MPI communicator
768: .  m - number of local rows (or PETSC_DECIDE to have calculated if M is given)
769:            This value should be the same as the local size used in creating the
770:            y vector for the matrix-vector product y = Ax.
771: .  n - This value should be the same as the local size used in creating the
772:        x vector for the matrix-vector product y = Ax. (or PETSC_DECIDE to have
773:        calculated if N is given) For square matrices n is almost always m.
774: .  M - number of global rows (or PETSC_DETERMINE to have calculated if m is given)
775: -  N - number of global columns (or PETSC_DETERMINE to have calculated if n is given)


778:    Output Parameter:
779: .  J - the matrix-free matrix

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


787:    Level: advanced

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

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

796: .vb
797:      F'(u)*a = [F(u+h*a) - F(u)]/h where
798:      h = error_rel*u'a/||a||^2                        if  |u'a| > umin*||a||_{1}
799:        = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2   otherwise
800:  where
801:      error_rel = square root of relative error in function evaluation
802:      umin = minimum iterate parameter
803: .ve

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

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

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

814:    Options Database Keys:
815: +  -mat_mffd_err <error_rel> - Sets error_rel
816: .  -mat_mffd_unim <umin> - Sets umin (for default PETSc routine that computes h only)
817: -  -mat_mffd_check_positivity

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

821: .seealso: MatDestroy(), MatMFFDSetFunctionError(), MatMFFDDSSetUmin(), MatMFFDSetFunction()
822:           MatMFFDSetHHistory(), MatMFFDResetHHistory(), MatCreateSNESMF(),
823:           MatMFFDGetH(), MatMFFDRegister(), MatMFFDComputeJacobian()

825: @*/
826: PetscErrorCode  MatCreateMFFD(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,Mat *J)
827: {

831:   MatCreate(comm,J);
832:   MatSetSizes(*J,m,n,M,N);
833:   MatSetType(*J,MATMFFD);
834:   MatSetUp(*J);
835:   return(0);
836: }


841: /*@
842:    MatMFFDGetH - Gets the last value that was used as the differencing
843:    parameter.

845:    Not Collective

847:    Input Parameters:
848: .  mat - the matrix obtained with MatCreateSNESMF()

850:    Output Paramter:
851: .  h - the differencing step size

853:    Level: advanced

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

857: .seealso: MatCreateSNESMF(),MatMFFDSetHHistory(), MatCreateMFFD(), MATMFFD, MatMFFDResetHHistory()
858: @*/
859: PetscErrorCode  MatMFFDGetH(Mat mat,PetscScalar *h)
860: {
861:   MatMFFD        ctx = (MatMFFD)mat->data;
863:   PetscBool      match;

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

869:   *h = ctx->currenth;
870:   return(0);
871: }

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

878:    Logically Collective on Mat

880:    Input Parameters:
881: +  mat - the matrix free matrix created via MatCreateSNESMF() or MatCreateMFFD()
882: .  func - the function to use
883: -  funcctx - optional function context passed to function

885:    Calling Sequence of func:
886: $     func (void *funcctx, Vec x, Vec f)

888: +  funcctx - user provided context
889: .  x - input vector
890: -  f - computed output function

892:    Level: advanced

894:    Notes:
895:     If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
896:     matrix inside your compute Jacobian routine

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

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

902: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD,
903:           MatMFFDSetHHistory(), MatMFFDResetHHistory(), SNESetFunction()
904: @*/
905: PetscErrorCode  MatMFFDSetFunction(Mat mat,PetscErrorCode (*func)(void*,Vec,Vec),void *funcctx)
906: {

910:   PetscTryMethod(mat,"MatMFFDSetFunction_C",(Mat,PetscErrorCode (*)(void*,Vec,Vec),void*),(mat,func,funcctx));
911:   return(0);
912: }

916: /*@C
917:    MatMFFDSetFunctioni - Sets the function for a single component

919:    Logically Collective on Mat

921:    Input Parameters:
922: +  mat - the matrix free matrix created via MatCreateSNESMF()
923: -  funci - the function to use

925:    Level: advanced

927:    Notes:
928:     If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
929:     matrix inside your compute Jacobian routine


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

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

936: @*/
937: PetscErrorCode  MatMFFDSetFunctioni(Mat mat,PetscErrorCode (*funci)(void*,PetscInt,Vec,PetscScalar*))
938: {

943:   PetscTryMethod(mat,"MatMFFDSetFunctioni_C",(Mat,PetscErrorCode (*)(void*,PetscInt,Vec,PetscScalar*)),(mat,funci));
944:   return(0);
945: }


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

953:    Logically Collective on Mat

955:    Input Parameters:
956: +  mat - the matrix free matrix created via MatCreateSNESMF()
957: -  func - the function to use

959:    Level: advanced

961:    Notes:
962:     If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
963:     matrix inside your compute Jacobian routine


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

968: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD
969:           MatMFFDSetHHistory(), MatMFFDResetHHistory(), SNESetFunction()
970: @*/
971: PetscErrorCode  MatMFFDSetFunctioniBase(Mat mat,PetscErrorCode (*func)(void*,Vec))
972: {

977:   PetscTryMethod(mat,"MatMFFDSetFunctioniBase_C",(Mat,PetscErrorCode (*)(void*,Vec)),(mat,func));
978:   return(0);
979: }

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

986:    Logically Collective on Mat

988:    Input Parameters:
989: +  mat - the matrix free matrix created via MatCreateSNESMF()
990: -  period - 1 for everytime, 2 for every second etc

992:    Options Database Keys:
993: +  -mat_mffd_period <period>

995:    Level: advanced


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

1000: .seealso: MatCreateSNESMF(),MatMFFDGetH(),
1001:           MatMFFDSetHHistory(), MatMFFDResetHHistory()
1002: @*/
1003: PetscErrorCode  MatMFFDSetPeriod(Mat mat,PetscInt period)
1004: {

1008:   PetscTryMethod(mat,"MatMFFDSetPeriod_C",(Mat,PetscInt),(mat,period));
1009:   return(0);
1010: }

1014: /*@
1015:    MatMFFDSetFunctionError - Sets the error_rel for the approximation of
1016:    matrix-vector products using finite differences.

1018:    Logically Collective on Mat

1020:    Input Parameters:
1021: +  mat - the matrix free matrix created via MatCreateMFFD() or MatCreateSNESMF()
1022: -  error_rel - relative error (should be set to the square root of
1023:                the relative error in the function evaluations)

1025:    Options Database Keys:
1026: +  -mat_mffd_err <error_rel> - Sets error_rel

1028:    Level: advanced

1030:    Notes:
1031:    The default matrix-free matrix-vector product routine computes
1032: .vb
1033:      F'(u)*a = [F(u+h*a) - F(u)]/h where
1034:      h = error_rel*u'a/||a||^2                        if  |u'a| > umin*||a||_{1}
1035:        = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2   else
1036: .ve

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

1040: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD
1041:           MatMFFDSetHHistory(), MatMFFDResetHHistory()
1042: @*/
1043: PetscErrorCode  MatMFFDSetFunctionError(Mat mat,PetscReal error)
1044: {

1048:   PetscTryMethod(mat,"MatMFFDSetFunctionError_C",(Mat,PetscReal),(mat,error));
1049:   return(0);
1050: }

1054: /*@
1055:    MatMFFDAddNullSpace - Provides a null space that an operator is
1056:    supposed to have.  Since roundoff will create a small component in
1057:    the null space, if you know the null space you may have it
1058:    automatically removed.

1060:    Logically Collective on Mat

1062:    Input Parameters:
1063: +  J - the matrix-free matrix context
1064: -  nullsp - object created with MatNullSpaceCreate()

1066:    Level: advanced

1068: .keywords: SNES, matrix-free, null space

1070: .seealso: MatNullSpaceCreate(), MatMFFDGetH(), MatCreateSNESMF(), MatCreateMFFD(), MATMFFD
1071:           MatMFFDSetHHistory(), MatMFFDResetHHistory()
1072: @*/
1073: PetscErrorCode  MatMFFDAddNullSpace(Mat J,MatNullSpace nullsp)
1074: {

1078:   PetscTryMethod(J,"MatMFFDAddNullSpace_C",(Mat,MatNullSpace),(J,nullsp));
1079:   return(0);
1080: }

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

1088:    Logically Collective on Mat

1090:    Input Parameters:
1091: +  J - the matrix-free matrix context
1092: .  histroy - space to hold the history
1093: -  nhistory - number of entries in history, if more entries are generated than
1094:               nhistory, then the later ones are discarded

1096:    Level: advanced

1098:    Notes:
1099:    Use MatMFFDResetHHistory() to reset the history counter and collect
1100:    a new batch of differencing parameters, h.

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

1104: .seealso: MatMFFDGetH(), MatCreateSNESMF(),
1105:           MatMFFDResetHHistory(), MatMFFDSetFunctionError()

1107: @*/
1108: PetscErrorCode  MatMFFDSetHHistory(Mat J,PetscScalar history[],PetscInt nhistory)
1109: {
1110:   MatMFFD        ctx = (MatMFFD)J->data;
1112:   PetscBool      match;

1115:   PetscObjectTypeCompare((PetscObject)J,MATMFFD,&match);
1116:   if (!match) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_WRONG,"Not a MFFD matrix");
1117:   ctx->historyh    = history;
1118:   ctx->maxcurrenth = nhistory;
1119:   ctx->currenth    = 0.;
1120:   return(0);
1121: }


1126: /*@
1127:    MatMFFDResetHHistory - Resets the counter to zero to begin
1128:    collecting a new set of differencing histories.

1130:    Logically Collective on Mat

1132:    Input Parameters:
1133: .  J - the matrix-free matrix context

1135:    Level: advanced

1137:    Notes:
1138:    Use MatMFFDSetHHistory() to create the original history counter.

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

1142: .seealso: MatMFFDGetH(), MatCreateSNESMF(),
1143:           MatMFFDSetHHistory(), MatMFFDSetFunctionError()

1145: @*/
1146: PetscErrorCode  MatMFFDResetHHistory(Mat J)
1147: {

1151:   PetscTryMethod(J,"MatMFFDResetHHistory_C",(Mat),(J));
1152:   return(0);
1153: }


1158: /*@
1159:     MatMFFDSetBase - Sets the vector U at which matrix vector products of the
1160:         Jacobian are computed

1162:     Logically Collective on Mat

1164:     Input Parameters:
1165: +   J - the MatMFFD matrix
1166: .   U - the vector
1167: -   F - (optional) vector that contains F(u) if it has been already computed

1169:     Notes: This is rarely used directly

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

1174:     Level: advanced

1176: @*/
1177: PetscErrorCode  MatMFFDSetBase(Mat J,Vec U,Vec F)
1178: {

1185:   PetscTryMethod(J,"MatMFFDSetBase_C",(Mat,Vec,Vec),(J,U,F));
1186:   return(0);
1187: }

1191: /*@C
1192:     MatMFFDSetCheckh - Sets a function that checks the computed h and adjusts
1193:         it to satisfy some criteria

1195:     Logically Collective on Mat

1197:     Input Parameters:
1198: +   J - the MatMFFD matrix
1199: .   fun - the function that checks h
1200: -   ctx - any context needed by the function

1202:     Options Database Keys:
1203: .   -mat_mffd_check_positivity

1205:     Level: advanced

1207:     Notes: For example, MatMFFDSetCheckPositivity() insures that all entries
1208:        of U + h*a are non-negative

1210: .seealso:  MatMFFDSetCheckPositivity()
1211: @*/
1212: PetscErrorCode  MatMFFDSetCheckh(Mat J,PetscErrorCode (*fun)(void*,Vec,Vec,PetscScalar*),void *ctx)
1213: {

1218:   PetscTryMethod(J,"MatMFFDSetCheckh_C",(Mat,PetscErrorCode (*)(void*,Vec,Vec,PetscScalar*),void*),(J,fun,ctx));
1219:   return(0);
1220: }

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

1228:     Logically Collective on Vec

1230:     Input Parameters:
1231: +   U - base vector that is added to
1232: .   a - vector that is added
1233: .   h - scaling factor on a
1234: -   dummy - context variable (unused)

1236:     Options Database Keys:
1237: .   -mat_mffd_check_positivity

1239:     Level: advanced

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

1244: .seealso:  MatMFFDSetCheckh()
1245: @*/
1246: PetscErrorCode  MatMFFDCheckPositivity(void *dummy,Vec U,Vec a,PetscScalar *h)
1247: {
1248:   PetscReal      val, minval;
1249:   PetscScalar    *u_vec, *a_vec;
1251:   PetscInt       i,n;
1252:   MPI_Comm       comm;

1255:   PetscObjectGetComm((PetscObject)U,&comm);
1256:   VecGetArray(U,&u_vec);
1257:   VecGetArray(a,&a_vec);
1258:   VecGetLocalSize(U,&n);
1259:   minval = PetscAbsScalar(*h*1.01);
1260:   for (i=0; i<n; i++) {
1261:     if (PetscRealPart(u_vec[i] + *h*a_vec[i]) <= 0.0) {
1262:       val = PetscAbsScalar(u_vec[i]/a_vec[i]);
1263:       if (val < minval) minval = val;
1264:     }
1265:   }
1266:   VecRestoreArray(U,&u_vec);
1267:   VecRestoreArray(a,&a_vec);
1268:   MPI_Allreduce(&minval,&val,1,MPIU_REAL,MPIU_MIN,comm);
1269:   if (val <= PetscAbsScalar(*h)) {
1270:     PetscInfo2(U,"Scaling back h from %g to %g\n",(double)PetscRealPart(*h),(double)(.99*val));
1271:     if (PetscRealPart(*h) > 0.0) *h =  0.99*val;
1272:     else                         *h = -0.99*val;
1273:   }
1274:   return(0);
1275: }