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           bfgsUpdates = 0;
 56:   PetscInt           needH;

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

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

 66:   KSPGetType(tao->ksp,&ksp_type);
 67:   PetscStrcmp(ksp_type,KSPNASH,&is_nash);
 68:   PetscStrcmp(ksp_type,KSPSTCG,&is_stcg);
 69:   PetscStrcmp(ksp_type,KSPGLTR,&is_gltr);

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

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

 94:   /* Check convergence criteria */
 95:   TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);
 96:   TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm);
 98:   needH = 1;

100:   tao->reason = TAO_CONTINUE_ITERATING;
101:   TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);
102:   TaoMonitor(tao,tao->niter,f,gnorm,0.0,1.0);
103:   (*tao->ops->convergencetest)(tao,tao->cnvP);
104:   if (tao->reason != TAO_CONTINUE_ITERATING) return 0;

106:   /*  Initialize trust-region radius */
107:   switch (tr->init_type) {
108:   case NTR_INIT_CONSTANT:
109:     /*  Use the initial radius specified */
110:     break;

112:   case NTR_INIT_INTERPOLATION:
113:     /*  Use the initial radius specified */
114:     max_radius = 0.0;

116:     for (j = 0; j < j_max; ++j) {
117:       fmin = f;
118:       sigma = 0.0;

120:       if (needH) {
121:         TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);
122:         needH = 0;
123:       }

125:       for (i = 0; i < i_max; ++i) {

127:         VecCopy(tao->solution, tr->W);
128:         VecAXPY(tr->W, -tao->trust/gnorm, tao->gradient);
129:         TaoComputeObjective(tao, tr->W, &ftrial);

131:         if (PetscIsInfOrNanReal(ftrial)) {
132:           tau = tr->gamma1_i;
133:         }
134:         else {
135:           if (ftrial < fmin) {
136:             fmin = ftrial;
137:             sigma = -tao->trust / gnorm;
138:           }

140:           MatMult(tao->hessian, tao->gradient, tao->stepdirection);
141:           VecDot(tao->gradient, tao->stepdirection, &prered);

143:           prered = tao->trust * (gnorm - 0.5 * tao->trust * prered / (gnorm * gnorm));
144:           actred = f - ftrial;
145:           if ((PetscAbsScalar(actred) <= tr->epsilon) &&
146:               (PetscAbsScalar(prered) <= tr->epsilon)) {
147:             kappa = 1.0;
148:           }
149:           else {
150:             kappa = actred / prered;
151:           }

153:           tau_1 = tr->theta_i * gnorm * tao->trust / (tr->theta_i * gnorm * tao->trust + (1.0 - tr->theta_i) * prered - actred);
154:           tau_2 = tr->theta_i * gnorm * tao->trust / (tr->theta_i * gnorm * tao->trust - (1.0 + tr->theta_i) * prered + actred);
155:           tau_min = PetscMin(tau_1, tau_2);
156:           tau_max = PetscMax(tau_1, tau_2);

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

162:             if (tau_max < 1.0) {
163:               tau = tr->gamma3_i;
164:             }
165:             else if (tau_max > tr->gamma4_i) {
166:               tau = tr->gamma4_i;
167:             }
168:             else {
169:               tau = tau_max;
170:             }
171:           }
172:           else if (PetscAbsScalar(kappa - (PetscReal)1.0) <= tr->mu2_i) {
173:             /*  Good agreement */
174:             max_radius = PetscMax(max_radius, tao->trust);

176:             if (tau_max < tr->gamma2_i) {
177:               tau = tr->gamma2_i;
178:             }
179:             else if (tau_max > tr->gamma3_i) {
180:               tau = tr->gamma3_i;
181:             }
182:             else {
183:               tau = tau_max;
184:             }
185:           }
186:           else {
187:             /*  Not good agreement */
188:             if (tau_min > 1.0) {
189:               tau = tr->gamma2_i;
190:             }
191:             else if (tau_max < tr->gamma1_i) {
192:               tau = tr->gamma1_i;
193:             }
194:             else if ((tau_min < tr->gamma1_i) && (tau_max >= 1.0)) {
195:               tau = tr->gamma1_i;
196:             }
197:             else if ((tau_1 >= tr->gamma1_i) && (tau_1 < 1.0) &&
198:                      ((tau_2 < tr->gamma1_i) || (tau_2 >= 1.0))) {
199:               tau = tau_1;
200:             }
201:             else if ((tau_2 >= tr->gamma1_i) && (tau_2 < 1.0) &&
202:                      ((tau_1 < tr->gamma1_i) || (tau_2 >= 1.0))) {
203:               tau = tau_2;
204:             }
205:             else {
206:               tau = tau_max;
207:             }
208:           }
209:         }
210:         tao->trust = tau * tao->trust;
211:       }

213:       if (fmin < f) {
214:         f = fmin;
215:         VecAXPY(tao->solution, sigma, tao->gradient);
216:         TaoComputeGradient(tao,tao->solution, tao->gradient);

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

222:         TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);
223:         TaoMonitor(tao,tao->niter,f,gnorm,0.0,1.0);
224:         (*tao->ops->convergencetest)(tao,tao->cnvP);
225:         if (tao->reason != TAO_CONTINUE_ITERATING) {
226:           return 0;
227:         }
228:       }
229:     }
230:     tao->trust = PetscMax(tao->trust, max_radius);

232:     /*  Modify the radius if it is too large or small */
233:     tao->trust = PetscMax(tao->trust, tr->min_radius);
234:     tao->trust = PetscMin(tao->trust, tr->max_radius);
235:     break;

237:   default:
238:     /*  Norm of the first direction will initialize radius */
239:     tao->trust = 0.0;
240:     break;
241:   }

243:   /* Have not converged; continue with Newton method */
244:   while (tao->reason == TAO_CONTINUE_ITERATING) {
245:     /* Call general purpose update function */
246:     if (tao->ops->update) {
247:       (*tao->ops->update)(tao, tao->niter, tao->user_update);
248:     }
249:     ++tao->niter;
250:     tao->ksp_its=0;
251:     /* Compute the Hessian */
252:     if (needH) {
253:       TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);
254:       needH = 0;
255:     }

257:     if (tr->bfgs_pre) {
258:       /* Update the limited memory preconditioner */
259:       MatLMVMUpdate(tr->M, tao->solution, tao->gradient);
260:       ++bfgsUpdates;
261:     }

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

266:       /* Solve the trust region subproblem */
267:       KSPCGSetRadius(tao->ksp,tao->trust);
268:       KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);
269:       KSPGetIterationNumber(tao->ksp,&its);
270:       tao->ksp_its+=its;
271:       tao->ksp_tot_its+=its;
272:       KSPCGGetNormD(tao->ksp, &norm_d);

274:       if (0.0 == tao->trust) {
275:         /* Radius was uninitialized; use the norm of the direction */
276:         if (norm_d > 0.0) {
277:           tao->trust = norm_d;

279:           /* Modify the radius if it is too large or small */
280:           tao->trust = PetscMax(tao->trust, tr->min_radius);
281:           tao->trust = PetscMin(tao->trust, tr->max_radius);
282:         }
283:         else {
284:           /* The direction was bad; set radius to default value and re-solve
285:              the trust-region subproblem to get a direction */
286:           tao->trust = tao->trust0;

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

292:           KSPCGSetRadius(tao->ksp,tao->trust);
293:           KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);
294:           KSPGetIterationNumber(tao->ksp,&its);
295:           tao->ksp_its+=its;
296:           tao->ksp_tot_its+=its;
297:           KSPCGGetNormD(tao->ksp, &norm_d);

300:         }
301:       }
302:       VecScale(tao->stepdirection, -1.0);
303:       KSPGetConvergedReason(tao->ksp, &ksp_reason);
304:       if ((KSP_DIVERGED_INDEFINITE_PC == ksp_reason) && (tr->bfgs_pre)) {
305:         /* Preconditioner is numerically indefinite; reset the
306:            approximate if using BFGS preconditioning. */
307:         MatLMVMReset(tr->M, PETSC_FALSE);
308:         MatLMVMUpdate(tr->M, tao->solution, tao->gradient);
309:         bfgsUpdates = 1;
310:       }

312:       if (NTR_UPDATE_REDUCTION == tr->update_type) {
313:         /* Get predicted reduction */
314:         KSPCGGetObjFcn(tao->ksp,&prered);
315:         if (prered >= 0.0) {
316:           /* The predicted reduction has the wrong sign.  This cannot
317:              happen in infinite precision arithmetic.  Step should
318:              be rejected! */
319:           tao->trust = tr->alpha1 * PetscMin(tao->trust, norm_d);
320:         }
321:         else {
322:           /* Compute trial step and function value */
323:           VecCopy(tao->solution,tr->W);
324:           VecAXPY(tr->W, 1.0, tao->stepdirection);
325:           TaoComputeObjective(tao, tr->W, &ftrial);

327:           if (PetscIsInfOrNanReal(ftrial)) {
328:             tao->trust = tr->alpha1 * PetscMin(tao->trust, norm_d);
329:           } else {
330:             /* Compute and actual reduction */
331:             actred = f - ftrial;
332:             prered = -prered;
333:             if ((PetscAbsScalar(actred) <= tr->epsilon) &&
334:                 (PetscAbsScalar(prered) <= tr->epsilon)) {
335:               kappa = 1.0;
336:             }
337:             else {
338:               kappa = actred / prered;
339:             }

341:             /* Accept or reject the step and update radius */
342:             if (kappa < tr->eta1) {
343:               /* Reject the step */
344:               tao->trust = tr->alpha1 * PetscMin(tao->trust, norm_d);
345:             }
346:             else {
347:               /* Accept the step */
348:               if (kappa < tr->eta2) {
349:                 /* Marginal bad step */
350:                 tao->trust = tr->alpha2 * PetscMin(tao->trust, norm_d);
351:               }
352:               else if (kappa < tr->eta3) {
353:                 /* Reasonable step */
354:                 tao->trust = tr->alpha3 * tao->trust;
355:               }
356:               else if (kappa < tr->eta4) {
357:                 /* Good step */
358:                 tao->trust = PetscMax(tr->alpha4 * norm_d, tao->trust);
359:               }
360:               else {
361:                 /* Very good step */
362:                 tao->trust = PetscMax(tr->alpha5 * norm_d, tao->trust);
363:               }
364:               break;
365:             }
366:           }
367:         }
368:       }
369:       else {
370:         /* Get predicted reduction */
371:         KSPCGGetObjFcn(tao->ksp,&prered);
372:         if (prered >= 0.0) {
373:           /* The predicted reduction has the wrong sign.  This cannot
374:              happen in infinite precision arithmetic.  Step should
375:              be rejected! */
376:           tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d);
377:         }
378:         else {
379:           VecCopy(tao->solution, tr->W);
380:           VecAXPY(tr->W, 1.0, tao->stepdirection);
381:           TaoComputeObjective(tao, tr->W, &ftrial);
382:           if (PetscIsInfOrNanReal(ftrial)) {
383:             tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d);
384:           }
385:           else {
386:             VecDot(tao->gradient, tao->stepdirection, &beta);
387:             actred = f - ftrial;
388:             prered = -prered;
389:             if ((PetscAbsScalar(actred) <= tr->epsilon) &&
390:                 (PetscAbsScalar(prered) <= tr->epsilon)) {
391:               kappa = 1.0;
392:             }
393:             else {
394:               kappa = actred / prered;
395:             }

397:             tau_1 = tr->theta * beta / (tr->theta * beta - (1.0 - tr->theta) * prered + actred);
398:             tau_2 = tr->theta * beta / (tr->theta * beta + (1.0 + tr->theta) * prered - actred);
399:             tau_min = PetscMin(tau_1, tau_2);
400:             tau_max = PetscMax(tau_1, tau_2);

402:             if (kappa >= 1.0 - tr->mu1) {
403:               /* Great agreement; accept step and update radius */
404:               if (tau_max < 1.0) {
405:                 tao->trust = PetscMax(tao->trust, tr->gamma3 * norm_d);
406:               }
407:               else if (tau_max > tr->gamma4) {
408:                 tao->trust = PetscMax(tao->trust, tr->gamma4 * norm_d);
409:               }
410:               else {
411:                 tao->trust = PetscMax(tao->trust, tau_max * norm_d);
412:               }
413:               break;
414:             }
415:             else if (kappa >= 1.0 - tr->mu2) {
416:               /* Good agreement */

418:               if (tau_max < tr->gamma2) {
419:                 tao->trust = tr->gamma2 * PetscMin(tao->trust, norm_d);
420:               }
421:               else if (tau_max > tr->gamma3) {
422:                 tao->trust = PetscMax(tao->trust, tr->gamma3 * norm_d);
423:               }
424:               else if (tau_max < 1.0) {
425:                 tao->trust = tau_max * PetscMin(tao->trust, norm_d);
426:               }
427:               else {
428:                 tao->trust = PetscMax(tao->trust, tau_max * norm_d);
429:               }
430:               break;
431:             }
432:             else {
433:               /* Not good agreement */
434:               if (tau_min > 1.0) {
435:                 tao->trust = tr->gamma2 * PetscMin(tao->trust, norm_d);
436:               }
437:               else if (tau_max < tr->gamma1) {
438:                 tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d);
439:               }
440:               else if ((tau_min < tr->gamma1) && (tau_max >= 1.0)) {
441:                 tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d);
442:               }
443:               else if ((tau_1 >= tr->gamma1) && (tau_1 < 1.0) &&
444:                        ((tau_2 < tr->gamma1) || (tau_2 >= 1.0))) {
445:                 tao->trust = tau_1 * PetscMin(tao->trust, norm_d);
446:               }
447:               else if ((tau_2 >= tr->gamma1) && (tau_2 < 1.0) &&
448:                        ((tau_1 < tr->gamma1) || (tau_2 >= 1.0))) {
449:                 tao->trust = tau_2 * PetscMin(tao->trust, norm_d);
450:               }
451:               else {
452:                 tao->trust = tau_max * PetscMin(tao->trust, norm_d);
453:               }
454:             }
455:           }
456:         }
457:       }

459:       /* The step computed was not good and the radius was decreased.
460:          Monitor the radius to terminate. */
461:       TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);
462:       TaoMonitor(tao,tao->niter,f,gnorm,0.0,tao->trust);
463:       (*tao->ops->convergencetest)(tao,tao->cnvP);
464:     }

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

469:     if (tao->reason == TAO_CONTINUE_ITERATING) {
470:       VecCopy(tr->W, tao->solution);
471:       f = ftrial;
472:       TaoComputeGradient(tao, tao->solution, tao->gradient);
473:       TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm);
475:       needH = 1;
476:       TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);
477:       TaoMonitor(tao,tao->niter,f,gnorm,0.0,tao->trust);
478:       (*tao->ops->convergencetest)(tao,tao->cnvP);
479:     }
480:   }
481:   return 0;
482: }

484: /*------------------------------------------------------------*/
485: static PetscErrorCode TaoSetUp_NTR(Tao tao)
486: {
487:   TAO_NTR *tr = (TAO_NTR *)tao->data;

489:   if (!tao->gradient) VecDuplicate(tao->solution, &tao->gradient);
490:   if (!tao->stepdirection) VecDuplicate(tao->solution, &tao->stepdirection);
491:   if (!tr->W) VecDuplicate(tao->solution, &tr->W);

493:   tr->bfgs_pre = NULL;
494:   tr->M = NULL;
495:   return 0;
496: }

498: /*------------------------------------------------------------*/
499: static PetscErrorCode TaoDestroy_NTR(Tao tao)
500: {
501:   TAO_NTR        *tr = (TAO_NTR *)tao->data;

503:   if (tao->setupcalled) {
504:     VecDestroy(&tr->W);
505:   }
506:   PetscFree(tao->data);
507:   return 0;
508: }

510: /*------------------------------------------------------------*/
511: static PetscErrorCode TaoSetFromOptions_NTR(PetscOptionItems *PetscOptionsObject,Tao tao)
512: {
513:   TAO_NTR        *tr = (TAO_NTR *)tao->data;

515:   PetscOptionsHead(PetscOptionsObject,"Newton trust region method for unconstrained optimization");
516:   PetscOptionsEList("-tao_ntr_init_type", "tao->trust initialization type", "", NTR_INIT, NTR_INIT_TYPES, NTR_INIT[tr->init_type], &tr->init_type,NULL);
517:   PetscOptionsEList("-tao_ntr_update_type", "radius update type", "", NTR_UPDATE, NTR_UPDATE_TYPES, NTR_UPDATE[tr->update_type], &tr->update_type,NULL);
518:   PetscOptionsReal("-tao_ntr_eta1", "step is unsuccessful if actual reduction < eta1 * predicted reduction", "", tr->eta1, &tr->eta1,NULL);
519:   PetscOptionsReal("-tao_ntr_eta2", "", "", tr->eta2, &tr->eta2,NULL);
520:   PetscOptionsReal("-tao_ntr_eta3", "", "", tr->eta3, &tr->eta3,NULL);
521:   PetscOptionsReal("-tao_ntr_eta4", "", "", tr->eta4, &tr->eta4,NULL);
522:   PetscOptionsReal("-tao_ntr_alpha1", "", "", tr->alpha1, &tr->alpha1,NULL);
523:   PetscOptionsReal("-tao_ntr_alpha2", "", "", tr->alpha2, &tr->alpha2,NULL);
524:   PetscOptionsReal("-tao_ntr_alpha3", "", "", tr->alpha3, &tr->alpha3,NULL);
525:   PetscOptionsReal("-tao_ntr_alpha4", "", "", tr->alpha4, &tr->alpha4,NULL);
526:   PetscOptionsReal("-tao_ntr_alpha5", "", "", tr->alpha5, &tr->alpha5,NULL);
527:   PetscOptionsReal("-tao_ntr_mu1", "", "", tr->mu1, &tr->mu1,NULL);
528:   PetscOptionsReal("-tao_ntr_mu2", "", "", tr->mu2, &tr->mu2,NULL);
529:   PetscOptionsReal("-tao_ntr_gamma1", "", "", tr->gamma1, &tr->gamma1,NULL);
530:   PetscOptionsReal("-tao_ntr_gamma2", "", "", tr->gamma2, &tr->gamma2,NULL);
531:   PetscOptionsReal("-tao_ntr_gamma3", "", "", tr->gamma3, &tr->gamma3,NULL);
532:   PetscOptionsReal("-tao_ntr_gamma4", "", "", tr->gamma4, &tr->gamma4,NULL);
533:   PetscOptionsReal("-tao_ntr_theta", "", "", tr->theta, &tr->theta,NULL);
534:   PetscOptionsReal("-tao_ntr_mu1_i", "", "", tr->mu1_i, &tr->mu1_i,NULL);
535:   PetscOptionsReal("-tao_ntr_mu2_i", "", "", tr->mu2_i, &tr->mu2_i,NULL);
536:   PetscOptionsReal("-tao_ntr_gamma1_i", "", "", tr->gamma1_i, &tr->gamma1_i,NULL);
537:   PetscOptionsReal("-tao_ntr_gamma2_i", "", "", tr->gamma2_i, &tr->gamma2_i,NULL);
538:   PetscOptionsReal("-tao_ntr_gamma3_i", "", "", tr->gamma3_i, &tr->gamma3_i,NULL);
539:   PetscOptionsReal("-tao_ntr_gamma4_i", "", "", tr->gamma4_i, &tr->gamma4_i,NULL);
540:   PetscOptionsReal("-tao_ntr_theta_i", "", "", tr->theta_i, &tr->theta_i,NULL);
541:   PetscOptionsReal("-tao_ntr_min_radius", "lower bound on initial trust-region radius", "", tr->min_radius, &tr->min_radius,NULL);
542:   PetscOptionsReal("-tao_ntr_max_radius", "upper bound on trust-region radius", "", tr->max_radius, &tr->max_radius,NULL);
543:   PetscOptionsReal("-tao_ntr_epsilon", "tolerance used when computing actual and predicted reduction", "", tr->epsilon, &tr->epsilon,NULL);
544:   PetscOptionsTail();
545:   KSPSetFromOptions(tao->ksp);
546:   return 0;
547: }

549: /*------------------------------------------------------------*/
550: /*MC
551:   TAONTR - Newton's method with trust region for unconstrained minimization.
552:   At each iteration, the Newton trust region method solves the system.
553:   NTR expects a KSP solver with a trust region radius.
554:             min_d  .5 dT Hk d + gkT d,  s.t.   ||d|| < Delta_k

556:   Options Database Keys:
557: + -tao_ntr_init_type - "constant","direction","interpolation"
558: . -tao_ntr_update_type - "reduction","interpolation"
559: . -tao_ntr_min_radius - lower bound on trust region radius
560: . -tao_ntr_max_radius - upper bound on trust region radius
561: . -tao_ntr_epsilon - tolerance for accepting actual / predicted reduction
562: . -tao_ntr_mu1_i - mu1 interpolation init factor
563: . -tao_ntr_mu2_i - mu2 interpolation init factor
564: . -tao_ntr_gamma1_i - gamma1 interpolation init factor
565: . -tao_ntr_gamma2_i - gamma2 interpolation init factor
566: . -tao_ntr_gamma3_i - gamma3 interpolation init factor
567: . -tao_ntr_gamma4_i - gamma4 interpolation init factor
568: . -tao_ntr_theta_i - theta1 interpolation init factor
569: . -tao_ntr_eta1 - eta1 reduction update factor
570: . -tao_ntr_eta2 - eta2 reduction update factor
571: . -tao_ntr_eta3 - eta3 reduction update factor
572: . -tao_ntr_eta4 - eta4 reduction update factor
573: . -tao_ntr_alpha1 - alpha1 reduction update factor
574: . -tao_ntr_alpha2 - alpha2 reduction update factor
575: . -tao_ntr_alpha3 - alpha3 reduction update factor
576: . -tao_ntr_alpha4 - alpha4 reduction update factor
577: . -tao_ntr_alpha4 - alpha4 reduction update factor
578: . -tao_ntr_mu1 - mu1 interpolation update
579: . -tao_ntr_mu2 - mu2 interpolation update
580: . -tao_ntr_gamma1 - gamma1 interpolcation update
581: . -tao_ntr_gamma2 - gamma2 interpolcation update
582: . -tao_ntr_gamma3 - gamma3 interpolcation update
583: . -tao_ntr_gamma4 - gamma4 interpolation update
584: - -tao_ntr_theta - theta interpolation update

586:   Level: beginner
587: M*/

589: PETSC_EXTERN PetscErrorCode TaoCreate_NTR(Tao tao)
590: {
591:   TAO_NTR *tr;


594:   PetscNewLog(tao,&tr);

596:   tao->ops->setup = TaoSetUp_NTR;
597:   tao->ops->solve = TaoSolve_NTR;
598:   tao->ops->setfromoptions = TaoSetFromOptions_NTR;
599:   tao->ops->destroy = TaoDestroy_NTR;

601:   /* Override default settings (unless already changed) */
602:   if (!tao->max_it_changed) tao->max_it = 50;
603:   if (!tao->trust0_changed) tao->trust0 = 100.0;
604:   tao->data = (void*)tr;

606:   /*  Standard trust region update parameters */
607:   tr->eta1 = 1.0e-4;
608:   tr->eta2 = 0.25;
609:   tr->eta3 = 0.50;
610:   tr->eta4 = 0.90;

612:   tr->alpha1 = 0.25;
613:   tr->alpha2 = 0.50;
614:   tr->alpha3 = 1.00;
615:   tr->alpha4 = 2.00;
616:   tr->alpha5 = 4.00;

618:   /*  Interpolation trust region update parameters */
619:   tr->mu1 = 0.10;
620:   tr->mu2 = 0.50;

622:   tr->gamma1 = 0.25;
623:   tr->gamma2 = 0.50;
624:   tr->gamma3 = 2.00;
625:   tr->gamma4 = 4.00;

627:   tr->theta = 0.05;

629:   /*  Interpolation parameters for initialization */
630:   tr->mu1_i = 0.35;
631:   tr->mu2_i = 0.50;

633:   tr->gamma1_i = 0.0625;
634:   tr->gamma2_i = 0.50;
635:   tr->gamma3_i = 2.00;
636:   tr->gamma4_i = 5.00;

638:   tr->theta_i = 0.25;

640:   tr->min_radius = 1.0e-10;
641:   tr->max_radius = 1.0e10;
642:   tr->epsilon    = 1.0e-6;

644:   tr->init_type       = NTR_INIT_INTERPOLATION;
645:   tr->update_type     = NTR_UPDATE_REDUCTION;

647:   /* Set linear solver to default for trust region */
648:   KSPCreate(((PetscObject)tao)->comm,&tao->ksp);
649:   PetscObjectIncrementTabLevel((PetscObject)tao->ksp,(PetscObject)tao,1);
650:   KSPSetOptionsPrefix(tao->ksp,tao->hdr.prefix);
651:   KSPAppendOptionsPrefix(tao->ksp,"tao_ntr_");
652:   KSPSetType(tao->ksp,KSPSTCG);
653:   return 0;
654: }