Actual source code: ntr.c

petsc-3.10.5 2019-03-28
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 KSPCGNASH, KSPCGSTCG,
 40:           or KSPCGGLTR.  Thus, we set KSPCGNASH, KSPCGSTCG, or KSPCGGLTR 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:     PetscPrintf(((PetscObject)tao)->comm,"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,KSPCGNASH,&is_nash);
 70:   PetscStrcmp(ksp_type,KSPCGSTCG,&is_stcg);
 71:   PetscStrcmp(ksp_type,KSPCGGLTR,&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 - 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 - 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:     ++tao->niter;
251:     tao->ksp_its=0;
252:     /* Compute the Hessian */
253:     if (needH) {
254:       TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);
255:       needH = 0;
256:     }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

492:   if (!tao->gradient) {VecDuplicate(tao->solution, &tao->gradient);}
493:   if (!tao->stepdirection) {VecDuplicate(tao->solution, &tao->stepdirection);}
494:   if (!tr->W) {VecDuplicate(tao->solution, &tr->W);}

496:   tr->bfgs_pre = 0;
497:   tr->M = 0;
498:   return(0);
499: }

501: /*------------------------------------------------------------*/
502: static PetscErrorCode TaoDestroy_NTR(Tao tao)
503: {
504:   TAO_NTR        *tr = (TAO_NTR *)tao->data;

508:   if (tao->setupcalled) {
509:     VecDestroy(&tr->W);
510:   }
511:   PetscFree(tao->data);
512:   return(0);
513: }

515: /*------------------------------------------------------------*/
516: static PetscErrorCode TaoSetFromOptions_NTR(PetscOptionItems *PetscOptionsObject,Tao tao)
517: {
518:   TAO_NTR        *tr = (TAO_NTR *)tao->data;

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

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

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

593:   Level: beginner
594: M*/

596: PETSC_EXTERN PetscErrorCode TaoCreate_NTR(Tao tao)
597: {
598:   TAO_NTR *tr;


603:   PetscNewLog(tao,&tr);

605:   tao->ops->setup = TaoSetUp_NTR;
606:   tao->ops->solve = TaoSolve_NTR;
607:   tao->ops->setfromoptions = TaoSetFromOptions_NTR;
608:   tao->ops->destroy = TaoDestroy_NTR;

610:   /* Override default settings (unless already changed) */
611:   if (!tao->max_it_changed) tao->max_it = 50;
612:   if (!tao->trust0_changed) tao->trust0 = 100.0;
613:   tao->data = (void*)tr;

615:   /*  Standard trust region update parameters */
616:   tr->eta1 = 1.0e-4;
617:   tr->eta2 = 0.25;
618:   tr->eta3 = 0.50;
619:   tr->eta4 = 0.90;

621:   tr->alpha1 = 0.25;
622:   tr->alpha2 = 0.50;
623:   tr->alpha3 = 1.00;
624:   tr->alpha4 = 2.00;
625:   tr->alpha5 = 4.00;

627:   /*  Interpolation trust region update parameters */
628:   tr->mu1 = 0.10;
629:   tr->mu2 = 0.50;

631:   tr->gamma1 = 0.25;
632:   tr->gamma2 = 0.50;
633:   tr->gamma3 = 2.00;
634:   tr->gamma4 = 4.00;

636:   tr->theta = 0.05;

638:   /*  Interpolation parameters for initialization */
639:   tr->mu1_i = 0.35;
640:   tr->mu2_i = 0.50;

642:   tr->gamma1_i = 0.0625;
643:   tr->gamma2_i = 0.50;
644:   tr->gamma3_i = 2.00;
645:   tr->gamma4_i = 5.00;

647:   tr->theta_i = 0.25;

649:   tr->min_radius = 1.0e-10;
650:   tr->max_radius = 1.0e10;
651:   tr->epsilon    = 1.0e-6;

653:   tr->init_type       = NTR_INIT_INTERPOLATION;
654:   tr->update_type     = NTR_UPDATE_REDUCTION;

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