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:   VecDestroy(&ctx->W1);
 11:   VecDestroy(&ctx->W2);
 12:   VecDestroy(&ctx->Gold);
 13:   VecDestroy(&ctx->x);
 14:   PetscFree(ls->data);
 15:   return 0;
 16: }

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

 23:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
 24:   if (isascii) PetscViewerASCIIPrintf(viewer, " GPCG Line search");
 25:   return 0;
 26: }

 28: /*------------------------------------------------------------*/
 29: static PetscErrorCode TaoLineSearchApply_GPCG(TaoLineSearch ls, Vec x, PetscReal *f, Vec g, Vec s)
 30: {
 31:   TaoLineSearch_GPCG *neP = (TaoLineSearch_GPCG *)ls->data;
 32:   PetscInt            i;
 33:   PetscBool           g_computed = PETSC_FALSE; /* to prevent extra gradient computation */
 34:   PetscReal           d1, finit, actred, prered, rho, gdx;

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

 45:   TaoLineSearchMonitor(ls, 0, *f, 0.0);

 47:   ls->reason = TAOLINESEARCH_CONTINUE_ITERATING;
 48:   ls->step   = ls->initstep;
 49:   if (!neP->W2) {
 50:     VecDuplicate(x, &neP->W2);
 51:     VecDuplicate(x, &neP->W1);
 52:     VecDuplicate(x, &neP->Gold);
 53:     neP->x = x;
 54:     PetscObjectReference((PetscObject)neP->x);
 55:   } else if (x != neP->x) {
 56:     VecDestroy(&neP->x);
 57:     VecDestroy(&neP->W1);
 58:     VecDestroy(&neP->W2);
 59:     VecDestroy(&neP->Gold);
 60:     VecDuplicate(x, &neP->W1);
 61:     VecDuplicate(x, &neP->W2);
 62:     VecDuplicate(x, &neP->Gold);
 63:     PetscObjectDereference((PetscObject)neP->x);
 64:     neP->x = x;
 65:     PetscObjectReference((PetscObject)neP->x);
 66:   }

 68:   VecDot(g, s, &gdx);
 69:   if (gdx > 0) {
 70:     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:     return 0;
 73:   }
 74:   VecCopy(x, neP->W2);
 75:   VecCopy(g, neP->Gold);
 76:   if (ls->bounded) {
 77:     /* Compute the smallest steplength that will make one nonbinding variable  equal the bound */
 78:     VecStepBoundInfo(x, s, ls->lower, ls->upper, &rho, &actred, &d1);
 79:     ls->step = PetscMin(ls->step, d1);
 80:   }
 81:   rho    = 0;
 82:   actred = 0;

 84:   if (ls->step < 0) {
 85:     PetscInfo(ls, "Line search error: initial step parameter %g< 0\n", (double)ls->step);
 86:     ls->reason = TAOLINESEARCH_HALTED_OTHER;
 87:     return 0;
 88:   }

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

 97:     VecWAXPY(neP->W2, ls->step, s, x);
 98:     if (ls->bounded) {
 99:       /* Make sure new vector is numerically within bounds */
100:       VecMedian(neP->W2, ls->lower, ls->upper, neP->W2);
101:     }

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:       TaoLineSearchComputeObjective(ls, neP->W2, f);
108:       g_computed = PETSC_FALSE;
109:     } else if (ls->usegts) {
110:       TaoLineSearchComputeObjectiveAndGTS(ls, neP->W2, f, &gdx);
111:       g_computed = PETSC_FALSE;
112:     } else {
113:       TaoLineSearchComputeObjectiveAndGradient(ls, neP->W2, f, g);
114:       g_computed = PETSC_TRUE;
115:     }

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

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

121:     actred = *f - finit;
122:     VecWAXPY(neP->W1, -1.0, x, neP->W2); /* W1 = W2 - X */
123:     VecDot(neP->W1, neP->Gold, &prered);

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

128:     /*
129:        If sufficient progress has been obtained, accept the
130:        point.  Otherwise, backtrack.
131:     */

133:     if (actred > 0) {
134:       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:     }

142:     /* Convergence testing */

144:     if (ls->step <= ls->stepmin || ls->step >= ls->stepmax) {
145:       ls->reason = TAOLINESEARCH_HALTED_OTHER;
146:       PetscInfo(ls, "Rounding errors may prevent further progress.  May not be a step satisfying\n");
147:       PetscInfo(ls, "sufficient decrease and curvature conditions. Tolerances may be too small.\n");
148:       break;
149:     }
150:     if (ls->step == ls->stepmax) {
151:       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:       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:       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:       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:   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:   VecCopy(neP->W2, x);
174:   if (ls->reason == TAOLINESEARCH_CONTINUE_ITERATING) ls->reason = TAOLINESEARCH_SUCCESS;
175:   if (!g_computed) TaoLineSearchComputeGradient(ls, x, g);
176:   return 0;
177: }

179: /* ---------------------------------------------------------- */

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.

185:    Level: developer

187: .keywords: Tao, linesearch
188: M*/
189: PETSC_EXTERN PetscErrorCode TaoLineSearchCreate_GPCG(TaoLineSearch ls)
190: {
191:   TaoLineSearch_GPCG *neP;

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;

202:   PetscNew(&neP);
203:   neP->bracket = 0;
204:   neP->infoc   = 1;
205:   ls->data     = (void *)neP;

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:   return 0;
215: }