Actual source code: gltr.c


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

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

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

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

 14:     Collective

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

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

 22:     Level: advanced

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

 34: /*@
 35:     KSPGLTRGetLambda - Get the multiplier on the trust-region constraint when using `KSPGLTR`
 36: t
 37:     Not Collective

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

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

 45:     Level: advanced

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

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

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

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

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

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

 84:   PetscFunctionBegin;
 85:   /***************************************************************************/
 86:   /* Check the arguments and parameters.                                     */
 87:   /***************************************************************************/

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

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

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

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

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

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

116:   cg->o_fcn = 0.0;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

315:   piv = 1.0;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

862:     if (0 == info) break;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1080:   cg->norm_d = norm_t;

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

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

1090:   sigma = -1;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1200:   PetscFunctionBegin;
1201:   /***************************************************************************/
1202:   /* Free memory allocated for the data.                                     */
1203:   /***************************************************************************/

1205:   PetscCall(PetscFree5(cg->diag, cg->offd, cg->alpha, cg->beta, cg->norm_r));
1206:   if (cg->alloced) PetscCall(PetscFree2(cg->rwork, cg->iwork));

1208:   /***************************************************************************/
1209:   /* Clear composed functions                                                */
1210:   /***************************************************************************/

1212:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetRadius_C", NULL));
1213:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetNormD_C", NULL));
1214:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", NULL));
1215:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGLTRGetMinEig_C", NULL));
1216:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGLTRGetLambda_C", NULL));

1218:   /***************************************************************************/
1219:   /* Destroy KSP object.                                                     */
1220:   /***************************************************************************/
1221:   PetscCall(KSPDestroyDefault(ksp));
1222:   PetscFunctionReturn(PETSC_SUCCESS);
1223: }

1225: static PetscErrorCode KSPCGSetRadius_GLTR(KSP ksp, PetscReal radius)
1226: {
1227:   KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;

1229:   PetscFunctionBegin;
1230:   cg->radius = radius;
1231:   PetscFunctionReturn(PETSC_SUCCESS);
1232: }

1234: static PetscErrorCode KSPCGGetNormD_GLTR(KSP ksp, PetscReal *norm_d)
1235: {
1236:   KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;

1238:   PetscFunctionBegin;
1239:   *norm_d = cg->norm_d;
1240:   PetscFunctionReturn(PETSC_SUCCESS);
1241: }

1243: static PetscErrorCode KSPCGGetObjFcn_GLTR(KSP ksp, PetscReal *o_fcn)
1244: {
1245:   KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;

1247:   PetscFunctionBegin;
1248:   *o_fcn = cg->o_fcn;
1249:   PetscFunctionReturn(PETSC_SUCCESS);
1250: }

1252: static PetscErrorCode KSPGLTRGetMinEig_GLTR(KSP ksp, PetscReal *e_min)
1253: {
1254:   KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;

1256:   PetscFunctionBegin;
1257:   *e_min = cg->e_min;
1258:   PetscFunctionReturn(PETSC_SUCCESS);
1259: }

1261: static PetscErrorCode KSPGLTRGetLambda_GLTR(KSP ksp, PetscReal *lambda)
1262: {
1263:   KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;

1265:   PetscFunctionBegin;
1266:   *lambda = cg->lambda;
1267:   PetscFunctionReturn(PETSC_SUCCESS);
1268: }

1270: static PetscErrorCode KSPCGSetFromOptions_GLTR(KSP ksp, PetscOptionItems *PetscOptionsObject)
1271: {
1272:   KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;

1274:   PetscFunctionBegin;
1275:   PetscOptionsHeadBegin(PetscOptionsObject, "KSP GLTR options");

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

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

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

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

1288:   PetscOptionsHeadEnd();
1289:   PetscFunctionReturn(PETSC_SUCCESS);
1290: }

1292: /*MC
1293:      KSPGLTR -   Code to run conjugate gradient method subject to a constraint
1294:          on the solution norm.

1296:    Options Database Keys:
1297: .      -ksp_cg_radius <r> - Trust Region Radius

1299:    Level: developer

1301:   Notes:
1302:   Uses preconditioned conjugate gradient to compute
1303:   an approximate minimizer of the quadratic function

1305:             q(s) = g^T * s + .5 * s^T * H * s

1307:    subject to the trust region constraint

1309:             || s || <= delta,

1311:    where

1313:      delta is the trust region radius,
1314:      g is the gradient vector,
1315:      H is the Hessian approximation,
1316:      M is the positive definite preconditioner matrix.

1318:    `KSPConvergedReason` may have the additional values
1319: .vb
1320:    KSP_CONVERGED_NEG_CURVE if convergence is reached along a negative curvature direction,
1321:    KSP_CONVERGED_STEP_LENGTH if convergence is reached along a constrained step.
1322: .ve

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

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

1328:   Reference:
1329: . * -  Gould, N. and Lucidi, S. and Roma, M. and Toint, P., Solving the Trust-Region Subproblem using the Lanczos Method,
1330:    SIAM Journal on Optimization, volume 9, number 2, 1999, 504-525

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

1335: PETSC_EXTERN PetscErrorCode KSPCreate_GLTR(KSP ksp)
1336: {
1337:   KSPCG_GLTR *cg;

1339:   PetscFunctionBegin;
1340:   PetscCall(PetscNew(&cg));
1341:   cg->radius = 0.0;
1342:   cg->dtype  = GLTR_UNPRECONDITIONED_DIRECTION;

1344:   cg->init_pert  = 1.0e-8;
1345:   cg->eigen_tol  = 1.0e-10;
1346:   cg->newton_tol = 1.0e-6;

1348:   cg->alloced    = 0;
1349:   cg->init_alloc = 1024;

1351:   cg->max_lanczos_its = 20;
1352:   cg->max_newton_its  = 10;

1354:   ksp->data = (void *)cg;
1355:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 3));
1356:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 2));
1357:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 2));
1358:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
1359:   PetscCall(KSPSetConvergedNegativeCurvature(ksp, PETSC_TRUE));

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

1366:   ksp->ops->setup          = KSPCGSetUp_GLTR;
1367:   ksp->ops->solve          = KSPCGSolve_GLTR;
1368:   ksp->ops->destroy        = KSPCGDestroy_GLTR;
1369:   ksp->ops->setfromoptions = KSPCGSetFromOptions_GLTR;
1370:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
1371:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
1372:   ksp->ops->view           = NULL;

1374:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetRadius_C", KSPCGSetRadius_GLTR));
1375:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetNormD_C", KSPCGGetNormD_GLTR));
1376:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", KSPCGGetObjFcn_GLTR));
1377:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGLTRGetMinEig_C", KSPGLTRGetMinEig_GLTR));
1378:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGLTRGetLambda_C", KSPGLTRGetLambda_GLTR));
1379:   PetscFunctionReturn(PETSC_SUCCESS);
1380: }