Actual source code: gltr.c

petsc-3.6.0 2015-06-09
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:   KSPCheckDot(ksp,rr);

226:   /***************************************************************************/
227:   /* Check the preconditioner for numerical problems and for positive        */
228:   /* definiteness.  The check for not-a-number and infinite values need be   */
229:   /* performed only once.                                                    */
230:   /***************************************************************************/

232:   KSP_PCApply(ksp, r, z);          /* z = inv(M) r      */
233:   VecDot(r, z, &rz);               /* rz = r^T inv(M) r */
234:   if (PetscIsInfOrNanScalar(rz)) {
235:     /*************************************************************************/
236:     /* The preconditioner contains not-a-number or an infinite value.        */
237:     /* Return the gradient direction intersected with the trust region.      */
238:     /*************************************************************************/

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

243:     if (cg->radius) {
244:       if (r2 >= rr) {
245:         alpha      = 1.0;
246:         cg->norm_d = PetscSqrtReal(rr);
247:       } else {
248:         alpha      = PetscSqrtReal(r2 / rr);
249:         cg->norm_d = cg->radius;
250:       }

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

254:       /***********************************************************************/
255:       /* Compute objective function.                                         */
256:       /***********************************************************************/

258:       KSP_MatMult(ksp, Qmat, d, z);
259:       VecAYPX(z, -0.5, ksp->vec_rhs);
260:       VecDot(d, z, &cg->o_fcn);
261:       cg->o_fcn = -cg->o_fcn;
262:       ++ksp->its;
263:     }
264:     return(0);
265:   }

267:   if (rz < 0.0) {
268:     /*************************************************************************/
269:     /* The preconditioner is indefinite.  Because this is the first          */
270:     /* and we do not have a direction yet, we use the gradient step.  Note   */
271:     /* that we cannot use the preconditioned norm when computing the step    */
272:     /* because the matrix is indefinite.                                     */
273:     /*************************************************************************/

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

278:     if (cg->radius) {
279:       if (r2 >= rr) {
280:         alpha      = 1.0;
281:         cg->norm_d = PetscSqrtReal(rr);
282:       } else {
283:         alpha      = PetscSqrtReal(r2 / rr);
284:         cg->norm_d = cg->radius;
285:       }

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

289:       /***********************************************************************/
290:       /* Compute objective function.                                         */
291:       /***********************************************************************/

293:       KSP_MatMult(ksp, Qmat, d, z);
294:       VecAYPX(z, -0.5, ksp->vec_rhs);
295:       VecDot(d, z, &cg->o_fcn);
296:       cg->o_fcn = -cg->o_fcn;
297:       ++ksp->its;
298:     }
299:     return(0);
300:   }

302:   /***************************************************************************/
303:   /* As far as we know, the preconditioner is positive semidefinite.         */
304:   /* Compute and log the residual.  Check convergence because this           */
305:   /* initializes things, but do not terminate until at least one conjugate   */
306:   /* gradient iteration has been performed.                                  */
307:   /***************************************************************************/

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

311:   switch (ksp->normtype) {
312:   case KSP_NORM_PRECONDITIONED:
313:     VecNorm(z, NORM_2, &norm_r);   /* norm_r = |z|      */
314:     break;

316:   case KSP_NORM_UNPRECONDITIONED:
317:     norm_r = PetscSqrtReal(rr);                                 /* norm_r = |r|      */
318:     break;

320:   case KSP_NORM_NATURAL:
321:     norm_r = cg->norm_r[0];                             /* norm_r = |r|_M    */
322:     break;

324:   default:
325:     norm_r = 0.0;
326:     break;
327:   }

329:   KSPLogResidualHistory(ksp, norm_r);
330:   KSPMonitor(ksp, ksp->its, norm_r);
331:   ksp->rnorm = norm_r;

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

335:   /***************************************************************************/
336:   /* Compute the first direction and update the iteration.                   */
337:   /***************************************************************************/

339:   VecCopy(z, p);                   /* p = z             */
340:   KSP_MatMult(ksp, Qmat, p, z);    /* z = Q * p         */
341:   ++ksp->its;

343:   /***************************************************************************/
344:   /* Check the matrix for numerical problems.                                */
345:   /***************************************************************************/

347:   VecDot(p, z, &kappa);            /* kappa = p^T Q p   */
348:   if (PetscIsInfOrNanScalar(kappa)) {
349:     /*************************************************************************/
350:     /* The matrix produced not-a-number or an infinite value.  In this case, */
351:     /* we must stop and use the gradient direction.  This condition need     */
352:     /* only be checked once.                                                 */
353:     /*************************************************************************/

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

358:     if (cg->radius) {
359:       if (r2 >= rr) {
360:         alpha      = 1.0;
361:         cg->norm_d = PetscSqrtReal(rr);
362:       } else {
363:         alpha      = PetscSqrtReal(r2 / rr);
364:         cg->norm_d = cg->radius;
365:       }

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

369:       /***********************************************************************/
370:       /* Compute objective function.                                         */
371:       /***********************************************************************/

373:       KSP_MatMult(ksp, Qmat, d, z);
374:       VecAYPX(z, -0.5, ksp->vec_rhs);
375:       VecDot(d, z, &cg->o_fcn);
376:       cg->o_fcn = -cg->o_fcn;
377:       ++ksp->its;
378:     }
379:     return(0);
380:   }

382:   /***************************************************************************/
383:   /* Initialize variables for calculating the norm of the direction and for  */
384:   /* the Lanczos tridiagonal matrix.  Note that we track the diagonal value  */
385:   /* of the Cholesky factorization of the Lanczos matrix in order to         */
386:   /* determine when negative curvature is encountered.                       */
387:   /***************************************************************************/

389:   dMp    = 0.0;
390:   norm_d = 0.0;
391:   switch (cg->dtype) {
392:   case GLTR_PRECONDITIONED_DIRECTION:
393:     norm_p = rz;
394:     break;

396:   default:
397:     VecDot(p, p, &norm_p);
398:     break;
399:   }

401:   cg->diag[t_size] = kappa / rz;
402:   cg->offd[t_size] = 0.0;
403:   ++t_size;

405:   piv = 1.0;

407:   /***************************************************************************/
408:   /* Check for breakdown of the conjugate gradient method; this occurs when  */
409:   /* kappa is zero.                                                          */
410:   /***************************************************************************/

412:   if (PetscAbsReal(kappa) <= 0.0) {
413:     /*************************************************************************/
414:     /* The curvature is zero.  In this case, we must stop and use follow     */
415:     /* the direction of negative curvature since the Lanczos matrix is zero. */
416:     /*************************************************************************/

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

421:     if (cg->radius && norm_p > 0.0) {
422:       /***********************************************************************/
423:       /* Follow direction of negative curvature to the boundary of the       */
424:       /* trust region.                                                       */
425:       /***********************************************************************/

427:       step       = PetscSqrtReal(r2 / norm_p);
428:       cg->norm_d = cg->radius;

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

432:       /***********************************************************************/
433:       /* Update objective function.                                          */
434:       /***********************************************************************/

436:       cg->o_fcn += step * (0.5 * step * kappa - rz);
437:     } else if (cg->radius) {
438:       /***********************************************************************/
439:       /* The norm of the preconditioned direction is zero; use the gradient  */
440:       /* step.                                                               */
441:       /***********************************************************************/

443:       if (r2 >= rr) {
444:         alpha      = 1.0;
445:         cg->norm_d = PetscSqrtReal(rr);
446:       } else {
447:         alpha      = PetscSqrtReal(r2 / rr);
448:         cg->norm_d = cg->radius;
449:       }

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

453:       /***********************************************************************/
454:       /* Compute objective function.                                         */
455:       /***********************************************************************/

457:       KSP_MatMult(ksp, Qmat, d, z);
458:       VecAYPX(z, -0.5, ksp->vec_rhs);
459:       VecDot(d, z, &cg->o_fcn);
460:       cg->o_fcn = -cg->o_fcn;
461:       ++ksp->its;
462:     }
463:     return(0);
464:   }

466:   /***************************************************************************/
467:   /* Begin the first part of the GLTR algorithm which runs the conjugate     */
468:   /* gradient method until either the problem is solved, we encounter the    */
469:   /* boundary of the trust region, or the conjugate gradient method breaks   */
470:   /* down.                                                                   */
471:   /***************************************************************************/

473:   while (1) {
474:     /*************************************************************************/
475:     /* Know that kappa is nonzero, because we have not broken down, so we    */
476:     /* can compute the steplength.                                           */
477:     /*************************************************************************/

479:     alpha             = rz / kappa;
480:     cg->alpha[l_size] = alpha;

482:     /*************************************************************************/
483:     /* Compute the diagonal value of the Cholesky factorization of the       */
484:     /* Lanczos matrix and check to see if the Lanczos matrix is indefinite.  */
485:     /* This indicates a direction of negative curvature.                     */
486:     /*************************************************************************/

488:     piv = cg->diag[l_size] - cg->offd[l_size]*cg->offd[l_size] / piv;
489:     if (piv <= 0.0) {
490:       /***********************************************************************/
491:       /* In this case, the matrix is indefinite and we have encountered      */
492:       /* a direction of negative curvature.  Follow the direction to the     */
493:       /* boundary of the trust region.                                       */
494:       /***********************************************************************/

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

499:       if (cg->radius && norm_p > 0.0) {
500:         /*********************************************************************/
501:         /* Follow direction of negative curvature to boundary.               */
502:         /*********************************************************************/

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

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

509:         /*********************************************************************/
510:         /* Update objective function.                                        */
511:         /*********************************************************************/

513:         cg->o_fcn += step * (0.5 * step * kappa - rz);
514:       } else if (cg->radius) {
515:         /*********************************************************************/
516:         /* The norm of the direction is zero; there is nothing to follow.    */
517:         /*********************************************************************/
518:       }
519:       break;
520:     }

522:     /*************************************************************************/
523:     /* Compute the steplength and check for intersection with the trust      */
524:     /* region.                                                               */
525:     /*************************************************************************/

527:     norm_dp1 = norm_d + alpha*(2.0*dMp + alpha*norm_p);
528:     if (cg->radius && norm_dp1 >= r2) {
529:       /***********************************************************************/
530:       /* In this case, the matrix is positive definite as far as we know.    */
531:       /* However, the full step goes beyond the trust region.                */
532:       /***********************************************************************/

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

537:       if (norm_p > 0.0) {
538:         /*********************************************************************/
539:         /* Follow the direction to the boundary of the trust region.         */
540:         /*********************************************************************/

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

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

547:         /*********************************************************************/
548:         /* Update objective function.                                        */
549:         /*********************************************************************/

551:         cg->o_fcn += step * (0.5 * step * kappa - rz);
552:       } else {
553:         /*********************************************************************/
554:         /* The norm of the direction is zero; there is nothing to follow.    */
555:         /*********************************************************************/
556:       }
557:       break;
558:     }

560:     /*************************************************************************/
561:     /* Now we can update the direction and residual.                         */
562:     /*************************************************************************/

564:     VecAXPY(d, alpha, p);          /* d = d + alpha p   */
565:     VecAXPY(r, -alpha, z);         /* r = r - alpha Q p */
566:     KSP_PCApply(ksp, r, z);        /* z = inv(M) r      */

568:     switch (cg->dtype) {
569:     case GLTR_PRECONDITIONED_DIRECTION:
570:       norm_d = norm_dp1;
571:       break;

573:     default:
574:       VecDot(d, d, &norm_d);
575:       break;
576:     }
577:     cg->norm_d = PetscSqrtReal(norm_d);

579:     /*************************************************************************/
580:     /* Update objective function.                                            */
581:     /*************************************************************************/

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

585:     /*************************************************************************/
586:     /* Check that the preconditioner appears positive semidefinite.          */
587:     /*************************************************************************/

589:     rzm1 = rz;
590:     VecDot(r, z, &rz);             /* rz = r^T z        */
591:     if (rz < 0.0) {
592:       /***********************************************************************/
593:       /* The preconditioner is indefinite.                                   */
594:       /***********************************************************************/

596:       ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
597:       PetscInfo1(ksp, "KSPSolve_GLTR: cg indefinite preconditioner: rz=%g\n", (double)rz);
598:       break;
599:     }

601:     /*************************************************************************/
602:     /* As far as we know, the preconditioner is positive semidefinite.       */
603:     /* Compute the residual and check for convergence.                       */
604:     /*************************************************************************/

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

608:     switch (ksp->normtype) {
609:     case KSP_NORM_PRECONDITIONED:
610:       VecNorm(z, NORM_2, &norm_r); /* norm_r = |z|      */
611:       break;

613:     case KSP_NORM_UNPRECONDITIONED:
614:       VecNorm(r, NORM_2, &norm_r); /* norm_r = |r|      */
615:       break;

617:     case KSP_NORM_NATURAL:
618:       norm_r = cg->norm_r[l_size+1];                    /* norm_r = |r|_M    */
619:       break;

621:     default:
622:       norm_r = 0.0;
623:       break;
624:     }

626:     KSPLogResidualHistory(ksp, norm_r);
627:     KSPMonitor(ksp, ksp->its, norm_r);
628:     ksp->rnorm = norm_r;

630:     (*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP);
631:     if (ksp->reason) {
632:       /***********************************************************************/
633:       /* The method has converged.                                           */
634:       /***********************************************************************/

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

640:     /*************************************************************************/
641:     /* We have not converged yet.  Check for breakdown.                      */
642:     /*************************************************************************/

644:     beta = rz / rzm1;
645:     if (PetscAbsReal(beta) <= 0.0) {
646:       /***********************************************************************/
647:       /* Conjugate gradients has broken down.                                */
648:       /***********************************************************************/

650:       ksp->reason = KSP_DIVERGED_BREAKDOWN;
651:       PetscInfo1(ksp, "KSPSolve_GLTR: breakdown: beta=%g\n", (double)beta);
652:       break;
653:     }

655:     /*************************************************************************/
656:     /* Check iteration limit.                                                */
657:     /*************************************************************************/

659:     if (ksp->its >= max_cg_its) {
660:       ksp->reason = KSP_DIVERGED_ITS;
661:       PetscInfo1(ksp, "KSPSolve_GLTR: iterlim: its=%D\n", ksp->its);
662:       break;
663:     }

665:     /*************************************************************************/
666:     /* Update p and the norms.                                               */
667:     /*************************************************************************/

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

672:     switch (cg->dtype) {
673:     case GLTR_PRECONDITIONED_DIRECTION:
674:       dMp    = beta*(dMp + alpha*norm_p);
675:       norm_p = beta*(rzm1 + beta*norm_p);
676:       break;

678:     default:
679:       VecDot(d, p, &dMp);
680:       VecDot(p, p, &norm_p);
681:       break;
682:     }

684:     /*************************************************************************/
685:     /* Compute the new direction and update the iteration.                   */
686:     /*************************************************************************/

688:     KSP_MatMult(ksp, Qmat, p, z);  /* z = Q * p         */
689:     VecDot(p, z, &kappa);          /* kappa = p^T Q p   */
690:     ++ksp->its;

692:     /*************************************************************************/
693:     /* Update the Lanczos tridiagonal matrix.                            */
694:     /*************************************************************************/

696:     ++l_size;
697:     cg->offd[t_size] = PetscSqrtReal(beta) / PetscAbsReal(alpha);
698:     cg->diag[t_size] = kappa / rz + beta / alpha;
699:     ++t_size;

701:     /*************************************************************************/
702:     /* Check for breakdown of the conjugate gradient method; this occurs     */
703:     /* when kappa is zero.                                                   */
704:     /*************************************************************************/

706:     if (PetscAbsReal(kappa) <= 0.0) {
707:       /***********************************************************************/
708:       /* The method breaks down; move along the direction as if the matrix   */
709:       /* were indefinite.                                                    */
710:       /***********************************************************************/

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

715:       if (cg->radius && norm_p > 0.0) {
716:         /*********************************************************************/
717:         /* Follow direction to boundary.                                     */
718:         /*********************************************************************/

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

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

725:         /*********************************************************************/
726:         /* Update objective function.                                        */
727:         /*********************************************************************/

729:         cg->o_fcn += step * (0.5 * step * kappa - rz);
730:       } else if (cg->radius) {
731:         /*********************************************************************/
732:         /* The norm of the direction is zero; there is nothing to follow.    */
733:         /*********************************************************************/
734:       }
735:       break;
736:     }
737:   }

739:   /***************************************************************************/
740:   /* Check to see if we need to continue with the Lanczos method.            */
741:   /***************************************************************************/

743:   if (!cg->radius) {
744:     /*************************************************************************/
745:     /* There is no radius.  Therefore, we cannot move along the boundary.    */
746:     /*************************************************************************/
747:     return(0);
748:   }

750:   if (KSP_CONVERGED_CG_NEG_CURVE != ksp->reason) {
751:     /*************************************************************************/
752:     /* The method either converged to an interior point, hit the boundary of */
753:     /* the trust-region without encountering a direction of negative         */
754:     /* curvature or the iteration limit was reached.                         */
755:     /*************************************************************************/
756:     return(0);
757:   }

759:   /***************************************************************************/
760:   /* Switch to contructing the Lanczos basis by way of the conjugate         */
761:   /* directions.                                                             */
762:   /***************************************************************************/

764:   for (i = 0; i < max_lanczos_its; ++i) {
765:     /*************************************************************************/
766:     /* Check for breakdown of the conjugate gradient method; this occurs     */
767:     /* when kappa is zero.                                                   */
768:     /*************************************************************************/

770:     if (PetscAbsReal(kappa) <= 0.0) {
771:       ksp->reason = KSP_DIVERGED_BREAKDOWN;
772:       PetscInfo1(ksp, "KSPSolve_GLTR: lanczos breakdown: kappa=%g\n", (double)kappa);
773:       break;
774:     }

776:     /*************************************************************************/
777:     /* Update the direction and residual.                                    */
778:     /*************************************************************************/

780:     alpha             = rz / kappa;
781:     cg->alpha[l_size] = alpha;

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

786:     /*************************************************************************/
787:     /* Check that the preconditioner appears positive semidefinite.          */
788:     /*************************************************************************/

790:     rzm1 = rz;
791:     VecDot(r, z, &rz);             /* rz = r^T z        */
792:     if (rz < 0.0) {
793:       /***********************************************************************/
794:       /* The preconditioner is indefinite.                                   */
795:       /***********************************************************************/

797:       ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
798:       PetscInfo1(ksp, "KSPSolve_GLTR: lanczos indefinite preconditioner: rz=%g\n", (double)rz);
799:       break;
800:     }

802:     /*************************************************************************/
803:     /* As far as we know, the preconditioner is positive definite.  Compute  */
804:     /* the residual.  Do NOT check for convergence.                          */
805:     /*************************************************************************/

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

809:     switch (ksp->normtype) {
810:     case KSP_NORM_PRECONDITIONED:
811:       VecNorm(z, NORM_2, &norm_r); /* norm_r = |z|      */
812:       break;

814:     case KSP_NORM_UNPRECONDITIONED:
815:       VecNorm(r, NORM_2, &norm_r); /* norm_r = |r|      */
816:       break;

818:     case KSP_NORM_NATURAL:
819:       norm_r = cg->norm_r[l_size+1];                    /* norm_r = |r|_M    */
820:       break;

822:     default:
823:       norm_r = 0.0;
824:       break;
825:     }

827:     KSPLogResidualHistory(ksp, norm_r);
828:     KSPMonitor(ksp, ksp->its, norm_r);
829:     ksp->rnorm = norm_r;

831:     /*************************************************************************/
832:     /* Check for breakdown.                                                  */
833:     /*************************************************************************/

835:     beta = rz / rzm1;
836:     if (PetscAbsReal(beta) <= 0.0) {
837:       /***********************************************************************/
838:       /* Conjugate gradients has broken down.                                */
839:       /***********************************************************************/

841:       ksp->reason = KSP_DIVERGED_BREAKDOWN;
842:       PetscInfo1(ksp, "KSPSolve_GLTR: breakdown: beta=%g\n",(double) beta);
843:       break;
844:     }

846:     /*************************************************************************/
847:     /* Update p and the norms.                                               */
848:     /*************************************************************************/

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

853:     /*************************************************************************/
854:     /* Compute the new direction and update the iteration.                   */
855:     /*************************************************************************/

857:     KSP_MatMult(ksp, Qmat, p, z);  /* z = Q * p         */
858:     VecDot(p, z, &kappa);          /* kappa = p^T Q p   */
859:     ++ksp->its;

861:     /*************************************************************************/
862:     /* Update the Lanczos tridiagonal matrix.                                */
863:     /*************************************************************************/

865:     ++l_size;
866:     cg->offd[t_size] = PetscSqrtReal(beta) / PetscAbsReal(alpha);
867:     cg->diag[t_size] = kappa / rz + beta / alpha;
868:     ++t_size;
869:   }

871:   /***************************************************************************/
872:   /* We have the Lanczos basis, solve the tridiagonal trust-region problem   */
873:   /* to obtain the Lanczos direction.  We know that the solution lies on     */
874:   /* the boundary of the trust region.  We start by checking that the        */
875:   /* workspace allocated is large enough.                                    */
876:   /***************************************************************************/
877:   /* Note that the current version only computes the solution by using the   */
878:   /* preconditioned direction.  Need to think about how to do the            */
879:   /* unpreconditioned direction calculation.                                 */
880:   /***************************************************************************/

882:   if (t_size > cg->alloced) {
883:     if (cg->alloced) {
884:       PetscFree(cg->rwork);
885:       PetscFree(cg->iwork);
886:       cg->alloced += cg->init_alloc;
887:     } else {
888:       cg->alloced = cg->init_alloc;
889:     }

891:     while (t_size > cg->alloced) {
892:       cg->alloced += cg->init_alloc;
893:     }

895:     cg->alloced = PetscMin(cg->alloced, t_size);
896:     PetscMalloc2(10*cg->alloced, &cg->rwork,5*cg->alloced, &cg->iwork);
897:   }

899:   /***************************************************************************/
900:   /* Set up the required vectors.                                            */
901:   /***************************************************************************/

903:   t_soln = cg->rwork + 0*t_size;                        /* Solution          */
904:   t_diag = cg->rwork + 1*t_size;                        /* Diagonal of T     */
905:   t_offd = cg->rwork + 2*t_size;                        /* Off-diagonal of T */
906:   e_valu = cg->rwork + 3*t_size;                        /* Eigenvalues of T  */
907:   e_vect = cg->rwork + 4*t_size;                        /* Eigenvector of T  */
908:   e_rwrk = cg->rwork + 5*t_size;                        /* Eigen workspace   */

910:   e_iblk = cg->iwork + 0*t_size;                        /* Eigen blocks      */
911:   e_splt = cg->iwork + 1*t_size;                        /* Eigen splits      */
912:   e_iwrk = cg->iwork + 2*t_size;                        /* Eigen workspace   */

914:   /***************************************************************************/
915:   /* Compute the minimum eigenvalue of T.                                    */
916:   /***************************************************************************/

918:   vl = 0.0;
919:   vu = 0.0;
920:   il = 1;
921:   iu = 1;

923: #if defined(PETSC_MISSING_LAPACK_STEBZ)
924:   SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"STEBZ - Lapack routine is unavailable.");
925: #else
926:   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));

928:   if ((0 != info) || (1 != e_valus)) {
929:     /*************************************************************************/
930:     /* Calculation of the minimum eigenvalue failed.  Return the             */
931:     /* Steihaug-Toint direction.                                             */
932:     /*************************************************************************/

934:     PetscInfo(ksp, "KSPSolve_GLTR: failed to compute eigenvalue.\n");
935:     ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
936:     return(0);
937:   }

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

941:   /***************************************************************************/
942:   /* Compute the initial value of lambda to make (T + lamba I) positive      */
943:   /* definite.                                                               */
944:   /***************************************************************************/

946:   pert = cg->init_pert;
947:   if (e_valu[0] < 0.0) cg->lambda = pert - e_valu[0];
948: #endif

950:   while (1) {
951:     for (i = 0; i < t_size; ++i) {
952:       t_diag[i] = cg->diag[i] + cg->lambda;
953:       t_offd[i] = cg->offd[i];
954:     }

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

961:     if (0 == info) break;

963:     pert      += pert;
964:     cg->lambda = cg->lambda * (1.0 + pert) + pert;
965: #endif
966:   }

968:   /***************************************************************************/
969:   /* Compute the initial step and its norm.                                  */
970:   /***************************************************************************/

972:   nrhs = 1;
973:   nldb = t_size;

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

978: #if defined(PETSC_MISSING_LAPACK_PTTRS)
979:   SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"PTTRS - Lapack routine is unavailable.");
980: #else
981:   PetscStackCallBLAS("LAPACKpttrs",LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, t_soln, &nldb, &info));
982: #endif

984:   if (0 != info) {
985:     /*************************************************************************/
986:     /* Calculation of the initial step failed; return the Steihaug-Toint     */
987:     /* direction.                                                            */
988:     /*************************************************************************/

990:     PetscInfo(ksp, "KSPSolve_GLTR: failed to compute step.\n");
991:     ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
992:     return(0);
993:   }

995:   norm_t = 0.;
996:   for (i = 0; i < t_size; ++i) norm_t += t_soln[i] * t_soln[i];
997:   norm_t = PetscSqrtReal(norm_t);

999:   /***************************************************************************/
1000:   /* Determine the case we are in.                                           */
1001:   /***************************************************************************/

1003:   if (norm_t <= cg->radius) {
1004:     /*************************************************************************/
1005:     /* The step is within the trust region; check if we are in the hard case */
1006:     /* and need to move to the boundary by following a direction of negative */
1007:     /* curvature.                                                            */
1008:     /*************************************************************************/

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

1016: #if defined(PETSC_MISSING_LAPACK_STEIN)
1017:       SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"STEIN - Lapack routine is unavailable.");
1018: #else
1019:       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));
1020: #endif

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

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

1033:       coef1 = 0.0;
1034:       coef2 = 0.0;
1035:       coef3 = -cg->radius * cg->radius;
1036:       for (i = 0; i < t_size; ++i) {
1037:         coef1 += e_vect[i] * e_vect[i];
1038:         coef2 += e_vect[i] * t_soln[i];
1039:         coef3 += t_soln[i] * t_soln[i];
1040:       }

1042:       coef3 = PetscSqrtReal(coef2 * coef2 - coef1 * coef3);
1043:       root1 = (-coef2 + coef3) / coef1;
1044:       root2 = (-coef2 - coef3) / coef1;

1046:       /***********************************************************************/
1047:       /* Compute objective value for (t_soln + root1 * e_vect)               */
1048:       /***********************************************************************/

1050:       for (i = 0; i < t_size; ++i) {
1051:         e_rwrk[i] = t_soln[i] + root1 * e_vect[i];
1052:       }

1054:       obj1 = e_rwrk[0]*(0.5*(cg->diag[0]*e_rwrk[0]+
1055:                              cg->offd[1]*e_rwrk[1])+cg->norm_r[0]);
1056:       for (i = 1; i < t_size - 1; ++i) {
1057:         obj1 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
1058:                                cg->diag[i]*e_rwrk[i]+
1059:                                cg->offd[i+1]*e_rwrk[i+1]);
1060:       }
1061:       obj1 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
1062:                              cg->diag[i]*e_rwrk[i]);

1064:       /***********************************************************************/
1065:       /* Compute objective value for (t_soln + root2 * e_vect)               */
1066:       /***********************************************************************/

1068:       for (i = 0; i < t_size; ++i) {
1069:         e_rwrk[i] = t_soln[i] + root2 * e_vect[i];
1070:       }

1072:       obj2 = e_rwrk[0]*(0.5*(cg->diag[0]*e_rwrk[0]+
1073:                              cg->offd[1]*e_rwrk[1])+cg->norm_r[0]);
1074:       for (i = 1; i < t_size - 1; ++i) {
1075:         obj2 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
1076:                                cg->diag[i]*e_rwrk[i]+
1077:                                cg->offd[i+1]*e_rwrk[i+1]);
1078:       }
1079:       obj2 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
1080:                              cg->diag[i]*e_rwrk[i]);

1082:       /***********************************************************************/
1083:       /* Choose the point with the best objective function value.            */
1084:       /***********************************************************************/

1086:       if (obj1 <= obj2) {
1087:         for (i = 0; i < t_size; ++i) {
1088:           t_soln[i] += root1 * e_vect[i];
1089:         }
1090:       }
1091:       else {
1092:         for (i = 0; i < t_size; ++i) {
1093:           t_soln[i] += root2 * e_vect[i];
1094:         }
1095:       }
1096:     } else {
1097:       /***********************************************************************/
1098:       /* The matrix is positive definite or there was no room to move; the   */
1099:       /* solution is already contained in t_soln.                            */
1100:       /***********************************************************************/
1101:     }
1102:   } else {
1103:     /*************************************************************************/
1104:     /* The step is outside the trust-region.  Compute the correct value for  */
1105:     /* lambda by performing Newton's method.                                 */
1106:     /*************************************************************************/

1108:     for (i = 0; i < max_newton_its; ++i) {
1109:       /***********************************************************************/
1110:       /* Check for convergence.                                              */
1111:       /***********************************************************************/

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

1115:       /***********************************************************************/
1116:       /* Compute the update.                                                 */
1117:       /***********************************************************************/

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

1121: #if defined(PETSC_MISSING_LAPACK_PTTRS)
1122:       SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"PTTRS - Lapack routine is unavailable.");
1123: #else
1124:       PetscStackCallBLAS("LAPACKpttrs",LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, e_rwrk, &nldb, &info));
1125: #endif

1127:       if (0 != info) {
1128:         /*********************************************************************/
1129:         /* Calculation of the step failed; return the Steihaug-Toint         */
1130:         /* direction.                                                        */
1131:         /*********************************************************************/

1133:         PetscInfo(ksp, "KSPSolve_GLTR: failed to compute step.\n");
1134:         ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1135:         return(0);
1136:       }

1138:       /***********************************************************************/
1139:       /* Modify lambda.                                                      */
1140:       /***********************************************************************/

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

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

1147:       /***********************************************************************/
1148:       /* Factor T + lambda I                                                 */
1149:       /***********************************************************************/

1151:       for (j = 0; j < t_size; ++j) {
1152:         t_diag[j] = cg->diag[j] + cg->lambda;
1153:         t_offd[j] = cg->offd[j];
1154:       }

1156: #if defined(PETSC_MISSING_LAPACK_PTTRF)
1157:       SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"PTTRF - Lapack routine is unavailable.");
1158: #else
1159:       PetscStackCallBLAS("LAPACKpttrf",LAPACKpttrf_(&t_size, t_diag, t_offd + 1, &info));
1160: #endif

1162:       if (0 != info) {
1163:         /*********************************************************************/
1164:         /* Calculation of factorization failed; return the Steihaug-Toint    */
1165:         /* direction.                                                        */
1166:         /*********************************************************************/

1168:         PetscInfo(ksp, "KSPSolve_GLTR: factorization failed.\n");
1169:         ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1170:         return(0);
1171:       }

1173:       /***********************************************************************/
1174:       /* Compute the new step and its norm.                                  */
1175:       /***********************************************************************/

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

1180: #if defined(PETSC_MISSING_LAPACK_PTTRS)
1181:       SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"PTTRS - Lapack routine is unavailable.");
1182: #else
1183:       PetscStackCallBLAS("LAPACKpttrs",LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, t_soln, &nldb, &info));
1184: #endif

1186:       if (0 != info) {
1187:         /*********************************************************************/
1188:         /* Calculation of the step failed; return the Steihaug-Toint         */
1189:         /* direction.                                                        */
1190:         /*********************************************************************/

1192:         PetscInfo(ksp, "KSPSolve_GLTR: failed to compute step.\n");
1193:         ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1194:         return(0);
1195:       }

1197:       norm_t = 0.;
1198:       for (j = 0; j < t_size; ++j) norm_t += t_soln[j] * t_soln[j];
1199:       norm_t = PetscSqrtReal(norm_t);
1200:     }

1202:     /*************************************************************************/
1203:     /* Check for convergence.                                                */
1204:     /*************************************************************************/

1206:     if (PetscAbsReal(norm_t - cg->radius) > cg->newton_tol * cg->radius) {
1207:       /***********************************************************************/
1208:       /* Newton method failed to converge in iteration limit.                */
1209:       /***********************************************************************/

1211:       PetscInfo(ksp, "KSPSolve_GLTR: failed to converge.\n");
1212:       ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1213:       return(0);
1214:     }
1215:   }

1217:   /***************************************************************************/
1218:   /* Recover the norm of the direction and objective function value.         */
1219:   /***************************************************************************/

1221:   cg->norm_d = norm_t;

1223:   cg->o_fcn = t_soln[0]*(0.5*(cg->diag[0]*t_soln[0]+cg->offd[1]*t_soln[1])+cg->norm_r[0]);
1224:   for (i = 1; i < t_size - 1; ++i) {
1225:     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]);
1226:   }
1227:   cg->o_fcn += 0.5*t_soln[i]*(cg->offd[i]*t_soln[i-1]+cg->diag[i]*t_soln[i]);

1229:   /***************************************************************************/
1230:   /* Recover the direction.                                                  */
1231:   /***************************************************************************/

1233:   sigma = -1;

1235:   /***************************************************************************/
1236:   /* Start conjugate gradient method from the beginning                      */
1237:   /***************************************************************************/

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

1242:   /***************************************************************************/
1243:   /* Accumulate Q * s                                                        */
1244:   /***************************************************************************/

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

1249:   /***************************************************************************/
1250:   /* Compute the first direction.                                            */
1251:   /***************************************************************************/

1253:   VecCopy(z, p);                   /* p = z             */
1254:   KSP_MatMult(ksp, Qmat, p, z);    /* z = Q * p         */
1255:   ++ksp->its;

1257:   for (i = 0; i < l_size - 1; ++i) {
1258:     /*************************************************************************/
1259:     /* Update the residual and direction.                                    */
1260:     /*************************************************************************/

1262:     alpha = cg->alpha[i];
1263:     if (alpha >= 0.0) sigma = -sigma;

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

1268:     /*************************************************************************/
1269:     /* Accumulate Q * s                                                      */
1270:     /*************************************************************************/

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

1274:     /*************************************************************************/
1275:     /* Update p.                                                             */
1276:     /*************************************************************************/

1278:     beta = cg->beta[i];
1279:     VecAYPX(p, beta, z);          /* p = z + beta p    */
1280:     KSP_MatMult(ksp, Qmat, p, z);  /* z = Q * p         */
1281:     ++ksp->its;
1282:   }

1284:   /***************************************************************************/
1285:   /* Update the residual and direction.                                      */
1286:   /***************************************************************************/

1288:   alpha = cg->alpha[i];
1289:   if (alpha >= 0.0) sigma = -sigma;

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

1294:   /***************************************************************************/
1295:   /* Accumulate Q * s                                                        */
1296:   /***************************************************************************/

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

1300:   /***************************************************************************/
1301:   /* Set the termination reason.                                             */
1302:   /***************************************************************************/

1304:   ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1305:   return(0);
1306: #endif
1307: }

1311: PetscErrorCode KSPSetUp_GLTR(KSP ksp)
1312: {
1313:   KSP_GLTR       *cg = (KSP_GLTR*)ksp->data;
1314:   PetscInt       max_its;

1318:   /***************************************************************************/
1319:   /* Determine the total maximum number of iterations.                       */
1320:   /***************************************************************************/

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

1324:   /***************************************************************************/
1325:   /* Set work vectors needed by conjugate gradient method and allocate       */
1326:   /* workspace for Lanczos matrix.                                           */
1327:   /***************************************************************************/

1329:   KSPSetWorkVecs(ksp, 3);

1331:   PetscCalloc5(max_its,&cg->diag,max_its,&cg->offd,max_its,&cg->alpha,max_its,&cg->beta,max_its,&cg->norm_r);
1332:   PetscLogObjectMemory((PetscObject)ksp, 5*max_its*sizeof(PetscReal));
1333:   return(0);
1334: }

1338: PetscErrorCode KSPDestroy_GLTR(KSP ksp)
1339: {
1340:   KSP_GLTR       *cg = (KSP_GLTR*)ksp->data;

1344:   /***************************************************************************/
1345:   /* Free memory allocated for the data.                                     */
1346:   /***************************************************************************/

1348:   PetscFree5(cg->diag,cg->offd,cg->alpha,cg->beta,cg->norm_r);
1349:   if (cg->alloced) {
1350:     PetscFree2(cg->rwork,cg->iwork);
1351:   }

1353:   /***************************************************************************/
1354:   /* Clear composed functions                                                */
1355:   /***************************************************************************/

1357:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRSetRadius_C",NULL);
1358:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetNormD_C",NULL);
1359:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetObjFcn_C",NULL);
1360:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetMinEig_C",NULL);
1361:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetLambda_C",NULL);

1363:   /***************************************************************************/
1364:   /* Destroy KSP object.                                                     */
1365:   /***************************************************************************/
1366:   KSPDestroyDefault(ksp);
1367:   return(0);
1368: }

1372: static PetscErrorCode  KSPGLTRSetRadius_GLTR(KSP ksp, PetscReal radius)
1373: {
1374:   KSP_GLTR *cg = (KSP_GLTR*)ksp->data;

1377:   cg->radius = radius;
1378:   return(0);
1379: }

1383: static PetscErrorCode  KSPGLTRGetNormD_GLTR(KSP ksp, PetscReal *norm_d)
1384: {
1385:   KSP_GLTR *cg = (KSP_GLTR*)ksp->data;

1388:   *norm_d = cg->norm_d;
1389:   return(0);
1390: }

1394: static PetscErrorCode  KSPGLTRGetObjFcn_GLTR(KSP ksp, PetscReal *o_fcn)
1395: {
1396:   KSP_GLTR *cg = (KSP_GLTR*)ksp->data;

1399:   *o_fcn = cg->o_fcn;
1400:   return(0);
1401: }

1405: static PetscErrorCode  KSPGLTRGetMinEig_GLTR(KSP ksp, PetscReal *e_min)
1406: {
1407:   KSP_GLTR *cg = (KSP_GLTR*)ksp->data;

1410:   *e_min = cg->e_min;
1411:   return(0);
1412: }

1416: static PetscErrorCode  KSPGLTRGetLambda_GLTR(KSP ksp, PetscReal *lambda)
1417: {
1418:   KSP_GLTR *cg = (KSP_GLTR*)ksp->data;

1421:   *lambda = cg->lambda;
1422:   return(0);
1423: }

1427: PetscErrorCode KSPSetFromOptions_GLTR(PetscOptions *PetscOptionsObject,KSP ksp)
1428: {
1430:   KSP_GLTR       *cg = (KSP_GLTR*)ksp->data;

1433:   PetscOptionsHead(PetscOptionsObject,"KSP GLTR options");

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

1437:   PetscOptionsReal("-ksp_gltr_init_pert", "Initial perturbation", "", cg->init_pert, &cg->init_pert, NULL);
1438:   PetscOptionsReal("-ksp_gltr_eigen_tol", "Eigenvalue tolerance", "", cg->eigen_tol, &cg->eigen_tol, NULL);
1439:   PetscOptionsReal("-ksp_gltr_newton_tol", "Newton tolerance", "", cg->newton_tol, &cg->newton_tol, NULL);

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

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

1446:   PetscOptionsTail();
1447:   return(0);
1448: }

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

1455:    Options Database Keys:
1456: .      -ksp_gltr_radius <r> - Trust Region Radius

1458:    Notes: This is rarely used directly

1460:   Use preconditioned conjugate gradient to compute
1461:   an approximate minimizer of the quadratic function

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

1465:    subject to the trust region constraint

1467:             || s || <= delta,

1469:    where

1471:      delta is the trust region radius,
1472:      g is the gradient vector,
1473:      H is the Hessian approximation,
1474:      M is the positive definite preconditioner matrix.

1476:    KSPConvergedReason may be
1477: $  KSP_CONVERGED_CG_NEG_CURVE if convergence is reached along a negative curvature direction,
1478: $  KSP_CONVERGED_CG_CONSTRAINED if convergence is reached along a constrained step,
1479: $  other KSP converged/diverged reasons

1481:   Notes:
1482:   The preconditioner supplied should be symmetric and positive definite.

1484:    Level: developer

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

1491: PETSC_EXTERN PetscErrorCode KSPCreate_GLTR(KSP ksp)
1492: {
1494:   KSP_GLTR       *cg;

1497:   PetscNewLog(ksp,&cg);
1498:   cg->radius = 0.0;
1499:   cg->dtype  = GLTR_UNPRECONDITIONED_DIRECTION;

1501:   cg->init_pert  = 1.0e-8;
1502:   cg->eigen_tol  = 1.0e-10;
1503:   cg->newton_tol = 1.0e-6;

1505:   cg->alloced    = 0;
1506:   cg->init_alloc = 1024;

1508:   cg->max_lanczos_its = 20;
1509:   cg->max_newton_its  = 10;

1511:   ksp->data = (void*) cg;
1512:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,3);
1513:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);
1514:   KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,2);

1516:   /***************************************************************************/
1517:   /* Sets the functions that are associated with this data structure         */
1518:   /* (in C++ this is the same as defining virtual functions).                */
1519:   /***************************************************************************/

1521:   ksp->ops->setup          = KSPSetUp_GLTR;
1522:   ksp->ops->solve          = KSPSolve_GLTR;
1523:   ksp->ops->destroy        = KSPDestroy_GLTR;
1524:   ksp->ops->setfromoptions = KSPSetFromOptions_GLTR;
1525:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
1526:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
1527:   ksp->ops->view           = 0;

1529:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRSetRadius_C",KSPGLTRSetRadius_GLTR);
1530:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetNormD_C", KSPGLTRGetNormD_GLTR);
1531:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetObjFcn_C",KSPGLTRGetObjFcn_GLTR);
1532:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetMinEig_C",KSPGLTRGetMinEig_GLTR);
1533:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetLambda_C",KSPGLTRGetLambda_GLTR);
1534:   return(0);
1535: }