Actual source code: gpcglinesearch.c

  1: #include <petsc/private/taolinesearchimpl.h>
  2: #include <../src/tao/linesearch/impls/gpcglinesearch/gpcglinesearch.h>

  4: /* ---------------------------------------------------------- */

  6: static PetscErrorCode TaoLineSearchDestroy_GPCG(TaoLineSearch ls)
  7: {
  8:   TaoLineSearch_GPCG *ctx = (TaoLineSearch_GPCG *)ls->data;

 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: }

 19: /*------------------------------------------------------------*/
 20: static PetscErrorCode TaoLineSearchView_GPCG(TaoLineSearch ls, PetscViewer viewer)
 21: {
 22:   PetscBool isascii;

 24:   PetscFunctionBegin;
 25:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
 26:   if (isascii) PetscCall(PetscViewerASCIIPrintf(viewer, " GPCG Line search"));
 27:   PetscFunctionReturn(PETSC_SUCCESS);
 28: }

 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;

 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 */

 48:   PetscCall(TaoLineSearchMonitor(ls, 0, *f, 0.0));

 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:   }

 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;

 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:   }

 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);

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:     }

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:     }

120:     PetscCall(TaoLineSearchMonitor(ls, i + 1, *f, ls->step));

122:     if (0 == i) ls->f_fullstep = *f;

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));

128:     if (PetscAbsReal(prered) < 1.0e-100) prered = 1.0e-12;
129:     rho = actred / prered;

131:     /*
132:        If sufficient progress has been obtained, accept the
133:        point.  Otherwise, backtrack.
134:     */

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:     }

145:     /* Convergence testing */

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: }

182: /* ---------------------------------------------------------- */

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.

188:    Level: developer

190: .keywords: Tao, linesearch
191: M*/
192: PETSC_EXTERN PetscErrorCode TaoLineSearchCreate_GPCG(TaoLineSearch ls)
193: {
194:   TaoLineSearch_GPCG *neP;

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;

206:   PetscCall(PetscNew(&neP));
207:   neP->bracket = 0;
208:   neP->infoc   = 1;
209:   ls->data     = (void *)neP;

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: }