13d443117SPeter Brune #include <petsc-private/linesearchimpl.h> /*I "petscsnes.h" I*/ 2b45d2f2cSJed Brown #include <petsc-private/snesimpl.h> 371b4ebd2SPeter Brune 471b4ebd2SPeter Brune typedef struct { 571b4ebd2SPeter Brune PetscReal alpha; /* sufficient decrease parameter */ 6f1c6b773SPeter Brune } SNESLineSearch_BT; 771b4ebd2SPeter Brune 871b4ebd2SPeter Brune #undef __FUNCT__ 92f4102e2SPeter Brune #define __FUNCT__ "SNESLineSearchBTSetAlpha" 102f4102e2SPeter Brune /*@ 112f4102e2SPeter Brune SNESLineSearchBTSetAlpha - Sets the descent parameter, alpha, in the BT linesearch variant. 122f4102e2SPeter Brune 132f4102e2SPeter Brune Input Parameters: 142f4102e2SPeter Brune + linesearch - linesearch context 152f4102e2SPeter Brune - alpha - The descent parameter 162f4102e2SPeter Brune 172f4102e2SPeter Brune Level: intermediate 182f4102e2SPeter Brune 191a4f838cSPeter Brune .seealso: SNESLineSearchSetLambda(), SNESLineSearchGetTolerances() SNESLINESEARCHBT 202f4102e2SPeter Brune @*/ 212f4102e2SPeter Brune PetscErrorCode SNESLineSearchBTSetAlpha(SNESLineSearch linesearch, PetscReal alpha) 222f4102e2SPeter Brune { 232f4102e2SPeter Brune SNESLineSearch_BT *bt; 242f4102e2SPeter Brune PetscFunctionBegin; 252f4102e2SPeter Brune PetscValidHeaderSpecific(linesearch,SNESLINESEARCH_CLASSID,1); 262f4102e2SPeter Brune bt = (SNESLineSearch_BT *)linesearch->data; 272f4102e2SPeter Brune bt->alpha = alpha; 282f4102e2SPeter Brune PetscFunctionReturn(0); 292f4102e2SPeter Brune } 302f4102e2SPeter Brune 312f4102e2SPeter Brune 322f4102e2SPeter Brune #undef __FUNCT__ 332f4102e2SPeter Brune #define __FUNCT__ "SNESLineSearchBTGetAlpha" 342f4102e2SPeter Brune /*@ 352f4102e2SPeter Brune SNESLineSearchBTGetAlpha - Gets the descent parameter, alpha, in the BT linesearch variant. 362f4102e2SPeter Brune 372f4102e2SPeter Brune Input Parameters: 382f4102e2SPeter Brune . linesearch - linesearch context 398e557f58SPeter Brune 408e557f58SPeter Brune Output Parameters: 412f4102e2SPeter Brune . alpha - The descent parameter 422f4102e2SPeter Brune 432f4102e2SPeter Brune Level: intermediate 442f4102e2SPeter Brune 451a4f838cSPeter Brune .seealso: SNESLineSearchGetLambda(), SNESLineSearchGetTolerances() SNESLINESEARCHBT 462f4102e2SPeter Brune @*/ 472f4102e2SPeter Brune PetscErrorCode SNESLineSearchBTGetAlpha(SNESLineSearch linesearch, PetscReal *alpha) 482f4102e2SPeter Brune { 492f4102e2SPeter Brune SNESLineSearch_BT *bt; 502f4102e2SPeter Brune PetscFunctionBegin; 512f4102e2SPeter Brune PetscValidHeaderSpecific(linesearch,SNESLINESEARCH_CLASSID,1); 522f4102e2SPeter Brune bt = (SNESLineSearch_BT *)linesearch->data; 532f4102e2SPeter Brune *alpha = bt->alpha; 542f4102e2SPeter Brune PetscFunctionReturn(0); 552f4102e2SPeter Brune } 562f4102e2SPeter Brune 572f4102e2SPeter Brune #undef __FUNCT__ 58f1c6b773SPeter Brune #define __FUNCT__ "SNESLineSearchApply_BT" 59f1c6b773SPeter Brune static PetscErrorCode SNESLineSearchApply_BT(SNESLineSearch linesearch) 6071b4ebd2SPeter Brune { 6171b4ebd2SPeter Brune PetscBool changed_y,changed_w; 6271b4ebd2SPeter Brune PetscErrorCode ierr; 6371b4ebd2SPeter Brune Vec X,F,Y,W,G; 6471b4ebd2SPeter Brune SNES snes; 6571b4ebd2SPeter Brune PetscReal fnorm, xnorm, ynorm, gnorm, gnormprev; 663b288129SPeter Brune PetscReal lambda,lambdatemp,lambdaprev,minlambda,maxstep,initslope,alpha,stol; 67dfcd3828SPeter Brune PetscReal t1,t2,a,b,d; 6871b4ebd2SPeter Brune #if defined(PETSC_USE_COMPLEX) 6971b4ebd2SPeter Brune PetscScalar cinitslope; 7071b4ebd2SPeter Brune #endif 7171b4ebd2SPeter Brune PetscBool domainerror; 7271b4ebd2SPeter Brune PetscViewer monitor; 7371b4ebd2SPeter Brune PetscInt max_its,count; 74f1c6b773SPeter Brune SNESLineSearch_BT *bt; 7571b4ebd2SPeter Brune Mat jac; 7671b4ebd2SPeter Brune 7771b4ebd2SPeter Brune 7871b4ebd2SPeter Brune PetscFunctionBegin; 7971b4ebd2SPeter Brune 80f1c6b773SPeter Brune ierr = SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, &G);CHKERRQ(ierr); 81f1c6b773SPeter Brune ierr = SNESLineSearchGetNorms(linesearch, &xnorm, &fnorm, &ynorm);CHKERRQ(ierr); 82f1c6b773SPeter Brune ierr = SNESLineSearchGetLambda(linesearch, &lambda);CHKERRQ(ierr); 83f1c6b773SPeter Brune ierr = SNESLineSearchGetSNES(linesearch, &snes);CHKERRQ(ierr); 84f1c6b773SPeter Brune ierr = SNESLineSearchGetMonitor(linesearch, &monitor);CHKERRQ(ierr); 85dfcd3828SPeter Brune ierr = SNESLineSearchGetTolerances(linesearch,&minlambda,&maxstep,PETSC_NULL,PETSC_NULL,PETSC_NULL,&max_its); 863b288129SPeter Brune ierr = SNESGetTolerances(snes,PETSC_NULL,PETSC_NULL,&stol,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 87f1c6b773SPeter Brune bt = (SNESLineSearch_BT *)linesearch->data; 8871b4ebd2SPeter Brune 8971b4ebd2SPeter Brune alpha = bt->alpha; 9071b4ebd2SPeter Brune 9171b4ebd2SPeter Brune ierr = SNESGetJacobian(snes, &jac, PETSC_NULL, PETSC_NULL, PETSC_NULL);CHKERRQ(ierr); 92c69d1a72SBarry Smith if (!jac) SETERRQ(((PetscObject)linesearch)->comm, PETSC_ERR_USER, "SNESLineSearchBT requires a Jacobian matrix"); 93c69d1a72SBarry Smith 9471b4ebd2SPeter Brune /* precheck */ 957b1df9c1SPeter Brune ierr = SNESLineSearchPreCheck(linesearch,X,Y,&changed_y);CHKERRQ(ierr); 96f1c6b773SPeter Brune ierr = SNESLineSearchSetSuccess(linesearch, PETSC_TRUE);CHKERRQ(ierr); 9771b4ebd2SPeter Brune 983b288129SPeter Brune ierr = VecNormBegin(Y, NORM_2, &ynorm);CHKERRQ(ierr); 993b288129SPeter Brune ierr = VecNormBegin(X, NORM_2, &xnorm);CHKERRQ(ierr); 1003b288129SPeter Brune ierr = VecNormEnd(Y, NORM_2, &ynorm);CHKERRQ(ierr); 1013b288129SPeter Brune ierr = VecNormEnd(X, NORM_2, &xnorm);CHKERRQ(ierr); 1023b288129SPeter Brune 10371b4ebd2SPeter Brune if (ynorm == 0.0) { 10471b4ebd2SPeter Brune if (monitor) { 105ea5d4fccSPeter Brune ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 10671b4ebd2SPeter Brune ierr = PetscViewerASCIIPrintf(monitor," Line search: Initial direction and size is 0\n");CHKERRQ(ierr); 107ea5d4fccSPeter Brune ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 10871b4ebd2SPeter Brune } 10971b4ebd2SPeter Brune ierr = VecCopy(X,W);CHKERRQ(ierr); 11071b4ebd2SPeter Brune ierr = VecCopy(F,G);CHKERRQ(ierr); 111f1c6b773SPeter Brune ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr); 11271b4ebd2SPeter Brune PetscFunctionReturn(0); 11371b4ebd2SPeter Brune } 11471b4ebd2SPeter Brune if (ynorm > maxstep) { /* Step too big, so scale back */ 11571b4ebd2SPeter Brune if (monitor) { 116ea5d4fccSPeter Brune ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 117c69d1a72SBarry Smith ierr = PetscViewerASCIIPrintf(monitor," Line search: Scaling step by %14.12e old ynorm %14.12e\n", (double)(maxstep/ynorm),(double)ynorm);CHKERRQ(ierr); 118ea5d4fccSPeter Brune ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 11971b4ebd2SPeter Brune } 12071b4ebd2SPeter Brune ierr = VecScale(Y,maxstep/(ynorm));CHKERRQ(ierr); 12171b4ebd2SPeter Brune ynorm = maxstep; 12271b4ebd2SPeter Brune } 12371b4ebd2SPeter Brune ierr = MatMult(jac,Y,W);CHKERRQ(ierr); 12471b4ebd2SPeter Brune #if defined(PETSC_USE_COMPLEX) 12571b4ebd2SPeter Brune ierr = VecDot(F,W,&cinitslope);CHKERRQ(ierr); 12671b4ebd2SPeter Brune initslope = PetscRealPart(cinitslope); 12771b4ebd2SPeter Brune #else 12871b4ebd2SPeter Brune ierr = VecDot(F,W,&initslope);CHKERRQ(ierr); 12971b4ebd2SPeter Brune #endif 13071b4ebd2SPeter Brune if (initslope > 0.0) initslope = -initslope; 13171b4ebd2SPeter Brune if (initslope == 0.0) initslope = -1.0; 13271b4ebd2SPeter Brune 133c427506bSPeter Brune ierr = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr); 1349bd66eb0SPeter Brune if (linesearch->ops->viproject) { 1359bd66eb0SPeter Brune ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr); 1369bd66eb0SPeter Brune } 13771b4ebd2SPeter Brune if (snes->nfuncs >= snes->max_funcs) { 13871b4ebd2SPeter Brune ierr = PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n");CHKERRQ(ierr); 13971b4ebd2SPeter Brune snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 140f1c6b773SPeter Brune ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr); 14171b4ebd2SPeter Brune PetscFunctionReturn(0); 14271b4ebd2SPeter Brune } 14371b4ebd2SPeter Brune ierr = SNESComputeFunction(snes,W,G);CHKERRQ(ierr); 144c427506bSPeter Brune ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr); 145c427506bSPeter Brune if (domainerror) { 146f1c6b773SPeter Brune ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr); 14771b4ebd2SPeter Brune PetscFunctionReturn(0); 14871b4ebd2SPeter Brune } 1499bd66eb0SPeter Brune if (linesearch->ops->vinorm) { 1509bd66eb0SPeter Brune gnorm = fnorm; 1519bd66eb0SPeter Brune ierr = (*linesearch->ops->vinorm)(snes, G, W, &gnorm);CHKERRQ(ierr); 1529bd66eb0SPeter Brune } else { 15371b4ebd2SPeter Brune ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr); 1549bd66eb0SPeter Brune } 1559bd66eb0SPeter Brune 15671b4ebd2SPeter Brune if (PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 157c69d1a72SBarry Smith ierr = PetscInfo2(snes,"Initial fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm);CHKERRQ(ierr); 158782c0b7fSBarry Smith if (.5*gnorm*gnorm <= .5*fnorm*fnorm + lambda*alpha*initslope) { /* Sufficient reduction or step tolerance convergence */ 15971b4ebd2SPeter Brune if (monitor) { 160ea5d4fccSPeter Brune ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 161c69d1a72SBarry Smith ierr = PetscViewerASCIIPrintf(monitor," Line search: Using full step: fnorm %14.12e gnorm %14.12e\n", (double)fnorm, (double)gnorm);CHKERRQ(ierr); 162ea5d4fccSPeter Brune ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 16371b4ebd2SPeter Brune } 16471b4ebd2SPeter Brune } else { 165*c21ba15cSPeter Brune /* Since the full step didn't work and the step is tiny, quit */ 166*c21ba15cSPeter Brune if (stol*xnorm > ynorm) { 167*c21ba15cSPeter Brune ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr); 168*c21ba15cSPeter Brune ierr = PetscInfo2(monitor,"Aborted due to ynorm < stol*xnorm (%14.12e < %14.12e) and inadequate full step.\n",ynorm,stol*xnorm);CHKERRQ(ierr); 169*c21ba15cSPeter Brune PetscFunctionReturn(0); 170*c21ba15cSPeter Brune } 17171b4ebd2SPeter Brune /* Fit points with quadratic */ 17271b4ebd2SPeter Brune lambdatemp = -initslope/(gnorm*gnorm - fnorm*fnorm - 2.0*lambda*initslope); 17371b4ebd2SPeter Brune lambdaprev = lambda; 17471b4ebd2SPeter Brune gnormprev = gnorm; 17571b4ebd2SPeter Brune if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 17671b4ebd2SPeter Brune if (lambdatemp <= .1*lambda) lambda = .1*lambda; 17771b4ebd2SPeter Brune else lambda = lambdatemp; 17871b4ebd2SPeter Brune 17971b4ebd2SPeter Brune ierr = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr); 1809bd66eb0SPeter Brune if (linesearch->ops->viproject) { 1819bd66eb0SPeter Brune ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr); 1829bd66eb0SPeter Brune } 18371b4ebd2SPeter Brune if (snes->nfuncs >= snes->max_funcs) { 18471b4ebd2SPeter Brune ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while attempting quadratic backtracking! %D \n",snes->nfuncs);CHKERRQ(ierr); 18571b4ebd2SPeter Brune snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 186f1c6b773SPeter Brune ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr); 18771b4ebd2SPeter Brune PetscFunctionReturn(0); 18871b4ebd2SPeter Brune } 18971b4ebd2SPeter Brune ierr = SNESComputeFunction(snes,W,G);CHKERRQ(ierr); 1906188f407SPeter Brune ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr); 1916188f407SPeter Brune if (domainerror) { 19271b4ebd2SPeter Brune PetscFunctionReturn(0); 19371b4ebd2SPeter Brune } 1949bd66eb0SPeter Brune if (linesearch->ops->vinorm) { 1959bd66eb0SPeter Brune gnorm = fnorm; 1969bd66eb0SPeter Brune ierr = (*linesearch->ops->vinorm)(snes, G, W, &gnorm);CHKERRQ(ierr); 1979bd66eb0SPeter Brune } else { 19871b4ebd2SPeter Brune ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr); 1999bd66eb0SPeter Brune } 20071b4ebd2SPeter Brune if (PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 20171b4ebd2SPeter Brune if (monitor) { 202ea5d4fccSPeter Brune ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 203c69d1a72SBarry Smith ierr = PetscViewerASCIIPrintf(monitor," Line search: gnorm after quadratic fit %14.12e\n",(double)gnorm);CHKERRQ(ierr); 204ea5d4fccSPeter Brune ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 20571b4ebd2SPeter Brune } 20671b4ebd2SPeter Brune if (.5*gnorm*gnorm < .5*fnorm*fnorm + lambda*alpha*initslope) { /* sufficient reduction */ 20771b4ebd2SPeter Brune if (monitor) { 208ea5d4fccSPeter Brune ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 20971b4ebd2SPeter Brune ierr = PetscViewerASCIIPrintf(monitor," Line search: Quadratically determined step, lambda=%18.16e\n",(double)lambda);CHKERRQ(ierr); 210ea5d4fccSPeter Brune ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 21171b4ebd2SPeter Brune } 21271b4ebd2SPeter Brune } else { 21371b4ebd2SPeter Brune /* Fit points with cubic */ 21471b4ebd2SPeter Brune for (count = 0; count < max_its; count++) { 21571b4ebd2SPeter Brune if (lambda <= minlambda) { 21671b4ebd2SPeter Brune if (monitor) { 217ea5d4fccSPeter Brune ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 21871b4ebd2SPeter Brune ierr = PetscViewerASCIIPrintf(monitor," Line search: unable to find good step length! After %D tries \n",count);CHKERRQ(ierr); 21971b4ebd2SPeter Brune ierr = PetscViewerASCIIPrintf(monitor, 22071b4ebd2SPeter Brune " Line search: fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, minlambda=%18.16e, lambda=%18.16e, initial slope=%18.16e\n", 221c69d1a72SBarry Smith (double)fnorm, (double)gnorm, (double)ynorm, (double)minlambda, (double)lambda, (double)initslope);CHKERRQ(ierr); 222ea5d4fccSPeter Brune ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 22371b4ebd2SPeter Brune } 224f1c6b773SPeter Brune ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr); 22571b4ebd2SPeter Brune PetscFunctionReturn(0); 22671b4ebd2SPeter Brune } 227b000cd8dSPeter Brune if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) { 22871b4ebd2SPeter Brune t1 = .5*(gnorm*gnorm - fnorm*fnorm) - lambda*initslope; 22971b4ebd2SPeter Brune t2 = .5*(gnormprev*gnormprev - fnorm*fnorm) - lambdaprev*initslope; 23071b4ebd2SPeter Brune a = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 23171b4ebd2SPeter Brune b = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 23271b4ebd2SPeter Brune d = b*b - 3*a*initslope; 23371b4ebd2SPeter Brune if (d < 0.0) d = 0.0; 23471b4ebd2SPeter Brune if (a == 0.0) { 23571b4ebd2SPeter Brune lambdatemp = -initslope/(2.0*b); 23671b4ebd2SPeter Brune } else { 23771b4ebd2SPeter Brune lambdatemp = (-b + PetscSqrtReal(d))/(3.0*a); 23871b4ebd2SPeter Brune } 239b000cd8dSPeter Brune } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) { 24071b4ebd2SPeter Brune lambdatemp = -initslope/(gnorm*gnorm - fnorm*fnorm - 2.0*initslope); 24159405d5eSPeter Brune } else { 24259405d5eSPeter Brune SETERRQ(((PetscObject)linesearch)->comm, PETSC_ERR_SUP, "unsupported line search order for type bt"); 24371b4ebd2SPeter Brune } 24471b4ebd2SPeter Brune lambdaprev = lambda; 24571b4ebd2SPeter Brune gnormprev = gnorm; 24671b4ebd2SPeter Brune if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 24771b4ebd2SPeter Brune if (lambdatemp <= .1*lambda) lambda = .1*lambda; 24871b4ebd2SPeter Brune else lambda = lambdatemp; 24971b4ebd2SPeter Brune ierr = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr); 250b7b2e573SPeter Brune if (linesearch->ops->viproject) { 251b7b2e573SPeter Brune ierr = (*linesearch->ops->viproject)(snes,W);CHKERRQ(ierr); 252b7b2e573SPeter Brune } 25371b4ebd2SPeter Brune if (snes->nfuncs >= snes->max_funcs) { 25471b4ebd2SPeter Brune ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);CHKERRQ(ierr); 25571b4ebd2SPeter Brune ierr = PetscInfo5(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lambda=%18.16e, initial slope=%18.16e\n", 256c69d1a72SBarry Smith (double)fnorm,(double)gnorm,(double)ynorm,(double)lambda,(double)initslope);CHKERRQ(ierr); 257f1c6b773SPeter Brune ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr); 25871b4ebd2SPeter Brune snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 25971b4ebd2SPeter Brune PetscFunctionReturn(0); 26071b4ebd2SPeter Brune } 26171b4ebd2SPeter Brune ierr = SNESComputeFunction(snes,W,G);CHKERRQ(ierr); 262c427506bSPeter Brune ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr); 263c427506bSPeter Brune if (domainerror) { 264f1c6b773SPeter Brune ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr); 26571b4ebd2SPeter Brune PetscFunctionReturn(0); 26671b4ebd2SPeter Brune } 2679bd66eb0SPeter Brune if (linesearch->ops->vinorm) { 2689bd66eb0SPeter Brune gnorm = fnorm; 2699bd66eb0SPeter Brune ierr = (*linesearch->ops->vinorm)(snes, G, W, &gnorm);CHKERRQ(ierr); 2709bd66eb0SPeter Brune } else { 27171b4ebd2SPeter Brune ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr); 2729bd66eb0SPeter Brune } 27371b4ebd2SPeter Brune if (PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 27471b4ebd2SPeter Brune if (.5*gnorm*gnorm < .5*fnorm*fnorm + lambda*alpha*initslope) { /* is reduction enough? */ 27571b4ebd2SPeter Brune if (monitor) { 276ea5d4fccSPeter Brune ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 277b000cd8dSPeter Brune if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) { 278c69d1a72SBarry Smith ierr = PetscViewerASCIIPrintf(monitor," Line search: Cubically determined step, current gnorm %14.12e lambda=%18.16e\n",(double)gnorm,(double)lambda);CHKERRQ(ierr); 27971b4ebd2SPeter Brune } else { 280c69d1a72SBarry Smith ierr = PetscViewerASCIIPrintf(monitor," Line search: Quadratically determined step, current gnorm %14.12e lambda=%18.16e\n",(double)gnorm,(double)lambda);CHKERRQ(ierr); 28171b4ebd2SPeter Brune } 282ea5d4fccSPeter Brune ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 28371b4ebd2SPeter Brune } 28471b4ebd2SPeter Brune break; 28571b4ebd2SPeter Brune } else { 28671b4ebd2SPeter Brune if (monitor) { 287ea5d4fccSPeter Brune ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 288b000cd8dSPeter Brune if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) { 289c69d1a72SBarry Smith ierr = PetscViewerASCIIPrintf(monitor," Line search: Cubic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n",(double)gnorm,(double)lambda);CHKERRQ(ierr); 29071b4ebd2SPeter Brune } else { 291c69d1a72SBarry Smith ierr = PetscViewerASCIIPrintf(monitor," Line search: Quadratic step no good, shrinking lambda, current gnorm %12.12e lambda=%18.16e\n",(double)gnorm,(double)lambda);CHKERRQ(ierr); 29271b4ebd2SPeter Brune } 293ea5d4fccSPeter Brune ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 29471b4ebd2SPeter Brune } 29571b4ebd2SPeter Brune } 29671b4ebd2SPeter Brune } 29771b4ebd2SPeter Brune } 29871b4ebd2SPeter Brune } 29971b4ebd2SPeter Brune 30071b4ebd2SPeter Brune /* postcheck */ 3017b1df9c1SPeter Brune ierr = SNESLineSearchPostCheck(linesearch,X,Y,W,&changed_y,&changed_w);CHKERRQ(ierr); 30271b4ebd2SPeter Brune if (changed_y) { 30371b4ebd2SPeter Brune ierr = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr); 3049bd66eb0SPeter Brune if (linesearch->ops->viproject) { 3059bd66eb0SPeter Brune ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr); 3069bd66eb0SPeter Brune } 30771b4ebd2SPeter Brune } 308c427506bSPeter Brune if (changed_y || changed_w) { /* recompute the function if the step has changed */ 309c427506bSPeter Brune ierr = SNESComputeFunction(snes,W,G);CHKERRQ(ierr); 31071b4ebd2SPeter Brune ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr); 31171b4ebd2SPeter Brune if (domainerror) { 312f1c6b773SPeter Brune ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr); 31371b4ebd2SPeter Brune PetscFunctionReturn(0); 31471b4ebd2SPeter Brune } 3159bd66eb0SPeter Brune if (linesearch->ops->vinorm) { 3169bd66eb0SPeter Brune gnorm = fnorm; 3179bd66eb0SPeter Brune ierr = (*linesearch->ops->vinorm)(snes, G, W, &gnorm);CHKERRQ(ierr); 3189bd66eb0SPeter Brune } else { 3199bd66eb0SPeter Brune ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr); 3209bd66eb0SPeter Brune } 3219bd66eb0SPeter Brune ierr = VecNorm(Y,NORM_2,&ynorm);CHKERRQ(ierr); 322c427506bSPeter Brune if (PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 3239bd66eb0SPeter Brune 324c427506bSPeter Brune } 32571b4ebd2SPeter Brune 32671b4ebd2SPeter Brune /* copy the solution over */ 32771b4ebd2SPeter Brune ierr = VecCopy(W, X);CHKERRQ(ierr); 328c427506bSPeter Brune ierr = VecCopy(G, F);CHKERRQ(ierr); 329c427506bSPeter Brune ierr = VecNorm(X, NORM_2, &xnorm);CHKERRQ(ierr); 330f1c6b773SPeter Brune ierr = SNESLineSearchSetLambda(linesearch, lambda);CHKERRQ(ierr); 331f1c6b773SPeter Brune ierr = SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm);CHKERRQ(ierr); 33271b4ebd2SPeter Brune PetscFunctionReturn(0); 33371b4ebd2SPeter Brune } 33471b4ebd2SPeter Brune 33571b4ebd2SPeter Brune #undef __FUNCT__ 3367f1410a3SPeter Brune #define __FUNCT__ "SNESLineSearchView_BT" 3377f1410a3SPeter Brune PetscErrorCode SNESLineSearchView_BT(SNESLineSearch linesearch, PetscViewer viewer) 3387f1410a3SPeter Brune { 3397f1410a3SPeter Brune PetscErrorCode ierr; 3407f1410a3SPeter Brune PetscBool iascii; 3417f1410a3SPeter Brune SNESLineSearch_BT *bt; 3427f1410a3SPeter Brune PetscFunctionBegin; 343251f4c67SDmitry Karpeev ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 3447f1410a3SPeter Brune bt = (SNESLineSearch_BT*)linesearch->data; 3457f1410a3SPeter Brune if (iascii) { 346b000cd8dSPeter Brune if (linesearch->order == SNES_LINESEARCH_ORDER_CUBIC) { 3477f1410a3SPeter Brune ierr = PetscViewerASCIIPrintf(viewer, " interpolation: cubic\n");CHKERRQ(ierr); 348b000cd8dSPeter Brune } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) { 3497f1410a3SPeter Brune ierr = PetscViewerASCIIPrintf(viewer, " interpolation: quadratic\n");CHKERRQ(ierr); 3507f1410a3SPeter Brune } 351007b6e36SPeter Brune ierr = PetscViewerASCIIPrintf(viewer, " alpha=%e\n", bt->alpha);CHKERRQ(ierr); 3527f1410a3SPeter Brune } 3537f1410a3SPeter Brune PetscFunctionReturn(0); 3547f1410a3SPeter Brune } 3557f1410a3SPeter Brune 3567f1410a3SPeter Brune 3577f1410a3SPeter Brune #undef __FUNCT__ 358f1c6b773SPeter Brune #define __FUNCT__ "SNESLineSearchDestroy_BT" 359f1c6b773SPeter Brune static PetscErrorCode SNESLineSearchDestroy_BT(SNESLineSearch linesearch) 360731c067cSJed Brown { 361731c067cSJed Brown PetscErrorCode ierr; 36271b4ebd2SPeter Brune 36371b4ebd2SPeter Brune PetscFunctionBegin; 364731c067cSJed Brown ierr = PetscFree(linesearch->data);CHKERRQ(ierr); 36571b4ebd2SPeter Brune PetscFunctionReturn(0); 36671b4ebd2SPeter Brune } 36771b4ebd2SPeter Brune 36871b4ebd2SPeter Brune 36971b4ebd2SPeter Brune #undef __FUNCT__ 370f1c6b773SPeter Brune #define __FUNCT__ "SNESLineSearchSetFromOptions_BT" 371f1c6b773SPeter Brune static PetscErrorCode SNESLineSearchSetFromOptions_BT(SNESLineSearch linesearch) 37271b4ebd2SPeter Brune { 37371b4ebd2SPeter Brune 37471b4ebd2SPeter Brune PetscErrorCode ierr; 375f1c6b773SPeter Brune SNESLineSearch_BT *bt; 37671b4ebd2SPeter Brune PetscFunctionBegin; 37771b4ebd2SPeter Brune 378f1c6b773SPeter Brune bt = (SNESLineSearch_BT*)linesearch->data; 37971b4ebd2SPeter Brune 380f1c6b773SPeter Brune ierr = PetscOptionsHead("SNESLineSearch BT options");CHKERRQ(ierr); 3817a35526eSPeter Brune ierr = PetscOptionsReal("-snes_linesearch_alpha", "Descent tolerance", "SNESLineSearchBT", bt->alpha, &bt->alpha, PETSC_NULL);CHKERRQ(ierr); 38271b4ebd2SPeter Brune 38371b4ebd2SPeter Brune ierr = PetscOptionsTail();CHKERRQ(ierr); 38471b4ebd2SPeter Brune PetscFunctionReturn(0); 38571b4ebd2SPeter Brune } 38671b4ebd2SPeter Brune 38771b4ebd2SPeter Brune 38871b4ebd2SPeter Brune #undef __FUNCT__ 389f1c6b773SPeter Brune #define __FUNCT__ "SNESLineSearchCreate_BT" 390954494b2SJed Brown /*MC 391db609ea7SPeter Brune SNESLINESEARCHBT - Backtracking line search. 392954494b2SJed Brown 393db609ea7SPeter Brune This line search finds the minimum of a polynomial fitting of the L2 norm of the 394db609ea7SPeter Brune function. If this fit does not satisfy the conditions for progress, the interval shrinks 395db609ea7SPeter Brune and the fit is reattempted at most max_it times or until lambda is below minlambda. 396954494b2SJed Brown 397cd7522eaSPeter Brune Options Database Keys: 398cd7522eaSPeter Brune + -snes_linesearch_alpha<1e-4> - slope descent parameter 399db609ea7SPeter Brune . -snes_linesearch_damping<1.0> - initial step length 400db609ea7SPeter Brune . -snes_linesearch_max_it<40> - maximum number of shrinking step 401db609ea7SPeter Brune . -snes_linesearch_minlambda<1e-12> - minimum step length allowed 402cd7522eaSPeter Brune - -snes_linesearch_order<cubic,quadratic> - order of the approximation 403cd7522eaSPeter Brune 404954494b2SJed Brown Level: advanced 405954494b2SJed Brown 406db609ea7SPeter Brune Notes: 407db609ea7SPeter Brune This line search is taken from "Numerical Methods for Unconstrained 408db609ea7SPeter Brune Optimization and Nonlinear Equations" by Dennis and Schnabel, page 325. 409db609ea7SPeter Brune 410f1c6b773SPeter Brune .keywords: SNES, SNESLineSearch, damping 411954494b2SJed Brown 412f1c6b773SPeter Brune .seealso: SNESLineSearchCreate(), SNESLineSearchSetType() 413954494b2SJed Brown M*/ 414f1c6b773SPeter Brune PETSC_EXTERN_C PetscErrorCode SNESLineSearchCreate_BT(SNESLineSearch linesearch) 41571b4ebd2SPeter Brune { 41671b4ebd2SPeter Brune 417f1c6b773SPeter Brune SNESLineSearch_BT *bt; 41871b4ebd2SPeter Brune PetscErrorCode ierr; 41971b4ebd2SPeter Brune 42071b4ebd2SPeter Brune PetscFunctionBegin; 421f1c6b773SPeter Brune linesearch->ops->apply = SNESLineSearchApply_BT; 422f1c6b773SPeter Brune linesearch->ops->destroy = SNESLineSearchDestroy_BT; 423f1c6b773SPeter Brune linesearch->ops->setfromoptions = SNESLineSearchSetFromOptions_BT; 42471b4ebd2SPeter Brune linesearch->ops->reset = PETSC_NULL; 4257f1410a3SPeter Brune linesearch->ops->view = SNESLineSearchView_BT; 42671b4ebd2SPeter Brune linesearch->ops->setup = PETSC_NULL; 42771b4ebd2SPeter Brune 428f1c6b773SPeter Brune ierr = PetscNewLog(linesearch, SNESLineSearch_BT, &bt);CHKERRQ(ierr); 42971b4ebd2SPeter Brune linesearch->data = (void *)bt; 43071b4ebd2SPeter Brune linesearch->max_its = 40; 431b000cd8dSPeter Brune linesearch->order = SNES_LINESEARCH_ORDER_CUBIC; 43271b4ebd2SPeter Brune bt->alpha = 1e-4; 43371b4ebd2SPeter Brune 43471b4ebd2SPeter Brune PetscFunctionReturn(0); 43571b4ebd2SPeter Brune } 436