1 #include <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_old; 16 PetscInt i, max_its; 17 18 PetscViewer monitor; 19 20 PetscFunctionBegin; 21 22 ierr = SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, PETSC_NULL);CHKERRQ(ierr); 23 ierr = SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm);CHKERRQ(ierr); 24 ierr = SNESLineSearchGetSNES(linesearch, &snes);CHKERRQ(ierr); 25 ierr = SNESLineSearchGetLambda(linesearch, &lambda);CHKERRQ(ierr); 26 ierr = SNESLineSearchGetTolerances(linesearch, &steptol, &maxstep, &rtol, &atol, <ol, &max_its);CHKERRQ(ierr); 27 ierr = SNESLineSearchSetSuccess(linesearch, PETSC_TRUE);CHKERRQ(ierr); 28 ierr = SNESLineSearchGetMonitor(linesearch, &monitor);CHKERRQ(ierr); 29 30 /* precheck */ 31 ierr = SNESLineSearchPreCheck(linesearch, &changed_y);CHKERRQ(ierr); 32 lambda_old = 0.0; 33 ierr = VecDot(F, Y, &fty_old);CHKERRQ(ierr); 34 for (i = 0; i < max_its; i++) { 35 36 /* compute the norm at lambda */ 37 ierr = VecCopy(X, W);CHKERRQ(ierr); 38 ierr = VecAXPY(W, -lambda, Y);CHKERRQ(ierr); 39 if (linesearch->ops->viproject) { 40 ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr); 41 } 42 ierr = SNESComputeFunction(snes, W, F);CHKERRQ(ierr); 43 44 ierr = VecDot(F, Y, &fty);CHKERRQ(ierr); 45 46 delLambda = lambda - lambda_old; 47 if (PetscAbsReal(delLambda) < steptol) break; 48 if (monitor) { 49 ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 50 ierr = PetscViewerASCIIPrintf(monitor," Line search: lambdas = [%g, %g], ftys = [%g, %g]\n", 51 lambda, lambda_old, PetscRealPart(fty), PetscRealPart(fty_old));CHKERRQ(ierr); 52 ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 53 } 54 55 /* compute the search direction */ 56 lambda_update = PetscRealPart((fty*lambda_old - fty_old*lambda) / (fty - fty_old)); 57 if (PetscIsInfOrNanScalar(lambda_update)) break; 58 if (lambda_update > maxstep) { 59 break; 60 } 61 62 /* compute the new state of the line search */ 63 lambda_old = lambda; 64 lambda = lambda_update; 65 fty_old = fty; 66 } 67 /* construct the solution */ 68 ierr = VecCopy(X, W);CHKERRQ(ierr); 69 ierr = VecAXPY(W, -lambda, Y);CHKERRQ(ierr); 70 if (linesearch->ops->viproject) { 71 ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr); 72 } 73 /* postcheck */ 74 ierr = SNESLineSearchPostCheck(linesearch, &changed_y, &changed_w);CHKERRQ(ierr); 75 if (changed_y) { 76 ierr = VecAXPY(X, -lambda, Y);CHKERRQ(ierr); 77 if (linesearch->ops->viproject) { 78 ierr = (*linesearch->ops->viproject)(snes, X);CHKERRQ(ierr); 79 } 80 } else { 81 ierr = VecCopy(W, X);CHKERRQ(ierr); 82 } 83 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); 84 ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr); 85 if (domainerror) { 86 ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr); 87 PetscFunctionReturn(0); 88 } 89 90 ierr = SNESLineSearchComputeNorms(linesearch);CHKERRQ(ierr); 91 ierr = SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm);CHKERRQ(ierr); 92 93 ierr = SNESLineSearchSetLambda(linesearch, lambda);CHKERRQ(ierr); 94 95 if (monitor) { 96 ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 97 ierr = PetscViewerASCIIPrintf(monitor," Line search terminated: lambda = %g, fnorms = %g\n", lambda, gnorm);CHKERRQ(ierr); 98 ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr); 99 } 100 if (lambda <= steptol) { 101 ierr = SNESLineSearchSetSuccess(linesearch, PETSC_FALSE);CHKERRQ(ierr); 102 } 103 PetscFunctionReturn(0); 104 } 105 106 #undef __FUNCT__ 107 #define __FUNCT__ "SNESLineSearchCreate_CP" 108 /*MC 109 SNES_LINESEARCH_CP - Critical point line search 110 111 Level: advanced 112 113 .keywords: SNES, SNESLineSearch, damping 114 115 .seealso: SNESLineSearchCreate(), SNESLineSearchSetType() 116 M*/ 117 PETSC_EXTERN_C PetscErrorCode SNESLineSearchCreate_CP(SNESLineSearch linesearch) 118 { 119 PetscFunctionBegin; 120 linesearch->ops->apply = SNESLineSearchApply_CP; 121 linesearch->ops->destroy = PETSC_NULL; 122 linesearch->ops->setfromoptions = PETSC_NULL; 123 linesearch->ops->reset = PETSC_NULL; 124 linesearch->ops->view = PETSC_NULL; 125 linesearch->ops->setup = PETSC_NULL; 126 PetscFunctionReturn(0); 127 } 128