1 #include <petsc-private/linesearchimpl.h> 2 #include <petscsnes.h> 3 4 #undef __FUNCT__ 5 #define __FUNCT__ "SNESLineSearchApply_CP" 6 static PetscErrorCode SNESLineSearchApply_CP(SNESLineSearch linesearch) 7 { 8 PetscBool changed_y, changed_w, domainerror; 9 PetscErrorCode ierr; 10 Vec X, Y, F, W; 11 SNES snes; 12 PetscReal xnorm, ynorm, gnorm, steptol, atol, rtol, ltol, maxstep; 13 14 PetscReal lambda, lambda_old, lambda_update, delLambda; 15 PetscScalar fty, fty_init, fty_old, fty_mid1, fty_mid2, s; 16 PetscInt i, max_its; 17 18 PetscViewer monitor; 19 20 PetscFunctionBegin; 21 ierr = SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, NULL);CHKERRQ(ierr); 22 ierr = SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm);CHKERRQ(ierr); 23 ierr = SNESLineSearchGetSNES(linesearch, &snes);CHKERRQ(ierr); 24 ierr = SNESLineSearchGetLambda(linesearch, &lambda);CHKERRQ(ierr); 25 ierr = SNESLineSearchGetTolerances(linesearch, &steptol, &maxstep, &rtol, &atol, <ol, &max_its);CHKERRQ(ierr); 26 ierr = SNESLineSearchSetSuccess(linesearch, PETSC_TRUE);CHKERRQ(ierr); 27 ierr = SNESLineSearchGetMonitor(linesearch, &monitor);CHKERRQ(ierr); 28 29 /* precheck */ 30 ierr = SNESLineSearchPreCheck(linesearch,X,Y,&changed_y);CHKERRQ(ierr); 31 lambda_old = 0.0; 32 33 34 for (i = 0; i < max_its; i++) { 35 if (i == 0) { 36 ierr = VecDotBegin(F,Y,&fty_old);CHKERRQ(ierr); 37 fty_init = fty_old; 38 } 39 /* compute the norm at lambda */ 40 ierr = VecCopy(X, W);CHKERRQ(ierr); 41 ierr = VecAXPY(W, -lambda, Y);CHKERRQ(ierr); 42 if (linesearch->ops->viproject) { 43 ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr); 44 } 45 ierr = (*linesearch->ops->snesfunc)(snes,W,F);CHKERRQ(ierr); 46 ierr = VecDotBegin(F,Y,&fty);CHKERRQ(ierr); 47 if (i == 0) { 48 ierr = VecDotEnd(F,Y,&fty_old);CHKERRQ(ierr); 49 } 50 ierr = VecDotEnd(F,Y,&fty);CHKERRQ(ierr); 51 52 delLambda = lambda - lambda_old; 53 54 /* check for convergence */ 55 if (PetscAbsReal(delLambda) < steptol*lambda) break; 56 if (PetscAbsScalar(fty) / PetscAbsScalar(fty_init) < rtol) break; 57 if (PetscAbsScalar(fty) < atol && i > 0) break; 58 if (monitor) { 59 ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 60 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); 61 ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 62 } 63 64 /* compute the search direction */ 65 if (linesearch->order == SNES_LINESEARCH_ORDER_LINEAR) { 66 s = (fty - fty_old) / delLambda; 67 } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) { 68 ierr = VecCopy(X, W);CHKERRQ(ierr); 69 ierr = VecAXPY(W, -0.5*(lambda + lambda_old), Y);CHKERRQ(ierr); 70 if (linesearch->ops->viproject) { 71 ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr); 72 } 73 ierr = (*linesearch->ops->snesfunc)(snes,W,F);CHKERRQ(ierr); 74 ierr = VecDot(F, Y, &fty_mid1);CHKERRQ(ierr); 75 s = (3.*fty - 4.*fty_mid1 + fty_old) / delLambda; 76 } else { 77 ierr = VecCopy(X, W);CHKERRQ(ierr); 78 ierr = VecAXPY(W, -0.5*(lambda + lambda_old), Y);CHKERRQ(ierr); 79 if (linesearch->ops->viproject) { 80 ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr); 81 } 82 ierr = (*linesearch->ops->snesfunc)(snes,W,F);CHKERRQ(ierr); 83 ierr = VecDotBegin(F, Y, &fty_mid1);CHKERRQ(ierr); 84 ierr = VecCopy(X, W);CHKERRQ(ierr); 85 ierr = VecAXPY(W, -(lambda + 0.5*(lambda - lambda_old)), Y);CHKERRQ(ierr); 86 if (linesearch->ops->viproject) { 87 ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr); 88 } 89 ierr = (*linesearch->ops->snesfunc)(snes, W, F);CHKERRQ(ierr); 90 ierr = VecDotBegin(F, Y, &fty_mid2);CHKERRQ(ierr); 91 ierr = VecDotEnd(F, Y, &fty_mid1);CHKERRQ(ierr); 92 ierr = VecDotEnd(F, Y, &fty_mid2);CHKERRQ(ierr); 93 s = (2.*fty_mid2 + 3.*fty - 6.*fty_mid1 + fty_old) / (3.*delLambda); 94 } 95 /* if the solve is going in the wrong direction, fix it */ 96 if (PetscRealPart(s) > 0.) s = -s; 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 (PetscIsInfOrNanScalar(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 ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr); 128 if (domainerror) { 129 ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr); 130 PetscFunctionReturn(0); 131 } 132 133 ierr = SNESLineSearchComputeNorms(linesearch);CHKERRQ(ierr); 134 ierr = SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm);CHKERRQ(ierr); 135 136 ierr = SNESLineSearchSetLambda(linesearch, lambda);CHKERRQ(ierr); 137 138 if (monitor) { 139 ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 140 ierr = PetscViewerASCIIPrintf(monitor," Line search terminated: lambda = %g, fnorms = %g\n", (double)lambda, (double)gnorm);CHKERRQ(ierr); 141 ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 142 } 143 if (lambda <= steptol) { 144 ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr); 145 } 146 PetscFunctionReturn(0); 147 } 148 149 #undef __FUNCT__ 150 #define __FUNCT__ "SNESLineSearchCreate_CP" 151 /*MC 152 SNESLINESEARCHCP - Critical point line search. This line search assumes that there exists some 153 artificial G(x) for which the SNESFunction F(x) = grad G(x). Therefore, this line search seeks 154 to find roots of dot(F, Y) via a secant method. 155 156 Options Database Keys: 157 + -snes_linesearch_minlambda - the minimum acceptable lambda 158 . -snes_linesearch_damping - initial trial step length 159 - -snes_linesearch_max_it - the maximum number of secant steps performed. 160 161 Notes: 162 This method is the preferred line search for SNESQN and SNESNCG. 163 164 Level: advanced 165 166 .keywords: SNES, SNESLineSearch, damping 167 168 .seealso: SNESLineSearchCreate(), SNESLineSearchSetType() 169 M*/ 170 PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_CP(SNESLineSearch linesearch) 171 { 172 PetscFunctionBegin; 173 linesearch->ops->apply = SNESLineSearchApply_CP; 174 linesearch->ops->destroy = NULL; 175 linesearch->ops->setfromoptions = NULL; 176 linesearch->ops->reset = NULL; 177 linesearch->ops->view = NULL; 178 linesearch->ops->setup = NULL; 179 linesearch->order = SNES_LINESEARCH_ORDER_LINEAR; 180 181 linesearch->max_its = 1; 182 PetscFunctionReturn(0); 183 } 184