Actual source code: ntl.c

  1: #include <../src/tao/unconstrained/impls/ntl/ntlimpl.h>

  3: #include <petscksp.h>

  5: #define NTL_INIT_CONSTANT      0
  6: #define NTL_INIT_DIRECTION     1
  7: #define NTL_INIT_INTERPOLATION 2
  8: #define NTL_INIT_TYPES         3

 10: #define NTL_UPDATE_REDUCTION     0
 11: #define NTL_UPDATE_INTERPOLATION 1
 12: #define NTL_UPDATE_TYPES         2

 14: static const char *NTL_INIT[64] = {"constant", "direction", "interpolation"};

 16: static const char *NTL_UPDATE[64] = {"reduction", "interpolation"};

 18: /* Implements Newton's Method with a trust-region, line-search approach for
 19:    solving unconstrained minimization problems.  A More'-Thuente line search
 20:    is used to guarantee that the bfgs preconditioner remains positive
 21:    definite. */

 23: #define NTL_NEWTON          0
 24: #define NTL_BFGS            1
 25: #define NTL_SCALED_GRADIENT 2
 26: #define NTL_GRADIENT        3

 28: static PetscErrorCode TaoSolve_NTL(Tao tao)
 29: {
 30:   TAO_NTL                     *tl = (TAO_NTL *)tao->data;
 31:   KSPType                      ksp_type;
 32:   PetscBool                    is_nash, is_stcg, is_gltr, is_bfgs, is_jacobi, is_symmetric, sym_set;
 33:   KSPConvergedReason           ksp_reason;
 34:   PC                           pc;
 35:   TaoLineSearchConvergedReason ls_reason;

 37:   PetscReal fmin, ftrial, prered, actred, kappa, sigma;
 38:   PetscReal tau, tau_1, tau_2, tau_max, tau_min, max_radius;
 39:   PetscReal f, fold, gdx, gnorm;
 40:   PetscReal step = 1.0;

 42:   PetscReal norm_d = 0.0;
 43:   PetscInt  stepType;
 44:   PetscInt  its;

 46:   PetscInt bfgsUpdates = 0;
 47:   PetscInt needH;

 49:   PetscInt i_max = 5;
 50:   PetscInt j_max = 1;
 51:   PetscInt i, j, n, N;

 53:   PetscInt tr_reject;

 55:   if (tao->XL || tao->XU || tao->ops->computebounds) PetscInfo(tao, "WARNING: Variable bounds have been set but will be ignored by ntl algorithm\n");

 57:   KSPGetType(tao->ksp, &ksp_type);
 58:   PetscStrcmp(ksp_type, KSPNASH, &is_nash);
 59:   PetscStrcmp(ksp_type, KSPSTCG, &is_stcg);
 60:   PetscStrcmp(ksp_type, KSPGLTR, &is_gltr);

 63:   /* Initialize the radius and modify if it is too large or small */
 64:   tao->trust = tao->trust0;
 65:   tao->trust = PetscMax(tao->trust, tl->min_radius);
 66:   tao->trust = PetscMin(tao->trust, tl->max_radius);

 68:   /* Allocate the vectors needed for the BFGS approximation */
 69:   KSPGetPC(tao->ksp, &pc);
 70:   PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &is_bfgs);
 71:   PetscObjectTypeCompare((PetscObject)pc, PCJACOBI, &is_jacobi);
 72:   if (is_bfgs) {
 73:     tl->bfgs_pre = pc;
 74:     PCLMVMGetMatLMVM(tl->bfgs_pre, &tl->M);
 75:     VecGetLocalSize(tao->solution, &n);
 76:     VecGetSize(tao->solution, &N);
 77:     MatSetSizes(tl->M, n, n, N, N);
 78:     MatLMVMAllocate(tl->M, tao->solution, tao->gradient);
 79:     MatIsSymmetricKnown(tl->M, &sym_set, &is_symmetric);
 81:   } else if (is_jacobi) PCJacobiSetUseAbs(pc, PETSC_TRUE);

 83:   /* Check convergence criteria */
 84:   TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);
 85:   VecNorm(tao->gradient, NORM_2, &gnorm);
 87:   needH = 1;

 89:   tao->reason = TAO_CONTINUE_ITERATING;
 90:   TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its);
 91:   TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step);
 92:   PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
 93:   if (tao->reason != TAO_CONTINUE_ITERATING) return 0;

 95:   /* Initialize trust-region radius */
 96:   switch (tl->init_type) {
 97:   case NTL_INIT_CONSTANT:
 98:     /* Use the initial radius specified */
 99:     break;

101:   case NTL_INIT_INTERPOLATION:
102:     /* Use the initial radius specified */
103:     max_radius = 0.0;

105:     for (j = 0; j < j_max; ++j) {
106:       fmin  = f;
107:       sigma = 0.0;

109:       if (needH) {
110:         TaoComputeHessian(tao, tao->solution, tao->hessian, tao->hessian_pre);
111:         needH = 0;
112:       }

114:       for (i = 0; i < i_max; ++i) {
115:         VecCopy(tao->solution, tl->W);
116:         VecAXPY(tl->W, -tao->trust / gnorm, tao->gradient);

118:         TaoComputeObjective(tao, tl->W, &ftrial);
119:         if (PetscIsInfOrNanReal(ftrial)) {
120:           tau = tl->gamma1_i;
121:         } else {
122:           if (ftrial < fmin) {
123:             fmin  = ftrial;
124:             sigma = -tao->trust / gnorm;
125:           }

127:           MatMult(tao->hessian, tao->gradient, tao->stepdirection);
128:           VecDot(tao->gradient, tao->stepdirection, &prered);

130:           prered = tao->trust * (gnorm - 0.5 * tao->trust * prered / (gnorm * gnorm));
131:           actred = f - ftrial;
132:           if ((PetscAbsScalar(actred) <= tl->epsilon) && (PetscAbsScalar(prered) <= tl->epsilon)) {
133:             kappa = 1.0;
134:           } else {
135:             kappa = actred / prered;
136:           }

138:           tau_1   = tl->theta_i * gnorm * tao->trust / (tl->theta_i * gnorm * tao->trust + (1.0 - tl->theta_i) * prered - actred);
139:           tau_2   = tl->theta_i * gnorm * tao->trust / (tl->theta_i * gnorm * tao->trust - (1.0 + tl->theta_i) * prered + actred);
140:           tau_min = PetscMin(tau_1, tau_2);
141:           tau_max = PetscMax(tau_1, tau_2);

143:           if (PetscAbsScalar(kappa - (PetscReal)1.0) <= tl->mu1_i) {
144:             /* Great agreement */
145:             max_radius = PetscMax(max_radius, tao->trust);

147:             if (tau_max < 1.0) {
148:               tau = tl->gamma3_i;
149:             } else if (tau_max > tl->gamma4_i) {
150:               tau = tl->gamma4_i;
151:             } else if (tau_1 >= 1.0 && tau_1 <= tl->gamma4_i && tau_2 < 1.0) {
152:               tau = tau_1;
153:             } else if (tau_2 >= 1.0 && tau_2 <= tl->gamma4_i && tau_1 < 1.0) {
154:               tau = tau_2;
155:             } else {
156:               tau = tau_max;
157:             }
158:           } else if (PetscAbsScalar(kappa - (PetscReal)1.0) <= tl->mu2_i) {
159:             /* Good agreement */
160:             max_radius = PetscMax(max_radius, tao->trust);

162:             if (tau_max < tl->gamma2_i) {
163:               tau = tl->gamma2_i;
164:             } else if (tau_max > tl->gamma3_i) {
165:               tau = tl->gamma3_i;
166:             } else {
167:               tau = tau_max;
168:             }
169:           } else {
170:             /* Not good agreement */
171:             if (tau_min > 1.0) {
172:               tau = tl->gamma2_i;
173:             } else if (tau_max < tl->gamma1_i) {
174:               tau = tl->gamma1_i;
175:             } else if ((tau_min < tl->gamma1_i) && (tau_max >= 1.0)) {
176:               tau = tl->gamma1_i;
177:             } else if ((tau_1 >= tl->gamma1_i) && (tau_1 < 1.0) && ((tau_2 < tl->gamma1_i) || (tau_2 >= 1.0))) {
178:               tau = tau_1;
179:             } else if ((tau_2 >= tl->gamma1_i) && (tau_2 < 1.0) && ((tau_1 < tl->gamma1_i) || (tau_2 >= 1.0))) {
180:               tau = tau_2;
181:             } else {
182:               tau = tau_max;
183:             }
184:           }
185:         }
186:         tao->trust = tau * tao->trust;
187:       }

189:       if (fmin < f) {
190:         f = fmin;
191:         VecAXPY(tao->solution, sigma, tao->gradient);
192:         TaoComputeGradient(tao, tao->solution, tao->gradient);

194:         VecNorm(tao->gradient, NORM_2, &gnorm);
196:         needH = 1;

198:         TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its);
199:         TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step);
200:         PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
201:         if (tao->reason != TAO_CONTINUE_ITERATING) return 0;
202:       }
203:     }
204:     tao->trust = PetscMax(tao->trust, max_radius);

206:     /* Modify the radius if it is too large or small */
207:     tao->trust = PetscMax(tao->trust, tl->min_radius);
208:     tao->trust = PetscMin(tao->trust, tl->max_radius);
209:     break;

211:   default:
212:     /* Norm of the first direction will initialize radius */
213:     tao->trust = 0.0;
214:     break;
215:   }

217:   /* Set counter for gradient/reset steps */
218:   tl->ntrust = 0;
219:   tl->newt   = 0;
220:   tl->bfgs   = 0;
221:   tl->grad   = 0;

223:   /* Have not converged; continue with Newton method */
224:   while (tao->reason == TAO_CONTINUE_ITERATING) {
225:     /* Call general purpose update function */
226:     PetscTryTypeMethod(tao, update, tao->niter, tao->user_update);
227:     ++tao->niter;
228:     tao->ksp_its = 0;
229:     /* Compute the Hessian */
230:     if (needH) TaoComputeHessian(tao, tao->solution, tao->hessian, tao->hessian_pre);

232:     if (tl->bfgs_pre) {
233:       /* Update the limited memory preconditioner */
234:       MatLMVMUpdate(tl->M, tao->solution, tao->gradient);
235:       ++bfgsUpdates;
236:     }
237:     KSPSetOperators(tao->ksp, tao->hessian, tao->hessian_pre);
238:     /* Solve the Newton system of equations */
239:     KSPCGSetRadius(tao->ksp, tl->max_radius);
240:     KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);
241:     KSPGetIterationNumber(tao->ksp, &its);
242:     tao->ksp_its += its;
243:     tao->ksp_tot_its += its;
244:     KSPCGGetNormD(tao->ksp, &norm_d);

246:     if (0.0 == tao->trust) {
247:       /* Radius was uninitialized; use the norm of the direction */
248:       if (norm_d > 0.0) {
249:         tao->trust = norm_d;

251:         /* Modify the radius if it is too large or small */
252:         tao->trust = PetscMax(tao->trust, tl->min_radius);
253:         tao->trust = PetscMin(tao->trust, tl->max_radius);
254:       } else {
255:         /* The direction was bad; set radius to default value and re-solve
256:            the trust-region subproblem to get a direction */
257:         tao->trust = tao->trust0;

259:         /* Modify the radius if it is too large or small */
260:         tao->trust = PetscMax(tao->trust, tl->min_radius);
261:         tao->trust = PetscMin(tao->trust, tl->max_radius);

263:         KSPCGSetRadius(tao->ksp, tl->max_radius);
264:         KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);
265:         KSPGetIterationNumber(tao->ksp, &its);
266:         tao->ksp_its += its;
267:         tao->ksp_tot_its += its;
268:         KSPCGGetNormD(tao->ksp, &norm_d);

271:       }
272:     }

274:     VecScale(tao->stepdirection, -1.0);
275:     KSPGetConvergedReason(tao->ksp, &ksp_reason);
276:     if ((KSP_DIVERGED_INDEFINITE_PC == ksp_reason) && (tl->bfgs_pre)) {
277:       /* Preconditioner is numerically indefinite; reset the
278:          approximate if using BFGS preconditioning. */
279:       MatLMVMReset(tl->M, PETSC_FALSE);
280:       MatLMVMUpdate(tl->M, tao->solution, tao->gradient);
281:       bfgsUpdates = 1;
282:     }

284:     /* Check trust-region reduction conditions */
285:     tr_reject = 0;
286:     if (NTL_UPDATE_REDUCTION == tl->update_type) {
287:       /* Get predicted reduction */
288:       KSPCGGetObjFcn(tao->ksp, &prered);
289:       if (prered >= 0.0) {
290:         /* The predicted reduction has the wrong sign.  This cannot
291:            happen in infinite precision arithmetic.  Step should
292:            be rejected! */
293:         tao->trust = tl->alpha1 * PetscMin(tao->trust, norm_d);
294:         tr_reject  = 1;
295:       } else {
296:         /* Compute trial step and function value */
297:         VecCopy(tao->solution, tl->W);
298:         VecAXPY(tl->W, 1.0, tao->stepdirection);
299:         TaoComputeObjective(tao, tl->W, &ftrial);

301:         if (PetscIsInfOrNanReal(ftrial)) {
302:           tao->trust = tl->alpha1 * PetscMin(tao->trust, norm_d);
303:           tr_reject  = 1;
304:         } else {
305:           /* Compute and actual reduction */
306:           actred = f - ftrial;
307:           prered = -prered;
308:           if ((PetscAbsScalar(actred) <= tl->epsilon) && (PetscAbsScalar(prered) <= tl->epsilon)) {
309:             kappa = 1.0;
310:           } else {
311:             kappa = actred / prered;
312:           }

314:           /* Accept of reject the step and update radius */
315:           if (kappa < tl->eta1) {
316:             /* Reject the step */
317:             tao->trust = tl->alpha1 * PetscMin(tao->trust, norm_d);
318:             tr_reject  = 1;
319:           } else {
320:             /* Accept the step */
321:             if (kappa < tl->eta2) {
322:               /* Marginal bad step */
323:               tao->trust = tl->alpha2 * PetscMin(tao->trust, norm_d);
324:             } else if (kappa < tl->eta3) {
325:               /* Reasonable step */
326:               tao->trust = tl->alpha3 * tao->trust;
327:             } else if (kappa < tl->eta4) {
328:               /* Good step */
329:               tao->trust = PetscMax(tl->alpha4 * norm_d, tao->trust);
330:             } else {
331:               /* Very good step */
332:               tao->trust = PetscMax(tl->alpha5 * norm_d, tao->trust);
333:             }
334:           }
335:         }
336:       }
337:     } else {
338:       /* Get predicted reduction */
339:       KSPCGGetObjFcn(tao->ksp, &prered);
340:       if (prered >= 0.0) {
341:         /* The predicted reduction has the wrong sign.  This cannot
342:            happen in infinite precision arithmetic.  Step should
343:            be rejected! */
344:         tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d);
345:         tr_reject  = 1;
346:       } else {
347:         VecCopy(tao->solution, tl->W);
348:         VecAXPY(tl->W, 1.0, tao->stepdirection);
349:         TaoComputeObjective(tao, tl->W, &ftrial);
350:         if (PetscIsInfOrNanReal(ftrial)) {
351:           tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d);
352:           tr_reject  = 1;
353:         } else {
354:           VecDot(tao->gradient, tao->stepdirection, &gdx);

356:           actred = f - ftrial;
357:           prered = -prered;
358:           if ((PetscAbsScalar(actred) <= tl->epsilon) && (PetscAbsScalar(prered) <= tl->epsilon)) {
359:             kappa = 1.0;
360:           } else {
361:             kappa = actred / prered;
362:           }

364:           tau_1   = tl->theta * gdx / (tl->theta * gdx - (1.0 - tl->theta) * prered + actred);
365:           tau_2   = tl->theta * gdx / (tl->theta * gdx + (1.0 + tl->theta) * prered - actred);
366:           tau_min = PetscMin(tau_1, tau_2);
367:           tau_max = PetscMax(tau_1, tau_2);

369:           if (kappa >= 1.0 - tl->mu1) {
370:             /* Great agreement; accept step and update radius */
371:             if (tau_max < 1.0) {
372:               tao->trust = PetscMax(tao->trust, tl->gamma3 * norm_d);
373:             } else if (tau_max > tl->gamma4) {
374:               tao->trust = PetscMax(tao->trust, tl->gamma4 * norm_d);
375:             } else {
376:               tao->trust = PetscMax(tao->trust, tau_max * norm_d);
377:             }
378:           } else if (kappa >= 1.0 - tl->mu2) {
379:             /* Good agreement */

381:             if (tau_max < tl->gamma2) {
382:               tao->trust = tl->gamma2 * PetscMin(tao->trust, norm_d);
383:             } else if (tau_max > tl->gamma3) {
384:               tao->trust = PetscMax(tao->trust, tl->gamma3 * norm_d);
385:             } else if (tau_max < 1.0) {
386:               tao->trust = tau_max * PetscMin(tao->trust, norm_d);
387:             } else {
388:               tao->trust = PetscMax(tao->trust, tau_max * norm_d);
389:             }
390:           } else {
391:             /* Not good agreement */
392:             if (tau_min > 1.0) {
393:               tao->trust = tl->gamma2 * PetscMin(tao->trust, norm_d);
394:             } else if (tau_max < tl->gamma1) {
395:               tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d);
396:             } else if ((tau_min < tl->gamma1) && (tau_max >= 1.0)) {
397:               tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d);
398:             } else if ((tau_1 >= tl->gamma1) && (tau_1 < 1.0) && ((tau_2 < tl->gamma1) || (tau_2 >= 1.0))) {
399:               tao->trust = tau_1 * PetscMin(tao->trust, norm_d);
400:             } else if ((tau_2 >= tl->gamma1) && (tau_2 < 1.0) && ((tau_1 < tl->gamma1) || (tau_2 >= 1.0))) {
401:               tao->trust = tau_2 * PetscMin(tao->trust, norm_d);
402:             } else {
403:               tao->trust = tau_max * PetscMin(tao->trust, norm_d);
404:             }
405:             tr_reject = 1;
406:           }
407:         }
408:       }
409:     }

411:     if (tr_reject) {
412:       /* The trust-region constraints rejected the step.  Apply a linesearch.
413:          Check for descent direction. */
414:       VecDot(tao->stepdirection, tao->gradient, &gdx);
415:       if ((gdx >= 0.0) || PetscIsInfOrNanReal(gdx)) {
416:         /* Newton step is not descent or direction produced Inf or NaN */

418:         if (!tl->bfgs_pre) {
419:           /* We don't have the bfgs matrix around and updated
420:              Must use gradient direction in this case */
421:           VecCopy(tao->gradient, tao->stepdirection);
422:           VecScale(tao->stepdirection, -1.0);
423:           ++tl->grad;
424:           stepType = NTL_GRADIENT;
425:         } else {
426:           /* Attempt to use the BFGS direction */
427:           MatSolve(tl->M, tao->gradient, tao->stepdirection);
428:           VecScale(tao->stepdirection, -1.0);

430:           /* Check for success (descent direction) */
431:           VecDot(tao->stepdirection, tao->gradient, &gdx);
432:           if ((gdx >= 0) || PetscIsInfOrNanReal(gdx)) {
433:             /* BFGS direction is not descent or direction produced not a number
434:                We can assert bfgsUpdates > 1 in this case because
435:                the first solve produces the scaled gradient direction,
436:                which is guaranteed to be descent */

438:             /* Use steepest descent direction (scaled) */
439:             MatLMVMReset(tl->M, PETSC_FALSE);
440:             MatLMVMUpdate(tl->M, tao->solution, tao->gradient);
441:             MatSolve(tl->M, tao->gradient, tao->stepdirection);
442:             VecScale(tao->stepdirection, -1.0);

444:             bfgsUpdates = 1;
445:             ++tl->grad;
446:             stepType = NTL_GRADIENT;
447:           } else {
448:             MatLMVMGetUpdateCount(tl->M, &bfgsUpdates);
449:             if (1 == bfgsUpdates) {
450:               /* The first BFGS direction is always the scaled gradient */
451:               ++tl->grad;
452:               stepType = NTL_GRADIENT;
453:             } else {
454:               ++tl->bfgs;
455:               stepType = NTL_BFGS;
456:             }
457:           }
458:         }
459:       } else {
460:         /* Computed Newton step is descent */
461:         ++tl->newt;
462:         stepType = NTL_NEWTON;
463:       }

465:       /* Perform the linesearch */
466:       fold = f;
467:       VecCopy(tao->solution, tl->Xold);
468:       VecCopy(tao->gradient, tl->Gold);

470:       step = 1.0;
471:       TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_reason);
472:       TaoAddLineSearchCounts(tao);

474:       while (ls_reason != TAOLINESEARCH_SUCCESS && ls_reason != TAOLINESEARCH_SUCCESS_USER && stepType != NTL_GRADIENT) { /* Linesearch failed */
475:         /* Linesearch failed */
476:         f = fold;
477:         VecCopy(tl->Xold, tao->solution);
478:         VecCopy(tl->Gold, tao->gradient);

480:         switch (stepType) {
481:         case NTL_NEWTON:
482:           /* Failed to obtain acceptable iterate with Newton step */

484:           if (tl->bfgs_pre) {
485:             /* We don't have the bfgs matrix around and being updated
486:                Must use gradient direction in this case */
487:             VecCopy(tao->gradient, tao->stepdirection);
488:             ++tl->grad;
489:             stepType = NTL_GRADIENT;
490:           } else {
491:             /* Attempt to use the BFGS direction */
492:             MatSolve(tl->M, tao->gradient, tao->stepdirection);

494:             /* Check for success (descent direction) */
495:             VecDot(tao->stepdirection, tao->gradient, &gdx);
496:             if ((gdx <= 0) || PetscIsInfOrNanReal(gdx)) {
497:               /* BFGS direction is not descent or direction produced
498:                  not a number.  We can assert bfgsUpdates > 1 in this case
499:                  Use steepest descent direction (scaled) */
500:               MatLMVMReset(tl->M, PETSC_FALSE);
501:               MatLMVMUpdate(tl->M, tao->solution, tao->gradient);
502:               MatSolve(tl->M, tao->gradient, tao->stepdirection);

504:               bfgsUpdates = 1;
505:               ++tl->grad;
506:               stepType = NTL_GRADIENT;
507:             } else {
508:               MatLMVMGetUpdateCount(tl->M, &bfgsUpdates);
509:               if (1 == bfgsUpdates) {
510:                 /* The first BFGS direction is always the scaled gradient */
511:                 ++tl->grad;
512:                 stepType = NTL_GRADIENT;
513:               } else {
514:                 ++tl->bfgs;
515:                 stepType = NTL_BFGS;
516:               }
517:             }
518:           }
519:           break;

521:         case NTL_BFGS:
522:           /* Can only enter if pc_type == NTL_PC_BFGS
523:              Failed to obtain acceptable iterate with BFGS step
524:              Attempt to use the scaled gradient direction */
525:           MatLMVMReset(tl->M, PETSC_FALSE);
526:           MatLMVMUpdate(tl->M, tao->solution, tao->gradient);
527:           MatSolve(tl->M, tao->gradient, tao->stepdirection);

529:           bfgsUpdates = 1;
530:           ++tl->grad;
531:           stepType = NTL_GRADIENT;
532:           break;
533:         }
534:         VecScale(tao->stepdirection, -1.0);

536:         /* This may be incorrect; linesearch has values for stepmax and stepmin
537:            that should be reset. */
538:         step = 1.0;
539:         TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_reason);
540:         TaoAddLineSearchCounts(tao);
541:       }

543:       if (ls_reason != TAOLINESEARCH_SUCCESS && ls_reason != TAOLINESEARCH_SUCCESS_USER) {
544:         /* Failed to find an improving point */
545:         f = fold;
546:         VecCopy(tl->Xold, tao->solution);
547:         VecCopy(tl->Gold, tao->gradient);
548:         tao->trust  = 0.0;
549:         step        = 0.0;
550:         tao->reason = TAO_DIVERGED_LS_FAILURE;
551:         break;
552:       } else if (stepType == NTL_NEWTON) {
553:         if (step < tl->nu1) {
554:           /* Very bad step taken; reduce radius */
555:           tao->trust = tl->omega1 * PetscMin(norm_d, tao->trust);
556:         } else if (step < tl->nu2) {
557:           /* Reasonably bad step taken; reduce radius */
558:           tao->trust = tl->omega2 * PetscMin(norm_d, tao->trust);
559:         } else if (step < tl->nu3) {
560:           /* Reasonable step was taken; leave radius alone */
561:           if (tl->omega3 < 1.0) {
562:             tao->trust = tl->omega3 * PetscMin(norm_d, tao->trust);
563:           } else if (tl->omega3 > 1.0) {
564:             tao->trust = PetscMax(tl->omega3 * norm_d, tao->trust);
565:           }
566:         } else if (step < tl->nu4) {
567:           /* Full step taken; increase the radius */
568:           tao->trust = PetscMax(tl->omega4 * norm_d, tao->trust);
569:         } else {
570:           /* More than full step taken; increase the radius */
571:           tao->trust = PetscMax(tl->omega5 * norm_d, tao->trust);
572:         }
573:       } else {
574:         /* Newton step was not good; reduce the radius */
575:         tao->trust = tl->omega1 * PetscMin(norm_d, tao->trust);
576:       }
577:     } else {
578:       /* Trust-region step is accepted */
579:       VecCopy(tl->W, tao->solution);
580:       f = ftrial;
581:       TaoComputeGradient(tao, tao->solution, tao->gradient);
582:       ++tl->ntrust;
583:     }

585:     /* The radius may have been increased; modify if it is too large */
586:     tao->trust = PetscMin(tao->trust, tl->max_radius);

588:     /* Check for converged */
589:     VecNorm(tao->gradient, NORM_2, &gnorm);
591:     needH = 1;

593:     TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its);
594:     TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step);
595:     PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
596:   }
597:   return 0;
598: }

600: /* ---------------------------------------------------------- */
601: static PetscErrorCode TaoSetUp_NTL(Tao tao)
602: {
603:   TAO_NTL *tl = (TAO_NTL *)tao->data;

605:   if (!tao->gradient) VecDuplicate(tao->solution, &tao->gradient);
606:   if (!tao->stepdirection) VecDuplicate(tao->solution, &tao->stepdirection);
607:   if (!tl->W) VecDuplicate(tao->solution, &tl->W);
608:   if (!tl->Xold) VecDuplicate(tao->solution, &tl->Xold);
609:   if (!tl->Gold) VecDuplicate(tao->solution, &tl->Gold);
610:   tl->bfgs_pre = NULL;
611:   tl->M        = NULL;
612:   return 0;
613: }

615: /*------------------------------------------------------------*/
616: static PetscErrorCode TaoDestroy_NTL(Tao tao)
617: {
618:   TAO_NTL *tl = (TAO_NTL *)tao->data;

620:   if (tao->setupcalled) {
621:     VecDestroy(&tl->W);
622:     VecDestroy(&tl->Xold);
623:     VecDestroy(&tl->Gold);
624:   }
625:   KSPDestroy(&tao->ksp);
626:   PetscFree(tao->data);
627:   return 0;
628: }

630: /*------------------------------------------------------------*/
631: static PetscErrorCode TaoSetFromOptions_NTL(Tao tao, PetscOptionItems *PetscOptionsObject)
632: {
633:   TAO_NTL *tl = (TAO_NTL *)tao->data;

635:   PetscOptionsHeadBegin(PetscOptionsObject, "Newton trust region with line search method for unconstrained optimization");
636:   PetscOptionsEList("-tao_ntl_init_type", "radius initialization type", "", NTL_INIT, NTL_INIT_TYPES, NTL_INIT[tl->init_type], &tl->init_type, NULL);
637:   PetscOptionsEList("-tao_ntl_update_type", "radius update type", "", NTL_UPDATE, NTL_UPDATE_TYPES, NTL_UPDATE[tl->update_type], &tl->update_type, NULL);
638:   PetscOptionsReal("-tao_ntl_eta1", "poor steplength; reduce radius", "", tl->eta1, &tl->eta1, NULL);
639:   PetscOptionsReal("-tao_ntl_eta2", "reasonable steplength; leave radius alone", "", tl->eta2, &tl->eta2, NULL);
640:   PetscOptionsReal("-tao_ntl_eta3", "good steplength; increase radius", "", tl->eta3, &tl->eta3, NULL);
641:   PetscOptionsReal("-tao_ntl_eta4", "excellent steplength; greatly increase radius", "", tl->eta4, &tl->eta4, NULL);
642:   PetscOptionsReal("-tao_ntl_alpha1", "", "", tl->alpha1, &tl->alpha1, NULL);
643:   PetscOptionsReal("-tao_ntl_alpha2", "", "", tl->alpha2, &tl->alpha2, NULL);
644:   PetscOptionsReal("-tao_ntl_alpha3", "", "", tl->alpha3, &tl->alpha3, NULL);
645:   PetscOptionsReal("-tao_ntl_alpha4", "", "", tl->alpha4, &tl->alpha4, NULL);
646:   PetscOptionsReal("-tao_ntl_alpha5", "", "", tl->alpha5, &tl->alpha5, NULL);
647:   PetscOptionsReal("-tao_ntl_nu1", "poor steplength; reduce radius", "", tl->nu1, &tl->nu1, NULL);
648:   PetscOptionsReal("-tao_ntl_nu2", "reasonable steplength; leave radius alone", "", tl->nu2, &tl->nu2, NULL);
649:   PetscOptionsReal("-tao_ntl_nu3", "good steplength; increase radius", "", tl->nu3, &tl->nu3, NULL);
650:   PetscOptionsReal("-tao_ntl_nu4", "excellent steplength; greatly increase radius", "", tl->nu4, &tl->nu4, NULL);
651:   PetscOptionsReal("-tao_ntl_omega1", "", "", tl->omega1, &tl->omega1, NULL);
652:   PetscOptionsReal("-tao_ntl_omega2", "", "", tl->omega2, &tl->omega2, NULL);
653:   PetscOptionsReal("-tao_ntl_omega3", "", "", tl->omega3, &tl->omega3, NULL);
654:   PetscOptionsReal("-tao_ntl_omega4", "", "", tl->omega4, &tl->omega4, NULL);
655:   PetscOptionsReal("-tao_ntl_omega5", "", "", tl->omega5, &tl->omega5, NULL);
656:   PetscOptionsReal("-tao_ntl_mu1_i", "", "", tl->mu1_i, &tl->mu1_i, NULL);
657:   PetscOptionsReal("-tao_ntl_mu2_i", "", "", tl->mu2_i, &tl->mu2_i, NULL);
658:   PetscOptionsReal("-tao_ntl_gamma1_i", "", "", tl->gamma1_i, &tl->gamma1_i, NULL);
659:   PetscOptionsReal("-tao_ntl_gamma2_i", "", "", tl->gamma2_i, &tl->gamma2_i, NULL);
660:   PetscOptionsReal("-tao_ntl_gamma3_i", "", "", tl->gamma3_i, &tl->gamma3_i, NULL);
661:   PetscOptionsReal("-tao_ntl_gamma4_i", "", "", tl->gamma4_i, &tl->gamma4_i, NULL);
662:   PetscOptionsReal("-tao_ntl_theta_i", "", "", tl->theta_i, &tl->theta_i, NULL);
663:   PetscOptionsReal("-tao_ntl_mu1", "", "", tl->mu1, &tl->mu1, NULL);
664:   PetscOptionsReal("-tao_ntl_mu2", "", "", tl->mu2, &tl->mu2, NULL);
665:   PetscOptionsReal("-tao_ntl_gamma1", "", "", tl->gamma1, &tl->gamma1, NULL);
666:   PetscOptionsReal("-tao_ntl_gamma2", "", "", tl->gamma2, &tl->gamma2, NULL);
667:   PetscOptionsReal("-tao_ntl_gamma3", "", "", tl->gamma3, &tl->gamma3, NULL);
668:   PetscOptionsReal("-tao_ntl_gamma4", "", "", tl->gamma4, &tl->gamma4, NULL);
669:   PetscOptionsReal("-tao_ntl_theta", "", "", tl->theta, &tl->theta, NULL);
670:   PetscOptionsReal("-tao_ntl_min_radius", "lower bound on initial radius", "", tl->min_radius, &tl->min_radius, NULL);
671:   PetscOptionsReal("-tao_ntl_max_radius", "upper bound on radius", "", tl->max_radius, &tl->max_radius, NULL);
672:   PetscOptionsReal("-tao_ntl_epsilon", "tolerance used when computing actual and predicted reduction", "", tl->epsilon, &tl->epsilon, NULL);
673:   PetscOptionsHeadEnd();
674:   TaoLineSearchSetFromOptions(tao->linesearch);
675:   KSPSetFromOptions(tao->ksp);
676:   return 0;
677: }

679: /*------------------------------------------------------------*/
680: static PetscErrorCode TaoView_NTL(Tao tao, PetscViewer viewer)
681: {
682:   TAO_NTL  *tl = (TAO_NTL *)tao->data;
683:   PetscBool isascii;

685:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
686:   if (isascii) {
687:     PetscViewerASCIIPushTab(viewer);
688:     PetscViewerASCIIPrintf(viewer, "Trust-region steps: %" PetscInt_FMT "\n", tl->ntrust);
689:     PetscViewerASCIIPrintf(viewer, "Newton search steps: %" PetscInt_FMT "\n", tl->newt);
690:     PetscViewerASCIIPrintf(viewer, "BFGS search steps: %" PetscInt_FMT "\n", tl->bfgs);
691:     PetscViewerASCIIPrintf(viewer, "Gradient search steps: %" PetscInt_FMT "\n", tl->grad);
692:     PetscViewerASCIIPopTab(viewer);
693:   }
694:   return 0;
695: }

697: /* ---------------------------------------------------------- */
698: /*MC
699:   TAONTL - Newton's method with trust region globalization and line search fallback.
700:   At each iteration, the Newton trust region method solves the system for d
701:   and performs a line search in the d direction:

703:             min_d  .5 dT Hk d + gkT d,  s.t.   ||d|| < Delta_k

705:   Options Database Keys:
706: + -tao_ntl_init_type - "constant","direction","interpolation"
707: . -tao_ntl_update_type - "reduction","interpolation"
708: . -tao_ntl_min_radius - lower bound on trust region radius
709: . -tao_ntl_max_radius - upper bound on trust region radius
710: . -tao_ntl_epsilon - tolerance for accepting actual / predicted reduction
711: . -tao_ntl_mu1_i - mu1 interpolation init factor
712: . -tao_ntl_mu2_i - mu2 interpolation init factor
713: . -tao_ntl_gamma1_i - gamma1 interpolation init factor
714: . -tao_ntl_gamma2_i - gamma2 interpolation init factor
715: . -tao_ntl_gamma3_i - gamma3 interpolation init factor
716: . -tao_ntl_gamma4_i - gamma4 interpolation init factor
717: . -tao_ntl_theta_i - theta1 interpolation init factor
718: . -tao_ntl_eta1 - eta1 reduction update factor
719: . -tao_ntl_eta2 - eta2 reduction update factor
720: . -tao_ntl_eta3 - eta3 reduction update factor
721: . -tao_ntl_eta4 - eta4 reduction update factor
722: . -tao_ntl_alpha1 - alpha1 reduction update factor
723: . -tao_ntl_alpha2 - alpha2 reduction update factor
724: . -tao_ntl_alpha3 - alpha3 reduction update factor
725: . -tao_ntl_alpha4 - alpha4 reduction update factor
726: . -tao_ntl_alpha4 - alpha4 reduction update factor
727: . -tao_ntl_mu1 - mu1 interpolation update
728: . -tao_ntl_mu2 - mu2 interpolation update
729: . -tao_ntl_gamma1 - gamma1 interpolcation update
730: . -tao_ntl_gamma2 - gamma2 interpolcation update
731: . -tao_ntl_gamma3 - gamma3 interpolcation update
732: . -tao_ntl_gamma4 - gamma4 interpolation update
733: - -tao_ntl_theta - theta1 interpolation update

735:   Level: beginner
736: M*/
737: PETSC_EXTERN PetscErrorCode TaoCreate_NTL(Tao tao)
738: {
739:   TAO_NTL    *tl;
740:   const char *morethuente_type = TAOLINESEARCHMT;

742:   PetscNew(&tl);
743:   tao->ops->setup          = TaoSetUp_NTL;
744:   tao->ops->solve          = TaoSolve_NTL;
745:   tao->ops->view           = TaoView_NTL;
746:   tao->ops->setfromoptions = TaoSetFromOptions_NTL;
747:   tao->ops->destroy        = TaoDestroy_NTL;

749:   /* Override default settings (unless already changed) */
750:   if (!tao->max_it_changed) tao->max_it = 50;
751:   if (!tao->trust0_changed) tao->trust0 = 100.0;

753:   tao->data = (void *)tl;

755:   /* Default values for trust-region radius update based on steplength */
756:   tl->nu1 = 0.25;
757:   tl->nu2 = 0.50;
758:   tl->nu3 = 1.00;
759:   tl->nu4 = 1.25;

761:   tl->omega1 = 0.25;
762:   tl->omega2 = 0.50;
763:   tl->omega3 = 1.00;
764:   tl->omega4 = 2.00;
765:   tl->omega5 = 4.00;

767:   /* Default values for trust-region radius update based on reduction */
768:   tl->eta1 = 1.0e-4;
769:   tl->eta2 = 0.25;
770:   tl->eta3 = 0.50;
771:   tl->eta4 = 0.90;

773:   tl->alpha1 = 0.25;
774:   tl->alpha2 = 0.50;
775:   tl->alpha3 = 1.00;
776:   tl->alpha4 = 2.00;
777:   tl->alpha5 = 4.00;

779:   /* Default values for trust-region radius update based on interpolation */
780:   tl->mu1 = 0.10;
781:   tl->mu2 = 0.50;

783:   tl->gamma1 = 0.25;
784:   tl->gamma2 = 0.50;
785:   tl->gamma3 = 2.00;
786:   tl->gamma4 = 4.00;

788:   tl->theta = 0.05;

790:   /* Default values for trust region initialization based on interpolation */
791:   tl->mu1_i = 0.35;
792:   tl->mu2_i = 0.50;

794:   tl->gamma1_i = 0.0625;
795:   tl->gamma2_i = 0.5;
796:   tl->gamma3_i = 2.0;
797:   tl->gamma4_i = 5.0;

799:   tl->theta_i = 0.25;

801:   /* Remaining parameters */
802:   tl->min_radius = 1.0e-10;
803:   tl->max_radius = 1.0e10;
804:   tl->epsilon    = 1.0e-6;

806:   tl->init_type   = NTL_INIT_INTERPOLATION;
807:   tl->update_type = NTL_UPDATE_REDUCTION;

809:   TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);
810:   PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1);
811:   TaoLineSearchSetType(tao->linesearch, morethuente_type);
812:   TaoLineSearchUseTaoRoutines(tao->linesearch, tao);
813:   TaoLineSearchSetOptionsPrefix(tao->linesearch, tao->hdr.prefix);
814:   KSPCreate(((PetscObject)tao)->comm, &tao->ksp);
815:   PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1);
816:   KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix);
817:   KSPAppendOptionsPrefix(tao->ksp, "tao_ntl_");
818:   KSPSetType(tao->ksp, KSPSTCG);
819:   return 0;
820: }