Actual source code: gltr.c

petsc-3.5.2 2014-09-08
Report Typos and Errors
  2: #include <petsc-private/kspimpl.h>             /*I "petscksp.h" I*/
  3: #include <../src/ksp/ksp/impls/cg/gltr/gltrimpl.h>
  4: #include <petscblaslapack.h>

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

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

 14: /*@
 15:     KSPGLTRSetRadius - Sets the radius of the trust region.

 17:     Logically Collective on KSP

 19:     Input Parameters:
 20: +   ksp    - the iterative context
 21: -   radius - the trust region radius (Infinity is the default)

 23:     Options Database Key:
 24: .   -ksp_gltr_radius <r>

 26:     Level: advanced

 28: .keywords: KSP, GLTR, set, trust region radius
 29: @*/
 30: PetscErrorCode  KSPGLTRSetRadius(KSP ksp, PetscReal radius)
 31: {

 36:   if (radius < 0.0) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE, "Radius negative");
 38:   PetscTryMethod(ksp,"KSPGLTRSetRadius_C",(KSP,PetscReal),(ksp,radius));
 39:   return(0);
 40: }

 44: /*@
 45:     KSPGLTRGetNormD - Get norm of the direction.

 47:     Collective

 49:     Input Parameters:
 50: +   ksp    - the iterative context
 51: -   norm_d - the norm of the direction

 53:     Level: advanced

 55: .keywords: KSP, GLTR, get, norm direction
 56: @*/
 57: PetscErrorCode  KSPGLTRGetNormD(KSP ksp, PetscReal *norm_d)
 58: {

 63:   PetscUseMethod(ksp,"KSPGLTRGetNormD_C",(KSP,PetscReal*),(ksp,norm_d));
 64:   return(0);
 65: }

 69: /*@
 70:     KSPGLTRGetObjFcn - Get objective function value.

 72:     Note Collective

 74:     Input Parameters:
 75: +   ksp   - the iterative context
 76: -   o_fcn - the objective function value

 78:     Level: advanced

 80: .keywords: KSP, GLTR, get, objective function
 81: @*/
 82: PetscErrorCode  KSPGLTRGetObjFcn(KSP ksp, PetscReal *o_fcn)
 83: {

 88:   PetscUseMethod(ksp,"KSPGLTRGetObjFcn_C",(KSP,PetscReal*),(ksp,o_fcn));
 89:   return(0);
 90: }

 94: /*@
 95:     KSPGLTRGetMinEig - Get minimum eigenvalue.

 97:     Collective on KSP

 99:     Input Parameters:
100: +   ksp   - the iterative context
101: -   e_min - the minimum eigenvalue

103:     Level: advanced

105: .keywords: KSP, GLTR, get, minimum eigenvalue
106: @*/
107: PetscErrorCode  KSPGLTRGetMinEig(KSP ksp, PetscReal *e_min)
108: {

113:   PetscUseMethod(ksp,"KSPGLTRGetMinEig_C",(KSP,PetscReal*),(ksp,e_min));
114:   return(0);
115: }

119: /*@
120:     KSPGLTRGetLambda - Get multiplier on trust-region constraint.

122:     Not Collective

124:     Input Parameters:
125: +   ksp    - the iterative context
126: -   lambda - the multiplier

128:     Level: advanced

130: .keywords: KSP, GLTR, get, multiplier
131: @*/
132: PetscErrorCode  KSPGLTRGetLambda(KSP ksp, PetscReal *lambda)
133: {

138:   PetscUseMethod(ksp,"KSPGLTRGetLambda_C",(KSP,PetscReal*),(ksp,lambda));
139:   return(0);
140: }

144: PetscErrorCode KSPSolve_GLTR(KSP ksp)
145: {
146: #if defined(PETSC_USE_COMPLEX)
147:   SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP, "GLTR is not available for complex systems");
148: #else
149:   KSP_GLTR     *cg = (KSP_GLTR*)ksp->data;
150:   PetscReal    *t_soln, *t_diag, *t_offd, *e_valu, *e_vect, *e_rwrk;
151:   PetscBLASInt *e_iblk, *e_splt, *e_iwrk;

154:   Mat            Qmat, Mmat;
155:   Vec            r, z, p, d;
156:   PC             pc;

158:   PetscReal norm_r, norm_d, norm_dp1, norm_p, dMp;
159:   PetscReal alpha, beta, kappa, rz, rzm1;
160:   PetscReal rr, r2, piv, step;
161:   PetscReal vl, vu;
162:   PetscReal coef1, coef2, coef3, root1, root2, obj1, obj2;
163:   PetscReal norm_t, norm_w, pert;

165:   PetscInt     i, j, max_cg_its, max_lanczos_its, max_newton_its, sigma;
166:   PetscBLASInt t_size = 0, l_size = 0, il, iu, info;
167:   PetscBLASInt nrhs, nldb;

169: #if !defined(PETSC_MISSING_LAPACK_STEBZ)
170:   PetscBLASInt e_valus, e_splts;
171: #endif
172:   PetscBool diagonalscale;

175:   /***************************************************************************/
176:   /* Check the arguments and parameters.                                     */
177:   /***************************************************************************/

179:   PCGetDiagonalScale(ksp->pc, &diagonalscale);
180:   if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
181:   if (cg->radius < 0.0) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE, "Input error: radius < 0");

183:   /***************************************************************************/
184:   /* Get the workspace vectors and initialize variables                      */
185:   /***************************************************************************/

187:   r2 = cg->radius * cg->radius;
188:   r  = ksp->work[0];
189:   z  = ksp->work[1];
190:   p  = ksp->work[2];
191:   d  = ksp->vec_sol;
192:   pc = ksp->pc;

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

196:   VecGetSize(d, &max_cg_its);
197:   max_cg_its      = PetscMin(max_cg_its, ksp->max_it);
198:   max_lanczos_its = cg->max_lanczos_its;
199:   max_newton_its  = cg->max_newton_its;
200:   ksp->its        = 0;

202:   /***************************************************************************/
203:   /* Initialize objective function direction, and minimum eigenvalue.        */
204:   /***************************************************************************/

206:   cg->o_fcn = 0.0;

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

211:   cg->e_min  = 0.0;
212:   cg->lambda = 0.0;

214:   /***************************************************************************/
215:   /* The first phase of GLTR performs a standard conjugate gradient method,  */
216:   /* but stores the values required for the Lanczos matrix.  We switch to    */
217:   /* the Lanczos when the conjugate gradient method breaks down.  Check the  */
218:   /* right-hand side for numerical problems.  The check for not-a-number and */
219:   /* infinite values need be performed only once.                            */
220:   /***************************************************************************/

222:   VecCopy(ksp->vec_rhs, r);        /* r = -grad         */
223:   VecDot(r, r, &rr);               /* rr = r^T r        */
224:   if (PetscIsInfOrNanScalar(rr)) {
225:     /*************************************************************************/
226:     /* The right-hand side contains not-a-number or an infinite value.       */
227:     /* The gradient step does not work; return a zero value for the step.    */
228:     /*************************************************************************/

230:     ksp->reason = KSP_DIVERGED_NANORINF;
231:     PetscInfo1(ksp, "KSPSolve_GLTR: bad right-hand side: rr=%g\n", (double)rr);
232:     return(0);
233:   }

235:   /***************************************************************************/
236:   /* Check the preconditioner for numerical problems and for positive        */
237:   /* definiteness.  The check for not-a-number and infinite values need be   */
238:   /* performed only once.                                                    */
239:   /***************************************************************************/

241:   KSP_PCApply(ksp, r, z);          /* z = inv(M) r      */
242:   VecDot(r, z, &rz);               /* rz = r^T inv(M) r */
243:   if (PetscIsInfOrNanScalar(rz)) {
244:     /*************************************************************************/
245:     /* The preconditioner contains not-a-number or an infinite value.        */
246:     /* Return the gradient direction intersected with the trust region.      */
247:     /*************************************************************************/

249:     ksp->reason = KSP_DIVERGED_NANORINF;
250:     PetscInfo1(ksp, "KSPSolve_GLTR: bad preconditioner: rz=%g\n", (double)rz);

252:     if (cg->radius) {
253:       if (r2 >= rr) {
254:         alpha      = 1.0;
255:         cg->norm_d = PetscSqrtReal(rr);
256:       } else {
257:         alpha      = PetscSqrtReal(r2 / rr);
258:         cg->norm_d = cg->radius;
259:       }

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

263:       /***********************************************************************/
264:       /* Compute objective function.                                         */
265:       /***********************************************************************/

267:       KSP_MatMult(ksp, Qmat, d, z);
268:       VecAYPX(z, -0.5, ksp->vec_rhs);
269:       VecDot(d, z, &cg->o_fcn);
270:       cg->o_fcn = -cg->o_fcn;
271:       ++ksp->its;
272:     }
273:     return(0);
274:   }

276:   if (rz < 0.0) {
277:     /*************************************************************************/
278:     /* The preconditioner is indefinite.  Because this is the first          */
279:     /* and we do not have a direction yet, we use the gradient step.  Note   */
280:     /* that we cannot use the preconditioned norm when computing the step    */
281:     /* because the matrix is indefinite.                                     */
282:     /*************************************************************************/

284:     ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
285:     PetscInfo1(ksp, "KSPSolve_GLTR: indefinite preconditioner: rz=%g\n", (double)rz);

287:     if (cg->radius) {
288:       if (r2 >= rr) {
289:         alpha      = 1.0;
290:         cg->norm_d = PetscSqrtReal(rr);
291:       } else {
292:         alpha      = PetscSqrtReal(r2 / rr);
293:         cg->norm_d = cg->radius;
294:       }

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

298:       /***********************************************************************/
299:       /* Compute objective function.                                         */
300:       /***********************************************************************/

302:       KSP_MatMult(ksp, Qmat, d, z);
303:       VecAYPX(z, -0.5, ksp->vec_rhs);
304:       VecDot(d, z, &cg->o_fcn);
305:       cg->o_fcn = -cg->o_fcn;
306:       ++ksp->its;
307:     }
308:     return(0);
309:   }

311:   /***************************************************************************/
312:   /* As far as we know, the preconditioner is positive semidefinite.         */
313:   /* Compute and log the residual.  Check convergence because this           */
314:   /* initializes things, but do not terminate until at least one conjugate   */
315:   /* gradient iteration has been performed.                                  */
316:   /***************************************************************************/

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

320:   switch (ksp->normtype) {
321:   case KSP_NORM_PRECONDITIONED:
322:     VecNorm(z, NORM_2, &norm_r);   /* norm_r = |z|      */
323:     break;

325:   case KSP_NORM_UNPRECONDITIONED:
326:     norm_r = PetscSqrtReal(rr);                                 /* norm_r = |r|      */
327:     break;

329:   case KSP_NORM_NATURAL:
330:     norm_r = cg->norm_r[0];                             /* norm_r = |r|_M    */
331:     break;

333:   default:
334:     norm_r = 0.0;
335:     break;
336:   }

338:   KSPLogResidualHistory(ksp, norm_r);
339:   KSPMonitor(ksp, ksp->its, norm_r);
340:   ksp->rnorm = norm_r;

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

344:   /***************************************************************************/
345:   /* Compute the first direction and update the iteration.                   */
346:   /***************************************************************************/

348:   VecCopy(z, p);                   /* p = z             */
349:   KSP_MatMult(ksp, Qmat, p, z);    /* z = Q * p         */
350:   ++ksp->its;

352:   /***************************************************************************/
353:   /* Check the matrix for numerical problems.                                */
354:   /***************************************************************************/

356:   VecDot(p, z, &kappa);            /* kappa = p^T Q p   */
357:   if (PetscIsInfOrNanScalar(kappa)) {
358:     /*************************************************************************/
359:     /* The matrix produced not-a-number or an infinite value.  In this case, */
360:     /* we must stop and use the gradient direction.  This condition need     */
361:     /* only be checked once.                                                 */
362:     /*************************************************************************/

364:     ksp->reason = KSP_DIVERGED_NANORINF;
365:     PetscInfo1(ksp, "KSPSolve_GLTR: bad matrix: kappa=%g\n", (double)kappa);

367:     if (cg->radius) {
368:       if (r2 >= rr) {
369:         alpha      = 1.0;
370:         cg->norm_d = PetscSqrtReal(rr);
371:       } else {
372:         alpha      = PetscSqrtReal(r2 / rr);
373:         cg->norm_d = cg->radius;
374:       }

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

378:       /***********************************************************************/
379:       /* Compute objective function.                                         */
380:       /***********************************************************************/

382:       KSP_MatMult(ksp, Qmat, d, z);
383:       VecAYPX(z, -0.5, ksp->vec_rhs);
384:       VecDot(d, z, &cg->o_fcn);
385:       cg->o_fcn = -cg->o_fcn;
386:       ++ksp->its;
387:     }
388:     return(0);
389:   }

391:   /***************************************************************************/
392:   /* Initialize variables for calculating the norm of the direction and for  */
393:   /* the Lanczos tridiagonal matrix.  Note that we track the diagonal value  */
394:   /* of the Cholesky factorization of the Lanczos matrix in order to         */
395:   /* determine when negative curvature is encountered.                       */
396:   /***************************************************************************/

398:   dMp    = 0.0;
399:   norm_d = 0.0;
400:   switch (cg->dtype) {
401:   case GLTR_PRECONDITIONED_DIRECTION:
402:     norm_p = rz;
403:     break;

405:   default:
406:     VecDot(p, p, &norm_p);
407:     break;
408:   }

410:   cg->diag[t_size] = kappa / rz;
411:   cg->offd[t_size] = 0.0;
412:   ++t_size;

414:   piv = 1.0;

416:   /***************************************************************************/
417:   /* Check for breakdown of the conjugate gradient method; this occurs when  */
418:   /* kappa is zero.                                                          */
419:   /***************************************************************************/

421:   if (PetscAbsReal(kappa) <= 0.0) {
422:     /*************************************************************************/
423:     /* The curvature is zero.  In this case, we must stop and use follow     */
424:     /* the direction of negative curvature since the Lanczos matrix is zero. */
425:     /*************************************************************************/

427:     ksp->reason = KSP_DIVERGED_BREAKDOWN;
428:     PetscInfo1(ksp, "KSPSolve_GLTR: breakdown: kappa=%g\n", (double)kappa);

430:     if (cg->radius && norm_p > 0.0) {
431:       /***********************************************************************/
432:       /* Follow direction of negative curvature to the boundary of the       */
433:       /* trust region.                                                       */
434:       /***********************************************************************/

436:       step       = PetscSqrtReal(r2 / norm_p);
437:       cg->norm_d = cg->radius;

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

441:       /***********************************************************************/
442:       /* Update objective function.                                          */
443:       /***********************************************************************/

445:       cg->o_fcn += step * (0.5 * step * kappa - rz);
446:     } else if (cg->radius) {
447:       /***********************************************************************/
448:       /* The norm of the preconditioned direction is zero; use the gradient  */
449:       /* step.                                                               */
450:       /***********************************************************************/

452:       if (r2 >= rr) {
453:         alpha      = 1.0;
454:         cg->norm_d = PetscSqrtReal(rr);
455:       } else {
456:         alpha      = PetscSqrtReal(r2 / rr);
457:         cg->norm_d = cg->radius;
458:       }

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

462:       /***********************************************************************/
463:       /* Compute objective function.                                         */
464:       /***********************************************************************/

466:       KSP_MatMult(ksp, Qmat, d, z);
467:       VecAYPX(z, -0.5, ksp->vec_rhs);
468:       VecDot(d, z, &cg->o_fcn);
469:       cg->o_fcn = -cg->o_fcn;
470:       ++ksp->its;
471:     }
472:     return(0);
473:   }

475:   /***************************************************************************/
476:   /* Begin the first part of the GLTR algorithm which runs the conjugate     */
477:   /* gradient method until either the problem is solved, we encounter the    */
478:   /* boundary of the trust region, or the conjugate gradient method breaks   */
479:   /* down.                                                                   */
480:   /***************************************************************************/

482:   while (1) {
483:     /*************************************************************************/
484:     /* Know that kappa is nonzero, because we have not broken down, so we    */
485:     /* can compute the steplength.                                           */
486:     /*************************************************************************/

488:     alpha             = rz / kappa;
489:     cg->alpha[l_size] = alpha;

491:     /*************************************************************************/
492:     /* Compute the diagonal value of the Cholesky factorization of the       */
493:     /* Lanczos matrix and check to see if the Lanczos matrix is indefinite.  */
494:     /* This indicates a direction of negative curvature.                     */
495:     /*************************************************************************/

497:     piv = cg->diag[l_size] - cg->offd[l_size]*cg->offd[l_size] / piv;
498:     if (piv <= 0.0) {
499:       /***********************************************************************/
500:       /* In this case, the matrix is indefinite and we have encountered      */
501:       /* a direction of negative curvature.  Follow the direction to the     */
502:       /* boundary of the trust region.                                       */
503:       /***********************************************************************/

505:       ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
506:       PetscInfo1(ksp, "KSPSolve_GLTR: negative curvature: kappa=%g\n", (double)kappa);

508:       if (cg->radius && norm_p > 0.0) {
509:         /*********************************************************************/
510:         /* Follow direction of negative curvature to boundary.               */
511:         /*********************************************************************/

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

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

518:         /*********************************************************************/
519:         /* Update objective function.                                        */
520:         /*********************************************************************/

522:         cg->o_fcn += step * (0.5 * step * kappa - rz);
523:       } else if (cg->radius) {
524:         /*********************************************************************/
525:         /* The norm of the direction is zero; there is nothing to follow.    */
526:         /*********************************************************************/
527:       }
528:       break;
529:     }

531:     /*************************************************************************/
532:     /* Compute the steplength and check for intersection with the trust      */
533:     /* region.                                                               */
534:     /*************************************************************************/

536:     norm_dp1 = norm_d + alpha*(2.0*dMp + alpha*norm_p);
537:     if (cg->radius && norm_dp1 >= r2) {
538:       /***********************************************************************/
539:       /* In this case, the matrix is positive definite as far as we know.    */
540:       /* However, the full step goes beyond the trust region.                */
541:       /***********************************************************************/

543:       ksp->reason = KSP_CONVERGED_CG_CONSTRAINED;
544:       PetscInfo1(ksp, "KSPSolve_GLTR: constrained step: radius=%g\n", (double)cg->radius);

546:       if (norm_p > 0.0) {
547:         /*********************************************************************/
548:         /* Follow the direction to the boundary of the trust region.         */
549:         /*********************************************************************/

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

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

556:         /*********************************************************************/
557:         /* Update objective function.                                        */
558:         /*********************************************************************/

560:         cg->o_fcn += step * (0.5 * step * kappa - rz);
561:       } else {
562:         /*********************************************************************/
563:         /* The norm of the direction is zero; there is nothing to follow.    */
564:         /*********************************************************************/
565:       }
566:       break;
567:     }

569:     /*************************************************************************/
570:     /* Now we can update the direction and residual.                         */
571:     /*************************************************************************/

573:     VecAXPY(d, alpha, p);          /* d = d + alpha p   */
574:     VecAXPY(r, -alpha, z);         /* r = r - alpha Q p */
575:     KSP_PCApply(ksp, r, z);        /* z = inv(M) r      */

577:     switch (cg->dtype) {
578:     case GLTR_PRECONDITIONED_DIRECTION:
579:       norm_d = norm_dp1;
580:       break;

582:     default:
583:       VecDot(d, d, &norm_d);
584:       break;
585:     }
586:     cg->norm_d = PetscSqrtReal(norm_d);

588:     /*************************************************************************/
589:     /* Update objective function.                                            */
590:     /*************************************************************************/

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

594:     /*************************************************************************/
595:     /* Check that the preconditioner appears positive semidefinite.          */
596:     /*************************************************************************/

598:     rzm1 = rz;
599:     VecDot(r, z, &rz);             /* rz = r^T z        */
600:     if (rz < 0.0) {
601:       /***********************************************************************/
602:       /* The preconditioner is indefinite.                                   */
603:       /***********************************************************************/

605:       ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
606:       PetscInfo1(ksp, "KSPSolve_GLTR: cg indefinite preconditioner: rz=%g\n", (double)rz);
607:       break;
608:     }

610:     /*************************************************************************/
611:     /* As far as we know, the preconditioner is positive semidefinite.       */
612:     /* Compute the residual and check for convergence.                       */
613:     /*************************************************************************/

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

617:     switch (ksp->normtype) {
618:     case KSP_NORM_PRECONDITIONED:
619:       VecNorm(z, NORM_2, &norm_r); /* norm_r = |z|      */
620:       break;

622:     case KSP_NORM_UNPRECONDITIONED:
623:       VecNorm(r, NORM_2, &norm_r); /* norm_r = |r|      */
624:       break;

626:     case KSP_NORM_NATURAL:
627:       norm_r = cg->norm_r[l_size+1];                    /* norm_r = |r|_M    */
628:       break;

630:     default:
631:       norm_r = 0.0;
632:       break;
633:     }

635:     KSPLogResidualHistory(ksp, norm_r);
636:     KSPMonitor(ksp, ksp->its, norm_r);
637:     ksp->rnorm = norm_r;

639:     (*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP);
640:     if (ksp->reason) {
641:       /***********************************************************************/
642:       /* The method has converged.                                           */
643:       /***********************************************************************/

645:       PetscInfo2(ksp, "KSPSolve_GLTR: cg truncated step: rnorm=%g, radius=%g\n", (double)norm_r, (double)cg->radius);
646:       break;
647:     }

649:     /*************************************************************************/
650:     /* We have not converged yet.  Check for breakdown.                      */
651:     /*************************************************************************/

653:     beta = rz / rzm1;
654:     if (PetscAbsReal(beta) <= 0.0) {
655:       /***********************************************************************/
656:       /* Conjugate gradients has broken down.                                */
657:       /***********************************************************************/

659:       ksp->reason = KSP_DIVERGED_BREAKDOWN;
660:       PetscInfo1(ksp, "KSPSolve_GLTR: breakdown: beta=%g\n", (double)beta);
661:       break;
662:     }

664:     /*************************************************************************/
665:     /* Check iteration limit.                                                */
666:     /*************************************************************************/

668:     if (ksp->its >= max_cg_its) {
669:       ksp->reason = KSP_DIVERGED_ITS;
670:       PetscInfo1(ksp, "KSPSolve_GLTR: iterlim: its=%D\n", ksp->its);
671:       break;
672:     }

674:     /*************************************************************************/
675:     /* Update p and the norms.                                               */
676:     /*************************************************************************/

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

681:     switch (cg->dtype) {
682:     case GLTR_PRECONDITIONED_DIRECTION:
683:       dMp    = beta*(dMp + alpha*norm_p);
684:       norm_p = beta*(rzm1 + beta*norm_p);
685:       break;

687:     default:
688:       VecDot(d, p, &dMp);
689:       VecDot(p, p, &norm_p);
690:       break;
691:     }

693:     /*************************************************************************/
694:     /* Compute the new direction and update the iteration.                   */
695:     /*************************************************************************/

697:     KSP_MatMult(ksp, Qmat, p, z);  /* z = Q * p         */
698:     VecDot(p, z, &kappa);          /* kappa = p^T Q p   */
699:     ++ksp->its;

701:     /*************************************************************************/
702:     /* Update the Lanczos tridiagonal matrix.                            */
703:     /*************************************************************************/

705:     ++l_size;
706:     cg->offd[t_size] = PetscSqrtReal(beta) / PetscAbsReal(alpha);
707:     cg->diag[t_size] = kappa / rz + beta / alpha;
708:     ++t_size;

710:     /*************************************************************************/
711:     /* Check for breakdown of the conjugate gradient method; this occurs     */
712:     /* when kappa is zero.                                                   */
713:     /*************************************************************************/

715:     if (PetscAbsReal(kappa) <= 0.0) {
716:       /***********************************************************************/
717:       /* The method breaks down; move along the direction as if the matrix   */
718:       /* were indefinite.                                                    */
719:       /***********************************************************************/

721:       ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
722:       PetscInfo1(ksp, "KSPSolve_GLTR: cg breakdown: kappa=%g\n", (double)kappa);

724:       if (cg->radius && norm_p > 0.0) {
725:         /*********************************************************************/
726:         /* Follow direction to boundary.                                     */
727:         /*********************************************************************/

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

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

734:         /*********************************************************************/
735:         /* Update objective function.                                        */
736:         /*********************************************************************/

738:         cg->o_fcn += step * (0.5 * step * kappa - rz);
739:       } else if (cg->radius) {
740:         /*********************************************************************/
741:         /* The norm of the direction is zero; there is nothing to follow.    */
742:         /*********************************************************************/
743:       }
744:       break;
745:     }
746:   }

748:   /***************************************************************************/
749:   /* Check to see if we need to continue with the Lanczos method.            */
750:   /***************************************************************************/

752:   if (!cg->radius) {
753:     /*************************************************************************/
754:     /* There is no radius.  Therefore, we cannot move along the boundary.    */
755:     /*************************************************************************/
756:     return(0);
757:   }

759:   if (KSP_CONVERGED_CG_NEG_CURVE != ksp->reason) {
760:     /*************************************************************************/
761:     /* The method either converged to an interior point, hit the boundary of */
762:     /* the trust-region without encountering a direction of negative         */
763:     /* curvature or the iteration limit was reached.                         */
764:     /*************************************************************************/
765:     return(0);
766:   }

768:   /***************************************************************************/
769:   /* Switch to contructing the Lanczos basis by way of the conjugate         */
770:   /* directions.                                                             */
771:   /***************************************************************************/

773:   for (i = 0; i < max_lanczos_its; ++i) {
774:     /*************************************************************************/
775:     /* Check for breakdown of the conjugate gradient method; this occurs     */
776:     /* when kappa is zero.                                                   */
777:     /*************************************************************************/

779:     if (PetscAbsReal(kappa) <= 0.0) {
780:       ksp->reason = KSP_DIVERGED_BREAKDOWN;
781:       PetscInfo1(ksp, "KSPSolve_GLTR: lanczos breakdown: kappa=%g\n", (double)kappa);
782:       break;
783:     }

785:     /*************************************************************************/
786:     /* Update the direction and residual.                                    */
787:     /*************************************************************************/

789:     alpha             = rz / kappa;
790:     cg->alpha[l_size] = alpha;

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

795:     /*************************************************************************/
796:     /* Check that the preconditioner appears positive semidefinite.          */
797:     /*************************************************************************/

799:     rzm1 = rz;
800:     VecDot(r, z, &rz);             /* rz = r^T z        */
801:     if (rz < 0.0) {
802:       /***********************************************************************/
803:       /* The preconditioner is indefinite.                                   */
804:       /***********************************************************************/

806:       ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
807:       PetscInfo1(ksp, "KSPSolve_GLTR: lanczos indefinite preconditioner: rz=%g\n", (double)rz);
808:       break;
809:     }

811:     /*************************************************************************/
812:     /* As far as we know, the preconditioner is positive definite.  Compute  */
813:     /* the residual.  Do NOT check for convergence.                          */
814:     /*************************************************************************/

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

818:     switch (ksp->normtype) {
819:     case KSP_NORM_PRECONDITIONED:
820:       VecNorm(z, NORM_2, &norm_r); /* norm_r = |z|      */
821:       break;

823:     case KSP_NORM_UNPRECONDITIONED:
824:       VecNorm(r, NORM_2, &norm_r); /* norm_r = |r|      */
825:       break;

827:     case KSP_NORM_NATURAL:
828:       norm_r = cg->norm_r[l_size+1];                    /* norm_r = |r|_M    */
829:       break;

831:     default:
832:       norm_r = 0.0;
833:       break;
834:     }

836:     KSPLogResidualHistory(ksp, norm_r);
837:     KSPMonitor(ksp, ksp->its, norm_r);
838:     ksp->rnorm = norm_r;

840:     /*************************************************************************/
841:     /* Check for breakdown.                                                  */
842:     /*************************************************************************/

844:     beta = rz / rzm1;
845:     if (PetscAbsReal(beta) <= 0.0) {
846:       /***********************************************************************/
847:       /* Conjugate gradients has broken down.                                */
848:       /***********************************************************************/

850:       ksp->reason = KSP_DIVERGED_BREAKDOWN;
851:       PetscInfo1(ksp, "KSPSolve_GLTR: breakdown: beta=%g\n",(double) beta);
852:       break;
853:     }

855:     /*************************************************************************/
856:     /* Update p and the norms.                                               */
857:     /*************************************************************************/

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

862:     /*************************************************************************/
863:     /* Compute the new direction and update the iteration.                   */
864:     /*************************************************************************/

866:     KSP_MatMult(ksp, Qmat, p, z);  /* z = Q * p         */
867:     VecDot(p, z, &kappa);          /* kappa = p^T Q p   */
868:     ++ksp->its;

870:     /*************************************************************************/
871:     /* Update the Lanczos tridiagonal matrix.                                */
872:     /*************************************************************************/

874:     ++l_size;
875:     cg->offd[t_size] = PetscSqrtReal(beta) / PetscAbsReal(alpha);
876:     cg->diag[t_size] = kappa / rz + beta / alpha;
877:     ++t_size;
878:   }

880:   /***************************************************************************/
881:   /* We have the Lanczos basis, solve the tridiagonal trust-region problem   */
882:   /* to obtain the Lanczos direction.  We know that the solution lies on     */
883:   /* the boundary of the trust region.  We start by checking that the        */
884:   /* workspace allocated is large enough.                                    */
885:   /***************************************************************************/
886:   /* Note that the current version only computes the solution by using the   */
887:   /* preconditioned direction.  Need to think about how to do the            */
888:   /* unpreconditioned direction calculation.                                 */
889:   /***************************************************************************/

891:   if (t_size > cg->alloced) {
892:     if (cg->alloced) {
893:       PetscFree(cg->rwork);
894:       PetscFree(cg->iwork);
895:       cg->alloced += cg->init_alloc;
896:     } else {
897:       cg->alloced = cg->init_alloc;
898:     }

900:     while (t_size > cg->alloced) {
901:       cg->alloced += cg->init_alloc;
902:     }

904:     cg->alloced = PetscMin(cg->alloced, t_size);
905:     PetscMalloc2(10*cg->alloced, &cg->rwork,5*cg->alloced, &cg->iwork);
906:   }

908:   /***************************************************************************/
909:   /* Set up the required vectors.                                            */
910:   /***************************************************************************/

912:   t_soln = cg->rwork + 0*t_size;                        /* Solution          */
913:   t_diag = cg->rwork + 1*t_size;                        /* Diagonal of T     */
914:   t_offd = cg->rwork + 2*t_size;                        /* Off-diagonal of T */
915:   e_valu = cg->rwork + 3*t_size;                        /* Eigenvalues of T  */
916:   e_vect = cg->rwork + 4*t_size;                        /* Eigenvector of T  */
917:   e_rwrk = cg->rwork + 5*t_size;                        /* Eigen workspace   */

919:   e_iblk = cg->iwork + 0*t_size;                        /* Eigen blocks      */
920:   e_splt = cg->iwork + 1*t_size;                        /* Eigen splits      */
921:   e_iwrk = cg->iwork + 2*t_size;                        /* Eigen workspace   */

923:   /***************************************************************************/
924:   /* Compute the minimum eigenvalue of T.                                    */
925:   /***************************************************************************/

927:   vl = 0.0;
928:   vu = 0.0;
929:   il = 1;
930:   iu = 1;

932: #if defined(PETSC_MISSING_LAPACK_STEBZ)
933:   SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"STEBZ - Lapack routine is unavailable.");
934: #else
935:   PetscStackCallBLAS("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));

937:   if ((0 != info) || (1 != e_valus)) {
938:     /*************************************************************************/
939:     /* Calculation of the minimum eigenvalue failed.  Return the             */
940:     /* Steihaug-Toint direction.                                             */
941:     /*************************************************************************/

943:     PetscInfo(ksp, "KSPSolve_GLTR: failed to compute eigenvalue.\n");
944:     ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
945:     return(0);
946:   }

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

950:   /***************************************************************************/
951:   /* Compute the initial value of lambda to make (T + lamba I) positive      */
952:   /* definite.                                                               */
953:   /***************************************************************************/

955:   pert = cg->init_pert;
956:   if (e_valu[0] < 0.0) cg->lambda = pert - e_valu[0];
957: #endif

959:   while (1) {
960:     for (i = 0; i < t_size; ++i) {
961:       t_diag[i] = cg->diag[i] + cg->lambda;
962:       t_offd[i] = cg->offd[i];
963:     }

965: #if defined(PETSC_MISSING_LAPACK_PTTRF)
966:     SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"PTTRF - Lapack routine is unavailable.");
967: #else
968:     PetscStackCallBLAS("LAPACKpttrf",LAPACKpttrf_(&t_size, t_diag, t_offd + 1, &info));

970:     if (0 == info) break;

972:     pert      += pert;
973:     cg->lambda = cg->lambda * (1.0 + pert) + pert;
974: #endif
975:   }

977:   /***************************************************************************/
978:   /* Compute the initial step and its norm.                                  */
979:   /***************************************************************************/

981:   nrhs = 1;
982:   nldb = t_size;

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

987: #if defined(PETSC_MISSING_LAPACK_PTTRS)
988:   SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"PTTRS - Lapack routine is unavailable.");
989: #else
990:   PetscStackCallBLAS("LAPACKpttrs",LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, t_soln, &nldb, &info));
991: #endif

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

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

1004:   norm_t = 0.;
1005:   for (i = 0; i < t_size; ++i) norm_t += t_soln[i] * t_soln[i];
1006:   norm_t = PetscSqrtReal(norm_t);

1008:   /***************************************************************************/
1009:   /* Determine the case we are in.                                           */
1010:   /***************************************************************************/

1012:   if (norm_t <= cg->radius) {
1013:     /*************************************************************************/
1014:     /* The step is within the trust region; check if we are in the hard case */
1015:     /* and need to move to the boundary by following a direction of negative */
1016:     /* curvature.                                                            */
1017:     /*************************************************************************/

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

1025: #if defined(PETSC_MISSING_LAPACK_STEIN)
1026:       SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"STEIN - Lapack routine is unavailable.");
1027: #else
1028:       PetscStackCallBLAS("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));
1029: #endif

1031:       if (0 != info) {
1032:         /*********************************************************************/
1033:         /* Calculation of the minimum eigenvalue failed.  Return the         */
1034:         /* Steihaug-Toint direction.                                         */
1035:         /*********************************************************************/

1037:         PetscInfo(ksp, "KSPSolve_GLTR: failed to compute eigenvector.\n");
1038:         ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1039:         return(0);
1040:       }

1042:       coef1 = 0.0;
1043:       coef2 = 0.0;
1044:       coef3 = -cg->radius * cg->radius;
1045:       for (i = 0; i < t_size; ++i) {
1046:         coef1 += e_vect[i] * e_vect[i];
1047:         coef2 += e_vect[i] * t_soln[i];
1048:         coef3 += t_soln[i] * t_soln[i];
1049:       }

1051:       coef3 = PetscSqrtReal(coef2 * coef2 - coef1 * coef3);
1052:       root1 = (-coef2 + coef3) / coef1;
1053:       root2 = (-coef2 - coef3) / coef1;

1055:       /***********************************************************************/
1056:       /* Compute objective value for (t_soln + root1 * e_vect)               */
1057:       /***********************************************************************/

1059:       for (i = 0; i < t_size; ++i) {
1060:         e_rwrk[i] = t_soln[i] + root1 * e_vect[i];
1061:       }

1063:       obj1 = e_rwrk[0]*(0.5*(cg->diag[0]*e_rwrk[0]+
1064:                              cg->offd[1]*e_rwrk[1])+cg->norm_r[0]);
1065:       for (i = 1; i < t_size - 1; ++i) {
1066:         obj1 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
1067:                                cg->diag[i]*e_rwrk[i]+
1068:                                cg->offd[i+1]*e_rwrk[i+1]);
1069:       }
1070:       obj1 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
1071:                              cg->diag[i]*e_rwrk[i]);

1073:       /***********************************************************************/
1074:       /* Compute objective value for (t_soln + root2 * e_vect)               */
1075:       /***********************************************************************/

1077:       for (i = 0; i < t_size; ++i) {
1078:         e_rwrk[i] = t_soln[i] + root2 * e_vect[i];
1079:       }

1081:       obj2 = e_rwrk[0]*(0.5*(cg->diag[0]*e_rwrk[0]+
1082:                              cg->offd[1]*e_rwrk[1])+cg->norm_r[0]);
1083:       for (i = 1; i < t_size - 1; ++i) {
1084:         obj2 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
1085:                                cg->diag[i]*e_rwrk[i]+
1086:                                cg->offd[i+1]*e_rwrk[i+1]);
1087:       }
1088:       obj2 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
1089:                              cg->diag[i]*e_rwrk[i]);

1091:       /***********************************************************************/
1092:       /* Choose the point with the best objective function value.            */
1093:       /***********************************************************************/

1095:       if (obj1 <= obj2) {
1096:         for (i = 0; i < t_size; ++i) {
1097:           t_soln[i] += root1 * e_vect[i];
1098:         }
1099:       }
1100:       else {
1101:         for (i = 0; i < t_size; ++i) {
1102:           t_soln[i] += root2 * e_vect[i];
1103:         }
1104:       }
1105:     } else {
1106:       /***********************************************************************/
1107:       /* The matrix is positive definite or there was no room to move; the   */
1108:       /* solution is already contained in t_soln.                            */
1109:       /***********************************************************************/
1110:     }
1111:   } else {
1112:     /*************************************************************************/
1113:     /* The step is outside the trust-region.  Compute the correct value for  */
1114:     /* lambda by performing Newton's method.                                 */
1115:     /*************************************************************************/

1117:     for (i = 0; i < max_newton_its; ++i) {
1118:       /***********************************************************************/
1119:       /* Check for convergence.                                              */
1120:       /***********************************************************************/

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

1124:       /***********************************************************************/
1125:       /* Compute the update.                                                 */
1126:       /***********************************************************************/

1128:       PetscMemcpy(e_rwrk, t_soln, sizeof(PetscReal)*t_size);

1130: #if defined(PETSC_MISSING_LAPACK_PTTRS)
1131:       SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"PTTRS - Lapack routine is unavailable.");
1132: #else
1133:       PetscStackCallBLAS("LAPACKpttrs",LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, e_rwrk, &nldb, &info));
1134: #endif

1136:       if (0 != info) {
1137:         /*********************************************************************/
1138:         /* Calculation of the step failed; return the Steihaug-Toint         */
1139:         /* direction.                                                        */
1140:         /*********************************************************************/

1142:         PetscInfo(ksp, "KSPSolve_GLTR: failed to compute step.\n");
1143:         ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1144:         return(0);
1145:       }

1147:       /***********************************************************************/
1148:       /* Modify lambda.                                                      */
1149:       /***********************************************************************/

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

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

1156:       /***********************************************************************/
1157:       /* Factor T + lambda I                                                 */
1158:       /***********************************************************************/

1160:       for (j = 0; j < t_size; ++j) {
1161:         t_diag[j] = cg->diag[j] + cg->lambda;
1162:         t_offd[j] = cg->offd[j];
1163:       }

1165: #if defined(PETSC_MISSING_LAPACK_PTTRF)
1166:       SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"PTTRF - Lapack routine is unavailable.");
1167: #else
1168:       PetscStackCallBLAS("LAPACKpttrf",LAPACKpttrf_(&t_size, t_diag, t_offd + 1, &info));
1169: #endif

1171:       if (0 != info) {
1172:         /*********************************************************************/
1173:         /* Calculation of factorization failed; return the Steihaug-Toint    */
1174:         /* direction.                                                        */
1175:         /*********************************************************************/

1177:         PetscInfo(ksp, "KSPSolve_GLTR: factorization failed.\n");
1178:         ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1179:         return(0);
1180:       }

1182:       /***********************************************************************/
1183:       /* Compute the new step and its norm.                                  */
1184:       /***********************************************************************/

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

1189: #if defined(PETSC_MISSING_LAPACK_PTTRS)
1190:       SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"PTTRS - Lapack routine is unavailable.");
1191: #else
1192:       PetscStackCallBLAS("LAPACKpttrs",LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, t_soln, &nldb, &info));
1193: #endif

1195:       if (0 != info) {
1196:         /*********************************************************************/
1197:         /* Calculation of the step failed; return the Steihaug-Toint         */
1198:         /* direction.                                                        */
1199:         /*********************************************************************/

1201:         PetscInfo(ksp, "KSPSolve_GLTR: failed to compute step.\n");
1202:         ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1203:         return(0);
1204:       }

1206:       norm_t = 0.;
1207:       for (j = 0; j < t_size; ++j) norm_t += t_soln[j] * t_soln[j];
1208:       norm_t = PetscSqrtReal(norm_t);
1209:     }

1211:     /*************************************************************************/
1212:     /* Check for convergence.                                                */
1213:     /*************************************************************************/

1215:     if (PetscAbsReal(norm_t - cg->radius) > cg->newton_tol * cg->radius) {
1216:       /***********************************************************************/
1217:       /* Newton method failed to converge in iteration limit.                */
1218:       /***********************************************************************/

1220:       PetscInfo(ksp, "KSPSolve_GLTR: failed to converge.\n");
1221:       ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1222:       return(0);
1223:     }
1224:   }

1226:   /***************************************************************************/
1227:   /* Recover the norm of the direction and objective function value.         */
1228:   /***************************************************************************/

1230:   cg->norm_d = norm_t;

1232:   cg->o_fcn = t_soln[0]*(0.5*(cg->diag[0]*t_soln[0]+cg->offd[1]*t_soln[1])+cg->norm_r[0]);
1233:   for (i = 1; i < t_size - 1; ++i) {
1234:     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]);
1235:   }
1236:   cg->o_fcn += 0.5*t_soln[i]*(cg->offd[i]*t_soln[i-1]+cg->diag[i]*t_soln[i]);

1238:   /***************************************************************************/
1239:   /* Recover the direction.                                                  */
1240:   /***************************************************************************/

1242:   sigma = -1;

1244:   /***************************************************************************/
1245:   /* Start conjugate gradient method from the beginning                      */
1246:   /***************************************************************************/

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

1251:   /***************************************************************************/
1252:   /* Accumulate Q * s                                                        */
1253:   /***************************************************************************/

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

1258:   /***************************************************************************/
1259:   /* Compute the first direction.                                            */
1260:   /***************************************************************************/

1262:   VecCopy(z, p);                   /* p = z             */
1263:   KSP_MatMult(ksp, Qmat, p, z);    /* z = Q * p         */
1264:   ++ksp->its;

1266:   for (i = 0; i < l_size - 1; ++i) {
1267:     /*************************************************************************/
1268:     /* Update the residual and direction.                                    */
1269:     /*************************************************************************/

1271:     alpha = cg->alpha[i];
1272:     if (alpha >= 0.0) sigma = -sigma;

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

1277:     /*************************************************************************/
1278:     /* Accumulate Q * s                                                      */
1279:     /*************************************************************************/

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

1283:     /*************************************************************************/
1284:     /* Update p.                                                             */
1285:     /*************************************************************************/

1287:     beta = cg->beta[i];
1288:     VecAYPX(p, beta, z);          /* p = z + beta p    */
1289:     KSP_MatMult(ksp, Qmat, p, z);  /* z = Q * p         */
1290:     ++ksp->its;
1291:   }

1293:   /***************************************************************************/
1294:   /* Update the residual and direction.                                      */
1295:   /***************************************************************************/

1297:   alpha = cg->alpha[i];
1298:   if (alpha >= 0.0) sigma = -sigma;

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

1303:   /***************************************************************************/
1304:   /* Accumulate Q * s                                                        */
1305:   /***************************************************************************/

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

1309:   /***************************************************************************/
1310:   /* Set the termination reason.                                             */
1311:   /***************************************************************************/

1313:   ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1314:   return(0);
1315: #endif
1316: }

1320: PetscErrorCode KSPSetUp_GLTR(KSP ksp)
1321: {
1322:   KSP_GLTR       *cg = (KSP_GLTR*)ksp->data;
1323:   PetscInt       max_its;

1327:   /***************************************************************************/
1328:   /* Determine the total maximum number of iterations.                       */
1329:   /***************************************************************************/

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

1333:   /***************************************************************************/
1334:   /* Set work vectors needed by conjugate gradient method and allocate       */
1335:   /* workspace for Lanczos matrix.                                           */
1336:   /***************************************************************************/

1338:   KSPSetWorkVecs(ksp, 3);

1340:   PetscCalloc5(max_its,&cg->diag,max_its,&cg->offd,max_its,&cg->alpha,max_its,&cg->beta,max_its,&cg->norm_r);
1341:   PetscLogObjectMemory((PetscObject)ksp, 5*max_its*sizeof(PetscReal));
1342:   return(0);
1343: }

1347: PetscErrorCode KSPDestroy_GLTR(KSP ksp)
1348: {
1349:   KSP_GLTR       *cg = (KSP_GLTR*)ksp->data;

1353:   /***************************************************************************/
1354:   /* Free memory allocated for the data.                                     */
1355:   /***************************************************************************/

1357:   PetscFree5(cg->diag,cg->offd,cg->alpha,cg->beta,cg->norm_r);
1358:   if (cg->alloced) {
1359:     PetscFree2(cg->rwork,cg->iwork);
1360:   }

1362:   /***************************************************************************/
1363:   /* Clear composed functions                                                */
1364:   /***************************************************************************/

1366:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRSetRadius_C",NULL);
1367:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetNormD_C",NULL);
1368:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetObjFcn_C",NULL);
1369:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetMinEig_C",NULL);
1370:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetLambda_C",NULL);

1372:   /***************************************************************************/
1373:   /* Destroy KSP object.                                                     */
1374:   /***************************************************************************/
1375:   KSPDestroyDefault(ksp);
1376:   return(0);
1377: }

1381: static PetscErrorCode  KSPGLTRSetRadius_GLTR(KSP ksp, PetscReal radius)
1382: {
1383:   KSP_GLTR *cg = (KSP_GLTR*)ksp->data;

1386:   cg->radius = radius;
1387:   return(0);
1388: }

1392: static PetscErrorCode  KSPGLTRGetNormD_GLTR(KSP ksp, PetscReal *norm_d)
1393: {
1394:   KSP_GLTR *cg = (KSP_GLTR*)ksp->data;

1397:   *norm_d = cg->norm_d;
1398:   return(0);
1399: }

1403: static PetscErrorCode  KSPGLTRGetObjFcn_GLTR(KSP ksp, PetscReal *o_fcn)
1404: {
1405:   KSP_GLTR *cg = (KSP_GLTR*)ksp->data;

1408:   *o_fcn = cg->o_fcn;
1409:   return(0);
1410: }

1414: static PetscErrorCode  KSPGLTRGetMinEig_GLTR(KSP ksp, PetscReal *e_min)
1415: {
1416:   KSP_GLTR *cg = (KSP_GLTR*)ksp->data;

1419:   *e_min = cg->e_min;
1420:   return(0);
1421: }

1425: static PetscErrorCode  KSPGLTRGetLambda_GLTR(KSP ksp, PetscReal *lambda)
1426: {
1427:   KSP_GLTR *cg = (KSP_GLTR*)ksp->data;

1430:   *lambda = cg->lambda;
1431:   return(0);
1432: }

1436: PetscErrorCode KSPSetFromOptions_GLTR(KSP ksp)
1437: {
1439:   KSP_GLTR       *cg = (KSP_GLTR*)ksp->data;

1442:   PetscOptionsHead("KSP GLTR options");

1444:   PetscOptionsReal("-ksp_gltr_radius", "Trust Region Radius", "KSPGLTRSetRadius", cg->radius, &cg->radius, NULL);

1446:   PetscOptionsReal("-ksp_gltr_init_pert", "Initial perturbation", "", cg->init_pert, &cg->init_pert, NULL);
1447:   PetscOptionsReal("-ksp_gltr_eigen_tol", "Eigenvalue tolerance", "", cg->eigen_tol, &cg->eigen_tol, NULL);
1448:   PetscOptionsReal("-ksp_gltr_newton_tol", "Newton tolerance", "", cg->newton_tol, &cg->newton_tol, NULL);

1450:   PetscOptionsInt("-ksp_gltr_max_lanczos_its", "Maximum Lanczos Iters", "", cg->max_lanczos_its, &cg->max_lanczos_its, NULL);
1451:   PetscOptionsInt("-ksp_gltr_max_newton_its", "Maximum Newton Iters", "", cg->max_newton_its, &cg->max_newton_its, NULL);

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

1455:   PetscOptionsTail();
1456:   return(0);
1457: }

1459: /*MC
1460:      KSPGLTR -   Code to run conjugate gradient method subject to a constraint
1461:          on the solution norm. This is used in Trust Region methods for
1462:          nonlinear equations, SNESNEWTONTR

1464:    Options Database Keys:
1465: .      -ksp_gltr_radius <r> - Trust Region Radius

1467:    Notes: This is rarely used directly

1469:   Use preconditioned conjugate gradient to compute
1470:   an approximate minimizer of the quadratic function

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

1474:    subject to the trust region constraint

1476:             || s || <= delta,

1478:    where

1480:      delta is the trust region radius,
1481:      g is the gradient vector,
1482:      H is the Hessian approximation,
1483:      M is the positive definite preconditioner matrix.

1485:    KSPConvergedReason may be
1486: $  KSP_CONVERGED_CG_NEG_CURVE if convergence is reached along a negative curvature direction,
1487: $  KSP_CONVERGED_CG_CONSTRAINED if convergence is reached along a constrained step,
1488: $  other KSP converged/diverged reasons

1490:   Notes:
1491:   The preconditioner supplied should be symmetric and positive definite.

1493:    Level: developer

1495: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPGLTRSetRadius(), KSPGLTRGetNormD(), KSPGLTRGetObjFcn(), KSPGLTRGetMinEig(), KSPGLTRGetLambda()
1496: M*/

1500: PETSC_EXTERN PetscErrorCode KSPCreate_GLTR(KSP ksp)
1501: {
1503:   KSP_GLTR       *cg;

1506:   PetscNewLog(ksp,&cg);
1507:   cg->radius = 0.0;
1508:   cg->dtype  = GLTR_UNPRECONDITIONED_DIRECTION;

1510:   cg->init_pert  = 1.0e-8;
1511:   cg->eigen_tol  = 1.0e-10;
1512:   cg->newton_tol = 1.0e-6;

1514:   cg->alloced    = 0;
1515:   cg->init_alloc = 1024;

1517:   cg->max_lanczos_its = 20;
1518:   cg->max_newton_its  = 10;

1520:   ksp->data = (void*) cg;
1521:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,3);
1522:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);
1523:   KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,2);

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

1530:   ksp->ops->setup          = KSPSetUp_GLTR;
1531:   ksp->ops->solve          = KSPSolve_GLTR;
1532:   ksp->ops->destroy        = KSPDestroy_GLTR;
1533:   ksp->ops->setfromoptions = KSPSetFromOptions_GLTR;
1534:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
1535:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
1536:   ksp->ops->view           = 0;

1538:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRSetRadius_C",KSPGLTRSetRadius_GLTR);
1539:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetNormD_C", KSPGLTRGetNormD_GLTR);
1540:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetObjFcn_C",KSPGLTRGetObjFcn_GLTR);
1541:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetMinEig_C",KSPGLTRGetMinEig_GLTR);
1542:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetLambda_C",KSPGLTRGetLambda_GLTR);
1543:   return(0);
1544: }