Actual source code: ntr.c

  1: /*$Id$*/

  3: #include "ntr.h"

  5: #ifdef TAO_USE_PETSC
  6: #include "petscksp.h"
  7: #include "petscpc.h"
  8: #include "src/petsctao/linearsolver/taolinearsolver_petsc.h"
  9: #include "src/petsctao/vector/taovec_petsc.h"

 11: #include "private/kspimpl.h"
 12: #include "private/pcimpl.h"

 14: #define NTR_KSP_NASH    0
 15: #define NTR_KSP_STCG    1
 16: #define NTR_KSP_GLTR    2
 17: #define NTR_KSP_TYPES   3

 19: #define NTR_PC_NONE        0
 20: #define NTR_PC_AHESS    1
 21: #define NTR_PC_BFGS     2
 22: #define NTR_PC_PETSC    3
 23: #define NTR_PC_TYPES    4

 25: #define BFGS_SCALE_AHESS   0
 26: #define BFGS_SCALE_BFGS    1
 27: #define BFGS_SCALE_TYPES   2

 29: #define NTR_INIT_CONSTANT          0
 30: #define NTR_INIT_DIRECTION          1
 31: #define NTR_INIT_INTERPOLATION          2
 32: #define NTR_INIT_TYPES                  3

 34: #define NTR_UPDATE_REDUCTION      0
 35: #define NTR_UPDATE_INTERPOLATION  1
 36: #define NTR_UPDATE_TYPES          2

 38: static const char *NTR_KSP[64] = {
 39:   "nash", "stcg", "gltr"
 40: };

 42: static const char *NTR_PC[64] = {
 43:   "none", "ahess", "bfgs", "petsc"
 44: };

 46: static const char *BFGS_SCALE[64] = {
 47:   "ahess", "bfgs"
 48: };

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

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

 58: // Routine for BFGS preconditioner

 62: static PetscErrorCode bfgs_apply(PC pc, Vec xin, Vec xout)
 63: {
 64:   TaoLMVMMat *M;
 65:   TaoVecPetsc Xin(xin);
 66:   TaoVecPetsc Xout(xout);
 67:   TaoTruth info2;
 68:   int info;

 71:   info = PCShellGetContext(pc,(void**)&M); CHKERRQ(info);
 72:   info = M->Solve(&Xin, &Xout, &info2); CHKERRQ(info);
 73:   return(0);
 74: }

 76: // TaoSolve_NTR - Implements Newton's Method with a trust region approach 
 77: // for solving unconstrained minimization problems.  
 78: //
 79: // The basic algorithm is taken from MINPACK-2 (dstrn).
 80: //
 81: // TaoSolve_NTR computes a local minimizer of a twice differentiable function
 82: // f by applying a trust region variant of Newton's method.  At each stage 
 83: // of the algorithm, we use the prconditioned conjugate gradient method to
 84: // determine an approximate minimizer of the quadratic equation
 85: //
 86: //      q(s) = <s, Hs + g>
 87: //
 88: // subject to the trust region constraint
 89: //
 90: //      || s ||_M <= radius,
 91: //
 92: // where radius is the trust region radius and M is a symmetric positive
 93: // definite matrix (the preconditioner).  Here g is the gradient and H 
 94: // is the Hessian matrix.
 95: //
 96: // Note:  TaoSolve_NTR MUST use the iterative solver KSPNASH, KSPSTCG,
 97: //        or KSPGLTR.  Thus, we set KSPNASH, KSPSTCG, or KSPGLTR in this 
 98: //        routine regardless of what the user may have previously specified.

102: static int TaoSolve_NTR(TAO_SOLVER tao, void *solver)
103: {
104:   TAO_NTR *tr = (TAO_NTR *)solver;
105:   TaoVec *X, *G = tr->G, *D = tr->D, *W = tr->W;
106:   TaoVec *Diag = tr->Diag;
107:   TaoMat *H;
108:   TaoLMVMMat *M = tr->M;

110:   TaoLinearSolver *tls;
111:   TaoLinearSolverPetsc *pls;

113:   KSP pksp;
114:   PC ppc;

116:   KSPConvergedReason ksp_reason;
117:   TaoTerminateReason reason;
118:   TaoTruth success;

120:   double fmin, ftrial, prered, actred, kappa, sigma, beta;
121:   double tau, tau_1, tau_2, tau_max, tau_min, max_radius;
122:   double f, gnorm;

124:   double delta;
125:   double radius, norm_d;
126:   int info;

128:   TaoInt iter = 0;
129:   TaoInt bfgsUpdates = 0;
130:   TaoInt needH;

132:   TaoInt i_max = 5;
133:   TaoInt j_max = 1;
134:   TaoInt i, j;

136:   TaoFunctionBegin;

138:   // Get the initial trust-region radius
139:   info = TaoGetInitialTrustRegionRadius(tao, &radius); CHKERRQ(info);
140:   if (radius < 0.0) {
141:     SETERRQ(1, "Initial radius negative");
142:   }

144:   // Modify the radius if it is too large or small
145:   radius = TaoMax(radius, tr->min_radius);
146:   radius = TaoMin(radius, tr->max_radius);

148:   // Get vectors we will need
149:   info = TaoGetSolution(tao, &X); CHKERRQ(info);
150:   info = TaoGetHessian(tao, &H); CHKERRQ(info);

152:   if (NTR_PC_BFGS == tr->pc_type && !M) {
153:     tr->M = new TaoLMVMMat(X);
154:     M = tr->M;
155:   }

157:   // Check convergence criteria
158:   info = TaoComputeFunctionGradient(tao, X, &f, G); CHKERRQ(info);
159:   info = G->Norm2(&gnorm); CHKERRQ(info);
160:   if (TaoInfOrNaN(f) || TaoInfOrNaN(gnorm)) {
161:     SETERRQ(1, "User provided compute function generated Inf or NaN");
162:   }
163:   needH = 1;

165:   info = TaoMonitor(tao, iter, f, gnorm, 0.0, 1.0, &reason); CHKERRQ(info);
166:   if (reason != TAO_CONTINUE_ITERATING) {
167:     TaoFunctionReturn(0);
168:   }

170:   // Create vectors for the limited memory preconditioner
171:   if ((NTR_PC_BFGS == tr->pc_type) && 
172:       (BFGS_SCALE_BFGS != tr->bfgs_scale_type)) {
173:     if (!Diag) {
174:       info = X->Clone(&tr->Diag); CHKERRQ(info);
175:       Diag = tr->Diag;
176:     }
177:   }

179:   // Modify the linear solver to a conjugate gradient method
180:   info = TaoGetLinearSolver(tao, &tls); CHKERRQ(info);
181:   pls  = dynamic_cast <TaoLinearSolverPetsc *> (tls);

183:   pksp = pls->GetKSP();
184:   switch(tr->ksp_type) {
185:   case NTR_KSP_NASH:
186:     info = KSPSetType(pksp, KSPNASH); CHKERRQ(info);
187:     if (pksp->ops->setfromoptions) {
188:       (*pksp->ops->setfromoptions)(pksp);
189:     }
190:     break;

192:   case NTR_KSP_STCG:
193:     info = KSPSetType(pksp, KSPSTCG); CHKERRQ(info);
194:     if (pksp->ops->setfromoptions) {
195:       (*pksp->ops->setfromoptions)(pksp);
196:     }
197:     break;

199:   default:
200:     info = KSPSetType(pksp, KSPGLTR); CHKERRQ(info);
201:     if (pksp->ops->setfromoptions) {
202:       (*pksp->ops->setfromoptions)(pksp);
203:     }
204:     break;
205:   }

207:   // Modify the preconditioner to use the bfgs approximation
208:   info = KSPGetPC(pksp, &ppc); CHKERRQ(info);
209:   switch(tr->pc_type) {
210:   case NTR_PC_NONE:
211:     info = PCSetType(ppc, PCNONE); CHKERRQ(info);
212:     if (ppc->ops->setfromoptions) {
213:       (*ppc->ops->setfromoptions)(ppc);
214:     }
215:     break;
216:  
217:   case NTR_PC_AHESS:
218:     info = PCSetType(ppc, PCJACOBI); CHKERRQ(info);
219:     if (ppc->ops->setfromoptions) {
220:       (*ppc->ops->setfromoptions)(ppc);
221:     }
222:     info = PCJacobiSetUseAbs(ppc); CHKERRQ(info);
223:     break;

225:   case NTR_PC_BFGS:
226:     info = PCSetType(ppc, PCSHELL); CHKERRQ(info);
227:     if (ppc->ops->setfromoptions) {
228:       (*ppc->ops->setfromoptions)(ppc);
229:     }
230:     info = PCShellSetName(ppc, "bfgs"); CHKERRQ(info);
231:     info = PCShellSetContext(ppc, M); CHKERRQ(info);
232:     info = PCShellSetApply(ppc, bfgs_apply); CHKERRQ(info);
233:     break;

235:   default:
236:     // Use the pc method set by pc_type
237:     break;
238:   }

240:   // Initialize trust-region radius
241:   switch(tr->init_type) {
242:   case NTR_INIT_CONSTANT:
243:     // Use the initial radius specified
244:     break;

246:   case NTR_INIT_INTERPOLATION:
247:     // Use the initial radius specified
248:     max_radius = 0.0;

250:     for (j = 0; j < j_max; ++j) {
251:       fmin = f;
252:       sigma = 0.0;

254:       if (needH) {
255:         info = TaoComputeHessian(tao, X, H); CHKERRQ(info);
256:         needH = 0;
257:       }

259:       for (i = 0; i < i_max; ++i) {
260:         info = W->Waxpby(1.0, X, -radius / gnorm, G); CHKERRQ(info);

262:         info = TaoComputeFunction(tao, W, &ftrial); CHKERRQ(info);
263:         if (TaoInfOrNaN(ftrial)) {
264:           tau = tr->gamma1_i;
265:         }
266:         else {
267:           if (ftrial < fmin) {
268:             fmin = ftrial;
269:             sigma = -radius / gnorm;
270:           }

272:           info = H->Multiply(G, D); CHKERRQ(info);
273:           info = D->Dot(G, &prered); CHKERRQ(info);

275:           prered = radius * (gnorm - 0.5 * radius * prered / (gnorm * gnorm));
276:           actred = f - ftrial;
277:           if ((fabs(actred) <= tr->epsilon) && 
278:               (fabs(prered) <= tr->epsilon)) {
279:             kappa = 1.0;
280:           }
281:           else {
282:             kappa = actred / prered;
283:           }

285:           tau_1 = tr->theta_i * gnorm * radius / (tr->theta_i * gnorm * radius + (1.0 - tr->theta_i) * prered - actred);
286:           tau_2 = tr->theta_i * gnorm * radius / (tr->theta_i * gnorm * radius - (1.0 + tr->theta_i) * prered + actred);
287:           tau_min = TaoMin(tau_1, tau_2);
288:           tau_max = TaoMax(tau_1, tau_2);

290:           if (fabs(kappa - 1.0) <= tr->mu1_i) {
291:             // Great agreement
292:             max_radius = TaoMax(max_radius, radius);

294:             if (tau_max < 1.0) {
295:               tau = tr->gamma3_i;
296:             }
297:             else if (tau_max > tr->gamma4_i) {
298:               tau = tr->gamma4_i;
299:             }
300:             else {
301:               tau = tau_max;
302:             }
303:           }
304:           else if (fabs(kappa - 1.0) <= tr->mu2_i) {
305:             // Good agreement
306:             max_radius = TaoMax(max_radius, radius);

308:             if (tau_max < tr->gamma2_i) {
309:               tau = tr->gamma2_i;
310:             }
311:             else if (tau_max > tr->gamma3_i) {
312:               tau = tr->gamma3_i;
313:             }
314:             else {
315:               tau = tau_max;
316:             }
317:           }
318:           else {
319:             // Not good agreement
320:             if (tau_min > 1.0) {
321:               tau = tr->gamma2_i;
322:             }
323:             else if (tau_max < tr->gamma1_i) {
324:               tau = tr->gamma1_i;
325:             }
326:             else if ((tau_min < tr->gamma1_i) && (tau_max >= 1.0)) {
327:               tau = tr->gamma1_i;
328:             }
329:             else if ((tau_1 >= tr->gamma1_i) && (tau_1 < 1.0) && 
330:                      ((tau_2 < tr->gamma1_i) || (tau_2 >= 1.0))) {
331:               tau = tau_1;
332:             }
333:             else if ((tau_2 >= tr->gamma1_i) && (tau_2 < 1.0) && 
334:                      ((tau_1 < tr->gamma1_i) || (tau_2 >= 1.0))) {
335:               tau = tau_2;
336:             }
337:             else {
338:               tau = tau_max;
339:             }
340:           }
341:         }
342:         radius = tau * radius;
343:       }
344:   
345:       if (fmin < f) {
346:         f = fmin;
347:         info = X->Axpy(sigma, G); CHKERRQ(info);
348:         info = TaoComputeGradient(tao, X, G); CHKERRQ(info);

350:         info = G->Norm2(&gnorm); CHKERRQ(info);
351:         if (TaoInfOrNaN(f) || TaoInfOrNaN(gnorm)) {
352:           SETERRQ(1, "User provided compute function generated Inf or NaN");
353:         }
354:         needH = 1;

356:         info = TaoMonitor(tao, iter, f, gnorm, 0.0, 1.0, &reason); CHKERRQ(info);
357:         if (reason != TAO_CONTINUE_ITERATING) {
358:           TaoFunctionReturn(0);
359:         }
360:       }
361:     }
362:     radius = TaoMax(radius, max_radius);

364:     // Modify the radius if it is too large or small
365:     radius = TaoMax(radius, tr->min_radius);
366:     radius = TaoMin(radius, tr->max_radius);
367:     break;

369:   default:
370:     // Norm of the first direction will initialize radius
371:     radius = 0.0;
372:     break;
373:   }

375:   // Set initial scaling for the BFGS preconditioner 
376:   // This step is done after computing the initial trust-region radius
377:   // since the function value may have decreased
378:   if (NTR_PC_BFGS == tr->pc_type) {
379:     if (f != 0.0) {
380:       delta = 2.0 * TaoAbsDouble(f) / (gnorm*gnorm);
381:     }
382:     else {
383:       delta = 2.0 / (gnorm*gnorm);
384:     }
385:     info = M->SetDelta(delta); CHKERRQ(info);
386:   }

388:   // Have not converged; continue with Newton method
389:   while (reason == TAO_CONTINUE_ITERATING) {
390:     ++iter;

392:     // Compute the Hessian
393:     if (needH) {
394:       info = TaoComputeHessian(tao, X, H); CHKERRQ(info);
395:       needH = 0;
396:     }

398:     if (NTR_PC_BFGS == tr->pc_type) {
399:       if (BFGS_SCALE_AHESS == tr->bfgs_scale_type) {
400:         // Obtain diagonal for the bfgs preconditioner
401:         info = H->GetDiagonal(Diag); CHKERRQ(info);
402:         info = Diag->AbsoluteValue(); CHKERRQ(info);
403:         info = Diag->Reciprocal(); CHKERRQ(info);
404:         info = M->SetScale(Diag); CHKERRQ(info);
405:       }

407:       // Update the limited memory preconditioner
408:       info = M->Update(X, G); CHKERRQ(info);
409:       ++bfgsUpdates;
410:     }

412:     while (reason == TAO_CONTINUE_ITERATING) {
413:       // Solve the trust region subproblem
414:       info = TaoPreLinearSolve(tao, H); CHKERRQ(info);
415:       info = TaoLinearSolveTrustRegion(tao, H, G, D, radius, &success); CHKERRQ(info);
416:       info = pls->GetNormDirection(&norm_d); CHKERRQ(info);
417:       if (0.0 == radius) {
418:         // Radius was uninitialized; use the norm of the direction
419:         if (norm_d > 0.0) {
420:           radius = norm_d;

422:           // Modify the radius if it is too large or small
423:           radius = TaoMax(radius, tr->min_radius);
424:           radius = TaoMin(radius, tr->max_radius);
425:         }
426:         else {
427:           // The direction was bad; set radius to default value and re-solve 
428:           // the trust-region subproblem to get a direction
429:           info = TaoGetInitialTrustRegionRadius(tao, &radius); CHKERRQ(info);

431:           // Modify the radius if it is too large or small
432:           radius = TaoMax(radius, tr->min_radius);
433:           radius = TaoMin(radius, tr->max_radius);

435:           info = TaoLinearSolveTrustRegion(tao, H, G, D, radius, &success); CHKERRQ(info);
436:           info = pls->GetNormDirection(&norm_d); CHKERRQ(info);
437:           if (norm_d == 0.0) {
438:             SETERRQ(1, "Initial direction zero");
439:           }
440:         }
441:       }
442:       info = D->Negate(); CHKERRQ(info);

444:       info = KSPGetConvergedReason(pksp, &ksp_reason); CHKERRQ(info);
445:       if ((KSP_DIVERGED_INDEFINITE_PC == ksp_reason) &&
446:           (NTR_PC_BFGS == tr->pc_type) && (bfgsUpdates > 1)) {
447:         // Preconditioner is numerically indefinite; reset the
448:         // approximate if using BFGS preconditioning.
449:   
450:         if (f != 0.0) {
451:           delta = 2.0 * TaoAbsDouble(f) / (gnorm*gnorm);
452:         }
453:         else {
454:           delta = 2.0 / (gnorm*gnorm);
455:         }
456:         info = M->SetDelta(delta); CHKERRQ(info);
457:         info = M->Reset(); CHKERRQ(info);
458:         info = M->Update(X, G); CHKERRQ(info);
459:         bfgsUpdates = 1;
460:       }

462:       if (NTR_UPDATE_REDUCTION == tr->update_type) {
463:         // Get predicted reduction
464:         info = pls->GetObjFcn(&prered); CHKERRQ(info);
465:         
466:         if (prered >= 0.0) {
467:           // The predicted reduction has the wrong sign.  This cannot
468:           // happen in infinite precision arithmetic.  Step should
469:           // be rejected!
470:           radius = tr->alpha1 * TaoMin(radius, norm_d);
471:         }
472:         else {
473:           // Compute trial step and function value
474:           info = W->Waxpby(1.0, X, 1.0, D); CHKERRQ(info);
475:           info = TaoComputeFunction(tao, W, &ftrial); CHKERRQ(info);
476:           if (TaoInfOrNaN(ftrial)) {
477:             radius = tr->alpha1 * TaoMin(radius, norm_d);
478:           }
479:           else {
480:             // Compute and actual reduction
481:             actred = f - ftrial;
482:             prered = -prered;
483:             if ((fabs(actred) <= tr->epsilon) && 
484:                 (fabs(prered) <= tr->epsilon)) {
485:               kappa = 1.0;
486:             }
487:             else {
488:               kappa = actred / prered;
489:             }
490:             
491:             // Accept or reject the step and update radius
492:             if (kappa < tr->eta1) {
493:               // Reject the step
494:               radius = tr->alpha1 * TaoMin(radius, norm_d);
495:             } 
496:             else {
497:               // Accept the step
498:               if (kappa < tr->eta2) { 
499:                 // Marginal bad step
500:                 radius = tr->alpha2 * TaoMin(radius, norm_d);
501:               }
502:               else if (kappa < tr->eta3) {
503:                 // Reasonable step
504:                 radius = tr->alpha3 * radius;
505:               }
506:               else if (kappa < tr->eta4) { 
507:                 // Good step
508:                 radius = TaoMax(tr->alpha4 * norm_d, radius);
509:               }
510:               else {
511:                 // Very good step
512:                 radius = TaoMax(tr->alpha5 * norm_d, radius);
513:               }
514:               break;
515:             }
516:           }
517:         } 
518:       }
519:       else {
520:         // Get predicted reduction
521:         info = pls->GetObjFcn(&prered); CHKERRQ(info);

523:         if (prered >= 0.0) {
524:           // The predicted reduction has the wrong sign.  This cannot
525:           // happen in infinite precision arithmetic.  Step should
526:           // be rejected!
527:           radius = tr->gamma1 * TaoMin(radius, norm_d);
528:         }
529:         else {
530:           info = W->Waxpby(1.0, X, 1.0, D); CHKERRQ(info);
531:           info = TaoComputeFunction(tao, W, &ftrial); CHKERRQ(info);
532:           if (TaoInfOrNaN(ftrial)) {
533:             radius = tr->gamma1 * TaoMin(radius, norm_d);
534:           }
535:           else {
536:             info = D->Dot(G, &beta); CHKERRQ(info);

538:             actred = f - ftrial;
539:             prered = -prered;
540:             if ((fabs(actred) <= tr->epsilon) && 
541:                 (fabs(prered) <= tr->epsilon)) {
542:               kappa = 1.0;
543:             }
544:             else {
545:               kappa = actred / prered;
546:             }

548:             tau_1 = tr->theta * beta / (tr->theta * beta - (1.0 - tr->theta) * prered + actred);
549:             tau_2 = tr->theta * beta / (tr->theta * beta + (1.0 + tr->theta) * prered - actred);
550:             tau_min = TaoMin(tau_1, tau_2);
551:             tau_max = TaoMax(tau_1, tau_2);

553:             if (kappa >= 1.0 - tr->mu1) {
554:               // Great agreement; accept step and update radius
555:               if (tau_max < 1.0) {
556:                 radius = TaoMax(radius, tr->gamma3 * norm_d);
557:               }
558:               else if (tau_max > tr->gamma4) {
559:                 radius = TaoMax(radius, tr->gamma4 * norm_d);
560:               }
561:               else {
562:                 radius = TaoMax(radius, tau_max * norm_d);
563:               }
564:               break;
565:             }
566:             else if (kappa >= 1.0 - tr->mu2) {
567:               // Good agreement

569:               if (tau_max < tr->gamma2) {
570:                 radius = tr->gamma2 * TaoMin(radius, norm_d);
571:               }
572:               else if (tau_max > tr->gamma3) {
573:                 radius = TaoMax(radius, tr->gamma3 * norm_d);
574:               }
575:               else if (tau_max < 1.0) {
576:                 radius = tau_max * TaoMin(radius, norm_d);
577:               }
578:               else {
579:                 radius = TaoMax(radius, tau_max * norm_d);
580:               }
581:               break;
582:             }
583:             else {
584:               // Not good agreement
585:               if (tau_min > 1.0) {
586:                 radius = tr->gamma2 * TaoMin(radius, norm_d);
587:               }
588:               else if (tau_max < tr->gamma1) {
589:                 radius = tr->gamma1 * TaoMin(radius, norm_d);
590:               }
591:               else if ((tau_min < tr->gamma1) && (tau_max >= 1.0)) {
592:                 radius = tr->gamma1 * TaoMin(radius, norm_d);
593:               }
594:               else if ((tau_1 >= tr->gamma1) && (tau_1 < 1.0) && 
595:                        ((tau_2 < tr->gamma1) || (tau_2 >= 1.0))) {
596:                 radius = tau_1 * TaoMin(radius, norm_d);
597:               }
598:               else if ((tau_2 >= tr->gamma1) && (tau_2 < 1.0) && 
599:                        ((tau_1 < tr->gamma1) || (tau_2 >= 1.0))) {
600:                 radius = tau_2 * TaoMin(radius, norm_d);
601:               }
602:               else {
603:                 radius = tau_max * TaoMin(radius, norm_d);
604:               }
605:             }
606:           }
607:         }
608:       }

610:       // The step computed was not good and the radius was decreased.
611:       // Monitor the radius to terminate.
612:       info = TaoMonitor(tao, iter, f, gnorm, 0.0, radius, &reason); CHKERRQ(info);
613:     }

615:     // The radius may have been increased; modify if it is too large
616:     radius = TaoMin(radius, tr->max_radius);

618:     if (reason == TAO_CONTINUE_ITERATING) {
619:       info = X->CopyFrom(W); CHKERRQ(info);
620:       f = ftrial;
621:       info = TaoComputeGradient(tao, X, G); CHKERRQ(info);
622:       info = G->Norm2(&gnorm); CHKERRQ(info);
623:       if (TaoInfOrNaN(f) || TaoInfOrNaN(gnorm)) {
624:         SETERRQ(1, "User provided compute function generated Inf or NaN");
625:       }
626:       needH = 1;

628:       info = TaoMonitor(tao, iter, f, gnorm, 0.0, radius, &reason); CHKERRQ(info);
629:     }
630:   }
631:   TaoFunctionReturn(0);
632: }

634: /*------------------------------------------------------------*/
637: static int TaoSetUp_NTR(TAO_SOLVER tao, void *solver)
638: {
639:   TAO_NTR *tr = (TAO_NTR *)solver;
640:   TaoVec *X;
641:   TaoMat *H;
642:   int info;

644:   TaoFunctionBegin;

646:   info = TaoGetSolution(tao, &X); CHKERRQ(info);
647:   info = X->Clone(&tr->G); CHKERRQ(info);
648:   info = X->Clone(&tr->D); CHKERRQ(info);
649:   info = X->Clone(&tr->W); CHKERRQ(info);

651:   tr->Diag = 0;
652:   tr->M = 0;

654:   info = TaoSetLagrangianGradientVector(tao, tr->G); CHKERRQ(info);
655:   info = TaoSetStepDirectionVector(tao, tr->D); CHKERRQ(info);

657:   // Set linear solver to default for trust region
658:   info = TaoGetHessian(tao, &H); CHKERRQ(info);
659:   info = TaoCreateLinearSolver(tao, H, 200, 0); CHKERRQ(info); 

661:   // Check sizes for compatability
662:   info = TaoCheckFGH(tao); CHKERRQ(info);
663:   TaoFunctionReturn(0);
664: }

666: /*------------------------------------------------------------*/
669: static int TaoDestroy_NTR(TAO_SOLVER tao, void *solver)
670: {
671:   TAO_NTR *tr = (TAO_NTR *)solver;
672:   int info;

674:   TaoFunctionBegin;
675:   info = TaoVecDestroy(tr->G); CHKERRQ(info);
676:   info = TaoVecDestroy(tr->D); CHKERRQ(info);
677:   info = TaoVecDestroy(tr->W); CHKERRQ(info);

679:   info = TaoSetLagrangianGradientVector(tao, 0); CHKERRQ(info);
680:   info = TaoSetStepDirectionVector(tao, 0); CHKERRQ(info);

682:   if (tr->Diag) {
683:     info = TaoVecDestroy(tr->Diag); CHKERRQ(info);
684:     tr->Diag = 0;
685:   }

687:   if (tr->M) {
688:     info = TaoMatDestroy(tr->M); CHKERRQ(info);
689:     tr->M = 0;
690:   }
691:   TaoFunctionReturn(0);
692: }

694: /*------------------------------------------------------------*/
697: static int TaoSetOptions_NTR(TAO_SOLVER tao, void*solver)
698: {
699:   TAO_NTR *tr = (TAO_NTR *)solver;
700:   int info;

702:   TaoFunctionBegin;
703:   info = TaoOptionsHead("Newton trust region method for unconstrained optimization"); CHKERRQ(info);
704:   info = TaoOptionList("-tao_ntr_ksp_type", "ksp type", "", NTR_KSP, NTR_KSP_TYPES, NTR_KSP[tr->ksp_type], &tr->ksp_type, 0); CHKERRQ(info);
705:   info = TaoOptionList("-tao_ntr_pc_type", "pc type", "", NTR_PC, NTR_PC_TYPES, NTR_PC[tr->pc_type], &tr->pc_type, 0); CHKERRQ(info);
706:   info = TaoOptionList("-tao_ntr_bfgs_scale_type", "bfgs scale type", "", BFGS_SCALE, BFGS_SCALE_TYPES, BFGS_SCALE[tr->bfgs_scale_type], &tr->bfgs_scale_type, 0); CHKERRQ(info);
707:   info = TaoOptionList("-tao_ntr_init_type", "radius initialization type", "", NTR_INIT, NTR_INIT_TYPES, NTR_INIT[tr->init_type], &tr->init_type, 0); CHKERRQ(info);
708:   info = TaoOptionList("-tao_ntr_update_type", "radius update type", "", NTR_UPDATE, NTR_UPDATE_TYPES, NTR_UPDATE[tr->update_type], &tr->update_type, 0); CHKERRQ(info);
709:   info = TaoOptionDouble("-tao_ntr_eta1", "step is unsuccessful if actual reduction < eta1 * predicted reduction", "", tr->eta1, &tr->eta1, 0); CHKERRQ(info);
710:   info = TaoOptionDouble("-tao_ntr_eta2", "", "", tr->eta2, &tr->eta2, 0); CHKERRQ(info);
711:   info = TaoOptionDouble("-tao_ntr_eta3", "", "", tr->eta3, &tr->eta3, 0); CHKERRQ(info);
712:   info = TaoOptionDouble("-tao_ntr_eta4", "", "", tr->eta4, &tr->eta4, 0); CHKERRQ(info);
713:   info = TaoOptionDouble("-tao_ntr_alpha1", "", "", tr->alpha1, &tr->alpha1, 0); CHKERRQ(info);
714:   info = TaoOptionDouble("-tao_ntr_alpha2", "", "", tr->alpha2, &tr->alpha2, 0); CHKERRQ(info);
715:   info = TaoOptionDouble("-tao_ntr_alpha3", "", "", tr->alpha3, &tr->alpha3, 0); CHKERRQ(info);
716:   info = TaoOptionDouble("-tao_ntr_alpha4", "", "", tr->alpha4, &tr->alpha4, 0); CHKERRQ(info);
717:   info = TaoOptionDouble("-tao_ntr_alpha5", "", "", tr->alpha5, &tr->alpha5, 0); CHKERRQ(info);
718:   info = TaoOptionDouble("-tao_ntr_mu1", "", "", tr->mu1, &tr->mu1, 0); CHKERRQ(info);
719:   info = TaoOptionDouble("-tao_ntr_mu2", "", "", tr->mu2, &tr->mu2, 0); CHKERRQ(info);
720:   info = TaoOptionDouble("-tao_ntr_gamma1", "", "", tr->gamma1, &tr->gamma1, 0); CHKERRQ(info);
721:   info = TaoOptionDouble("-tao_ntr_gamma2", "", "", tr->gamma2, &tr->gamma2, 0); CHKERRQ(info);
722:   info = TaoOptionDouble("-tao_ntr_gamma3", "", "", tr->gamma3, &tr->gamma3, 0); CHKERRQ(info);
723:   info = TaoOptionDouble("-tao_ntr_gamma4", "", "", tr->gamma4, &tr->gamma4, 0); CHKERRQ(info);
724:   info = TaoOptionDouble("-tao_ntr_theta", "", "", tr->theta, &tr->theta, 0); CHKERRQ(info);
725:   info = TaoOptionDouble("-tao_ntr_mu1_i", "", "", tr->mu1_i, &tr->mu1_i, 0); CHKERRQ(info);
726:   info = TaoOptionDouble("-tao_ntr_mu2_i", "", "", tr->mu2_i, &tr->mu2_i, 0); CHKERRQ(info);
727:   info = TaoOptionDouble("-tao_ntr_gamma1_i", "", "", tr->gamma1_i, &tr->gamma1_i, 0); CHKERRQ(info);
728:   info = TaoOptionDouble("-tao_ntr_gamma2_i", "", "", tr->gamma2_i, &tr->gamma2_i, 0); CHKERRQ(info);
729:   info = TaoOptionDouble("-tao_ntr_gamma3_i", "", "", tr->gamma3_i, &tr->gamma3_i, 0); CHKERRQ(info);
730:   info = TaoOptionDouble("-tao_ntr_gamma4_i", "", "", tr->gamma4_i, &tr->gamma4_i, 0); CHKERRQ(info);
731:   info = TaoOptionDouble("-tao_ntr_theta_i", "", "", tr->theta_i, &tr->theta_i, 0); CHKERRQ(info);
732:   info = TaoOptionDouble("-tao_ntr_min_radius", "lower bound on initial trust-region radius", "", tr->min_radius, &tr->min_radius, 0); CHKERRQ(info);
733:   info = TaoOptionDouble("-tao_ntr_max_radius", "upper bound on trust-region radius", "", tr->max_radius, &tr->max_radius, 0); CHKERRQ(info);
734:   info = TaoOptionDouble("-tao_ntr_epsilon", "tolerance used when computing actual and predicted reduction", "", tr->epsilon, &tr->epsilon, 0); CHKERRQ(info);
735:   info = TaoOptionsTail(); CHKERRQ(info);
736:   TaoFunctionReturn(0);
737: }

739: /*------------------------------------------------------------*/
742: static int TaoView_NTR(TAO_SOLVER tao,void*solver)
743: {
744:   TAO_NTR *tr = (TAO_NTR *)solver;
745:   int info;

747:   TaoFunctionBegin;
748:   if (NTR_PC_BFGS == tr->pc_type && tr->M) {
749:     info = TaoPrintInt(tao, "  Rejected matrix updates: %d\n", tr->M->GetRejects()); CHKERRQ(info);
750:   }
751:   TaoFunctionReturn(0);
752: }

754: /*------------------------------------------------------------*/
758: int TaoCreate_NTR(TAO_SOLVER tao)
759: {
760:   TAO_NTR *tr;
761:   int info;

763:   TaoFunctionBegin;

765:   info = TaoNew(TAO_NTR, &tr); CHKERRQ(info);
766:   info = PetscLogObjectMemory(tao, sizeof(TAO_NTR)); CHKERRQ(info);

768:   info = TaoSetTaoSolveRoutine(tao, TaoSolve_NTR, (void *)tr); CHKERRQ(info);
769:   info = TaoSetTaoSetUpDownRoutines(tao, TaoSetUp_NTR, TaoDestroy_NTR); CHKERRQ(info);
770:   info = TaoSetTaoOptionsRoutine(tao, TaoSetOptions_NTR); CHKERRQ(info);
771:   info = TaoSetTaoViewRoutine(tao, TaoView_NTR); CHKERRQ(info);

773:   info = TaoSetMaximumIterates(tao, 50); CHKERRQ(info);
774:   info = TaoSetTolerances(tao, 1e-10, 1e-10, 0, 0); CHKERRQ(info);

776:   info = TaoSetTrustRegionRadius(tao, 100.0); CHKERRQ(info);
777:   info = TaoSetTrustRegionTolerance(tao, 1.0e-12); CHKERRQ(info);

779:   // Standard trust region update parameters
780:   tr->eta1 = 1.0e-4;
781:   tr->eta2 = 0.25;
782:   tr->eta3 = 0.50;
783:   tr->eta4 = 0.90;

785:   tr->alpha1 = 0.25;
786:   tr->alpha2 = 0.50;
787:   tr->alpha3 = 1.00;
788:   tr->alpha4 = 2.00;
789:   tr->alpha5 = 4.00;

791:   // Interpolation parameters
792:   tr->mu1_i = 0.35;
793:   tr->mu2_i = 0.50;

795:   tr->gamma1_i = 0.0625;
796:   tr->gamma2_i = 0.50;
797:   tr->gamma3_i = 2.00;
798:   tr->gamma4_i = 5.00;

800:   tr->theta_i = 0.25;

802:   // Interpolation trust region update parameters
803:   tr->mu1 = 0.10;
804:   tr->mu2 = 0.50;

806:   tr->gamma1 = 0.25;
807:   tr->gamma2 = 0.50;
808:   tr->gamma3 = 2.00;
809:   tr->gamma4 = 4.00;

811:   tr->theta = 0.05;

813:   tr->min_radius = 1.0e-10;
814:   tr->max_radius = 1.0e10;
815:   tr->epsilon = 1.0e-6;

817:   tr->ksp_type        = NTR_KSP_STCG;
818:   tr->pc_type         = NTR_PC_BFGS;
819:   tr->bfgs_scale_type = BFGS_SCALE_AHESS;
820:   tr->init_type              = NTR_INIT_INTERPOLATION;
821:   tr->update_type     = NTR_UPDATE_REDUCTION;
822:   TaoFunctionReturn(0);
823: }

826: #endif