Actual source code: gltr.c

      1: #include <../src/ksp/ksp/impls/cg/gltr/gltrimpl.h>
      2: #include <petscblaslapack.h>

      4: #define GLTR_PRECONDITIONED_DIRECTION   0
      5: #define GLTR_UNPRECONDITIONED_DIRECTION 1
      6: #define GLTR_DIRECTION_TYPES            2

      8: static const char *DType_Table[64] = {"preconditioned", "unpreconditioned"};

     10: /*@
     11:   KSPGLTRGetMinEig - Get minimum eigenvalue computed by `KSPGLTR`

     13:   Collective

     15:   Input Parameter:
     16: . ksp - the iterative context

     18:   Output Parameter:
     19: . e_min - the minimum eigenvalue

     21:   Level: advanced

     23: .seealso: [](ch_ksp), `KSP`, `KSPGLTR`, `KSPGLTRGetLambda()`
     24: @*/
     25: PetscErrorCode KSPGLTRGetMinEig(KSP ksp, PetscReal *e_min)
     26: {
     27:   PetscFunctionBegin;
     29:   PetscUseMethod(ksp, "KSPGLTRGetMinEig_C", (KSP, PetscReal *), (ksp, e_min));
     30:   PetscFunctionReturn(PETSC_SUCCESS);
     31: }

     33: /*@
     34:   KSPGLTRGetLambda - Get the multiplier on the trust-region constraint when using `KSPGLTR`

     36:   Not Collective

     38:   Input Parameter:
     39: . ksp - the iterative context

     41:   Output Parameter:
     42: . lambda - the multiplier

     44:   Level: advanced

     46: .seealso: [](ch_ksp), `KSP`, `KSPGLTR`, `KSPGLTRGetMinEig()`
     47: @*/
     48: PetscErrorCode KSPGLTRGetLambda(KSP ksp, PetscReal *lambda)
     49: {
     50:   PetscFunctionBegin;
     52:   PetscUseMethod(ksp, "KSPGLTRGetLambda_C", (KSP, PetscReal *), (ksp, lambda));
     53:   PetscFunctionReturn(PETSC_SUCCESS);
     54: }

     56: static PetscErrorCode KSPCGSolve_GLTR(KSP ksp)
     57: {
     58: #if defined(PETSC_USE_COMPLEX)
     59:   SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "GLTR is not available for complex systems");
     60: #else
     61:   KSPCG_GLTR   *cg = (KSPCG_GLTR *)ksp->data;
     62:   PetscReal    *t_soln, *t_diag, *t_offd, *e_valu, *e_vect, *e_rwrk;
     63:   PetscBLASInt *e_iblk, *e_splt, *e_iwrk;

     65:   Mat Qmat, Mmat;
     66:   Vec r, z, p, d;
     67:   PC  pc;

     69:   PetscReal norm_r, norm_d, norm_dp1, norm_p, dMp;
     70:   PetscReal alpha, beta, kappa, rz, rzm1;
     71:   PetscReal rr, r2, piv, step;
     72:   PetscReal vl, vu;
     73:   PetscReal coef1, coef2, coef3, root1, root2, obj1, obj2;
     74:   PetscReal norm_t, norm_w, pert;

     76:   PetscInt     i, j, max_cg_its, max_lanczos_its, max_newton_its, sigma;
     77:   PetscBLASInt t_size = 0, l_size = 0, il, iu, info;
     78:   PetscBLASInt nrhs, nldb;

     80:   PetscBLASInt e_valus = 0, e_splts;
     81:   PetscBool    diagonalscale;

     83:   PetscFunctionBegin;
     84:   /* Check the arguments and parameters.                                     */
     85:   PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
     86:   PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
     87:   PetscCheck(cg->radius >= 0.0, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "Input error: radius < 0");

     89:   /* Get the workspace vectors and initialize variables                      */
     90:   r2 = cg->radius * cg->radius;
     91:   r  = ksp->work[0];
     92:   z  = ksp->work[1];
     93:   p  = ksp->work[2];
     94:   d  = ksp->vec_sol;
     95:   pc = ksp->pc;

     97:   PetscCall(PCGetOperators(pc, &Qmat, &Mmat));

     99:   PetscCall(VecGetSize(d, &max_cg_its));
    100:   max_cg_its      = PetscMin(max_cg_its, ksp->max_it);
    101:   max_lanczos_its = cg->max_lanczos_its;
    102:   max_newton_its  = cg->max_newton_its;
    103:   ksp->its        = 0;

    105:   /* Initialize objective function direction, and minimum eigenvalue.        */
    106:   cg->o_fcn = 0.0;

    108:   PetscCall(VecSet(d, 0.0)); /* d = 0             */
    109:   cg->norm_d = 0.0;

    111:   cg->e_min  = 0.0;
    112:   cg->lambda = 0.0;

    114:   /*
    115:     The first phase of GLTR performs a standard conjugate gradient method,
    116:     but stores the values required for the Lanczos matrix.  We switch to
    117:     the Lanczos when the conjugate gradient method breaks down.  Check the
    118:     right-hand side for numerical problems.  The check for not-a-number and
    119:     infinite values need be performed only once.
    120:   */
    121:   PetscCall(VecCopy(ksp->vec_rhs, r)); /* r = -grad         */
    122:   PetscCall(VecDot(r, r, &rr));        /* rr = r^T r        */
    123:   KSPCheckDot(ksp, rr);

    125:   /*
    126:     Check the preconditioner for numerical problems and for positive
    127:     definiteness.  The check for not-a-number and infinite values need be
    128:     performed only once.
    129:   */
    130:   PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r      */
    131:   PetscCall(VecDot(r, z, &rz));      /* rz = r^T inv(M) r */
    132:   if (PetscIsInfOrNanScalar(rz)) {
    133:     /*
    134:       The preconditioner contains not-a-number or an infinite value.
    135:       Return the gradient direction intersected with the trust region.
    136:     */
    137:     ksp->reason = KSP_DIVERGED_NANORINF;
    138:     PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: bad preconditioner: rz=%g\n", (double)rz));

    140:     if (cg->radius) {
    141:       if (r2 >= rr) {
    142:         alpha      = 1.0;
    143:         cg->norm_d = PetscSqrtReal(rr);
    144:       } else {
    145:         alpha      = PetscSqrtReal(r2 / rr);
    146:         cg->norm_d = cg->radius;
    147:       }

    149:       PetscCall(VecAXPY(d, alpha, r)); /* d = d + alpha r   */

    151:       /* Compute objective function.                                         */
    152:       PetscCall(KSP_MatMult(ksp, Qmat, d, z));
    153:       PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
    154:       PetscCall(VecDot(d, z, &cg->o_fcn));
    155:       cg->o_fcn = -cg->o_fcn;
    156:       ++ksp->its;
    157:     }
    158:     PetscFunctionReturn(PETSC_SUCCESS);
    159:   }

    161:   if (rz < 0.0) {
    162:     /*
    163:       The preconditioner is indefinite.  Because this is the first
    164:       and we do not have a direction yet, we use the gradient step.  Note
    165:       that we cannot use the preconditioned norm when computing the step
    166:       because the matrix is indefinite.
    167:     */
    168:     ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
    169:     PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: indefinite preconditioner: rz=%g\n", (double)rz));

    171:     if (cg->radius) {
    172:       if (r2 >= rr) {
    173:         alpha      = 1.0;
    174:         cg->norm_d = PetscSqrtReal(rr);
    175:       } else {
    176:         alpha      = PetscSqrtReal(r2 / rr);
    177:         cg->norm_d = cg->radius;
    178:       }

    180:       PetscCall(VecAXPY(d, alpha, r)); /* d = d + alpha r   */

    182:       /* Compute objective function.                                         */
    183:       PetscCall(KSP_MatMult(ksp, Qmat, d, z));
    184:       PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
    185:       PetscCall(VecDot(d, z, &cg->o_fcn));
    186:       cg->o_fcn = -cg->o_fcn;
    187:       ++ksp->its;
    188:     }
    189:     PetscFunctionReturn(PETSC_SUCCESS);
    190:   }

    192:   /*
    193:     As far as we know, the preconditioner is positive semidefinite.
    194:    Compute and log the residual.  Check convergence because this
    195:    initializes things, but do not terminate until at least one conjugate
    196:    gradient iteration has been performed.
    197:   */
    198:   cg->norm_r[0] = PetscSqrtReal(rz); /* norm_r = |r|_M    */

    200:   switch (ksp->normtype) {
    201:   case KSP_NORM_PRECONDITIONED:
    202:     PetscCall(VecNorm(z, NORM_2, &norm_r)); /* norm_r = |z|      */
    203:     break;

    205:   case KSP_NORM_UNPRECONDITIONED:
    206:     norm_r = PetscSqrtReal(rr); /* norm_r = |r|      */
    207:     break;

    209:   case KSP_NORM_NATURAL:
    210:     norm_r = cg->norm_r[0]; /* norm_r = |r|_M    */
    211:     break;

    213:   default:
    214:     norm_r = 0.0;
    215:     break;
    216:   }

    218:   PetscCall(KSPLogResidualHistory(ksp, norm_r));
    219:   PetscCall(KSPMonitor(ksp, ksp->its, norm_r));
    220:   ksp->rnorm = norm_r;

    222:   PetscCall((*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP));

    224:   /* Compute the first direction and update the iteration.                   */
    225:   PetscCall(VecCopy(z, p));                /* p = z             */
    226:   PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p         */
    227:   ++ksp->its;

    229:   /* Check the matrix for numerical problems.                                */
    230:   PetscCall(VecDot(p, z, &kappa)); /* kappa = p^T Q p   */
    231:   if (PetscIsInfOrNanScalar(kappa)) {
    232:     /*
    233:       The matrix produced not-a-number or an infinite value.  In this case
    234:       we must stop and use the gradient direction.  This condition need
    235:       only be checked once.
    236:       */
    237:     ksp->reason = KSP_DIVERGED_NANORINF;
    238:     PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: bad matrix: kappa=%g\n", (double)kappa));

    240:     if (cg->radius) {
    241:       if (r2 >= rr) {
    242:         alpha      = 1.0;
    243:         cg->norm_d = PetscSqrtReal(rr);
    244:       } else {
    245:         alpha      = PetscSqrtReal(r2 / rr);
    246:         cg->norm_d = cg->radius;
    247:       }

    249:       PetscCall(VecAXPY(d, alpha, r)); /* d = d + alpha r   */

    251:       /* Compute objective function.                                         */
    252:       PetscCall(KSP_MatMult(ksp, Qmat, d, z));
    253:       PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
    254:       PetscCall(VecDot(d, z, &cg->o_fcn));
    255:       cg->o_fcn = -cg->o_fcn;
    256:       ++ksp->its;
    257:     }
    258:     PetscFunctionReturn(PETSC_SUCCESS);
    259:   }

    261:   /*
    262:     Initialize variables for calculating the norm of the direction and for
    263:     the Lanczos tridiagonal matrix.  Note that we track the diagonal value
    264:     of the Cholesky factorization of the Lanczos matrix in order to
    265:     determine when negative curvature is encountered.
    266:   */
    267:   dMp    = 0.0;
    268:   norm_d = 0.0;
    269:   switch (cg->dtype) {
    270:   case GLTR_PRECONDITIONED_DIRECTION:
    271:     norm_p = rz;
    272:     break;

    274:   default:
    275:     PetscCall(VecDot(p, p, &norm_p));
    276:     break;
    277:   }

    279:   cg->diag[t_size] = kappa / rz;
    280:   cg->offd[t_size] = 0.0;
    281:   ++t_size;

    283:   piv = 1.0;

    285:   /*
    286:      Check for breakdown of the conjugate gradient method; this occurs when
    287:      kappa is zero.
    288:    */
    289:   if (PetscAbsReal(kappa) <= 0.0) {
    290:     /* The curvature is zero.  In this case, we must stop and use follow
    291:        the direction of negative curvature since the Lanczos matrix is zero. */
    292:     ksp->reason = KSP_DIVERGED_BREAKDOWN;
    293:     PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: breakdown: kappa=%g\n", (double)kappa));

    295:     if (cg->radius && norm_p > 0.0) {
    296:       /* Follow direction of negative curvature to the boundary of the
    297:          trust region.                                                       */
    298:       step       = PetscSqrtReal(r2 / norm_p);
    299:       cg->norm_d = cg->radius;

    301:       PetscCall(VecAXPY(d, step, p)); /* d = d + step p    */

    303:       /* Update objective function.                                          */
    304:       cg->o_fcn += step * (0.5 * step * kappa - rz);
    305:     } else if (cg->radius) {
    306:       /* The norm of the preconditioned direction is zero; use the gradient
    307:          step.                                                               */
    308:       if (r2 >= rr) {
    309:         alpha      = 1.0;
    310:         cg->norm_d = PetscSqrtReal(rr);
    311:       } else {
    312:         alpha      = PetscSqrtReal(r2 / rr);
    313:         cg->norm_d = cg->radius;
    314:       }

    316:       PetscCall(VecAXPY(d, alpha, r)); /* d = d + alpha r   */

    318:       /* Compute objective function.                                         */
    319:       PetscCall(KSP_MatMult(ksp, Qmat, d, z));
    320:       PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
    321:       PetscCall(VecDot(d, z, &cg->o_fcn));
    322:       cg->o_fcn = -cg->o_fcn;
    323:       ++ksp->its;
    324:     }
    325:     PetscFunctionReturn(PETSC_SUCCESS);
    326:   }

    328:   /*
    329:      Begin the first part of the GLTR algorithm which runs the conjugate
    330:      gradient method until either the problem is solved, we encounter the
    331:      boundary of the trust region, or the conjugate gradient method breaks
    332:      down.
    333:   */
    334:   while (1) {
    335:     /* Know that kappa is nonzero, because we have not broken down, so we    */
    336:     /* can compute the steplength.                                           */
    337:     alpha             = rz / kappa;
    338:     cg->alpha[l_size] = alpha;

    340:     /* Compute the diagonal value of the Cholesky factorization of the       */
    341:     /* Lanczos matrix and check to see if the Lanczos matrix is indefinite.  */
    342:     /* This indicates a direction of negative curvature.                     */
    343:     piv = cg->diag[l_size] - cg->offd[l_size] * cg->offd[l_size] / piv;
    344:     if (piv <= 0.0) {
    345:       /* In this case, the matrix is indefinite and we have encountered      */
    346:       /* a direction of negative curvature.  Follow the direction to the     */
    347:       /* boundary of the trust region.                                       */
    348:       ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
    349:       PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: negative curvature: kappa=%g\n", (double)kappa));

    351:       if (cg->radius && norm_p > 0.0) {
    352:         /* Follow direction of negative curvature to boundary.               */
    353:         step       = (PetscSqrtReal(dMp * dMp + norm_p * (r2 - norm_d)) - dMp) / norm_p;
    354:         cg->norm_d = cg->radius;

    356:         PetscCall(VecAXPY(d, step, p)); /* d = d + step p    */

    358:         /* Update objective function.                                        */
    359:         cg->o_fcn += step * (0.5 * step * kappa - rz);
    360:       } else if (cg->radius) {
    361:         /* The norm of the direction is zero; there is nothing to follow.    */
    362:       }
    363:       break;
    364:     }

    366:     /* Compute the steplength and check for intersection with the trust      */
    367:     /* region.                                                               */
    368:     norm_dp1 = norm_d + alpha * (2.0 * dMp + alpha * norm_p);
    369:     if (cg->radius && norm_dp1 >= r2) {
    370:       /* In this case, the matrix is positive definite as far as we know.    */
    371:       /* However, the full step goes beyond the trust region.                */
    372:       ksp->reason = KSP_CONVERGED_STEP_LENGTH;
    373:       PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: constrained step: radius=%g\n", (double)cg->radius));

    375:       if (norm_p > 0.0) {
    376:         /* Follow the direction to the boundary of the trust region.         */

    378:         step       = (PetscSqrtReal(dMp * dMp + norm_p * (r2 - norm_d)) - dMp) / norm_p;
    379:         cg->norm_d = cg->radius;

    381:         PetscCall(VecAXPY(d, step, p)); /* d = d + step p    */

    383:         /* Update objective function.                                        */
    384:         cg->o_fcn += step * (0.5 * step * kappa - rz);
    385:       } else {
    386:         /* The norm of the direction is zero; there is nothing to follow.    */
    387:       }
    388:       break;
    389:     }

    391:     /* Now we can update the direction and residual.                         */
    392:     PetscCall(VecAXPY(d, alpha, p));   /* d = d + alpha p   */
    393:     PetscCall(VecAXPY(r, -alpha, z));  /* r = r - alpha Q p */
    394:     PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r      */

    396:     switch (cg->dtype) {
    397:     case GLTR_PRECONDITIONED_DIRECTION:
    398:       norm_d = norm_dp1;
    399:       break;

    401:     default:
    402:       PetscCall(VecDot(d, d, &norm_d));
    403:       break;
    404:     }
    405:     cg->norm_d = PetscSqrtReal(norm_d);

    407:     /* Update objective function.                                            */
    408:     cg->o_fcn -= 0.5 * alpha * rz;

    410:     /* Check that the preconditioner appears positive semidefinite.          */
    411:     rzm1 = rz;
    412:     PetscCall(VecDot(r, z, &rz)); /* rz = r^T z        */
    413:     if (rz < 0.0) {
    414:       /* The preconditioner is indefinite.                                   */
    415:       ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
    416:       PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: cg indefinite preconditioner: rz=%g\n", (double)rz));
    417:       break;
    418:     }

    420:     /* As far as we know, the preconditioner is positive semidefinite.       */
    421:     /* Compute the residual and check for convergence.                       */
    422:     cg->norm_r[l_size + 1] = PetscSqrtReal(rz); /* norm_r = |r|_M   */

    424:     switch (ksp->normtype) {
    425:     case KSP_NORM_PRECONDITIONED:
    426:       PetscCall(VecNorm(z, NORM_2, &norm_r)); /* norm_r = |z|      */
    427:       break;

    429:     case KSP_NORM_UNPRECONDITIONED:
    430:       PetscCall(VecNorm(r, NORM_2, &norm_r)); /* norm_r = |r|      */
    431:       break;

    433:     case KSP_NORM_NATURAL:
    434:       norm_r = cg->norm_r[l_size + 1]; /* norm_r = |r|_M    */
    435:       break;

    437:     default:
    438:       norm_r = 0.0;
    439:       break;
    440:     }

    442:     PetscCall(KSPLogResidualHistory(ksp, norm_r));
    443:     PetscCall(KSPMonitor(ksp, ksp->its, norm_r));
    444:     ksp->rnorm = norm_r;

    446:     PetscCall((*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP));
    447:     if (ksp->reason) {
    448:       /* The method has converged.                                           */
    449:       PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: cg truncated step: rnorm=%g, radius=%g\n", (double)norm_r, (double)cg->radius));
    450:       break;
    451:     }

    453:     /* We have not converged yet.  Check for breakdown.                      */
    454:     beta = rz / rzm1;
    455:     if (PetscAbsReal(beta) <= 0.0) {
    456:       /* Conjugate gradients has broken down.                                */
    457:       ksp->reason = KSP_DIVERGED_BREAKDOWN;
    458:       PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: breakdown: beta=%g\n", (double)beta));
    459:       break;
    460:     }

    462:     /* Check iteration limit.                                                */
    463:     if (ksp->its >= max_cg_its) {
    464:       ksp->reason = KSP_DIVERGED_ITS;
    465:       PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: iterlim: its=%" PetscInt_FMT "\n", ksp->its));
    466:       break;
    467:     }

    469:     /* Update p and the norms.                                               */
    470:     cg->beta[l_size] = beta;
    471:     PetscCall(VecAYPX(p, beta, z)); /* p = z + beta p    */

    473:     switch (cg->dtype) {
    474:     case GLTR_PRECONDITIONED_DIRECTION:
    475:       dMp    = beta * (dMp + alpha * norm_p);
    476:       norm_p = beta * (rzm1 + beta * norm_p);
    477:       break;

    479:     default:
    480:       PetscCall(VecDot(d, p, &dMp));
    481:       PetscCall(VecDot(p, p, &norm_p));
    482:       break;
    483:     }

    485:     /* Compute the new direction and update the iteration.                   */
    486:     PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p         */
    487:     PetscCall(VecDot(p, z, &kappa));         /* kappa = p^T Q p   */
    488:     ++ksp->its;

    490:     /* Update the Lanczos tridiagonal matrix.                            */
    491:     ++l_size;
    492:     cg->offd[t_size] = PetscSqrtReal(beta) / PetscAbsReal(alpha);
    493:     cg->diag[t_size] = kappa / rz + beta / alpha;
    494:     ++t_size;

    496:     /* Check for breakdown of the conjugate gradient method; this occurs     */
    497:     /* when kappa is zero.                                                   */
    498:     if (PetscAbsReal(kappa) <= 0.0) {
    499:       /* The method breaks down; move along the direction as if the matrix   */
    500:       /* were indefinite.                                                    */
    501:       ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
    502:       PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: cg breakdown: kappa=%g\n", (double)kappa));

    504:       if (cg->radius && norm_p > 0.0) {
    505:         /* Follow direction to boundary.                                     */
    506:         step       = (PetscSqrtReal(dMp * dMp + norm_p * (r2 - norm_d)) - dMp) / norm_p;
    507:         cg->norm_d = cg->radius;

    509:         PetscCall(VecAXPY(d, step, p)); /* d = d + step p    */

    511:         /* Update objective function.                                        */
    512:         cg->o_fcn += step * (0.5 * step * kappa - rz);
    513:       } else if (cg->radius) {
    514:         /* The norm of the direction is zero; there is nothing to follow.    */
    515:       }
    516:       break;
    517:     }
    518:   }

    520:   /* Check to see if we need to continue with the Lanczos method.            */
    521:   if (!cg->radius) {
    522:     /* There is no radius.  Therefore, we cannot move along the boundary.    */
    523:     PetscFunctionReturn(PETSC_SUCCESS);
    524:   }

    526:   if (KSP_CONVERGED_NEG_CURVE != ksp->reason) {
    527:     /* The method either converged to an interior point, hit the boundary of */
    528:     /* the trust-region without encountering a direction of negative         */
    529:     /* curvature or the iteration limit was reached.                         */
    530:     PetscFunctionReturn(PETSC_SUCCESS);
    531:   }

    533:   /* Switch to constructing the Lanczos basis by way of the conjugate        */
    534:   /* directions.                                                             */
    535:   for (i = 0; i < max_lanczos_its; ++i) {
    536:     /* Check for breakdown of the conjugate gradient method; this occurs     */
    537:     /* when kappa is zero.                                                   */
    538:     if (PetscAbsReal(kappa) <= 0.0) {
    539:       ksp->reason = KSP_DIVERGED_BREAKDOWN;
    540:       PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: lanczos breakdown: kappa=%g\n", (double)kappa));
    541:       break;
    542:     }

    544:     /* Update the direction and residual.                                    */
    545:     alpha             = rz / kappa;
    546:     cg->alpha[l_size] = alpha;

    548:     PetscCall(VecAXPY(r, -alpha, z));  /* r = r - alpha Q p */
    549:     PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r      */

    551:     /* Check that the preconditioner appears positive semidefinite.          */
    552:     rzm1 = rz;
    553:     PetscCall(VecDot(r, z, &rz)); /* rz = r^T z        */
    554:     if (rz < 0.0) {
    555:       /* The preconditioner is indefinite.                                   */
    556:       ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
    557:       PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: lanczos indefinite preconditioner: rz=%g\n", (double)rz));
    558:       break;
    559:     }

    561:     /* As far as we know, the preconditioner is positive definite.  Compute  */
    562:     /* the residual.  Do NOT check for convergence.                          */
    563:     cg->norm_r[l_size + 1] = PetscSqrtReal(rz); /* norm_r = |r|_M    */

    565:     switch (ksp->normtype) {
    566:     case KSP_NORM_PRECONDITIONED:
    567:       PetscCall(VecNorm(z, NORM_2, &norm_r)); /* norm_r = |z|      */
    568:       break;

    570:     case KSP_NORM_UNPRECONDITIONED:
    571:       PetscCall(VecNorm(r, NORM_2, &norm_r)); /* norm_r = |r|      */
    572:       break;

    574:     case KSP_NORM_NATURAL:
    575:       norm_r = cg->norm_r[l_size + 1]; /* norm_r = |r|_M    */
    576:       break;

    578:     default:
    579:       norm_r = 0.0;
    580:       break;
    581:     }

    583:     PetscCall(KSPLogResidualHistory(ksp, norm_r));
    584:     PetscCall(KSPMonitor(ksp, ksp->its, norm_r));
    585:     ksp->rnorm = norm_r;

    587:     /* Check for breakdown.                                                  */
    588:     beta = rz / rzm1;
    589:     if (PetscAbsReal(beta) <= 0.0) {
    590:       /* Conjugate gradients has broken down.                                */
    591:       ksp->reason = KSP_DIVERGED_BREAKDOWN;
    592:       PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: breakdown: beta=%g\n", (double)beta));
    593:       break;
    594:     }

    596:     /* Update p and the norms.                                               */
    597:     cg->beta[l_size] = beta;
    598:     PetscCall(VecAYPX(p, beta, z)); /* p = z + beta p    */

    600:     /* Compute the new direction and update the iteration.                   */
    601:     PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p         */
    602:     PetscCall(VecDot(p, z, &kappa));         /* kappa = p^T Q p   */
    603:     ++ksp->its;

    605:     /* Update the Lanczos tridiagonal matrix.                                */
    606:     ++l_size;
    607:     cg->offd[t_size] = PetscSqrtReal(beta) / PetscAbsReal(alpha);
    608:     cg->diag[t_size] = kappa / rz + beta / alpha;
    609:     ++t_size;
    610:   }

    612:   /*
    613:     We have the Lanczos basis, solve the tridiagonal trust-region problem
    614:     to obtain the Lanczos direction.  We know that the solution lies on
    615:     the boundary of the trust region.  We start by checking that the
    616:     workspace allocated is large enough.

    618:     Note that the current version only computes the solution by using the
    619:     preconditioned direction.  Need to think about how to do the
    620:     unpreconditioned direction calculation.
    621:   */

    623:   if (t_size > cg->alloced) {
    624:     if (cg->alloced) {
    625:       PetscCall(PetscFree(cg->rwork));
    626:       PetscCall(PetscFree(cg->iwork));
    627:       cg->alloced += cg->init_alloc;
    628:     } else {
    629:       cg->alloced = cg->init_alloc;
    630:     }

    632:     while (t_size > cg->alloced) cg->alloced += cg->init_alloc;

    634:     cg->alloced = PetscMin(cg->alloced, t_size);
    635:     PetscCall(PetscMalloc2(10 * cg->alloced, &cg->rwork, 5 * cg->alloced, &cg->iwork));
    636:   }

    638:   /* Set up the required vectors.                                            */
    639:   t_soln = cg->rwork + 0 * t_size; /* Solution          */
    640:   t_diag = cg->rwork + 1 * t_size; /* Diagonal of T     */
    641:   t_offd = cg->rwork + 2 * t_size; /* Off-diagonal of T */
    642:   e_valu = cg->rwork + 3 * t_size; /* Eigenvalues of T  */
    643:   e_vect = cg->rwork + 4 * t_size; /* Eigenvector of T  */
    644:   e_rwrk = cg->rwork + 5 * t_size; /* Eigen workspace   */

    646:   e_iblk = cg->iwork + 0 * t_size; /* Eigen blocks      */
    647:   e_splt = cg->iwork + 1 * t_size; /* Eigen splits      */
    648:   e_iwrk = cg->iwork + 2 * t_size; /* Eigen workspace   */

    650:   /* Compute the minimum eigenvalue of T.                                    */
    651:   vl = 0.0;
    652:   vu = 0.0;
    653:   il = 1;
    654:   iu = 1;

    656:   PetscCallBLAS("LAPACKstebz", LAPACKstebz_("I", "E", &t_size, &vl, &vu, &il, &iu, &cg->eigen_tol, cg->diag, cg->offd + 1, &e_valus, &e_splts, e_valu, e_iblk, e_splt, e_rwrk, e_iwrk, &info));

    658:   if ((0 != info) || (1 != e_valus)) {
    659:     /* Calculation of the minimum eigenvalue failed.  Return the             */
    660:     /* Steihaug-Toint direction.                                             */
    661:     PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute eigenvalue.\n"));
    662:     ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
    663:     PetscFunctionReturn(PETSC_SUCCESS);
    664:   }

    666:   cg->e_min = e_valu[0];

    668:   /* Compute the initial value of lambda to make (T + lambda I) positive      */
    669:   /* definite.                                                               */
    670:   pert = cg->init_pert;
    671:   if (e_valu[0] < 0.0) cg->lambda = pert - e_valu[0];

    673:   while (1) {
    674:     for (i = 0; i < t_size; ++i) {
    675:       t_diag[i] = cg->diag[i] + cg->lambda;
    676:       t_offd[i] = cg->offd[i];
    677:     }

    679:     PetscCallBLAS("LAPACKpttrf", LAPACKpttrf_(&t_size, t_diag, t_offd + 1, &info));

    681:     if (0 == info) break;

    683:     pert += pert;
    684:     cg->lambda = cg->lambda * (1.0 + pert) + pert;
    685:   }

    687:   /* Compute the initial step and its norm.                                  */
    688:   nrhs = 1;
    689:   nldb = t_size;

    691:   t_soln[0] = -cg->norm_r[0];
    692:   for (i = 1; i < t_size; ++i) t_soln[i] = 0.0;

    694:   PetscCallBLAS("LAPACKpttrs", LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, t_soln, &nldb, &info));

    696:   if (0 != info) {
    697:     /* Calculation of the initial step failed; return the Steihaug-Toint     */
    698:     /* direction.                                                            */
    699:     PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute step.\n"));
    700:     ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
    701:     PetscFunctionReturn(PETSC_SUCCESS);
    702:   }

    704:   norm_t = 0.;
    705:   for (i = 0; i < t_size; ++i) norm_t += t_soln[i] * t_soln[i];
    706:   norm_t = PetscSqrtReal(norm_t);

    708:   /* Determine the case we are in.                                           */
    709:   if (norm_t <= cg->radius) {
    710:     /* The step is within the trust region; check if we are in the hard case */
    711:     /* and need to move to the boundary by following a direction of negative */
    712:     /* curvature.                                                            */
    713:     if ((e_valu[0] <= 0.0) && (norm_t < cg->radius)) {
    714:       /* This is the hard case; compute the eigenvector associated with the  */
    715:       /* minimum eigenvalue and move along this direction to the boundary.   */
    716:       PetscCallBLAS("LAPACKstein", LAPACKstein_(&t_size, cg->diag, cg->offd + 1, &e_valus, e_valu, e_iblk, e_splt, e_vect, &nldb, e_rwrk, e_iwrk, e_iwrk + t_size, &info));

    718:       if (0 != info) {
    719:         /* Calculation of the minimum eigenvalue failed.  Return the         */
    720:         /* Steihaug-Toint direction.                                         */
    721:         PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute eigenvector.\n"));
    722:         ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
    723:         PetscFunctionReturn(PETSC_SUCCESS);
    724:       }

    726:       coef1 = 0.0;
    727:       coef2 = 0.0;
    728:       coef3 = -cg->radius * cg->radius;
    729:       for (i = 0; i < t_size; ++i) {
    730:         coef1 += e_vect[i] * e_vect[i];
    731:         coef2 += e_vect[i] * t_soln[i];
    732:         coef3 += t_soln[i] * t_soln[i];
    733:       }

    735:       coef3 = PetscSqrtReal(coef2 * coef2 - coef1 * coef3);
    736:       root1 = (-coef2 + coef3) / coef1;
    737:       root2 = (-coef2 - coef3) / coef1;

    739:       /* Compute objective value for (t_soln + root1 * e_vect)               */
    740:       for (i = 0; i < t_size; ++i) e_rwrk[i] = t_soln[i] + root1 * e_vect[i];

    742:       obj1 = e_rwrk[0] * (0.5 * (cg->diag[0] * e_rwrk[0] + cg->offd[1] * e_rwrk[1]) + cg->norm_r[0]);
    743:       for (i = 1; i < t_size - 1; ++i) obj1 += 0.5 * e_rwrk[i] * (cg->offd[i] * e_rwrk[i - 1] + cg->diag[i] * e_rwrk[i] + cg->offd[i + 1] * e_rwrk[i + 1]);
    744:       obj1 += 0.5 * e_rwrk[i] * (cg->offd[i] * e_rwrk[i - 1] + cg->diag[i] * e_rwrk[i]);

    746:       /* Compute objective value for (t_soln + root2 * e_vect)               */
    747:       for (i = 0; i < t_size; ++i) e_rwrk[i] = t_soln[i] + root2 * e_vect[i];

    749:       obj2 = e_rwrk[0] * (0.5 * (cg->diag[0] * e_rwrk[0] + cg->offd[1] * e_rwrk[1]) + cg->norm_r[0]);
    750:       for (i = 1; i < t_size - 1; ++i) obj2 += 0.5 * e_rwrk[i] * (cg->offd[i] * e_rwrk[i - 1] + cg->diag[i] * e_rwrk[i] + cg->offd[i + 1] * e_rwrk[i + 1]);
    751:       obj2 += 0.5 * e_rwrk[i] * (cg->offd[i] * e_rwrk[i - 1] + cg->diag[i] * e_rwrk[i]);

    753:       /* Choose the point with the best objective function value.            */
    754:       if (obj1 <= obj2) {
    755:         for (i = 0; i < t_size; ++i) t_soln[i] += root1 * e_vect[i];
    756:       } else {
    757:         for (i = 0; i < t_size; ++i) t_soln[i] += root2 * e_vect[i];
    758:       }
    759:     } else {
    760:       /* The matrix is positive definite or there was no room to move; the   */
    761:       /* solution is already contained in t_soln.                            */
    762:     }
    763:   } else {
    764:     /* The step is outside the trust-region.  Compute the correct value for  */
    765:     /* lambda by performing Newton's method.                                 */

    767:     for (i = 0; i < max_newton_its; ++i) {
    768:       /* Check for convergence.                                              */
    769:       if (PetscAbsReal(norm_t - cg->radius) <= cg->newton_tol * cg->radius) break;

    771:       /* Compute the update.                                                 */
    772:       PetscCall(PetscArraycpy(e_rwrk, t_soln, t_size));

    774:       PetscCallBLAS("LAPACKpttrs", LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, e_rwrk, &nldb, &info));

    776:       if (0 != info) {
    777:         /* Calculation of the step failed; return the Steihaug-Toint         */
    778:         /* direction.                                                        */
    779:         PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute step.\n"));
    780:         ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
    781:         PetscFunctionReturn(PETSC_SUCCESS);
    782:       }

    784:       /* Modify lambda.                                                      */
    785:       norm_w = 0.;
    786:       for (j = 0; j < t_size; ++j) norm_w += t_soln[j] * e_rwrk[j];

    788:       cg->lambda += (norm_t - cg->radius) / cg->radius * (norm_t * norm_t) / norm_w;

    790:       /* Factor T + lambda I                                                 */
    791:       for (j = 0; j < t_size; ++j) {
    792:         t_diag[j] = cg->diag[j] + cg->lambda;
    793:         t_offd[j] = cg->offd[j];
    794:       }

    796:       PetscCallBLAS("LAPACKpttrf", LAPACKpttrf_(&t_size, t_diag, t_offd + 1, &info));

    798:       if (0 != info) {
    799:         /* Calculation of factorization failed; return the Steihaug-Toint    */
    800:         /* direction.                                                        */
    801:         PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: factorization failed.\n"));
    802:         ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
    803:         PetscFunctionReturn(PETSC_SUCCESS);
    804:       }

    806:       /* Compute the new step and its norm.                                  */
    807:       t_soln[0] = -cg->norm_r[0];
    808:       for (j = 1; j < t_size; ++j) t_soln[j] = 0.0;

    810:       PetscCallBLAS("LAPACKpttrs", LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, t_soln, &nldb, &info));

    812:       if (0 != info) {
    813:         /* Calculation of the step failed; return the Steihaug-Toint         */
    814:         /* direction.                                                        */
    815:         PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute step.\n"));
    816:         ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
    817:         PetscFunctionReturn(PETSC_SUCCESS);
    818:       }

    820:       norm_t = 0.;
    821:       for (j = 0; j < t_size; ++j) norm_t += t_soln[j] * t_soln[j];
    822:       norm_t = PetscSqrtReal(norm_t);
    823:     }

    825:     /* Check for convergence.                                                */
    826:     if (PetscAbsReal(norm_t - cg->radius) > cg->newton_tol * cg->radius) {
    827:       /* Newton method failed to converge in iteration limit.                */
    828:       PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to converge.\n"));
    829:       ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
    830:       PetscFunctionReturn(PETSC_SUCCESS);
    831:     }
    832:   }

    834:   /* Recover the norm of the direction and objective function value.         */
    835:   cg->norm_d = norm_t;

    837:   cg->o_fcn = t_soln[0] * (0.5 * (cg->diag[0] * t_soln[0] + cg->offd[1] * t_soln[1]) + cg->norm_r[0]);
    838:   for (i = 1; i < t_size - 1; ++i) cg->o_fcn += 0.5 * t_soln[i] * (cg->offd[i] * t_soln[i - 1] + cg->diag[i] * t_soln[i] + cg->offd[i + 1] * t_soln[i + 1]);
    839:   cg->o_fcn += 0.5 * t_soln[i] * (cg->offd[i] * t_soln[i - 1] + cg->diag[i] * t_soln[i]);

    841:   /* Recover the direction.                                                  */
    842:   sigma = -1;

    844:   /* Start conjugate gradient method from the beginning                      */
    845:   PetscCall(VecCopy(ksp->vec_rhs, r)); /* r = -grad         */
    846:   PetscCall(KSP_PCApply(ksp, r, z));   /* z = inv(M) r      */

    848:   /* Accumulate Q * s                                                        */
    849:   PetscCall(VecCopy(z, d));
    850:   PetscCall(VecScale(d, sigma * t_soln[0] / cg->norm_r[0]));

    852:   /* Compute the first direction.                                            */
    853:   PetscCall(VecCopy(z, p));                /* p = z             */
    854:   PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p         */
    855:   ++ksp->its;

    857:   for (i = 0; i < l_size - 1; ++i) {
    858:     /* Update the residual and direction.                                    */
    859:     alpha = cg->alpha[i];
    860:     if (alpha >= 0.0) sigma = -sigma;

    862:     PetscCall(VecAXPY(r, -alpha, z));  /* r = r - alpha Q p */
    863:     PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r      */

    865:     /* Accumulate Q * s                                                      */
    866:     PetscCall(VecAXPY(d, sigma * t_soln[i + 1] / cg->norm_r[i + 1], z));

    868:     /* Update p.                                                             */
    869:     beta = cg->beta[i];
    870:     PetscCall(VecAYPX(p, beta, z));          /* p = z + beta p    */
    871:     PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p         */
    872:     ++ksp->its;
    873:   }

    875:   /* Update the residual and direction.                                      */
    876:   alpha = cg->alpha[i];
    877:   if (alpha >= 0.0) sigma = -sigma;

    879:   PetscCall(VecAXPY(r, -alpha, z));  /* r = r - alpha Q p */
    880:   PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r      */

    882:   /* Accumulate Q * s                                                        */
    883:   PetscCall(VecAXPY(d, sigma * t_soln[i + 1] / cg->norm_r[i + 1], z));

    885:   /* Set the termination reason.                                             */
    886:   ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
    887:   PetscFunctionReturn(PETSC_SUCCESS);
    888: #endif
    889: }

    891: static PetscErrorCode KSPCGSetUp_GLTR(KSP ksp)
    892: {
    893:   KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
    894:   PetscInt    max_its;

    896:   PetscFunctionBegin;
    897:   /* Determine the total maximum number of iterations.                       */
    898:   max_its = ksp->max_it + cg->max_lanczos_its + 1;

    900:   /* Set work vectors needed by conjugate gradient method and allocate       */
    901:   /* workspace for Lanczos matrix.                                           */
    902:   PetscCall(KSPSetWorkVecs(ksp, 3));
    903:   if (cg->diag) {
    904:     PetscCall(PetscArrayzero(cg->diag, max_its));
    905:     PetscCall(PetscArrayzero(cg->offd, max_its));
    906:     PetscCall(PetscArrayzero(cg->alpha, max_its));
    907:     PetscCall(PetscArrayzero(cg->beta, max_its));
    908:     PetscCall(PetscArrayzero(cg->norm_r, max_its));
    909:   } else {
    910:     PetscCall(PetscCalloc5(max_its, &cg->diag, max_its, &cg->offd, max_its, &cg->alpha, max_its, &cg->beta, max_its, &cg->norm_r));
    911:   }
    912:   PetscFunctionReturn(PETSC_SUCCESS);
    913: }

    915: static PetscErrorCode KSPCGDestroy_GLTR(KSP ksp)
    916: {
    917:   KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;

    919:   PetscFunctionBegin;
    920:   PetscCall(PetscFree5(cg->diag, cg->offd, cg->alpha, cg->beta, cg->norm_r));
    921:   if (cg->alloced) PetscCall(PetscFree2(cg->rwork, cg->iwork));
    922:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetRadius_C", NULL));
    923:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetNormD_C", NULL));
    924:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", NULL));
    925:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGLTRGetMinEig_C", NULL));
    926:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGLTRGetLambda_C", NULL));
    927:   PetscCall(KSPDestroyDefault(ksp));
    928:   PetscFunctionReturn(PETSC_SUCCESS);
    929: }

    931: static PetscErrorCode KSPCGSetRadius_GLTR(KSP ksp, PetscReal radius)
    932: {
    933:   KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;

    935:   PetscFunctionBegin;
    936:   cg->radius = radius;
    937:   PetscFunctionReturn(PETSC_SUCCESS);
    938: }

    940: static PetscErrorCode KSPCGGetNormD_GLTR(KSP ksp, PetscReal *norm_d)
    941: {
    942:   KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;

    944:   PetscFunctionBegin;
    945:   *norm_d = cg->norm_d;
    946:   PetscFunctionReturn(PETSC_SUCCESS);
    947: }

    949: static PetscErrorCode KSPCGGetObjFcn_GLTR(KSP ksp, PetscReal *o_fcn)
    950: {
    951:   KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;

    953:   PetscFunctionBegin;
    954:   *o_fcn = cg->o_fcn;
    955:   PetscFunctionReturn(PETSC_SUCCESS);
    956: }

    958: static PetscErrorCode KSPGLTRGetMinEig_GLTR(KSP ksp, PetscReal *e_min)
    959: {
    960:   KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;

    962:   PetscFunctionBegin;
    963:   *e_min = cg->e_min;
    964:   PetscFunctionReturn(PETSC_SUCCESS);
    965: }

    967: static PetscErrorCode KSPGLTRGetLambda_GLTR(KSP ksp, PetscReal *lambda)
    968: {
    969:   KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;

    971:   PetscFunctionBegin;
    972:   *lambda = cg->lambda;
    973:   PetscFunctionReturn(PETSC_SUCCESS);
    974: }

    976: static PetscErrorCode KSPCGSetFromOptions_GLTR(KSP ksp, PetscOptionItems *PetscOptionsObject)
    977: {
    978:   KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;

    980:   PetscFunctionBegin;
    981:   PetscOptionsHeadBegin(PetscOptionsObject, "KSP GLTR options");

    983:   PetscCall(PetscOptionsReal("-ksp_cg_radius", "Trust Region Radius", "KSPCGSetRadius", cg->radius, &cg->radius, NULL));

    985:   PetscCall(PetscOptionsEList("-ksp_cg_dtype", "Norm used for direction", "", DType_Table, GLTR_DIRECTION_TYPES, DType_Table[cg->dtype], &cg->dtype, NULL));

    987:   PetscCall(PetscOptionsReal("-ksp_cg_gltr_init_pert", "Initial perturbation", "", cg->init_pert, &cg->init_pert, NULL));
    988:   PetscCall(PetscOptionsReal("-ksp_cg_gltr_eigen_tol", "Eigenvalue tolerance", "", cg->eigen_tol, &cg->eigen_tol, NULL));
    989:   PetscCall(PetscOptionsReal("-ksp_cg_gltr_newton_tol", "Newton tolerance", "", cg->newton_tol, &cg->newton_tol, NULL));

    991:   PetscCall(PetscOptionsInt("-ksp_cg_gltr_max_lanczos_its", "Maximum Lanczos Iters", "", cg->max_lanczos_its, &cg->max_lanczos_its, NULL));
    992:   PetscCall(PetscOptionsInt("-ksp_cg_gltr_max_newton_its", "Maximum Newton Iters", "", cg->max_newton_its, &cg->max_newton_its, NULL));

    994:   PetscOptionsHeadEnd();
    995:   PetscFunctionReturn(PETSC_SUCCESS);
    996: }

    998: /*MC
    999:    KSPGLTR -   Code to run conjugate gradient method subject to a constraint on the solution norm, used within trust region methods {cite}`gould1999solving`

    1001:    Options Database Key:
    1002: .  -ksp_cg_radius <r> - Trust Region Radius

    1004:    Level: developer

    1006:    Notes:
    1007:    Uses preconditioned conjugate gradient to compute  an approximate minimizer of the quadratic function

    1009:    $$
    1010:    q(s) = g^T * s + .5 * s^T * H * s
    1011:    $$

    1013:    subject to the trust region constraint

    1015:    $$
    1016:    || s || \le delta,
    1017:    $$

    1019:    where
    1020: .vb
    1021:      delta is the trust region radius,
    1022:      g is the gradient vector,
    1023:      H is the Hessian approximation,
    1024:      M is the positive definite preconditioner matrix.
    1025: .ve

    1027:    `KSPConvergedReason` may have the additional values
    1028: +  `KSP_CONVERGED_NEG_CURVE` - if convergence is reached along a negative curvature direction,
    1029: -  `KSP_CONVERGED_STEP_LENGTH` - if convergence is reached along a constrained step.

    1031:   The operator and the preconditioner supplied must be symmetric and positive definite.

    1033:   This is rarely used directly, it is used in Trust Region methods for nonlinear equations, `SNESNEWTONTR`

    1035: .seealso: [](ch_ksp), `KSPQCG`, `KSPNASH`, `KSPSTCG`, `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPCGSetRadius()`, `KSPCGGetNormD()`, `KSPCGGetObjFcn()`, `KSPGLTRGetMinEig()`, `KSPGLTRGetLambda()`, `KSPCG`
    1036: M*/

    1038: PETSC_EXTERN PetscErrorCode KSPCreate_GLTR(KSP ksp)
    1039: {
    1040:   KSPCG_GLTR *cg;

    1042:   PetscFunctionBegin;
    1043:   PetscCall(PetscNew(&cg));
    1044:   cg->radius = 0.0;
    1045:   cg->dtype  = GLTR_UNPRECONDITIONED_DIRECTION;

    1047:   cg->init_pert  = 1.0e-8;
    1048:   cg->eigen_tol  = 1.0e-10;
    1049:   cg->newton_tol = 1.0e-6;

    1051:   cg->alloced    = 0;
    1052:   cg->init_alloc = 1024;

    1054:   cg->max_lanczos_its = 20;
    1055:   cg->max_newton_its  = 10;

    1057:   ksp->data = (void *)cg;
    1058:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 3));
    1059:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 2));
    1060:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 2));
    1061:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
    1062:   PetscCall(KSPSetConvergedNegativeCurvature(ksp, PETSC_TRUE));

    1064:   /* Sets the functions that are associated with this data structure         */
    1065:   /* (in C++ this is the same as defining virtual functions).                */

    1067:   ksp->ops->setup          = KSPCGSetUp_GLTR;
    1068:   ksp->ops->solve          = KSPCGSolve_GLTR;
    1069:   ksp->ops->destroy        = KSPCGDestroy_GLTR;
    1070:   ksp->ops->setfromoptions = KSPCGSetFromOptions_GLTR;
    1071:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
    1072:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
    1073:   ksp->ops->view           = NULL;

    1075:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetRadius_C", KSPCGSetRadius_GLTR));
    1076:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetNormD_C", KSPCGGetNormD_GLTR));
    1077:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", KSPCGGetObjFcn_GLTR));
    1078:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGLTRGetMinEig_C", KSPGLTRGetMinEig_GLTR));
    1079:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGLTRGetLambda_C", KSPGLTRGetLambda_GLTR));
    1080:   PetscFunctionReturn(PETSC_SUCCESS);
    1081: }