Actual source code: snesmfj.c

  1: /*$Id: snesmfj.c,v 1.131 2001/09/05 18:45:40 bsmith Exp $*/

 3:  #include src/mat/matimpl.h
 4:  #include src/snes/mf/snesmfj.h

  6: PetscFList      MatSNESMPetscFList              = 0;
  7: PetscTruth MatSNESMFRegisterAllCalled = PETSC_FALSE;

  9: #undef __FUNCT__  
 11: /*@C
 12:     MatSNESMFSetType - Sets the method that is used to compute the 
 13:     differencing parameter for finite differene matrix-free formulations. 

 15:     Input Parameters:
 16: +   mat - the "matrix-free" matrix created via MatCreateSNESMF(), or MatCreateMF()
 17:           or MatSetType(mat,MATMFFD);
 18: -   ftype - the type requested

 20:     Level: advanced

 22:     Notes:
 23:     For example, such routines can compute h for use in
 24:     Jacobian-vector products of the form

 26:                         F(x+ha) - F(x)
 27:           F'(u)a  ~=  ----------------
 28:                               h

 30: .seealso: MatCreateSNESMF(), MatSNESMFRegisterDynamic)
 31: @*/
 32: int MatSNESMFSetType(Mat mat,MatSNESMFType ftype)
 33: {
 34:   int          ierr,(*r)(MatSNESMFCtx);
 35:   MatSNESMFCtx ctx = (MatSNESMFCtx)mat->data;
 36:   PetscTruth   match;
 37: 

 42:   /* already set, so just return */
 43:   PetscTypeCompare((PetscObject)ctx,ftype,&match);
 44:   if (match) return(0);

 46:   /* destroy the old one if it exists */
 47:   if (ctx->ops->destroy) {
 48:     (*ctx->ops->destroy)(ctx);
 49:   }

 51:   /* Get the function pointers for the requrested method */
 52:   if (!MatSNESMFRegisterAllCalled) {MatSNESMFRegisterAll(PETSC_NULL);}

 54:    PetscFListFind(ctx->comm,MatSNESMPetscFList,ftype,(void (**)(void)) &r);

 56:   if (!r) SETERRQ(1,"Unknown MatSNESMF type given");

 58:   (*r)(ctx);

 60:   PetscObjectChangeTypeName((PetscObject)ctx,ftype);

 62:   return(0);
 63: }

 65: EXTERN_C_BEGIN
 66: #undef __FUNCT__  
 68: int MatSNESMFSetFunctioniBase_FD(Mat mat,int (*func)(Vec,void *))
 69: {
 70:   MatSNESMFCtx ctx = (MatSNESMFCtx)mat->data;

 73:   ctx->funcisetbase = func;
 74:   return(0);
 75: }
 76: EXTERN_C_END

 78: EXTERN_C_BEGIN
 79: #undef __FUNCT__  
 81: int MatSNESMFSetFunctioni_FD(Mat mat,int (*funci)(int,Vec,PetscScalar*,void *))
 82: {
 83:   MatSNESMFCtx ctx = (MatSNESMFCtx)mat->data;

 86:   ctx->funci = funci;
 87:   return(0);
 88: }
 89: EXTERN_C_END

 91: /*MC
 92:    MatSNESMFRegisterDynamic - Adds a method to the MatSNESMF registry.

 94:    Synopsis:
 95:    int MatSNESMFRegisterDynamic(char *name_solver,char *path,char *name_create,int (*routine_create)(MatSNESMF))

 97:    Not Collective

 99:    Input Parameters:
100: +  name_solver - name of a new user-defined compute-h module
101: .  path - path (either absolute or relative) the library containing this solver
102: .  name_create - name of routine to create method context
103: -  routine_create - routine to create method context

105:    Level: developer

107:    Notes:
108:    MatSNESMFRegisterDynamic) may be called multiple times to add several user-defined solvers.

110:    If dynamic libraries are used, then the fourth input argument (routine_create)
111:    is ignored.

113:    Sample usage:
114: .vb
115:    MatSNESMFRegisterDynamic"my_h",/home/username/my_lib/lib/libO/solaris/mylib.a,
116:                "MyHCreate",MyHCreate);
117: .ve

119:    Then, your solver can be chosen with the procedural interface via
120: $     MatSNESMFSetType(mfctx,"my_h")
121:    or at runtime via the option
122: $     -snes_mf_type my_h

124: .keywords: MatSNESMF, register

126: .seealso: MatSNESMFRegisterAll(), MatSNESMFRegisterDestroy()
127: M*/

129: #undef __FUNCT__  
131: int MatSNESMFRegister(char *sname,char *path,char *name,int (*function)(MatSNESMFCtx))
132: {
134:   char fullname[256];

137:   PetscFListConcat(path,name,fullname);
138:   PetscFListAdd(&MatSNESMPetscFList,sname,fullname,(void (*)(void))function);
139:   return(0);
140: }


143: #undef __FUNCT__  
145: /*@C
146:    MatSNESMFRegisterDestroy - Frees the list of MatSNESMF methods that were
147:    registered by MatSNESMFRegisterDynamic).

149:    Not Collective

151:    Level: developer

153: .keywords: MatSNESMF, register, destroy

155: .seealso: MatSNESMFRegisterDynamic), MatSNESMFRegisterAll()
156: @*/
157: int MatSNESMFRegisterDestroy(void)
158: {

162:   if (MatSNESMPetscFList) {
163:     PetscFListDestroy(&MatSNESMPetscFList);
164:     MatSNESMPetscFList = 0;
165:   }
166:   MatSNESMFRegisterAllCalled = PETSC_FALSE;
167:   return(0);
168: }

170: /* ----------------------------------------------------------------------------------------*/
171: #undef __FUNCT__  
173: int MatDestroy_MFFD(Mat mat)
174: {
175:   int          ierr;
176:   MatSNESMFCtx ctx = (MatSNESMFCtx)mat->data;

179:   if (ctx->w != PETSC_NULL) {
180:     VecDestroy(ctx->w);
181:   }
182:   if (ctx->ops->destroy) {(*ctx->ops->destroy)(ctx);}
183:   if (ctx->sp) {MatNullSpaceDestroy(ctx->sp);}
184:   PetscHeaderDestroy(ctx);
185:   return(0);
186: }

188: #undef __FUNCT__  
190: /*
191:    MatSNESMFView_MFFD - Views matrix-free parameters.

193: */
194: int MatView_MFFD(Mat J,PetscViewer viewer)
195: {
196:   int          ierr;
197:   MatSNESMFCtx ctx = (MatSNESMFCtx)J->data;
198:   PetscTruth   isascii;

201:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);
202:   if (isascii) {
203:      PetscViewerASCIIPrintf(viewer,"  SNES matrix-free approximation:n");
204:      PetscViewerASCIIPrintf(viewer,"    err=%g (relative error in function evaluation)n",ctx->error_rel);
205:      if (!ctx->type_name) {
206:        PetscViewerASCIIPrintf(viewer,"    The compute h routine has not yet been setn");
207:      } else {
208:        PetscViewerASCIIPrintf(viewer,"    Using %s compute h routinen",ctx->type_name);
209:      }
210:      if (ctx->ops->view) {
211:        (*ctx->ops->view)(ctx,viewer);
212:      }
213:   } else {
214:     SETERRQ1(1,"Viewer type %s not supported for SNES matrix free matrix",((PetscObject)viewer)->type_name);
215:   }
216:   return(0);
217: }

219: #undef __FUNCT__  
221: /*
222:    MatSNESMFAssemblyEnd_Private - Resets the ctx->ncurrenth to zero. This 
223:    allows the user to indicate the beginning of a new linear solve by calling
224:    MatAssemblyXXX() on the matrix free matrix. This then allows the 
225:    MatSNESMFCreate_WP() to properly compute ||U|| only the first time
226:    in the linear solver rather than every time.
227: */
228: int MatAssemblyEnd_MFFD(Mat J,MatAssemblyType mt)
229: {
230:   int             ierr;
231:   MatSNESMFCtx    j = (MatSNESMFCtx)J->data;

234:   MatSNESMFResetHHistory(J);
235:   if (j->usesnes) {
236:     SNESGetSolution(j->snes,&j->current_u);
237:     SNESGetFunction(j->snes,&j->current_f,PETSC_NULL,PETSC_NULL);
238:     if (j->w == PETSC_NULL) {
239:       VecDuplicate(j->current_u, &j->w);
240:     }
241:   }
242:   j->vshift = 0.0;
243:   j->vscale = 1.0;
244:   return(0);
245: }

247: #undef __FUNCT__  
249: /*
250:   MatMult_MFFD - Default matrix-free form for Jacobian-vector product, y = F'(u)*a:

252:         y ~= (F(u + ha) - F(u))/h, 
253:   where F = nonlinear function, as set by SNESSetFunction()
254:         u = current iterate
255:         h = difference interval
256: */
257: int MatMult_MFFD(Mat mat,Vec a,Vec y)
258: {
259:   MatSNESMFCtx    ctx = (MatSNESMFCtx)mat->data;
260:   SNES            snes;
261:   PetscScalar     h,mone = -1.0;
262:   Vec             w,U,F;
263:   int             ierr,(*eval_fct)(SNES,Vec,Vec)=0;

266:   /* We log matrix-free matrix-vector products separately, so that we can
267:      separate the performance monitoring from the cases that use conventional
268:      storage.  We may eventually modify event logging to associate events
269:      with particular objects, hence alleviating the more general problem. */
270:   PetscLogEventBegin(MAT_MultMatrixFree,a,y,0,0);

272:   snes = ctx->snes;
273:   w    = ctx->w;
274:   U    = ctx->current_u;

276:   /* 
277:       Compute differencing parameter 
278:   */
279:   if (!ctx->ops->compute) {
280:     MatSNESMFSetType(mat,MATSNESMF_DEFAULT);
281:     MatSNESMFSetFromOptions(mat);
282:   }
283:   (*ctx->ops->compute)(ctx,U,a,&h);

285:   if (ctx->checkh) {
286:     (*ctx->checkh)(U,a,&h,ctx->checkhctx);
287:   }

289:   /* keep a record of the current differencing parameter h */
290:   ctx->currenth = h;
291: #if defined(PETSC_USE_COMPLEX)
292:   PetscLogInfo(mat,"MatMult_MFFD:Current differencing parameter: %g + %g in",PetscRealPart(h),PetscImaginaryPart(h));
293: #else
294:   PetscLogInfo(mat,"MatMult_MFFD:Current differencing parameter: %15.12en",h);
295: #endif
296:   if (ctx->historyh && ctx->ncurrenth < ctx->maxcurrenth) {
297:     ctx->historyh[ctx->ncurrenth] = h;
298:   }
299:   ctx->ncurrenth++;

301:   /* w = u + ha */
302:   VecWAXPY(&h,a,U,w);

304:   if (ctx->usesnes) {
305:     eval_fct = SNESComputeFunction;
306:     F    = ctx->current_f;
307:     if (!F) SETERRQ(1,"You must call MatAssembly() even on matrix-free matrices");
308:     (*eval_fct)(snes,w,y);
309:   } else {
310:     F = ctx->funcvec;
311:     /* compute func(U) as base for differencing */
312:     if (ctx->ncurrenth == 1) {
313:       (*ctx->func)(snes,U,F,ctx->funcctx);
314:     }
315:     (*ctx->func)(snes,w,y,ctx->funcctx);
316:   }

318:   VecAXPY(&mone,F,y);
319:   h    = 1.0/h;
320:   VecScale(&h,y);


323:   if (ctx->vshift != 0.0 && ctx->vscale != 1.0) {
324:     VecAXPBY(&ctx->vshift,&ctx->vscale,a,y);
325:   } else if (ctx->vscale != 1.0) {
326:     VecScale(&ctx->vscale,y);
327:   } else if (ctx->vshift != 0.0) {
328:     VecAXPY(&ctx->vshift,a,y);
329:   }

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

333:   PetscLogEventEnd(MAT_MultMatrixFree,a,y,0,0);
334:   return(0);
335: }

337: #undef __FUNCT__  
339: /*
340:   MatGetDiagonal_MFFD - Gets the diagonal for a matrix free matrix

342:         y ~= (F(u + ha) - F(u))/h, 
343:   where F = nonlinear function, as set by SNESSetFunction()
344:         u = current iterate
345:         h = difference interval
346: */
347: int MatGetDiagonal_MFFD(Mat mat,Vec a)
348: {
349:   MatSNESMFCtx ctx = (MatSNESMFCtx)mat->data;
350:   PetscScalar  h,*aa,*ww,v;
351:   PetscReal    epsilon = PETSC_SQRT_MACHINE_EPSILON,umin = 100.0*PETSC_SQRT_MACHINE_EPSILON;
352:   Vec          w,U;
353:   int          i,ierr,rstart,rend;

356:   if (!ctx->funci) {
357:     SETERRQ(1,"Requirers calling MatSNESMFSetFunctioni() first");
358:   }

360:   w    = ctx->w;
361:   U    = ctx->current_u;
362:   (*ctx->func)(0,U,a,ctx->funcctx);
363:   (*ctx->funcisetbase)(U,ctx->funcctx);
364:   VecCopy(U,w);

366:   VecGetOwnershipRange(a,&rstart,&rend);
367:   VecGetArray(a,&aa);
368:   for (i=rstart; i<rend; i++) {
369:     VecGetArray(w,&ww);
370:     h  = ww[i-rstart];
371:     if (h == 0.0) h = 1.0;
372: #if !defined(PETSC_USE_COMPLEX)
373:     if (h < umin && h >= 0.0)      h = umin;
374:     else if (h < 0.0 && h > -umin) h = -umin;
375: #else
376:     if (PetscAbsScalar(h) < umin && PetscRealPart(h) >= 0.0)     h = umin;
377:     else if (PetscRealPart(h) < 0.0 && PetscAbsScalar(h) < umin) h = -umin;
378: #endif
379:     h     *= epsilon;
380: 
381:     ww[i-rstart] += h;
382:     VecRestoreArray(w,&ww);
383:     ierr          = (*ctx->funci)(i,w,&v,ctx->funcctx);
384:     aa[i-rstart]  = (v - aa[i-rstart])/h;

386:     /* possibly shift and scale result */
387:     aa[i - rstart] = ctx->vshift + ctx->vscale*aa[i-rstart];

389:     VecGetArray(w,&ww);
390:     ww[i-rstart] -= h;
391:     VecRestoreArray(w,&ww);
392:   }
393:   VecRestoreArray(a,&aa);
394:   return(0);
395: }

397: #undef __FUNCT__  
399: int MatShift_MFFD(PetscScalar *a,Mat Y)
400: {
401:   MatSNESMFCtx shell = (MatSNESMFCtx)Y->data;
403:   shell->vshift += *a;
404:   return(0);
405: }

407: #undef __FUNCT__  
409: int MatScale_MFFD(PetscScalar *a,Mat Y)
410: {
411:   MatSNESMFCtx shell = (MatSNESMFCtx)Y->data;
413:   shell->vscale *= *a;
414:   return(0);
415: }


418: #undef __FUNCT__  
420: /*@C
421:    MatCreateSNESMF - Creates a matrix-free matrix context for use with
422:    a SNES solver.  This matrix can be used as the Jacobian argument for
423:    the routine SNESSetJacobian().

425:    Collective on SNES and Vec

427:    Input Parameters:
428: +  snes - the SNES context
429: -  x - vector where SNES solution is to be stored.

431:    Output Parameter:
432: .  J - the matrix-free matrix

434:    Level: advanced

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

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

443: .vb
444:      F'(u)*a = [F(u+h*a) - F(u)]/h where
445:      h = error_rel*u'a/||a||^2                        if  |u'a| > umin*||a||_{1}
446:        = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2   otherwise
447:  where
448:      error_rel = square root of relative error in function evaluation
449:      umin = minimum iterate parameter
450: .ve

452:    The user can set the error_rel via MatSNESMFSetFunctionError() and 
453:    umin via MatSNESMFDefaultSetUmin(); see the nonlinear solvers chapter
454:    of the users manual for details.

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

459:    Options Database Keys:
460: +  -snes_mf_err <error_rel> - Sets error_rel
461: .  -snes_mf_unim <umin> - Sets umin (for default PETSc routine that computes h only)
462: -  -snes_mf_ksp_monitor - KSP monitor routine that prints differencing h

464: .keywords: SNES, default, matrix-free, create, matrix

466: .seealso: MatDestroy(), MatSNESMFSetFunctionError(), MatSNESMFDefaultSetUmin()
467:           MatSNESMFSetHHistory(), MatSNESMFResetHHistory(), MatCreateMF(),
468:           MatSNESMFGetH(),MatSNESMFKSPMonitor(), MatSNESMFRegisterDynamic), MatSNESMFComputeJacobian()
469:  
470: @*/
471: int MatCreateSNESMF(SNES snes,Vec x,Mat *J)
472: {
473:   MatSNESMFCtx mfctx;
474:   int          ierr;

477:   MatCreateMF(x,J);

479:   mfctx          = (MatSNESMFCtx)(*J)->data;
480:   mfctx->snes    = snes;
481:   mfctx->usesnes = PETSC_TRUE;
482:   PetscLogObjectParent(snes,*J);
483:   return(0);
484: }

486: EXTERN_C_BEGIN
487: #undef __FUNCT__  
489: int MatSNESMFSetBase_FD(Mat J,Vec U)
490: {
491:   int          ierr;
492:   MatSNESMFCtx ctx = (MatSNESMFCtx)J->data;

495:   MatSNESMFResetHHistory(J);
496:   ctx->current_u = U;
497:   ctx->usesnes   = PETSC_FALSE;
498:   if (ctx->w == PETSC_NULL) {
499:     VecDuplicate(ctx->current_u, &ctx->w);
500:   }
501:   return(0);
502: }
503: EXTERN_C_END

505: EXTERN_C_BEGIN
506: #undef __FUNCT__  
508: int MatSNESMFSetCheckh_FD(Mat J,int (*fun)(Vec,Vec,PetscScalar*,void*),void*ectx)
509: {
510:   MatSNESMFCtx ctx = (MatSNESMFCtx)J->data;

513:   ctx->checkh    = fun;
514:   ctx->checkhctx = ectx;
515:   return(0);
516: }
517: EXTERN_C_END

519: #undef __FUNCT__  
521: /*@
522:    MatSNESMFSetFromOptions - Sets the MatSNESMF options from the command line
523:    parameter.

525:    Collective on Mat

527:    Input Parameters:
528: .  mat - the matrix obtained with MatCreateSNESMF()

530:    Options Database Keys:
531: +  -snes_mf_type - <default,wp>
532: -  -snes_mf_err - square root of estimated relative error in function evaluation
533: -  -snes_mf_period - how often h is recomputed, defaults to 1, everytime

535:    Level: advanced

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

539: .seealso: MatCreateSNESMF(),MatSNESMFSetHHistory(), 
540:           MatSNESMFResetHHistory(), MatSNESMFKSPMonitor()
541: @*/
542: int MatSNESMFSetFromOptions(Mat mat)
543: {
544:   MatSNESMFCtx mfctx = (MatSNESMFCtx)mat->data;
545:   int          ierr;
546:   PetscTruth   flg;
547:   char         ftype[256];

550:   if (!MatSNESMFRegisterAllCalled) {MatSNESMFRegisterAll(PETSC_NULL);}
551: 
552:   PetscOptionsBegin(mfctx->comm,mfctx->prefix,"Set matrix free computation parameters","MatSNESMF");
553:   PetscOptionsList("-snes_mf_type","Matrix free type","MatSNESMFSetType",MatSNESMPetscFList,mfctx->type_name,ftype,256,&flg);
554:   if (flg) {
555:     MatSNESMFSetType(mat,ftype);
556:   }

558:   PetscOptionsReal("-snes_mf_err","set sqrt relative error in function","MatSNESMFSetFunctionError",mfctx->error_rel,&mfctx->error_rel,0);
559:   PetscOptionsInt("-snes_mf_period","how often h is recomputed","MatSNESMFSetPeriod",mfctx->recomputeperiod,&mfctx->recomputeperiod,0);
560:   if (mfctx->snes) {
561:     PetscOptionsName("-snes_mf_ksp_monitor","Monitor matrix-free parameters","MatSNESMFKSPMonitor",&flg);
562:     if (flg) {
563:       SLES sles;
564:       KSP  ksp;
565:       SNESGetSLES(mfctx->snes,&sles);
566:       SLESGetKSP(sles,&ksp);
567:       KSPSetMonitor(ksp,MatSNESMFKSPMonitor,PETSC_NULL,0);
568:     }
569:   }
570:   PetscOptionsName("-snes_mf_check_positivity","Insure that U + h*a is nonnegative","MatSNESMFSetCheckh",&flg);
571:   if (flg) {
572:     MatSNESMFSetCheckh(mat,MatSNESMFCheckPositivity,0);
573:   }
574:   if (mfctx->ops->setfromoptions) {
575:     (*mfctx->ops->setfromoptions)(mfctx);
576:   }
577:   PetscOptionsEnd();
578:   return(0);
579: }

581: #undef __FUNCT__  
583: EXTERN_C_BEGIN
584: int MatCreate_MFFD(Mat A)
585: {
586:   MatSNESMFCtx mfctx;
587:   int          ierr;

590: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
591:   SNESInitializePackage(PETSC_NULL);
592: #endif

594:   PetscHeaderCreate(mfctx,_p_MatSNESMFCtx,struct _MFOps,MATSNESMFCTX_COOKIE,0,"SNESMF",A->comm,MatDestroy_MFFD,MatView_MFFD);
595:   PetscLogObjectCreate(mfctx);
596:   mfctx->sp              = 0;
597:   mfctx->snes            = 0;
598:   mfctx->error_rel       = PETSC_SQRT_MACHINE_EPSILON;
599:   mfctx->recomputeperiod = 1;
600:   mfctx->count           = 0;
601:   mfctx->currenth        = 0.0;
602:   mfctx->historyh        = PETSC_NULL;
603:   mfctx->ncurrenth       = 0;
604:   mfctx->maxcurrenth     = 0;
605:   mfctx->type_name       = 0;
606:   mfctx->usesnes         = PETSC_FALSE;

608:   mfctx->vshift          = 0.0;
609:   mfctx->vscale          = 1.0;

611:   /* 
612:      Create the empty data structure to contain compute-h routines.
613:      These will be filled in below from the command line options or 
614:      a later call with MatSNESMFSetType() or if that is not called 
615:      then it will default in the first use of MatMult_MFFD()
616:   */
617:   mfctx->ops->compute        = 0;
618:   mfctx->ops->destroy        = 0;
619:   mfctx->ops->view           = 0;
620:   mfctx->ops->setfromoptions = 0;
621:   mfctx->hctx                = 0;

623:   mfctx->func                = 0;
624:   mfctx->funcctx             = 0;
625:   mfctx->funcvec             = 0;
626:   mfctx->w                   = PETSC_NULL;

628:   A->data                = mfctx;

630:   A->ops->mult           = MatMult_MFFD;
631:   A->ops->destroy        = MatDestroy_MFFD;
632:   A->ops->view           = MatView_MFFD;
633:   A->ops->assemblyend    = MatAssemblyEnd_MFFD;
634:   A->ops->getdiagonal    = MatGetDiagonal_MFFD;
635:   A->ops->scale          = MatScale_MFFD;
636:   A->ops->shift          = MatShift_MFFD;
637:   A->ops->setfromoptions = MatSNESMFSetFromOptions;

639:   PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSNESMFSetBase_C","MatSNESMFSetBase_FD",MatSNESMFSetBase_FD);
640:   PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSNESMFSetFunctioniBase_C","MatSNESMFSetFunctioniBase_FD",MatSNESMFSetFunctioniBase_FD);
641:   PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSNESMFSetFunctioni_C","MatSNESMFSetFunctioni_FD",MatSNESMFSetFunctioni_FD);
642:   PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSNESMFSetCheckh_C","MatSNESMFSetCheckh_FD",MatSNESMFSetCheckh_FD);
643:   mfctx->mat = A;

645:   return(0);
646: }

648: EXTERN_C_END

650: #undef __FUNCT__  
652: /*@C
653:    MatCreateMF - Creates a matrix-free matrix. See also MatCreateSNESMF() 

655:    Collective on Vec

657:    Input Parameters:
658: .  x - vector that defines layout of the vectors and matrices

660:    Output Parameter:
661: .  J - the matrix-free matrix

663:    Level: advanced

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

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

672: .vb
673:      F'(u)*a = [F(u+h*a) - F(u)]/h where
674:      h = error_rel*u'a/||a||^2                        if  |u'a| > umin*||a||_{1}
675:        = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2   otherwise
676:  where
677:      error_rel = square root of relative error in function evaluation
678:      umin = minimum iterate parameter
679: .ve

681:    The user can set the error_rel via MatSNESMFSetFunctionError() and 
682:    umin via MatSNESMFDefaultSetUmin(); see the nonlinear solvers chapter
683:    of the users manual for details.

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

688:    Options Database Keys:
689: +  -snes_mf_err <error_rel> - Sets error_rel
690: .  -snes_mf_unim <umin> - Sets umin (for default PETSc routine that computes h only)
691: .  -snes_mf_ksp_monitor - KSP monitor routine that prints differencing h
692: -  -snes_mf_check_positivity

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

696: .seealso: MatDestroy(), MatSNESMFSetFunctionError(), MatSNESMFDefaultSetUmin()
697:           MatSNESMFSetHHistory(), MatSNESMFResetHHistory(), MatCreateSNESMF(),
698:           MatSNESMFGetH(),MatSNESMFKSPMonitor(), MatSNESMFRegisterDynamic),, MatSNESMFComputeJacobian()
699:  
700: @*/
701: int MatCreateMF(Vec x,Mat *J)
702: {
703:   MPI_Comm     comm;
704:   int          n,nloc,ierr;

707:   PetscObjectGetComm((PetscObject)x,&comm);
708:   VecGetSize(x,&n);
709:   VecGetLocalSize(x,&nloc);
710:   MatCreate(comm,nloc,nloc,n,n,J);
711:   MatRegister(MATMFFD,0,"MatCreate_MFFD",MatCreate_MFFD);
712:   MatSetType(*J,MATMFFD);
713:   return(0);
714: }


717: #undef __FUNCT__  
719: /*@
720:    MatSNESMFGetH - Gets the last value that was used as the differencing 
721:    parameter.

723:    Not Collective

725:    Input Parameters:
726: .  mat - the matrix obtained with MatCreateSNESMF()

728:    Output Paramter:
729: .  h - the differencing step size

731:    Level: advanced

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

735: .seealso: MatCreateSNESMF(),MatSNESMFSetHHistory(), 
736:           MatSNESMFResetHHistory(),MatSNESMFKSPMonitor()
737: @*/
738: int MatSNESMFGetH(Mat mat,PetscScalar *h)
739: {
740:   MatSNESMFCtx ctx = (MatSNESMFCtx)mat->data;

743:   *h = ctx->currenth;
744:   return(0);
745: }

747: #undef __FUNCT__  
749: /*
750:    MatSNESMFKSPMonitor - A KSP monitor for use with the default PETSc
751:    SNES matrix free routines. Prints the differencing parameter used at 
752:    each step.
753: */
754: int MatSNESMFKSPMonitor(KSP ksp,int n,PetscReal rnorm,void *dummy)
755: {
756:   PC             pc;
757:   MatSNESMFCtx   ctx;
758:   int            ierr;
759:   Mat            mat;
760:   MPI_Comm       comm;
761:   PetscTruth     nonzeroinitialguess;

764:   PetscObjectGetComm((PetscObject)ksp,&comm);
765:   KSPGetPC(ksp,&pc);
766:   KSPGetInitialGuessNonzero(ksp,&nonzeroinitialguess);
767:   PCGetOperators(pc,&mat,PETSC_NULL,PETSC_NULL);
768:   ctx  = (MatSNESMFCtx)mat->data;

770:   if (n > 0 || nonzeroinitialguess) {
771: #if defined(PETSC_USE_COMPLEX)
772:     PetscPrintf(comm,"%d KSP Residual norm %14.12e h %g + %g in",n,rnorm,
773:                 PetscRealPart(ctx->currenth),PetscImaginaryPart(ctx->currenth));
774: #else
775:     PetscPrintf(comm,"%d KSP Residual norm %14.12e h %g n",n,rnorm,ctx->currenth);
776: #endif
777:   } else {
778:     PetscPrintf(comm,"%d KSP Residual norm %14.12en",n,rnorm);
779:   }
780:   return(0);
781: }

783: #undef __FUNCT__  
785: /*@C
786:    MatSNESMFSetFunction - Sets the function used in applying the matrix free.

788:    Collective on Mat

790:    Input Parameters:
791: +  mat - the matrix free matrix created via MatCreateSNESMF()
792: .  v   - workspace vector
793: .  func - the function to use
794: -  funcctx - optional function context passed to function

796:    Level: advanced

798:    Notes:
799:     If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
800:     matrix inside your compute Jacobian routine

802:     If this is not set then it will use the function set with SNESSetFunction()

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

806: .seealso: MatCreateSNESMF(),MatSNESMFGetH(),
807:           MatSNESMFSetHHistory(), MatSNESMFResetHHistory(),
808:           MatSNESMFKSPMonitor(), SNESetFunction()
809: @*/
810: int MatSNESMFSetFunction(Mat mat,Vec v,int (*func)(SNES,Vec,Vec,void *),void *funcctx)
811: {
812:   MatSNESMFCtx ctx = (MatSNESMFCtx)mat->data;

815:   ctx->func    = func;
816:   ctx->funcctx = funcctx;
817:   ctx->funcvec = v;
818:   return(0);
819: }

821: #undef __FUNCT__  
823: /*@C
824:    MatSNESMFSetFunctioni - Sets the function for a single component

826:    Collective on Mat

828:    Input Parameters:
829: +  mat - the matrix free matrix created via MatCreateSNESMF()
830: -  funci - the function to use

832:    Level: advanced

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


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

841: .seealso: MatCreateSNESMF(),MatSNESMFGetH(),
842:           MatSNESMFSetHHistory(), MatSNESMFResetHHistory(),
843:           MatSNESMFKSPMonitor(), SNESetFunction()
844: @*/
845: int MatSNESMFSetFunctioni(Mat mat,int (*funci)(int,Vec,PetscScalar*,void *))
846: {
847:   int  ierr,(*f)(Mat,int (*)(int,Vec,PetscScalar*,void *));

851:   PetscObjectQueryFunction((PetscObject)mat,"MatSNESMFSetFunctioni_C",(void (**)(void))&f);
852:   if (f) {
853:     (*f)(mat,funci);
854:   }
855:   return(0);
856: }


859: #undef __FUNCT__  
861: /*@C
862:    MatSNESMFSetFunctioniBase - Sets the base vector for a single component function evaluation

864:    Collective on Mat

866:    Input Parameters:
867: +  mat - the matrix free matrix created via MatCreateSNESMF()
868: -  func - the function to use

870:    Level: advanced

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


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

879: .seealso: MatCreateSNESMF(),MatSNESMFGetH(),
880:           MatSNESMFSetHHistory(), MatSNESMFResetHHistory(),
881:           MatSNESMFKSPMonitor(), SNESetFunction()
882: @*/
883: int MatSNESMFSetFunctioniBase(Mat mat,int (*func)(Vec,void *))
884: {
885:   int  ierr,(*f)(Mat,int (*)(Vec,void *));

889:   PetscObjectQueryFunction((PetscObject)mat,"MatSNESMFSetFunctioniBase_C",(void (**)(void))&f);
890:   if (f) {
891:     (*f)(mat,func);
892:   }
893:   return(0);
894: }


897: #undef __FUNCT__  
899: /*@
900:    MatSNESMFSetPeriod - Sets how often h is recomputed, by default it is everytime

902:    Collective on Mat

904:    Input Parameters:
905: +  mat - the matrix free matrix created via MatCreateSNESMF()
906: -  period - 1 for everytime, 2 for every second etc

908:    Options Database Keys:
909: +  -snes_mf_period <period>

911:    Level: advanced


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

916: .seealso: MatCreateSNESMF(),MatSNESMFGetH(),
917:           MatSNESMFSetHHistory(), MatSNESMFResetHHistory(),
918:           MatSNESMFKSPMonitor()
919: @*/
920: int MatSNESMFSetPeriod(Mat mat,int period)
921: {
922:   MatSNESMFCtx ctx = (MatSNESMFCtx)mat->data;

925:   ctx->recomputeperiod = period;
926:   return(0);
927: }

929: #undef __FUNCT__  
931: /*@
932:    MatSNESMFSetFunctionError - Sets the error_rel for the approximation of
933:    matrix-vector products using finite differences.

935:    Collective on Mat

937:    Input Parameters:
938: +  mat - the matrix free matrix created via MatCreateSNESMF()
939: -  error_rel - relative error (should be set to the square root of
940:                the relative error in the function evaluations)

942:    Options Database Keys:
943: +  -snes_mf_err <error_rel> - Sets error_rel

945:    Level: advanced

947:    Notes:
948:    The default matrix-free matrix-vector product routine computes
949: .vb
950:      F'(u)*a = [F(u+h*a) - F(u)]/h where
951:      h = error_rel*u'a/||a||^2                        if  |u'a| > umin*||a||_{1}
952:        = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2   else
953: .ve

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

957: .seealso: MatCreateSNESMF(),MatSNESMFGetH(),
958:           MatSNESMFSetHHistory(), MatSNESMFResetHHistory(),
959:           MatSNESMFKSPMonitor()
960: @*/
961: int MatSNESMFSetFunctionError(Mat mat,PetscReal error)
962: {
963:   MatSNESMFCtx ctx = (MatSNESMFCtx)mat->data;

966:   if (error != PETSC_DEFAULT) ctx->error_rel = error;
967:   return(0);
968: }

970: #undef __FUNCT__  
972: /*@
973:    MatSNESMFAddNullSpace - Provides a null space that an operator is
974:    supposed to have.  Since roundoff will create a small component in
975:    the null space, if you know the null space you may have it
976:    automatically removed.

978:    Collective on Mat 

980:    Input Parameters:
981: +  J - the matrix-free matrix context
982: -  nullsp - object created with MatNullSpaceCreate()

984:    Level: advanced

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

988: .seealso: MatNullSpaceCreate(), MatSNESMFGetH(), MatCreateSNESMF(),
989:           MatSNESMFSetHHistory(), MatSNESMFResetHHistory(),
990:           MatSNESMFKSPMonitor(), MatSNESMFErrorRel()
991: @*/
992: int MatSNESMFAddNullSpace(Mat J,MatNullSpace nullsp)
993: {
994:   int          ierr;
995:   MatSNESMFCtx ctx = (MatSNESMFCtx)J->data;
996:   MPI_Comm     comm;

999:   PetscObjectGetComm((PetscObject)J,&comm);

1001:   ctx->sp = nullsp;
1002:   ierr    = PetscObjectReference((PetscObject)nullsp);
1003:   return(0);
1004: }

1006: #undef __FUNCT__  
1008: /*@
1009:    MatSNESMFSetHHistory - Sets an array to collect a history of the
1010:    differencing values (h) computed for the matrix-free product.

1012:    Collective on Mat 

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

1020:    Level: advanced

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

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

1028: .seealso: MatSNESMFGetH(), MatCreateSNESMF(),
1029:           MatSNESMFResetHHistory(),
1030:           MatSNESMFKSPMonitor(), MatSNESMFSetFunctionError()

1032: @*/
1033: int MatSNESMFSetHHistory(Mat J,PetscScalar *history,int nhistory)
1034: {
1035:   MatSNESMFCtx ctx = (MatSNESMFCtx)J->data;

1038:   ctx->historyh    = history;
1039:   ctx->maxcurrenth = nhistory;
1040:   ctx->currenth    = 0;
1041:   return(0);
1042: }

1044: #undef __FUNCT__  
1046: /*@
1047:    MatSNESMFResetHHistory - Resets the counter to zero to begin 
1048:    collecting a new set of differencing histories.

1050:    Collective on Mat 

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

1055:    Level: advanced

1057:    Notes:
1058:    Use MatSNESMFSetHHistory() to create the original history counter.

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

1062: .seealso: MatSNESMFGetH(), MatCreateSNESMF(),
1063:           MatSNESMFSetHHistory(),
1064:           MatSNESMFKSPMonitor(), MatSNESMFSetFunctionError()

1066: @*/
1067: int MatSNESMFResetHHistory(Mat J)
1068: {
1069:   MatSNESMFCtx ctx = (MatSNESMFCtx)J->data;

1072:   ctx->ncurrenth    = 0;
1073:   return(0);
1074: }

1076: #undef __FUNCT__  
1078: int MatSNESMFComputeJacobian(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure *flag,void *dummy)
1079: {
1082:   MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
1083:   MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
1084:   return(0);
1085: }

1087: #undef __FUNCT__  
1089: /*@
1090:     MatSNESMFSetBase - Sets the vector U at which matrix vector products of the 
1091:         Jacobian are computed

1093:     Collective on Mat

1095:     Input Parameters:
1096: +   J - the MatSNESMF matrix
1097: -   U - the vector

1099:     Notes: This is rarely used directly

1101:     Level: advanced

1103: @*/
1104: int MatSNESMFSetBase(Mat J,Vec U)
1105: {
1106:   int  ierr,(*f)(Mat,Vec);

1111:   PetscObjectQueryFunction((PetscObject)J,"MatSNESMFSetBase_C",(void (**)(void))&f);
1112:   if (f) {
1113:     (*f)(J,U);
1114:   }
1115:   return(0);
1116: }

1118: #undef __FUNCT__  
1120: /*@C
1121:     MatSNESMFSetCheckh - Sets a function that checks the computed h and adjusts
1122:         it to satisfy some criteria

1124:     Collective on Mat

1126:     Input Parameters:
1127: +   J - the MatSNESMF matrix
1128: .   fun - the function that checks h
1129: -   ctx - any context needed by the function

1131:     Options Database Keys:
1132: .   -snes_mf_check_positivity

1134:     Level: advanced

1136:     Notes: For example, MatSNESMFSetCheckPositivity() insures that all entries
1137:        of U + h*a are non-negative

1139: .seealso:  MatSNESMFSetCheckPositivity()
1140: @*/
1141: int MatSNESMFSetCheckh(Mat J,int (*fun)(Vec,Vec,PetscScalar*,void*),void* ctx)
1142: {
1143:   int  ierr,(*f)(Mat,int (*)(Vec,Vec,PetscScalar*,void*),void*);

1147:   PetscObjectQueryFunction((PetscObject)J,"MatSNESMFSetCheckh_C",(void (**)(void))&f);
1148:   if (f) {
1149:     (*f)(J,fun,ctx);
1150:   }
1151:   return(0);
1152: }

1154: #undef __FUNCT__  
1156: /*@
1157:     MatSNESMFCheckPositivity - Checks that all entries in U + h*a are positive or
1158:         zero, decreases h until this is satisfied.

1160:     Collective on Vec

1162:     Input Parameters:
1163: +   U - base vector that is added to
1164: .   a - vector that is added
1165: .   h - scaling factor on a
1166: -   dummy - context variable (unused)

1168:     Options Database Keys:
1169: .   -snes_mf_check_positivity

1171:     Level: advanced

1173:     Notes: This is rarely used directly, rather it is passed as an argument to 
1174:            MatSNESMFSetCheckh()

1176: .seealso:  MatSNESMFSetCheckh()
1177: @*/
1178: int MatSNESMFCheckPositivity(Vec U,Vec a,PetscScalar *h,void *dummy)
1179: {
1180:   PetscReal     val, minval;
1181:   PetscScalar   *u_vec, *a_vec;
1182:   int           ierr, i, size;
1183:   MPI_Comm      comm;

1186:   PetscObjectGetComm((PetscObject)U,&comm);
1187:   VecGetArray(U,&u_vec);
1188:   VecGetArray(a,&a_vec);
1189:   VecGetLocalSize(U,&size);
1190:   minval = PetscAbsScalar(*h*1.01);
1191:   for(i=0;i<size;i++) {
1192:     if (PetscRealPart(u_vec[i] + *h*a_vec[i]) <= 0.0) {
1193:       val = PetscAbsScalar(u_vec[i]/a_vec[i]);
1194:       if (val < minval) minval = val;
1195:     }
1196:   }
1197:   VecRestoreArray(U,&u_vec);
1198:   VecRestoreArray(a,&a_vec);
1199:   PetscGlobalMin(&minval,&val,comm);
1200:   if (val <= PetscAbsScalar(*h)) {
1201:     PetscLogInfo(U,"MatSNESMFCheckPositivity: Scaling back h from %g to %gn",PetscRealPart(*h),.99*val);
1202:     if (PetscRealPart(*h) > 0.0) *h =  0.99*val;
1203:     else                         *h = -0.99*val;
1204:   }
1205:   return(0);
1206: }