Actual source code: gltr.c

petsc-3.3-p7 2013-05-11
  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(((PetscObject)ksp)->comm,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: #ifdef PETSC_USE_COMPLEX
147:   SETERRQ(((PetscObject)ksp)->comm,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:   MatStructure   pflag;
155:   Mat            Qmat, Mmat;
156:   Vec            r, z, p, d;
157:   PC             pc;

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

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

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

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

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

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

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

195:   PCGetOperators(pc, &Qmat, &Mmat, &pflag);

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

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

207:   cg->o_fcn = 0.0;

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

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

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

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

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

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

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

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

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

263:       VecAXPY(d, alpha, r);        /* d = d + alpha r   */
264: 
265:       /***********************************************************************/
266:       /* Compute objective function.                                         */
267:       /***********************************************************************/

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

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

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

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

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

301:       /***********************************************************************/
302:       /* Compute objective function.                                         */
303:       /***********************************************************************/

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

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

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

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

328:   case KSP_NORM_UNPRECONDITIONED:
329:     norm_r = PetscSqrtReal(rr);                                        /* norm_r = |r|      */
330:     break;

332:   case KSP_NORM_NATURAL:
333:     norm_r = cg->norm_r[0];                                /* norm_r = |r|_M    */
334:     break;

336:   default:
337:     norm_r = 0.0;
338:     break;
339:   }

341:   KSPLogResidualHistory(ksp, norm_r);
342:   KSPMonitor(ksp, ksp->its, norm_r);
343:   ksp->rnorm = norm_r;

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

347:   /***************************************************************************/
348:   /* Compute the first direction and update the iteration.                   */
349:   /***************************************************************************/

351:   VecCopy(z, p);                        /* p = z             */
352:   KSP_MatMult(ksp, Qmat, p, z);        /* z = Q * p         */
353:   ++ksp->its;

355:   /***************************************************************************/
356:   /* Check the matrix for numerical problems.                                */
357:   /***************************************************************************/

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

367:     ksp->reason = KSP_DIVERGED_NAN;
368:     PetscInfo1(ksp, "KSPSolve_GLTR: bad matrix: kappa=%g\n", kappa);

370:     if (cg->radius) {
371:       if (r2 >= rr) {
372:         alpha = 1.0;
373:         cg->norm_d = PetscSqrtReal(rr);
374:       }
375:       else {
376:         alpha = PetscSqrtReal(r2 / rr);
377:         cg->norm_d = cg->radius;
378:       }

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

382:       /***********************************************************************/
383:       /* Compute objective function.                                         */
384:       /***********************************************************************/

386:       KSP_MatMult(ksp, Qmat, d, z);
387:       VecAYPX(z, -0.5, ksp->vec_rhs);
388:       VecDot(d, z, &cg->o_fcn);
389:       cg->o_fcn = -cg->o_fcn;
390:       ++ksp->its;
391:     }
392:     return(0);
393:   }

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

402:   dMp = 0.0;
403:   norm_d = 0.0;
404:   switch(cg->dtype) {
405:   case GLTR_PRECONDITIONED_DIRECTION:
406:     norm_p = rz;
407:     break;

409:   default:
410:     VecDot(p, p, &norm_p);
411:     break;
412:   }

414:   cg->diag[t_size] = kappa / rz;
415:   cg->offd[t_size] = 0.0;
416:   ++t_size;

418:   piv = 1.0;

420:   /***************************************************************************/
421:   /* Check for breakdown of the conjugate gradient method; this occurs when  */
422:   /* kappa is zero.                                                          */
423:   /***************************************************************************/

425:   if (fabs(kappa) <= 0.0) {
426:     /*************************************************************************/
427:     /* The curvature is zero.  In this case, we must stop and use follow     */
428:     /* the direction of negative curvature since the Lanczos matrix is zero. */
429:     /*************************************************************************/

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

434:     if (cg->radius && norm_p > 0.0) {
435:       /***********************************************************************/
436:       /* Follow direction of negative curvature to the boundary of the       */
437:       /* trust region.                                                       */
438:       /***********************************************************************/

440:       step = PetscSqrtReal(r2 / norm_p);
441:       cg->norm_d = cg->radius;

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

445:       /***********************************************************************/
446:       /* Update objective function.                                          */
447:       /***********************************************************************/

449:       cg->o_fcn += step * (0.5 * step * kappa - rz);
450:     }
451:     else if (cg->radius) {
452:       /***********************************************************************/
453:       /* The norm of the preconditioned direction is zero; use the gradient  */
454:       /* step.                                                               */
455:       /***********************************************************************/

457:       if (r2 >= rr) {
458:         alpha = 1.0;
459:         cg->norm_d = PetscSqrtReal(rr);
460:       }
461:       else {
462:         alpha = PetscSqrtReal(r2 / rr);
463:         cg->norm_d = cg->radius;
464:       }

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

468:       /***********************************************************************/
469:       /* Compute objective function.                                         */
470:       /***********************************************************************/

472:       KSP_MatMult(ksp, Qmat, d, z);
473:       VecAYPX(z, -0.5, ksp->vec_rhs);
474:       VecDot(d, z, &cg->o_fcn);
475:       cg->o_fcn = -cg->o_fcn;
476:       ++ksp->its;
477:     }
478:     return(0);
479:   }

481:   /***************************************************************************/
482:   /* Begin the first part of the GLTR algorithm which runs the conjugate     */
483:   /* gradient method until either the problem is solved, we encounter the    */
484:   /* boundary of the trust region, or the conjugate gradient method breaks   */
485:   /* down.                                                                   */
486:   /***************************************************************************/

488:   while(1) {
489:     /*************************************************************************/
490:     /* Know that kappa is nonzero, because we have not broken down, so we    */
491:     /* can compute the steplength.                                           */
492:     /*************************************************************************/

494:     alpha = rz / kappa;
495:     cg->alpha[l_size] = alpha;

497:     /*************************************************************************/
498:     /* Compute the diagonal value of the Cholesky factorization of the       */
499:     /* Lanczos matrix and check to see if the Lanczos matrix is indefinite.  */
500:     /* This indicates a direction of negative curvature.                     */
501:     /*************************************************************************/

503:     piv = cg->diag[l_size] - cg->offd[l_size]*cg->offd[l_size] / piv;
504:     if (piv <= 0.0) {
505:       /***********************************************************************/
506:       /* In this case, the matrix is indefinite and we have encountered      */
507:       /* a direction of negative curvature.  Follow the direction to the     */
508:       /* boundary of the trust region.                                       */
509:       /***********************************************************************/

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

514:       if (cg->radius && norm_p > 0.0) {
515:         /*********************************************************************/
516:         /* Follow direction of negative curvature to boundary.               */
517:         /*********************************************************************/

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

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

524:         /*********************************************************************/
525:         /* Update objective function.                                        */
526:         /*********************************************************************/

528:         cg->o_fcn += step * (0.5 * step * kappa - rz);
529:       }
530:       else if (cg->radius) {
531:         /*********************************************************************/
532:         /* The norm of the direction is zero; there is nothing to follow.    */
533:         /*********************************************************************/
534:       }
535:       break;
536:     }

538:     /*************************************************************************/
539:     /* Compute the steplength and check for intersection with the trust      */
540:     /* region.                                                               */
541:     /*************************************************************************/

543:     norm_dp1 = norm_d + alpha*(2.0*dMp + alpha*norm_p);
544:     if (cg->radius && norm_dp1 >= r2) {
545:       /***********************************************************************/
546:       /* In this case, the matrix is positive definite as far as we know.    */
547:       /* However, the full step goes beyond the trust region.                */
548:       /***********************************************************************/

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

553:       if (norm_p > 0.0) {
554:         /*********************************************************************/
555:         /* Follow the direction to the boundary of the trust region.         */
556:         /*********************************************************************/

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

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

563:         /*********************************************************************/
564:         /* Update objective function.                                        */
565:         /*********************************************************************/

567:         cg->o_fcn += step * (0.5 * step * kappa - rz);
568:       }
569:       else {
570:         /*********************************************************************/
571:         /* The norm of the direction is zero; there is nothing to follow.    */
572:         /*********************************************************************/
573:       }
574:       break;
575:     }

577:     /*************************************************************************/
578:     /* Now we can update the direction and residual.                         */
579:     /*************************************************************************/

581:     VecAXPY(d, alpha, p);                /* d = d + alpha p   */
582:     VecAXPY(r, -alpha, z);                        /* r = r - alpha Q p */
583:     KSP_PCApply(ksp, r, z);        /* z = inv(M) r      */

585:     switch(cg->dtype) {
586:     case GLTR_PRECONDITIONED_DIRECTION:
587:       norm_d = norm_dp1;
588:       break;

590:     default:
591:       VecDot(d, d, &norm_d);
592:       break;
593:     }
594:     cg->norm_d = PetscSqrtReal(norm_d);

596:     /*************************************************************************/
597:     /* Update objective function.                                            */
598:     /*************************************************************************/

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

602:     /*************************************************************************/
603:     /* Check that the preconditioner appears positive semidefinite.          */
604:     /*************************************************************************/

606:     rzm1 = rz;
607:     VecDot(r, z, &rz);                /* rz = r^T z        */
608:     if (rz < 0.0) {
609:       /***********************************************************************/
610:       /* The preconditioner is indefinite.                                   */
611:       /***********************************************************************/

613:       ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
614:       PetscInfo1(ksp, "KSPSolve_GLTR: cg indefinite preconditioner: rz=%g\n", rz);
615:       break;
616:     }

618:     /*************************************************************************/
619:     /* As far as we know, the preconditioner is positive semidefinite.       */
620:     /* Compute the residual and check for convergence.                       */
621:     /*************************************************************************/

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

625:     switch(ksp->normtype) {
626:     case KSP_NORM_PRECONDITIONED:
627:       VecNorm(z, NORM_2, &norm_r);/* norm_r = |z|      */
628:       break;

630:     case KSP_NORM_UNPRECONDITIONED:
631:       VecNorm(r, NORM_2, &norm_r);/* norm_r = |r|      */
632:       break;

634:     case KSP_NORM_NATURAL:
635:       norm_r = cg->norm_r[l_size+1];                        /* norm_r = |r|_M    */
636:       break;

638:     default:
639:       norm_r = 0.0;
640:       break;
641:     }

643:     KSPLogResidualHistory(ksp, norm_r);
644:     KSPMonitor(ksp, ksp->its, norm_r);
645:     ksp->rnorm = norm_r;
646: 
647:     (*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP);
648:     if (ksp->reason) {
649:       /***********************************************************************/
650:       /* The method has converged.                                           */
651:       /***********************************************************************/

653:       PetscInfo2(ksp, "KSPSolve_GLTR: cg truncated step: rnorm=%g, radius=%g\n", norm_r, cg->radius);
654:       break;
655:     }

657:     /*************************************************************************/
658:     /* We have not converged yet.  Check for breakdown.                      */
659:     /*************************************************************************/

661:     beta = rz / rzm1;
662:     if (fabs(beta) <= 0.0) {
663:       /***********************************************************************/
664:       /* Conjugate gradients has broken down.                                */
665:       /***********************************************************************/

667:       ksp->reason = KSP_DIVERGED_BREAKDOWN;
668:       PetscInfo1(ksp, "KSPSolve_GLTR: breakdown: beta=%g\n", beta);
669:       break;
670:     }

672:     /*************************************************************************/
673:     /* Check iteration limit.                                                */
674:     /*************************************************************************/

676:     if (ksp->its >= max_cg_its) {
677:       ksp->reason = KSP_DIVERGED_ITS;
678:       PetscInfo1(ksp, "KSPSolve_GLTR: iterlim: its=%d\n", ksp->its);
679:       break;
680:     }

682:     /*************************************************************************/
683:     /* Update p and the norms.                                               */
684:     /*************************************************************************/

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

689:     switch(cg->dtype) {
690:     case GLTR_PRECONDITIONED_DIRECTION:
691:       dMp = beta*(dMp + alpha*norm_p);
692:       norm_p = beta*(rzm1 + beta*norm_p);
693:       break;

695:     default:
696:       VecDot(d, p, &dMp);
697:       VecDot(p, p, &norm_p);
698:       break;
699:     }

701:     /*************************************************************************/
702:     /* Compute the new direction and update the iteration.                   */
703:     /*************************************************************************/

705:     KSP_MatMult(ksp, Qmat, p, z);        /* z = Q * p         */
706:     VecDot(p, z, &kappa);                /* kappa = p^T Q p   */
707:     ++ksp->its;

709:     /*************************************************************************/
710:     /* Update the the Lanczos tridiagonal matrix.                            */
711:     /*************************************************************************/

713:     ++l_size;
714:     cg->offd[t_size] = PetscSqrtReal(beta) / fabs(alpha);
715:     cg->diag[t_size] = kappa / rz + beta / alpha;
716:     ++t_size;

718:     /*************************************************************************/
719:     /* Check for breakdown of the conjugate gradient method; this occurs     */
720:     /* when kappa is zero.                                                   */
721:     /*************************************************************************/

723:     if (fabs(kappa) <= 0.0) {
724:       /***********************************************************************/
725:       /* The method breaks down; move along the direction as if the matrix   */
726:       /* were indefinite.                                                    */
727:       /***********************************************************************/

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

732:       if (cg->radius && norm_p > 0.0) {
733:         /*********************************************************************/
734:         /* Follow direction to boundary.                                     */
735:         /*********************************************************************/

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

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

742:         /*********************************************************************/
743:         /* Update objective function.                                        */
744:         /*********************************************************************/

746:         cg->o_fcn += step * (0.5 * step * kappa - rz);
747:       }
748:       else if (cg->radius) {
749:         /*********************************************************************/
750:         /* The norm of the direction is zero; there is nothing to follow.    */
751:         /*********************************************************************/
752:       }
753:       break;
754:     }
755:   }

757:   /***************************************************************************/
758:   /* Check to see if we need to continue with the Lanczos method.            */
759:   /***************************************************************************/

761:   if (!cg->radius) {
762:     /*************************************************************************/
763:     /* There is no radius.  Therefore, we cannot move along the boundary.    */
764:     /*************************************************************************/

766:     return(0);
767:   }

769:   if (KSP_CONVERGED_CG_NEG_CURVE != ksp->reason) {
770:     /*************************************************************************/
771:     /* The method either converged to an interior point, hit the boundary of */
772:     /* the trust-region without encountering a direction of negative         */
773:     /* curvature or the iteration limit was reached.                         */
774:     /*************************************************************************/

776:     return(0);
777:   }

779:   /***************************************************************************/
780:   /* Switch to contructing the Lanczos basis by way of the conjugate         */
781:   /* directions.                                                             */
782:   /***************************************************************************/

784:   for (i = 0; i < max_lanczos_its; ++i) {
785:     /*************************************************************************/
786:     /* Check for breakdown of the conjugate gradient method; this occurs     */
787:     /* when kappa is zero.                                                   */
788:     /*************************************************************************/

790:     if (fabs(kappa) <= 0.0) {
791:       ksp->reason = KSP_DIVERGED_BREAKDOWN;
792:       PetscInfo1(ksp, "KSPSolve_GLTR: lanczos breakdown: kappa=%g\n", kappa);
793:       break;
794:     }

796:     /*************************************************************************/
797:     /* Update the direction and residual.                                    */
798:     /*************************************************************************/
799: 
800:     alpha = rz / kappa;
801:     cg->alpha[l_size] = alpha;

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

806:     /*************************************************************************/
807:     /* Check that the preconditioner appears positive semidefinite.          */
808:     /*************************************************************************/

810:     rzm1 = rz;
811:     VecDot(r, z, &rz);                /* rz = r^T z        */
812:     if (rz < 0.0) {
813:       /***********************************************************************/
814:       /* The preconditioner is indefinite.                                   */
815:       /***********************************************************************/

817:       ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
818:       PetscInfo1(ksp, "KSPSolve_GLTR: lanczos indefinite preconditioner: rz=%g\n", rz);
819:       break;
820:     }

822:     /*************************************************************************/
823:     /* As far as we know, the preconditioner is positive definite.  Compute  */
824:     /* the residual.  Do NOT check for convergence.                          */
825:     /*************************************************************************/

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

829:     switch(ksp->normtype) {
830:     case KSP_NORM_PRECONDITIONED:
831:       VecNorm(z, NORM_2, &norm_r);/* norm_r = |z|      */
832:       break;

834:     case KSP_NORM_UNPRECONDITIONED:
835:       VecNorm(r, NORM_2, &norm_r);/* norm_r = |r|      */
836:       break;

838:     case KSP_NORM_NATURAL:
839:       norm_r = cg->norm_r[l_size+1];                        /* norm_r = |r|_M    */
840:       break;

842:     default:
843:       norm_r = 0.0;
844:       break;
845:     }

847:     KSPLogResidualHistory(ksp, norm_r);
848:     KSPMonitor(ksp, ksp->its, norm_r);
849:     ksp->rnorm = norm_r;
850: 
851:     /*************************************************************************/
852:     /* Check for breakdown.                                                  */
853:     /*************************************************************************/

855:     beta = rz / rzm1;
856:     if (fabs(beta) <= 0.0) {
857:       /***********************************************************************/
858:       /* Conjugate gradients has broken down.                                */
859:       /***********************************************************************/

861:       ksp->reason = KSP_DIVERGED_BREAKDOWN;
862:       PetscInfo1(ksp, "KSPSolve_GLTR: breakdown: beta=%g\n", beta);
863:       break;
864:     }

866:     /*************************************************************************/
867:     /* Update p and the norms.                                               */
868:     /*************************************************************************/

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

873:     /*************************************************************************/
874:     /* Compute the new direction and update the iteration.                   */
875:     /*************************************************************************/

877:     KSP_MatMult(ksp, Qmat, p, z);        /* z = Q * p         */
878:     VecDot(p, z, &kappa);                /* kappa = p^T Q p   */
879:     ++ksp->its;

881:     /*************************************************************************/
882:     /* Update the Lanczos tridiagonal matrix.                                */
883:     /*************************************************************************/

885:     ++l_size;
886:     cg->offd[t_size] = PetscSqrtReal(beta) / fabs(alpha);
887:     cg->diag[t_size] = kappa / rz + beta / alpha;
888:     ++t_size;
889:   }

891:   /***************************************************************************/
892:   /* We have the Lanczos basis, solve the tridiagonal trust-region problem   */
893:   /* to obtain the Lanczos direction.  We know that the solution lies on     */
894:   /* the boundary of the trust region.  We start by checking that the        */
895:   /* workspace allocated is large enough.                                    */
896:   /***************************************************************************/
897:   /* Note that the current version only computes the solution by using the   */
898:   /* preconditioned direction.  Need to think about how to do the            */
899:   /* unpreconditioned direction calculation.                                 */
900:   /***************************************************************************/

902:   if (t_size > cg->alloced) {
903:     if (cg->alloced) {
904:       PetscFree(cg->rwork);
905:       PetscFree(cg->iwork);
906:       cg->alloced += cg->init_alloc;
907:     } else {
908:       cg->alloced = cg->init_alloc;
909:     }

911:     while (t_size > cg->alloced) {
912:       cg->alloced += cg->init_alloc;
913:     }
914: 
915:     cg->alloced = PetscMin(cg->alloced, t_size);
916:     PetscMalloc2(10*cg->alloced,PetscReal, &cg->rwork,5*cg->alloced,PetscBLASInt, &cg->iwork);
917:   }

919:   /***************************************************************************/
920:   /* Set up the required vectors.                                            */
921:   /***************************************************************************/

923:   t_soln = cg->rwork + 0*t_size;                        /* Solution          */
924:   t_diag = cg->rwork + 1*t_size;                        /* Diagonal of T     */
925:   t_offd = cg->rwork + 2*t_size;                        /* Off-diagonal of T */
926:   e_valu = cg->rwork + 3*t_size;                        /* Eigenvalues of T  */
927:   e_vect = cg->rwork + 4*t_size;                        /* Eigenvector of T  */
928:   e_rwrk = cg->rwork + 5*t_size;                        /* Eigen workspace   */
929: 
930:   e_iblk = cg->iwork + 0*t_size;                        /* Eigen blocks      */
931:   e_splt = cg->iwork + 1*t_size;                        /* Eigen splits      */
932:   e_iwrk = cg->iwork + 2*t_size;                        /* Eigen workspace   */

934:   /***************************************************************************/
935:   /* Compute the minimum eigenvalue of T.                                    */
936:   /***************************************************************************/

938:   vl = 0.0;
939:   vu = 0.0;
940:   il = 1;
941:   iu = 1;

943: #if defined(PETSC_MISSING_LAPACK_STEBZ)
944:   SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_SUP,"STEBZ - Lapack routine is unavailable.");
945: #else
946:   LAPACKstebz_("I", "E", &t_size, &vl, &vu, &il, &iu, &cg->eigen_tol,
947:                cg->diag, cg->offd + 1, &e_valus, &e_splts, e_valu,
948:                e_iblk, e_splt, e_rwrk, e_iwrk, &info);

950:   if ((0 != info) || (1 != e_valus)) {
951:     /*************************************************************************/
952:     /* Calculation of the minimum eigenvalue failed.  Return the             */
953:     /* Steihaug-Toint direction.                                             */
954:     /*************************************************************************/

956:     PetscInfo(ksp, "KSPSolve_GLTR: failed to compute eigenvalue.\n");
957:     ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
958:     return(0);
959:   }

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

963:   /***************************************************************************/
964:   /* Compute the initial value of lambda to make (T + lamba I) positive      */
965:   /* definite.                                                               */
966:   /***************************************************************************/

968:   pert = cg->init_pert;
969:   if (e_valu[0] < 0.0) {
970:     cg->lambda = pert - e_valu[0];
971:   }
972: #endif

974:   while(1) {
975:     for (i = 0; i < t_size; ++i) {
976:       t_diag[i] = cg->diag[i] + cg->lambda;
977:       t_offd[i] = cg->offd[i];
978:     }

980: #if defined(PETSC_MISSING_LAPACK_PTTRF)
981:     SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_SUP,"PTTRF - Lapack routine is unavailable.");
982: #else
983:     LAPACKpttrf_(&t_size, t_diag, t_offd + 1, &info);

985:     if (0 == info) {
986:       break;
987:     }

989:     pert += pert;
990:     cg->lambda = cg->lambda * (1.0 + pert) + pert;
991: #endif
992:   }

994:   /***************************************************************************/
995:   /* Compute the initial step and its norm.                                  */
996:   /***************************************************************************/

998:   nrhs = 1;
999:   nldb = t_size;

1001:   t_soln[0] = -cg->norm_r[0];
1002:   for (i = 1; i < t_size; ++i) {
1003:     t_soln[i] = 0.0;
1004:   }

1006: #if defined(PETSC_MISSING_LAPACK_PTTRS)
1007:   SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_SUP,"PTTRS - Lapack routine is unavailable.");
1008: #else
1009:   LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, t_soln, &nldb, &info);
1010: #endif

1012:   if (0 != info) {
1013:     /*************************************************************************/
1014:     /* Calculation of the initial step failed; return the Steihaug-Toint     */
1015:     /* direction.                                                            */
1016:     /*************************************************************************/

1018:     PetscInfo(ksp, "KSPSolve_GLTR: failed to compute step.\n");
1019:     ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1020:     return(0);
1021:   }

1023:   norm_t = 0.;
1024:   for (i = 0; i < t_size; ++i) {
1025:     norm_t += t_soln[i] * t_soln[i];
1026:   }
1027:   norm_t = PetscSqrtReal(norm_t);

1029:   /***************************************************************************/
1030:   /* Determine the case we are in.                                           */
1031:   /***************************************************************************/

1033:   if (norm_t <= cg->radius) {
1034:     /*************************************************************************/
1035:     /* The step is within the trust region; check if we are in the hard case */
1036:     /* and need to move to the boundary by following a direction of negative */
1037:     /* curvature.                                                            */
1038:     /*************************************************************************/
1039: 
1040:     if ((e_valu[0] <= 0.0) && (norm_t < cg->radius)) {
1041:       /***********************************************************************/
1042:       /* This is the hard case; compute the eigenvector associated with the  */
1043:       /* minimum eigenvalue and move along this direction to the boundary.   */
1044:       /***********************************************************************/

1046: #if defined(PETSC_MISSING_LAPACK_STEIN)
1047:   SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_SUP,"STEIN - Lapack routine is unavailable.");
1048: #else
1049:       LAPACKstein_(&t_size, cg->diag, cg->offd + 1, &e_valus, e_valu,
1050:                    e_iblk, e_splt, e_vect, &nldb,
1051:                    e_rwrk, e_iwrk, e_iwrk + t_size, &info);
1052: #endif

1054:       if (0 != info) {
1055:         /*********************************************************************/
1056:         /* Calculation of the minimum eigenvalue failed.  Return the         */
1057:         /* Steihaug-Toint direction.                                         */
1058:         /*********************************************************************/
1059: 
1060:         PetscInfo(ksp, "KSPSolve_GLTR: failed to compute eigenvector.\n");
1061:         ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1062:         return(0);
1063:       }
1064: 
1065:       coef1 = 0.0;
1066:       coef2 = 0.0;
1067:       coef3 = -cg->radius * cg->radius;
1068:       for (i = 0; i < t_size; ++i) {
1069:         coef1 += e_vect[i] * e_vect[i];
1070:         coef2 += e_vect[i] * t_soln[i];
1071:         coef3 += t_soln[i] * t_soln[i];
1072:       }
1073: 
1074:       coef3 = PetscSqrtReal(coef2 * coef2 - coef1 * coef3);
1075:       root1 = (-coef2 + coef3) / coef1;
1076:       root2 = (-coef2 - coef3) / coef1;

1078:       /***********************************************************************/
1079:       /* Compute objective value for (t_soln + root1 * e_vect)               */
1080:       /***********************************************************************/

1082:       for (i = 0; i < t_size; ++i) {
1083:         e_rwrk[i] = t_soln[i] + root1 * e_vect[i];
1084:       }
1085: 
1086:       obj1 = e_rwrk[0]*(0.5*(cg->diag[0]*e_rwrk[0]+
1087:                              cg->offd[1]*e_rwrk[1])+cg->norm_r[0]);
1088:       for (i = 1; i < t_size - 1; ++i) {
1089:         obj1 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
1090:                                cg->diag[i]*e_rwrk[i]+
1091:                                cg->offd[i+1]*e_rwrk[i+1]);
1092:       }
1093:       obj1 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
1094:                              cg->diag[i]*e_rwrk[i]);

1096:       /***********************************************************************/
1097:       /* Compute objective value for (t_soln + root2 * e_vect)               */
1098:       /***********************************************************************/

1100:       for (i = 0; i < t_size; ++i) {
1101:         e_rwrk[i] = t_soln[i] + root2 * e_vect[i];
1102:       }
1103: 
1104:       obj2 = e_rwrk[0]*(0.5*(cg->diag[0]*e_rwrk[0]+
1105:                              cg->offd[1]*e_rwrk[1])+cg->norm_r[0]);
1106:       for (i = 1; i < t_size - 1; ++i) {
1107:         obj2 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
1108:                                cg->diag[i]*e_rwrk[i]+
1109:                                cg->offd[i+1]*e_rwrk[i+1]);
1110:       }
1111:       obj2 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
1112:                              cg->diag[i]*e_rwrk[i]);
1113: 
1114:       /***********************************************************************/
1115:       /* Choose the point with the best objective function value.            */
1116:       /***********************************************************************/

1118:       if (obj1 <= obj2) {
1119:         for (i = 0; i < t_size; ++i) {
1120:           t_soln[i] += root1 * e_vect[i];
1121:         }
1122:       }
1123:       else {
1124:         for (i = 0; i < t_size; ++i) {
1125:           t_soln[i] += root2 * e_vect[i];
1126:         }
1127:       }
1128:     }
1129:     else {
1130:       /***********************************************************************/
1131:       /* The matrix is positive definite or there was no room to move; the   */
1132:       /* solution is already contained in t_soln.                            */
1133:       /***********************************************************************/
1134:     }
1135:   }
1136:   else {
1137:     /*************************************************************************/
1138:     /* The step is outside the trust-region.  Compute the correct value for  */
1139:     /* lambda by performing Newton's method.                                 */
1140:     /*************************************************************************/

1142:     for (i = 0; i < max_newton_its; ++i) {
1143:       /***********************************************************************/
1144:       /* Check for convergence.                                              */
1145:       /***********************************************************************/

1147:       if (fabs(norm_t - cg->radius) <= cg->newton_tol * cg->radius) {
1148:         break;
1149:       }

1151:       /***********************************************************************/
1152:       /* Compute the update.                                                 */
1153:       /***********************************************************************/

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

1157: #if defined(PETSC_MISSING_LAPACK_PTTRS)
1158:   SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_SUP,"PTTRS - Lapack routine is unavailable.");
1159: #else      
1160:       LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, e_rwrk, &nldb, &info);
1161: #endif

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

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

1174:       /***********************************************************************/
1175:       /* Modify lambda.                                                      */
1176:       /***********************************************************************/

1178:       norm_w = 0.;
1179:       for (j = 0; j < t_size; ++j) {
1180:         norm_w += t_soln[j] * e_rwrk[j];
1181:       }
1182: 
1183:       cg->lambda += (norm_t - cg->radius)/cg->radius * (norm_t * norm_t) / norm_w;

1185:       /***********************************************************************/
1186:       /* Factor T + lambda I                                                 */
1187:       /***********************************************************************/
1188: 
1189:       for (j = 0; j < t_size; ++j) {
1190:         t_diag[j] = cg->diag[j] + cg->lambda;
1191:         t_offd[j] = cg->offd[j];
1192:       }

1194: #if defined(PETSC_MISSING_LAPACK_PTTRF)
1195:   SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_SUP,"PTTRF - Lapack routine is unavailable.");
1196: #else
1197:       LAPACKpttrf_(&t_size, t_diag, t_offd + 1, &info);
1198: #endif

1200:       if (0 != info) {
1201:         /*********************************************************************/
1202:         /* Calculation of factorization failed; return the Steihaug-Toint    */
1203:         /* direction.                                                        */
1204:         /*********************************************************************/

1206:         PetscInfo(ksp, "KSPSolve_GLTR: factorization failed.\n");
1207:         ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1208:         return(0);
1209:       }

1211:       /***********************************************************************/
1212:       /* Compute the new step and its norm.                                  */
1213:       /***********************************************************************/

1215:       t_soln[0] = -cg->norm_r[0];
1216:       for (j = 1; j < t_size; ++j) {
1217:         t_soln[j] = 0.0;
1218:       }

1220: #if defined(PETSC_MISSING_LAPACK_PTTRS)
1221:   SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_SUP,"PTTRS - Lapack routine is unavailable.");
1222: #else
1223:       LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, t_soln, &nldb, &info);
1224: #endif

1226:       if (0 != info) {
1227:         /*********************************************************************/
1228:         /* Calculation of the step failed; return the Steihaug-Toint         */
1229:         /* direction.                                                        */
1230:         /*********************************************************************/
1231: 
1232:         PetscInfo(ksp, "KSPSolve_GLTR: failed to compute step.\n");
1233:         ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1234:         return(0);
1235:       }

1237:       norm_t = 0.;
1238:       for (j = 0; j < t_size; ++j) {
1239:         norm_t += t_soln[j] * t_soln[j];
1240:       }
1241:       norm_t = PetscSqrtReal(norm_t);
1242:     }

1244:     /*************************************************************************/
1245:     /* Check for convergence.                                                */
1246:     /*************************************************************************/

1248:     if (fabs(norm_t - cg->radius) > cg->newton_tol * cg->radius) {
1249:       /***********************************************************************/
1250:       /* Newton method failed to converge in iteration limit.                */
1251:       /***********************************************************************/

1253:       PetscInfo(ksp, "KSPSolve_GLTR: failed to converge.\n");
1254:       ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1255:       return(0);
1256:     }
1257:   }

1259:   /***************************************************************************/
1260:   /* Recover the norm of the direction and objective function value.         */
1261:   /***************************************************************************/

1263:   cg->norm_d = norm_t;

1265:   cg->o_fcn = t_soln[0]*(0.5*(cg->diag[0]*t_soln[0]+cg->offd[1]*t_soln[1])+cg->norm_r[0]);
1266:   for (i = 1; i < t_size - 1; ++i) {
1267:     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]);
1268:   }
1269:   cg->o_fcn += 0.5*t_soln[i]*(cg->offd[i]*t_soln[i-1]+cg->diag[i]*t_soln[i]);

1271:   /***************************************************************************/
1272:   /* Recover the direction.                                                  */
1273:   /***************************************************************************/

1275:   sigma = -1;

1277:   /***************************************************************************/
1278:   /* Start conjugate gradient method from the beginning                      */
1279:   /***************************************************************************/

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

1284:   /***************************************************************************/
1285:   /* Accumulate Q * s                                                        */
1286:   /***************************************************************************/

1288:   VecCopy(z, d);
1289:   VecScale(d, sigma * t_soln[0] / cg->norm_r[0]);
1290: 
1291:   /***************************************************************************/
1292:   /* Compute the first direction.                                            */
1293:   /***************************************************************************/

1295:   VecCopy(z, p);                        /* p = z             */
1296:   KSP_MatMult(ksp, Qmat, p, z);        /* z = Q * p         */
1297:   ++ksp->its;

1299:   for (i = 0; i < l_size - 1; ++i) {
1300:     /*************************************************************************/
1301:     /* Update the residual and direction.                                    */
1302:     /*************************************************************************/

1304:     alpha = cg->alpha[i];
1305:     if (alpha >= 0.0) {
1306:       sigma = -sigma;
1307:     }

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

1312:     /*************************************************************************/
1313:     /* Accumulate Q * s                                                      */
1314:     /*************************************************************************/

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

1318:     /*************************************************************************/
1319:     /* Update p.                                                             */
1320:     /*************************************************************************/

1322:     beta = cg->beta[i];
1323:     VecAYPX(p, beta, z);          /* p = z + beta p    */
1324:     KSP_MatMult(ksp, Qmat, p, z);        /* z = Q * p         */
1325:     ++ksp->its;
1326:   }

1328:   /***************************************************************************/
1329:   /* Update the residual and direction.                                      */
1330:   /***************************************************************************/

1332:   alpha = cg->alpha[i];
1333:   if (alpha >= 0.0) {
1334:     sigma = -sigma;
1335:   }

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

1340:   /***************************************************************************/
1341:   /* Accumulate Q * s                                                        */
1342:   /***************************************************************************/

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

1346:   /***************************************************************************/
1347:   /* Set the termination reason.                                             */
1348:   /***************************************************************************/

1350:   ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1351:   return(0);
1352: #endif
1353: }

1357: PetscErrorCode KSPSetUp_GLTR(KSP ksp)
1358: {
1359:   KSP_GLTR       *cg = (KSP_GLTR *)ksp->data;
1360:   PetscInt       max_its;

1364:   /***************************************************************************/
1365:   /* Determine the total maximum number of iterations.                       */
1366:   /***************************************************************************/

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

1370:   /***************************************************************************/
1371:   /* Set work vectors needed by conjugate gradient method and allocate       */
1372:   /* workspace for Lanczos matrix.                                           */
1373:   /***************************************************************************/

1375:   KSPDefaultGetWork(ksp, 3);

1377:   PetscMalloc5(max_its,PetscReal,&cg->diag,max_its,PetscReal,&cg->offd,max_its,PetscReal,&cg->alpha,max_its,PetscReal,&cg->beta,max_its,PetscReal,&cg->norm_r);
1378:   PetscMemzero(cg->diag, max_its*sizeof(PetscReal));
1379:   PetscMemzero(cg->offd, max_its*sizeof(PetscReal));
1380:   PetscMemzero(cg->alpha, max_its*sizeof(PetscReal));
1381:   PetscMemzero(cg->beta, max_its*sizeof(PetscReal));
1382:   PetscMemzero(cg->norm_r, max_its*sizeof(PetscReal));
1383:   PetscLogObjectMemory(ksp, 5*max_its*sizeof(PetscReal));
1384:   return(0);
1385: }

1389: PetscErrorCode KSPDestroy_GLTR(KSP ksp)
1390: {
1391:   KSP_GLTR       *cg = (KSP_GLTR *)ksp->data;

1395:   /***************************************************************************/
1396:   /* Free memory allocated for the data.                                     */
1397:   /***************************************************************************/

1399:   PetscFree5(cg->diag,cg->offd,cg->alpha,cg->beta,cg->norm_r);
1400:   if (cg->alloced) {
1401:     PetscFree2(cg->rwork,cg->iwork);
1402:   }

1404:   /***************************************************************************/
1405:   /* Clear composed functions                                                */
1406:   /***************************************************************************/

1408:   PetscObjectComposeFunctionDynamic((PetscObject)ksp,"KSPGLTRSetRadius_C","",PETSC_NULL);
1409:   PetscObjectComposeFunctionDynamic((PetscObject)ksp,"KSPGLTRGetNormD_C","",PETSC_NULL);
1410:   PetscObjectComposeFunctionDynamic((PetscObject)ksp,"KSPGLTRGetObjFcn_C","",PETSC_NULL);
1411:   PetscObjectComposeFunctionDynamic((PetscObject)ksp,"KSPGLTRGetMinEig_C","",PETSC_NULL);
1412:   PetscObjectComposeFunctionDynamic((PetscObject)ksp,"KSPGLTRGetLambda_C","",PETSC_NULL);

1414:   /***************************************************************************/
1415:   /* Destroy KSP object.                                                     */
1416:   /***************************************************************************/
1417:   KSPDefaultDestroy(ksp);
1418:   return(0);
1419: }

1421: EXTERN_C_BEGIN
1424: PetscErrorCode  KSPGLTRSetRadius_GLTR(KSP ksp, PetscReal radius)
1425: {
1426:   KSP_GLTR *cg = (KSP_GLTR *)ksp->data;

1429:   cg->radius = radius;
1430:   return(0);
1431: }

1435: PetscErrorCode  KSPGLTRGetNormD_GLTR(KSP ksp, PetscReal *norm_d)
1436: {
1437:   KSP_GLTR *cg = (KSP_GLTR *)ksp->data;

1440:   *norm_d = cg->norm_d;
1441:   return(0);
1442: }

1446: PetscErrorCode  KSPGLTRGetObjFcn_GLTR(KSP ksp, PetscReal *o_fcn)
1447: {
1448:   KSP_GLTR *cg = (KSP_GLTR *)ksp->data;

1451:   *o_fcn = cg->o_fcn;
1452:   return(0);
1453: }

1457: PetscErrorCode  KSPGLTRGetMinEig_GLTR(KSP ksp, PetscReal *e_min)
1458: {
1459:   KSP_GLTR *cg = (KSP_GLTR *)ksp->data;

1462:   *e_min = cg->e_min;
1463:   return(0);
1464: }

1468: PetscErrorCode  KSPGLTRGetLambda_GLTR(KSP ksp, PetscReal *lambda)
1469: {
1470:   KSP_GLTR *cg = (KSP_GLTR *)ksp->data;

1473:   *lambda = cg->lambda;
1474:   return(0);
1475: }
1476: EXTERN_C_END

1480: PetscErrorCode KSPSetFromOptions_GLTR(KSP ksp)
1481: {
1483:   KSP_GLTR       *cg = (KSP_GLTR *)ksp->data;

1486:   PetscOptionsHead("KSP GLTR options");

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

1490:   PetscOptionsReal("-ksp_gltr_init_pert", "Initial perturbation", "", cg->init_pert, &cg->init_pert, PETSC_NULL);
1491:   PetscOptionsReal("-ksp_gltr_eigen_tol", "Eigenvalue tolerance", "", cg->eigen_tol, &cg->eigen_tol, PETSC_NULL);
1492:   PetscOptionsReal("-ksp_gltr_newton_tol", "Newton tolerance", "", cg->newton_tol, &cg->newton_tol, PETSC_NULL);

1494:   PetscOptionsInt("-ksp_gltr_max_lanczos_its", "Maximum Lanczos Iters", "", cg->max_lanczos_its, &cg->max_lanczos_its, PETSC_NULL);
1495:   PetscOptionsInt("-ksp_gltr_max_newton_its", "Maximum Newton Iters", "", cg->max_newton_its, &cg->max_newton_its, PETSC_NULL);

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

1499:   PetscOptionsTail();
1500:   return(0);
1501: }

1503: /*MC
1504:      KSPGLTR -   Code to run conjugate gradient method subject to a constraint
1505:          on the solution norm. This is used in Trust Region methods for
1506:          nonlinear equations, SNESTR

1508:    Options Database Keys:
1509: .      -ksp_gltr_radius <r> - Trust Region Radius

1511:    Notes: This is rarely used directly

1513:   Use preconditioned conjugate gradient to compute
1514:   an approximate minimizer of the quadratic function

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

1518:    subject to the trust region constraint

1520:             || s || <= delta,

1522:    where

1524:      delta is the trust region radius,
1525:      g is the gradient vector,
1526:      H is the Hessian approximation,
1527:      M is the positive definite preconditioner matrix.

1529:    KSPConvergedReason may be
1530: $  KSP_CONVERGED_CG_NEG_CURVE if convergence is reached along a negative curvature direction,
1531: $  KSP_CONVERGED_CG_CONSTRAINED if convergence is reached along a constrained step,
1532: $  other KSP converged/diverged reasons

1534:   Notes:
1535:   The preconditioner supplied should be symmetric and positive definite.

1537:    Level: developer

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

1542: EXTERN_C_BEGIN
1545: PetscErrorCode  KSPCreate_GLTR(KSP ksp)
1546: {
1548:   KSP_GLTR       *cg;

1551:   PetscNewLog(ksp, KSP_GLTR, &cg);
1552:   cg->radius = 0.0;
1553:   cg->dtype = GLTR_UNPRECONDITIONED_DIRECTION;

1555:   cg->init_pert = 1.0e-8;
1556:   cg->eigen_tol = 1.0e-10;
1557:   cg->newton_tol = 1.0e-6;

1559:   cg->alloced = 0;
1560:   cg->init_alloc = 1024;

1562:   cg->max_lanczos_its = 20;
1563:   cg->max_newton_its  = 10;

1565:   ksp->data = (void *) cg;
1566:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);
1567:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,1);
1568:   KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,1);

1570:   /***************************************************************************/
1571:   /* Sets the functions that are associated with this data structure         */
1572:   /* (in C++ this is the same as defining virtual functions).                */
1573:   /***************************************************************************/

1575:   ksp->ops->setup          = KSPSetUp_GLTR;
1576:   ksp->ops->solve          = KSPSolve_GLTR;
1577:   ksp->ops->destroy        = KSPDestroy_GLTR;
1578:   ksp->ops->setfromoptions = KSPSetFromOptions_GLTR;
1579:   ksp->ops->buildsolution  = KSPDefaultBuildSolution;
1580:   ksp->ops->buildresidual  = KSPDefaultBuildResidual;
1581:   ksp->ops->view           = 0;

1583:   PetscObjectComposeFunctionDynamic((PetscObject)ksp,
1584:                                            "KSPGLTRSetRadius_C",
1585:                                            "KSPGLTRSetRadius_GLTR",
1586:                                             KSPGLTRSetRadius_GLTR);
1587:   PetscObjectComposeFunctionDynamic((PetscObject)ksp,
1588:                                            "KSPGLTRGetNormD_C",
1589:                                            "KSPGLTRGetNormD_GLTR",
1590:                                             KSPGLTRGetNormD_GLTR);
1591:   PetscObjectComposeFunctionDynamic((PetscObject)ksp,
1592:                                            "KSPGLTRGetObjFcn_C",
1593:                                            "KSPGLTRGetObjFcn_GLTR",
1594:                                             KSPGLTRGetObjFcn_GLTR);
1595:   PetscObjectComposeFunctionDynamic((PetscObject)ksp,
1596:                                            "KSPGLTRGetMinEig_C",
1597:                                            "KSPGLTRGetMinEig_GLTR",
1598:                                             KSPGLTRGetMinEig_GLTR);
1599:   PetscObjectComposeFunctionDynamic((PetscObject)ksp,
1600:                                            "KSPGLTRGetLambda_C",
1601:                                            "KSPGLTRGetLambda_GLTR",
1602:                                             KSPGLTRGetLambda_GLTR);
1603:   return(0);
1604: }
1605: EXTERN_C_END