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: }