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