Actual source code: ntr.c

petsc-3.12.5 2020-03-29
Report Typos and Errors
  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:   PetscErrorCode     ierr;
 56:   PetscInt           bfgsUpdates = 0;
 57:   PetscInt           needH;

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

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

 68:   KSPGetType(tao->ksp,&ksp_type);
 69:   PetscStrcmp(ksp_type,KSPNASH,&is_nash);
 70:   PetscStrcmp(ksp_type,KSPSTCG,&is_stcg);
 71:   PetscStrcmp(ksp_type,KSPGLTR,&is_gltr);
 72:   if (!is_nash && !is_stcg && !is_gltr) {
 73:     SETERRQ(PETSC_COMM_SELF,1,"TAO_NTR requires nash, stcg, or gltr for the KSP");
 74:   }

 76:   /* Initialize the radius and modify if it is too large or small */
 77:   tao->trust = tao->trust0;
 78:   tao->trust = PetscMax(tao->trust, tr->min_radius);
 79:   tao->trust = PetscMin(tao->trust, tr->max_radius);

 81:   /* Allocate the vectors needed for the BFGS approximation */
 82:   KSPGetPC(tao->ksp, &pc);
 83:   PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &is_bfgs);
 84:   PetscObjectTypeCompare((PetscObject)pc, PCJACOBI, &is_jacobi);
 85:   if (is_bfgs) {
 86:     tr->bfgs_pre = pc;
 87:     PCLMVMGetMatLMVM(tr->bfgs_pre, &tr->M);
 88:     VecGetLocalSize(tao->solution, &n);
 89:     VecGetSize(tao->solution, &N);
 90:     MatSetSizes(tr->M, n, n, N, N);
 91:     MatLMVMAllocate(tr->M, tao->solution, tao->gradient);
 92:     MatIsSymmetricKnown(tr->M, &sym_set, &is_symmetric);
 93:     if (!sym_set || !is_symmetric) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM matrix in the LMVM preconditioner must be symmetric.");
 94:   } else if (is_jacobi) {
 95:     PCJacobiSetUseAbs(pc,PETSC_TRUE);
 96:   }

 98:   /* Check convergence criteria */
 99:   TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);
100:   TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm);
101:   if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1,"User provided compute function generated Inf or NaN");
102:   needH = 1;

104:   tao->reason = TAO_CONTINUE_ITERATING;
105:   TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);
106:   TaoMonitor(tao,tao->niter,f,gnorm,0.0,1.0);
107:   (*tao->ops->convergencetest)(tao,tao->cnvP);
108:   if (tao->reason != TAO_CONTINUE_ITERATING) return(0);

110:   /*  Initialize trust-region radius */
111:   switch (tr->init_type) {
112:   case NTR_INIT_CONSTANT:
113:     /*  Use the initial radius specified */
114:     break;

116:   case NTR_INIT_INTERPOLATION:
117:     /*  Use the initial radius specified */
118:     max_radius = 0.0;

120:     for (j = 0; j < j_max; ++j) {
121:       fmin = f;
122:       sigma = 0.0;

124:       if (needH) {
125:         TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);
126:         needH = 0;
127:       }

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

131:         VecCopy(tao->solution, tr->W);
132:         VecAXPY(tr->W, -tao->trust/gnorm, tao->gradient);
133:         TaoComputeObjective(tao, tr->W, &ftrial);

135:         if (PetscIsInfOrNanReal(ftrial)) {
136:           tau = tr->gamma1_i;
137:         }
138:         else {
139:           if (ftrial < fmin) {
140:             fmin = ftrial;
141:             sigma = -tao->trust / gnorm;
142:           }

144:           MatMult(tao->hessian, tao->gradient, tao->stepdirection);
145:           VecDot(tao->gradient, tao->stepdirection, &prered);

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

157:           tau_1 = tr->theta_i * gnorm * tao->trust / (tr->theta_i * gnorm * tao->trust + (1.0 - tr->theta_i) * prered - actred);
158:           tau_2 = tr->theta_i * gnorm * tao->trust / (tr->theta_i * gnorm * tao->trust - (1.0 + tr->theta_i) * prered + actred);
159:           tau_min = PetscMin(tau_1, tau_2);
160:           tau_max = PetscMax(tau_1, tau_2);

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

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

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

217:       if (fmin < f) {
218:         f = fmin;
219:         VecAXPY(tao->solution, sigma, tao->gradient);
220:         TaoComputeGradient(tao,tao->solution, tao->gradient);

222:         TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm);

224:         if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");
225:         needH = 1;

227:         TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);
228:         TaoMonitor(tao,tao->niter,f,gnorm,0.0,1.0);
229:         (*tao->ops->convergencetest)(tao,tao->cnvP);
230:         if (tao->reason != TAO_CONTINUE_ITERATING) {
231:           return(0);
232:         }
233:       }
234:     }
235:     tao->trust = PetscMax(tao->trust, max_radius);

237:     /*  Modify the radius if it is too large or small */
238:     tao->trust = PetscMax(tao->trust, tr->min_radius);
239:     tao->trust = PetscMin(tao->trust, tr->max_radius);
240:     break;

242:   default:
243:     /*  Norm of the first direction will initialize radius */
244:     tao->trust = 0.0;
245:     break;
246:   }

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

262:     if (tr->bfgs_pre) {
263:       /* Update the limited memory preconditioner */
264:       MatLMVMUpdate(tr->M, tao->solution, tao->gradient);
265:       ++bfgsUpdates;
266:     }

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

271:       /* Solve the trust region subproblem */
272:       KSPCGSetRadius(tao->ksp,tao->trust);
273:       KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);
274:       KSPGetIterationNumber(tao->ksp,&its);
275:       tao->ksp_its+=its;
276:       tao->ksp_tot_its+=its;
277:       KSPCGGetNormD(tao->ksp, &norm_d);

279:       if (0.0 == tao->trust) {
280:         /* Radius was uninitialized; use the norm of the direction */
281:         if (norm_d > 0.0) {
282:           tao->trust = norm_d;

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

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

297:           KSPCGSetRadius(tao->ksp,tao->trust);
298:           KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);
299:           KSPGetIterationNumber(tao->ksp,&its);
300:           tao->ksp_its+=its;
301:           tao->ksp_tot_its+=its;
302:           KSPCGGetNormD(tao->ksp, &norm_d);

304:           if (norm_d == 0.0) SETERRQ(PETSC_COMM_SELF,1, "Initial direction zero");
305:         }
306:       }
307:       VecScale(tao->stepdirection, -1.0);
308:       KSPGetConvergedReason(tao->ksp, &ksp_reason);
309:       if ((KSP_DIVERGED_INDEFINITE_PC == ksp_reason) && (tr->bfgs_pre)) {
310:         /* Preconditioner is numerically indefinite; reset the
311:            approximate if using BFGS preconditioning. */
312:         MatLMVMReset(tr->M, PETSC_FALSE);
313:         MatLMVMUpdate(tr->M, tao->solution, tao->gradient);
314:         bfgsUpdates = 1;
315:       }

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

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

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

402:             tau_1 = tr->theta * beta / (tr->theta * beta - (1.0 - tr->theta) * prered + actred);
403:             tau_2 = tr->theta * beta / (tr->theta * beta + (1.0 + tr->theta) * prered - actred);
404:             tau_min = PetscMin(tau_1, tau_2);
405:             tau_max = PetscMax(tau_1, tau_2);

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

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

464:       /* The step computed was not good and the radius was decreased.
465:          Monitor the radius to terminate. */
466:       TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);
467:       TaoMonitor(tao,tao->niter,f,gnorm,0.0,tao->trust);
468:       (*tao->ops->convergencetest)(tao,tao->cnvP);
469:     }

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

474:     if (tao->reason == TAO_CONTINUE_ITERATING) {
475:       VecCopy(tr->W, tao->solution);
476:       f = ftrial;
477:       TaoComputeGradient(tao, tao->solution, tao->gradient);
478:       TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm);
479:       if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");
480:       needH = 1;
481:       TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);
482:       TaoMonitor(tao,tao->niter,f,gnorm,0.0,tao->trust);
483:       (*tao->ops->convergencetest)(tao,tao->cnvP);
484:     }
485:   }
486:   return(0);
487: }

489: /*------------------------------------------------------------*/
490: static PetscErrorCode TaoSetUp_NTR(Tao tao)
491: {
492:   TAO_NTR *tr = (TAO_NTR *)tao->data;

496:   if (!tao->gradient) {VecDuplicate(tao->solution, &tao->gradient);}
497:   if (!tao->stepdirection) {VecDuplicate(tao->solution, &tao->stepdirection);}
498:   if (!tr->W) {VecDuplicate(tao->solution, &tr->W);}

500:   tr->bfgs_pre = 0;
501:   tr->M = 0;
502:   return(0);
503: }

505: /*------------------------------------------------------------*/
506: static PetscErrorCode TaoDestroy_NTR(Tao tao)
507: {
508:   TAO_NTR        *tr = (TAO_NTR *)tao->data;

512:   if (tao->setupcalled) {
513:     VecDestroy(&tr->W);
514:   }
515:   PetscFree(tao->data);
516:   return(0);
517: }

519: /*------------------------------------------------------------*/
520: static PetscErrorCode TaoSetFromOptions_NTR(PetscOptionItems *PetscOptionsObject,Tao tao)
521: {
522:   TAO_NTR        *tr = (TAO_NTR *)tao->data;

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

560: /*------------------------------------------------------------*/
561: /*MC
562:   TAONTR - Newton's method with trust region for unconstrained minimization.
563:   At each iteration, the Newton trust region method solves the system.
564:   NTR expects a KSP solver with a trust region radius.
565:             min_d  .5 dT Hk d + gkT d,  s.t.   ||d|| < Delta_k

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

597:   Level: beginner
598: M*/

600: PETSC_EXTERN PetscErrorCode TaoCreate_NTR(Tao tao)
601: {
602:   TAO_NTR *tr;


607:   PetscNewLog(tao,&tr);

609:   tao->ops->setup = TaoSetUp_NTR;
610:   tao->ops->solve = TaoSolve_NTR;
611:   tao->ops->setfromoptions = TaoSetFromOptions_NTR;
612:   tao->ops->destroy = TaoDestroy_NTR;

614:   /* Override default settings (unless already changed) */
615:   if (!tao->max_it_changed) tao->max_it = 50;
616:   if (!tao->trust0_changed) tao->trust0 = 100.0;
617:   tao->data = (void*)tr;

619:   /*  Standard trust region update parameters */
620:   tr->eta1 = 1.0e-4;
621:   tr->eta2 = 0.25;
622:   tr->eta3 = 0.50;
623:   tr->eta4 = 0.90;

625:   tr->alpha1 = 0.25;
626:   tr->alpha2 = 0.50;
627:   tr->alpha3 = 1.00;
628:   tr->alpha4 = 2.00;
629:   tr->alpha5 = 4.00;

631:   /*  Interpolation trust region update parameters */
632:   tr->mu1 = 0.10;
633:   tr->mu2 = 0.50;

635:   tr->gamma1 = 0.25;
636:   tr->gamma2 = 0.50;
637:   tr->gamma3 = 2.00;
638:   tr->gamma4 = 4.00;

640:   tr->theta = 0.05;

642:   /*  Interpolation parameters for initialization */
643:   tr->mu1_i = 0.35;
644:   tr->mu2_i = 0.50;

646:   tr->gamma1_i = 0.0625;
647:   tr->gamma2_i = 0.50;
648:   tr->gamma3_i = 2.00;
649:   tr->gamma4_i = 5.00;

651:   tr->theta_i = 0.25;

653:   tr->min_radius = 1.0e-10;
654:   tr->max_radius = 1.0e10;
655:   tr->epsilon    = 1.0e-6;

657:   tr->init_type       = NTR_INIT_INTERPOLATION;
658:   tr->update_type     = NTR_UPDATE_REDUCTION;

660:   /* Set linear solver to default for trust region */
661:   KSPCreate(((PetscObject)tao)->comm,&tao->ksp);
662:   PetscObjectIncrementTabLevel((PetscObject)tao->ksp,(PetscObject)tao,1);
663:   KSPSetOptionsPrefix(tao->ksp,tao->hdr.prefix);
664:   KSPAppendOptionsPrefix(tao->ksp,"tao_ntr_");
665:   KSPSetType(tao->ksp,KSPSTCG);
666:   return(0);
667: }