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: [](chapter_ksp), `KSPGLTR`, `KSPGLTRGetLambda()`
 25: @*/
 26: PetscErrorCode KSPGLTRGetMinEig(KSP ksp, PetscReal *e_min)
 27: {
 29:   PetscUseMethod(ksp, "KSPGLTRGetMinEig_C", (KSP, PetscReal *), (ksp, e_min));
 30:   return 0;
 31: }

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

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

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

 44:     Level: advanced

 46: .seealso: [](chapter_ksp), `KSPGLTR`, `KSPGLTRGetMinEig()`
 47: @*/
 48: PetscErrorCode KSPGLTRGetLambda(KSP ksp, PetscReal *lambda)
 49: {
 51:   PetscUseMethod(ksp, "KSPGLTRGetLambda_C", (KSP, PetscReal *), (ksp, lambda));
 52:   return 0;
 53: }

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

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

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

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

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

 82:   /***************************************************************************/
 83:   /* Check the arguments and parameters.                                     */
 84:   /***************************************************************************/

 86:   PCGetDiagonalScale(ksp->pc, &diagonalscale);

 90:   /***************************************************************************/
 91:   /* Get the workspace vectors and initialize variables                      */
 92:   /***************************************************************************/

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

101:   PCGetOperators(pc, &Qmat, &Mmat);

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

109:   /***************************************************************************/
110:   /* Initialize objective function direction, and minimum eigenvalue.        */
111:   /***************************************************************************/

113:   cg->o_fcn = 0.0;

115:   VecSet(d, 0.0); /* d = 0             */
116:   cg->norm_d = 0.0;

118:   cg->e_min  = 0.0;
119:   cg->lambda = 0.0;

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

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

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

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

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

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

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

161:       /***********************************************************************/
162:       /* Compute objective function.                                         */
163:       /***********************************************************************/

165:       KSP_MatMult(ksp, Qmat, d, z);
166:       VecAYPX(z, -0.5, ksp->vec_rhs);
167:       VecDot(d, z, &cg->o_fcn);
168:       cg->o_fcn = -cg->o_fcn;
169:       ++ksp->its;
170:     }
171:     return 0;
172:   }

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

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

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

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

196:       /***********************************************************************/
197:       /* Compute objective function.                                         */
198:       /***********************************************************************/

200:       KSP_MatMult(ksp, Qmat, d, z);
201:       VecAYPX(z, -0.5, ksp->vec_rhs);
202:       VecDot(d, z, &cg->o_fcn);
203:       cg->o_fcn = -cg->o_fcn;
204:       ++ksp->its;
205:     }
206:     return 0;
207:   }

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

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

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

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

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

231:   default:
232:     norm_r = 0.0;
233:     break;
234:   }

236:   KSPLogResidualHistory(ksp, norm_r);
237:   KSPMonitor(ksp, ksp->its, norm_r);
238:   ksp->rnorm = norm_r;

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

242:   /***************************************************************************/
243:   /* Compute the first direction and update the iteration.                   */
244:   /***************************************************************************/

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

250:   /***************************************************************************/
251:   /* Check the matrix for numerical problems.                                */
252:   /***************************************************************************/

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

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

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

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

276:       /***********************************************************************/
277:       /* Compute objective function.                                         */
278:       /***********************************************************************/

280:       KSP_MatMult(ksp, Qmat, d, z);
281:       VecAYPX(z, -0.5, ksp->vec_rhs);
282:       VecDot(d, z, &cg->o_fcn);
283:       cg->o_fcn = -cg->o_fcn;
284:       ++ksp->its;
285:     }
286:     return 0;
287:   }

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

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

303:   default:
304:     VecDot(p, p, &norm_p);
305:     break;
306:   }

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

312:   piv = 1.0;

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

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

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

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

334:       step       = PetscSqrtReal(r2 / norm_p);
335:       cg->norm_d = cg->radius;

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

339:       /***********************************************************************/
340:       /* Update objective function.                                          */
341:       /***********************************************************************/

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

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

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

360:       /***********************************************************************/
361:       /* Compute objective function.                                         */
362:       /***********************************************************************/

364:       KSP_MatMult(ksp, Qmat, d, z);
365:       VecAYPX(z, -0.5, ksp->vec_rhs);
366:       VecDot(d, z, &cg->o_fcn);
367:       cg->o_fcn = -cg->o_fcn;
368:       ++ksp->its;
369:     }
370:     return 0;
371:   }

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

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

386:     alpha             = rz / kappa;
387:     cg->alpha[l_size] = alpha;

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

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

403:       ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
404:       PetscInfo(ksp, "KSPCGSolve_GLTR: negative curvature: kappa=%g\n", (double)kappa);

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

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

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

416:         /*********************************************************************/
417:         /* Update objective function.                                        */
418:         /*********************************************************************/

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

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

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

441:       ksp->reason = KSP_CONVERGED_CG_CONSTRAINED;
442:       PetscInfo(ksp, "KSPCGSolve_GLTR: constrained step: radius=%g\n", (double)cg->radius);

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

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

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

454:         /*********************************************************************/
455:         /* Update objective function.                                        */
456:         /*********************************************************************/

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

467:     /*************************************************************************/
468:     /* Now we can update the direction and residual.                         */
469:     /*************************************************************************/

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

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

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

486:     /*************************************************************************/
487:     /* Update objective function.                                            */
488:     /*************************************************************************/

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

492:     /*************************************************************************/
493:     /* Check that the preconditioner appears positive semidefinite.          */
494:     /*************************************************************************/

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

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

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

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

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

520:     case KSP_NORM_UNPRECONDITIONED:
521:       VecNorm(r, NORM_2, &norm_r); /* norm_r = |r|      */
522:       break;

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

528:     default:
529:       norm_r = 0.0;
530:       break;
531:     }

533:     KSPLogResidualHistory(ksp, norm_r);
534:     KSPMonitor(ksp, ksp->its, norm_r);
535:     ksp->rnorm = norm_r;

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

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

547:     /*************************************************************************/
548:     /* We have not converged yet.  Check for breakdown.                      */
549:     /*************************************************************************/

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

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

562:     /*************************************************************************/
563:     /* Check iteration limit.                                                */
564:     /*************************************************************************/

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

572:     /*************************************************************************/
573:     /* Update p and the norms.                                               */
574:     /*************************************************************************/

576:     cg->beta[l_size] = beta;
577:     VecAYPX(p, beta, z); /* p = z + beta p    */

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

585:     default:
586:       VecDot(d, p, &dMp);
587:       VecDot(p, p, &norm_p);
588:       break;
589:     }

591:     /*************************************************************************/
592:     /* Compute the new direction and update the iteration.                   */
593:     /*************************************************************************/

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

599:     /*************************************************************************/
600:     /* Update the Lanczos tridiagonal matrix.                            */
601:     /*************************************************************************/

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

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

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

619:       ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
620:       PetscInfo(ksp, "KSPCGSolve_GLTR: cg breakdown: kappa=%g\n", (double)kappa);

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

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

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

632:         /*********************************************************************/
633:         /* Update objective function.                                        */
634:         /*********************************************************************/

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

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

650:   if (!cg->radius) {
651:     /*************************************************************************/
652:     /* There is no radius.  Therefore, we cannot move along the boundary.    */
653:     /*************************************************************************/
654:     return 0;
655:   }

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

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

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

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

683:     /*************************************************************************/
684:     /* Update the direction and residual.                                    */
685:     /*************************************************************************/

687:     alpha             = rz / kappa;
688:     cg->alpha[l_size] = alpha;

690:     VecAXPY(r, -alpha, z);  /* r = r - alpha Q p */
691:     KSP_PCApply(ksp, r, z); /* z = inv(M) r      */

693:     /*************************************************************************/
694:     /* Check that the preconditioner appears positive semidefinite.          */
695:     /*************************************************************************/

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

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

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

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

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

721:     case KSP_NORM_UNPRECONDITIONED:
722:       VecNorm(r, NORM_2, &norm_r); /* norm_r = |r|      */
723:       break;

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

729:     default:
730:       norm_r = 0.0;
731:       break;
732:     }

734:     KSPLogResidualHistory(ksp, norm_r);
735:     KSPMonitor(ksp, ksp->its, norm_r);
736:     ksp->rnorm = norm_r;

738:     /*************************************************************************/
739:     /* Check for breakdown.                                                  */
740:     /*************************************************************************/

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

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

753:     /*************************************************************************/
754:     /* Update p and the norms.                                               */
755:     /*************************************************************************/

757:     cg->beta[l_size] = beta;
758:     VecAYPX(p, beta, z); /* p = z + beta p    */

760:     /*************************************************************************/
761:     /* Compute the new direction and update the iteration.                   */
762:     /*************************************************************************/

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

768:     /*************************************************************************/
769:     /* Update the Lanczos tridiagonal matrix.                                */
770:     /*************************************************************************/

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

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

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

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

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

804:   /***************************************************************************/
805:   /* Set up the required vectors.                                            */
806:   /***************************************************************************/

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

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

819:   /***************************************************************************/
820:   /* Compute the minimum eigenvalue of T.                                    */
821:   /***************************************************************************/

823:   vl = 0.0;
824:   vu = 0.0;
825:   il = 1;
826:   iu = 1;

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

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

836:     PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute eigenvalue.\n");
837:     ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
838:     return 0;
839:   }

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

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

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

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

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

859:     if (0 == info) break;

861:     pert += pert;
862:     cg->lambda = cg->lambda * (1.0 + pert) + pert;
863:   }

865:   /***************************************************************************/
866:   /* Compute the initial step and its norm.                                  */
867:   /***************************************************************************/

869:   nrhs = 1;
870:   nldb = t_size;

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

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

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

883:     PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute step.\n");
884:     ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
885:     return 0;
886:   }

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

892:   /***************************************************************************/
893:   /* Determine the case we are in.                                           */
894:   /***************************************************************************/

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

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

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

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

917:         PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute eigenvector.\n");
918:         ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
919:         return 0;
920:       }

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

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

935:       /***********************************************************************/
936:       /* Compute objective value for (t_soln + root1 * e_vect)               */
937:       /***********************************************************************/

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

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

945:       /***********************************************************************/
946:       /* Compute objective value for (t_soln + root2 * e_vect)               */
947:       /***********************************************************************/

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

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

955:       /***********************************************************************/
956:       /* Choose the point with the best objective function value.            */
957:       /***********************************************************************/

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

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

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

983:       /***********************************************************************/
984:       /* Compute the update.                                                 */
985:       /***********************************************************************/

987:       PetscArraycpy(e_rwrk, t_soln, t_size);

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

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

997:         PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute step.\n");
998:         ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
999:         return 0;
1000:       }

1002:       /***********************************************************************/
1003:       /* Modify lambda.                                                      */
1004:       /***********************************************************************/

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

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

1011:       /***********************************************************************/
1012:       /* Factor T + lambda I                                                 */
1013:       /***********************************************************************/

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

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

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

1028:         PetscInfo(ksp, "KSPCGSolve_GLTR: factorization failed.\n");
1029:         ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1030:         return 0;
1031:       }

1033:       /***********************************************************************/
1034:       /* Compute the new step and its norm.                                  */
1035:       /***********************************************************************/

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

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

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

1048:         PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute step.\n");
1049:         ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1050:         return 0;
1051:       }

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

1058:     /*************************************************************************/
1059:     /* Check for convergence.                                                */
1060:     /*************************************************************************/

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

1067:       PetscInfo(ksp, "KSPCGSolve_GLTR: failed to converge.\n");
1068:       ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1069:       return 0;
1070:     }
1071:   }

1073:   /***************************************************************************/
1074:   /* Recover the norm of the direction and objective function value.         */
1075:   /***************************************************************************/

1077:   cg->norm_d = norm_t;

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

1083:   /***************************************************************************/
1084:   /* Recover the direction.                                                  */
1085:   /***************************************************************************/

1087:   sigma = -1;

1089:   /***************************************************************************/
1090:   /* Start conjugate gradient method from the beginning                      */
1091:   /***************************************************************************/

1093:   VecCopy(ksp->vec_rhs, r); /* r = -grad         */
1094:   KSP_PCApply(ksp, r, z);   /* z = inv(M) r      */

1096:   /***************************************************************************/
1097:   /* Accumulate Q * s                                                        */
1098:   /***************************************************************************/

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

1103:   /***************************************************************************/
1104:   /* Compute the first direction.                                            */
1105:   /***************************************************************************/

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

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

1116:     alpha = cg->alpha[i];
1117:     if (alpha >= 0.0) sigma = -sigma;

1119:     VecAXPY(r, -alpha, z);  /* r = r - alpha Q p */
1120:     KSP_PCApply(ksp, r, z); /* z = inv(M) r      */

1122:     /*************************************************************************/
1123:     /* Accumulate Q * s                                                      */
1124:     /*************************************************************************/

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

1128:     /*************************************************************************/
1129:     /* Update p.                                                             */
1130:     /*************************************************************************/

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

1138:   /***************************************************************************/
1139:   /* Update the residual and direction.                                      */
1140:   /***************************************************************************/

1142:   alpha = cg->alpha[i];
1143:   if (alpha >= 0.0) sigma = -sigma;

1145:   VecAXPY(r, -alpha, z);  /* r = r - alpha Q p */
1146:   KSP_PCApply(ksp, r, z); /* z = inv(M) r      */

1148:   /***************************************************************************/
1149:   /* Accumulate Q * s                                                        */
1150:   /***************************************************************************/

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

1154:   /***************************************************************************/
1155:   /* Set the termination reason.                                             */
1156:   /***************************************************************************/

1158:   ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1159:   return 0;
1160: #endif
1161: }

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

1168:   /***************************************************************************/
1169:   /* Determine the total maximum number of iterations.                       */
1170:   /***************************************************************************/

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

1174:   /***************************************************************************/
1175:   /* Set work vectors needed by conjugate gradient method and allocate       */
1176:   /* workspace for Lanczos matrix.                                           */
1177:   /***************************************************************************/

1179:   KSPSetWorkVecs(ksp, 3);
1180:   if (cg->diag) {
1181:     PetscArrayzero(cg->diag, max_its);
1182:     PetscArrayzero(cg->offd, max_its);
1183:     PetscArrayzero(cg->alpha, max_its);
1184:     PetscArrayzero(cg->beta, max_its);
1185:     PetscArrayzero(cg->norm_r, max_its);
1186:   } else {
1187:     PetscCalloc5(max_its, &cg->diag, max_its, &cg->offd, max_its, &cg->alpha, max_its, &cg->beta, max_its, &cg->norm_r);
1188:   }
1189:   return 0;
1190: }

1192: static PetscErrorCode KSPCGDestroy_GLTR(KSP ksp)
1193: {
1194:   KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;

1196:   /***************************************************************************/
1197:   /* Free memory allocated for the data.                                     */
1198:   /***************************************************************************/

1200:   PetscFree5(cg->diag, cg->offd, cg->alpha, cg->beta, cg->norm_r);
1201:   if (cg->alloced) PetscFree2(cg->rwork, cg->iwork);

1203:   /***************************************************************************/
1204:   /* Clear composed functions                                                */
1205:   /***************************************************************************/

1207:   PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetRadius_C", NULL);
1208:   PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetNormD_C", NULL);
1209:   PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", NULL);
1210:   PetscObjectComposeFunction((PetscObject)ksp, "KSPGLTRGetMinEig_C", NULL);
1211:   PetscObjectComposeFunction((PetscObject)ksp, "KSPGLTRGetLambda_C", NULL);

1213:   /***************************************************************************/
1214:   /* Destroy KSP object.                                                     */
1215:   /***************************************************************************/
1216:   KSPDestroyDefault(ksp);
1217:   return 0;
1218: }

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

1224:   cg->radius = radius;
1225:   return 0;
1226: }

1228: static PetscErrorCode KSPCGGetNormD_GLTR(KSP ksp, PetscReal *norm_d)
1229: {
1230:   KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;

1232:   *norm_d = cg->norm_d;
1233:   return 0;
1234: }

1236: static PetscErrorCode KSPCGGetObjFcn_GLTR(KSP ksp, PetscReal *o_fcn)
1237: {
1238:   KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;

1240:   *o_fcn = cg->o_fcn;
1241:   return 0;
1242: }

1244: static PetscErrorCode KSPGLTRGetMinEig_GLTR(KSP ksp, PetscReal *e_min)
1245: {
1246:   KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;

1248:   *e_min = cg->e_min;
1249:   return 0;
1250: }

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

1256:   *lambda = cg->lambda;
1257:   return 0;
1258: }

1260: static PetscErrorCode KSPCGSetFromOptions_GLTR(KSP ksp, PetscOptionItems *PetscOptionsObject)
1261: {
1262:   KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;

1264:   PetscOptionsHeadBegin(PetscOptionsObject, "KSP GLTR options");

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

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

1270:   PetscOptionsReal("-ksp_cg_gltr_init_pert", "Initial perturbation", "", cg->init_pert, &cg->init_pert, NULL);
1271:   PetscOptionsReal("-ksp_cg_gltr_eigen_tol", "Eigenvalue tolerance", "", cg->eigen_tol, &cg->eigen_tol, NULL);
1272:   PetscOptionsReal("-ksp_cg_gltr_newton_tol", "Newton tolerance", "", cg->newton_tol, &cg->newton_tol, NULL);

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

1277:   PetscOptionsHeadEnd();
1278:   return 0;
1279: }

1281: /*MC
1282:      KSPGLTR -   Code to run conjugate gradient method subject to a constraint
1283:          on the solution norm.

1285:    Options Database Keys:
1286: .      -ksp_cg_radius <r> - Trust Region Radius

1288:    Level: developer

1290:   Notes:
1291:   Uses preconditioned conjugate gradient to compute
1292:   an approximate minimizer of the quadratic function

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

1296:    subject to the trust region constraint

1298:             || s || <= delta,

1300:    where

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

1307:    `KSPConvergedReason` may have the additional values
1308: .vb
1309:    KSP_CONVERGED_CG_NEG_CURVE if convergence is reached along a negative curvature direction,
1310:    KSP_CONVERGED_CG_CONSTRAINED if convergence is reached along a constrained step.
1311: .ve

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

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

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

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

1324: PETSC_EXTERN PetscErrorCode KSPCreate_GLTR(KSP ksp)
1325: {
1326:   KSPCG_GLTR *cg;

1328:   PetscNew(&cg);
1329:   cg->radius = 0.0;
1330:   cg->dtype  = GLTR_UNPRECONDITIONED_DIRECTION;

1332:   cg->init_pert  = 1.0e-8;
1333:   cg->eigen_tol  = 1.0e-10;
1334:   cg->newton_tol = 1.0e-6;

1336:   cg->alloced    = 0;
1337:   cg->init_alloc = 1024;

1339:   cg->max_lanczos_its = 20;
1340:   cg->max_newton_its  = 10;

1342:   ksp->data = (void *)cg;
1343:   KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 3);
1344:   KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 2);
1345:   KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 2);
1346:   KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1);

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

1353:   ksp->ops->setup          = KSPCGSetUp_GLTR;
1354:   ksp->ops->solve          = KSPCGSolve_GLTR;
1355:   ksp->ops->destroy        = KSPCGDestroy_GLTR;
1356:   ksp->ops->setfromoptions = KSPCGSetFromOptions_GLTR;
1357:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
1358:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
1359:   ksp->ops->view           = NULL;

1361:   PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetRadius_C", KSPCGSetRadius_GLTR);
1362:   PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetNormD_C", KSPCGGetNormD_GLTR);
1363:   PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", KSPCGGetObjFcn_GLTR);
1364:   PetscObjectComposeFunction((PetscObject)ksp, "KSPGLTRGetMinEig_C", KSPGLTRGetMinEig_GLTR);
1365:   PetscObjectComposeFunction((PetscObject)ksp, "KSPGLTRGetLambda_C", KSPGLTRGetLambda_GLTR);
1366:   return 0;
1367: }