Actual source code: snes.c

  1: /*$Id: snes.c,v 1.235 2001/08/21 21:03:49 bsmith Exp $*/

 3:  #include src/snes/snesimpl.h

  5: PetscTruth SNESRegisterAllCalled = PETSC_FALSE;
  6: PetscFList SNESList              = PETSC_NULL;

  8: /* Logging support */
  9: int SNES_COOKIE;
 10: int SNES_Solve, SNES_LineSearch, SNES_FunctionEval, SNES_JacobianEval;

 12: #undef __FUNCT__  
 14: /*@C
 15:    SNESView - Prints the SNES data structure.

 17:    Collective on SNES

 19:    Input Parameters:
 20: +  SNES - the SNES context
 21: -  viewer - visualization context

 23:    Options Database Key:
 24: .  -snes_view - Calls SNESView() at end of SNESSolve()

 26:    Notes:
 27:    The available visualization contexts include
 28: +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
 29: -     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
 30:          output where only the first processor opens
 31:          the file.  All other processors send their 
 32:          data to the first processor to print. 

 34:    The user can open an alternative visualization context with
 35:    PetscViewerASCIIOpen() - output to a specified file.

 37:    Level: beginner

 39: .keywords: SNES, view

 41: .seealso: PetscViewerASCIIOpen()
 42: @*/
 43: int SNESView(SNES snes,PetscViewer viewer)
 44: {
 45:   SNES_KSP_EW_ConvCtx *kctx;
 46:   int                 ierr;
 47:   SLES                sles;
 48:   char                *type;
 49:   PetscTruth          isascii,isstring;

 53:   if (!viewer) viewer = PETSC_VIEWER_STDOUT_(snes->comm);

 57:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);
 58:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_STRING,&isstring);
 59:   if (isascii) {
 60:     if (snes->prefix) {
 61:       PetscViewerASCIIPrintf(viewer,"SNES Object:(%s)n",snes->prefix);
 62:     } else {
 63:       PetscViewerASCIIPrintf(viewer,"SNES Object:n");
 64:     }
 65:     SNESGetType(snes,&type);
 66:     if (type) {
 67:       PetscViewerASCIIPrintf(viewer,"  type: %sn",type);
 68:     } else {
 69:       PetscViewerASCIIPrintf(viewer,"  type: not set yetn");
 70:     }
 71:     if (snes->view) {
 72:       PetscViewerASCIIPushTab(viewer);
 73:       (*snes->view)(snes,viewer);
 74:       PetscViewerASCIIPopTab(viewer);
 75:     }
 76:     PetscViewerASCIIPrintf(viewer,"  maximum iterations=%d, maximum function evaluations=%dn",snes->max_its,snes->max_funcs);
 77:     PetscViewerASCIIPrintf(viewer,"  tolerances: relative=%g, absolute=%g, solution=%gn",
 78:                  snes->rtol,snes->atol,snes->xtol);
 79:     PetscViewerASCIIPrintf(viewer,"  total number of linear solver iterations=%dn",snes->linear_its);
 80:     PetscViewerASCIIPrintf(viewer,"  total number of function evaluations=%dn",snes->nfuncs);
 81:     if (snes->ksp_ewconv) {
 82:       kctx = (SNES_KSP_EW_ConvCtx *)snes->kspconvctx;
 83:       if (kctx) {
 84:         PetscViewerASCIIPrintf(viewer,"  Eisenstat-Walker computation of KSP relative tolerance (version %d)n",kctx->version);
 85:         PetscViewerASCIIPrintf(viewer,"    rtol_0=%g, rtol_max=%g, threshold=%gn",kctx->rtol_0,kctx->rtol_max,kctx->threshold);
 86:         PetscViewerASCIIPrintf(viewer,"    gamma=%g, alpha=%g, alpha2=%gn",kctx->gamma,kctx->alpha,kctx->alpha2);
 87:       }
 88:     }
 89:   } else if (isstring) {
 90:     SNESGetType(snes,&type);
 91:     PetscViewerStringSPrintf(viewer," %-3.3s",type);
 92:   }
 93:   SNESGetSLES(snes,&sles);
 94:   PetscViewerASCIIPushTab(viewer);
 95:   SLESView(sles,viewer);
 96:   PetscViewerASCIIPopTab(viewer);
 97:   return(0);
 98: }

100: /*
101:   We retain a list of functions that also take SNES command 
102:   line options. These are called at the end SNESSetFromOptions()
103: */
104: #define MAXSETFROMOPTIONS 5
105: static int numberofsetfromoptions;
106: static int (*othersetfromoptions[MAXSETFROMOPTIONS])(SNES);

108: #undef __FUNCT__  
110: /*@
111:   SNESAddOptionsChecker - Adds an additional function to check for SNES options.

113:   Not Collective

115:   Input Parameter:
116: . snescheck - function that checks for options

118:   Level: developer

120: .seealso: SNESSetFromOptions()
121: @*/
122: int SNESAddOptionsChecker(int (*snescheck)(SNES))
123: {
125:   if (numberofsetfromoptions >= MAXSETFROMOPTIONS) {
126:     SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE, "Too many options checkers, only %d allowed", MAXSETFROMOPTIONS);
127:   }

129:   othersetfromoptions[numberofsetfromoptions++] = snescheck;
130:   return(0);
131: }

133: #undef __FUNCT__  
135: /*@
136:    SNESSetFromOptions - Sets various SNES and SLES parameters from user options.

138:    Collective on SNES

140:    Input Parameter:
141: .  snes - the SNES context

143:    Options Database Keys:
144: +  -snes_type <type> - ls, tr, umls, umtr, test
145: .  -snes_stol - convergence tolerance in terms of the norm
146:                 of the change in the solution between steps
147: .  -snes_atol <atol> - absolute tolerance of residual norm
148: .  -snes_rtol <rtol> - relative decrease in tolerance norm from initial
149: .  -snes_max_it <max_it> - maximum number of iterations
150: .  -snes_max_funcs <max_funcs> - maximum number of function evaluations
151: .  -snes_max_fail <max_fail> - maximum number of failures
152: .  -snes_trtol <trtol> - trust region tolerance
153: .  -snes_no_convergence_test - skip convergence test in nonlinear or minimization 
154:                                solver; hence iterations will continue until max_it
155:                                or some other criterion is reached. Saves expense
156:                                of convergence test
157: .  -snes_monitor - prints residual norm at each iteration 
158: .  -snes_vecmonitor - plots solution at each iteration
159: .  -snes_vecmonitor_update - plots update to solution at each iteration 
160: .  -snes_xmonitor - plots residual norm at each iteration 
161: .  -snes_fd - use finite differences to compute Jacobian; very slow, only for testing
162: -  -snes_mf_ksp_monitor - if using matrix-free multiply then print h at each KSP iteration

164:     Options Database for Eisenstat-Walker method:
165: +  -snes_ksp_ew_conv - use Eisenstat-Walker method for determining linear system convergence
166: .  -snes_ksp_ew_version ver - version of  Eisenstat-Walker method
167: .  -snes_ksp_ew_rtol0 <rtol0> - Sets rtol0
168: .  -snes_ksp_ew_rtolmax <rtolmax> - Sets rtolmax
169: .  -snes_ksp_ew_gamma <gamma> - Sets gamma
170: .  -snes_ksp_ew_alpha <alpha> - Sets alpha
171: .  -snes_ksp_ew_alpha2 <alpha2> - Sets alpha2 
172: -  -snes_ksp_ew_threshold <threshold> - Sets threshold

174:    Notes:
175:    To see all options, run your program with the -help option or consult
176:    the users manual.

178:    Level: beginner

180: .keywords: SNES, nonlinear, set, options, database

182: .seealso: SNESSetOptionsPrefix()
183: @*/
184: int SNESSetFromOptions(SNES snes)
185: {
186:   SLES                sles;
187:   SNES_KSP_EW_ConvCtx *kctx = (SNES_KSP_EW_ConvCtx *)snes->kspconvctx;
188:   PetscTruth          flg;
189:   int                 ierr, i;
190:   char                *deft,type[256];


195:   PetscOptionsBegin(snes->comm,snes->prefix,"Nonlinear solver (SNES) options","SNES");
196:     if (snes->type_name) {
197:       deft = snes->type_name;
198:     } else {
199:       deft = SNESLS;
200:     }

202:     if (!SNESRegisterAllCalled) {SNESRegisterAll(PETSC_NULL);}
203:     PetscOptionsList("-snes_type","Nonlinear solver method","SNESSetType",SNESList,deft,type,256,&flg);
204:     if (flg) {
205:       SNESSetType(snes,type);
206:     } else if (!snes->type_name) {
207:       SNESSetType(snes,deft);
208:     }
209:     PetscOptionsName("-snes_view","Print detailed information on solver used","SNESView",0);

211:     PetscOptionsReal("-snes_stol","Stop if step length less then","SNESSetTolerances",snes->xtol,&snes->xtol,0);
212:     PetscOptionsReal("-snes_atol","Stop if function norm less then","SNESSetTolerances",snes->atol,&snes->atol,0);

214:     PetscOptionsReal("-snes_rtol","Stop if decrease in function norm less then","SNESSetTolerances",snes->rtol,&snes->rtol,0);
215:     PetscOptionsInt("-snes_max_it","Maximum iterations","SNESSetTolerances",snes->max_its,&snes->max_its,PETSC_NULL);
216:     PetscOptionsInt("-snes_max_funcs","Maximum function evaluations","SNESSetTolerances",snes->max_funcs,&snes->max_funcs,PETSC_NULL);
217:     PetscOptionsInt("-snes_max_fail","Maximum failures","SNESSetTolerances",snes->maxFailures,&snes->maxFailures,PETSC_NULL);

219:     PetscOptionsName("-snes_ksp_ew_conv","Use Eisentat-Walker linear system convergence test","SNES_KSP_SetParametersEW",&snes->ksp_ewconv);

221:     PetscOptionsInt("-snes_ksp_ew_version","Version 1 or 2","SNES_KSP_SetParametersEW",kctx->version,&kctx->version,0);
222:     PetscOptionsReal("-snes_ksp_ew_rtol0","0 <= rtol0 < 1","SNES_KSP_SetParametersEW",kctx->rtol_0,&kctx->rtol_0,0);
223:     PetscOptionsReal("-snes_ksp_ew_rtolmax","0 <= rtolmax < 1","SNES_KSP_SetParametersEW",kctx->rtol_max,&kctx->rtol_max,0);
224:     PetscOptionsReal("-snes_ksp_ew_gamma","0 <= gamma <= 1","SNES_KSP_SetParametersEW",kctx->gamma,&kctx->gamma,0);
225:     PetscOptionsReal("-snes_ksp_ew_alpha","1 < alpha <= 2","SNES_KSP_SetParametersEW",kctx->alpha,&kctx->alpha,0);
226:     PetscOptionsReal("-snes_ksp_ew_alpha2","alpha2","SNES_KSP_SetParametersEW",kctx->alpha2,&kctx->alpha2,0);
227:     PetscOptionsReal("-snes_ksp_ew_threshold","0 < threshold < 1","SNES_KSP_SetParametersEW",kctx->threshold,&kctx->threshold,0);

229:     PetscOptionsName("-snes_no_convergence_test","Don't test for convergence","None",&flg);
230:     if (flg) {snes->converged = 0;}
231:     PetscOptionsName("-snes_cancelmonitors","Remove all monitors","SNESClearMonitor",&flg);
232:     if (flg) {SNESClearMonitor(snes);}
233:     PetscOptionsName("-snes_monitor","Monitor norm of function","SNESDefaultMonitor",&flg);
234:     if (flg) {SNESSetMonitor(snes,SNESDefaultMonitor,0,0);}
235:     PetscOptionsName("-snes_ratiomonitor","Monitor norm of function","SNESSetRatioMonitor",&flg);
236:     if (flg) {SNESSetRatioMonitor(snes);}
237:     PetscOptionsName("-snes_smonitor","Monitor norm of function (fewer digits)","SNESDefaultSMonitor",&flg);
238:     if (flg) {SNESSetMonitor(snes,SNESDefaultSMonitor,0,0);}
239:     PetscOptionsName("-snes_vecmonitor","Plot solution at each iteration","SNESVecViewMonitor",&flg);
240:     if (flg) {SNESSetMonitor(snes,SNESVecViewMonitor,0,0);}
241:     PetscOptionsName("-snes_vecmonitor_update","Plot correction at each iteration","SNESVecViewUpdateMonitor",&flg);
242:     if (flg) {SNESSetMonitor(snes,SNESVecViewUpdateMonitor,0,0);}
243:     PetscOptionsName("-snes_vecmonitor_residual","Plot residual at each iteration","SNESVecViewResidualMonitor",&flg);
244:     if (flg) {SNESSetMonitor(snes,SNESVecViewResidualMonitor,0,0);}
245:     PetscOptionsName("-snes_xmonitor","Plot function norm at each iteration","SNESLGMonitor",&flg);
246:     if (flg) {SNESSetMonitor(snes,SNESLGMonitor,PETSC_NULL,PETSC_NULL);}

248:     PetscOptionsName("-snes_fd","Use finite differences (slow) to compute Jacobian","SNESDefaultComputeJacobian",&flg);
249:     if (flg) {
250:       SNESSetJacobian(snes,snes->jacobian,snes->jacobian_pre,SNESDefaultComputeJacobian,snes->funP);
251:       PetscLogInfo(snes,"SNESSetFromOptions: Setting default finite difference Jacobian matrixn");
252:     }

254:     for(i = 0; i < numberofsetfromoptions; i++) {
255:       (*othersetfromoptions[i])(snes);
256:     }

258:     if (snes->setfromoptions) {
259:       (*snes->setfromoptions)(snes);
260:     }

262:   PetscOptionsEnd();

264:   SNESGetSLES(snes,&sles);
265:   SLESSetFromOptions(sles);

267:   return(0);
268: }


271: #undef __FUNCT__  
273: /*@
274:    SNESSetApplicationContext - Sets the optional user-defined context for 
275:    the nonlinear solvers.  

277:    Collective on SNES

279:    Input Parameters:
280: +  snes - the SNES context
281: -  usrP - optional user context

283:    Level: intermediate

285: .keywords: SNES, nonlinear, set, application, context

287: .seealso: SNESGetApplicationContext()
288: @*/
289: int SNESSetApplicationContext(SNES snes,void *usrP)
290: {
293:   snes->user                = usrP;
294:   return(0);
295: }

297: #undef __FUNCT__  
299: /*@C
300:    SNESGetApplicationContext - Gets the user-defined context for the 
301:    nonlinear solvers.  

303:    Not Collective

305:    Input Parameter:
306: .  snes - SNES context

308:    Output Parameter:
309: .  usrP - user context

311:    Level: intermediate

313: .keywords: SNES, nonlinear, get, application, context

315: .seealso: SNESSetApplicationContext()
316: @*/
317: int SNESGetApplicationContext(SNES snes,void **usrP)
318: {
321:   *usrP = snes->user;
322:   return(0);
323: }

325: #undef __FUNCT__  
327: /*@
328:    SNESGetIterationNumber - Gets the number of nonlinear iterations completed
329:    at this time.

331:    Not Collective

333:    Input Parameter:
334: .  snes - SNES context

336:    Output Parameter:
337: .  iter - iteration number

339:    Notes:
340:    For example, during the computation of iteration 2 this would return 1.

342:    This is useful for using lagged Jacobians (where one does not recompute the 
343:    Jacobian at each SNES iteration). For example, the code
344: .vb
345:       SNESGetIterationNumber(snes,&it);
346:       if (!(it % 2)) {
347:         [compute Jacobian here]
348:       }
349: .ve
350:    can be used in your ComputeJacobian() function to cause the Jacobian to be
351:    recomputed every second SNES iteration.

353:    Level: intermediate

355: .keywords: SNES, nonlinear, get, iteration, number
356: @*/
357: int SNESGetIterationNumber(SNES snes,int* iter)
358: {
362:   *iter = snes->iter;
363:   return(0);
364: }

366: #undef __FUNCT__  
368: /*@
369:    SNESGetFunctionNorm - Gets the norm of the current function that was set
370:    with SNESSSetFunction().

372:    Collective on SNES

374:    Input Parameter:
375: .  snes - SNES context

377:    Output Parameter:
378: .  fnorm - 2-norm of function

380:    Level: intermediate

382: .keywords: SNES, nonlinear, get, function, norm

384: .seealso: SNESGetFunction()
385: @*/
386: int SNESGetFunctionNorm(SNES snes,PetscScalar *fnorm)
387: {
391:   *fnorm = snes->norm;
392:   return(0);
393: }

395: #undef __FUNCT__  
397: /*@
398:    SNESGetNumberUnsuccessfulSteps - Gets the number of unsuccessful steps
399:    attempted by the nonlinear solver.

401:    Not Collective

403:    Input Parameter:
404: .  snes - SNES context

406:    Output Parameter:
407: .  nfails - number of unsuccessful steps attempted

409:    Notes:
410:    This counter is reset to zero for each successive call to SNESSolve().

412:    Level: intermediate

414: .keywords: SNES, nonlinear, get, number, unsuccessful, steps
415: @*/
416: int SNESGetNumberUnsuccessfulSteps(SNES snes,int* nfails)
417: {
421:   *nfails = snes->numFailures;
422:   return(0);
423: }

425: #undef __FUNCT__  
427: /*@
428:    SNESSetMaximumUnsuccessfulSteps - Sets the maximum number of unsuccessful steps
429:    attempted by the nonlinear solver before it gives up.

431:    Not Collective

433:    Input Parameters:
434: +  snes     - SNES context
435: -  maxFails - maximum of unsuccessful steps

437:    Level: intermediate

439: .keywords: SNES, nonlinear, set, maximum, unsuccessful, steps
440: @*/
441: int SNESSetMaximumUnsuccessfulSteps(SNES snes, int maxFails)
442: {
445:   snes->maxFailures = maxFails;
446:   return(0);
447: }

449: #undef __FUNCT__  
451: /*@
452:    SNESGetMaximumUnsuccessfulSteps - Gets the maximum number of unsuccessful steps
453:    attempted by the nonlinear solver before it gives up.

455:    Not Collective

457:    Input Parameter:
458: .  snes     - SNES context

460:    Output Parameter:
461: .  maxFails - maximum of unsuccessful steps

463:    Level: intermediate

465: .keywords: SNES, nonlinear, get, maximum, unsuccessful, steps
466: @*/
467: int SNESGetMaximumUnsuccessfulSteps(SNES snes, int *maxFails)
468: {
472:   *maxFails = snes->maxFailures;
473:   return(0);
474: }

476: #undef __FUNCT__  
478: /*@
479:    SNESGetNumberLinearIterations - Gets the total number of linear iterations
480:    used by the nonlinear solver.

482:    Not Collective

484:    Input Parameter:
485: .  snes - SNES context

487:    Output Parameter:
488: .  lits - number of linear iterations

490:    Notes:
491:    This counter is reset to zero for each successive call to SNESSolve().

493:    Level: intermediate

495: .keywords: SNES, nonlinear, get, number, linear, iterations
496: @*/
497: int SNESGetNumberLinearIterations(SNES snes,int* lits)
498: {
502:   *lits = snes->linear_its;
503:   return(0);
504: }

506: #undef __FUNCT__  
508: /*@C
509:    SNESGetSLES - Returns the SLES context for a SNES solver.

511:    Not Collective, but if SNES object is parallel, then SLES object is parallel

513:    Input Parameter:
514: .  snes - the SNES context

516:    Output Parameter:
517: .  sles - the SLES context

519:    Notes:
520:    The user can then directly manipulate the SLES context to set various
521:    options, etc.  Likewise, the user can then extract and manipulate the 
522:    KSP and PC contexts as well.

524:    Level: beginner

526: .keywords: SNES, nonlinear, get, SLES, context

528: .seealso: SLESGetPC(), SLESGetKSP()
529: @*/
530: int SNESGetSLES(SNES snes,SLES *sles)
531: {
534:   *sles = snes->sles;
535:   return(0);
536: }

538: #undef __FUNCT__  
540: static int SNESPublish_Petsc(PetscObject obj)
541: {
542: #if defined(PETSC_HAVE_AMS)
543:   SNES          v = (SNES) obj;
544:   int          ierr;
545: #endif


549: #if defined(PETSC_HAVE_AMS)
550:   /* if it is already published then return */
551:   if (v->amem >=0) return(0);

553:   PetscObjectPublishBaseBegin(obj);
554:   AMS_Memory_add_field((AMS_Memory)v->amem,"Iteration",&v->iter,1,AMS_INT,AMS_READ,
555:                                 AMS_COMMON,AMS_REDUCT_UNDEF);
556:   AMS_Memory_add_field((AMS_Memory)v->amem,"Residual",&v->norm,1,AMS_DOUBLE,AMS_READ,
557:                                 AMS_COMMON,AMS_REDUCT_UNDEF);
558:   PetscObjectPublishBaseEnd(obj);
559: #endif
560:   return(0);
561: }

563: /* -----------------------------------------------------------*/
564: #undef __FUNCT__  
566: /*@C
567:    SNESCreate - Creates a nonlinear solver context.

569:    Collective on MPI_Comm

571:    Input Parameters:
572: +  comm - MPI communicator

574:    Output Parameter:
575: .  outsnes - the new SNES context

577:    Options Database Keys:
578: +   -snes_mf - Activates default matrix-free Jacobian-vector products,
579:                and no preconditioning matrix
580: .   -snes_mf_operator - Activates default matrix-free Jacobian-vector
581:                products, and a user-provided preconditioning matrix
582:                as set by SNESSetJacobian()
583: -   -snes_fd - Uses (slow!) finite differences to compute Jacobian

585:    Level: beginner

587: .keywords: SNES, nonlinear, create, context

589: .seealso: SNESSolve(), SNESDestroy(), SNES
590: @*/
591: int SNESCreate(MPI_Comm comm,SNES *outsnes)
592: {
593:   int                 ierr;
594:   SNES                snes;
595:   SNES_KSP_EW_ConvCtx *kctx;

599:   *outsnes = PETSC_NULL;
600: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
601:   SNESInitializePackage(PETSC_NULL);
602: #endif

604:   PetscHeaderCreate(snes,_p_SNES,int,SNES_COOKIE,0,"SNES",comm,SNESDestroy,SNESView);
605:   PetscLogObjectCreate(snes);
606:   snes->bops->publish     = SNESPublish_Petsc;
607:   snes->max_its           = 50;
608:   snes->max_funcs          = 10000;
609:   snes->norm                  = 0.0;
610:   snes->rtol                  = 1.e-8;
611:   snes->ttol              = 0.0;
612:   snes->atol                  = 1.e-50;
613:   snes->xtol                  = 1.e-8;
614:   snes->deltatol          = 1.e-12;
615:   snes->nfuncs            = 0;
616:   snes->numFailures       = 0;
617:   snes->maxFailures       = 1;
618:   snes->linear_its        = 0;
619:   snes->numbermonitors    = 0;
620:   snes->data              = 0;
621:   snes->view              = 0;
622:   snes->setupcalled       = 0;
623:   snes->ksp_ewconv        = PETSC_FALSE;
624:   snes->vwork             = 0;
625:   snes->nwork             = 0;
626:   snes->conv_hist_len     = 0;
627:   snes->conv_hist_max     = 0;
628:   snes->conv_hist         = PETSC_NULL;
629:   snes->conv_hist_its     = PETSC_NULL;
630:   snes->conv_hist_reset   = PETSC_TRUE;
631:   snes->reason            = SNES_CONVERGED_ITERATING;

633:   /* Create context to compute Eisenstat-Walker relative tolerance for KSP */
634:   PetscNew(SNES_KSP_EW_ConvCtx,&kctx);
635:   PetscLogObjectMemory(snes,sizeof(SNES_KSP_EW_ConvCtx));
636:   snes->kspconvctx  = (void*)kctx;
637:   kctx->version     = 2;
638:   kctx->rtol_0      = .3; /* Eisenstat and Walker suggest rtol_0=.5, but 
639:                              this was too large for some test cases */
640:   kctx->rtol_last   = 0;
641:   kctx->rtol_max    = .9;
642:   kctx->gamma       = 1.0;
643:   kctx->alpha2      = .5*(1.0 + sqrt(5.0));
644:   kctx->alpha       = kctx->alpha2;
645:   kctx->threshold   = .1;
646:   kctx->lresid_last = 0;
647:   kctx->norm_last   = 0;

649:   SLESCreate(comm,&snes->sles);
650:   PetscLogObjectParent(snes,snes->sles)

652:   *outsnes = snes;
653:   PetscPublishAll(snes);
654:   return(0);
655: }

657: /* --------------------------------------------------------------- */
658: #undef __FUNCT__  
660: /*@C
661:    SNESSetFunction - Sets the function evaluation routine and function 
662:    vector for use by the SNES routines in solving systems of nonlinear
663:    equations.

665:    Collective on SNES

667:    Input Parameters:
668: +  snes - the SNES context
669: .  func - function evaluation routine
670: .  r - vector to store function value
671: -  ctx - [optional] user-defined context for private data for the 
672:          function evaluation routine (may be PETSC_NULL)

674:    Calling sequence of func:
675: $    func (SNES snes,Vec x,Vec f,void *ctx);

677: .  f - function vector
678: -  ctx - optional user-defined function context 

680:    Notes:
681:    The Newton-like methods typically solve linear systems of the form
682: $      f'(x) x = -f(x),
683:    where f'(x) denotes the Jacobian matrix and f(x) is the function.

685:    Level: beginner

687: .keywords: SNES, nonlinear, set, function

689: .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian()
690: @*/
691: int SNESSetFunction(SNES snes,Vec r,int (*func)(SNES,Vec,Vec,void*),void *ctx)
692: {

698:   snes->computefunction     = func;
699:   snes->vec_func            = snes->vec_func_always = r;
700:   snes->funP                = ctx;
701:   return(0);
702: }

704: #undef __FUNCT__  
706: /*@
707:    SNESComputeFunction - Calls the function that has been set with
708:                          SNESSetFunction().  

710:    Collective on SNES

712:    Input Parameters:
713: +  snes - the SNES context
714: -  x - input vector

716:    Output Parameter:
717: .  y - function vector, as set by SNESSetFunction()

719:    Notes:
720:    SNESComputeFunction() is typically used within nonlinear solvers
721:    implementations, so most users would not generally call this routine
722:    themselves.

724:    Level: developer

726: .keywords: SNES, nonlinear, compute, function

728: .seealso: SNESSetFunction(), SNESGetFunction()
729: @*/
730: int SNESComputeFunction(SNES snes,Vec x,Vec y)
731: {
732:   int    ierr;


741:   PetscLogEventBegin(SNES_FunctionEval,snes,x,y,0);
742:   PetscStackPush("SNES user function");
743:   (*snes->computefunction)(snes,x,y,snes->funP);
744:   PetscStackPop;
745:   snes->nfuncs++;
746:   PetscLogEventEnd(SNES_FunctionEval,snes,x,y,0);
747:   return(0);
748: }

750: #undef __FUNCT__  
752: /*@
753:    SNESComputeJacobian - Computes the Jacobian matrix that has been
754:    set with SNESSetJacobian().

756:    Collective on SNES and Mat

758:    Input Parameters:
759: +  snes - the SNES context
760: -  x - input vector

762:    Output Parameters:
763: +  A - Jacobian matrix
764: .  B - optional preconditioning matrix
765: -  flag - flag indicating matrix structure

767:    Notes: 
768:    Most users should not need to explicitly call this routine, as it 
769:    is used internally within the nonlinear solvers. 

771:    See SLESSetOperators() for important information about setting the
772:    flag parameter.

774:    Level: developer

776: .keywords: SNES, compute, Jacobian, matrix

778: .seealso:  SNESSetJacobian(), SLESSetOperators()
779: @*/
780: int SNESComputeJacobian(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *flg)
781: {
782:   int    ierr;

788:   if (!snes->computejacobian) return(0);
789:   PetscLogEventBegin(SNES_JacobianEval,snes,X,*A,*B);
790:   *flg = DIFFERENT_NONZERO_PATTERN;
791:   PetscStackPush("SNES user Jacobian function");
792:   (*snes->computejacobian)(snes,X,A,B,flg,snes->jacP);
793:   PetscStackPop;
794:   PetscLogEventEnd(SNES_JacobianEval,snes,X,*A,*B);
795:   /* make sure user returned a correct Jacobian and preconditioner */
798:   return(0);
799: }

801: #undef __FUNCT__  
803: /*@C
804:    SNESSetJacobian - Sets the function to compute Jacobian as well as the
805:    location to store the matrix.

807:    Collective on SNES and Mat

809:    Input Parameters:
810: +  snes - the SNES context
811: .  A - Jacobian matrix
812: .  B - preconditioner matrix (usually same as the Jacobian)
813: .  func - Jacobian evaluation routine
814: -  ctx - [optional] user-defined context for private data for the 
815:          Jacobian evaluation routine (may be PETSC_NULL)

817:    Calling sequence of func:
818: $     func (SNES snes,Vec x,Mat *A,Mat *B,int *flag,void *ctx);

820: +  x - input vector
821: .  A - Jacobian matrix
822: .  B - preconditioner matrix, usually the same as A
823: .  flag - flag indicating information about the preconditioner matrix
824:    structure (same as flag in SLESSetOperators())
825: -  ctx - [optional] user-defined Jacobian context

827:    Notes: 
828:    See SLESSetOperators() for important information about setting the flag
829:    output parameter in the routine func().  Be sure to read this information!

831:    The routine func() takes Mat * as the matrix arguments rather than Mat.  
832:    This allows the Jacobian evaluation routine to replace A and/or B with a 
833:    completely new new matrix structure (not just different matrix elements)
834:    when appropriate, for instance, if the nonzero structure is changing
835:    throughout the global iterations.

837:    Level: beginner

839: .keywords: SNES, nonlinear, set, Jacobian, matrix

841: .seealso: SLESSetOperators(), SNESSetFunction()
842: @*/
843: int SNESSetJacobian(SNES snes,Mat A,Mat B,int (*func)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void *ctx)
844: {

853:   if (func) snes->computejacobian = func;
854:   if (ctx)  snes->jacP            = ctx;
855:   if (A) {
856:     if (snes->jacobian) {MatDestroy(snes->jacobian);}
857:     snes->jacobian = A;
858:     ierr           = PetscObjectReference((PetscObject)A);
859:   }
860:   if (B) {
861:     if (snes->jacobian_pre) {MatDestroy(snes->jacobian_pre);}
862:     snes->jacobian_pre = B;
863:     ierr               = PetscObjectReference((PetscObject)B);
864:   }
865:   return(0);
866: }

868: #undef __FUNCT__  
870: /*@C
871:    SNESGetJacobian - Returns the Jacobian matrix and optionally the user 
872:    provided context for evaluating the Jacobian.

874:    Not Collective, but Mat object will be parallel if SNES object is

876:    Input Parameter:
877: .  snes - the nonlinear solver context

879:    Output Parameters:
880: +  A - location to stash Jacobian matrix (or PETSC_NULL)
881: .  B - location to stash preconditioner matrix (or PETSC_NULL)
882: .  ctx - location to stash Jacobian ctx (or PETSC_NULL)
883: -  func - location to put Jacobian function (or PETSC_NULL)

885:    Level: advanced

887: .seealso: SNESSetJacobian(), SNESComputeJacobian()
888: @*/
889: int SNESGetJacobian(SNES snes,Mat *A,Mat *B,void **ctx,int (**func)(SNES,Vec,Mat*,Mat*,MatStructure*,void*))
890: {
893:   if (A)    *A    = snes->jacobian;
894:   if (B)    *B    = snes->jacobian_pre;
895:   if (ctx)  *ctx  = snes->jacP;
896:   if (func) *func = snes->computejacobian;
897:   return(0);
898: }

900: /* ----- Routines to initialize and destroy a nonlinear solver ---- */
901: extern int SNESDefaultMatrixFreeCreate2(SNES,Vec,Mat*);

903: #undef __FUNCT__  
905: /*@
906:    SNESSetUp - Sets up the internal data structures for the later use
907:    of a nonlinear solver.

909:    Collective on SNES

911:    Input Parameters:
912: +  snes - the SNES context
913: -  x - the solution vector

915:    Notes:
916:    For basic use of the SNES solvers the user need not explicitly call
917:    SNESSetUp(), since these actions will automatically occur during
918:    the call to SNESSolve().  However, if one wishes to control this
919:    phase separately, SNESSetUp() should be called after SNESCreate()
920:    and optional routines of the form SNESSetXXX(), but before SNESSolve().  

922:    Level: advanced

924: .keywords: SNES, nonlinear, setup

926: .seealso: SNESCreate(), SNESSolve(), SNESDestroy()
927: @*/
928: int SNESSetUp(SNES snes,Vec x)
929: {
930:   int        ierr;
931:   PetscTruth flg, iseqtr;

937:   snes->vec_sol = snes->vec_sol_always = x;

939:   PetscOptionsHasName(snes->prefix,"-snes_mf_operator",&flg);
940:   /*
941:       This version replaces the user provided Jacobian matrix with a
942:       matrix-free version but still employs the user-provided preconditioner matrix
943:   */
944:   if (flg) {
945:     Mat J;
946:     MatCreateSNESMF(snes,snes->vec_sol,&J);
947:     MatSNESMFSetFromOptions(J);
948:     PetscLogInfo(snes,"SNESSetUp: Setting default matrix-free operator routinesn");
949:     SNESSetJacobian(snes,J,0,0,0);
950:     MatDestroy(J);
951:   }

953: #if !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE)
954:   PetscOptionsHasName(snes->prefix,"-snes_mf_operator2",&flg);
955:   if (flg) {
956:     Mat J;
957:     SNESDefaultMatrixFreeCreate2(snes,snes->vec_sol,&J);
958:     SNESSetJacobian(snes,J,0,0,0);
959:     MatDestroy(J);
960:   }
961: #endif

963:   PetscOptionsHasName(snes->prefix,"-snes_mf",&flg);
964:   /*
965:       This version replaces both the user-provided Jacobian and the user-
966:       provided preconditioner matrix with the default matrix free version.
967:    */
968:   if (flg) {
969:     Mat  J;
970:     SLES sles;
971:     PC   pc;

973:     MatCreateSNESMF(snes,snes->vec_sol,&J);
974:     MatSNESMFSetFromOptions(J);
975:     PetscLogInfo(snes,"SNESSetUp: Setting default matrix-free operator and preconditioner routinesn");
976:     SNESSetJacobian(snes,J,J,MatSNESMFComputeJacobian,snes->funP);
977:     MatDestroy(J);

979:     /* force no preconditioner */
980:     SNESGetSLES(snes,&sles);
981:     SLESGetPC(sles,&pc);
982:     PetscTypeCompare((PetscObject)pc,PCSHELL,&flg);
983:     if (!flg) {
984:       PCSetType(pc,PCNONE);
985:     }
986:   }

988:   if (!snes->vec_func) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Must call SNESSetFunction() first");
989:   if (!snes->computefunction) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Must call SNESSetFunction() first");
990:   if (!snes->jacobian) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Must call SNESSetJacobian() first n or use -snes_mf option");
991:   if (snes->vec_func == snes->vec_sol) {
992:     SETERRQ(PETSC_ERR_ARG_IDN,"Solution vector cannot be function vector");
993:   }

995:   /* Set the KSP stopping criterion to use the Eisenstat-Walker method */
996:   PetscTypeCompare((PetscObject)snes,SNESTR,&iseqtr);
997:   if (snes->ksp_ewconv && !iseqtr) {
998:     SLES sles; KSP ksp;
999:     SNESGetSLES(snes,&sles);
1000:     SLESGetKSP(sles,&ksp);
1001:     KSPSetConvergenceTest(ksp,SNES_KSP_EW_Converged_Private,snes);
1002:   }

1004:   if (snes->setup) {(*snes->setup)(snes);}
1005:   snes->setupcalled = 1;
1006:   return(0);
1007: }

1009: #undef __FUNCT__  
1011: /*@C
1012:    SNESDestroy - Destroys the nonlinear solver context that was created
1013:    with SNESCreate().

1015:    Collective on SNES

1017:    Input Parameter:
1018: .  snes - the SNES context

1020:    Level: beginner

1022: .keywords: SNES, nonlinear, destroy

1024: .seealso: SNESCreate(), SNESSolve()
1025: @*/
1026: int SNESDestroy(SNES snes)
1027: {
1028:   int i,ierr;

1032:   if (--snes->refct > 0) return(0);

1034:   /* if memory was published with AMS then destroy it */
1035:   PetscObjectDepublish(snes);

1037:   if (snes->destroy) {(*(snes)->destroy)(snes);}
1038:   if (snes->kspconvctx) {PetscFree(snes->kspconvctx);}
1039:   if (snes->jacobian) {MatDestroy(snes->jacobian);}
1040:   if (snes->jacobian_pre) {MatDestroy(snes->jacobian_pre);}
1041:   SLESDestroy(snes->sles);
1042:   if (snes->vwork) {VecDestroyVecs(snes->vwork,snes->nvwork);}
1043:   for (i=0; i<snes->numbermonitors; i++) {
1044:     if (snes->monitordestroy[i]) {
1045:       (*snes->monitordestroy[i])(snes->monitorcontext[i]);
1046:     }
1047:   }
1048:   PetscLogObjectDestroy((PetscObject)snes);
1049:   PetscHeaderDestroy((PetscObject)snes);
1050:   return(0);
1051: }

1053: /* ----------- Routines to set solver parameters ---------- */

1055: #undef __FUNCT__  
1057: /*@
1058:    SNESSetTolerances - Sets various parameters used in convergence tests.

1060:    Collective on SNES

1062:    Input Parameters:
1063: +  snes - the SNES context
1064: .  atol - absolute convergence tolerance
1065: .  rtol - relative convergence tolerance
1066: .  stol -  convergence tolerance in terms of the norm
1067:            of the change in the solution between steps
1068: .  maxit - maximum number of iterations
1069: -  maxf - maximum number of function evaluations

1071:    Options Database Keys: 
1072: +    -snes_atol <atol> - Sets atol
1073: .    -snes_rtol <rtol> - Sets rtol
1074: .    -snes_stol <stol> - Sets stol
1075: .    -snes_max_it <maxit> - Sets maxit
1076: -    -snes_max_funcs <maxf> - Sets maxf

1078:    Notes:
1079:    The default maximum number of iterations is 50.
1080:    The default maximum number of function evaluations is 1000.

1082:    Level: intermediate

1084: .keywords: SNES, nonlinear, set, convergence, tolerances

1086: .seealso: SNESSetTrustRegionTolerance(), SNESSetMinimizationFunctionTolerance()
1087: @*/
1088: int SNESSetTolerances(SNES snes,PetscReal atol,PetscReal rtol,PetscReal stol,int maxit,int maxf)
1089: {
1092:   if (atol != PETSC_DEFAULT)  snes->atol      = atol;
1093:   if (rtol != PETSC_DEFAULT)  snes->rtol      = rtol;
1094:   if (stol != PETSC_DEFAULT)  snes->xtol      = stol;
1095:   if (maxit != PETSC_DEFAULT) snes->max_its   = maxit;
1096:   if (maxf != PETSC_DEFAULT)  snes->max_funcs = maxf;
1097:   return(0);
1098: }

1100: #undef __FUNCT__  
1102: /*@
1103:    SNESGetTolerances - Gets various parameters used in convergence tests.

1105:    Not Collective

1107:    Input Parameters:
1108: +  snes - the SNES context
1109: .  atol - absolute convergence tolerance
1110: .  rtol - relative convergence tolerance
1111: .  stol -  convergence tolerance in terms of the norm
1112:            of the change in the solution between steps
1113: .  maxit - maximum number of iterations
1114: -  maxf - maximum number of function evaluations

1116:    Notes:
1117:    The user can specify PETSC_NULL for any parameter that is not needed.

1119:    Level: intermediate

1121: .keywords: SNES, nonlinear, get, convergence, tolerances

1123: .seealso: SNESSetTolerances()
1124: @*/
1125: int SNESGetTolerances(SNES snes,PetscReal *atol,PetscReal *rtol,PetscReal *stol,int *maxit,int *maxf)
1126: {
1129:   if (atol)  *atol  = snes->atol;
1130:   if (rtol)  *rtol  = snes->rtol;
1131:   if (stol)  *stol  = snes->xtol;
1132:   if (maxit) *maxit = snes->max_its;
1133:   if (maxf)  *maxf  = snes->max_funcs;
1134:   return(0);
1135: }

1137: #undef __FUNCT__  
1139: /*@
1140:    SNESSetTrustRegionTolerance - Sets the trust region parameter tolerance.  

1142:    Collective on SNES

1144:    Input Parameters:
1145: +  snes - the SNES context
1146: -  tol - tolerance
1147:    
1148:    Options Database Key: 
1149: .  -snes_trtol <tol> - Sets tol

1151:    Level: intermediate

1153: .keywords: SNES, nonlinear, set, trust region, tolerance

1155: .seealso: SNESSetTolerances(), SNESSetMinimizationFunctionTolerance()
1156: @*/
1157: int SNESSetTrustRegionTolerance(SNES snes,PetscReal tol)
1158: {
1161:   snes->deltatol = tol;
1162:   return(0);
1163: }

1165: /* 
1166:    Duplicate the lg monitors for SNES from KSP; for some reason with 
1167:    dynamic libraries things don't work under Sun4 if we just use 
1168:    macros instead of functions
1169: */
1170: #undef __FUNCT__  
1172: int SNESLGMonitor(SNES snes,int it,PetscReal norm,void *ctx)
1173: {

1178:   KSPLGMonitor((KSP)snes,it,norm,ctx);
1179:   return(0);
1180: }

1182: #undef __FUNCT__  
1184: int SNESLGMonitorCreate(char *host,char *label,int x,int y,int m,int n,PetscDrawLG *draw)
1185: {

1189:   KSPLGMonitorCreate(host,label,x,y,m,n,draw);
1190:   return(0);
1191: }

1193: #undef __FUNCT__  
1195: int SNESLGMonitorDestroy(PetscDrawLG draw)
1196: {

1200:   KSPLGMonitorDestroy(draw);
1201:   return(0);
1202: }

1204: /* ------------ Routines to set performance monitoring options ----------- */

1206: #undef __FUNCT__  
1208: /*@C
1209:    SNESSetMonitor - Sets an ADDITIONAL function that is to be used at every
1210:    iteration of the nonlinear solver to display the iteration's 
1211:    progress.   

1213:    Collective on SNES

1215:    Input Parameters:
1216: +  snes - the SNES context
1217: .  func - monitoring routine
1218: .  mctx - [optional] user-defined context for private data for the 
1219:           monitor routine (use PETSC_NULL if no context is desitre)
1220: -  monitordestroy - [optional] routine that frees monitor context
1221:           (may be PETSC_NULL)

1223:    Calling sequence of func:
1224: $     int func(SNES snes,int its, PetscReal norm,void *mctx)

1226: +    snes - the SNES context
1227: .    its - iteration number
1228: .    norm - 2-norm function value (may be estimated)
1229: -    mctx - [optional] monitoring context

1231:    Options Database Keys:
1232: +    -snes_monitor        - sets SNESDefaultMonitor()
1233: .    -snes_xmonitor       - sets line graph monitor,
1234:                             uses SNESLGMonitorCreate()
1235: _    -snes_cancelmonitors - cancels all monitors that have
1236:                             been hardwired into a code by 
1237:                             calls to SNESSetMonitor(), but
1238:                             does not cancel those set via
1239:                             the options database.

1241:    Notes: 
1242:    Several different monitoring routines may be set by calling
1243:    SNESSetMonitor() multiple times; all will be called in the 
1244:    order in which they were set.

1246:    Level: intermediate

1248: .keywords: SNES, nonlinear, set, monitor

1250: .seealso: SNESDefaultMonitor(), SNESClearMonitor()
1251: @*/
1252: int SNESSetMonitor(SNES snes,int (*func)(SNES,int,PetscReal,void*),void *mctx,int (*monitordestroy)(void *))
1253: {
1256:   if (snes->numbermonitors >= MAXSNESMONITORS) {
1257:     SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Too many monitors set");
1258:   }

1260:   snes->monitor[snes->numbermonitors]           = func;
1261:   snes->monitordestroy[snes->numbermonitors]    = monitordestroy;
1262:   snes->monitorcontext[snes->numbermonitors++]  = (void*)mctx;
1263:   return(0);
1264: }

1266: #undef __FUNCT__  
1268: /*@C
1269:    SNESClearMonitor - Clears all the monitor functions for a SNES object.

1271:    Collective on SNES

1273:    Input Parameters:
1274: .  snes - the SNES context

1276:    Options Database Key:
1277: .  -snes_cancelmonitors - cancels all monitors that have been hardwired
1278:     into a code by calls to SNESSetMonitor(), but does not cancel those 
1279:     set via the options database

1281:    Notes: 
1282:    There is no way to clear one specific monitor from a SNES object.

1284:    Level: intermediate

1286: .keywords: SNES, nonlinear, set, monitor

1288: .seealso: SNESDefaultMonitor(), SNESSetMonitor()
1289: @*/
1290: int SNESClearMonitor(SNES snes)
1291: {
1294:   snes->numbermonitors = 0;
1295:   return(0);
1296: }

1298: #undef __FUNCT__  
1300: /*@C
1301:    SNESSetConvergenceTest - Sets the function that is to be used 
1302:    to test for convergence of the nonlinear iterative solution.   

1304:    Collective on SNES

1306:    Input Parameters:
1307: +  snes - the SNES context
1308: .  func - routine to test for convergence
1309: -  cctx - [optional] context for private data for the convergence routine 
1310:           (may be PETSC_NULL)

1312:    Calling sequence of func:
1313: $     int func (SNES snes,PetscReal xnorm,PetscReal gnorm,PetscReal f,SNESConvergedReason *reason,void *cctx)

1315: +    snes - the SNES context
1316: .    cctx - [optional] convergence context
1317: .    reason - reason for convergence/divergence
1318: .    xnorm - 2-norm of current iterate
1319: .    gnorm - 2-norm of current step
1320: -    f - 2-norm of function

1322:    Level: advanced

1324: .keywords: SNES, nonlinear, set, convergence, test

1326: .seealso: SNESConverged_LS(), SNESConverged_TR()
1327: @*/
1328: int SNESSetConvergenceTest(SNES snes,int (*func)(SNES,PetscReal,PetscReal,PetscReal,SNESConvergedReason*,void*),void *cctx)
1329: {
1332:   (snes)->converged = func;
1333:   (snes)->cnvP      = cctx;
1334:   return(0);
1335: }

1337: #undef __FUNCT__  
1339: /*@C
1340:    SNESGetConvergedReason - Gets the reason the SNES iteration was stopped.

1342:    Not Collective

1344:    Input Parameter:
1345: .  snes - the SNES context

1347:    Output Parameter:
1348: .  reason - negative value indicates diverged, positive value converged, see petscsnes.h or the 
1349:             manual pages for the individual convergence tests for complete lists

1351:    Level: intermediate

1353:    Notes: Can only be called after the call the SNESSolve() is complete.

1355: .keywords: SNES, nonlinear, set, convergence, test

1357: .seealso: SNESSetConvergenceTest(), SNESConverged_LS(), SNESConverged_TR(), SNESConvergedReason
1358: @*/
1359: int SNESGetConvergedReason(SNES snes,SNESConvergedReason *reason)
1360: {
1363:   *reason = snes->reason;
1364:   return(0);
1365: }

1367: #undef __FUNCT__  
1369: /*@
1370:    SNESSetConvergenceHistory - Sets the array used to hold the convergence history.

1372:    Collective on SNES

1374:    Input Parameters:
1375: +  snes - iterative context obtained from SNESCreate()
1376: .  a   - array to hold history
1377: .  its - integer array holds the number of linear iterations for each solve.
1378: .  na  - size of a and its
1379: -  reset - PETSC_TRUE indicates each new nonlinear solve resets the history counter to zero,
1380:            else it continues storing new values for new nonlinear solves after the old ones

1382:    Notes:
1383:    If set, this array will contain the function norms computed
1384:    at each step.

1386:    This routine is useful, e.g., when running a code for purposes
1387:    of accurate performance monitoring, when no I/O should be done
1388:    during the section of code that is being timed.

1390:    Level: intermediate

1392: .keywords: SNES, set, convergence, history

1394: .seealso: SNESGetConvergenceHistory()

1396: @*/
1397: int SNESSetConvergenceHistory(SNES snes,PetscReal *a,int *its,int na,PetscTruth reset)
1398: {
1402:   snes->conv_hist       = a;
1403:   snes->conv_hist_its   = its;
1404:   snes->conv_hist_max   = na;
1405:   snes->conv_hist_reset = reset;
1406:   return(0);
1407: }

1409: #undef __FUNCT__  
1411: /*@C
1412:    SNESGetConvergenceHistory - Gets the array used to hold the convergence history.

1414:    Collective on SNES

1416:    Input Parameter:
1417: .  snes - iterative context obtained from SNESCreate()

1419:    Output Parameters:
1420: .  a   - array to hold history
1421: .  its - integer array holds the number of linear iterations (or
1422:          negative if not converged) for each solve.
1423: -  na  - size of a and its

1425:    Notes:
1426:     The calling sequence for this routine in Fortran is
1427: $   call SNESGetConvergenceHistory(SNES snes, integer na, integer ierr)

1429:    This routine is useful, e.g., when running a code for purposes
1430:    of accurate performance monitoring, when no I/O should be done
1431:    during the section of code that is being timed.

1433:    Level: intermediate

1435: .keywords: SNES, get, convergence, history

1437: .seealso: SNESSetConvergencHistory()

1439: @*/
1440: int SNESGetConvergenceHistory(SNES snes,PetscReal **a,int **its,int *na)
1441: {
1444:   if (a)   *a   = snes->conv_hist;
1445:   if (its) *its = snes->conv_hist_its;
1446:   if (na) *na   = snes->conv_hist_len;
1447:   return(0);
1448: }

1450: #undef __FUNCT__  
1452: /*@
1453:   SNESSetRhsBC - Sets the function which applies boundary conditions
1454:   to the Rhs of each system.

1456:   Collective on SNES

1458:   Input Parameters:
1459: . snes - The nonlinear solver context
1460: . func - The function

1462:   Calling sequence of func:
1463: . func (SNES snes, Vec rhs, void *ctx);

1465: . rhs - The current rhs vector
1466: . ctx - The user-context

1468:   Level: intermediate

1470: .keywords: SNES, Rhs, boundary conditions
1471: .seealso SNESDefaultRhsBC(), SNESSetSolutionBC(), SNESSetUpdate()
1472: @*/
1473: int SNESSetRhsBC(SNES snes, int (*func)(SNES, Vec, void *))
1474: {
1477:   snes->applyrhsbc = func;
1478:   return(0);
1479: }

1481: #undef __FUNCT__  
1483: /*@
1484:   SNESDefaultRhsBC - The default boundary condition function which does nothing.

1486:   Not collective

1488:   Input Parameters:
1489: . snes - The nonlinear solver context
1490: . rhs  - The Rhs
1491: . ctx  - The user-context

1493:   Level: beginner

1495: .keywords: SNES, Rhs, boundary conditions
1496: .seealso SNESSetRhsBC(), SNESDefaultSolutionBC(), SNESDefaultUpdate()
1497: @*/
1498: int SNESDefaultRhsBC(SNES snes, Vec rhs, void *ctx)
1499: {
1501:   return(0);
1502: }

1504: #undef __FUNCT__  
1506: /*@
1507:   SNESSetSolutionBC - Sets the function which applies boundary conditions
1508:   to the solution of each system.

1510:   Collective on SNES

1512:   Input Parameters:
1513: . snes - The nonlinear solver context
1514: . func - The function

1516:   Calling sequence of func:
1517: . func (SNES snes, Vec rsol, void *ctx);

1519: . sol - The current solution vector
1520: . ctx - The user-context

1522:   Level: intermediate

1524: .keywords: SNES, solution, boundary conditions
1525: .seealso SNESDefaultSolutionBC(), SNESSetRhsBC(), SNESSetUpdate()
1526: @*/
1527: int SNESSetSolutionBC(SNES snes, int (*func)(SNES, Vec, void *))
1528: {
1531:   snes->applysolbc = func;
1532:   return(0);
1533: }

1535: #undef __FUNCT__  
1537: /*@
1538:   SNESDefaultSolutionBC - The default boundary condition function which does nothing.

1540:   Not collective

1542:   Input Parameters:
1543: . snes - The nonlinear solver context
1544: . sol  - The solution
1545: . ctx  - The user-context

1547:   Level: beginner

1549: .keywords: SNES, solution, boundary conditions
1550: .seealso SNESSetSolutionBC(), SNESDefaultRhsBC(), SNESDefaultUpdate()
1551: @*/
1552: int SNESDefaultSolutionBC(SNES snes, Vec sol, void *ctx)
1553: {
1555:   return(0);
1556: }

1558: #undef __FUNCT__  
1560: /*@
1561:   SNESSetUpdate - Sets the general-purpose update function called
1562:   at the beginning of every step of the iteration.

1564:   Collective on SNES

1566:   Input Parameters:
1567: . snes - The nonlinear solver context
1568: . func - The function

1570:   Calling sequence of func:
1571: . func (TS ts, int step);

1573: . step - The current step of the iteration

1575:   Level: intermediate

1577: .keywords: SNES, update
1578: .seealso SNESDefaultUpdate(), SNESSetRhsBC(), SNESSetSolutionBC()
1579: @*/
1580: int SNESSetUpdate(SNES snes, int (*func)(SNES, int))
1581: {
1584:   snes->update = func;
1585:   return(0);
1586: }

1588: #undef __FUNCT__  
1590: /*@
1591:   SNESDefaultUpdate - The default update function which does nothing.

1593:   Not collective

1595:   Input Parameters:
1596: . snes - The nonlinear solver context
1597: . step - The current step of the iteration

1599:   Level: intermediate

1601: .keywords: SNES, update
1602: .seealso SNESSetUpdate(), SNESDefaultRhsBC(), SNESDefaultSolutionBC()
1603: @*/
1604: int SNESDefaultUpdate(SNES snes, int step)
1605: {
1607:   return(0);
1608: }

1610: #undef __FUNCT__  
1612: /*
1613:    SNESScaleStep_Private - Scales a step so that its length is less than the
1614:    positive parameter delta.

1616:     Input Parameters:
1617: +   snes - the SNES context
1618: .   y - approximate solution of linear system
1619: .   fnorm - 2-norm of current function
1620: -   delta - trust region size

1622:     Output Parameters:
1623: +   gpnorm - predicted function norm at the new point, assuming local 
1624:     linearization.  The value is zero if the step lies within the trust 
1625:     region, and exceeds zero otherwise.
1626: -   ynorm - 2-norm of the step

1628:     Note:
1629:     For non-trust region methods such as SNESLS, the parameter delta 
1630:     is set to be the maximum allowable step size.  

1632: .keywords: SNES, nonlinear, scale, step
1633: */
1634: int SNESScaleStep_Private(SNES snes,Vec y,PetscReal *fnorm,PetscReal *delta,PetscReal *gpnorm,PetscReal *ynorm)
1635: {
1636:   PetscReal   nrm;
1637:   PetscScalar cnorm;
1638:   int         ierr;


1645:   VecNorm(y,NORM_2,&nrm);
1646:   if (nrm > *delta) {
1647:      nrm = *delta/nrm;
1648:      *gpnorm = (1.0 - nrm)*(*fnorm);
1649:      cnorm = nrm;
1650:      VecScale(&cnorm,y);
1651:      *ynorm = *delta;
1652:   } else {
1653:      *gpnorm = 0.0;
1654:      *ynorm = nrm;
1655:   }
1656:   return(0);
1657: }

1659: #undef __FUNCT__  
1661: /*@
1662:    SNESSolve - Solves a nonlinear system.  Call SNESSolve after calling 
1663:    SNESCreate() and optional routines of the form SNESSetXXX().

1665:    Collective on SNES

1667:    Input Parameters:
1668: +  snes - the SNES context
1669: -  x - the solution vector

1671:    Output Parameter:
1672: .  its - number of iterations until termination

1674:    Notes:
1675:    The user should initialize the vector,x, with the initial guess
1676:    for the nonlinear solve prior to calling SNESSolve.  In particular,
1677:    to employ an initial guess of zero, the user should explicitly set
1678:    this vector to zero by calling VecSet().

1680:    Level: beginner

1682: .keywords: SNES, nonlinear, solve

1684: .seealso: SNESCreate(), SNESDestroy()
1685: @*/
1686: int SNESSolve(SNES snes,Vec x,int *its)
1687: {
1688:   int        ierr;
1689:   PetscTruth flg;

1696:   if (!snes->solve) SETERRQ(1,"SNESSetType() or SNESSetFromOptions() must be called before SNESSolve()");

1698:   if (!snes->setupcalled) {SNESSetUp(snes,x);}
1699:   else {snes->vec_sol = snes->vec_sol_always = x;}
1700:   if (snes->conv_hist_reset == PETSC_TRUE) snes->conv_hist_len = 0;
1701:   PetscLogEventBegin(SNES_Solve,snes,0,0,0);
1702:   snes->nfuncs = 0; snes->linear_its = 0; snes->numFailures = 0;
1703:   (*(snes)->solve)(snes,its);
1704:   PetscLogEventEnd(SNES_Solve,snes,0,0,0);
1705:   PetscOptionsHasName(snes->prefix,"-snes_view",&flg);
1706:   if (flg && !PetscPreLoadingOn) { SNESView(snes,PETSC_VIEWER_STDOUT_(snes->comm)); }
1707:   PetscOptionsHasName(snes->prefix,"-snes_test_local_min",&flg);
1708:   if (flg && !PetscPreLoadingOn) { SNESTestLocalMin(snes); }
1709:   return(0);
1710: }

1712: /* --------- Internal routines for SNES Package --------- */

1714: #undef __FUNCT__  
1716: /*@C
1717:    SNESSetType - Sets the method for the nonlinear solver.  

1719:    Collective on SNES

1721:    Input Parameters:
1722: +  snes - the SNES context
1723: -  type - a known method

1725:    Options Database Key:
1726: .  -snes_type <type> - Sets the method; use -help for a list
1727:    of available methods (for instance, ls or tr)

1729:    Notes:
1730:    See "petsc/include/petscsnes.h" for available methods (for instance)
1731: +    SNESLS - Newton's method with line search
1732:      (systems of nonlinear equations)
1733: .    SNESTR - Newton's method with trust region
1734:      (systems of nonlinear equations)

1736:   Normally, it is best to use the SNESSetFromOptions() command and then
1737:   set the SNES solver type from the options database rather than by using
1738:   this routine.  Using the options database provides the user with
1739:   maximum flexibility in evaluating the many nonlinear solvers.
1740:   The SNESSetType() routine is provided for those situations where it
1741:   is necessary to set the nonlinear solver independently of the command
1742:   line or options database.  This might be the case, for example, when
1743:   the choice of solver changes during the execution of the program,
1744:   and the user's application is taking responsibility for choosing the
1745:   appropriate method.

1747:   Level: intermediate

1749: .keywords: SNES, set, type

1751: .seealso: SNESType, SNESCreate()

1753: @*/
1754: int SNESSetType(SNES snes,SNESType type)
1755: {
1756:   int        ierr,(*r)(SNES);
1757:   PetscTruth match;


1763:   PetscTypeCompare((PetscObject)snes,type,&match);
1764:   if (match) return(0);

1766:   if (snes->setupcalled) {
1767:     ierr       = (*(snes)->destroy)(snes);
1768:     snes->data = 0;
1769:   }

1771:   /* Get the function pointers for the iterative method requested */
1772:   if (!SNESRegisterAllCalled) {SNESRegisterAll(PETSC_NULL);}

1774:    PetscFListFind(snes->comm,SNESList,type,(void (**)(void)) &r);

1776:   if (!r) SETERRQ1(1,"Unable to find requested SNES type %s",type);

1778:   if (snes->data) {PetscFree(snes->data);}
1779:   snes->data = 0;
1780:   (*r)(snes);
1781:   PetscObjectChangeTypeName((PetscObject)snes,type);

1783:   return(0);
1784: }


1787: /* --------------------------------------------------------------------- */
1788: #undef __FUNCT__  
1790: /*@C
1791:    SNESRegisterDestroy - Frees the list of nonlinear solvers that were
1792:    registered by SNESRegisterDynamic().

1794:    Not Collective

1796:    Level: advanced

1798: .keywords: SNES, nonlinear, register, destroy

1800: .seealso: SNESRegisterAll(), SNESRegisterAll()
1801: @*/
1802: int SNESRegisterDestroy(void)
1803: {

1807:   if (SNESList) {
1808:     PetscFListDestroy(&SNESList);
1809:     SNESList = 0;
1810:   }
1811:   SNESRegisterAllCalled = PETSC_FALSE;
1812:   return(0);
1813: }

1815: #undef __FUNCT__  
1817: /*@C
1818:    SNESGetType - Gets the SNES method type and name (as a string).

1820:    Not Collective

1822:    Input Parameter:
1823: .  snes - nonlinear solver context

1825:    Output Parameter:
1826: .  type - SNES method (a character string)

1828:    Level: intermediate

1830: .keywords: SNES, nonlinear, get, type, name
1831: @*/
1832: int SNESGetType(SNES snes,SNESType *type)
1833: {
1836:   *type = snes->type_name;
1837:   return(0);
1838: }

1840: #undef __FUNCT__  
1842: /*@C
1843:    SNESGetSolution - Returns the vector where the approximate solution is
1844:    stored.

1846:    Not Collective, but Vec is parallel if SNES is parallel

1848:    Input Parameter:
1849: .  snes - the SNES context

1851:    Output Parameter:
1852: .  x - the solution

1854:    Level: advanced

1856: .keywords: SNES, nonlinear, get, solution

1858: .seealso: SNESGetFunction(), SNESGetSolutionUpdate()
1859: @*/
1860: int SNESGetSolution(SNES snes,Vec *x)
1861: {
1864:   *x = snes->vec_sol_always;
1865:   return(0);
1866: }

1868: #undef __FUNCT__  
1870: /*@C
1871:    SNESGetSolutionUpdate - Returns the vector where the solution update is
1872:    stored. 

1874:    Not Collective, but Vec is parallel if SNES is parallel

1876:    Input Parameter:
1877: .  snes - the SNES context

1879:    Output Parameter:
1880: .  x - the solution update

1882:    Level: advanced

1884: .keywords: SNES, nonlinear, get, solution, update

1886: .seealso: SNESGetSolution(), SNESGetFunction
1887: @*/
1888: int SNESGetSolutionUpdate(SNES snes,Vec *x)
1889: {
1892:   *x = snes->vec_sol_update_always;
1893:   return(0);
1894: }

1896: #undef __FUNCT__  
1898: /*@C
1899:    SNESGetFunction - Returns the vector where the function is stored.

1901:    Not Collective, but Vec is parallel if SNES is parallel

1903:    Input Parameter:
1904: .  snes - the SNES context

1906:    Output Parameter:
1907: +  r - the function (or PETSC_NULL)
1908: .  ctx - the function context (or PETSC_NULL)
1909: -  func - the function (or PETSC_NULL)

1911:    Level: advanced

1913: .keywords: SNES, nonlinear, get, function

1915: .seealso: SNESSetFunction(), SNESGetSolution()
1916: @*/
1917: int SNESGetFunction(SNES snes,Vec *r,void **ctx,int (**func)(SNES,Vec,Vec,void*))
1918: {
1921:   if (r)    *r    = snes->vec_func_always;
1922:   if (ctx)  *ctx  = snes->funP;
1923:   if (func) *func = snes->computefunction;
1924:   return(0);
1925: }

1927: #undef __FUNCT__  
1929: /*@C
1930:    SNESSetOptionsPrefix - Sets the prefix used for searching for all 
1931:    SNES options in the database.

1933:    Collective on SNES

1935:    Input Parameter:
1936: +  snes - the SNES context
1937: -  prefix - the prefix to prepend to all option names

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

1943:    Level: advanced

1945: .keywords: SNES, set, options, prefix, database

1947: .seealso: SNESSetFromOptions()
1948: @*/
1949: int SNESSetOptionsPrefix(SNES snes,char *prefix)
1950: {

1955:   PetscObjectSetOptionsPrefix((PetscObject)snes,prefix);
1956:   SLESSetOptionsPrefix(snes->sles,prefix);
1957:   return(0);
1958: }

1960: #undef __FUNCT__  
1962: /*@C
1963:    SNESAppendOptionsPrefix - Appends to the prefix used for searching for all 
1964:    SNES options in the database.

1966:    Collective on SNES

1968:    Input Parameters:
1969: +  snes - the SNES context
1970: -  prefix - the prefix to prepend to all option names

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

1976:    Level: advanced

1978: .keywords: SNES, append, options, prefix, database

1980: .seealso: SNESGetOptionsPrefix()
1981: @*/
1982: int SNESAppendOptionsPrefix(SNES snes,char *prefix)
1983: {
1985: 
1988:   PetscObjectAppendOptionsPrefix((PetscObject)snes,prefix);
1989:   SLESAppendOptionsPrefix(snes->sles,prefix);
1990:   return(0);
1991: }

1993: #undef __FUNCT__  
1995: /*@C
1996:    SNESGetOptionsPrefix - Sets the prefix used for searching for all 
1997:    SNES options in the database.

1999:    Not Collective

2001:    Input Parameter:
2002: .  snes - the SNES context

2004:    Output Parameter:
2005: .  prefix - pointer to the prefix string used

2007:    Notes: On the fortran side, the user should pass in a string 'prifix' of
2008:    sufficient length to hold the prefix.

2010:    Level: advanced

2012: .keywords: SNES, get, options, prefix, database

2014: .seealso: SNESAppendOptionsPrefix()
2015: @*/
2016: int SNESGetOptionsPrefix(SNES snes,char **prefix)
2017: {

2022:   PetscObjectGetOptionsPrefix((PetscObject)snes,prefix);
2023:   return(0);
2024: }

2026: /*MC
2027:    SNESRegisterDynamic - Adds a method to the nonlinear solver package.

2029:    Synopsis:
2030:    int SNESRegisterDynamic(char *name_solver,char *path,char *name_create,int (*routine_create)(SNES))

2032:    Not collective

2034:    Input Parameters:
2035: +  name_solver - name of a new user-defined solver
2036: .  path - path (either absolute or relative) the library containing this solver
2037: .  name_create - name of routine to create method context
2038: -  routine_create - routine to create method context

2040:    Notes:
2041:    SNESRegisterDynamic() may be called multiple times to add several user-defined solvers.

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

2046:    Environmental variables such as ${PETSC_ARCH}, ${PETSC_DIR}, ${PETSC_LIB_DIR}, ${BOPT},
2047:    and others of the form ${any_environmental_variable} occuring in pathname will be 
2048:    replaced with appropriate values.

2050:    Sample usage:
2051: .vb
2052:    SNESRegisterDynamic("my_solver",/home/username/my_lib/lib/libg/solaris/mylib.a,
2053:                 "MySolverCreate",MySolverCreate);
2054: .ve

2056:    Then, your solver can be chosen with the procedural interface via
2057: $     SNESSetType(snes,"my_solver")
2058:    or at runtime via the option
2059: $     -snes_type my_solver

2061:    Level: advanced

2063:     Note: If your function is not being put into a shared library then use SNESRegister() instead

2065: .keywords: SNES, nonlinear, register

2067: .seealso: SNESRegisterAll(), SNESRegisterDestroy()
2068: M*/

2070: #undef __FUNCT__  
2072: /*@C
2073:   SNESRegister - See SNESRegisterDynamic()

2075:   Level: advanced
2076: @*/
2077: int SNESRegister(char *sname,char *path,char *name,int (*function)(SNES))
2078: {
2079:   char fullname[256];
2080:   int  ierr;

2083:   PetscFListConcat(path,name,fullname);
2084:   PetscFListAdd(&SNESList,sname,fullname,(void (*)(void))function);
2085:   return(0);
2086: }

2088: #undef __FUNCT__  
2090: int SNESTestLocalMin(SNES snes)
2091: {
2092:   int         ierr,N,i,j;
2093:   Vec         u,uh,fh;
2094:   PetscScalar value;
2095:   PetscReal   norm;

2098:   SNESGetSolution(snes,&u);
2099:   VecDuplicate(u,&uh);
2100:   VecDuplicate(u,&fh);

2102:   /* currently only works for sequential */
2103:   PetscPrintf(PETSC_COMM_WORLD,"Testing FormFunction() for local minn");
2104:   VecGetSize(u,&N);
2105:   for (i=0; i<N; i++) {
2106:     VecCopy(u,uh);
2107:     ierr  = PetscPrintf(PETSC_COMM_WORLD,"i = %dn",i);
2108:     for (j=-10; j<11; j++) {
2109:       value = PetscSign(j)*exp(PetscAbs(j)-10.0);
2110:       ierr  = VecSetValue(uh,i,value,ADD_VALUES);
2111:       ierr  = (*snes->computefunction)(snes,uh,fh,snes->funP);
2112:       ierr  = VecNorm(fh,NORM_2,&norm);
2113:       ierr  = PetscPrintf(PETSC_COMM_WORLD,"       j norm %d %18.16en",j,norm);
2114:       value = -value;
2115:       ierr  = VecSetValue(uh,i,value,ADD_VALUES);
2116:     }
2117:   }
2118:   VecDestroy(uh);
2119:   VecDestroy(fh);
2120:   return(0);
2121: }