Actual source code: ntl.c

  1: #include "src/matrix/lmvmmat.h"
  2: #include "ntl.h"

  4: #include "petscksp.h"
  5: #include "petscpc.h"
  6: #include "private/kspimpl.h"
  7: #include "private/pcimpl.h"

  9: #define NTL_KSP_NASH        0
 10: #define NTL_KSP_STCG        1
 11: #define NTL_KSP_GLTR        2
 12: #define NTL_KSP_TYPES        3

 14: #define NTL_PC_NONE        0
 15: #define NTL_PC_AHESS        1
 16: #define NTL_PC_BFGS        2
 17: #define NTL_PC_PETSC        3
 18: #define NTL_PC_TYPES        4

 20: #define BFGS_SCALE_AHESS        0
 21: #define BFGS_SCALE_BFGS                1
 22: #define BFGS_SCALE_TYPES        2

 24: #define NTL_INIT_CONSTANT         0
 25: #define NTL_INIT_DIRECTION        1
 26: #define NTL_INIT_INTERPOLATION    2
 27: #define NTL_INIT_TYPES            3

 29: #define NTL_UPDATE_REDUCTION      0
 30: #define NTL_UPDATE_INTERPOLATION  1
 31: #define NTL_UPDATE_TYPES          2

 33: static const char *NTL_KSP[64] = {
 34:   "nash", "stcg", "gltr"
 35: };

 37: static const char *NTL_PC[64] = {
 38:   "none", "ahess", "bfgs", "petsc"
 39: };

 41: static const char *BFGS_SCALE[64] = {
 42:   "ahess", "bfgs"
 43: };

 45: static const char *NTL_INIT[64] = {
 46:   "constant", "direction", "interpolation"
 47: };

 49: static const char *NTL_UPDATE[64] = {
 50:   "reduction", "interpolation"
 51: };

 53: /* Routine for BFGS preconditioner */

 57: static PetscErrorCode MatLMVMSolveShell(PC pc, Vec b, Vec x)
 58: {
 60:   Mat M;

 66:   PCShellGetContext(pc,(void**)&M); 
 67:   MatLMVMSolve(M, b, x); 
 68:   return(0);
 69: }

 71: /* Implements Newton's Method with a trust-region, line-search approach for 
 72:    solving unconstrained minimization problems.  A More'-Thuente line search 
 73:    is used to guarantee that the bfgs preconditioner remains positive
 74:    definite. */

 76: #define NTL_NEWTON                 0
 77: #define NTL_BFGS                 1
 78: #define NTL_SCALED_GRADIENT         2
 79: #define NTL_GRADIENT                 3

 83: static PetscErrorCode TaoSolve_NTL(TaoSolver tao)
 84: {
 85:   TAO_NTL *tl = (TAO_NTL *)tao->data;

 87:   PC pc;
 88:   KSPConvergedReason ksp_reason;
 89:   TaoSolverTerminationReason reason;
 90:   TaoLineSearchTerminationReason ls_reason;

 92:   PetscReal fmin, ftrial, prered, actred, kappa, sigma;
 93:   PetscReal tau, tau_1, tau_2, tau_max, tau_min, max_radius;
 94:   PetscReal f, fold, gdx, gnorm;
 95:   PetscReal step = 1.0;

 97:   PetscReal delta;
 98:   PetscReal norm_d = 0.0;
 99:   MatStructure matflag;
101:   PetscInt stepType;
102:   PetscInt iter = 0,its;

104:   PetscInt bfgsUpdates = 0;
105:   PetscInt needH;

107:   PetscInt i_max = 5;
108:   PetscInt j_max = 1;
109:   PetscInt i, j, n, N;

111:   PetscInt tr_reject;


115:   if (tao->XL || tao->XU || tao->ops->computebounds) {
116:     PetscPrintf(((PetscObject)tao)->comm,"WARNING: Variable bounds have been set but will be ignored by ntl algorithm\n"); 
117:   }

119:   /* Initialize trust-region radius */
120:   tao->trust = tao->trust0;

122:   /* Modify the radius if it is too large or small */
123:   tao->trust = PetscMax(tao->trust, tl->min_radius);
124:   tao->trust = PetscMin(tao->trust, tl->max_radius);

126:   if (NTL_PC_BFGS == tl->pc_type && !tl->M) {
127:     VecGetLocalSize(tao->solution,&n); 
128:     VecGetSize(tao->solution,&N); 
129:     MatCreateLMVM(((PetscObject)tao)->comm,n,N,&tl->M); 
130:     MatLMVMAllocateVectors(tl->M,tao->solution); 
131:   }

133:   /* Check convergence criteria */
134:   TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient); 
135:   VecNorm(tao->gradient, NORM_2, &gnorm); 
136:   if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) {
137:     SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");
138:   }
139:   needH = 1;

141:   TaoMonitor(tao, iter, f, gnorm, 0.0, 1.0, &reason); 
142:   if (reason != TAO_CONTINUE_ITERATING) {
143:     return(0);
144:   }

146:   /* Create vectors for the limited memory preconditioner */
147:   if ((NTL_PC_BFGS == tl->pc_type) && 
148:       (BFGS_SCALE_BFGS != tl->bfgs_scale_type)) {
149:     if (!tl->Diag) {
150:       VecDuplicate(tao->solution, &tl->Diag); 
151:     }
152:   }

154:   /* Modify the linear solver to a conjugate gradient method */
155:   switch(tl->ksp_type) {
156:   case NTL_KSP_NASH:
157:     KSPSetType(tao->ksp, KSPNASH); 
158:     if (tao->ksp->ops->setfromoptions) {
159:       (*tao->ksp->ops->setfromoptions)(tao->ksp);
160:     }
161:     break;

163:   case NTL_KSP_STCG:
164:     KSPSetType(tao->ksp, KSPSTCG); 
165:     if (tao->ksp->ops->setfromoptions) {
166:       (*tao->ksp->ops->setfromoptions)(tao->ksp);
167:     }
168:     break;

170:   default:
171:     KSPSetType(tao->ksp, KSPGLTR); 
172:     if (tao->ksp->ops->setfromoptions) {
173:       (*tao->ksp->ops->setfromoptions)(tao->ksp);
174:     }
175:     break;
176:   }

178:   /* Modify the preconditioner to use the bfgs approximation */
179:   KSPGetPC(tao->ksp, &pc); 
180:   switch(tl->pc_type) {
181:   case NTL_PC_NONE:
182:     PCSetType(pc, PCNONE); 
183:     if (pc->ops->setfromoptions) {
184:       (*pc->ops->setfromoptions)(pc);
185:     }
186:     break;

188:   case NTL_PC_AHESS:
189:     PCSetType(pc, PCJACOBI); 
190:     if (pc->ops->setfromoptions) {
191:       (*pc->ops->setfromoptions)(pc);
192:     }
193:     PCJacobiSetUseAbs(pc); 
194:     break;

196:   case NTL_PC_BFGS:
197:     PCSetType(pc, PCSHELL); 
198:     if (pc->ops->setfromoptions) {
199:       (*pc->ops->setfromoptions)(pc);
200:     }
201:     PCShellSetName(pc, "bfgs"); 
202:     PCShellSetContext(pc, tl->M); 
203:     PCShellSetApply(pc, MatLMVMSolveShell); 
204:     break;

206:   default:
207:     /* Use the pc method set by pc_type */
208:     break;
209:   }

211:   /* Initialize trust-region radius.  The initialization is only performed 
212:      when we are using Steihaug-Toint or the Generalized Lanczos method. */
213:   switch(tl->init_type) {
214:   case NTL_INIT_CONSTANT:
215:     /* Use the initial radius specified */
216:     break;

218:   case NTL_INIT_INTERPOLATION:
219:     /* Use the initial radius specified */
220:     max_radius = 0.0;
221:   
222:     for (j = 0; j < j_max; ++j) {
223:       fmin = f;
224:       sigma = 0.0;
225:   
226:       if (needH) {
227:         TaoComputeHessian(tao, tao->solution, &tao->hessian, &tao->hessian_pre, &matflag); 
228:         needH = 0;
229:       }
230:   
231:       for (i = 0; i < i_max; ++i) {
232:         VecCopy(tao->solution, tl->W); 
233:         VecAXPY(tl->W, -tao->trust/gnorm, tao->gradient); 

235:         TaoComputeObjective(tao, tl->W, &ftrial); 
236:         if (PetscIsInfOrNanReal(ftrial)) {
237:           tau = tl->gamma1_i;
238:         }
239:         else {
240:           if (ftrial < fmin) {
241:             fmin = ftrial;
242:             sigma = -tao->trust / gnorm;
243:           }

245:           MatMult(tao->hessian, tao->gradient, tao->stepdirection); 
246:           VecDot(tao->gradient, tao->stepdirection, &prered); 

248:           prered = tao->trust * (gnorm - 0.5 * tao->trust * prered / (gnorm * gnorm));
249:           actred = f - ftrial;
250:           if ((PetscAbsScalar(actred) <= tl->epsilon) && 
251:               (PetscAbsScalar(prered) <= tl->epsilon)) {
252:             kappa = 1.0;
253:           }
254:           else {
255:             kappa = actred / prered;
256:           }

258:           tau_1 = tl->theta_i * gnorm * tao->trust / (tl->theta_i * gnorm * tao->trust + (1.0 - tl->theta_i) * prered - actred);
259:           tau_2 = tl->theta_i * gnorm * tao->trust / (tl->theta_i * gnorm * tao->trust - (1.0 + tl->theta_i) * prered + actred);
260:           tau_min = PetscMin(tau_1, tau_2);
261:           tau_max = PetscMax(tau_1, tau_2);

263:           if (PetscAbsScalar(kappa - 1.0) <= tl->mu1_i) {
264:             /* Great agreement */
265:             max_radius = PetscMax(max_radius, tao->trust);

267:             if (tau_max < 1.0) {
268:               tau = tl->gamma3_i;
269:             }
270:             else if (tau_max > tl->gamma4_i) {
271:               tau = tl->gamma4_i;
272:             }
273:             else if (tau_1 >= 1.0 && tau_1 <= tl->gamma4_i && tau_2 < 1.0) {
274:               tau = tau_1;
275:             }
276:             else if (tau_2 >= 1.0 && tau_2 <= tl->gamma4_i && tau_1 < 1.0) {
277:               tau = tau_2;
278:             }
279:             else {
280:               tau = tau_max;
281:             }
282:           }
283:           else if (PetscAbsScalar(kappa - 1.0) <= tl->mu2_i) {
284:             /* Good agreement */
285:             max_radius = PetscMax(max_radius, tao->trust);

287:             if (tau_max < tl->gamma2_i) {
288:               tau = tl->gamma2_i;
289:             }
290:             else if (tau_max > tl->gamma3_i) {
291:               tau = tl->gamma3_i;
292:             }
293:             else {
294:               tau = tau_max;
295:             }
296:           }
297:           else {
298:             /* Not good agreement */
299:             if (tau_min > 1.0) {
300:               tau = tl->gamma2_i;
301:             }
302:             else if (tau_max < tl->gamma1_i) {
303:               tau = tl->gamma1_i;
304:             }
305:             else if ((tau_min < tl->gamma1_i) && (tau_max >= 1.0)) {
306:               tau = tl->gamma1_i;
307:             }
308:             else if ((tau_1 >= tl->gamma1_i) && (tau_1 < 1.0) &&
309:                      ((tau_2 < tl->gamma1_i) || (tau_2 >= 1.0))) {
310:               tau = tau_1;
311:             }
312:             else if ((tau_2 >= tl->gamma1_i) && (tau_2 < 1.0) &&
313:                      ((tau_1 < tl->gamma1_i) || (tau_2 >= 1.0))) {
314:               tau = tau_2;
315:             }
316:             else {
317:               tau = tau_max;
318:             }
319:           }
320:         }
321:         tao->trust = tau * tao->trust;
322:       }
323:   
324:       if (fmin < f) {
325:         f = fmin;
326:         VecAXPY(tao->solution, sigma, tao->gradient); 
327:         TaoComputeGradient(tao, tao->solution, tao->gradient); 

329:         VecNorm(tao->gradient, NORM_2, &gnorm); 
330:         if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) {
331:           SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");
332:         }
333:         needH = 1;
334:   
335:         TaoMonitor(tao, iter, f, gnorm, 0.0, 1.0, &reason); 
336:         if (reason != TAO_CONTINUE_ITERATING) {
337:           return(0);
338:         }
339:       }
340:     }
341:     tao->trust = PetscMax(tao->trust, max_radius);

343:     /* Modify the radius if it is too large or small */
344:     tao->trust = PetscMax(tao->trust, tl->min_radius);
345:     tao->trust = PetscMin(tao->trust, tl->max_radius);
346:     break;

348:   default:
349:     /* Norm of the first direction will initialize radius */
350:     tao->trust = 0.0;
351:     break;
352:   }

354:   /* Set initial scaling for the BFGS preconditioner
355:      This step is done after computing the initial trust-region radius
356:      since the function value may have decreased */
357:   if (NTL_PC_BFGS == tl->pc_type) {
358:     if (f != 0.0) {
359:       delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
360:     }
361:     else {
362:       delta = 2.0 / (gnorm*gnorm);
363:     }
364:     MatLMVMSetDelta(tl->M, delta); 
365:   }

367:   /* Set counter for gradient/reset steps */
368:   tl->ntrust = 0;
369:   tl->newt = 0;
370:   tl->bfgs = 0;
371:   tl->sgrad = 0;
372:   tl->grad = 0;

374:   /* Have not converged; continue with Newton method */
375:   while (reason == TAO_CONTINUE_ITERATING) {
376:     ++iter;

378:     /* Compute the Hessian */
379:     if (needH) {
380:       TaoComputeHessian(tao, tao->solution, &tao->hessian, &tao->hessian_pre, &matflag); 
381:       needH = 0;
382:     }

384:     if (NTL_PC_BFGS == tl->pc_type) {
385:       if (BFGS_SCALE_AHESS == tl->bfgs_scale_type) {
386:         /* Obtain diagonal for the bfgs preconditioner */
387:         MatGetDiagonal(tao->hessian, tl->Diag); 
388:         VecAbs(tl->Diag); 
389:         VecReciprocal(tl->Diag); 
390:         MatLMVMSetScale(tl->M, tl->Diag); 
391:       }

393:       /* Update the limited memory preconditioner */
394:       MatLMVMUpdate(tl->M,tao->solution, tao->gradient); 
395:       ++bfgsUpdates;
396:     }
397:     KSPSetOperators(tao->ksp, tao->hessian, tao->hessian_pre, matflag); 
398:     /* Solve the Newton system of equations */
399:     if (NTL_KSP_NASH == tl->ksp_type) {
400:       KSPNASHSetRadius(tao->ksp,tl->max_radius); 
401:       KSPSolve(tao->ksp, tao->gradient, tao->stepdirection); 
402:       KSPGetIterationNumber(tao->ksp,&its); 
403:       tao->ksp_its+=its;
404:       KSPNASHGetNormD(tao->ksp, &norm_d); 
405:     } else if (NTL_KSP_STCG == tl->ksp_type) {
406:       KSPSTCGSetRadius(tao->ksp,tl->max_radius); 
407:       KSPSolve(tao->ksp, tao->gradient, tao->stepdirection); 
408:       KSPGetIterationNumber(tao->ksp,&its); 
409:       tao->ksp_its+=its;
410:       KSPSTCGGetNormD(tao->ksp, &norm_d); 
411:     } else { /* NTL_KSP_GLTR */
412:       KSPGLTRSetRadius(tao->ksp,tl->max_radius); 
413:       KSPSolve(tao->ksp, tao->gradient, tao->stepdirection); 
414:       KSPGetIterationNumber(tao->ksp,&its); 
415:       tao->ksp_its+=its;
416:       KSPGLTRGetNormD(tao->ksp, &norm_d); 
417:     }

419:     if (0.0 == tao->trust) {
420:       /* Radius was uninitialized; use the norm of the direction */
421:       if (norm_d > 0.0) {
422:         tao->trust = norm_d;

424:         /* Modify the radius if it is too large or small */
425:         tao->trust = PetscMax(tao->trust, tl->min_radius);
426:         tao->trust = PetscMin(tao->trust, tl->max_radius);
427:       }
428:       else {
429:         /* The direction was bad; set radius to default value and re-solve 
430:            the trust-region subproblem to get a direction */
431:         tao->trust = tao->trust0;

433:         /* Modify the radius if it is too large or small */
434:         tao->trust = PetscMax(tao->trust, tl->min_radius);
435:         tao->trust = PetscMin(tao->trust, tl->max_radius);

437:         if (NTL_KSP_NASH == tl->ksp_type) {
438:           KSPNASHSetRadius(tao->ksp,tl->max_radius); 
439:           KSPSolve(tao->ksp, tao->gradient, tao->stepdirection); 
440:           KSPGetIterationNumber(tao->ksp,&its); 
441:           tao->ksp_its+=its;
442:           KSPNASHGetNormD(tao->ksp, &norm_d); 
443:         } else if (NTL_KSP_STCG == tl->ksp_type) {
444:           KSPSTCGSetRadius(tao->ksp,tl->max_radius); 
445:           KSPSolve(tao->ksp, tao->gradient, tao->stepdirection); 
446:           KSPGetIterationNumber(tao->ksp,&its); 
447:           tao->ksp_its+=its;
448:           KSPSTCGGetNormD(tao->ksp, &norm_d); 
449:         } else { /* NTL_KSP_GLTR */
450:           KSPGLTRSetRadius(tao->ksp,tl->max_radius); 
451:           KSPSolve(tao->ksp, tao->gradient, tao->stepdirection); 
452:           KSPGetIterationNumber(tao->ksp,&its); 
453:           tao->ksp_its+=its;
454:           KSPGLTRGetNormD(tao->ksp, &norm_d); 
455:         }


458:         if (norm_d == 0.0) {
459:           SETERRQ(PETSC_COMM_SELF,1, "Initial direction zero");
460:         }
461:       }
462:     }

464:     VecScale(tao->stepdirection, -1.0); 
465:     KSPGetConvergedReason(tao->ksp, &ksp_reason); 
466:     if ((KSP_DIVERGED_INDEFINITE_PC == ksp_reason) &&
467:         (NTL_PC_BFGS == tl->pc_type) && (bfgsUpdates > 1)) {
468:       /* Preconditioner is numerically indefinite; reset the
469:          approximate if using BFGS preconditioning. */

471:       if (f != 0.0) {
472:         delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
473:       }
474:       else {
475:         delta = 2.0 / (gnorm*gnorm);
476:       }
477:       MatLMVMSetDelta(tl->M, delta); 
478:       MatLMVMReset(tl->M); 
479:       MatLMVMUpdate(tl->M, tao->solution, tao->gradient); 
480:       bfgsUpdates = 1;
481:     }

483:     /* Check trust-region reduction conditions */
484:     tr_reject = 0;
485:     if (NTL_UPDATE_REDUCTION == tl->update_type) {
486:       /* Get predicted reduction */
487:       if (NTL_KSP_NASH == tl->ksp_type) {
488:         KSPNASHGetObjFcn(tao->ksp,&prered); 
489:       } else if (NTL_KSP_STCG == tl->ksp_type) {
490:         KSPSTCGGetObjFcn(tao->ksp,&prered); 
491:       } else { /* gltr */
492:         KSPGLTRGetObjFcn(tao->ksp,&prered); 
493:       }

495:       if (prered >= 0.0) {
496:         /* The predicted reduction has the wrong sign.  This cannot
497:            happen in infinite precision arithmetic.  Step should
498:            be rejected! */
499:         tao->trust = tl->alpha1 * PetscMin(tao->trust, norm_d);
500:         tr_reject = 1;
501:       }
502:       else {
503:         /* Compute trial step and function value */
504:         VecCopy(tao->solution, tl->W); 
505:         VecAXPY(tl->W, 1.0, tao->stepdirection); 
506:         TaoComputeObjective(tao, tl->W, &ftrial); 

508:         if (PetscIsInfOrNanReal(ftrial)) {
509:           tao->trust = tl->alpha1 * PetscMin(tao->trust, norm_d);
510:           tr_reject = 1;
511:         }
512:         else {
513:           /* Compute and actual reduction */
514:           actred = f - ftrial;
515:           prered = -prered;
516:           if ((PetscAbsScalar(actred) <= tl->epsilon) &&
517:               (PetscAbsScalar(prered) <= tl->epsilon)) {
518:             kappa = 1.0;
519:           }
520:           else {
521:             kappa = actred / prered;
522:           }

524:           /* Accept of reject the step and update radius */
525:           if (kappa < tl->eta1) {
526:             /* Reject the step */
527:             tao->trust = tl->alpha1 * PetscMin(tao->trust, norm_d);
528:             tr_reject = 1;
529:           }
530:           else {
531:             /* Accept the step */
532:             if (kappa < tl->eta2) {
533:               /* Marginal bad step */
534:               tao->trust = tl->alpha2 * PetscMin(tao->trust, norm_d);
535:             }
536:             else if (kappa < tl->eta3) {
537:               /* Reasonable step */
538:               tao->trust = tl->alpha3 * tao->trust;
539:             }
540:             else if (kappa < tl->eta4) {
541:               /* Good step */
542:               tao->trust = PetscMax(tl->alpha4 * norm_d, tao->trust);
543:             }
544:             else {
545:               /* Very good step */
546:               tao->trust = PetscMax(tl->alpha5 * norm_d, tao->trust);
547:             }
548:           }
549:         }
550:       }
551:     }
552:     else {
553:       /* Get predicted reduction */
554:       if (NTL_KSP_NASH == tl->ksp_type) {
555:         KSPNASHGetObjFcn(tao->ksp,&prered); 
556:       } else if (NTL_KSP_STCG == tl->ksp_type) {
557:         KSPSTCGGetObjFcn(tao->ksp,&prered); 
558:       } else { /* gltr */
559:         KSPGLTRGetObjFcn(tao->ksp,&prered); 
560:       }

562:       if (prered >= 0.0) {
563:         /* The predicted reduction has the wrong sign.  This cannot
564:            happen in infinite precision arithmetic.  Step should
565:            be rejected! */
566:         tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d);
567:         tr_reject = 1;
568:       }
569:       else {
570:         VecCopy(tao->solution, tl->W); 
571:         VecAXPY(tl->W, 1.0, tao->stepdirection); 
572:         TaoComputeObjective(tao, tl->W, &ftrial); 
573:         if (PetscIsInfOrNanReal(ftrial)) {
574:           tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d);
575:           tr_reject = 1;
576:         }
577:         else {
578:           VecDot(tao->gradient, tao->stepdirection, &gdx); 

580:           actred = f - ftrial;
581:           prered = -prered;
582:           if ((PetscAbsScalar(actred) <= tl->epsilon) &&
583:               (PetscAbsScalar(prered) <= tl->epsilon)) {
584:             kappa = 1.0;
585:           }
586:           else {
587:             kappa = actred / prered;
588:           }

590:           tau_1 = tl->theta * gdx / (tl->theta * gdx - (1.0 - tl->theta) * prered + actred);
591:           tau_2 = tl->theta * gdx / (tl->theta * gdx + (1.0 + tl->theta) * prered - actred);
592:           tau_min = PetscMin(tau_1, tau_2);
593:           tau_max = PetscMax(tau_1, tau_2);

595:           if (kappa >= 1.0 - tl->mu1) {
596:             /* Great agreement; accept step and update radius */
597:             if (tau_max < 1.0) {
598:               tao->trust = PetscMax(tao->trust, tl->gamma3 * norm_d);
599:             }
600:             else if (tau_max > tl->gamma4) {
601:               tao->trust = PetscMax(tao->trust, tl->gamma4 * norm_d);
602:             }
603:             else {
604:               tao->trust = PetscMax(tao->trust, tau_max * norm_d);
605:             }
606:           }
607:           else if (kappa >= 1.0 - tl->mu2) {
608:             /* Good agreement */

610:             if (tau_max < tl->gamma2) {
611:               tao->trust = tl->gamma2 * PetscMin(tao->trust, norm_d);
612:             }
613:             else if (tau_max > tl->gamma3) {
614:               tao->trust = PetscMax(tao->trust, tl->gamma3 * norm_d);
615:             }              else if (tau_max < 1.0) {
616:               tao->trust = tau_max * PetscMin(tao->trust, norm_d);
617:             }
618:             else {
619:               tao->trust = PetscMax(tao->trust, tau_max * norm_d);
620:             }
621:           }
622:           else {
623:             /* Not good agreement */
624:             if (tau_min > 1.0) {
625:               tao->trust = tl->gamma2 * PetscMin(tao->trust, norm_d);
626:             }
627:             else if (tau_max < tl->gamma1) {
628:               tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d);
629:             }
630:             else if ((tau_min < tl->gamma1) && (tau_max >= 1.0)) {
631:               tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d);
632:             }
633:             else if ((tau_1 >= tl->gamma1) && (tau_1 < 1.0) &&
634:                      ((tau_2 < tl->gamma1) || (tau_2 >= 1.0))) {
635:               tao->trust = tau_1 * PetscMin(tao->trust, norm_d);
636:             }
637:             else if ((tau_2 >= tl->gamma1) && (tau_2 < 1.0) &&
638:                      ((tau_1 < tl->gamma1) || (tau_2 >= 1.0))) {
639:               tao->trust = tau_2 * PetscMin(tao->trust, norm_d);
640:             }
641:             else {
642:               tao->trust = tau_max * PetscMin(tao->trust, norm_d);
643:             }
644:             tr_reject = 1;
645:           }
646:         }
647:       }
648:     }

650:     if (tr_reject) {
651:       /* The trust-region constraints rejected the step.  Apply a linesearch.
652:          Check for descent direction. */
653:       VecDot(tao->stepdirection, tao->gradient, &gdx); 
654:       if ((gdx >= 0.0) || PetscIsInfOrNanReal(gdx)) {
655:         /* Newton step is not descent or direction produced Inf or NaN */
656:         
657:         if (NTL_PC_BFGS != tl->pc_type) {
658:           /* We don't have the bfgs matrix around and updated
659:              Must use gradient direction in this case */
660:           VecCopy(tao->gradient, tao->stepdirection); 
661:           VecScale(tao->stepdirection, -1.0); 
662:           ++tl->grad;
663:           stepType = NTL_GRADIENT;
664:         }
665:         else {
666:           /* Attempt to use the BFGS direction */
667:           MatLMVMSolve(tl->M, tao->gradient, tao->stepdirection); 
668:           VecScale(tao->stepdirection, -1.0); 
669:           
670:           /* Check for success (descent direction) */
671:           VecDot(tao->stepdirection, tao->gradient, &gdx); 
672:           if ((gdx >= 0) || PetscIsInfOrNanReal(gdx)) {
673:             /* BFGS direction is not descent or direction produced not a number
674:                We can assert bfgsUpdates > 1 in this case because
675:                the first solve produces the scaled gradient direction,
676:                which is guaranteed to be descent */
677:             
678:             /* Use steepest descent direction (scaled) */
679:             if (f != 0.0) {
680:               delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
681:             }
682:             else {
683:               delta = 2.0 / (gnorm*gnorm);
684:             }
685:             MatLMVMSetDelta(tl->M, delta); 
686:             MatLMVMReset(tl->M); 
687:             MatLMVMUpdate(tl->M, tao->solution, tao->gradient); 
688:             MatLMVMSolve(tl->M, tao->gradient, tao->stepdirection); 
689:             VecScale(tao->stepdirection, -1.0); 
690:             
691:             bfgsUpdates = 1;
692:             ++tl->sgrad;
693:             stepType = NTL_SCALED_GRADIENT;
694:           }
695:           else {
696:             if (1 == bfgsUpdates) {
697:               /* The first BFGS direction is always the scaled gradient */
698:               ++tl->sgrad;
699:               stepType = NTL_SCALED_GRADIENT;
700:             }
701:             else {
702:               ++tl->bfgs;
703:               stepType = NTL_BFGS;
704:             }
705:           }
706:         }
707:       }
708:       else {
709:         /* Computed Newton step is descent */
710:         ++tl->newt;
711:         stepType = NTL_NEWTON;
712:       }
713:       
714:       /* Perform the linesearch */
715:       fold = f;
716:       VecCopy(tao->solution, tl->Xold); 
717:       VecCopy(tao->gradient, tl->Gold); 

719:       step = 1.0;
720:       TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_reason); 
721:       TaoAddLineSearchCounts(tao); 

723:       
724:       while (ls_reason < 0 && stepType != NTL_GRADIENT) {
725:         /* Linesearch failed */
726:         f = fold;
727:         VecCopy(tl->Xold, tao->solution); 
728:         VecCopy(tl->Gold, tao->gradient); 
729:         
730:         switch(stepType) {
731:         case NTL_NEWTON:
732:           /* Failed to obtain acceptable iterate with Newton step */

734:           if (NTL_PC_BFGS != tl->pc_type) {
735:             /* We don't have the bfgs matrix around and being updated
736:                Must use gradient direction in this case */
737:             VecCopy(tao->gradient, tao->stepdirection); 
738:             ++tl->grad;
739:             stepType = NTL_GRADIENT;
740:           }
741:           else {
742:             /* Attempt to use the BFGS direction */
743:             MatLMVMSolve(tl->M, tao->gradient, tao->stepdirection); 

745:             
746:             /* Check for success (descent direction) */
747:             VecDot(tao->stepdirection, tao->gradient, &gdx); 
748:             if ((gdx <= 0) || PetscIsInfOrNanReal(gdx)) {
749:               /* BFGS direction is not descent or direction produced 
750:                  not a number.  We can assert bfgsUpdates > 1 in this case
751:                  Use steepest descent direction (scaled) */
752:     
753:               if (f != 0.0) {
754:                 delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
755:               }
756:               else {
757:                 delta = 2.0 / (gnorm*gnorm);
758:               }
759:               MatLMVMSetDelta(tl->M, delta); 
760:               MatLMVMReset(tl->M); 
761:               MatLMVMUpdate(tl->M, tao->solution, tao->gradient); 
762:               MatLMVMSolve(tl->M, tao->gradient, tao->stepdirection); 
763:               
764:               bfgsUpdates = 1;
765:               ++tl->sgrad;
766:               stepType = NTL_SCALED_GRADIENT;
767:             }
768:             else {
769:               if (1 == bfgsUpdates) {
770:                 /* The first BFGS direction is always the scaled gradient */
771:                 ++tl->sgrad;
772:                 stepType = NTL_SCALED_GRADIENT;
773:               }
774:               else {
775:                 ++tl->bfgs;
776:                 stepType = NTL_BFGS;
777:               }
778:             }
779:           }
780:           break;

782:         case NTL_BFGS:
783:           /* Can only enter if pc_type == NTL_PC_BFGS
784:              Failed to obtain acceptable iterate with BFGS step
785:              Attempt to use the scaled gradient direction */
786:           
787:           if (f != 0.0) {
788:             delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
789:           }
790:           else {
791:             delta = 2.0 / (gnorm*gnorm);
792:           }
793:           MatLMVMSetDelta(tl->M, delta); 
794:           MatLMVMReset(tl->M); 
795:           MatLMVMUpdate(tl->M, tao->solution, tao->gradient); 
796:           MatLMVMSolve(tl->M, tao->gradient, tao->stepdirection); 

798:           bfgsUpdates = 1;
799:           ++tl->sgrad;
800:           stepType = NTL_SCALED_GRADIENT;
801:           break;
802:           
803:         case NTL_SCALED_GRADIENT:
804:           /* Can only enter if pc_type == NTL_PC_BFGS
805:              The scaled gradient step did not produce a new iterate;
806:              attemp to use the gradient direction.
807:              Need to make sure we are not using a different diagonal scaling */
808:           MatLMVMSetScale(tl->M, tl->Diag); 
809:           MatLMVMSetDelta(tl->M, 1.0); 
810:           MatLMVMReset(tl->M); 
811:           MatLMVMUpdate(tl->M, tao->solution, tao->gradient); 
812:           MatLMVMSolve(tl->M, tao->gradient, tao->stepdirection); 
813:           
814:           bfgsUpdates = 1;
815:           ++tl->grad;
816:           stepType = NTL_GRADIENT;
817:           break;
818:         }
819:         VecScale(tao->stepdirection, -1.0); 

821:         /* This may be incorrect; linesearch has values for stepmax and stepmin
822:            that should be reset. */
823:         step = 1.0;
824:         TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_reason); 
825:         TaoAddLineSearchCounts(tao); 
826:       }

828:       if (ls_reason < 0) {
829:         /* Failed to find an improving point */
830:         f = fold;
831:         VecCopy(tl->Xold, tao->solution); 
832:         VecCopy(tl->Gold, tao->gradient); 
833:         tao->trust = 0.0;
834:       }
835:       else if (stepType == NTL_NEWTON) {
836:         if (step < tl->nu1) {
837:           /* Very bad step taken; reduce radius */
838:           tao->trust = tl->omega1 * PetscMin(norm_d, tao->trust);
839:         }
840:         else if (step < tl->nu2) {
841:           /* Reasonably bad step taken; reduce radius */
842:           tao->trust = tl->omega2 * PetscMin(norm_d, tao->trust);
843:         }
844:         else if (step < tl->nu3) {
845:           /* Reasonable step was taken; leave radius alone */
846:           if (tl->omega3 < 1.0) {
847:             tao->trust = tl->omega3 * PetscMin(norm_d, tao->trust);
848:           }
849:           else if (tl->omega3 > 1.0) {
850:             tao->trust = PetscMax(tl->omega3 * norm_d, tao->trust);
851:           }
852:         }
853:         else if (step < tl->nu4) {
854:           /* Full step taken; increase the radius */
855:           tao->trust = PetscMax(tl->omega4 * norm_d, tao->trust);
856:         }
857:         else {
858:           /* More than full step taken; increase the radius */
859:           tao->trust = PetscMax(tl->omega5 * norm_d, tao->trust);
860:         }
861:       }
862:       else {
863:         /* Newton step was not good; reduce the radius */
864:         tao->trust = tl->omega1 * PetscMin(norm_d, tao->trust);
865:       }
866:     }
867:     else {
868:       /* Trust-region step is accepted */
869:       VecCopy(tl->W, tao->solution); 
870:       f = ftrial;
871:       TaoComputeGradient(tao, tao->solution, tao->gradient); 
872:       ++tl->ntrust;
873:     }

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

878:     /* Check for termination */
879:     VecNorm(tao->gradient, NORM_2, &gnorm); 
880:     if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) {
881:       SETERRQ(PETSC_COMM_SELF,1,"User provided compute function generated Not-a-Number");
882:     }
883:     needH = 1;
884:     
885:     TaoMonitor(tao, iter, f, gnorm, 0.0, tao->trust, &reason); 
886:   }
887:   return(0);
888: }

890: /* ---------------------------------------------------------- */
893: static PetscErrorCode TaoSetUp_NTL(TaoSolver tao)
894: {
895:   TAO_NTL *tl = (TAO_NTL *)tao->data;

899:   if (!tao->gradient) {VecDuplicate(tao->solution, &tao->gradient);  }
900:   if (!tao->stepdirection) {VecDuplicate(tao->solution, &tao->stepdirection); }
901:   if (!tl->W) { VecDuplicate(tao->solution, &tl->W); }
902:   if (!tl->Xold) { VecDuplicate(tao->solution, &tl->Xold); }
903:   if (!tl->Gold) { VecDuplicate(tao->solution, &tl->Gold); }

905:   tl->Diag = 0;
906:   tl->M = 0;


909:   return(0);
910: }

912: /*------------------------------------------------------------*/
915: static PetscErrorCode TaoDestroy_NTL(TaoSolver tao)
916: {
917:   TAO_NTL *tl = (TAO_NTL *)tao->data;

921:   if (tao->setupcalled) {
922:     VecDestroy(&tl->W); 
923:     VecDestroy(&tl->Xold); 
924:     VecDestroy(&tl->Gold); 
925:   }
926:   if (tl->Diag) {
927:     VecDestroy(&tl->Diag); 
928:     tl->Diag = PETSC_NULL;
929:   }
930:   if (tl->M) {
931:     MatDestroy(&tl->M); 
932:     tl->M = PETSC_NULL;
933:   }

935:   PetscFree(tao->data); 
936:   tao->data = PETSC_NULL;
937:   
938:   return(0);
939: }

941: /*------------------------------------------------------------*/
944: static PetscErrorCode TaoSetFromOptions_NTL(TaoSolver tao)
945: {
946:   TAO_NTL *tl = (TAO_NTL *)tao->data;

950:   PetscOptionsHead("Newton line search method for unconstrained optimization"); 
951:   PetscOptionsEList("-tao_ntl_ksp_type", "ksp type", "", NTL_KSP, NTL_KSP_TYPES, NTL_KSP[tl->ksp_type], &tl->ksp_type, 0); 
952:   PetscOptionsEList("-tao_ntl_pc_type", "pc type", "", NTL_PC, NTL_PC_TYPES, NTL_PC[tl->pc_type], &tl->pc_type, 0); 
953:   PetscOptionsEList("-tao_ntl_bfgs_scale_type", "bfgs scale type", "", BFGS_SCALE, BFGS_SCALE_TYPES, BFGS_SCALE[tl->bfgs_scale_type], &tl->bfgs_scale_type, 0); 
954:   PetscOptionsEList("-tao_ntl_init_type", "radius initialization type", "", NTL_INIT, NTL_INIT_TYPES, NTL_INIT[tl->init_type], &tl->init_type, 0); 
955:   PetscOptionsEList("-tao_ntl_update_type", "radius update type", "", NTL_UPDATE, NTL_UPDATE_TYPES, NTL_UPDATE[tl->update_type], &tl->update_type, 0); 
956:   PetscOptionsReal("-tao_ntl_eta1", "poor steplength; reduce radius", "", tl->eta1, &tl->eta1, 0); 
957:   PetscOptionsReal("-tao_ntl_eta2", "reasonable steplength; leave radius alone", "", tl->eta2, &tl->eta2, 0); 
958:   PetscOptionsReal("-tao_ntl_eta3", "good steplength; increase radius", "", tl->eta3, &tl->eta3, 0); 
959:   PetscOptionsReal("-tao_ntl_eta4", "excellent steplength; greatly increase radius", "", tl->eta4, &tl->eta4, 0); 
960:   PetscOptionsReal("-tao_ntl_alpha1", "", "", tl->alpha1, &tl->alpha1, 0); 
961:   PetscOptionsReal("-tao_ntl_alpha2", "", "", tl->alpha2, &tl->alpha2, 0); 
962:   PetscOptionsReal("-tao_ntl_alpha3", "", "", tl->alpha3, &tl->alpha3, 0); 
963:   PetscOptionsReal("-tao_ntl_alpha4", "", "", tl->alpha4, &tl->alpha4, 0); 
964:   PetscOptionsReal("-tao_ntl_alpha5", "", "", tl->alpha5, &tl->alpha5, 0); 
965:   PetscOptionsReal("-tao_ntl_nu1", "poor steplength; reduce radius", "", tl->nu1, &tl->nu1, 0); 
966:   PetscOptionsReal("-tao_ntl_nu2", "reasonable steplength; leave radius alone", "", tl->nu2, &tl->nu2, 0); 
967:   PetscOptionsReal("-tao_ntl_nu3", "good steplength; increase radius", "", tl->nu3, &tl->nu3, 0); 
968:   PetscOptionsReal("-tao_ntl_nu4", "excellent steplength; greatly increase radius", "", tl->nu4, &tl->nu4, 0); 
969:   PetscOptionsReal("-tao_ntl_omega1", "", "", tl->omega1, &tl->omega1, 0); 
970:   PetscOptionsReal("-tao_ntl_omega2", "", "", tl->omega2, &tl->omega2, 0); 
971:   PetscOptionsReal("-tao_ntl_omega3", "", "", tl->omega3, &tl->omega3, 0); 
972:   PetscOptionsReal("-tao_ntl_omega4", "", "", tl->omega4, &tl->omega4, 0); 
973:   PetscOptionsReal("-tao_ntl_omega5", "", "", tl->omega5, &tl->omega5, 0); 
974:   PetscOptionsReal("-tao_ntl_mu1_i", "", "", tl->mu1_i, &tl->mu1_i, 0); 
975:   PetscOptionsReal("-tao_ntl_mu2_i", "", "", tl->mu2_i, &tl->mu2_i, 0); 
976:   PetscOptionsReal("-tao_ntl_gamma1_i", "", "", tl->gamma1_i, &tl->gamma1_i, 0); 
977:   PetscOptionsReal("-tao_ntl_gamma2_i", "", "", tl->gamma2_i, &tl->gamma2_i, 0); 
978:   PetscOptionsReal("-tao_ntl_gamma3_i", "", "", tl->gamma3_i, &tl->gamma3_i, 0); 
979:   PetscOptionsReal("-tao_ntl_gamma4_i", "", "", tl->gamma4_i, &tl->gamma4_i, 0); 
980:   PetscOptionsReal("-tao_ntl_theta_i", "", "", tl->theta_i, &tl->theta_i, 0); 
981:   PetscOptionsReal("-tao_ntl_mu1", "", "", tl->mu1, &tl->mu1, 0); 
982:   PetscOptionsReal("-tao_ntl_mu2", "", "", tl->mu2, &tl->mu2, 0); 
983:   PetscOptionsReal("-tao_ntl_gamma1", "", "", tl->gamma1, &tl->gamma1, 0); 
984:   PetscOptionsReal("-tao_ntl_gamma2", "", "", tl->gamma2, &tl->gamma2, 0); 
985:   PetscOptionsReal("-tao_ntl_gamma3", "", "", tl->gamma3, &tl->gamma3, 0); 
986:   PetscOptionsReal("-tao_ntl_gamma4", "", "", tl->gamma4, &tl->gamma4, 0); 
987:   PetscOptionsReal("-tao_ntl_theta", "", "", tl->theta, &tl->theta, 0); 
988:   PetscOptionsReal("-tao_ntl_min_radius", "lower bound on initial radius", "", tl->min_radius, &tl->min_radius, 0); 
989:   PetscOptionsReal("-tao_ntl_max_radius", "upper bound on radius", "", tl->max_radius, &tl->max_radius, 0); 
990:   PetscOptionsReal("-tao_ntl_epsilon", "tolerance used when computing actual and predicted reduction", "", tl->epsilon, &tl->epsilon, 0); 
991:   PetscOptionsTail(); 
992:   TaoLineSearchSetFromOptions(tao->linesearch); 
993:   KSPSetFromOptions(tao->ksp); 
994:   return(0);
995: }


998: /*------------------------------------------------------------*/
1001: static PetscErrorCode TaoView_NTL(TaoSolver tao, PetscViewer viewer)
1002: {
1003:   TAO_NTL *tl = (TAO_NTL *)tao->data;
1004:   PetscInt nrejects;
1005:   PetscBool isascii;

1009:   PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
1010:   if (isascii) {
1011:     PetscViewerASCIIPushTab(viewer); 
1012:     if (NTL_PC_BFGS == tl->pc_type && tl->M) {
1013:       MatLMVMGetRejects(tl->M, &nrejects); 
1014:       PetscViewerASCIIPrintf(viewer, "Rejected matrix updates: %D\n", &nrejects); 

1016:     }

1018:     PetscViewerASCIIPrintf(viewer, "Trust-region steps: %D\n", tl->ntrust); 
1019:     PetscViewerASCIIPrintf(viewer, "Newton search steps: %D\n", tl->newt); 
1020:     PetscViewerASCIIPrintf(viewer, "BFGS search steps: %D\n", tl->bfgs); 
1021:     PetscViewerASCIIPrintf(viewer, "Scaled gradient search steps: %D\n", tl->sgrad); 
1022:     PetscViewerASCIIPrintf(viewer, "Gradient search steps: %D\n", tl->grad); 
1023:     PetscViewerASCIIPopTab(viewer); 
1024:   } else {
1025:     SETERRQ1(((PetscObject)tao)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for TAO NTL",((PetscObject)viewer)->type_name);
1026:   }
1027:   return(0);
1028: }

1030: /* ---------------------------------------------------------- */
1034: PetscErrorCode TaoCreate_NTL(TaoSolver tao)
1035: {
1036:   TAO_NTL *tl;
1038:   const char *morethuente_type = TAOLINESEARCH_MT;
1040:   PetscNewLog(tao, TAO_NTL, &tl); 
1041:   
1042:   tao->ops->setup = TaoSetUp_NTL;
1043:   tao->ops->solve = TaoSolve_NTL;
1044:   tao->ops->view = TaoView_NTL;
1045:   tao->ops->setfromoptions = TaoSetFromOptions_NTL;
1046:   tao->ops->destroy = TaoDestroy_NTL;
1047:   
1048:   tao->max_it = 50;
1049:   tao->fatol = 1e-10;
1050:   tao->frtol = 1e-10;
1051:   tao->data = (void*)tl;
1052:   
1053:   tao->trust0 = 100.0;


1056:   /* Default values for trust-region radius update based on steplength */
1057:   tl->nu1 = 0.25;
1058:   tl->nu2 = 0.50;
1059:   tl->nu3 = 1.00;
1060:   tl->nu4 = 1.25;

1062:   tl->omega1 = 0.25;
1063:   tl->omega2 = 0.50;
1064:   tl->omega3 = 1.00;
1065:   tl->omega4 = 2.00;
1066:   tl->omega5 = 4.00;

1068:   /* Default values for trust-region radius update based on reduction */
1069:   tl->eta1 = 1.0e-4;
1070:   tl->eta2 = 0.25;
1071:   tl->eta3 = 0.50;
1072:   tl->eta4 = 0.90;

1074:   tl->alpha1 = 0.25;
1075:   tl->alpha2 = 0.50;
1076:   tl->alpha3 = 1.00;
1077:   tl->alpha4 = 2.00;
1078:   tl->alpha5 = 4.00;

1080:   /* Default values for trust-region radius update based on interpolation */
1081:   tl->mu1 = 0.10;
1082:   tl->mu2 = 0.50;

1084:   tl->gamma1 = 0.25;
1085:   tl->gamma2 = 0.50;
1086:   tl->gamma3 = 2.00;
1087:   tl->gamma4 = 4.00;

1089:   tl->theta = 0.05;

1091:   /* Default values for trust region initialization based on interpolation */
1092:   tl->mu1_i = 0.35;
1093:   tl->mu2_i = 0.50;

1095:   tl->gamma1_i = 0.0625;
1096:   tl->gamma2_i = 0.5;
1097:   tl->gamma3_i = 2.0;
1098:   tl->gamma4_i = 5.0;
1099:   
1100:   tl->theta_i = 0.25;

1102:   /* Remaining parameters */
1103:   tl->min_radius = 1.0e-10;
1104:   tl->max_radius = 1.0e10;
1105:   tl->epsilon = 1.0e-6;

1107:   tl->ksp_type        = NTL_KSP_STCG;
1108:   tl->pc_type         = NTL_PC_BFGS;
1109:   tl->bfgs_scale_type = BFGS_SCALE_AHESS;
1110:   tl->init_type       = NTL_INIT_INTERPOLATION;
1111:   tl->update_type     = NTL_UPDATE_REDUCTION;

1113:   TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch); 
1114:   TaoLineSearchSetType(tao->linesearch, morethuente_type); 
1115:   TaoLineSearchUseTaoSolverRoutines(tao->linesearch, tao); 

1117:   KSPCreate(((PetscObject)tao)->comm, &tao->ksp); 

1119:   return(0);
1120: }