1 #include <petsc/private/taolinesearchimpl.h> 2 #include <../src/tao/linesearch/impls/gpcglinesearch/gpcglinesearch.h> 3 4 /* ---------------------------------------------------------- */ 5 6 static PetscErrorCode TaoLineSearchDestroy_GPCG(TaoLineSearch ls) 7 { 8 TaoLineSearch_GPCG *ctx = (TaoLineSearch_GPCG *)ls->data; 9 10 PetscFunctionBegin; 11 PetscCall(VecDestroy(&ctx->W1)); 12 PetscCall(VecDestroy(&ctx->W2)); 13 PetscCall(VecDestroy(&ctx->Gold)); 14 PetscCall(VecDestroy(&ctx->x)); 15 PetscCall(PetscFree(ls->data)); 16 PetscFunctionReturn(0); 17 } 18 19 /*------------------------------------------------------------*/ 20 static PetscErrorCode TaoLineSearchView_GPCG(TaoLineSearch ls, PetscViewer viewer) 21 { 22 PetscBool isascii; 23 24 PetscFunctionBegin; 25 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii)); 26 if (isascii) { 27 PetscCall(PetscViewerASCIIPrintf(viewer," GPCG Line search")); 28 } 29 PetscFunctionReturn(0); 30 } 31 32 /*------------------------------------------------------------*/ 33 static PetscErrorCode TaoLineSearchApply_GPCG(TaoLineSearch ls, Vec x, PetscReal *f, Vec g, Vec s) 34 { 35 TaoLineSearch_GPCG *neP = (TaoLineSearch_GPCG *)ls->data; 36 PetscInt i; 37 PetscBool g_computed=PETSC_FALSE; /* to prevent extra gradient computation */ 38 PetscReal d1,finit,actred,prered,rho, gdx; 39 40 PetscFunctionBegin; 41 /* ls->stepmin - lower bound for step */ 42 /* ls->stepmax - upper bound for step */ 43 /* ls->rtol - relative tolerance for an acceptable step */ 44 /* ls->ftol - tolerance for sufficient decrease condition */ 45 /* ls->gtol - tolerance for curvature condition */ 46 /* ls->nfeval - number of function evaluations */ 47 /* ls->nfeval - number of function/gradient evaluations */ 48 /* ls->max_funcs - maximum number of function evaluations */ 49 50 PetscCall(TaoLineSearchMonitor(ls, 0, *f, 0.0)); 51 52 ls->reason = TAOLINESEARCH_CONTINUE_ITERATING; 53 ls->step = ls->initstep; 54 if (!neP->W2) { 55 PetscCall(VecDuplicate(x,&neP->W2)); 56 PetscCall(VecDuplicate(x,&neP->W1)); 57 PetscCall(VecDuplicate(x,&neP->Gold)); 58 neP->x = x; 59 PetscCall(PetscObjectReference((PetscObject)neP->x)); 60 } else if (x != neP->x) { 61 PetscCall(VecDestroy(&neP->x)); 62 PetscCall(VecDestroy(&neP->W1)); 63 PetscCall(VecDestroy(&neP->W2)); 64 PetscCall(VecDestroy(&neP->Gold)); 65 PetscCall(VecDuplicate(x,&neP->W1)); 66 PetscCall(VecDuplicate(x,&neP->W2)); 67 PetscCall(VecDuplicate(x,&neP->Gold)); 68 PetscCall(PetscObjectDereference((PetscObject)neP->x)); 69 neP->x = x; 70 PetscCall(PetscObjectReference((PetscObject)neP->x)); 71 } 72 73 PetscCall(VecDot(g,s,&gdx)); 74 if (gdx > 0) { 75 PetscCall(PetscInfo(ls,"Line search error: search direction is not descent direction. dot(g,s) = %g\n",(double)gdx)); 76 ls->reason = TAOLINESEARCH_FAILED_ASCENT; 77 PetscFunctionReturn(0); 78 } 79 PetscCall(VecCopy(x,neP->W2)); 80 PetscCall(VecCopy(g,neP->Gold)); 81 if (ls->bounded) { 82 /* Compute the smallest steplength that will make one nonbinding variable equal the bound */ 83 PetscCall(VecStepBoundInfo(x,s,ls->lower,ls->upper,&rho,&actred,&d1)); 84 ls->step = PetscMin(ls->step,d1); 85 } 86 rho=0; actred=0; 87 88 if (ls->step < 0) { 89 PetscCall(PetscInfo(ls,"Line search error: initial step parameter %g< 0\n",(double)ls->step)); 90 ls->reason = TAOLINESEARCH_HALTED_OTHER; 91 PetscFunctionReturn(0); 92 } 93 94 /* Initialization */ 95 finit = *f; 96 for (i=0; i< ls->max_funcs; i++) { 97 /* Force the step to be within the bounds */ 98 ls->step = PetscMax(ls->step,ls->stepmin); 99 ls->step = PetscMin(ls->step,ls->stepmax); 100 101 PetscCall(VecCopy(x,neP->W2)); 102 PetscCall(VecAXPY(neP->W2,ls->step,s)); 103 if (ls->bounded) { 104 /* Make sure new vector is numerically within bounds */ 105 PetscCall(VecMedian(neP->W2,ls->lower,ls->upper,neP->W2)); 106 } 107 108 /* Gradient is not needed here. Unless there is a separate 109 gradient routine, compute it here anyway to prevent recomputing at 110 the end of the line search */ 111 if (ls->hasobjective) { 112 PetscCall(TaoLineSearchComputeObjective(ls,neP->W2,f)); 113 g_computed=PETSC_FALSE; 114 } else if (ls->usegts) { 115 PetscCall(TaoLineSearchComputeObjectiveAndGTS(ls,neP->W2,f,&gdx)); 116 g_computed=PETSC_FALSE; 117 } else { 118 PetscCall(TaoLineSearchComputeObjectiveAndGradient(ls,neP->W2,f,g)); 119 g_computed=PETSC_TRUE; 120 } 121 122 PetscCall(TaoLineSearchMonitor(ls, i+1, *f, ls->step)); 123 124 if (0 == i) { 125 ls->f_fullstep = *f; 126 } 127 128 actred = *f - finit; 129 PetscCall(VecCopy(neP->W2,neP->W1)); 130 PetscCall(VecAXPY(neP->W1,-1.0,x)); /* W1 = W2 - X */ 131 PetscCall(VecDot(neP->W1,neP->Gold,&prered)); 132 133 if (PetscAbsReal(prered)<1.0e-100) prered=1.0e-12; 134 rho = actred/prered; 135 136 /* 137 If sufficient progress has been obtained, accept the 138 point. Otherwise, backtrack. 139 */ 140 141 if (actred > 0) { 142 PetscCall(PetscInfo(ls,"Step resulted in ascent, rejecting.\n")); 143 ls->step = (ls->step)/2; 144 } else if (rho > ls->ftol) { 145 break; 146 } else{ 147 ls->step = (ls->step)/2; 148 } 149 150 /* Convergence testing */ 151 152 if (ls->step <= ls->stepmin || ls->step >= ls->stepmax) { 153 ls->reason = TAOLINESEARCH_HALTED_OTHER; 154 PetscCall(PetscInfo(ls,"Rounding errors may prevent further progress. May not be a step satisfying\n")); 155 PetscCall(PetscInfo(ls,"sufficient decrease and curvature conditions. Tolerances may be too small.\n")); 156 break; 157 } 158 if (ls->step == ls->stepmax) { 159 PetscCall(PetscInfo(ls,"Step is at the upper bound, stepmax (%g)\n",(double)ls->stepmax)); 160 ls->reason = TAOLINESEARCH_HALTED_UPPERBOUND; 161 break; 162 } 163 if (ls->step == ls->stepmin) { 164 PetscCall(PetscInfo(ls,"Step is at the lower bound, stepmin (%g)\n",(double)ls->stepmin)); 165 ls->reason = TAOLINESEARCH_HALTED_LOWERBOUND; 166 break; 167 } 168 if ((ls->nfeval+ls->nfgeval) >= ls->max_funcs) { 169 PetscCall(PetscInfo(ls,"Number of line search function evals (%D) > maximum (%D)\n",ls->nfeval+ls->nfgeval,ls->max_funcs)); 170 ls->reason = TAOLINESEARCH_HALTED_MAXFCN; 171 break; 172 } 173 if ((neP->bracket) && (ls->stepmax - ls->stepmin <= ls->rtol*ls->stepmax)) { 174 PetscCall(PetscInfo(ls,"Relative width of interval of uncertainty is at most rtol (%g)\n",(double)ls->rtol)); 175 ls->reason = TAOLINESEARCH_HALTED_RTOL; 176 break; 177 } 178 } 179 PetscCall(PetscInfo(ls,"%D function evals in line search, step = %g\n",ls->nfeval+ls->nfgeval,(double)ls->step)); 180 /* set new solution vector and compute gradient if necessary */ 181 PetscCall(VecCopy(neP->W2, x)); 182 if (ls->reason == TAOLINESEARCH_CONTINUE_ITERATING) { 183 ls->reason = TAOLINESEARCH_SUCCESS; 184 } 185 if (!g_computed) { 186 PetscCall(TaoLineSearchComputeGradient(ls,x,g)); 187 } 188 PetscFunctionReturn(0); 189 } 190 191 /* ---------------------------------------------------------- */ 192 193 /*MC 194 TAOLINESEARCHGPCG - Special line-search method for the Gradient-Projected Conjugate Gradient (TAOGPCG) algorithm. 195 Should not be used with any other algorithm. 196 197 Level: developer 198 199 .keywords: Tao, linesearch 200 M*/ 201 PETSC_EXTERN PetscErrorCode TaoLineSearchCreate_GPCG(TaoLineSearch ls) 202 { 203 TaoLineSearch_GPCG *neP; 204 205 PetscFunctionBegin; 206 ls->ftol = 0.05; 207 ls->rtol = 0.0; 208 ls->gtol = 0.0; 209 ls->stepmin = 1.0e-20; 210 ls->stepmax = 1.0e+20; 211 ls->nfeval = 0; 212 ls->max_funcs = 30; 213 ls->step = 1.0; 214 215 PetscCall(PetscNewLog(ls,&neP)); 216 neP->bracket = 0; 217 neP->infoc = 1; 218 ls->data = (void*)neP; 219 220 ls->ops->setup = NULL; 221 ls->ops->reset = NULL; 222 ls->ops->apply = TaoLineSearchApply_GPCG; 223 ls->ops->view = TaoLineSearchView_GPCG; 224 ls->ops->destroy = TaoLineSearchDestroy_GPCG; 225 ls->ops->setfromoptions = NULL; 226 ls->ops->monitor = NULL; 227 PetscFunctionReturn(0); 228 } 229