Actual source code: ntl.c

      1: #include <../src/tao/unconstrained/impls/ntl/ntlimpl.h>

      3: #include <petscksp.h>

      5: #define NTL_INIT_CONSTANT      0
      6: #define NTL_INIT_DIRECTION     1
      7: #define NTL_INIT_INTERPOLATION 2
      8: #define NTL_INIT_TYPES         3

     10: #define NTL_UPDATE_REDUCTION     0
     11: #define NTL_UPDATE_INTERPOLATION 1
     12: #define NTL_UPDATE_TYPES         2

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

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

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

     23: #define NTL_NEWTON          0
     24: #define NTL_BFGS            1
     25: #define NTL_SCALED_GRADIENT 2
     26: #define NTL_GRADIENT        3

     28: static PetscErrorCode TaoSolve_NTL(Tao tao)
     29: {
     30:   TAO_NTL                     *tl = (TAO_NTL *)tao->data;
     31:   KSPType                      ksp_type;
     32:   PetscBool                    is_nash, is_stcg, is_gltr, is_bfgs, is_jacobi, is_symmetric, sym_set;
     33:   KSPConvergedReason           ksp_reason;
     34:   PC                           pc;
     35:   TaoLineSearchConvergedReason ls_reason;

     37:   PetscReal fmin, ftrial, prered, actred, kappa, sigma;
     38:   PetscReal tau, tau_1, tau_2, tau_max, tau_min, max_radius;
     39:   PetscReal f, fold, gdx, gnorm;
     40:   PetscReal step = 1.0;

     42:   PetscReal norm_d = 0.0;
     43:   PetscInt  stepType;
     44:   PetscInt  its;

     46:   PetscInt bfgsUpdates = 0;
     47:   PetscInt needH;

     49:   PetscInt i_max = 5;
     50:   PetscInt j_max = 1;
     51:   PetscInt i, j, n, N;

     53:   PetscInt tr_reject;

     55:   PetscFunctionBegin;
     56:   if (tao->XL || tao->XU || tao->ops->computebounds) PetscCall(PetscInfo(tao, "WARNING: Variable bounds have been set but will be ignored by ntl algorithm\n"));

     58:   PetscCall(KSPGetType(tao->ksp, &ksp_type));
     59:   PetscCall(PetscStrcmp(ksp_type, KSPNASH, &is_nash));
     60:   PetscCall(PetscStrcmp(ksp_type, KSPSTCG, &is_stcg));
     61:   PetscCall(PetscStrcmp(ksp_type, KSPGLTR, &is_gltr));
     62:   PetscCheck(is_nash || is_stcg || is_gltr, PetscObjectComm((PetscObject)tao), PETSC_ERR_SUP, "TAO_NTR requires nash, stcg, or gltr for the KSP");

     64:   /* Initialize the radius and modify if it is too large or small */
     65:   tao->trust = tao->trust0;
     66:   tao->trust = PetscMax(tao->trust, tl->min_radius);
     67:   tao->trust = PetscMin(tao->trust, tl->max_radius);

     69:   /* Allocate the vectors needed for the BFGS approximation */
     70:   PetscCall(KSPGetPC(tao->ksp, &pc));
     71:   PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &is_bfgs));
     72:   PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCJACOBI, &is_jacobi));
     73:   if (is_bfgs) {
     74:     tl->bfgs_pre = pc;
     75:     PetscCall(PCLMVMGetMatLMVM(tl->bfgs_pre, &tl->M));
     76:     PetscCall(VecGetLocalSize(tao->solution, &n));
     77:     PetscCall(VecGetSize(tao->solution, &N));
     78:     PetscCall(MatSetSizes(tl->M, n, n, N, N));
     79:     PetscCall(MatLMVMAllocate(tl->M, tao->solution, tao->gradient));
     80:     PetscCall(MatIsSymmetricKnown(tl->M, &sym_set, &is_symmetric));
     81:     PetscCheck(sym_set && is_symmetric, PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM matrix in the LMVM preconditioner must be symmetric.");
     82:   } else if (is_jacobi) PetscCall(PCJacobiSetUseAbs(pc, PETSC_TRUE));

     84:   /* Check convergence criteria */
     85:   PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient));
     86:   PetscCall(VecNorm(tao->gradient, NORM_2, &gnorm));
     87:   PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated Inf or NaN");
     88:   needH = 1;

     90:   tao->reason = TAO_CONTINUE_ITERATING;
     91:   PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its));
     92:   PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step));
     93:   PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
     94:   if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(PETSC_SUCCESS);

     96:   /* Initialize trust-region radius */
     97:   switch (tl->init_type) {
     98:   case NTL_INIT_CONSTANT:
     99:     /* Use the initial radius specified */
    100:     break;

    102:   case NTL_INIT_INTERPOLATION:
    103:     /* Use the initial radius specified */
    104:     max_radius = 0.0;

    106:     for (j = 0; j < j_max; ++j) {
    107:       fmin  = f;
    108:       sigma = 0.0;

    110:       if (needH) {
    111:         PetscCall(TaoComputeHessian(tao, tao->solution, tao->hessian, tao->hessian_pre));
    112:         needH = 0;
    113:       }

    115:       for (i = 0; i < i_max; ++i) {
    116:         PetscCall(VecCopy(tao->solution, tl->W));
    117:         PetscCall(VecAXPY(tl->W, -tao->trust / gnorm, tao->gradient));

    119:         PetscCall(TaoComputeObjective(tao, tl->W, &ftrial));
    120:         if (PetscIsInfOrNanReal(ftrial)) {
    121:           tau = tl->gamma1_i;
    122:         } else {
    123:           if (ftrial < fmin) {
    124:             fmin  = ftrial;
    125:             sigma = -tao->trust / gnorm;
    126:           }

    128:           PetscCall(MatMult(tao->hessian, tao->gradient, tao->stepdirection));
    129:           PetscCall(VecDot(tao->gradient, tao->stepdirection, &prered));

    131:           prered = tao->trust * (gnorm - 0.5 * tao->trust * prered / (gnorm * gnorm));
    132:           actred = f - ftrial;
    133:           if ((PetscAbsScalar(actred) <= tl->epsilon) && (PetscAbsScalar(prered) <= tl->epsilon)) {
    134:             kappa = 1.0;
    135:           } else {
    136:             kappa = actred / prered;
    137:           }

    139:           tau_1   = tl->theta_i * gnorm * tao->trust / (tl->theta_i * gnorm * tao->trust + (1.0 - tl->theta_i) * prered - actred);
    140:           tau_2   = tl->theta_i * gnorm * tao->trust / (tl->theta_i * gnorm * tao->trust - (1.0 + tl->theta_i) * prered + actred);
    141:           tau_min = PetscMin(tau_1, tau_2);
    142:           tau_max = PetscMax(tau_1, tau_2);

    144:           if (PetscAbsScalar(kappa - (PetscReal)1.0) <= tl->mu1_i) {
    145:             /* Great agreement */
    146:             max_radius = PetscMax(max_radius, tao->trust);

    148:             if (tau_max < 1.0) {
    149:               tau = tl->gamma3_i;
    150:             } else if (tau_max > tl->gamma4_i) {
    151:               tau = tl->gamma4_i;
    152:             } else if (tau_1 >= 1.0 && tau_1 <= tl->gamma4_i && tau_2 < 1.0) {
    153:               tau = tau_1;
    154:             } else if (tau_2 >= 1.0 && tau_2 <= tl->gamma4_i && tau_1 < 1.0) {
    155:               tau = tau_2;
    156:             } else {
    157:               tau = tau_max;
    158:             }
    159:           } else if (PetscAbsScalar(kappa - (PetscReal)1.0) <= tl->mu2_i) {
    160:             /* Good agreement */
    161:             max_radius = PetscMax(max_radius, tao->trust);

    163:             if (tau_max < tl->gamma2_i) {
    164:               tau = tl->gamma2_i;
    165:             } else if (tau_max > tl->gamma3_i) {
    166:               tau = tl->gamma3_i;
    167:             } else {
    168:               tau = tau_max;
    169:             }
    170:           } else {
    171:             /* Not good agreement */
    172:             if (tau_min > 1.0) {
    173:               tau = tl->gamma2_i;
    174:             } else if (tau_max < tl->gamma1_i) {
    175:               tau = tl->gamma1_i;
    176:             } else if ((tau_min < tl->gamma1_i) && (tau_max >= 1.0)) {
    177:               tau = tl->gamma1_i;
    178:             } else if ((tau_1 >= tl->gamma1_i) && (tau_1 < 1.0) && ((tau_2 < tl->gamma1_i) || (tau_2 >= 1.0))) {
    179:               tau = tau_1;
    180:             } else if ((tau_2 >= tl->gamma1_i) && (tau_2 < 1.0) && ((tau_1 < tl->gamma1_i) || (tau_2 >= 1.0))) {
    181:               tau = tau_2;
    182:             } else {
    183:               tau = tau_max;
    184:             }
    185:           }
    186:         }
    187:         tao->trust = tau * tao->trust;
    188:       }

    190:       if (fmin < f) {
    191:         f = fmin;
    192:         PetscCall(VecAXPY(tao->solution, sigma, tao->gradient));
    193:         PetscCall(TaoComputeGradient(tao, tao->solution, tao->gradient));

    195:         PetscCall(VecNorm(tao->gradient, NORM_2, &gnorm));
    196:         PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated Inf or NaN");
    197:         needH = 1;

    199:         PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its));
    200:         PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step));
    201:         PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
    202:         if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(PETSC_SUCCESS);
    203:       }
    204:     }
    205:     tao->trust = PetscMax(tao->trust, max_radius);

    207:     /* Modify the radius if it is too large or small */
    208:     tao->trust = PetscMax(tao->trust, tl->min_radius);
    209:     tao->trust = PetscMin(tao->trust, tl->max_radius);
    210:     break;

    212:   default:
    213:     /* Norm of the first direction will initialize radius */
    214:     tao->trust = 0.0;
    215:     break;
    216:   }

    218:   /* Set counter for gradient/reset steps */
    219:   tl->ntrust = 0;
    220:   tl->newt   = 0;
    221:   tl->bfgs   = 0;
    222:   tl->grad   = 0;

    224:   /* Have not converged; continue with Newton method */
    225:   while (tao->reason == TAO_CONTINUE_ITERATING) {
    226:     /* Call general purpose update function */
    227:     if (tao->ops->update) {
    228:       PetscUseTypeMethod(tao, update, tao->niter, tao->user_update);
    229:       PetscCall(TaoComputeObjective(tao, tao->solution, &f));
    230:     }
    231:     ++tao->niter;
    232:     tao->ksp_its = 0;
    233:     /* Compute the Hessian */
    234:     if (needH) PetscCall(TaoComputeHessian(tao, tao->solution, tao->hessian, tao->hessian_pre));

    236:     if (tl->bfgs_pre) {
    237:       /* Update the limited memory preconditioner */
    238:       PetscCall(MatLMVMUpdate(tl->M, tao->solution, tao->gradient));
    239:       ++bfgsUpdates;
    240:     }
    241:     PetscCall(KSPSetOperators(tao->ksp, tao->hessian, tao->hessian_pre));
    242:     /* Solve the Newton system of equations */
    243:     PetscCall(KSPCGSetRadius(tao->ksp, tl->max_radius));
    244:     PetscCall(KSPSolve(tao->ksp, tao->gradient, tao->stepdirection));
    245:     PetscCall(KSPGetIterationNumber(tao->ksp, &its));
    246:     tao->ksp_its += its;
    247:     tao->ksp_tot_its += its;
    248:     PetscCall(KSPCGGetNormD(tao->ksp, &norm_d));

    250:     if (0.0 == tao->trust) {
    251:       /* Radius was uninitialized; use the norm of the direction */
    252:       if (norm_d > 0.0) {
    253:         tao->trust = norm_d;

    255:         /* Modify the radius if it is too large or small */
    256:         tao->trust = PetscMax(tao->trust, tl->min_radius);
    257:         tao->trust = PetscMin(tao->trust, tl->max_radius);
    258:       } else {
    259:         /* The direction was bad; set radius to default value and re-solve
    260:            the trust-region subproblem to get a direction */
    261:         tao->trust = tao->trust0;

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

    267:         PetscCall(KSPCGSetRadius(tao->ksp, tl->max_radius));
    268:         PetscCall(KSPSolve(tao->ksp, tao->gradient, tao->stepdirection));
    269:         PetscCall(KSPGetIterationNumber(tao->ksp, &its));
    270:         tao->ksp_its += its;
    271:         tao->ksp_tot_its += its;
    272:         PetscCall(KSPCGGetNormD(tao->ksp, &norm_d));

    274:         PetscCheck(norm_d != 0.0, PetscObjectComm((PetscObject)tao), PETSC_ERR_PLIB, "Initial direction zero");
    275:       }
    276:     }

    278:     PetscCall(VecScale(tao->stepdirection, -1.0));
    279:     PetscCall(KSPGetConvergedReason(tao->ksp, &ksp_reason));
    280:     if ((KSP_DIVERGED_INDEFINITE_PC == ksp_reason) && (tl->bfgs_pre)) {
    281:       /* Preconditioner is numerically indefinite; reset the
    282:          approximate if using BFGS preconditioning. */
    283:       PetscCall(MatLMVMReset(tl->M, PETSC_FALSE));
    284:       PetscCall(MatLMVMUpdate(tl->M, tao->solution, tao->gradient));
    285:       bfgsUpdates = 1;
    286:     }

    288:     /* Check trust-region reduction conditions */
    289:     tr_reject = 0;
    290:     if (NTL_UPDATE_REDUCTION == tl->update_type) {
    291:       /* Get predicted reduction */
    292:       PetscCall(KSPCGGetObjFcn(tao->ksp, &prered));
    293:       if (prered >= 0.0) {
    294:         /* The predicted reduction has the wrong sign.  This cannot
    295:            happen in infinite precision arithmetic.  Step should
    296:            be rejected! */
    297:         tao->trust = tl->alpha1 * PetscMin(tao->trust, norm_d);
    298:         tr_reject  = 1;
    299:       } else {
    300:         /* Compute trial step and function value */
    301:         PetscCall(VecCopy(tao->solution, tl->W));
    302:         PetscCall(VecAXPY(tl->W, 1.0, tao->stepdirection));
    303:         PetscCall(TaoComputeObjective(tao, tl->W, &ftrial));

    305:         if (PetscIsInfOrNanReal(ftrial)) {
    306:           tao->trust = tl->alpha1 * PetscMin(tao->trust, norm_d);
    307:           tr_reject  = 1;
    308:         } else {
    309:           /* Compute and actual reduction */
    310:           actred = f - ftrial;
    311:           prered = -prered;
    312:           if ((PetscAbsScalar(actred) <= tl->epsilon) && (PetscAbsScalar(prered) <= tl->epsilon)) {
    313:             kappa = 1.0;
    314:           } else {
    315:             kappa = actred / prered;
    316:           }

    318:           /* Accept of reject the step and update radius */
    319:           if (kappa < tl->eta1) {
    320:             /* Reject the step */
    321:             tao->trust = tl->alpha1 * PetscMin(tao->trust, norm_d);
    322:             tr_reject  = 1;
    323:           } else {
    324:             /* Accept the step */
    325:             if (kappa < tl->eta2) {
    326:               /* Marginal bad step */
    327:               tao->trust = tl->alpha2 * PetscMin(tao->trust, norm_d);
    328:             } else if (kappa < tl->eta3) {
    329:               /* Reasonable step */
    330:               tao->trust = tl->alpha3 * tao->trust;
    331:             } else if (kappa < tl->eta4) {
    332:               /* Good step */
    333:               tao->trust = PetscMax(tl->alpha4 * norm_d, tao->trust);
    334:             } else {
    335:               /* Very good step */
    336:               tao->trust = PetscMax(tl->alpha5 * norm_d, tao->trust);
    337:             }
    338:           }
    339:         }
    340:       }
    341:     } else {
    342:       /* Get predicted reduction */
    343:       PetscCall(KSPCGGetObjFcn(tao->ksp, &prered));
    344:       if (prered >= 0.0) {
    345:         /* The predicted reduction has the wrong sign.  This cannot
    346:            happen in infinite precision arithmetic.  Step should
    347:            be rejected! */
    348:         tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d);
    349:         tr_reject  = 1;
    350:       } else {
    351:         PetscCall(VecCopy(tao->solution, tl->W));
    352:         PetscCall(VecAXPY(tl->W, 1.0, tao->stepdirection));
    353:         PetscCall(TaoComputeObjective(tao, tl->W, &ftrial));
    354:         if (PetscIsInfOrNanReal(ftrial)) {
    355:           tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d);
    356:           tr_reject  = 1;
    357:         } else {
    358:           PetscCall(VecDot(tao->gradient, tao->stepdirection, &gdx));

    360:           actred = f - ftrial;
    361:           prered = -prered;
    362:           if ((PetscAbsScalar(actred) <= tl->epsilon) && (PetscAbsScalar(prered) <= tl->epsilon)) {
    363:             kappa = 1.0;
    364:           } else {
    365:             kappa = actred / prered;
    366:           }

    368:           tau_1   = tl->theta * gdx / (tl->theta * gdx - (1.0 - tl->theta) * prered + actred);
    369:           tau_2   = tl->theta * gdx / (tl->theta * gdx + (1.0 + tl->theta) * prered - actred);
    370:           tau_min = PetscMin(tau_1, tau_2);
    371:           tau_max = PetscMax(tau_1, tau_2);

    373:           if (kappa >= 1.0 - tl->mu1) {
    374:             /* Great agreement; accept step and update radius */
    375:             if (tau_max < 1.0) {
    376:               tao->trust = PetscMax(tao->trust, tl->gamma3 * norm_d);
    377:             } else if (tau_max > tl->gamma4) {
    378:               tao->trust = PetscMax(tao->trust, tl->gamma4 * norm_d);
    379:             } else {
    380:               tao->trust = PetscMax(tao->trust, tau_max * norm_d);
    381:             }
    382:           } else if (kappa >= 1.0 - tl->mu2) {
    383:             /* Good agreement */

    385:             if (tau_max < tl->gamma2) {
    386:               tao->trust = tl->gamma2 * PetscMin(tao->trust, norm_d);
    387:             } else if (tau_max > tl->gamma3) {
    388:               tao->trust = PetscMax(tao->trust, tl->gamma3 * norm_d);
    389:             } else if (tau_max < 1.0) {
    390:               tao->trust = tau_max * PetscMin(tao->trust, norm_d);
    391:             } else {
    392:               tao->trust = PetscMax(tao->trust, tau_max * norm_d);
    393:             }
    394:           } else {
    395:             /* Not good agreement */
    396:             if (tau_min > 1.0) {
    397:               tao->trust = tl->gamma2 * PetscMin(tao->trust, norm_d);
    398:             } else if (tau_max < tl->gamma1) {
    399:               tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d);
    400:             } else if ((tau_min < tl->gamma1) && (tau_max >= 1.0)) {
    401:               tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d);
    402:             } else if ((tau_1 >= tl->gamma1) && (tau_1 < 1.0) && ((tau_2 < tl->gamma1) || (tau_2 >= 1.0))) {
    403:               tao->trust = tau_1 * PetscMin(tao->trust, norm_d);
    404:             } else if ((tau_2 >= tl->gamma1) && (tau_2 < 1.0) && ((tau_1 < tl->gamma1) || (tau_2 >= 1.0))) {
    405:               tao->trust = tau_2 * PetscMin(tao->trust, norm_d);
    406:             } else {
    407:               tao->trust = tau_max * PetscMin(tao->trust, norm_d);
    408:             }
    409:             tr_reject = 1;
    410:           }
    411:         }
    412:       }
    413:     }

    415:     if (tr_reject) {
    416:       /* The trust-region constraints rejected the step.  Apply a linesearch.
    417:          Check for descent direction. */
    418:       PetscCall(VecDot(tao->stepdirection, tao->gradient, &gdx));
    419:       if ((gdx >= 0.0) || PetscIsInfOrNanReal(gdx)) {
    420:         /* Newton step is not descent or direction produced Inf or NaN */

    422:         if (!tl->bfgs_pre) {
    423:           /* We don't have the bfgs matrix around and updated
    424:              Must use gradient direction in this case */
    425:           PetscCall(VecCopy(tao->gradient, tao->stepdirection));
    426:           PetscCall(VecScale(tao->stepdirection, -1.0));
    427:           ++tl->grad;
    428:           stepType = NTL_GRADIENT;
    429:         } else {
    430:           /* Attempt to use the BFGS direction */
    431:           PetscCall(MatSolve(tl->M, tao->gradient, tao->stepdirection));
    432:           PetscCall(VecScale(tao->stepdirection, -1.0));

    434:           /* Check for success (descent direction) */
    435:           PetscCall(VecDot(tao->stepdirection, tao->gradient, &gdx));
    436:           if ((gdx >= 0) || PetscIsInfOrNanReal(gdx)) {
    437:             /* BFGS direction is not descent or direction produced not a number
    438:                We can assert bfgsUpdates > 1 in this case because
    439:                the first solve produces the scaled gradient direction,
    440:                which is guaranteed to be descent */

    442:             /* Use steepest descent direction (scaled) */
    443:             PetscCall(MatLMVMReset(tl->M, PETSC_FALSE));
    444:             PetscCall(MatLMVMUpdate(tl->M, tao->solution, tao->gradient));
    445:             PetscCall(MatSolve(tl->M, tao->gradient, tao->stepdirection));
    446:             PetscCall(VecScale(tao->stepdirection, -1.0));

    448:             bfgsUpdates = 1;
    449:             ++tl->grad;
    450:             stepType = NTL_GRADIENT;
    451:           } else {
    452:             PetscCall(MatLMVMGetUpdateCount(tl->M, &bfgsUpdates));
    453:             if (1 == bfgsUpdates) {
    454:               /* The first BFGS direction is always the scaled gradient */
    455:               ++tl->grad;
    456:               stepType = NTL_GRADIENT;
    457:             } else {
    458:               ++tl->bfgs;
    459:               stepType = NTL_BFGS;
    460:             }
    461:           }
    462:         }
    463:       } else {
    464:         /* Computed Newton step is descent */
    465:         ++tl->newt;
    466:         stepType = NTL_NEWTON;
    467:       }

    469:       /* Perform the linesearch */
    470:       fold = f;
    471:       PetscCall(VecCopy(tao->solution, tl->Xold));
    472:       PetscCall(VecCopy(tao->gradient, tl->Gold));

    474:       step = 1.0;
    475:       PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_reason));
    476:       PetscCall(TaoAddLineSearchCounts(tao));

    478:       while (ls_reason != TAOLINESEARCH_SUCCESS && ls_reason != TAOLINESEARCH_SUCCESS_USER && stepType != NTL_GRADIENT) { /* Linesearch failed */
    479:         /* Linesearch failed */
    480:         f = fold;
    481:         PetscCall(VecCopy(tl->Xold, tao->solution));
    482:         PetscCall(VecCopy(tl->Gold, tao->gradient));

    484:         switch (stepType) {
    485:         case NTL_NEWTON:
    486:           /* Failed to obtain acceptable iterate with Newton step */

    488:           if (tl->bfgs_pre) {
    489:             /* We don't have the bfgs matrix around and being updated
    490:                Must use gradient direction in this case */
    491:             PetscCall(VecCopy(tao->gradient, tao->stepdirection));
    492:             ++tl->grad;
    493:             stepType = NTL_GRADIENT;
    494:           } else {
    495:             /* Attempt to use the BFGS direction */
    496:             PetscCall(MatSolve(tl->M, tao->gradient, tao->stepdirection));

    498:             /* Check for success (descent direction) */
    499:             PetscCall(VecDot(tao->stepdirection, tao->gradient, &gdx));
    500:             if ((gdx <= 0) || PetscIsInfOrNanReal(gdx)) {
    501:               /* BFGS direction is not descent or direction produced
    502:                  not a number.  We can assert bfgsUpdates > 1 in this case
    503:                  Use steepest descent direction (scaled) */
    504:               PetscCall(MatLMVMReset(tl->M, PETSC_FALSE));
    505:               PetscCall(MatLMVMUpdate(tl->M, tao->solution, tao->gradient));
    506:               PetscCall(MatSolve(tl->M, tao->gradient, tao->stepdirection));

    508:               bfgsUpdates = 1;
    509:               ++tl->grad;
    510:               stepType = NTL_GRADIENT;
    511:             } else {
    512:               PetscCall(MatLMVMGetUpdateCount(tl->M, &bfgsUpdates));
    513:               if (1 == bfgsUpdates) {
    514:                 /* The first BFGS direction is always the scaled gradient */
    515:                 ++tl->grad;
    516:                 stepType = NTL_GRADIENT;
    517:               } else {
    518:                 ++tl->bfgs;
    519:                 stepType = NTL_BFGS;
    520:               }
    521:             }
    522:           }
    523:           break;

    525:         case NTL_BFGS:
    526:           /* Can only enter if pc_type == NTL_PC_BFGS
    527:              Failed to obtain acceptable iterate with BFGS step
    528:              Attempt to use the scaled gradient direction */
    529:           PetscCall(MatLMVMReset(tl->M, PETSC_FALSE));
    530:           PetscCall(MatLMVMUpdate(tl->M, tao->solution, tao->gradient));
    531:           PetscCall(MatSolve(tl->M, tao->gradient, tao->stepdirection));

    533:           bfgsUpdates = 1;
    534:           ++tl->grad;
    535:           stepType = NTL_GRADIENT;
    536:           break;
    537:         }
    538:         PetscCall(VecScale(tao->stepdirection, -1.0));

    540:         /* This may be incorrect; linesearch has values for stepmax and stepmin
    541:            that should be reset. */
    542:         step = 1.0;
    543:         PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_reason));
    544:         PetscCall(TaoAddLineSearchCounts(tao));
    545:       }

    547:       if (ls_reason != TAOLINESEARCH_SUCCESS && ls_reason != TAOLINESEARCH_SUCCESS_USER) {
    548:         /* Failed to find an improving point */
    549:         f = fold;
    550:         PetscCall(VecCopy(tl->Xold, tao->solution));
    551:         PetscCall(VecCopy(tl->Gold, tao->gradient));
    552:         tao->trust  = 0.0;
    553:         step        = 0.0;
    554:         tao->reason = TAO_DIVERGED_LS_FAILURE;
    555:         break;
    556:       } else if (stepType == NTL_NEWTON) {
    557:         if (step < tl->nu1) {
    558:           /* Very bad step taken; reduce radius */
    559:           tao->trust = tl->omega1 * PetscMin(norm_d, tao->trust);
    560:         } else if (step < tl->nu2) {
    561:           /* Reasonably bad step taken; reduce radius */
    562:           tao->trust = tl->omega2 * PetscMin(norm_d, tao->trust);
    563:         } else if (step < tl->nu3) {
    564:           /* Reasonable step was taken; leave radius alone */
    565:           if (tl->omega3 < 1.0) {
    566:             tao->trust = tl->omega3 * PetscMin(norm_d, tao->trust);
    567:           } else if (tl->omega3 > 1.0) {
    568:             tao->trust = PetscMax(tl->omega3 * norm_d, tao->trust);
    569:           }
    570:         } else if (step < tl->nu4) {
    571:           /* Full step taken; increase the radius */
    572:           tao->trust = PetscMax(tl->omega4 * norm_d, tao->trust);
    573:         } else {
    574:           /* More than full step taken; increase the radius */
    575:           tao->trust = PetscMax(tl->omega5 * norm_d, tao->trust);
    576:         }
    577:       } else {
    578:         /* Newton step was not good; reduce the radius */
    579:         tao->trust = tl->omega1 * PetscMin(norm_d, tao->trust);
    580:       }
    581:     } else {
    582:       /* Trust-region step is accepted */
    583:       PetscCall(VecCopy(tl->W, tao->solution));
    584:       f = ftrial;
    585:       PetscCall(TaoComputeGradient(tao, tao->solution, tao->gradient));
    586:       ++tl->ntrust;
    587:     }

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

    592:     /* Check for converged */
    593:     PetscCall(VecNorm(tao->gradient, NORM_2, &gnorm));
    594:     PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated Not-a-Number");
    595:     needH = 1;

    597:     PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its));
    598:     PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step));
    599:     PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
    600:   }
    601:   PetscFunctionReturn(PETSC_SUCCESS);
    602: }

    604: /* ---------------------------------------------------------- */
    605: static PetscErrorCode TaoSetUp_NTL(Tao tao)
    606: {
    607:   TAO_NTL *tl = (TAO_NTL *)tao->data;

    609:   PetscFunctionBegin;
    610:   if (!tao->gradient) PetscCall(VecDuplicate(tao->solution, &tao->gradient));
    611:   if (!tao->stepdirection) PetscCall(VecDuplicate(tao->solution, &tao->stepdirection));
    612:   if (!tl->W) PetscCall(VecDuplicate(tao->solution, &tl->W));
    613:   if (!tl->Xold) PetscCall(VecDuplicate(tao->solution, &tl->Xold));
    614:   if (!tl->Gold) PetscCall(VecDuplicate(tao->solution, &tl->Gold));
    615:   tl->bfgs_pre = NULL;
    616:   tl->M        = NULL;
    617:   PetscFunctionReturn(PETSC_SUCCESS);
    618: }

    620: /*------------------------------------------------------------*/
    621: static PetscErrorCode TaoDestroy_NTL(Tao tao)
    622: {
    623:   TAO_NTL *tl = (TAO_NTL *)tao->data;

    625:   PetscFunctionBegin;
    626:   if (tao->setupcalled) {
    627:     PetscCall(VecDestroy(&tl->W));
    628:     PetscCall(VecDestroy(&tl->Xold));
    629:     PetscCall(VecDestroy(&tl->Gold));
    630:   }
    631:   PetscCall(KSPDestroy(&tao->ksp));
    632:   PetscCall(PetscFree(tao->data));
    633:   PetscFunctionReturn(PETSC_SUCCESS);
    634: }

    636: /*------------------------------------------------------------*/
    637: static PetscErrorCode TaoSetFromOptions_NTL(Tao tao, PetscOptionItems *PetscOptionsObject)
    638: {
    639:   TAO_NTL *tl = (TAO_NTL *)tao->data;

    641:   PetscFunctionBegin;
    642:   PetscOptionsHeadBegin(PetscOptionsObject, "Newton trust region with line search method for unconstrained optimization");
    643:   PetscCall(PetscOptionsEList("-tao_ntl_init_type", "radius initialization type", "", NTL_INIT, NTL_INIT_TYPES, NTL_INIT[tl->init_type], &tl->init_type, NULL));
    644:   PetscCall(PetscOptionsEList("-tao_ntl_update_type", "radius update type", "", NTL_UPDATE, NTL_UPDATE_TYPES, NTL_UPDATE[tl->update_type], &tl->update_type, NULL));
    645:   PetscCall(PetscOptionsReal("-tao_ntl_eta1", "poor steplength; reduce radius", "", tl->eta1, &tl->eta1, NULL));
    646:   PetscCall(PetscOptionsReal("-tao_ntl_eta2", "reasonable steplength; leave radius alone", "", tl->eta2, &tl->eta2, NULL));
    647:   PetscCall(PetscOptionsReal("-tao_ntl_eta3", "good steplength; increase radius", "", tl->eta3, &tl->eta3, NULL));
    648:   PetscCall(PetscOptionsReal("-tao_ntl_eta4", "excellent steplength; greatly increase radius", "", tl->eta4, &tl->eta4, NULL));
    649:   PetscCall(PetscOptionsReal("-tao_ntl_alpha1", "", "", tl->alpha1, &tl->alpha1, NULL));
    650:   PetscCall(PetscOptionsReal("-tao_ntl_alpha2", "", "", tl->alpha2, &tl->alpha2, NULL));
    651:   PetscCall(PetscOptionsReal("-tao_ntl_alpha3", "", "", tl->alpha3, &tl->alpha3, NULL));
    652:   PetscCall(PetscOptionsReal("-tao_ntl_alpha4", "", "", tl->alpha4, &tl->alpha4, NULL));
    653:   PetscCall(PetscOptionsReal("-tao_ntl_alpha5", "", "", tl->alpha5, &tl->alpha5, NULL));
    654:   PetscCall(PetscOptionsReal("-tao_ntl_nu1", "poor steplength; reduce radius", "", tl->nu1, &tl->nu1, NULL));
    655:   PetscCall(PetscOptionsReal("-tao_ntl_nu2", "reasonable steplength; leave radius alone", "", tl->nu2, &tl->nu2, NULL));
    656:   PetscCall(PetscOptionsReal("-tao_ntl_nu3", "good steplength; increase radius", "", tl->nu3, &tl->nu3, NULL));
    657:   PetscCall(PetscOptionsReal("-tao_ntl_nu4", "excellent steplength; greatly increase radius", "", tl->nu4, &tl->nu4, NULL));
    658:   PetscCall(PetscOptionsReal("-tao_ntl_omega1", "", "", tl->omega1, &tl->omega1, NULL));
    659:   PetscCall(PetscOptionsReal("-tao_ntl_omega2", "", "", tl->omega2, &tl->omega2, NULL));
    660:   PetscCall(PetscOptionsReal("-tao_ntl_omega3", "", "", tl->omega3, &tl->omega3, NULL));
    661:   PetscCall(PetscOptionsReal("-tao_ntl_omega4", "", "", tl->omega4, &tl->omega4, NULL));
    662:   PetscCall(PetscOptionsReal("-tao_ntl_omega5", "", "", tl->omega5, &tl->omega5, NULL));
    663:   PetscCall(PetscOptionsReal("-tao_ntl_mu1_i", "", "", tl->mu1_i, &tl->mu1_i, NULL));
    664:   PetscCall(PetscOptionsReal("-tao_ntl_mu2_i", "", "", tl->mu2_i, &tl->mu2_i, NULL));
    665:   PetscCall(PetscOptionsReal("-tao_ntl_gamma1_i", "", "", tl->gamma1_i, &tl->gamma1_i, NULL));
    666:   PetscCall(PetscOptionsReal("-tao_ntl_gamma2_i", "", "", tl->gamma2_i, &tl->gamma2_i, NULL));
    667:   PetscCall(PetscOptionsReal("-tao_ntl_gamma3_i", "", "", tl->gamma3_i, &tl->gamma3_i, NULL));
    668:   PetscCall(PetscOptionsReal("-tao_ntl_gamma4_i", "", "", tl->gamma4_i, &tl->gamma4_i, NULL));
    669:   PetscCall(PetscOptionsReal("-tao_ntl_theta_i", "", "", tl->theta_i, &tl->theta_i, NULL));
    670:   PetscCall(PetscOptionsReal("-tao_ntl_mu1", "", "", tl->mu1, &tl->mu1, NULL));
    671:   PetscCall(PetscOptionsReal("-tao_ntl_mu2", "", "", tl->mu2, &tl->mu2, NULL));
    672:   PetscCall(PetscOptionsReal("-tao_ntl_gamma1", "", "", tl->gamma1, &tl->gamma1, NULL));
    673:   PetscCall(PetscOptionsReal("-tao_ntl_gamma2", "", "", tl->gamma2, &tl->gamma2, NULL));
    674:   PetscCall(PetscOptionsReal("-tao_ntl_gamma3", "", "", tl->gamma3, &tl->gamma3, NULL));
    675:   PetscCall(PetscOptionsReal("-tao_ntl_gamma4", "", "", tl->gamma4, &tl->gamma4, NULL));
    676:   PetscCall(PetscOptionsReal("-tao_ntl_theta", "", "", tl->theta, &tl->theta, NULL));
    677:   PetscCall(PetscOptionsReal("-tao_ntl_min_radius", "lower bound on initial radius", "", tl->min_radius, &tl->min_radius, NULL));
    678:   PetscCall(PetscOptionsReal("-tao_ntl_max_radius", "upper bound on radius", "", tl->max_radius, &tl->max_radius, NULL));
    679:   PetscCall(PetscOptionsReal("-tao_ntl_epsilon", "tolerance used when computing actual and predicted reduction", "", tl->epsilon, &tl->epsilon, NULL));
    680:   PetscOptionsHeadEnd();
    681:   PetscCall(TaoLineSearchSetFromOptions(tao->linesearch));
    682:   PetscCall(KSPSetFromOptions(tao->ksp));
    683:   PetscFunctionReturn(PETSC_SUCCESS);
    684: }

    686: /*------------------------------------------------------------*/
    687: static PetscErrorCode TaoView_NTL(Tao tao, PetscViewer viewer)
    688: {
    689:   TAO_NTL  *tl = (TAO_NTL *)tao->data;
    690:   PetscBool isascii;

    692:   PetscFunctionBegin;
    693:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    694:   if (isascii) {
    695:     PetscCall(PetscViewerASCIIPushTab(viewer));
    696:     PetscCall(PetscViewerASCIIPrintf(viewer, "Trust-region steps: %" PetscInt_FMT "\n", tl->ntrust));
    697:     PetscCall(PetscViewerASCIIPrintf(viewer, "Newton search steps: %" PetscInt_FMT "\n", tl->newt));
    698:     PetscCall(PetscViewerASCIIPrintf(viewer, "BFGS search steps: %" PetscInt_FMT "\n", tl->bfgs));
    699:     PetscCall(PetscViewerASCIIPrintf(viewer, "Gradient search steps: %" PetscInt_FMT "\n", tl->grad));
    700:     PetscCall(PetscViewerASCIIPopTab(viewer));
    701:   }
    702:   PetscFunctionReturn(PETSC_SUCCESS);
    703: }

    705: /* ---------------------------------------------------------- */
    706: /*MC
    707:   TAONTL - Newton's method with trust region globalization and line search fallback.
    708:   At each iteration, the Newton trust region method solves the system for d
    709:   and performs a line search in the d direction:

    711:             min_d  .5 dT Hk d + gkT d,  s.t.   ||d|| < Delta_k

    713:   Options Database Keys:
    714: + -tao_ntl_init_type - "constant","direction","interpolation"
    715: . -tao_ntl_update_type - "reduction","interpolation"
    716: . -tao_ntl_min_radius - lower bound on trust region radius
    717: . -tao_ntl_max_radius - upper bound on trust region radius
    718: . -tao_ntl_epsilon - tolerance for accepting actual / predicted reduction
    719: . -tao_ntl_mu1_i - mu1 interpolation init factor
    720: . -tao_ntl_mu2_i - mu2 interpolation init factor
    721: . -tao_ntl_gamma1_i - gamma1 interpolation init factor
    722: . -tao_ntl_gamma2_i - gamma2 interpolation init factor
    723: . -tao_ntl_gamma3_i - gamma3 interpolation init factor
    724: . -tao_ntl_gamma4_i - gamma4 interpolation init factor
    725: . -tao_ntl_theta_i - theta1 interpolation init factor
    726: . -tao_ntl_eta1 - eta1 reduction update factor
    727: . -tao_ntl_eta2 - eta2 reduction update factor
    728: . -tao_ntl_eta3 - eta3 reduction update factor
    729: . -tao_ntl_eta4 - eta4 reduction update factor
    730: . -tao_ntl_alpha1 - alpha1 reduction update factor
    731: . -tao_ntl_alpha2 - alpha2 reduction update factor
    732: . -tao_ntl_alpha3 - alpha3 reduction update factor
    733: . -tao_ntl_alpha4 - alpha4 reduction update factor
    734: . -tao_ntl_alpha4 - alpha4 reduction update factor
    735: . -tao_ntl_mu1 - mu1 interpolation update
    736: . -tao_ntl_mu2 - mu2 interpolation update
    737: . -tao_ntl_gamma1 - gamma1 interpolcation update
    738: . -tao_ntl_gamma2 - gamma2 interpolcation update
    739: . -tao_ntl_gamma3 - gamma3 interpolcation update
    740: . -tao_ntl_gamma4 - gamma4 interpolation update
    741: - -tao_ntl_theta - theta1 interpolation update

    743:   Level: beginner
    744: M*/
    745: PETSC_EXTERN PetscErrorCode TaoCreate_NTL(Tao tao)
    746: {
    747:   TAO_NTL    *tl;
    748:   const char *morethuente_type = TAOLINESEARCHMT;

    750:   PetscFunctionBegin;
    751:   PetscCall(PetscNew(&tl));
    752:   tao->ops->setup          = TaoSetUp_NTL;
    753:   tao->ops->solve          = TaoSolve_NTL;
    754:   tao->ops->view           = TaoView_NTL;
    755:   tao->ops->setfromoptions = TaoSetFromOptions_NTL;
    756:   tao->ops->destroy        = TaoDestroy_NTL;

    758:   /* Override default settings (unless already changed) */
    759:   PetscCall(TaoParametersInitialize(tao));
    760:   PetscObjectParameterSetDefault(tao, max_it, 50);
    761:   PetscObjectParameterSetDefault(tao, trust0, 100.0);

    763:   tao->data = (void *)tl;

    765:   /* Default values for trust-region radius update based on steplength */
    766:   tl->nu1 = 0.25;
    767:   tl->nu2 = 0.50;
    768:   tl->nu3 = 1.00;
    769:   tl->nu4 = 1.25;

    771:   tl->omega1 = 0.25;
    772:   tl->omega2 = 0.50;
    773:   tl->omega3 = 1.00;
    774:   tl->omega4 = 2.00;
    775:   tl->omega5 = 4.00;

    777:   /* Default values for trust-region radius update based on reduction */
    778:   tl->eta1 = 1.0e-4;
    779:   tl->eta2 = 0.25;
    780:   tl->eta3 = 0.50;
    781:   tl->eta4 = 0.90;

    783:   tl->alpha1 = 0.25;
    784:   tl->alpha2 = 0.50;
    785:   tl->alpha3 = 1.00;
    786:   tl->alpha4 = 2.00;
    787:   tl->alpha5 = 4.00;

    789:   /* Default values for trust-region radius update based on interpolation */
    790:   tl->mu1 = 0.10;
    791:   tl->mu2 = 0.50;

    793:   tl->gamma1 = 0.25;
    794:   tl->gamma2 = 0.50;
    795:   tl->gamma3 = 2.00;
    796:   tl->gamma4 = 4.00;

    798:   tl->theta = 0.05;

    800:   /* Default values for trust region initialization based on interpolation */
    801:   tl->mu1_i = 0.35;
    802:   tl->mu2_i = 0.50;

    804:   tl->gamma1_i = 0.0625;
    805:   tl->gamma2_i = 0.5;
    806:   tl->gamma3_i = 2.0;
    807:   tl->gamma4_i = 5.0;

    809:   tl->theta_i = 0.25;

    811:   /* Remaining parameters */
    812:   tl->min_radius = 1.0e-10;
    813:   tl->max_radius = 1.0e10;
    814:   tl->epsilon    = 1.0e-6;

    816:   tl->init_type   = NTL_INIT_INTERPOLATION;
    817:   tl->update_type = NTL_UPDATE_REDUCTION;

    819:   PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch));
    820:   PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1));
    821:   PetscCall(TaoLineSearchSetType(tao->linesearch, morethuente_type));
    822:   PetscCall(TaoLineSearchUseTaoRoutines(tao->linesearch, tao));
    823:   PetscCall(TaoLineSearchSetOptionsPrefix(tao->linesearch, tao->hdr.prefix));
    824:   PetscCall(KSPCreate(((PetscObject)tao)->comm, &tao->ksp));
    825:   PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1));
    826:   PetscCall(KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix));
    827:   PetscCall(KSPAppendOptionsPrefix(tao->ksp, "tao_ntl_"));
    828:   PetscCall(KSPSetType(tao->ksp, KSPSTCG));
    829:   PetscFunctionReturn(PETSC_SUCCESS);
    830: }