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:   /***************************************************************************/
 85:   /* Check the arguments and parameters.                                     */
 86:   /***************************************************************************/

 88:   PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
 89:   PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
 90:   PetscCheck(cg->radius >= 0.0, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "Input error: radius < 0");

 92:   /***************************************************************************/
 93:   /* Get the workspace vectors and initialize variables                      */
 94:   /***************************************************************************/

 96:   r2 = cg->radius * cg->radius;
 97:   r  = ksp->work[0];
 98:   z  = ksp->work[1];
 99:   p  = ksp->work[2];
100:   d  = ksp->vec_sol;
101:   pc = ksp->pc;

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

105:   PetscCall(VecGetSize(d, &max_cg_its));
106:   max_cg_its      = PetscMin(max_cg_its, ksp->max_it);
107:   max_lanczos_its = cg->max_lanczos_its;
108:   max_newton_its  = cg->max_newton_its;
109:   ksp->its        = 0;

111:   /***************************************************************************/
112:   /* Initialize objective function direction, and minimum eigenvalue.        */
113:   /***************************************************************************/

115:   cg->o_fcn = 0.0;

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

120:   cg->e_min  = 0.0;
121:   cg->lambda = 0.0;

123:   /***************************************************************************/
124:   /* The first phase of GLTR performs a standard conjugate gradient method,  */
125:   /* but stores the values required for the Lanczos matrix.  We switch to    */
126:   /* the Lanczos when the conjugate gradient method breaks down.  Check the  */
127:   /* right-hand side for numerical problems.  The check for not-a-number and */
128:   /* infinite values need be performed only once.                            */
129:   /***************************************************************************/

131:   PetscCall(VecCopy(ksp->vec_rhs, r)); /* r = -grad         */
132:   PetscCall(VecDot(r, r, &rr));        /* rr = r^T r        */
133:   KSPCheckDot(ksp, rr);

135:   /***************************************************************************/
136:   /* Check the preconditioner for numerical problems and for positive        */
137:   /* definiteness.  The check for not-a-number and infinite values need be   */
138:   /* performed only once.                                                    */
139:   /***************************************************************************/

141:   PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r      */
142:   PetscCall(VecDot(r, z, &rz));      /* rz = r^T inv(M) r */
143:   if (PetscIsInfOrNanScalar(rz)) {
144:     /*************************************************************************/
145:     /* The preconditioner contains not-a-number or an infinite value.        */
146:     /* Return the gradient direction intersected with the trust region.      */
147:     /*************************************************************************/

149:     ksp->reason = KSP_DIVERGED_NANORINF;
150:     PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: bad preconditioner: rz=%g\n", (double)rz));

152:     if (cg->radius) {
153:       if (r2 >= rr) {
154:         alpha      = 1.0;
155:         cg->norm_d = PetscSqrtReal(rr);
156:       } else {
157:         alpha      = PetscSqrtReal(r2 / rr);
158:         cg->norm_d = cg->radius;
159:       }

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

163:       /***********************************************************************/
164:       /* Compute objective function.                                         */
165:       /***********************************************************************/

167:       PetscCall(KSP_MatMult(ksp, Qmat, d, z));
168:       PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
169:       PetscCall(VecDot(d, z, &cg->o_fcn));
170:       cg->o_fcn = -cg->o_fcn;
171:       ++ksp->its;
172:     }
173:     PetscFunctionReturn(PETSC_SUCCESS);
174:   }

176:   if (rz < 0.0) {
177:     /*************************************************************************/
178:     /* The preconditioner is indefinite.  Because this is the first          */
179:     /* and we do not have a direction yet, we use the gradient step.  Note   */
180:     /* that we cannot use the preconditioned norm when computing the step    */
181:     /* because the matrix is indefinite.                                     */
182:     /*************************************************************************/

184:     ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
185:     PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: indefinite preconditioner: rz=%g\n", (double)rz));

187:     if (cg->radius) {
188:       if (r2 >= rr) {
189:         alpha      = 1.0;
190:         cg->norm_d = PetscSqrtReal(rr);
191:       } else {
192:         alpha      = PetscSqrtReal(r2 / rr);
193:         cg->norm_d = cg->radius;
194:       }

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

198:       /***********************************************************************/
199:       /* Compute objective function.                                         */
200:       /***********************************************************************/

202:       PetscCall(KSP_MatMult(ksp, Qmat, d, z));
203:       PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
204:       PetscCall(VecDot(d, z, &cg->o_fcn));
205:       cg->o_fcn = -cg->o_fcn;
206:       ++ksp->its;
207:     }
208:     PetscFunctionReturn(PETSC_SUCCESS);
209:   }

211:   /***************************************************************************/
212:   /* As far as we know, the preconditioner is positive semidefinite.         */
213:   /* Compute and log the residual.  Check convergence because this           */
214:   /* initializes things, but do not terminate until at least one conjugate   */
215:   /* gradient iteration has been performed.                                  */
216:   /***************************************************************************/

218:   cg->norm_r[0] = PetscSqrtReal(rz); /* norm_r = |r|_M    */

220:   switch (ksp->normtype) {
221:   case KSP_NORM_PRECONDITIONED:
222:     PetscCall(VecNorm(z, NORM_2, &norm_r)); /* norm_r = |z|      */
223:     break;

225:   case KSP_NORM_UNPRECONDITIONED:
226:     norm_r = PetscSqrtReal(rr); /* norm_r = |r|      */
227:     break;

229:   case KSP_NORM_NATURAL:
230:     norm_r = cg->norm_r[0]; /* norm_r = |r|_M    */
231:     break;

233:   default:
234:     norm_r = 0.0;
235:     break;
236:   }

238:   PetscCall(KSPLogResidualHistory(ksp, norm_r));
239:   PetscCall(KSPMonitor(ksp, ksp->its, norm_r));
240:   ksp->rnorm = norm_r;

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

244:   /***************************************************************************/
245:   /* Compute the first direction and update the iteration.                   */
246:   /***************************************************************************/

248:   PetscCall(VecCopy(z, p));                /* p = z             */
249:   PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p         */
250:   ++ksp->its;

252:   /***************************************************************************/
253:   /* Check the matrix for numerical problems.                                */
254:   /***************************************************************************/

256:   PetscCall(VecDot(p, z, &kappa)); /* kappa = p^T Q p   */
257:   if (PetscIsInfOrNanScalar(kappa)) {
258:     /*************************************************************************/
259:     /* The matrix produced not-a-number or an infinite value.  In this case, */
260:     /* we must stop and use the gradient direction.  This condition need     */
261:     /* only be checked once.                                                 */
262:     /*************************************************************************/

264:     ksp->reason = KSP_DIVERGED_NANORINF;
265:     PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: bad matrix: kappa=%g\n", (double)kappa));

267:     if (cg->radius) {
268:       if (r2 >= rr) {
269:         alpha      = 1.0;
270:         cg->norm_d = PetscSqrtReal(rr);
271:       } else {
272:         alpha      = PetscSqrtReal(r2 / rr);
273:         cg->norm_d = cg->radius;
274:       }

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

278:       /***********************************************************************/
279:       /* Compute objective function.                                         */
280:       /***********************************************************************/

282:       PetscCall(KSP_MatMult(ksp, Qmat, d, z));
283:       PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
284:       PetscCall(VecDot(d, z, &cg->o_fcn));
285:       cg->o_fcn = -cg->o_fcn;
286:       ++ksp->its;
287:     }
288:     PetscFunctionReturn(PETSC_SUCCESS);
289:   }

291:   /***************************************************************************/
292:   /* Initialize variables for calculating the norm of the direction and for  */
293:   /* the Lanczos tridiagonal matrix.  Note that we track the diagonal value  */
294:   /* of the Cholesky factorization of the Lanczos matrix in order to         */
295:   /* determine when negative curvature is encountered.                       */
296:   /***************************************************************************/

298:   dMp    = 0.0;
299:   norm_d = 0.0;
300:   switch (cg->dtype) {
301:   case GLTR_PRECONDITIONED_DIRECTION:
302:     norm_p = rz;
303:     break;

305:   default:
306:     PetscCall(VecDot(p, p, &norm_p));
307:     break;
308:   }

310:   cg->diag[t_size] = kappa / rz;
311:   cg->offd[t_size] = 0.0;
312:   ++t_size;

314:   piv = 1.0;

316:   /***************************************************************************/
317:   /* Check for breakdown of the conjugate gradient method; this occurs when  */
318:   /* kappa is zero.                                                          */
319:   /***************************************************************************/

321:   if (PetscAbsReal(kappa) <= 0.0) {
322:     /*************************************************************************/
323:     /* The curvature is zero.  In this case, we must stop and use follow     */
324:     /* the direction of negative curvature since the Lanczos matrix is zero. */
325:     /*************************************************************************/

327:     ksp->reason = KSP_DIVERGED_BREAKDOWN;
328:     PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: breakdown: kappa=%g\n", (double)kappa));

330:     if (cg->radius && norm_p > 0.0) {
331:       /***********************************************************************/
332:       /* Follow direction of negative curvature to the boundary of the       */
333:       /* trust region.                                                       */
334:       /***********************************************************************/

336:       step       = PetscSqrtReal(r2 / norm_p);
337:       cg->norm_d = cg->radius;

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

341:       /***********************************************************************/
342:       /* Update objective function.                                          */
343:       /***********************************************************************/

345:       cg->o_fcn += step * (0.5 * step * kappa - rz);
346:     } else if (cg->radius) {
347:       /***********************************************************************/
348:       /* The norm of the preconditioned direction is zero; use the gradient  */
349:       /* step.                                                               */
350:       /***********************************************************************/

352:       if (r2 >= rr) {
353:         alpha      = 1.0;
354:         cg->norm_d = PetscSqrtReal(rr);
355:       } else {
356:         alpha      = PetscSqrtReal(r2 / rr);
357:         cg->norm_d = cg->radius;
358:       }

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

362:       /***********************************************************************/
363:       /* Compute objective function.                                         */
364:       /***********************************************************************/

366:       PetscCall(KSP_MatMult(ksp, Qmat, d, z));
367:       PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
368:       PetscCall(VecDot(d, z, &cg->o_fcn));
369:       cg->o_fcn = -cg->o_fcn;
370:       ++ksp->its;
371:     }
372:     PetscFunctionReturn(PETSC_SUCCESS);
373:   }

375:   /***************************************************************************/
376:   /* Begin the first part of the GLTR algorithm which runs the conjugate     */
377:   /* gradient method until either the problem is solved, we encounter the    */
378:   /* boundary of the trust region, or the conjugate gradient method breaks   */
379:   /* down.                                                                   */
380:   /***************************************************************************/

382:   while (1) {
383:     /*************************************************************************/
384:     /* Know that kappa is nonzero, because we have not broken down, so we    */
385:     /* can compute the steplength.                                           */
386:     /*************************************************************************/

388:     alpha             = rz / kappa;
389:     cg->alpha[l_size] = alpha;

391:     /*************************************************************************/
392:     /* Compute the diagonal value of the Cholesky factorization of the       */
393:     /* Lanczos matrix and check to see if the Lanczos matrix is indefinite.  */
394:     /* This indicates a direction of negative curvature.                     */
395:     /*************************************************************************/

397:     piv = cg->diag[l_size] - cg->offd[l_size] * cg->offd[l_size] / piv;
398:     if (piv <= 0.0) {
399:       /***********************************************************************/
400:       /* In this case, the matrix is indefinite and we have encountered      */
401:       /* a direction of negative curvature.  Follow the direction to the     */
402:       /* boundary of the trust region.                                       */
403:       /***********************************************************************/

405:       ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
406:       PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: negative curvature: kappa=%g\n", (double)kappa));

408:       if (cg->radius && norm_p > 0.0) {
409:         /*********************************************************************/
410:         /* Follow direction of negative curvature to boundary.               */
411:         /*********************************************************************/

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

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

418:         /*********************************************************************/
419:         /* Update objective function.                                        */
420:         /*********************************************************************/

422:         cg->o_fcn += step * (0.5 * step * kappa - rz);
423:       } else if (cg->radius) {
424:         /*********************************************************************/
425:         /* The norm of the direction is zero; there is nothing to follow.    */
426:         /*********************************************************************/
427:       }
428:       break;
429:     }

431:     /*************************************************************************/
432:     /* Compute the steplength and check for intersection with the trust      */
433:     /* region.                                                               */
434:     /*************************************************************************/

436:     norm_dp1 = norm_d + alpha * (2.0 * dMp + alpha * norm_p);
437:     if (cg->radius && norm_dp1 >= r2) {
438:       /***********************************************************************/
439:       /* In this case, the matrix is positive definite as far as we know.    */
440:       /* However, the full step goes beyond the trust region.                */
441:       /***********************************************************************/

443:       ksp->reason = KSP_CONVERGED_STEP_LENGTH;
444:       PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: constrained step: radius=%g\n", (double)cg->radius));

446:       if (norm_p > 0.0) {
447:         /*********************************************************************/
448:         /* Follow the direction to the boundary of the trust region.         */
449:         /*********************************************************************/

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

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

456:         /*********************************************************************/
457:         /* Update objective function.                                        */
458:         /*********************************************************************/

460:         cg->o_fcn += step * (0.5 * step * kappa - rz);
461:       } else {
462:         /*********************************************************************/
463:         /* The norm of the direction is zero; there is nothing to follow.    */
464:         /*********************************************************************/
465:       }
466:       break;
467:     }

469:     /*************************************************************************/
470:     /* Now we can update the direction and residual.                         */
471:     /*************************************************************************/

473:     PetscCall(VecAXPY(d, alpha, p));   /* d = d + alpha p   */
474:     PetscCall(VecAXPY(r, -alpha, z));  /* r = r - alpha Q p */
475:     PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r      */

477:     switch (cg->dtype) {
478:     case GLTR_PRECONDITIONED_DIRECTION:
479:       norm_d = norm_dp1;
480:       break;

482:     default:
483:       PetscCall(VecDot(d, d, &norm_d));
484:       break;
485:     }
486:     cg->norm_d = PetscSqrtReal(norm_d);

488:     /*************************************************************************/
489:     /* Update objective function.                                            */
490:     /*************************************************************************/

492:     cg->o_fcn -= 0.5 * alpha * rz;

494:     /*************************************************************************/
495:     /* Check that the preconditioner appears positive semidefinite.          */
496:     /*************************************************************************/

498:     rzm1 = rz;
499:     PetscCall(VecDot(r, z, &rz)); /* rz = r^T z        */
500:     if (rz < 0.0) {
501:       /***********************************************************************/
502:       /* The preconditioner is indefinite.                                   */
503:       /***********************************************************************/

505:       ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
506:       PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: cg indefinite preconditioner: rz=%g\n", (double)rz));
507:       break;
508:     }

510:     /*************************************************************************/
511:     /* As far as we know, the preconditioner is positive semidefinite.       */
512:     /* Compute the residual and check for convergence.                       */
513:     /*************************************************************************/

515:     cg->norm_r[l_size + 1] = PetscSqrtReal(rz); /* norm_r = |r|_M   */

517:     switch (ksp->normtype) {
518:     case KSP_NORM_PRECONDITIONED:
519:       PetscCall(VecNorm(z, NORM_2, &norm_r)); /* norm_r = |z|      */
520:       break;

522:     case KSP_NORM_UNPRECONDITIONED:
523:       PetscCall(VecNorm(r, NORM_2, &norm_r)); /* norm_r = |r|      */
524:       break;

526:     case KSP_NORM_NATURAL:
527:       norm_r = cg->norm_r[l_size + 1]; /* norm_r = |r|_M    */
528:       break;

530:     default:
531:       norm_r = 0.0;
532:       break;
533:     }

535:     PetscCall(KSPLogResidualHistory(ksp, norm_r));
536:     PetscCall(KSPMonitor(ksp, ksp->its, norm_r));
537:     ksp->rnorm = norm_r;

539:     PetscCall((*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP));
540:     if (ksp->reason) {
541:       /***********************************************************************/
542:       /* The method has converged.                                           */
543:       /***********************************************************************/

545:       PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: cg truncated step: rnorm=%g, radius=%g\n", (double)norm_r, (double)cg->radius));
546:       break;
547:     }

549:     /*************************************************************************/
550:     /* We have not converged yet.  Check for breakdown.                      */
551:     /*************************************************************************/

553:     beta = rz / rzm1;
554:     if (PetscAbsReal(beta) <= 0.0) {
555:       /***********************************************************************/
556:       /* Conjugate gradients has broken down.                                */
557:       /***********************************************************************/

559:       ksp->reason = KSP_DIVERGED_BREAKDOWN;
560:       PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: breakdown: beta=%g\n", (double)beta));
561:       break;
562:     }

564:     /*************************************************************************/
565:     /* Check iteration limit.                                                */
566:     /*************************************************************************/

568:     if (ksp->its >= max_cg_its) {
569:       ksp->reason = KSP_DIVERGED_ITS;
570:       PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: iterlim: its=%" PetscInt_FMT "\n", ksp->its));
571:       break;
572:     }

574:     /*************************************************************************/
575:     /* Update p and the norms.                                               */
576:     /*************************************************************************/

578:     cg->beta[l_size] = beta;
579:     PetscCall(VecAYPX(p, beta, z)); /* p = z + beta p    */

581:     switch (cg->dtype) {
582:     case GLTR_PRECONDITIONED_DIRECTION:
583:       dMp    = beta * (dMp + alpha * norm_p);
584:       norm_p = beta * (rzm1 + beta * norm_p);
585:       break;

587:     default:
588:       PetscCall(VecDot(d, p, &dMp));
589:       PetscCall(VecDot(p, p, &norm_p));
590:       break;
591:     }

593:     /*************************************************************************/
594:     /* Compute the new direction and update the iteration.                   */
595:     /*************************************************************************/

597:     PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p         */
598:     PetscCall(VecDot(p, z, &kappa));         /* kappa = p^T Q p   */
599:     ++ksp->its;

601:     /*************************************************************************/
602:     /* Update the Lanczos tridiagonal matrix.                            */
603:     /*************************************************************************/

605:     ++l_size;
606:     cg->offd[t_size] = PetscSqrtReal(beta) / PetscAbsReal(alpha);
607:     cg->diag[t_size] = kappa / rz + beta / alpha;
608:     ++t_size;

610:     /*************************************************************************/
611:     /* Check for breakdown of the conjugate gradient method; this occurs     */
612:     /* when kappa is zero.                                                   */
613:     /*************************************************************************/

615:     if (PetscAbsReal(kappa) <= 0.0) {
616:       /***********************************************************************/
617:       /* The method breaks down; move along the direction as if the matrix   */
618:       /* were indefinite.                                                    */
619:       /***********************************************************************/

621:       ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
622:       PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: cg breakdown: kappa=%g\n", (double)kappa));

624:       if (cg->radius && norm_p > 0.0) {
625:         /*********************************************************************/
626:         /* Follow direction to boundary.                                     */
627:         /*********************************************************************/

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

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

634:         /*********************************************************************/
635:         /* Update objective function.                                        */
636:         /*********************************************************************/

638:         cg->o_fcn += step * (0.5 * step * kappa - rz);
639:       } else if (cg->radius) {
640:         /*********************************************************************/
641:         /* The norm of the direction is zero; there is nothing to follow.    */
642:         /*********************************************************************/
643:       }
644:       break;
645:     }
646:   }

648:   /***************************************************************************/
649:   /* Check to see if we need to continue with the Lanczos method.            */
650:   /***************************************************************************/

652:   if (!cg->radius) {
653:     /*************************************************************************/
654:     /* There is no radius.  Therefore, we cannot move along the boundary.    */
655:     /*************************************************************************/
656:     PetscFunctionReturn(PETSC_SUCCESS);
657:   }

659:   if (KSP_CONVERGED_NEG_CURVE != ksp->reason) {
660:     /*************************************************************************/
661:     /* The method either converged to an interior point, hit the boundary of */
662:     /* the trust-region without encountering a direction of negative         */
663:     /* curvature or the iteration limit was reached.                         */
664:     /*************************************************************************/
665:     PetscFunctionReturn(PETSC_SUCCESS);
666:   }

668:   /***************************************************************************/
669:   /* Switch to constructing the Lanczos basis by way of the conjugate        */
670:   /* directions.                                                             */
671:   /***************************************************************************/

673:   for (i = 0; i < max_lanczos_its; ++i) {
674:     /*************************************************************************/
675:     /* Check for breakdown of the conjugate gradient method; this occurs     */
676:     /* when kappa is zero.                                                   */
677:     /*************************************************************************/

679:     if (PetscAbsReal(kappa) <= 0.0) {
680:       ksp->reason = KSP_DIVERGED_BREAKDOWN;
681:       PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: lanczos breakdown: kappa=%g\n", (double)kappa));
682:       break;
683:     }

685:     /*************************************************************************/
686:     /* Update the direction and residual.                                    */
687:     /*************************************************************************/

689:     alpha             = rz / kappa;
690:     cg->alpha[l_size] = alpha;

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

695:     /*************************************************************************/
696:     /* Check that the preconditioner appears positive semidefinite.          */
697:     /*************************************************************************/

699:     rzm1 = rz;
700:     PetscCall(VecDot(r, z, &rz)); /* rz = r^T z        */
701:     if (rz < 0.0) {
702:       /***********************************************************************/
703:       /* The preconditioner is indefinite.                                   */
704:       /***********************************************************************/

706:       ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
707:       PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: lanczos indefinite preconditioner: rz=%g\n", (double)rz));
708:       break;
709:     }

711:     /*************************************************************************/
712:     /* As far as we know, the preconditioner is positive definite.  Compute  */
713:     /* the residual.  Do NOT check for convergence.                          */
714:     /*************************************************************************/

716:     cg->norm_r[l_size + 1] = PetscSqrtReal(rz); /* norm_r = |r|_M    */

718:     switch (ksp->normtype) {
719:     case KSP_NORM_PRECONDITIONED:
720:       PetscCall(VecNorm(z, NORM_2, &norm_r)); /* norm_r = |z|      */
721:       break;

723:     case KSP_NORM_UNPRECONDITIONED:
724:       PetscCall(VecNorm(r, NORM_2, &norm_r)); /* norm_r = |r|      */
725:       break;

727:     case KSP_NORM_NATURAL:
728:       norm_r = cg->norm_r[l_size + 1]; /* norm_r = |r|_M    */
729:       break;

731:     default:
732:       norm_r = 0.0;
733:       break;
734:     }

736:     PetscCall(KSPLogResidualHistory(ksp, norm_r));
737:     PetscCall(KSPMonitor(ksp, ksp->its, norm_r));
738:     ksp->rnorm = norm_r;

740:     /*************************************************************************/
741:     /* Check for breakdown.                                                  */
742:     /*************************************************************************/

744:     beta = rz / rzm1;
745:     if (PetscAbsReal(beta) <= 0.0) {
746:       /***********************************************************************/
747:       /* Conjugate gradients has broken down.                                */
748:       /***********************************************************************/

750:       ksp->reason = KSP_DIVERGED_BREAKDOWN;
751:       PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: breakdown: beta=%g\n", (double)beta));
752:       break;
753:     }

755:     /*************************************************************************/
756:     /* Update p and the norms.                                               */
757:     /*************************************************************************/

759:     cg->beta[l_size] = beta;
760:     PetscCall(VecAYPX(p, beta, z)); /* p = z + beta p    */

762:     /*************************************************************************/
763:     /* Compute the new direction and update the iteration.                   */
764:     /*************************************************************************/

766:     PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p         */
767:     PetscCall(VecDot(p, z, &kappa));         /* kappa = p^T Q p   */
768:     ++ksp->its;

770:     /*************************************************************************/
771:     /* Update the Lanczos tridiagonal matrix.                                */
772:     /*************************************************************************/

774:     ++l_size;
775:     cg->offd[t_size] = PetscSqrtReal(beta) / PetscAbsReal(alpha);
776:     cg->diag[t_size] = kappa / rz + beta / alpha;
777:     ++t_size;
778:   }

780:   /***************************************************************************/
781:   /* We have the Lanczos basis, solve the tridiagonal trust-region problem   */
782:   /* to obtain the Lanczos direction.  We know that the solution lies on     */
783:   /* the boundary of the trust region.  We start by checking that the        */
784:   /* workspace allocated is large enough.                                    */
785:   /***************************************************************************/
786:   /* Note that the current version only computes the solution by using the   */
787:   /* preconditioned direction.  Need to think about how to do the            */
788:   /* unpreconditioned direction calculation.                                 */
789:   /***************************************************************************/

791:   if (t_size > cg->alloced) {
792:     if (cg->alloced) {
793:       PetscCall(PetscFree(cg->rwork));
794:       PetscCall(PetscFree(cg->iwork));
795:       cg->alloced += cg->init_alloc;
796:     } else {
797:       cg->alloced = cg->init_alloc;
798:     }

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

802:     cg->alloced = PetscMin(cg->alloced, t_size);
803:     PetscCall(PetscMalloc2(10 * cg->alloced, &cg->rwork, 5 * cg->alloced, &cg->iwork));
804:   }

806:   /***************************************************************************/
807:   /* Set up the required vectors.                                            */
808:   /***************************************************************************/

810:   t_soln = cg->rwork + 0 * t_size; /* Solution          */
811:   t_diag = cg->rwork + 1 * t_size; /* Diagonal of T     */
812:   t_offd = cg->rwork + 2 * t_size; /* Off-diagonal of T */
813:   e_valu = cg->rwork + 3 * t_size; /* Eigenvalues of T  */
814:   e_vect = cg->rwork + 4 * t_size; /* Eigenvector of T  */
815:   e_rwrk = cg->rwork + 5 * t_size; /* Eigen workspace   */

817:   e_iblk = cg->iwork + 0 * t_size; /* Eigen blocks      */
818:   e_splt = cg->iwork + 1 * t_size; /* Eigen splits      */
819:   e_iwrk = cg->iwork + 2 * t_size; /* Eigen workspace   */

821:   /***************************************************************************/
822:   /* Compute the minimum eigenvalue of T.                                    */
823:   /***************************************************************************/

825:   vl = 0.0;
826:   vu = 0.0;
827:   il = 1;
828:   iu = 1;

830:   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));

832:   if ((0 != info) || (1 != e_valus)) {
833:     /*************************************************************************/
834:     /* Calculation of the minimum eigenvalue failed.  Return the             */
835:     /* Steihaug-Toint direction.                                             */
836:     /*************************************************************************/

838:     PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute eigenvalue.\n"));
839:     ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
840:     PetscFunctionReturn(PETSC_SUCCESS);
841:   }

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

845:   /***************************************************************************/
846:   /* Compute the initial value of lambda to make (T + lambda I) positive      */
847:   /* definite.                                                               */
848:   /***************************************************************************/

850:   pert = cg->init_pert;
851:   if (e_valu[0] < 0.0) cg->lambda = pert - e_valu[0];

853:   while (1) {
854:     for (i = 0; i < t_size; ++i) {
855:       t_diag[i] = cg->diag[i] + cg->lambda;
856:       t_offd[i] = cg->offd[i];
857:     }

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

861:     if (0 == info) break;

863:     pert += pert;
864:     cg->lambda = cg->lambda * (1.0 + pert) + pert;
865:   }

867:   /***************************************************************************/
868:   /* Compute the initial step and its norm.                                  */
869:   /***************************************************************************/

871:   nrhs = 1;
872:   nldb = t_size;

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

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

879:   if (0 != info) {
880:     /*************************************************************************/
881:     /* Calculation of the initial step failed; return the Steihaug-Toint     */
882:     /* direction.                                                            */
883:     /*************************************************************************/

885:     PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute step.\n"));
886:     ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
887:     PetscFunctionReturn(PETSC_SUCCESS);
888:   }

890:   norm_t = 0.;
891:   for (i = 0; i < t_size; ++i) norm_t += t_soln[i] * t_soln[i];
892:   norm_t = PetscSqrtReal(norm_t);

894:   /***************************************************************************/
895:   /* Determine the case we are in.                                           */
896:   /***************************************************************************/

898:   if (norm_t <= cg->radius) {
899:     /*************************************************************************/
900:     /* The step is within the trust region; check if we are in the hard case */
901:     /* and need to move to the boundary by following a direction of negative */
902:     /* curvature.                                                            */
903:     /*************************************************************************/

905:     if ((e_valu[0] <= 0.0) && (norm_t < cg->radius)) {
906:       /***********************************************************************/
907:       /* This is the hard case; compute the eigenvector associated with the  */
908:       /* minimum eigenvalue and move along this direction to the boundary.   */
909:       /***********************************************************************/

911:       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));

913:       if (0 != info) {
914:         /*********************************************************************/
915:         /* Calculation of the minimum eigenvalue failed.  Return the         */
916:         /* Steihaug-Toint direction.                                         */
917:         /*********************************************************************/

919:         PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute eigenvector.\n"));
920:         ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
921:         PetscFunctionReturn(PETSC_SUCCESS);
922:       }

924:       coef1 = 0.0;
925:       coef2 = 0.0;
926:       coef3 = -cg->radius * cg->radius;
927:       for (i = 0; i < t_size; ++i) {
928:         coef1 += e_vect[i] * e_vect[i];
929:         coef2 += e_vect[i] * t_soln[i];
930:         coef3 += t_soln[i] * t_soln[i];
931:       }

933:       coef3 = PetscSqrtReal(coef2 * coef2 - coef1 * coef3);
934:       root1 = (-coef2 + coef3) / coef1;
935:       root2 = (-coef2 - coef3) / coef1;

937:       /***********************************************************************/
938:       /* Compute objective value for (t_soln + root1 * e_vect)               */
939:       /***********************************************************************/

941:       for (i = 0; i < t_size; ++i) e_rwrk[i] = t_soln[i] + root1 * e_vect[i];

943:       obj1 = e_rwrk[0] * (0.5 * (cg->diag[0] * e_rwrk[0] + cg->offd[1] * e_rwrk[1]) + cg->norm_r[0]);
944:       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]);
945:       obj1 += 0.5 * e_rwrk[i] * (cg->offd[i] * e_rwrk[i - 1] + cg->diag[i] * e_rwrk[i]);

947:       /***********************************************************************/
948:       /* Compute objective value for (t_soln + root2 * e_vect)               */
949:       /***********************************************************************/

951:       for (i = 0; i < t_size; ++i) e_rwrk[i] = t_soln[i] + root2 * e_vect[i];

953:       obj2 = e_rwrk[0] * (0.5 * (cg->diag[0] * e_rwrk[0] + cg->offd[1] * e_rwrk[1]) + cg->norm_r[0]);
954:       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]);
955:       obj2 += 0.5 * e_rwrk[i] * (cg->offd[i] * e_rwrk[i - 1] + cg->diag[i] * e_rwrk[i]);

957:       /***********************************************************************/
958:       /* Choose the point with the best objective function value.            */
959:       /***********************************************************************/

961:       if (obj1 <= obj2) {
962:         for (i = 0; i < t_size; ++i) t_soln[i] += root1 * e_vect[i];
963:       } else {
964:         for (i = 0; i < t_size; ++i) t_soln[i] += root2 * e_vect[i];
965:       }
966:     } else {
967:       /***********************************************************************/
968:       /* The matrix is positive definite or there was no room to move; the   */
969:       /* solution is already contained in t_soln.                            */
970:       /***********************************************************************/
971:     }
972:   } else {
973:     /*************************************************************************/
974:     /* The step is outside the trust-region.  Compute the correct value for  */
975:     /* lambda by performing Newton's method.                                 */
976:     /*************************************************************************/

978:     for (i = 0; i < max_newton_its; ++i) {
979:       /***********************************************************************/
980:       /* Check for convergence.                                              */
981:       /***********************************************************************/

983:       if (PetscAbsReal(norm_t - cg->radius) <= cg->newton_tol * cg->radius) break;

985:       /***********************************************************************/
986:       /* Compute the update.                                                 */
987:       /***********************************************************************/

989:       PetscCall(PetscArraycpy(e_rwrk, t_soln, t_size));

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

993:       if (0 != info) {
994:         /*********************************************************************/
995:         /* Calculation of the step failed; return the Steihaug-Toint         */
996:         /* direction.                                                        */
997:         /*********************************************************************/

999:         PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute step.\n"));
1000:         ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
1001:         PetscFunctionReturn(PETSC_SUCCESS);
1002:       }

1004:       /***********************************************************************/
1005:       /* Modify lambda.                                                      */
1006:       /***********************************************************************/

1008:       norm_w = 0.;
1009:       for (j = 0; j < t_size; ++j) norm_w += t_soln[j] * e_rwrk[j];

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

1013:       /***********************************************************************/
1014:       /* Factor T + lambda I                                                 */
1015:       /***********************************************************************/

1017:       for (j = 0; j < t_size; ++j) {
1018:         t_diag[j] = cg->diag[j] + cg->lambda;
1019:         t_offd[j] = cg->offd[j];
1020:       }

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

1024:       if (0 != info) {
1025:         /*********************************************************************/
1026:         /* Calculation of factorization failed; return the Steihaug-Toint    */
1027:         /* direction.                                                        */
1028:         /*********************************************************************/

1030:         PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: factorization failed.\n"));
1031:         ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
1032:         PetscFunctionReturn(PETSC_SUCCESS);
1033:       }

1035:       /***********************************************************************/
1036:       /* Compute the new step and its norm.                                  */
1037:       /***********************************************************************/

1039:       t_soln[0] = -cg->norm_r[0];
1040:       for (j = 1; j < t_size; ++j) t_soln[j] = 0.0;

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

1044:       if (0 != info) {
1045:         /*********************************************************************/
1046:         /* Calculation of the step failed; return the Steihaug-Toint         */
1047:         /* direction.                                                        */
1048:         /*********************************************************************/

1050:         PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute step.\n"));
1051:         ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
1052:         PetscFunctionReturn(PETSC_SUCCESS);
1053:       }

1055:       norm_t = 0.;
1056:       for (j = 0; j < t_size; ++j) norm_t += t_soln[j] * t_soln[j];
1057:       norm_t = PetscSqrtReal(norm_t);
1058:     }

1060:     /*************************************************************************/
1061:     /* Check for convergence.                                                */
1062:     /*************************************************************************/

1064:     if (PetscAbsReal(norm_t - cg->radius) > cg->newton_tol * cg->radius) {
1065:       /***********************************************************************/
1066:       /* Newton method failed to converge in iteration limit.                */
1067:       /***********************************************************************/

1069:       PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to converge.\n"));
1070:       ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
1071:       PetscFunctionReturn(PETSC_SUCCESS);
1072:     }
1073:   }

1075:   /***************************************************************************/
1076:   /* Recover the norm of the direction and objective function value.         */
1077:   /***************************************************************************/

1079:   cg->norm_d = norm_t;

1081:   cg->o_fcn = t_soln[0] * (0.5 * (cg->diag[0] * t_soln[0] + cg->offd[1] * t_soln[1]) + cg->norm_r[0]);
1082:   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]);
1083:   cg->o_fcn += 0.5 * t_soln[i] * (cg->offd[i] * t_soln[i - 1] + cg->diag[i] * t_soln[i]);

1085:   /***************************************************************************/
1086:   /* Recover the direction.                                                  */
1087:   /***************************************************************************/

1089:   sigma = -1;

1091:   /***************************************************************************/
1092:   /* Start conjugate gradient method from the beginning                      */
1093:   /***************************************************************************/

1095:   PetscCall(VecCopy(ksp->vec_rhs, r)); /* r = -grad         */
1096:   PetscCall(KSP_PCApply(ksp, r, z));   /* z = inv(M) r      */

1098:   /***************************************************************************/
1099:   /* Accumulate Q * s                                                        */
1100:   /***************************************************************************/

1102:   PetscCall(VecCopy(z, d));
1103:   PetscCall(VecScale(d, sigma * t_soln[0] / cg->norm_r[0]));

1105:   /***************************************************************************/
1106:   /* Compute the first direction.                                            */
1107:   /***************************************************************************/

1109:   PetscCall(VecCopy(z, p));                /* p = z             */
1110:   PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p         */
1111:   ++ksp->its;

1113:   for (i = 0; i < l_size - 1; ++i) {
1114:     /*************************************************************************/
1115:     /* Update the residual and direction.                                    */
1116:     /*************************************************************************/

1118:     alpha = cg->alpha[i];
1119:     if (alpha >= 0.0) sigma = -sigma;

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

1124:     /*************************************************************************/
1125:     /* Accumulate Q * s                                                      */
1126:     /*************************************************************************/

1128:     PetscCall(VecAXPY(d, sigma * t_soln[i + 1] / cg->norm_r[i + 1], z));

1130:     /*************************************************************************/
1131:     /* Update p.                                                             */
1132:     /*************************************************************************/

1134:     beta = cg->beta[i];
1135:     PetscCall(VecAYPX(p, beta, z));          /* p = z + beta p    */
1136:     PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p         */
1137:     ++ksp->its;
1138:   }

1140:   /***************************************************************************/
1141:   /* Update the residual and direction.                                      */
1142:   /***************************************************************************/

1144:   alpha = cg->alpha[i];
1145:   if (alpha >= 0.0) sigma = -sigma;

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

1150:   /***************************************************************************/
1151:   /* Accumulate Q * s                                                        */
1152:   /***************************************************************************/

1154:   PetscCall(VecAXPY(d, sigma * t_soln[i + 1] / cg->norm_r[i + 1], z));

1156:   /***************************************************************************/
1157:   /* Set the termination reason.                                             */
1158:   /***************************************************************************/

1160:   ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
1161:   PetscFunctionReturn(PETSC_SUCCESS);
1162: #endif
1163: }

1165: static PetscErrorCode KSPCGSetUp_GLTR(KSP ksp)
1166: {
1167:   KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
1168:   PetscInt    max_its;

1170:   PetscFunctionBegin;
1171:   /***************************************************************************/
1172:   /* Determine the total maximum number of iterations.                       */
1173:   /***************************************************************************/

1175:   max_its = ksp->max_it + cg->max_lanczos_its + 1;

1177:   /***************************************************************************/
1178:   /* Set work vectors needed by conjugate gradient method and allocate       */
1179:   /* workspace for Lanczos matrix.                                           */
1180:   /***************************************************************************/

1182:   PetscCall(KSPSetWorkVecs(ksp, 3));
1183:   if (cg->diag) {
1184:     PetscCall(PetscArrayzero(cg->diag, max_its));
1185:     PetscCall(PetscArrayzero(cg->offd, max_its));
1186:     PetscCall(PetscArrayzero(cg->alpha, max_its));
1187:     PetscCall(PetscArrayzero(cg->beta, max_its));
1188:     PetscCall(PetscArrayzero(cg->norm_r, max_its));
1189:   } else {
1190:     PetscCall(PetscCalloc5(max_its, &cg->diag, max_its, &cg->offd, max_its, &cg->alpha, max_its, &cg->beta, max_its, &cg->norm_r));
1191:   }
1192:   PetscFunctionReturn(PETSC_SUCCESS);
1193: }

1195: static PetscErrorCode KSPCGDestroy_GLTR(KSP ksp)
1196: {
1197:   KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;

1199:   PetscFunctionBegin;
1200:   PetscCall(PetscFree5(cg->diag, cg->offd, cg->alpha, cg->beta, cg->norm_r));
1201:   if (cg->alloced) PetscCall(PetscFree2(cg->rwork, cg->iwork));
1202:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetRadius_C", NULL));
1203:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetNormD_C", NULL));
1204:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", NULL));
1205:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGLTRGetMinEig_C", NULL));
1206:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGLTRGetLambda_C", NULL));
1207:   PetscCall(KSPDestroyDefault(ksp));
1208:   PetscFunctionReturn(PETSC_SUCCESS);
1209: }

1211: static PetscErrorCode KSPCGSetRadius_GLTR(KSP ksp, PetscReal radius)
1212: {
1213:   KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;

1215:   PetscFunctionBegin;
1216:   cg->radius = radius;
1217:   PetscFunctionReturn(PETSC_SUCCESS);
1218: }

1220: static PetscErrorCode KSPCGGetNormD_GLTR(KSP ksp, PetscReal *norm_d)
1221: {
1222:   KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;

1224:   PetscFunctionBegin;
1225:   *norm_d = cg->norm_d;
1226:   PetscFunctionReturn(PETSC_SUCCESS);
1227: }

1229: static PetscErrorCode KSPCGGetObjFcn_GLTR(KSP ksp, PetscReal *o_fcn)
1230: {
1231:   KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;

1233:   PetscFunctionBegin;
1234:   *o_fcn = cg->o_fcn;
1235:   PetscFunctionReturn(PETSC_SUCCESS);
1236: }

1238: static PetscErrorCode KSPGLTRGetMinEig_GLTR(KSP ksp, PetscReal *e_min)
1239: {
1240:   KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;

1242:   PetscFunctionBegin;
1243:   *e_min = cg->e_min;
1244:   PetscFunctionReturn(PETSC_SUCCESS);
1245: }

1247: static PetscErrorCode KSPGLTRGetLambda_GLTR(KSP ksp, PetscReal *lambda)
1248: {
1249:   KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;

1251:   PetscFunctionBegin;
1252:   *lambda = cg->lambda;
1253:   PetscFunctionReturn(PETSC_SUCCESS);
1254: }

1256: static PetscErrorCode KSPCGSetFromOptions_GLTR(KSP ksp, PetscOptionItems *PetscOptionsObject)
1257: {
1258:   KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;

1260:   PetscFunctionBegin;
1261:   PetscOptionsHeadBegin(PetscOptionsObject, "KSP GLTR options");

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

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

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

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

1274:   PetscOptionsHeadEnd();
1275:   PetscFunctionReturn(PETSC_SUCCESS);
1276: }

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

1281:    Options Database Key:
1282: .  -ksp_cg_radius <r> - Trust Region Radius

1284:    Level: developer

1286:    Notes:
1287:    Uses preconditioned conjugate gradient to compute  an approximate minimizer of the quadratic function

1289:    $$
1290:    q(s) = g^T * s + .5 * s^T * H * s
1291:    $$

1293:    subject to the trust region constraint

1295:    $$
1296:    || s || \le delta,
1297:    $$

1299:    where
1300: .vb
1301:      delta is the trust region radius,
1302:      g is the gradient vector,
1303:      H is the Hessian approximation,
1304:      M is the positive definite preconditioner matrix.
1305: .ve

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

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

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

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

1318: PETSC_EXTERN PetscErrorCode KSPCreate_GLTR(KSP ksp)
1319: {
1320:   KSPCG_GLTR *cg;

1322:   PetscFunctionBegin;
1323:   PetscCall(PetscNew(&cg));
1324:   cg->radius = 0.0;
1325:   cg->dtype  = GLTR_UNPRECONDITIONED_DIRECTION;

1327:   cg->init_pert  = 1.0e-8;
1328:   cg->eigen_tol  = 1.0e-10;
1329:   cg->newton_tol = 1.0e-6;

1331:   cg->alloced    = 0;
1332:   cg->init_alloc = 1024;

1334:   cg->max_lanczos_its = 20;
1335:   cg->max_newton_its  = 10;

1337:   ksp->data = (void *)cg;
1338:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 3));
1339:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 2));
1340:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 2));
1341:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
1342:   PetscCall(KSPSetConvergedNegativeCurvature(ksp, PETSC_TRUE));

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

1349:   ksp->ops->setup          = KSPCGSetUp_GLTR;
1350:   ksp->ops->solve          = KSPCGSolve_GLTR;
1351:   ksp->ops->destroy        = KSPCGDestroy_GLTR;
1352:   ksp->ops->setfromoptions = KSPCGSetFromOptions_GLTR;
1353:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
1354:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
1355:   ksp->ops->view           = NULL;

1357:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetRadius_C", KSPCGSetRadius_GLTR));
1358:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetNormD_C", KSPCGGetNormD_GLTR));
1359:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", KSPCGGetObjFcn_GLTR));
1360:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGLTRGetMinEig_C", KSPGLTRGetMinEig_GLTR));
1361:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGLTRGetLambda_C", KSPGLTRGetLambda_GLTR));
1362:   PetscFunctionReturn(PETSC_SUCCESS);
1363: }