Actual source code: gltr.c

petsc-master 2020-07-09
Report Typos and Errors

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

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

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

 11: /*@
 12:     KSPGLTRGetMinEig - Get minimum eigenvalue.

 14:     Collective on ksp

 16:     Input Parameters:
 17: +   ksp   - the iterative context
 18: -   e_min - the minimum eigenvalue

 20:     Level: advanced

 22: @*/
 23: PetscErrorCode  KSPGLTRGetMinEig(KSP ksp, PetscReal *e_min)
 24: {

 29:   PetscUseMethod(ksp,"KSPGLTRGetMinEig_C",(KSP,PetscReal*),(ksp,e_min));
 30:   return(0);
 31: }

 33: /*@
 34:     KSPGLTRGetLambda - Get multiplier on trust-region constraint.

 36:     Not Collective

 38:     Input Parameters:
 39: +   ksp    - the iterative context
 40: -   lambda - the multiplier

 42:     Level: advanced

 44: @*/
 45: PetscErrorCode  KSPGLTRGetLambda(KSP ksp, PetscReal *lambda)
 46: {

 51:   PetscUseMethod(ksp,"KSPGLTRGetLambda_C",(KSP,PetscReal*),(ksp,lambda));
 52:   return(0);
 53: }

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

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

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

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

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

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

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

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

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

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

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

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

115:   cg->o_fcn = 0.0;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

314:   piv = 1.0;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

405:       ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
406:       PetscInfo1(ksp, "KSPCGSolve_GLTR: negative curvature: kappa=%g\n", (double)kappa);

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

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

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

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

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

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

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

443:       ksp->reason = KSP_CONVERGED_CG_CONSTRAINED;
444:       PetscInfo1(ksp, "KSPCGSolve_GLTR: constrained step: radius=%g\n", (double)cg->radius);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

568:     if (ksp->its >= max_cg_its) {
569:       ksp->reason = KSP_DIVERGED_ITS;
570:       PetscInfo1(ksp, "KSPCGSolve_GLTR: iterlim: its=%D\n", ksp->its);
571:       break;
572:     }

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

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

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

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

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

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

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

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

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

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

621:       ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
622:       PetscInfo1(ksp, "KSPCGSolve_GLTR: cg breakdown: kappa=%g\n", (double)kappa);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

800:     while (t_size > cg->alloced) {
801:       cg->alloced += cg->init_alloc;
802:     }

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

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

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

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

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

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

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

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

840:     PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute eigenvalue.\n");
841:     ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
842:     return(0);
843:   }

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

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

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

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

861:     PetscStackCallBLAS("LAPACKpttrf",LAPACKpttrf_(&t_size, t_diag, t_offd + 1, &info));

863:     if (0 == info) break;

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

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

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

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

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

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

887:     PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute step.\n");
888:     ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
889:     return(0);
890:   }

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

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

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

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

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

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

921:         PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute eigenvector.\n");
922:         ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
923:         return(0);
924:       }

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

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

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

943:       for (i = 0; i < t_size; ++i) {
944:         e_rwrk[i] = t_soln[i] + root1 * e_vect[i];
945:       }

947:       obj1 = e_rwrk[0]*(0.5*(cg->diag[0]*e_rwrk[0]+
948:                              cg->offd[1]*e_rwrk[1])+cg->norm_r[0]);
949:       for (i = 1; i < t_size - 1; ++i) {
950:         obj1 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
951:                                cg->diag[i]*e_rwrk[i]+
952:                                cg->offd[i+1]*e_rwrk[i+1]);
953:       }
954:       obj1 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
955:                              cg->diag[i]*e_rwrk[i]);

957:       /***********************************************************************/
958:       /* Compute objective value for (t_soln + root2 * e_vect)               */
959:       /***********************************************************************/

961:       for (i = 0; i < t_size; ++i) {
962:         e_rwrk[i] = t_soln[i] + root2 * e_vect[i];
963:       }

965:       obj2 = e_rwrk[0]*(0.5*(cg->diag[0]*e_rwrk[0]+
966:                              cg->offd[1]*e_rwrk[1])+cg->norm_r[0]);
967:       for (i = 1; i < t_size - 1; ++i) {
968:         obj2 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
969:                                cg->diag[i]*e_rwrk[i]+
970:                                cg->offd[i+1]*e_rwrk[i+1]);
971:       }
972:       obj2 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
973:                              cg->diag[i]*e_rwrk[i]);

975:       /***********************************************************************/
976:       /* Choose the point with the best objective function value.            */
977:       /***********************************************************************/

979:       if (obj1 <= obj2) {
980:         for (i = 0; i < t_size; ++i) {
981:           t_soln[i] += root1 * e_vect[i];
982:         }
983:       }
984:       else {
985:         for (i = 0; i < t_size; ++i) {
986:           t_soln[i] += root2 * e_vect[i];
987:         }
988:       }
989:     } else {
990:       /***********************************************************************/
991:       /* The matrix is positive definite or there was no room to move; the   */
992:       /* solution is already contained in t_soln.                            */
993:       /***********************************************************************/
994:     }
995:   } else {
996:     /*************************************************************************/
997:     /* The step is outside the trust-region.  Compute the correct value for  */
998:     /* lambda by performing Newton's method.                                 */
999:     /*************************************************************************/

1001:     for (i = 0; i < max_newton_its; ++i) {
1002:       /***********************************************************************/
1003:       /* Check for convergence.                                              */
1004:       /***********************************************************************/

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

1008:       /***********************************************************************/
1009:       /* Compute the update.                                                 */
1010:       /***********************************************************************/

1012:       PetscArraycpy(e_rwrk, t_soln, t_size);

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

1016:       if (0 != info) {
1017:         /*********************************************************************/
1018:         /* Calculation of the step failed; return the Steihaug-Toint         */
1019:         /* direction.                                                        */
1020:         /*********************************************************************/

1022:         PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute step.\n");
1023:         ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1024:         return(0);
1025:       }

1027:       /***********************************************************************/
1028:       /* Modify lambda.                                                      */
1029:       /***********************************************************************/

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

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

1036:       /***********************************************************************/
1037:       /* Factor T + lambda I                                                 */
1038:       /***********************************************************************/

1040:       for (j = 0; j < t_size; ++j) {
1041:         t_diag[j] = cg->diag[j] + cg->lambda;
1042:         t_offd[j] = cg->offd[j];
1043:       }

1045:       PetscStackCallBLAS("LAPACKpttrf",LAPACKpttrf_(&t_size, t_diag, t_offd + 1, &info));

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

1053:         PetscInfo(ksp, "KSPCGSolve_GLTR: factorization failed.\n");
1054:         ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1055:         return(0);
1056:       }

1058:       /***********************************************************************/
1059:       /* Compute the new step and its norm.                                  */
1060:       /***********************************************************************/

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

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

1067:       if (0 != info) {
1068:         /*********************************************************************/
1069:         /* Calculation of the step failed; return the Steihaug-Toint         */
1070:         /* direction.                                                        */
1071:         /*********************************************************************/

1073:         PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute step.\n");
1074:         ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1075:         return(0);
1076:       }

1078:       norm_t = 0.;
1079:       for (j = 0; j < t_size; ++j) norm_t += t_soln[j] * t_soln[j];
1080:       norm_t = PetscSqrtReal(norm_t);
1081:     }

1083:     /*************************************************************************/
1084:     /* Check for convergence.                                                */
1085:     /*************************************************************************/

1087:     if (PetscAbsReal(norm_t - cg->radius) > cg->newton_tol * cg->radius) {
1088:       /***********************************************************************/
1089:       /* Newton method failed to converge in iteration limit.                */
1090:       /***********************************************************************/

1092:       PetscInfo(ksp, "KSPCGSolve_GLTR: failed to converge.\n");
1093:       ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1094:       return(0);
1095:     }
1096:   }

1098:   /***************************************************************************/
1099:   /* Recover the norm of the direction and objective function value.         */
1100:   /***************************************************************************/

1102:   cg->norm_d = norm_t;

1104:   cg->o_fcn = t_soln[0]*(0.5*(cg->diag[0]*t_soln[0]+cg->offd[1]*t_soln[1])+cg->norm_r[0]);
1105:   for (i = 1; i < t_size - 1; ++i) {
1106:     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]);
1107:   }
1108:   cg->o_fcn += 0.5*t_soln[i]*(cg->offd[i]*t_soln[i-1]+cg->diag[i]*t_soln[i]);

1110:   /***************************************************************************/
1111:   /* Recover the direction.                                                  */
1112:   /***************************************************************************/

1114:   sigma = -1;

1116:   /***************************************************************************/
1117:   /* Start conjugate gradient method from the beginning                      */
1118:   /***************************************************************************/

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

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

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

1130:   /***************************************************************************/
1131:   /* Compute the first direction.                                            */
1132:   /***************************************************************************/

1134:   VecCopy(z, p);                   /* p = z             */
1135:   KSP_MatMult(ksp, Qmat, p, z);    /* z = Q * p         */
1136:   ++ksp->its;

1138:   for (i = 0; i < l_size - 1; ++i) {
1139:     /*************************************************************************/
1140:     /* Update the residual and direction.                                    */
1141:     /*************************************************************************/

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

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

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

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

1155:     /*************************************************************************/
1156:     /* Update p.                                                             */
1157:     /*************************************************************************/

1159:     beta = cg->beta[i];
1160:     VecAYPX(p, beta, z);          /* p = z + beta p    */
1161:     KSP_MatMult(ksp, Qmat, p, z);  /* z = Q * p         */
1162:     ++ksp->its;
1163:   }

1165:   /***************************************************************************/
1166:   /* Update the residual and direction.                                      */
1167:   /***************************************************************************/

1169:   alpha = cg->alpha[i];
1170:   if (alpha >= 0.0) sigma = -sigma;

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

1175:   /***************************************************************************/
1176:   /* Accumulate Q * s                                                        */
1177:   /***************************************************************************/

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

1181:   /***************************************************************************/
1182:   /* Set the termination reason.                                             */
1183:   /***************************************************************************/

1185:   ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1186:   return(0);
1187: #endif
1188: }

1190: static PetscErrorCode KSPCGSetUp_GLTR(KSP ksp)
1191: {
1192:   KSPCG_GLTR     *cg = (KSPCG_GLTR*)ksp->data;
1193:   PetscInt       max_its;

1197:   /***************************************************************************/
1198:   /* Determine the total maximum number of iterations.                       */
1199:   /***************************************************************************/

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

1203:   /***************************************************************************/
1204:   /* Set work vectors needed by conjugate gradient method and allocate       */
1205:   /* workspace for Lanczos matrix.                                           */
1206:   /***************************************************************************/

1208:   KSPSetWorkVecs(ksp, 3);
1209:   if (cg->diag) {
1210:     PetscArrayzero(cg->diag, max_its);
1211:     PetscArrayzero(cg->offd, max_its);
1212:     PetscArrayzero(cg->alpha, max_its);
1213:     PetscArrayzero(cg->beta, max_its);
1214:     PetscArrayzero(cg->norm_r, max_its);
1215:   } else {
1216:     PetscCalloc5(max_its,&cg->diag,max_its,&cg->offd,max_its,&cg->alpha,max_its,&cg->beta,max_its,&cg->norm_r);
1217:     PetscLogObjectMemory((PetscObject)ksp, 5*max_its*sizeof(PetscReal));
1218:   }
1219:   return(0);
1220: }

1222: static PetscErrorCode KSPCGDestroy_GLTR(KSP ksp)
1223: {
1224:   KSPCG_GLTR     *cg = (KSPCG_GLTR*)ksp->data;

1228:   /***************************************************************************/
1229:   /* Free memory allocated for the data.                                     */
1230:   /***************************************************************************/

1232:   PetscFree5(cg->diag,cg->offd,cg->alpha,cg->beta,cg->norm_r);
1233:   if (cg->alloced) {
1234:     PetscFree2(cg->rwork,cg->iwork);
1235:   }

1237:   /***************************************************************************/
1238:   /* Clear composed functions                                                */
1239:   /***************************************************************************/

1241:   PetscObjectComposeFunction((PetscObject)ksp,"KSPCGSetRadius_C",NULL);
1242:   PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetNormD_C",NULL);
1243:   PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetObjFcn_C",NULL);
1244:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetMinEig_C",NULL);
1245:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetLambda_C",NULL);

1247:   /***************************************************************************/
1248:   /* Destroy KSP object.                                                     */
1249:   /***************************************************************************/
1250:   KSPDestroyDefault(ksp);
1251:   return(0);
1252: }

1254: static PetscErrorCode  KSPCGSetRadius_GLTR(KSP ksp, PetscReal radius)
1255: {
1256:   KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;

1259:   cg->radius = radius;
1260:   return(0);
1261: }

1263: static PetscErrorCode  KSPCGGetNormD_GLTR(KSP ksp, PetscReal *norm_d)
1264: {
1265:   KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;

1268:   *norm_d = cg->norm_d;
1269:   return(0);
1270: }

1272: static PetscErrorCode  KSPCGGetObjFcn_GLTR(KSP ksp, PetscReal *o_fcn)
1273: {
1274:   KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;

1277:   *o_fcn = cg->o_fcn;
1278:   return(0);
1279: }

1281: static PetscErrorCode  KSPGLTRGetMinEig_GLTR(KSP ksp, PetscReal *e_min)
1282: {
1283:   KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;

1286:   *e_min = cg->e_min;
1287:   return(0);
1288: }

1290: static PetscErrorCode  KSPGLTRGetLambda_GLTR(KSP ksp, PetscReal *lambda)
1291: {
1292:   KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;

1295:   *lambda = cg->lambda;
1296:   return(0);
1297: }

1299: static PetscErrorCode KSPCGSetFromOptions_GLTR(PetscOptionItems *PetscOptionsObject,KSP ksp)
1300: {
1302:   KSPCG_GLTR       *cg = (KSPCG_GLTR*)ksp->data;

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

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

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

1311:   PetscOptionsReal("-ksp_cg_gltr_init_pert", "Initial perturbation", "", cg->init_pert, &cg->init_pert, NULL);
1312:   PetscOptionsReal("-ksp_cg_gltr_eigen_tol", "Eigenvalue tolerance", "", cg->eigen_tol, &cg->eigen_tol, NULL);
1313:   PetscOptionsReal("-ksp_cg_gltr_newton_tol", "Newton tolerance", "", cg->newton_tol, &cg->newton_tol, NULL);

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

1318:   PetscOptionsTail();
1319:   return(0);
1320: }

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

1327:    Options Database Keys:
1328: .      -ksp_cg_radius <r> - Trust Region Radius

1330:    Notes:
1331:     This is rarely used directly

1333:   Use preconditioned conjugate gradient to compute
1334:   an approximate minimizer of the quadratic function

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

1338:    subject to the trust region constraint

1340:             || s || <= delta,

1342:    where

1344:      delta is the trust region radius,
1345:      g is the gradient vector,
1346:      H is the Hessian approximation,
1347:      M is the positive definite preconditioner matrix.

1349:    KSPConvergedReason may be
1350: $  KSP_CONVERGED_CG_NEG_CURVE if convergence is reached along a negative curvature direction,
1351: $  KSP_CONVERGED_CG_CONSTRAINED if convergence is reached along a constrained step,
1352: $  other KSP converged/diverged reasons

1354:   Notes:
1355:   The preconditioner supplied should be symmetric and positive definite.

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

1361:    Level: developer

1363: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPCGSetRadius(), KSPCGGetNormD(), KSPCGGetObjFcn(), KSPGLTRGetMinEig(), KSPGLTRGetLambda()
1364: M*/

1366: PETSC_EXTERN PetscErrorCode KSPCreate_GLTR(KSP ksp)
1367: {
1369:   KSPCG_GLTR       *cg;

1372:   PetscNewLog(ksp,&cg);
1373:   cg->radius = 0.0;
1374:   cg->dtype  = GLTR_UNPRECONDITIONED_DIRECTION;

1376:   cg->init_pert  = 1.0e-8;
1377:   cg->eigen_tol  = 1.0e-10;
1378:   cg->newton_tol = 1.0e-6;

1380:   cg->alloced    = 0;
1381:   cg->init_alloc = 1024;

1383:   cg->max_lanczos_its = 20;
1384:   cg->max_newton_its  = 10;

1386:   ksp->data = (void*) cg;
1387:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,3);
1388:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);
1389:   KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,2);
1390:   KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);

1392:   /***************************************************************************/
1393:   /* Sets the functions that are associated with this data structure         */
1394:   /* (in C++ this is the same as defining virtual functions).                */
1395:   /***************************************************************************/

1397:   ksp->ops->setup          = KSPCGSetUp_GLTR;
1398:   ksp->ops->solve          = KSPCGSolve_GLTR;
1399:   ksp->ops->destroy        = KSPCGDestroy_GLTR;
1400:   ksp->ops->setfromoptions = KSPCGSetFromOptions_GLTR;
1401:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
1402:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
1403:   ksp->ops->view           = NULL;

1405:   PetscObjectComposeFunction((PetscObject)ksp,"KSPCGSetRadius_C",KSPCGSetRadius_GLTR);
1406:   PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetNormD_C", KSPCGGetNormD_GLTR);
1407:   PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetObjFcn_C",KSPCGGetObjFcn_GLTR);
1408:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetMinEig_C",KSPGLTRGetMinEig_GLTR);
1409:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetLambda_C",KSPGLTRGetLambda_GLTR);
1410:   return(0);
1411: }