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