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