#define PETSCSNES_DLL #include "src/snes/impls/ls/ls.h" /* Checks if J^T F = 0 which implies we've found a local minimum of the function, but not a zero. In the case when one cannot compute J^T F we use the fact that 0 = (J^T F)^T W = F^T J W iff W not in the null space of J. Thanks for Jorge More for this trick. */ #undef __FUNCT__ #define __FUNCT__ "SNESLSCheckLocalMin_Private" PetscErrorCode SNESLSCheckLocalMin_Private(Mat A,Vec F,Vec W,PetscReal fnorm,PetscTruth *ismin) { PetscReal a1; PetscErrorCode ierr; PetscTruth hastranspose; PetscFunctionBegin; *ismin = PETSC_FALSE; ierr = MatHasOperation(A,MATOP_MULT_TRANSPOSE,&hastranspose);CHKERRQ(ierr); if (hastranspose) { /* Compute || J^T F|| */ ierr = MatMultTranspose(A,F,W);CHKERRQ(ierr); ierr = VecNorm(W,NORM_2,&a1);CHKERRQ(ierr); ierr = PetscInfo1(0,"|| J^T F|| %G near zero implies found a local minimum\n",a1/fnorm);CHKERRQ(ierr); if (a1/fnorm < 1.e-4) *ismin = PETSC_TRUE; } else { Vec work; PetscScalar result; PetscReal wnorm; ierr = VecSetRandom(W,PETSC_NULL);CHKERRQ(ierr); ierr = VecNorm(W,NORM_2,&wnorm);CHKERRQ(ierr); ierr = VecDuplicate(W,&work);CHKERRQ(ierr); ierr = MatMult(A,W,work);CHKERRQ(ierr); ierr = VecDot(F,work,&result);CHKERRQ(ierr); ierr = VecDestroy(work);CHKERRQ(ierr); a1 = PetscAbsScalar(result)/(fnorm*wnorm); ierr = PetscInfo1(0,"(F^T J random)/(|| F ||*||random|| %G near zero implies found a local minimum\n",a1);CHKERRQ(ierr); if (a1 < 1.e-4) *ismin = PETSC_TRUE; } PetscFunctionReturn(0); } /* Checks if J^T(F - J*X) = 0 */ #undef __FUNCT__ #define __FUNCT__ "SNESLSCheckResidual_Private" PetscErrorCode SNESLSCheckResidual_Private(Mat A,Vec F,Vec X,Vec W1,Vec W2) { PetscReal a1,a2; PetscErrorCode ierr; PetscTruth hastranspose; PetscFunctionBegin; ierr = MatHasOperation(A,MATOP_MULT_TRANSPOSE,&hastranspose);CHKERRQ(ierr); if (hastranspose) { ierr = MatMult(A,X,W1);CHKERRQ(ierr); ierr = VecAXPY(W1,-1.0,F);CHKERRQ(ierr); /* Compute || J^T W|| */ ierr = MatMultTranspose(A,W1,W2);CHKERRQ(ierr); ierr = VecNorm(W1,NORM_2,&a1);CHKERRQ(ierr); ierr = VecNorm(W2,NORM_2,&a2);CHKERRQ(ierr); if (a1) { ierr = PetscInfo1(0,"||J^T(F-Ax)||/||F-AX|| %G near zero implies inconsistent rhs\n",a2/a1);CHKERRQ(ierr); } } PetscFunctionReturn(0); } /* -------------------------------------------------------------------- This file implements a truncated Newton method with a line search, for solving a system of nonlinear equations, using the KSP, Vec, and Mat interfaces for linear solvers, vectors, and matrices, respectively. The following basic routines are required for each nonlinear solver: SNESCreate_XXX() - Creates a nonlinear solver context SNESSetFromOptions_XXX() - Sets runtime options SNESSolve_XXX() - Solves the nonlinear system SNESDestroy_XXX() - Destroys the nonlinear solver context The suffix "_XXX" denotes a particular implementation, in this case we use _LS (e.g., SNESCreate_LS, SNESSolve_LS) for solving systems of nonlinear equations with a line search (LS) method. These routines are actually called via the common user interface routines SNESCreate(), SNESSetFromOptions(), SNESSolve(), and SNESDestroy(), so the application code interface remains identical for all nonlinear solvers. Another key routine is: SNESSetUp_XXX() - Prepares for the use of a nonlinear solver by setting data structures and options. The interface routine SNESSetUp() is not usually called directly by the user, but instead is called by SNESSolve() if necessary. Additional basic routines are: SNESView_XXX() - Prints details of runtime options that have actually been used. These are called by application codes via the interface routines SNESView(). The various types of solvers (preconditioners, Krylov subspace methods, nonlinear solvers, timesteppers) are all organized similarly, so the above description applies to these categories also. -------------------------------------------------------------------- */ /* SNESSolve_LS - Solves a nonlinear system with a truncated Newton method with a line search. Input Parameters: . snes - the SNES context Output Parameter: . outits - number of iterations until termination Application Interface Routine: SNESSolve() Notes: This implements essentially a truncated Newton method with a line search. By default a cubic backtracking line search is employed, as described in the text "Numerical Methods for Unconstrained Optimization and Nonlinear Equations" by Dennis and Schnabel. */ #undef __FUNCT__ #define __FUNCT__ "SNESSolve_LS" PetscErrorCode SNESSolve_LS(SNES snes) { SNES_LS *neP = (SNES_LS*)snes->data; PetscErrorCode ierr; PetscInt maxits,i,lits; PetscTruth lssucceed; MatStructure flg = DIFFERENT_NONZERO_PATTERN; PetscReal fnorm,gnorm,xnorm,ynorm; Vec Y,X,F,G,W,TMP; KSPConvergedReason kspreason; PetscFunctionBegin; snes->numFailures = 0; snes->numLinearSolveFailures = 0; snes->reason = SNES_CONVERGED_ITERATING; maxits = snes->max_its; /* maximum number of iterations */ X = snes->vec_sol; /* solution vector */ F = snes->vec_func; /* residual vector */ Y = snes->work[0]; /* work vectors */ G = snes->work[1]; W = snes->work[2]; ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); snes->iter = 0; ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); ierr = VecNorm(F,NORM_2,&fnorm);CHKERRQ(ierr); /* fnorm <- ||F|| */ if (fnorm != fnorm) SETERRQ(PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); snes->norm = fnorm; ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); SNESLogConvHistory(snes,fnorm,0); SNESMonitor(snes,0,fnorm); /* set parameter for default relative tolerance convergence test */ snes->ttol = fnorm*snes->rtol; if (snes->ops->converged) { ierr = (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); } if (snes->reason) PetscFunctionReturn(0); for (i=0; iops->update) { ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); } /* Solve J Y = F, where J is Jacobian matrix */ ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); ierr = KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre,flg);CHKERRQ(ierr); ierr = SNES_KSPSolve(snes,snes->ksp,F,Y);CHKERRQ(ierr); ierr = KSPGetConvergedReason(snes->ksp,&kspreason);CHKERRQ(ierr); if (kspreason < 0) { if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) { snes->reason = SNES_DIVERGED_LINEAR_SOLVE; PetscFunctionReturn(0); } } ierr = KSPGetIterationNumber(snes->ksp,&lits);CHKERRQ(ierr); snes->linear_its += lits; ierr = PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);CHKERRQ(ierr); if (neP->precheckstep) { PetscTruth changed_y = PETSC_FALSE; ierr = (*neP->precheckstep)(snes,X,Y,neP->precheck,&changed_y);CHKERRQ(ierr); } if (PetscLogPrintInfo){ ierr = SNESLSCheckResidual_Private(snes->jacobian,F,Y,G,W);CHKERRQ(ierr); } /* Compute a (scaled) negative update in the line search routine: Y <- X - lambda*Y and evaluate G = function(Y) (depends on the line search). */ ierr = VecCopy(Y,snes->vec_sol_update_always);CHKERRQ(ierr); ynorm = 1; gnorm = fnorm; ierr = (*neP->LineSearch)(snes,neP->lsP,X,F,G,Y,W,fnorm,&ynorm,&gnorm,&lssucceed);CHKERRQ(ierr); ierr = PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",fnorm,gnorm,ynorm,(int)lssucceed);CHKERRQ(ierr); TMP = F; F = G; snes->vec_func_always = F; G = TMP; TMP = X; X = W; snes->vec_sol_always = X; W = TMP; fnorm = gnorm; if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break; ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); snes->iter = i+1; snes->norm = fnorm; ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); SNESLogConvHistory(snes,fnorm,lits); SNESMonitor(snes,i+1,fnorm); if (!lssucceed) { PetscTruth ismin; if (++snes->numFailures >= snes->maxFailures) { snes->reason = SNES_DIVERGED_LS_FAILURE; ierr = SNESLSCheckLocalMin_Private(snes->jacobian,F,W,fnorm,&ismin);CHKERRQ(ierr); if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN; break; } } /* Test for convergence */ if (snes->ops->converged) { ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); /* xnorm = || X || */ ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); if (snes->reason) break; } } if (X != snes->vec_sol) { ierr = VecCopy(X,snes->vec_sol);CHKERRQ(ierr); } if (F != snes->vec_func) { ierr = VecCopy(F,snes->vec_func);CHKERRQ(ierr); } snes->vec_sol_always = snes->vec_sol; snes->vec_func_always = snes->vec_func; if (i == maxits) { ierr = PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);CHKERRQ(ierr); snes->reason = SNES_DIVERGED_MAX_IT; } PetscFunctionReturn(0); } /* -------------------------------------------------------------------------- */ /* SNESSetUp_LS - Sets up the internal data structures for the later use of the SNESLS nonlinear solver. Input Parameter: . snes - the SNES context . x - the solution vector Application Interface Routine: SNESSetUp() Notes: For basic use of the SNES solvers, the user need not explicitly call SNESSetUp(), since these actions will automatically occur during the call to SNESSolve(). */ #undef __FUNCT__ #define __FUNCT__ "SNESSetUp_LS" PetscErrorCode SNESSetUp_LS(SNES snes) { PetscErrorCode ierr; PetscFunctionBegin; if (!snes->work) { snes->nwork = 4; ierr = VecDuplicateVecs(snes->vec_sol,snes->nwork,&snes->work);CHKERRQ(ierr); ierr = PetscLogObjectParents(snes,snes->nwork,snes->work);CHKERRQ(ierr); snes->vec_sol_update_always = snes->work[3]; } PetscFunctionReturn(0); } /* -------------------------------------------------------------------------- */ /* SNESDestroy_LS - Destroys the private SNES_LS context that was created with SNESCreate_LS(). Input Parameter: . snes - the SNES context Application Interface Routine: SNESDestroy() */ #undef __FUNCT__ #define __FUNCT__ "SNESDestroy_LS" PetscErrorCode SNESDestroy_LS(SNES snes) { PetscErrorCode ierr; PetscFunctionBegin; if (snes->nwork) { ierr = VecDestroyVecs(snes->work,snes->nwork);CHKERRQ(ierr); } ierr = PetscFree(snes->data);CHKERRQ(ierr); /* clear composed functions */ ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSet_C","",PETSC_NULL);CHKERRQ(ierr); ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetPostCheck_C","",PETSC_NULL);CHKERRQ(ierr); ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetPreCheck_C","",PETSC_NULL);CHKERRQ(ierr); PetscFunctionReturn(0); } /* -------------------------------------------------------------------------- */ #undef __FUNCT__ #define __FUNCT__ "SNESLineSearchNo" /*@C SNESLineSearchNo - This routine is not a line search at all; it simply uses the full Newton step. Thus, this routine is intended to serve as a template and is not recommended for general use. Collective on SNES and Vec Input Parameters: + snes - nonlinear context . lsctx - optional context for line search (not used here) . x - current iterate . f - residual evaluated at x . y - search direction . w - work vector - fnorm - 2-norm of f Output Parameters: + g - residual evaluated at new iterate y . w - new iterate . gnorm - 2-norm of g . ynorm - 2-norm of search length - flag - PETSC_TRUE on success, PETSC_FALSE on failure Options Database Key: . -snes_ls basic - Activates SNESLineSearchNo() Level: advanced .keywords: SNES, nonlinear, line search, cubic .seealso: SNESLineSearchCubic(), SNESLineSearchQuadratic(), SNESLineSearchSet(), SNESLineSearchNoNorms() @*/ PetscErrorCode PETSCSNES_DLLEXPORT SNESLineSearchNo(SNES snes,void *lsctx,Vec x,Vec f,Vec g,Vec y,Vec w,PetscReal fnorm,PetscReal *ynorm,PetscReal *gnorm,PetscTruth *flag) { PetscErrorCode ierr; SNES_LS *neP = (SNES_LS*)snes->data; PetscTruth changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; PetscFunctionBegin; *flag = PETSC_TRUE; ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr); /* ynorm = || y || */ ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ if (neP->postcheckstep) { ierr = (*neP->postcheckstep)(snes,x,y,w,neP->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); } if (changed_y) { ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ } ierr = SNESComputeFunction(snes,w,g); if (PetscExceptionValue(ierr)) { PetscErrorCode pierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(pierr); } CHKERRQ(ierr); ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); /* gnorm = || g || */ if (*gnorm != *gnorm) SETERRQ(PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); PetscFunctionReturn(0); } /* -------------------------------------------------------------------------- */ #undef __FUNCT__ #define __FUNCT__ "SNESLineSearchNoNorms" /*@C SNESLineSearchNoNorms - This routine is not a line search at all; it simply uses the full Newton step. This version does not even compute the norm of the function or search direction; this is intended only when you know the full step is fine and are not checking for convergence of the nonlinear iteration (for example, you are running always for a fixed number of Newton steps). Collective on SNES and Vec Input Parameters: + snes - nonlinear context . lsctx - optional context for line search (not used here) . x - current iterate . f - residual evaluated at x . y - search direction . w - work vector - fnorm - 2-norm of f Output Parameters: + g - residual evaluated at new iterate y . w - new iterate . gnorm - not changed . ynorm - not changed - flag - set to PETSC_TRUE indicating a successful line search Options Database Key: . -snes_ls basicnonorms - Activates SNESLineSearchNoNorms() Notes: SNESLineSearchNoNorms() must be used in conjunction with either the options $ -snes_no_convergence_test -snes_max_it or alternatively a user-defined custom test set via SNESSetConvergenceTest(); or a -snes_max_it of 1, otherwise, the SNES solver will generate an error. During the final iteration this will not evaluate the function at the solution point. This is to save a function evaluation while using pseudo-timestepping. The residual norms printed by monitoring routines such as SNESMonitorDefault() (as activated via -snes_monitor) will not be correct, since they are not computed. Level: advanced .keywords: SNES, nonlinear, line search, cubic .seealso: SNESLineSearchCubic(), SNESLineSearchQuadratic(), SNESLineSearchSet(), SNESLineSearchNo() @*/ PetscErrorCode PETSCSNES_DLLEXPORT SNESLineSearchNoNorms(SNES snes,void *lsctx,Vec x,Vec f,Vec g,Vec y,Vec w,PetscReal fnorm,PetscReal *ynorm,PetscReal *gnorm,PetscTruth *flag) { PetscErrorCode ierr; SNES_LS *neP = (SNES_LS*)snes->data; PetscTruth changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; PetscFunctionBegin; *flag = PETSC_TRUE; ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ if (neP->postcheckstep) { ierr = (*neP->postcheckstep)(snes,x,y,w,neP->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); } if (changed_y) { ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ } /* don't evaluate function the last time through */ if (snes->iter < snes->max_its-1) { ierr = SNESComputeFunction(snes,w,g); if (PetscExceptionValue(ierr)) { PetscErrorCode pierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(pierr); } CHKERRQ(ierr); } ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); PetscFunctionReturn(0); } /* -------------------------------------------------------------------------- */ #undef __FUNCT__ #define __FUNCT__ "SNESLineSearchCubic" /*@C SNESLineSearchCubic - Performs a cubic line search (default line search method). Collective on SNES Input Parameters: + snes - nonlinear context . lsctx - optional context for line search (not used here) . x - current iterate . f - residual evaluated at x . y - search direction . w - work vector - fnorm - 2-norm of f Output Parameters: + g - residual evaluated at new iterate y . w - new iterate . gnorm - 2-norm of g . ynorm - 2-norm of search length - flag - PETSC_TRUE if line search succeeds; PETSC_FALSE on failure. Options Database Key: $ -snes_ls cubic - Activates SNESLineSearchCubic() Notes: This line search is taken from "Numerical Methods for Unconstrained Optimization and Nonlinear Equations" by Dennis and Schnabel, page 325. Level: advanced .keywords: SNES, nonlinear, line search, cubic .seealso: SNESLineSearchQuadratic(), SNESLineSearchNo(), SNESLineSearchSet(), SNESLineSearchNoNorms() @*/ PetscErrorCode PETSCSNES_DLLEXPORT SNESLineSearchCubic(SNES snes,void *lsctx,Vec x,Vec f,Vec g,Vec y,Vec w,PetscReal fnorm,PetscReal *ynorm,PetscReal *gnorm,PetscTruth *flag) { /* Note that for line search purposes we work with with the related minimization problem: min z(x): R^n -> R, where z(x) = .5 * fnorm*fnorm, and fnorm = || f ||_2. */ PetscReal steptol,initslope,lambdaprev,gnormprev,a,b,d,t1,t2,rellength; PetscReal maxstep,minlambda,alpha,lambda,lambdatemp; #if defined(PETSC_USE_COMPLEX) PetscScalar cinitslope,clambda; #endif PetscErrorCode ierr; PetscInt count; SNES_LS *neP = (SNES_LS*)snes->data; PetscScalar scale; PetscTruth changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; PetscFunctionBegin; ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); *flag = PETSC_TRUE; alpha = neP->alpha; maxstep = neP->maxstep; steptol = neP->steptol; ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr); if (!*ynorm) { ierr = PetscInfo(snes,"Search direction and size is 0\n");CHKERRQ(ierr); *gnorm = fnorm; ierr = VecCopy(x,w);CHKERRQ(ierr); ierr = VecCopy(f,g);CHKERRQ(ierr); *flag = PETSC_FALSE; goto theend1; } if (*ynorm > maxstep) { /* Step too big, so scale back */ scale = maxstep/(*ynorm); ierr = PetscInfo2(snes,"Scaling step by %G old ynorm %G\n",PetscRealPart(scale),*ynorm);CHKERRQ(ierr); ierr = VecScale(y,scale);CHKERRQ(ierr); *ynorm = maxstep; } ierr = VecMaxPointwiseDivide(y,x,&rellength);CHKERRQ(ierr); minlambda = steptol/rellength; ierr = MatMult(snes->jacobian,y,w);CHKERRQ(ierr); #if defined(PETSC_USE_COMPLEX) ierr = VecDot(f,w,&cinitslope);CHKERRQ(ierr); initslope = PetscRealPart(cinitslope); #else ierr = VecDot(f,w,&initslope);CHKERRQ(ierr); #endif if (initslope > 0.0) initslope = -initslope; if (initslope == 0.0) initslope = -1.0; ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); if (snes->nfuncs >= snes->max_funcs) { ierr = PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n");CHKERRQ(ierr); *flag = PETSC_FALSE; snes->reason = SNES_DIVERGED_FUNCTION_COUNT; goto theend1; } ierr = SNESComputeFunction(snes,w,g); if (PetscExceptionValue(ierr)) { PetscErrorCode pierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(pierr); } CHKERRQ(ierr); ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); if (*gnorm != *gnorm) SETERRQ(PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); if (.5*(*gnorm)*(*gnorm) <= .5*fnorm*fnorm + alpha*initslope) { /* Sufficient reduction */ ierr = PetscInfo(snes,"Using full step\n");CHKERRQ(ierr); goto theend1; } /* Fit points with quadratic */ lambda = 1.0; lambdatemp = -initslope/((*gnorm)*(*gnorm) - fnorm*fnorm - 2.0*initslope); lambdaprev = lambda; gnormprev = *gnorm; if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; if (lambdatemp <= .1*lambda) lambda = .1*lambda; else lambda = lambdatemp; #if defined(PETSC_USE_COMPLEX) clambda = lambda; ierr = VecWAXPY(w,-clambda,y,x);CHKERRQ(ierr); #else ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); #endif if (snes->nfuncs >= snes->max_funcs) { ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while attempting quadratic backtracking! %D \n",snes->nfuncs);CHKERRQ(ierr); *flag = PETSC_FALSE; snes->reason = SNES_DIVERGED_FUNCTION_COUNT; goto theend1; } ierr = SNESComputeFunction(snes,w,g); if (PetscExceptionValue(ierr)) { PetscErrorCode pierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(pierr); } CHKERRQ(ierr); ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); if (*gnorm != *gnorm) SETERRQ(PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*alpha*initslope) { /* sufficient reduction */ ierr = PetscInfo1(snes,"Quadratically determined step, lambda=%18.16e\n",lambda);CHKERRQ(ierr); goto theend1; } /* Fit points with cubic */ count = 1; while (count < 10000) { if (lambda <= minlambda) { /* bad luck; use full step */ ierr = PetscInfo1(snes,"Unable to find good step length! %D \n",count);CHKERRQ(ierr); ierr = PetscInfo5(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",fnorm,*gnorm,*ynorm,lambda,initslope);CHKERRQ(ierr); *flag = PETSC_FALSE; break; } t1 = .5*((*gnorm)*(*gnorm) - fnorm*fnorm) - lambda*initslope; t2 = .5*(gnormprev*gnormprev - fnorm*fnorm) - lambdaprev*initslope; a = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); b = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); d = b*b - 3*a*initslope; if (d < 0.0) d = 0.0; if (a == 0.0) { lambdatemp = -initslope/(2.0*b); } else { lambdatemp = (-b + sqrt(d))/(3.0*a); } lambdaprev = lambda; gnormprev = *gnorm; if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; if (lambdatemp <= .1*lambda) lambda = .1*lambda; else lambda = lambdatemp; #if defined(PETSC_USE_COMPLEX) clambda = lambda; ierr = VecWAXPY(w,-clambda,y,x);CHKERRQ(ierr); #else ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); #endif if (snes->nfuncs >= snes->max_funcs) { ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);CHKERRQ(ierr); ierr = PetscInfo5(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",fnorm,*gnorm,*ynorm,lambda,initslope);CHKERRQ(ierr); *flag = PETSC_FALSE; snes->reason = SNES_DIVERGED_FUNCTION_COUNT; break; } ierr = SNESComputeFunction(snes,w,g); if (PetscExceptionValue(ierr)) { PetscErrorCode pierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(pierr); } CHKERRQ(ierr); ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); if (*gnorm != *gnorm) SETERRQ(PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*alpha*initslope) { /* is reduction enough? */ ierr = PetscInfo1(snes,"Cubically determined step, lambda=%18.16e\n",lambda);CHKERRQ(ierr); break; } else { ierr = PetscInfo1(snes,"Cubic step no good, shrinking lambda, lambda=%18.16e\n",lambda);CHKERRQ(ierr); } count++; } if (count >= 10000) { SETERRQ(PETSC_ERR_LIB, "Lambda was decreased more than 10,000 times, so something is probably wrong with the function evaluation"); } theend1: /* Optional user-defined check for line search step validity */ if (neP->postcheckstep && *flag) { ierr = (*neP->postcheckstep)(snes,x,y,w,neP->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); if (changed_y) { ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); } if (changed_y || changed_w) { /* recompute the function if the step has changed */ ierr = SNESComputeFunction(snes,w,g); if (PetscExceptionValue(ierr)) { PetscErrorCode pierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(pierr); } CHKERRQ(ierr); ierr = VecNormBegin(g,NORM_2,gnorm);CHKERRQ(ierr); if (*gnorm != *gnorm) SETERRQ(PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); ierr = VecNormBegin(y,NORM_2,ynorm);CHKERRQ(ierr); ierr = VecNormEnd(g,NORM_2,gnorm);CHKERRQ(ierr); ierr = VecNormEnd(y,NORM_2,ynorm);CHKERRQ(ierr); } } ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); PetscFunctionReturn(0); } /* -------------------------------------------------------------------------- */ #undef __FUNCT__ #define __FUNCT__ "SNESLineSearchQuadratic" /*@C SNESLineSearchQuadratic - Performs a quadratic line search. Collective on SNES and Vec Input Parameters: + snes - the SNES context . lsctx - optional context for line search (not used here) . x - current iterate . f - residual evaluated at x . y - search direction . w - work vector - fnorm - 2-norm of f Output Parameters: + g - residual evaluated at new iterate w . w - new iterate (x + alpha*y) . gnorm - 2-norm of g . ynorm - 2-norm of search length - flag - PETSC_TRUE if line search succeeds; PETSC_FALSE on failure. Options Database Keys: + -snes_ls quadratic - Activates SNESLineSearchQuadratic() . -snes_ls_alpha - Sets alpha . -snes_ls_maxstep - Sets maxstep - -snes_ls_steptol - Sets steptol, this is the minimum step size that the line search code will accept; min p[i]/x[i] < steptol. The -snes_stol is the minimum step length the default convergence test will use and is based on 2-norm(p) < stol*2-norm(x) Notes: Use SNESLineSearchSet() to set this routine within the SNESLS method. Level: advanced .keywords: SNES, nonlinear, quadratic, line search .seealso: SNESLineSearchCubic(), SNESLineSearchNo(), SNESLineSearchSet(), SNESLineSearchNoNorms() @*/ PetscErrorCode PETSCSNES_DLLEXPORT SNESLineSearchQuadratic(SNES snes,void *lsctx,Vec x,Vec f,Vec g,Vec y,Vec w,PetscReal fnorm,PetscReal *ynorm,PetscReal *gnorm,PetscTruth *flag) { /* Note that for line search purposes we work with with the related minimization problem: min z(x): R^n -> R, where z(x) = .5 * fnorm*fnorm,and fnorm = || f ||_2. */ PetscReal steptol,initslope,maxstep,minlambda,alpha,lambda,lambdatemp,rellength; #if defined(PETSC_USE_COMPLEX) PetscScalar cinitslope,clambda; #endif PetscErrorCode ierr; PetscInt count; SNES_LS *neP = (SNES_LS*)snes->data; PetscScalar scale; PetscTruth changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; PetscFunctionBegin; ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); *flag = PETSC_TRUE; alpha = neP->alpha; maxstep = neP->maxstep; steptol = neP->steptol; ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr); if (*ynorm == 0.0) { ierr = PetscInfo(snes,"Search direction and size is 0\n");CHKERRQ(ierr); *gnorm = fnorm; ierr = VecCopy(x,w);CHKERRQ(ierr); ierr = VecCopy(f,g);CHKERRQ(ierr); *flag = PETSC_FALSE; goto theend2; } if (*ynorm > maxstep) { /* Step too big, so scale back */ scale = maxstep/(*ynorm); ierr = VecScale(y,scale);CHKERRQ(ierr); *ynorm = maxstep; } ierr = VecMaxPointwiseDivide(y,x,&rellength);CHKERRQ(ierr); minlambda = steptol/rellength; ierr = MatMult(snes->jacobian,y,w);CHKERRQ(ierr); #if defined(PETSC_USE_COMPLEX) ierr = VecDot(f,w,&cinitslope);CHKERRQ(ierr); initslope = PetscRealPart(cinitslope); #else ierr = VecDot(f,w,&initslope);CHKERRQ(ierr); #endif if (initslope > 0.0) initslope = -initslope; if (initslope == 0.0) initslope = -1.0; ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); if (snes->nfuncs >= snes->max_funcs) { ierr = PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n");CHKERRQ(ierr); *flag = PETSC_FALSE; snes->reason = SNES_DIVERGED_FUNCTION_COUNT; goto theend2; } ierr = SNESComputeFunction(snes,w,g); if (PetscExceptionValue(ierr)) { PetscErrorCode pierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(pierr); } CHKERRQ(ierr); ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); if (*gnorm != *gnorm) SETERRQ(PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); if (.5*(*gnorm)*(*gnorm) <= .5*fnorm*fnorm + alpha*initslope) { /* Sufficient reduction */ ierr = PetscInfo(snes,"Using full step\n");CHKERRQ(ierr); goto theend2; } /* Fit points with quadratic */ lambda = 1.0; count = 1; while (PETSC_TRUE) { if (lambda <= minlambda) { /* bad luck; use full step */ ierr = PetscInfo1(snes,"Unable to find good step length! %D \n",count);CHKERRQ(ierr); ierr = PetscInfo5(snes,"fnorm=%G, gnorm=%G, ynorm=%G, lambda=%G, initial slope=%G\n",fnorm,*gnorm,*ynorm,lambda,initslope);CHKERRQ(ierr); ierr = VecCopy(x,w);CHKERRQ(ierr); *flag = PETSC_FALSE; break; } lambdatemp = -initslope/((*gnorm)*(*gnorm) - fnorm*fnorm - 2.0*initslope); if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; if (lambdatemp <= .1*lambda) lambda = .1*lambda; else lambda = lambdatemp; #if defined(PETSC_USE_COMPLEX) clambda = lambda; ierr = VecWAXPY(w,-clambda,y,x);CHKERRQ(ierr); #else ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); #endif if (snes->nfuncs >= snes->max_funcs) { ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);CHKERRQ(ierr); ierr = PetscInfo5(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",fnorm,*gnorm,*ynorm,lambda,initslope);CHKERRQ(ierr); *flag = PETSC_FALSE; snes->reason = SNES_DIVERGED_FUNCTION_COUNT; break; } ierr = SNESComputeFunction(snes,w,g); if (PetscExceptionValue(ierr)) { PetscErrorCode pierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(pierr); } CHKERRQ(ierr); ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); if (*gnorm != *gnorm) SETERRQ(PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*alpha*initslope) { /* sufficient reduction */ ierr = PetscInfo1(snes,"Quadratically determined step, lambda=%G\n",lambda);CHKERRQ(ierr); break; } count++; } theend2: /* Optional user-defined check for line search step validity */ if (neP->postcheckstep) { ierr = (*neP->postcheckstep)(snes,x,y,w,neP->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); if (changed_y) { ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); } if (changed_y || changed_w) { /* recompute the function if the step has changed */ ierr = SNESComputeFunction(snes,w,g); if (PetscExceptionValue(ierr)) { PetscErrorCode pierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(pierr); } CHKERRQ(ierr); ierr = VecNormBegin(g,NORM_2,gnorm);CHKERRQ(ierr); ierr = VecNormBegin(y,NORM_2,ynorm);CHKERRQ(ierr); ierr = VecNormEnd(g,NORM_2,gnorm);CHKERRQ(ierr); ierr = VecNormEnd(y,NORM_2,ynorm);CHKERRQ(ierr); if (*gnorm != *gnorm) SETERRQ(PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); } } ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); PetscFunctionReturn(0); } /* -------------------------------------------------------------------------- */ #undef __FUNCT__ #define __FUNCT__ "SNESLineSearchSet" /*@C SNESLineSearchSet - Sets the line search routine to be used by the method SNESLS. Input Parameters: + snes - nonlinear context obtained from SNESCreate() . lsctx - optional user-defined context for use by line search - func - pointer to int function Collective on SNES Available Routines: + SNESLineSearchCubic() - default line search . SNESLineSearchQuadratic() - quadratic line search . SNESLineSearchNo() - the full Newton step (actually not a line search) - SNESLineSearchNoNorms() - the full Newton step (calculating no norms; faster in parallel) Options Database Keys: + -snes_ls [cubic,quadratic,basic,basicnonorms] - Selects line search . -snes_ls_alpha - Sets alpha . -snes_ls_maxstep - Sets maxstep - -snes_ls_steptol - Sets steptol, this is the minimum step size that the line search code will accept; min p[i]/x[i] < steptol. The -snes_stol is the minimum step length the default convergence test will use and is based on 2-norm(p) < stol*2-norm(x) Calling sequence of func: .vb func (SNES snes,void *lsctx,Vec x,Vec f,Vec g,Vec y,Vec w, PetscReal fnorm,PetscReal *ynorm,PetscReal *gnorm,PetscTruth *flag) .ve Input parameters for func: + snes - nonlinear context . lsctx - optional user-defined context for line search . x - current iterate . f - residual evaluated at x . y - search direction . w - work vector - fnorm - 2-norm of f Output parameters for func: + g - residual evaluated at new iterate y . w - new iterate . gnorm - 2-norm of g . ynorm - 2-norm of search length - flag - set to PETSC_TRUE if the line search succeeds; PETSC_FALSE on failure. Level: advanced .keywords: SNES, nonlinear, set, line search, routine .seealso: SNESLineSearchCubic(), SNESLineSearchQuadratic(), SNESLineSearchNo(), SNESLineSearchNoNorms(), SNESLineSearchSetPostCheck(), SNESLineSearchSetParams(), SNESLineSearchGetParams(), SNESLineSearchSetPreCheck() @*/ PetscErrorCode PETSCSNES_DLLEXPORT SNESLineSearchSet(SNES snes,PetscErrorCode (*func)(SNES,void*,Vec,Vec,Vec,Vec,Vec,PetscReal,PetscReal*,PetscReal*,PetscTruth*),void *lsctx) { PetscErrorCode ierr,(*f)(SNES,PetscErrorCode (*)(SNES,void*,Vec,Vec,Vec,Vec,Vec,PetscReal,PetscReal*,PetscReal*,PetscTruth*),void*); PetscFunctionBegin; ierr = PetscObjectQueryFunction((PetscObject)snes,"SNESLineSearchSet_C",(void (**)(void))&f);CHKERRQ(ierr); if (f) { ierr = (*f)(snes,func,lsctx);CHKERRQ(ierr); } PetscFunctionReturn(0); } typedef PetscErrorCode (*FCN2)(SNES,void*,Vec,Vec,Vec,Vec,Vec,PetscReal,PetscReal*,PetscReal*,PetscTruth*); /* force argument to next function to not be extern C*/ /* -------------------------------------------------------------------------- */ EXTERN_C_BEGIN #undef __FUNCT__ #define __FUNCT__ "SNESLineSearchSet_LS" PetscErrorCode PETSCSNES_DLLEXPORT SNESLineSearchSet_LS(SNES snes,FCN2 func,void *lsctx) { PetscFunctionBegin; ((SNES_LS *)(snes->data))->LineSearch = func; ((SNES_LS *)(snes->data))->lsP = lsctx; PetscFunctionReturn(0); } EXTERN_C_END /* -------------------------------------------------------------------------- */ #undef __FUNCT__ #define __FUNCT__ "SNESLineSearchSetPostCheck" /*@C SNESLineSearchSetPostCheck - Sets a routine to check the validity of new iterate computed by the line search routine in the Newton-based method SNESLS. Input Parameters: + snes - nonlinear context obtained from SNESCreate() . func - pointer to function - checkctx - optional user-defined context for use by step checking routine Collective on SNES Calling sequence of func: .vb int func (SNES snes, Vec x,Vec y,Vec w,void *checkctx, PetscTruth *changed_y,PetscTruth *changed_w) .ve where func returns an error code of 0 on success and a nonzero on failure. Input parameters for func: + snes - nonlinear context . checkctx - optional user-defined context for use by step checking routine . x - previous iterate . y - new search direction and length - w - current candidate iterate Output parameters for func: + y - search direction (possibly changed) . w - current iterate (possibly modified) . changed_y - indicates search direction was changed by this routine - changed_w - indicates current iterate was changed by this routine Level: advanced Notes: All line searches accept the new iterate computed by the line search checking routine. Only one of changed_y and changed_w can be PETSC_TRUE On input w = x + y SNESLineSearchNo() and SNESLineSearchNoNorms() (1) compute a candidate iterate u_{i+1}, (2) pass control to the checking routine, and then (3) compute the corresponding nonlinear function f(u_{i+1}) with the (possibly altered) iterate u_{i+1}. SNESLineSearchQuadratic() and SNESLineSearchCubic() (1) compute a candidate iterate u_{i+1} as well as a candidate nonlinear function f(u_{i+1}), (2) pass control to the checking routine, and then (3) force a re-evaluation of f(u_{i+1}) if any changes were made to the candidate iterate in the checking routine (as indicated by flag=PETSC_TRUE). The overhead of this extra function re-evaluation can be very costly, so use this feature with caution! .keywords: SNES, nonlinear, set, line search check, step check, routine .seealso: SNESLineSearchSet(), SNESLineSearchSetPreCheck() @*/ PetscErrorCode PETSCSNES_DLLEXPORT SNESLineSearchSetPostCheck(SNES snes,PetscErrorCode (*func)(SNES,Vec,Vec,Vec,void*,PetscTruth*,PetscTruth*),void *checkctx) { PetscErrorCode ierr,(*f)(SNES,PetscErrorCode (*)(SNES,Vec,Vec,Vec,void*,PetscTruth*,PetscTruth*),void*); PetscFunctionBegin; ierr = PetscObjectQueryFunction((PetscObject)snes,"SNESLineSearchSetPostCheck_C",(void (**)(void))&f);CHKERRQ(ierr); if (f) { ierr = (*f)(snes,func,checkctx);CHKERRQ(ierr); } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "SNESLineSearchSetPreCheck" /*@C SNESLineSearchSetPreCheck - Sets a routine to check the validity of a new direction given by the linear solve before the line search is called. Input Parameters: + snes - nonlinear context obtained from SNESCreate() . func - pointer to function - checkctx - optional user-defined context for use by step checking routine Collective on SNES Calling sequence of func: .vb int func (SNES snes, Vec x,Vec y,void *checkctx, PetscTruth *changed_y) .ve where func returns an error code of 0 on success and a nonzero on failure. Input parameters for func: + snes - nonlinear context . checkctx - optional user-defined context for use by step checking routine . x - previous iterate - y - new search direction and length Output parameters for func: + y - search direction (possibly changed) - changed_y - indicates search direction was changed by this routine Level: advanced Notes: All line searches accept the new iterate computed by the line search checking routine. .keywords: SNES, nonlinear, set, line search check, step check, routine .seealso: SNESLineSearchSet(), SNESLineSearchSetPostCheck(), SNESSetUpdate() @*/ PetscErrorCode PETSCSNES_DLLEXPORT SNESLineSearchSetPreCheck(SNES snes,PetscErrorCode (*func)(SNES,Vec,Vec,void*,PetscTruth*),void *checkctx) { PetscErrorCode ierr,(*f)(SNES,PetscErrorCode (*)(SNES,Vec,Vec,void*,PetscTruth*),void*); PetscFunctionBegin; ierr = PetscObjectQueryFunction((PetscObject)snes,"SNESLineSearchSetPreCheck_C",(void (**)(void))&f);CHKERRQ(ierr); if (f) { ierr = (*f)(snes,func,checkctx);CHKERRQ(ierr); } PetscFunctionReturn(0); } /* -------------------------------------------------------------------------- */ typedef PetscErrorCode (*FCN1)(SNES,Vec,Vec,Vec,void*,PetscTruth*,PetscTruth*); /* force argument to next function to not be extern C*/ typedef PetscErrorCode (*FCN3)(SNES,Vec,Vec,void*,PetscTruth*); /* force argument to next function to not be extern C*/ EXTERN_C_BEGIN #undef __FUNCT__ #define __FUNCT__ "SNESLineSearchSetPostCheck_LS" PetscErrorCode PETSCSNES_DLLEXPORT SNESLineSearchSetPostCheck_LS(SNES snes,FCN1 func,void *checkctx) { PetscFunctionBegin; ((SNES_LS *)(snes->data))->postcheckstep = func; ((SNES_LS *)(snes->data))->postcheck = checkctx; PetscFunctionReturn(0); } EXTERN_C_END EXTERN_C_BEGIN #undef __FUNCT__ #define __FUNCT__ "SNESLineSearchSetPreCheck_LS" PetscErrorCode PETSCSNES_DLLEXPORT SNESLineSearchSetPreCheck_LS(SNES snes,FCN3 func,void *checkctx) { PetscFunctionBegin; ((SNES_LS *)(snes->data))->precheckstep = func; ((SNES_LS *)(snes->data))->precheck = checkctx; PetscFunctionReturn(0); } EXTERN_C_END /* SNESView_LS - Prints info from the SNESLS data structure. Input Parameters: . SNES - the SNES context . viewer - visualization context Application Interface Routine: SNESView() */ #undef __FUNCT__ #define __FUNCT__ "SNESView_LS" static PetscErrorCode SNESView_LS(SNES snes,PetscViewer viewer) { SNES_LS *ls = (SNES_LS *)snes->data; const char *cstr; PetscErrorCode ierr; PetscTruth iascii; PetscFunctionBegin; ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); if (iascii) { if (ls->LineSearch == SNESLineSearchNo) cstr = "SNESLineSearchNo"; else if (ls->LineSearch == SNESLineSearchQuadratic) cstr = "SNESLineSearchQuadratic"; else if (ls->LineSearch == SNESLineSearchCubic) cstr = "SNESLineSearchCubic"; else cstr = "unknown"; ierr = PetscViewerASCIIPrintf(viewer," line search variant: %s\n",cstr);CHKERRQ(ierr); ierr = PetscViewerASCIIPrintf(viewer," alpha=%G, maxstep=%G, steptol=%G\n",ls->alpha,ls->maxstep,ls->steptol);CHKERRQ(ierr); } else { SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for SNES EQ LS",((PetscObject)viewer)->type_name); } PetscFunctionReturn(0); } /* -------------------------------------------------------------------------- */ /* SNESSetFromOptions_LS - Sets various parameters for the SNESLS method. Input Parameter: . snes - the SNES context Application Interface Routine: SNESSetFromOptions() */ #undef __FUNCT__ #define __FUNCT__ "SNESSetFromOptions_LS" static PetscErrorCode SNESSetFromOptions_LS(SNES snes) { SNES_LS *ls = (SNES_LS *)snes->data; const char *lses[] = {"basic","basicnonorms","quadratic","cubic"}; PetscErrorCode ierr; PetscInt indx; PetscTruth flg; PetscFunctionBegin; ierr = PetscOptionsHead("SNES Line search options");CHKERRQ(ierr); ierr = PetscOptionsReal("-snes_ls_alpha","Function norm must decrease by","None",ls->alpha,&ls->alpha,0);CHKERRQ(ierr); ierr = PetscOptionsReal("-snes_ls_maxstep","Step must be less than","None",ls->maxstep,&ls->maxstep,0);CHKERRQ(ierr); ierr = PetscOptionsReal("-snes_ls_steptol","Step must be greater than","None",ls->steptol,&ls->steptol,0);CHKERRQ(ierr); ierr = PetscOptionsEList("-snes_ls","Line search used","SNESLineSearchSet",lses,4,"cubic",&indx,&flg);CHKERRQ(ierr); if (flg) { switch (indx) { case 0: ierr = SNESLineSearchSet(snes,SNESLineSearchNo,PETSC_NULL);CHKERRQ(ierr); break; case 1: ierr = SNESLineSearchSet(snes,SNESLineSearchNoNorms,PETSC_NULL);CHKERRQ(ierr); break; case 2: ierr = SNESLineSearchSet(snes,SNESLineSearchQuadratic,PETSC_NULL);CHKERRQ(ierr); break; case 3: ierr = SNESLineSearchSet(snes,SNESLineSearchCubic,PETSC_NULL);CHKERRQ(ierr); break; } } ierr = PetscOptionsTail();CHKERRQ(ierr); PetscFunctionReturn(0); } /* -------------------------------------------------------------------------- */ #undef __FUNCT__ #define __FUNCT__ "SNESConverged_LS" /*@C SNESConverged_LS - Monitors the convergence of the line search method SNESLS for solving systems of nonlinear equations (default). Collective on SNES Input Parameters: + snes - the SNES context . it - the iteration (0 indicates before any Newton steps) . xnorm - 2-norm of current iterate . pnorm - 2-norm of current step . fnorm - 2-norm of function at current iterate - dummy - unused context Output Parameter: . reason - one of $ SNES_CONVERGED_FNORM_ABS - (fnorm < abstol), $ SNES_CONVERGED_PNORM_RELATIVE - (pnorm < xtol*xnorm), $ SNES_CONVERGED_FNORM_RELATIVE - (fnorm < rtol*fnorm0), $ SNES_DIVERGED_FUNCTION_COUNT - (nfct > maxf), $ SNES_DIVERGED_FNORM_NAN - (fnorm == NaN), $ SNES_CONVERGED_ITERATING - (otherwise), where + maxf - maximum number of function evaluations, set with SNESSetTolerances() . nfct - number of function evaluations, . abstol - absolute function norm tolerance, set with SNESSetTolerances() - rtol - relative function norm tolerance, set with SNESSetTolerances() Level: intermediate .keywords: SNES, nonlinear, default, converged, convergence .seealso: SNESSetConvergenceTest() @*/ PetscErrorCode PETSCSNES_DLLEXPORT SNESConverged_LS(SNES snes,PetscInt it,PetscReal xnorm,PetscReal pnorm,PetscReal fnorm,SNESConvergedReason *reason,void *dummy) { PetscErrorCode ierr; PetscFunctionBegin; PetscValidHeaderSpecific(snes,SNES_COOKIE,1); PetscValidType(snes,1); PetscValidPointer(reason,6); ierr = SNESDefaultConverged(snes,it,xnorm,pnorm,fnorm,reason,dummy);CHKERRQ(ierr); PetscFunctionReturn(0); } /* -------------------------------------------------------------------------- */ /*MC SNESLS - Newton based nonlinear solver that uses a line search Options Database: + -snes_ls [cubic,quadratic,basic,basicnonorms] - Selects line search . -snes_ls_alpha - Sets alpha . -snes_ls_maxstep - Sets maxstep - -snes_ls_steptol - Sets steptol, this is the minimum step size that the line search code will accept; min p[i]/x[i] < steptol. The -snes_stol is the minimum step length the default convergence test will use and is based on 2-norm(p) < stol*2-norm(x) Notes: This is the default nonlinear solver in SNES Level: beginner .seealso: SNESCreate(), SNES, SNESSetType(), SNESTR, SNESLineSearchSet(), SNESLineSearchSetPostCheck(), SNESLineSearchNo(), SNESLineSearchCubic(), SNESLineSearchQuadratic(), SNESLineSearchSet(), SNESLineSearchNoNorms(), SNESLineSearchSetPreCheck() M*/ EXTERN_C_BEGIN #undef __FUNCT__ #define __FUNCT__ "SNESCreate_LS" PetscErrorCode PETSCSNES_DLLEXPORT SNESCreate_LS(SNES snes) { PetscErrorCode ierr; SNES_LS *neP; PetscFunctionBegin; snes->ops->setup = SNESSetUp_LS; snes->ops->solve = SNESSolve_LS; snes->ops->destroy = SNESDestroy_LS; snes->ops->converged = SNESConverged_LS; snes->ops->setfromoptions = SNESSetFromOptions_LS; snes->ops->view = SNESView_LS; snes->nwork = 0; ierr = PetscNew(SNES_LS,&neP);CHKERRQ(ierr); ierr = PetscLogObjectMemory(snes,sizeof(SNES_LS));CHKERRQ(ierr); snes->data = (void*)neP; neP->alpha = 1.e-4; neP->maxstep = 1.e8; neP->steptol = 1.e-12; neP->LineSearch = SNESLineSearchCubic; neP->lsP = PETSC_NULL; neP->postcheckstep = PETSC_NULL; neP->postcheck = PETSC_NULL; neP->precheckstep = PETSC_NULL; neP->precheck = PETSC_NULL; ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSet_C", "SNESLineSearchSet_LS", SNESLineSearchSet_LS);CHKERRQ(ierr); ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetPostCheck_C", "SNESLineSearchSetPostCheck_LS", SNESLineSearchSetPostCheck_LS);CHKERRQ(ierr); ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetPreCheck_C", "SNESLineSearchSetPreCheck_LS", SNESLineSearchSetPreCheck_LS);CHKERRQ(ierr); PetscFunctionReturn(0); } EXTERN_C_END