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 lambda_update = lambda - PetscRealPart(fty / s); 97 98 /* switch directions if we stepped out of bounds */ 99 if (lambda_update < steptol) lambda_update = lambda + PetscRealPart(fty / s); 100 101 if (PetscIsInfOrNanReal(lambda_update)) break; 102 if (lambda_update > maxstep) break; 103 104 /* compute the new state of the line search */ 105 lambda_old = lambda; 106 lambda = lambda_update; 107 fty_old = fty; 108 } 109 /* construct the solution */ 110 ierr = VecCopy(X, W);CHKERRQ(ierr); 111 ierr = VecAXPY(W, -lambda, Y);CHKERRQ(ierr); 112 if (linesearch->ops->viproject) { 113 ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr); 114 } 115 /* postcheck */ 116 ierr = SNESLineSearchPostCheck(linesearch,X,Y,W,&changed_y,&changed_w);CHKERRQ(ierr); 117 if (changed_y) { 118 ierr = VecAXPY(X, -lambda, Y);CHKERRQ(ierr); 119 if (linesearch->ops->viproject) { 120 ierr = (*linesearch->ops->viproject)(snes, X);CHKERRQ(ierr); 121 } 122 } else { 123 ierr = VecCopy(W, X);CHKERRQ(ierr); 124 } 125 ierr = (*linesearch->ops->snesfunc)(snes,X,F);CHKERRQ(ierr); 126 127 ierr = SNESLineSearchComputeNorms(linesearch);CHKERRQ(ierr); 128 ierr = SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm);CHKERRQ(ierr); 129 130 ierr = SNESLineSearchSetLambda(linesearch, lambda);CHKERRQ(ierr); 131 132 if (monitor) { 133 ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 134 ierr = PetscViewerASCIIPrintf(monitor," Line search terminated: lambda = %g, fnorms = %g\n", (double)lambda, (double)gnorm);CHKERRQ(ierr); 135 ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 136 } 137 if (lambda <= steptol) { 138 ierr = SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT);CHKERRQ(ierr); 139 } 140 PetscFunctionReturn(0); 141 } 142 143 /*MC 144 SNESLINESEARCHCP - Critical point line search. This line search assumes that there exists some 145 artificial G(x) for which the SNESFunction F(x) = grad G(x). Therefore, this line search seeks 146 to find roots of dot(F, Y) via a secant method. 147 148 Options Database Keys: 149 + -snes_linesearch_minlambda <minlambda> - the minimum acceptable lambda 150 . -snes_linesearch_maxstep <length> - the algorithm insures that a step length is never longer than this value 151 . -snes_linesearch_damping <damping> - initial trial step length is scaled by this factor, default is 1.0 152 - -snes_linesearch_max_it <max_it> - the maximum number of secant steps performed. 153 154 Notes: 155 This method does NOT use the objective function if it is provided with SNESSetObjective(). 156 157 This method is the preferred line search for SNESQN and SNESNCG. 158 159 Level: advanced 160 161 .keywords: SNES, SNESLineSearch, damping 162 163 .seealso: SNESLineSearchCreate(), SNESLineSearchSetType() 164 M*/ 165 PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_CP(SNESLineSearch linesearch) 166 { 167 PetscFunctionBegin; 168 linesearch->ops->apply = SNESLineSearchApply_CP; 169 linesearch->ops->destroy = NULL; 170 linesearch->ops->setfromoptions = NULL; 171 linesearch->ops->reset = NULL; 172 linesearch->ops->view = NULL; 173 linesearch->ops->setup = NULL; 174 linesearch->order = SNES_LINESEARCH_ORDER_LINEAR; 175 176 linesearch->max_its = 1; 177 PetscFunctionReturn(0); 178 } 179