1 #include <petsc/private/linesearchimpl.h> 2 #include <petscsnes.h> 3 4 static PetscErrorCode SNESLineSearchApply_CP(SNESLineSearch linesearch) 5 { 6 PetscBool changed_y, changed_w; 7 Vec X, Y, F, W; 8 SNES snes; 9 PetscReal xnorm, ynorm, gnorm, steptol, atol, rtol, ltol, maxstep; 10 PetscReal lambda, lambda_old, lambda_update, delLambda; 11 PetscScalar fty, fty_init, fty_old, fty_mid1, fty_mid2, s; 12 PetscInt i, max_its; 13 PetscViewer monitor; 14 15 PetscFunctionBegin; 16 PetscCall(SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, NULL)); 17 PetscCall(SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm)); 18 PetscCall(SNESLineSearchGetSNES(linesearch, &snes)); 19 PetscCall(SNESLineSearchGetLambda(linesearch, &lambda)); 20 PetscCall(SNESLineSearchGetTolerances(linesearch, &steptol, &maxstep, &rtol, &atol, <ol, &max_its)); 21 PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED)); 22 PetscCall(SNESLineSearchGetDefaultMonitor(linesearch, &monitor)); 23 24 /* precheck */ 25 PetscCall(SNESLineSearchPreCheck(linesearch,X,Y,&changed_y)); 26 lambda_old = 0.0; 27 28 PetscCall(VecDot(F,Y,&fty_old)); 29 if (PetscAbsScalar(fty_old) < atol * ynorm) { 30 if (monitor) { 31 PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel)); 32 PetscCall(PetscViewerASCIIPrintf(monitor," Line search terminated at initial point because dot(F,Y) = %g < atol*||y|| = %g\n",(double)PetscAbsScalar(fty_old), (double)(atol*ynorm))); 33 PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel)); 34 } 35 PetscCall(SNESSetConvergedReason(linesearch->snes,SNES_CONVERGED_FNORM_ABS)); 36 PetscFunctionReturn(0); 37 } 38 39 fty_init = fty_old; 40 41 for (i = 0; i < max_its; i++) { 42 /* compute the norm at lambda */ 43 PetscCall(VecWAXPY(W, -lambda, Y, X)); 44 if (linesearch->ops->viproject) { 45 PetscCall((*linesearch->ops->viproject)(snes, W)); 46 } 47 PetscCall((*linesearch->ops->snesfunc)(snes,W,F)); 48 PetscCall(VecDot(F,Y,&fty)); 49 50 delLambda = lambda - lambda_old; 51 52 /* check for convergence */ 53 if (PetscAbsReal(delLambda) < steptol*lambda) break; 54 if (PetscAbsScalar(fty) / PetscAbsScalar(fty_init) < rtol) break; 55 if (PetscAbsScalar(fty) < atol * ynorm && i > 0) break; 56 if (monitor) { 57 PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel)); 58 PetscCall(PetscViewerASCIIPrintf(monitor," Line search: lambdas = [%g, %g], ftys = [%g, %g]\n",(double)lambda, (double)lambda_old, (double)PetscRealPart(fty), (double)PetscRealPart(fty_old))); 59 PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel)); 60 } 61 62 /* compute the search direction */ 63 if (linesearch->order == SNES_LINESEARCH_ORDER_LINEAR) { 64 s = (fty - fty_old) / delLambda; 65 } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) { 66 PetscCall(VecWAXPY(W, -0.5*(lambda + lambda_old), Y, X)); 67 if (linesearch->ops->viproject) { 68 PetscCall((*linesearch->ops->viproject)(snes, W)); 69 } 70 PetscCall((*linesearch->ops->snesfunc)(snes,W,F)); 71 PetscCall(VecDot(F, Y, &fty_mid1)); 72 s = (3.*fty - 4.*fty_mid1 + fty_old) / delLambda; 73 } else { 74 PetscCall(VecWAXPY(W, -0.5*(lambda + lambda_old), Y, X)); 75 if (linesearch->ops->viproject) { 76 PetscCall((*linesearch->ops->viproject)(snes, W)); 77 } 78 PetscCall((*linesearch->ops->snesfunc)(snes,W,F)); 79 PetscCall(VecDot(F, Y, &fty_mid1)); 80 PetscCall(VecWAXPY(W, -(lambda + 0.5*(lambda - lambda_old)), Y, X)); 81 if (linesearch->ops->viproject) { 82 PetscCall((*linesearch->ops->viproject)(snes, W)); 83 } 84 PetscCall((*linesearch->ops->snesfunc)(snes, W, F)); 85 PetscCall(VecDot(F, Y, &fty_mid2)); 86 s = (2.*fty_mid2 + 3.*fty - 6.*fty_mid1 + fty_old) / (3.*delLambda); 87 } 88 /* if the solve is going in the wrong direction, fix it */ 89 if (PetscRealPart(s) > 0.) s = -s; 90 if (s == 0.0) break; 91 lambda_update = lambda - PetscRealPart(fty / s); 92 93 /* switch directions if we stepped out of bounds */ 94 if (lambda_update < steptol) lambda_update = lambda + PetscRealPart(fty / s); 95 96 if (PetscIsInfOrNanReal(lambda_update)) break; 97 if (lambda_update > maxstep) break; 98 99 /* compute the new state of the line search */ 100 lambda_old = lambda; 101 lambda = lambda_update; 102 fty_old = fty; 103 } 104 /* construct the solution */ 105 PetscCall(VecWAXPY(W, -lambda, Y, X)); 106 if (linesearch->ops->viproject) { 107 PetscCall((*linesearch->ops->viproject)(snes, W)); 108 } 109 /* postcheck */ 110 PetscCall(SNESLineSearchPostCheck(linesearch,X,Y,W,&changed_y,&changed_w)); 111 if (changed_y) { 112 PetscCall(VecAXPY(X, -lambda, Y)); 113 if (linesearch->ops->viproject) { 114 PetscCall((*linesearch->ops->viproject)(snes, X)); 115 } 116 } else { 117 PetscCall(VecCopy(W, X)); 118 } 119 PetscCall((*linesearch->ops->snesfunc)(snes,X,F)); 120 121 PetscCall(SNESLineSearchComputeNorms(linesearch)); 122 PetscCall(SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm)); 123 124 PetscCall(SNESLineSearchSetLambda(linesearch, lambda)); 125 126 if (monitor) { 127 PetscCall(PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel)); 128 PetscCall(PetscViewerASCIIPrintf(monitor," Line search terminated: lambda = %g, fnorms = %g\n", (double)lambda, (double)gnorm)); 129 PetscCall(PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel)); 130 } 131 if (lambda <= steptol) { 132 PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT)); 133 } 134 PetscFunctionReturn(0); 135 } 136 137 /*MC 138 SNESLINESEARCHCP - Critical point line search. This line search assumes that there exists some 139 artificial G(x) for which the SNESFunction F(x) = grad G(x). Therefore, this line search seeks 140 to find roots of dot(F, Y) via a secant method. 141 142 Options Database Keys: 143 + -snes_linesearch_minlambda <minlambda> - the minimum acceptable lambda 144 . -snes_linesearch_maxstep <length> - the algorithm insures that a step length is never longer than this value 145 . -snes_linesearch_damping <damping> - initial trial step length is scaled by this factor, default is 1.0 146 - -snes_linesearch_max_it <max_it> - the maximum number of secant steps performed. 147 148 Notes: 149 This method does NOT use the objective function if it is provided with SNESSetObjective(). 150 151 This method is the preferred line search for SNESQN and SNESNCG. 152 153 Level: advanced 154 155 .seealso: SNESLineSearchCreate(), SNESLineSearchSetType() 156 M*/ 157 PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_CP(SNESLineSearch linesearch) 158 { 159 PetscFunctionBegin; 160 linesearch->ops->apply = SNESLineSearchApply_CP; 161 linesearch->ops->destroy = NULL; 162 linesearch->ops->setfromoptions = NULL; 163 linesearch->ops->reset = NULL; 164 linesearch->ops->view = NULL; 165 linesearch->ops->setup = NULL; 166 linesearch->order = SNES_LINESEARCH_ORDER_LINEAR; 167 168 linesearch->max_its = 1; 169 PetscFunctionReturn(0); 170 } 171