Actual source code: ntr.c

  1: #include <../src/tao/unconstrained/impls/ntr/ntrimpl.h>

  3: #include <petscksp.h>

  5: #define NTR_INIT_CONSTANT      0
  6: #define NTR_INIT_DIRECTION     1
  7: #define NTR_INIT_INTERPOLATION 2
  8: #define NTR_INIT_TYPES         3

 10: #define NTR_UPDATE_REDUCTION     0
 11: #define NTR_UPDATE_INTERPOLATION 1
 12: #define NTR_UPDATE_TYPES         2

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

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

 18: /*
 19:    TaoSolve_NTR - Implements Newton's Method with a trust region approach
 20:    for solving unconstrained minimization problems.

 22:    The basic algorithm is taken from MINPACK-2 (dstrn).

 24:    TaoSolve_NTR computes a local minimizer of a twice differentiable function
 25:    f by applying a trust region variant of Newton's method.  At each stage
 26:    of the algorithm, we use the prconditioned conjugate gradient method to
 27:    determine an approximate minimizer of the quadratic equation

 29:         q(s) = <s, Hs + g>

 31:    subject to the trust region constraint

 33:         || s ||_M <= radius,

 35:    where radius is the trust region radius and M is a symmetric positive
 36:    definite matrix (the preconditioner).  Here g is the gradient and H
 37:    is the Hessian matrix.

 39:    Note:  TaoSolve_NTR MUST use the iterative solver KSPNASH, KSPSTCG,
 40:           or KSPGLTR.  Thus, we set KSPNASH, KSPSTCG, or KSPGLTR in this
 41:           routine regardless of what the user may have previously specified.
 42: */
 43: static PetscErrorCode TaoSolve_NTR(Tao tao)
 44: {
 45:   TAO_NTR           *tr = (TAO_NTR *)tao->data;
 46:   KSPType            ksp_type;
 47:   PetscBool          is_nash, is_stcg, is_gltr, is_bfgs, is_jacobi, is_symmetric, sym_set;
 48:   KSPConvergedReason ksp_reason;
 49:   PC                 pc;
 50:   PetscReal          fmin, ftrial, prered, actred, kappa, sigma, beta;
 51:   PetscReal          tau, tau_1, tau_2, tau_max, tau_min, max_radius;
 52:   PetscReal          f, gnorm;

 54:   PetscReal norm_d;
 55:   PetscInt  needH;

 57:   PetscInt i_max = 5;
 58:   PetscInt j_max = 1;
 59:   PetscInt i, j, N, n, its;

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

 63:   KSPGetType(tao->ksp, &ksp_type);
 64:   PetscStrcmp(ksp_type, KSPNASH, &is_nash);
 65:   PetscStrcmp(ksp_type, KSPSTCG, &is_stcg);
 66:   PetscStrcmp(ksp_type, KSPGLTR, &is_gltr);

 69:   /* Initialize the radius and modify if it is too large or small */
 70:   tao->trust = tao->trust0;
 71:   tao->trust = PetscMax(tao->trust, tr->min_radius);
 72:   tao->trust = PetscMin(tao->trust, tr->max_radius);

 74:   /* Allocate the vectors needed for the BFGS approximation */
 75:   KSPGetPC(tao->ksp, &pc);
 76:   PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &is_bfgs);
 77:   PetscObjectTypeCompare((PetscObject)pc, PCJACOBI, &is_jacobi);
 78:   if (is_bfgs) {
 79:     tr->bfgs_pre = pc;
 80:     PCLMVMGetMatLMVM(tr->bfgs_pre, &tr->M);
 81:     VecGetLocalSize(tao->solution, &n);
 82:     VecGetSize(tao->solution, &N);
 83:     MatSetSizes(tr->M, n, n, N, N);
 84:     MatLMVMAllocate(tr->M, tao->solution, tao->gradient);
 85:     MatIsSymmetricKnown(tr->M, &sym_set, &is_symmetric);
 87:   } else if (is_jacobi) PCJacobiSetUseAbs(pc, PETSC_TRUE);

 89:   /* Check convergence criteria */
 90:   TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);
 91:   TaoGradientNorm(tao, tao->gradient, NORM_2, &gnorm);
 93:   needH = 1;

 95:   tao->reason = TAO_CONTINUE_ITERATING;
 96:   TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its);
 97:   TaoMonitor(tao, tao->niter, f, gnorm, 0.0, 1.0);
 98:   PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
 99:   if (tao->reason != TAO_CONTINUE_ITERATING) return 0;

101:   /*  Initialize trust-region radius */
102:   switch (tr->init_type) {
103:   case NTR_INIT_CONSTANT:
104:     /*  Use the initial radius specified */
105:     break;

107:   case NTR_INIT_INTERPOLATION:
108:     /*  Use the initial radius specified */
109:     max_radius = 0.0;

111:     for (j = 0; j < j_max; ++j) {
112:       fmin  = f;
113:       sigma = 0.0;

115:       if (needH) {
116:         TaoComputeHessian(tao, tao->solution, tao->hessian, tao->hessian_pre);
117:         needH = 0;
118:       }

120:       for (i = 0; i < i_max; ++i) {
121:         VecCopy(tao->solution, tr->W);
122:         VecAXPY(tr->W, -tao->trust / gnorm, tao->gradient);
123:         TaoComputeObjective(tao, tr->W, &ftrial);

125:         if (PetscIsInfOrNanReal(ftrial)) {
126:           tau = tr->gamma1_i;
127:         } else {
128:           if (ftrial < fmin) {
129:             fmin  = ftrial;
130:             sigma = -tao->trust / gnorm;
131:           }

133:           MatMult(tao->hessian, tao->gradient, tao->stepdirection);
134:           VecDot(tao->gradient, tao->stepdirection, &prered);

136:           prered = tao->trust * (gnorm - 0.5 * tao->trust * prered / (gnorm * gnorm));
137:           actred = f - ftrial;
138:           if ((PetscAbsScalar(actred) <= tr->epsilon) && (PetscAbsScalar(prered) <= tr->epsilon)) {
139:             kappa = 1.0;
140:           } else {
141:             kappa = actred / prered;
142:           }

144:           tau_1   = tr->theta_i * gnorm * tao->trust / (tr->theta_i * gnorm * tao->trust + (1.0 - tr->theta_i) * prered - actred);
145:           tau_2   = tr->theta_i * gnorm * tao->trust / (tr->theta_i * gnorm * tao->trust - (1.0 + tr->theta_i) * prered + actred);
146:           tau_min = PetscMin(tau_1, tau_2);
147:           tau_max = PetscMax(tau_1, tau_2);

149:           if (PetscAbsScalar(kappa - (PetscReal)1.0) <= tr->mu1_i) {
150:             /*  Great agreement */
151:             max_radius = PetscMax(max_radius, tao->trust);

153:             if (tau_max < 1.0) {
154:               tau = tr->gamma3_i;
155:             } else if (tau_max > tr->gamma4_i) {
156:               tau = tr->gamma4_i;
157:             } else {
158:               tau = tau_max;
159:             }
160:           } else if (PetscAbsScalar(kappa - (PetscReal)1.0) <= tr->mu2_i) {
161:             /*  Good agreement */
162:             max_radius = PetscMax(max_radius, tao->trust);

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

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

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

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

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

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

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

231:     if (tr->bfgs_pre) {
232:       /* Update the limited memory preconditioner */
233:       MatLMVMUpdate(tr->M, tao->solution, tao->gradient);
234:     }

236:     while (tao->reason == TAO_CONTINUE_ITERATING) {
237:       KSPSetOperators(tao->ksp, tao->hessian, tao->hessian_pre);

239:       /* Solve the trust region subproblem */
240:       KSPCGSetRadius(tao->ksp, tao->trust);
241:       KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);
242:       KSPGetIterationNumber(tao->ksp, &its);
243:       tao->ksp_its += its;
244:       tao->ksp_tot_its += its;
245:       KSPCGGetNormD(tao->ksp, &norm_d);

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

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

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

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

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

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

297:           if (PetscIsInfOrNanReal(ftrial)) {
298:             tao->trust = tr->alpha1 * PetscMin(tao->trust, norm_d);
299:           } else {
300:             /* Compute and actual reduction */
301:             actred = f - ftrial;
302:             prered = -prered;
303:             if ((PetscAbsScalar(actred) <= tr->epsilon) && (PetscAbsScalar(prered) <= tr->epsilon)) {
304:               kappa = 1.0;
305:             } else {
306:               kappa = actred / prered;
307:             }

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

356:             tau_1   = tr->theta * beta / (tr->theta * beta - (1.0 - tr->theta) * prered + actred);
357:             tau_2   = tr->theta * beta / (tr->theta * beta + (1.0 + tr->theta) * prered - actred);
358:             tau_min = PetscMin(tau_1, tau_2);
359:             tau_max = PetscMax(tau_1, tau_2);

361:             if (kappa >= 1.0 - tr->mu1) {
362:               /* Great agreement; accept step and update radius */
363:               if (tau_max < 1.0) {
364:                 tao->trust = PetscMax(tao->trust, tr->gamma3 * norm_d);
365:               } else if (tau_max > tr->gamma4) {
366:                 tao->trust = PetscMax(tao->trust, tr->gamma4 * norm_d);
367:               } else {
368:                 tao->trust = PetscMax(tao->trust, tau_max * norm_d);
369:               }
370:               break;
371:             } else if (kappa >= 1.0 - tr->mu2) {
372:               /* Good agreement */

374:               if (tau_max < tr->gamma2) {
375:                 tao->trust = tr->gamma2 * PetscMin(tao->trust, norm_d);
376:               } else if (tau_max > tr->gamma3) {
377:                 tao->trust = PetscMax(tao->trust, tr->gamma3 * norm_d);
378:               } else if (tau_max < 1.0) {
379:                 tao->trust = tau_max * PetscMin(tao->trust, norm_d);
380:               } else {
381:                 tao->trust = PetscMax(tao->trust, tau_max * norm_d);
382:               }
383:               break;
384:             } else {
385:               /* Not good agreement */
386:               if (tau_min > 1.0) {
387:                 tao->trust = tr->gamma2 * PetscMin(tao->trust, norm_d);
388:               } else if (tau_max < tr->gamma1) {
389:                 tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d);
390:               } else if ((tau_min < tr->gamma1) && (tau_max >= 1.0)) {
391:                 tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d);
392:               } else if ((tau_1 >= tr->gamma1) && (tau_1 < 1.0) && ((tau_2 < tr->gamma1) || (tau_2 >= 1.0))) {
393:                 tao->trust = tau_1 * PetscMin(tao->trust, norm_d);
394:               } else if ((tau_2 >= tr->gamma1) && (tau_2 < 1.0) && ((tau_1 < tr->gamma1) || (tau_2 >= 1.0))) {
395:                 tao->trust = tau_2 * PetscMin(tao->trust, norm_d);
396:               } else {
397:                 tao->trust = tau_max * PetscMin(tao->trust, norm_d);
398:               }
399:             }
400:           }
401:         }
402:       }

404:       /* The step computed was not good and the radius was decreased.
405:          Monitor the radius to terminate. */
406:       TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its);
407:       TaoMonitor(tao, tao->niter, f, gnorm, 0.0, tao->trust);
408:       PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
409:     }

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

414:     if (tao->reason == TAO_CONTINUE_ITERATING) {
415:       VecCopy(tr->W, tao->solution);
416:       f = ftrial;
417:       TaoComputeGradient(tao, tao->solution, tao->gradient);
418:       TaoGradientNorm(tao, tao->gradient, NORM_2, &gnorm);
420:       needH = 1;
421:       TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its);
422:       TaoMonitor(tao, tao->niter, f, gnorm, 0.0, tao->trust);
423:       PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
424:     }
425:   }
426:   return 0;
427: }

429: /*------------------------------------------------------------*/
430: static PetscErrorCode TaoSetUp_NTR(Tao tao)
431: {
432:   TAO_NTR *tr = (TAO_NTR *)tao->data;

434:   if (!tao->gradient) VecDuplicate(tao->solution, &tao->gradient);
435:   if (!tao->stepdirection) VecDuplicate(tao->solution, &tao->stepdirection);
436:   if (!tr->W) VecDuplicate(tao->solution, &tr->W);

438:   tr->bfgs_pre = NULL;
439:   tr->M        = NULL;
440:   return 0;
441: }

443: /*------------------------------------------------------------*/
444: static PetscErrorCode TaoDestroy_NTR(Tao tao)
445: {
446:   TAO_NTR *tr = (TAO_NTR *)tao->data;

448:   if (tao->setupcalled) VecDestroy(&tr->W);
449:   KSPDestroy(&tao->ksp);
450:   PetscFree(tao->data);
451:   return 0;
452: }

454: /*------------------------------------------------------------*/
455: static PetscErrorCode TaoSetFromOptions_NTR(Tao tao, PetscOptionItems *PetscOptionsObject)
456: {
457:   TAO_NTR *tr = (TAO_NTR *)tao->data;

459:   PetscOptionsHeadBegin(PetscOptionsObject, "Newton trust region method for unconstrained optimization");
460:   PetscOptionsEList("-tao_ntr_init_type", "tao->trust initialization type", "", NTR_INIT, NTR_INIT_TYPES, NTR_INIT[tr->init_type], &tr->init_type, NULL);
461:   PetscOptionsEList("-tao_ntr_update_type", "radius update type", "", NTR_UPDATE, NTR_UPDATE_TYPES, NTR_UPDATE[tr->update_type], &tr->update_type, NULL);
462:   PetscOptionsReal("-tao_ntr_eta1", "step is unsuccessful if actual reduction < eta1 * predicted reduction", "", tr->eta1, &tr->eta1, NULL);
463:   PetscOptionsReal("-tao_ntr_eta2", "", "", tr->eta2, &tr->eta2, NULL);
464:   PetscOptionsReal("-tao_ntr_eta3", "", "", tr->eta3, &tr->eta3, NULL);
465:   PetscOptionsReal("-tao_ntr_eta4", "", "", tr->eta4, &tr->eta4, NULL);
466:   PetscOptionsReal("-tao_ntr_alpha1", "", "", tr->alpha1, &tr->alpha1, NULL);
467:   PetscOptionsReal("-tao_ntr_alpha2", "", "", tr->alpha2, &tr->alpha2, NULL);
468:   PetscOptionsReal("-tao_ntr_alpha3", "", "", tr->alpha3, &tr->alpha3, NULL);
469:   PetscOptionsReal("-tao_ntr_alpha4", "", "", tr->alpha4, &tr->alpha4, NULL);
470:   PetscOptionsReal("-tao_ntr_alpha5", "", "", tr->alpha5, &tr->alpha5, NULL);
471:   PetscOptionsReal("-tao_ntr_mu1", "", "", tr->mu1, &tr->mu1, NULL);
472:   PetscOptionsReal("-tao_ntr_mu2", "", "", tr->mu2, &tr->mu2, NULL);
473:   PetscOptionsReal("-tao_ntr_gamma1", "", "", tr->gamma1, &tr->gamma1, NULL);
474:   PetscOptionsReal("-tao_ntr_gamma2", "", "", tr->gamma2, &tr->gamma2, NULL);
475:   PetscOptionsReal("-tao_ntr_gamma3", "", "", tr->gamma3, &tr->gamma3, NULL);
476:   PetscOptionsReal("-tao_ntr_gamma4", "", "", tr->gamma4, &tr->gamma4, NULL);
477:   PetscOptionsReal("-tao_ntr_theta", "", "", tr->theta, &tr->theta, NULL);
478:   PetscOptionsReal("-tao_ntr_mu1_i", "", "", tr->mu1_i, &tr->mu1_i, NULL);
479:   PetscOptionsReal("-tao_ntr_mu2_i", "", "", tr->mu2_i, &tr->mu2_i, NULL);
480:   PetscOptionsReal("-tao_ntr_gamma1_i", "", "", tr->gamma1_i, &tr->gamma1_i, NULL);
481:   PetscOptionsReal("-tao_ntr_gamma2_i", "", "", tr->gamma2_i, &tr->gamma2_i, NULL);
482:   PetscOptionsReal("-tao_ntr_gamma3_i", "", "", tr->gamma3_i, &tr->gamma3_i, NULL);
483:   PetscOptionsReal("-tao_ntr_gamma4_i", "", "", tr->gamma4_i, &tr->gamma4_i, NULL);
484:   PetscOptionsReal("-tao_ntr_theta_i", "", "", tr->theta_i, &tr->theta_i, NULL);
485:   PetscOptionsReal("-tao_ntr_min_radius", "lower bound on initial trust-region radius", "", tr->min_radius, &tr->min_radius, NULL);
486:   PetscOptionsReal("-tao_ntr_max_radius", "upper bound on trust-region radius", "", tr->max_radius, &tr->max_radius, NULL);
487:   PetscOptionsReal("-tao_ntr_epsilon", "tolerance used when computing actual and predicted reduction", "", tr->epsilon, &tr->epsilon, NULL);
488:   PetscOptionsHeadEnd();
489:   KSPSetFromOptions(tao->ksp);
490:   return 0;
491: }

493: /*------------------------------------------------------------*/
494: /*MC
495:   TAONTR - Newton's method with trust region for unconstrained minimization.
496:   At each iteration, the Newton trust region method solves the system.
497:   NTR expects a KSP solver with a trust region radius.
498:             min_d  .5 dT Hk d + gkT d,  s.t.   ||d|| < Delta_k

500:   Options Database Keys:
501: + -tao_ntr_init_type - "constant","direction","interpolation"
502: . -tao_ntr_update_type - "reduction","interpolation"
503: . -tao_ntr_min_radius - lower bound on trust region radius
504: . -tao_ntr_max_radius - upper bound on trust region radius
505: . -tao_ntr_epsilon - tolerance for accepting actual / predicted reduction
506: . -tao_ntr_mu1_i - mu1 interpolation init factor
507: . -tao_ntr_mu2_i - mu2 interpolation init factor
508: . -tao_ntr_gamma1_i - gamma1 interpolation init factor
509: . -tao_ntr_gamma2_i - gamma2 interpolation init factor
510: . -tao_ntr_gamma3_i - gamma3 interpolation init factor
511: . -tao_ntr_gamma4_i - gamma4 interpolation init factor
512: . -tao_ntr_theta_i - theta1 interpolation init factor
513: . -tao_ntr_eta1 - eta1 reduction update factor
514: . -tao_ntr_eta2 - eta2 reduction update factor
515: . -tao_ntr_eta3 - eta3 reduction update factor
516: . -tao_ntr_eta4 - eta4 reduction update factor
517: . -tao_ntr_alpha1 - alpha1 reduction update factor
518: . -tao_ntr_alpha2 - alpha2 reduction update factor
519: . -tao_ntr_alpha3 - alpha3 reduction update factor
520: . -tao_ntr_alpha4 - alpha4 reduction update factor
521: . -tao_ntr_alpha4 - alpha4 reduction update factor
522: . -tao_ntr_mu1 - mu1 interpolation update
523: . -tao_ntr_mu2 - mu2 interpolation update
524: . -tao_ntr_gamma1 - gamma1 interpolcation update
525: . -tao_ntr_gamma2 - gamma2 interpolcation update
526: . -tao_ntr_gamma3 - gamma3 interpolcation update
527: . -tao_ntr_gamma4 - gamma4 interpolation update
528: - -tao_ntr_theta - theta interpolation update

530:   Level: beginner
531: M*/

533: PETSC_EXTERN PetscErrorCode TaoCreate_NTR(Tao tao)
534: {
535:   TAO_NTR *tr;


538:   PetscNew(&tr);

540:   tao->ops->setup          = TaoSetUp_NTR;
541:   tao->ops->solve          = TaoSolve_NTR;
542:   tao->ops->setfromoptions = TaoSetFromOptions_NTR;
543:   tao->ops->destroy        = TaoDestroy_NTR;

545:   /* Override default settings (unless already changed) */
546:   if (!tao->max_it_changed) tao->max_it = 50;
547:   if (!tao->trust0_changed) tao->trust0 = 100.0;
548:   tao->data = (void *)tr;

550:   /*  Standard trust region update parameters */
551:   tr->eta1 = 1.0e-4;
552:   tr->eta2 = 0.25;
553:   tr->eta3 = 0.50;
554:   tr->eta4 = 0.90;

556:   tr->alpha1 = 0.25;
557:   tr->alpha2 = 0.50;
558:   tr->alpha3 = 1.00;
559:   tr->alpha4 = 2.00;
560:   tr->alpha5 = 4.00;

562:   /*  Interpolation trust region update parameters */
563:   tr->mu1 = 0.10;
564:   tr->mu2 = 0.50;

566:   tr->gamma1 = 0.25;
567:   tr->gamma2 = 0.50;
568:   tr->gamma3 = 2.00;
569:   tr->gamma4 = 4.00;

571:   tr->theta = 0.05;

573:   /*  Interpolation parameters for initialization */
574:   tr->mu1_i = 0.35;
575:   tr->mu2_i = 0.50;

577:   tr->gamma1_i = 0.0625;
578:   tr->gamma2_i = 0.50;
579:   tr->gamma3_i = 2.00;
580:   tr->gamma4_i = 5.00;

582:   tr->theta_i = 0.25;

584:   tr->min_radius = 1.0e-10;
585:   tr->max_radius = 1.0e10;
586:   tr->epsilon    = 1.0e-6;

588:   tr->init_type   = NTR_INIT_INTERPOLATION;
589:   tr->update_type = NTR_UPDATE_REDUCTION;

591:   /* Set linear solver to default for trust region */
592:   KSPCreate(((PetscObject)tao)->comm, &tao->ksp);
593:   PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1);
594:   KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix);
595:   KSPAppendOptionsPrefix(tao->ksp, "tao_ntr_");
596:   KSPSetType(tao->ksp, KSPSTCG);
597:   return 0;
598: }