Actual source code: gltr.c

petsc-3.9.2 2018-05-20
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:     KSPCGGLTRGetMinEig - 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: .keywords: KSP, GLTR, get, minimum eigenvalue
 23: @*/
 24: PetscErrorCode  KSPCGGLTRGetMinEig(KSP ksp, PetscReal *e_min)
 25: {

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

 34: /*@
 35:     KSPCGGLTRGetLambda - Get multiplier on trust-region constraint.

 37:     Not Collective

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

 43:     Level: advanced

 45: .keywords: KSP, GLTR, get, multiplier
 46: @*/
 47: PetscErrorCode  KSPCGGLTRGetLambda(KSP ksp, PetscReal *lambda)
 48: {

 53:   PetscUseMethod(ksp,"KSPCGGLTRGetLambda_C",(KSP,PetscReal*),(ksp,lambda));
 54:   return(0);
 55: }

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

 67:   Mat            Qmat, Mmat;
 68:   Vec            r, z, p, d;
 69:   PC             pc;

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

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

 82: #if !defined(PETSC_MISSING_LAPACK_STEBZ)
 83:   PetscBLASInt e_valus, e_splts;
 84: #endif
 85:   PetscBool diagonalscale;

 88:   /***************************************************************************/
 89:   /* Check the arguments and parameters.                                     */
 90:   /***************************************************************************/

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

 96:   /***************************************************************************/
 97:   /* Get the workspace vectors and initialize variables                      */
 98:   /***************************************************************************/

100:   r2 = cg->radius * cg->radius;
101:   r  = ksp->work[0];
102:   z  = ksp->work[1];
103:   p  = ksp->work[2];
104:   d  = ksp->vec_sol;
105:   pc = ksp->pc;

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

109:   VecGetSize(d, &max_cg_its);
110:   max_cg_its      = PetscMin(max_cg_its, ksp->max_it);
111:   max_lanczos_its = cg->max_lanczos_its;
112:   max_newton_its  = cg->max_newton_its;
113:   ksp->its        = 0;

115:   /***************************************************************************/
116:   /* Initialize objective function direction, and minimum eigenvalue.        */
117:   /***************************************************************************/

119:   cg->o_fcn = 0.0;

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

124:   cg->e_min  = 0.0;
125:   cg->lambda = 0.0;

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

135:   VecCopy(ksp->vec_rhs, r);        /* r = -grad         */
136:   VecDot(r, r, &rr);               /* rr = r^T r        */
137:   KSPCheckDot(ksp,rr);

139:   /***************************************************************************/
140:   /* Check the preconditioner for numerical problems and for positive        */
141:   /* definiteness.  The check for not-a-number and infinite values need be   */
142:   /* performed only once.                                                    */
143:   /***************************************************************************/

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

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

156:     if (cg->radius) {
157:       if (r2 >= rr) {
158:         alpha      = 1.0;
159:         cg->norm_d = PetscSqrtReal(rr);
160:       } else {
161:         alpha      = PetscSqrtReal(r2 / rr);
162:         cg->norm_d = cg->radius;
163:       }

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

167:       /***********************************************************************/
168:       /* Compute objective function.                                         */
169:       /***********************************************************************/

171:       KSP_MatMult(ksp, Qmat, d, z);
172:       VecAYPX(z, -0.5, ksp->vec_rhs);
173:       VecDot(d, z, &cg->o_fcn);
174:       cg->o_fcn = -cg->o_fcn;
175:       ++ksp->its;
176:     }
177:     return(0);
178:   }

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

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

191:     if (cg->radius) {
192:       if (r2 >= rr) {
193:         alpha      = 1.0;
194:         cg->norm_d = PetscSqrtReal(rr);
195:       } else {
196:         alpha      = PetscSqrtReal(r2 / rr);
197:         cg->norm_d = cg->radius;
198:       }

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

202:       /***********************************************************************/
203:       /* Compute objective function.                                         */
204:       /***********************************************************************/

206:       KSP_MatMult(ksp, Qmat, d, z);
207:       VecAYPX(z, -0.5, ksp->vec_rhs);
208:       VecDot(d, z, &cg->o_fcn);
209:       cg->o_fcn = -cg->o_fcn;
210:       ++ksp->its;
211:     }
212:     return(0);
213:   }

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

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

224:   switch (ksp->normtype) {
225:   case KSP_NORM_PRECONDITIONED:
226:     VecNorm(z, NORM_2, &norm_r);   /* norm_r = |z|      */
227:     break;

229:   case KSP_NORM_UNPRECONDITIONED:
230:     norm_r = PetscSqrtReal(rr);                                 /* norm_r = |r|      */
231:     break;

233:   case KSP_NORM_NATURAL:
234:     norm_r = cg->norm_r[0];                             /* norm_r = |r|_M    */
235:     break;

237:   default:
238:     norm_r = 0.0;
239:     break;
240:   }

242:   KSPLogResidualHistory(ksp, norm_r);
243:   KSPMonitor(ksp, ksp->its, norm_r);
244:   ksp->rnorm = norm_r;

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

248:   /***************************************************************************/
249:   /* Compute the first direction and update the iteration.                   */
250:   /***************************************************************************/

252:   VecCopy(z, p);                   /* p = z             */
253:   KSP_MatMult(ksp, Qmat, p, z);    /* z = Q * p         */
254:   ++ksp->its;

256:   /***************************************************************************/
257:   /* Check the matrix for numerical problems.                                */
258:   /***************************************************************************/

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

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

271:     if (cg->radius) {
272:       if (r2 >= rr) {
273:         alpha      = 1.0;
274:         cg->norm_d = PetscSqrtReal(rr);
275:       } else {
276:         alpha      = PetscSqrtReal(r2 / rr);
277:         cg->norm_d = cg->radius;
278:       }

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

282:       /***********************************************************************/
283:       /* Compute objective function.                                         */
284:       /***********************************************************************/

286:       KSP_MatMult(ksp, Qmat, d, z);
287:       VecAYPX(z, -0.5, ksp->vec_rhs);
288:       VecDot(d, z, &cg->o_fcn);
289:       cg->o_fcn = -cg->o_fcn;
290:       ++ksp->its;
291:     }
292:     return(0);
293:   }

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

302:   dMp    = 0.0;
303:   norm_d = 0.0;
304:   switch (cg->dtype) {
305:   case GLTR_PRECONDITIONED_DIRECTION:
306:     norm_p = rz;
307:     break;

309:   default:
310:     VecDot(p, p, &norm_p);
311:     break;
312:   }

314:   cg->diag[t_size] = kappa / rz;
315:   cg->offd[t_size] = 0.0;
316:   ++t_size;

318:   piv = 1.0;

320:   /***************************************************************************/
321:   /* Check for breakdown of the conjugate gradient method; this occurs when  */
322:   /* kappa is zero.                                                          */
323:   /***************************************************************************/

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

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

334:     if (cg->radius && norm_p > 0.0) {
335:       /***********************************************************************/
336:       /* Follow direction of negative curvature to the boundary of the       */
337:       /* trust region.                                                       */
338:       /***********************************************************************/

340:       step       = PetscSqrtReal(r2 / norm_p);
341:       cg->norm_d = cg->radius;

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

345:       /***********************************************************************/
346:       /* Update objective function.                                          */
347:       /***********************************************************************/

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

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

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

366:       /***********************************************************************/
367:       /* Compute objective function.                                         */
368:       /***********************************************************************/

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

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

386:   while (1) {
387:     /*************************************************************************/
388:     /* Know that kappa is nonzero, because we have not broken down, so we    */
389:     /* can compute the steplength.                                           */
390:     /*************************************************************************/

392:     alpha             = rz / kappa;
393:     cg->alpha[l_size] = alpha;

395:     /*************************************************************************/
396:     /* Compute the diagonal value of the Cholesky factorization of the       */
397:     /* Lanczos matrix and check to see if the Lanczos matrix is indefinite.  */
398:     /* This indicates a direction of negative curvature.                     */
399:     /*************************************************************************/

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

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

412:       if (cg->radius && norm_p > 0.0) {
413:         /*********************************************************************/
414:         /* Follow direction of negative curvature to boundary.               */
415:         /*********************************************************************/

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

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

422:         /*********************************************************************/
423:         /* Update objective function.                                        */
424:         /*********************************************************************/

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

435:     /*************************************************************************/
436:     /* Compute the steplength and check for intersection with the trust      */
437:     /* region.                                                               */
438:     /*************************************************************************/

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

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

450:       if (norm_p > 0.0) {
451:         /*********************************************************************/
452:         /* Follow the direction to the boundary of the trust region.         */
453:         /*********************************************************************/

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

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

460:         /*********************************************************************/
461:         /* Update objective function.                                        */
462:         /*********************************************************************/

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

473:     /*************************************************************************/
474:     /* Now we can update the direction and residual.                         */
475:     /*************************************************************************/

477:     VecAXPY(d, alpha, p);          /* d = d + alpha p   */
478:     VecAXPY(r, -alpha, z);         /* r = r - alpha Q p */
479:     KSP_PCApply(ksp, r, z);        /* z = inv(M) r      */

481:     switch (cg->dtype) {
482:     case GLTR_PRECONDITIONED_DIRECTION:
483:       norm_d = norm_dp1;
484:       break;

486:     default:
487:       VecDot(d, d, &norm_d);
488:       break;
489:     }
490:     cg->norm_d = PetscSqrtReal(norm_d);

492:     /*************************************************************************/
493:     /* Update objective function.                                            */
494:     /*************************************************************************/

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

498:     /*************************************************************************/
499:     /* Check that the preconditioner appears positive semidefinite.          */
500:     /*************************************************************************/

502:     rzm1 = rz;
503:     VecDot(r, z, &rz);             /* rz = r^T z        */
504:     if (rz < 0.0) {
505:       /***********************************************************************/
506:       /* The preconditioner is indefinite.                                   */
507:       /***********************************************************************/

509:       ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
510:       PetscInfo1(ksp, "KSPCGSolve_GLTR: cg indefinite preconditioner: rz=%g\n", (double)rz);
511:       break;
512:     }

514:     /*************************************************************************/
515:     /* As far as we know, the preconditioner is positive semidefinite.       */
516:     /* Compute the residual and check for convergence.                       */
517:     /*************************************************************************/

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

521:     switch (ksp->normtype) {
522:     case KSP_NORM_PRECONDITIONED:
523:       VecNorm(z, NORM_2, &norm_r); /* norm_r = |z|      */
524:       break;

526:     case KSP_NORM_UNPRECONDITIONED:
527:       VecNorm(r, NORM_2, &norm_r); /* norm_r = |r|      */
528:       break;

530:     case KSP_NORM_NATURAL:
531:       norm_r = cg->norm_r[l_size+1];                    /* norm_r = |r|_M    */
532:       break;

534:     default:
535:       norm_r = 0.0;
536:       break;
537:     }

539:     KSPLogResidualHistory(ksp, norm_r);
540:     KSPMonitor(ksp, ksp->its, norm_r);
541:     ksp->rnorm = norm_r;

543:     (*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP);
544:     if (ksp->reason) {
545:       /***********************************************************************/
546:       /* The method has converged.                                           */
547:       /***********************************************************************/

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

553:     /*************************************************************************/
554:     /* We have not converged yet.  Check for breakdown.                      */
555:     /*************************************************************************/

557:     beta = rz / rzm1;
558:     if (PetscAbsReal(beta) <= 0.0) {
559:       /***********************************************************************/
560:       /* Conjugate gradients has broken down.                                */
561:       /***********************************************************************/

563:       ksp->reason = KSP_DIVERGED_BREAKDOWN;
564:       PetscInfo1(ksp, "KSPCGSolve_GLTR: breakdown: beta=%g\n", (double)beta);
565:       break;
566:     }

568:     /*************************************************************************/
569:     /* Check iteration limit.                                                */
570:     /*************************************************************************/

572:     if (ksp->its >= max_cg_its) {
573:       ksp->reason = KSP_DIVERGED_ITS;
574:       PetscInfo1(ksp, "KSPCGSolve_GLTR: iterlim: its=%D\n", ksp->its);
575:       break;
576:     }

578:     /*************************************************************************/
579:     /* Update p and the norms.                                               */
580:     /*************************************************************************/

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

585:     switch (cg->dtype) {
586:     case GLTR_PRECONDITIONED_DIRECTION:
587:       dMp    = beta*(dMp + alpha*norm_p);
588:       norm_p = beta*(rzm1 + beta*norm_p);
589:       break;

591:     default:
592:       VecDot(d, p, &dMp);
593:       VecDot(p, p, &norm_p);
594:       break;
595:     }

597:     /*************************************************************************/
598:     /* Compute the new direction and update the iteration.                   */
599:     /*************************************************************************/

601:     KSP_MatMult(ksp, Qmat, p, z);  /* z = Q * p         */
602:     VecDot(p, z, &kappa);          /* kappa = p^T Q p   */
603:     ++ksp->its;

605:     /*************************************************************************/
606:     /* Update the Lanczos tridiagonal matrix.                            */
607:     /*************************************************************************/

609:     ++l_size;
610:     cg->offd[t_size] = PetscSqrtReal(beta) / PetscAbsReal(alpha);
611:     cg->diag[t_size] = kappa / rz + beta / alpha;
612:     ++t_size;

614:     /*************************************************************************/
615:     /* Check for breakdown of the conjugate gradient method; this occurs     */
616:     /* when kappa is zero.                                                   */
617:     /*************************************************************************/

619:     if (PetscAbsReal(kappa) <= 0.0) {
620:       /***********************************************************************/
621:       /* The method breaks down; move along the direction as if the matrix   */
622:       /* were indefinite.                                                    */
623:       /***********************************************************************/

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

628:       if (cg->radius && norm_p > 0.0) {
629:         /*********************************************************************/
630:         /* Follow direction to boundary.                                     */
631:         /*********************************************************************/

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

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

638:         /*********************************************************************/
639:         /* Update objective function.                                        */
640:         /*********************************************************************/

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

652:   /***************************************************************************/
653:   /* Check to see if we need to continue with the Lanczos method.            */
654:   /***************************************************************************/

656:   if (!cg->radius) {
657:     /*************************************************************************/
658:     /* There is no radius.  Therefore, we cannot move along the boundary.    */
659:     /*************************************************************************/
660:     return(0);
661:   }

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

672:   /***************************************************************************/
673:   /* Switch to contructing the Lanczos basis by way of the conjugate         */
674:   /* directions.                                                             */
675:   /***************************************************************************/

677:   for (i = 0; i < max_lanczos_its; ++i) {
678:     /*************************************************************************/
679:     /* Check for breakdown of the conjugate gradient method; this occurs     */
680:     /* when kappa is zero.                                                   */
681:     /*************************************************************************/

683:     if (PetscAbsReal(kappa) <= 0.0) {
684:       ksp->reason = KSP_DIVERGED_BREAKDOWN;
685:       PetscInfo1(ksp, "KSPCGSolve_GLTR: lanczos breakdown: kappa=%g\n", (double)kappa);
686:       break;
687:     }

689:     /*************************************************************************/
690:     /* Update the direction and residual.                                    */
691:     /*************************************************************************/

693:     alpha             = rz / kappa;
694:     cg->alpha[l_size] = alpha;

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

699:     /*************************************************************************/
700:     /* Check that the preconditioner appears positive semidefinite.          */
701:     /*************************************************************************/

703:     rzm1 = rz;
704:     VecDot(r, z, &rz);             /* rz = r^T z        */
705:     if (rz < 0.0) {
706:       /***********************************************************************/
707:       /* The preconditioner is indefinite.                                   */
708:       /***********************************************************************/

710:       ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
711:       PetscInfo1(ksp, "KSPCGSolve_GLTR: lanczos indefinite preconditioner: rz=%g\n", (double)rz);
712:       break;
713:     }

715:     /*************************************************************************/
716:     /* As far as we know, the preconditioner is positive definite.  Compute  */
717:     /* the residual.  Do NOT check for convergence.                          */
718:     /*************************************************************************/

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

722:     switch (ksp->normtype) {
723:     case KSP_NORM_PRECONDITIONED:
724:       VecNorm(z, NORM_2, &norm_r); /* norm_r = |z|      */
725:       break;

727:     case KSP_NORM_UNPRECONDITIONED:
728:       VecNorm(r, NORM_2, &norm_r); /* norm_r = |r|      */
729:       break;

731:     case KSP_NORM_NATURAL:
732:       norm_r = cg->norm_r[l_size+1];                    /* norm_r = |r|_M    */
733:       break;

735:     default:
736:       norm_r = 0.0;
737:       break;
738:     }

740:     KSPLogResidualHistory(ksp, norm_r);
741:     KSPMonitor(ksp, ksp->its, norm_r);
742:     ksp->rnorm = norm_r;

744:     /*************************************************************************/
745:     /* Check for breakdown.                                                  */
746:     /*************************************************************************/

748:     beta = rz / rzm1;
749:     if (PetscAbsReal(beta) <= 0.0) {
750:       /***********************************************************************/
751:       /* Conjugate gradients has broken down.                                */
752:       /***********************************************************************/

754:       ksp->reason = KSP_DIVERGED_BREAKDOWN;
755:       PetscInfo1(ksp, "KSPCGSolve_GLTR: breakdown: beta=%g\n",(double) beta);
756:       break;
757:     }

759:     /*************************************************************************/
760:     /* Update p and the norms.                                               */
761:     /*************************************************************************/

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

766:     /*************************************************************************/
767:     /* Compute the new direction and update the iteration.                   */
768:     /*************************************************************************/

770:     KSP_MatMult(ksp, Qmat, p, z);  /* z = Q * p         */
771:     VecDot(p, z, &kappa);          /* kappa = p^T Q p   */
772:     ++ksp->its;

774:     /*************************************************************************/
775:     /* Update the Lanczos tridiagonal matrix.                                */
776:     /*************************************************************************/

778:     ++l_size;
779:     cg->offd[t_size] = PetscSqrtReal(beta) / PetscAbsReal(alpha);
780:     cg->diag[t_size] = kappa / rz + beta / alpha;
781:     ++t_size;
782:   }

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

795:   if (t_size > cg->alloced) {
796:     if (cg->alloced) {
797:       PetscFree(cg->rwork);
798:       PetscFree(cg->iwork);
799:       cg->alloced += cg->init_alloc;
800:     } else {
801:       cg->alloced = cg->init_alloc;
802:     }

804:     while (t_size > cg->alloced) {
805:       cg->alloced += cg->init_alloc;
806:     }

808:     cg->alloced = PetscMin(cg->alloced, t_size);
809:     PetscMalloc2(10*cg->alloced, &cg->rwork,5*cg->alloced, &cg->iwork);
810:   }

812:   /***************************************************************************/
813:   /* Set up the required vectors.                                            */
814:   /***************************************************************************/

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

823:   e_iblk = cg->iwork + 0*t_size;                        /* Eigen blocks      */
824:   e_splt = cg->iwork + 1*t_size;                        /* Eigen splits      */
825:   e_iwrk = cg->iwork + 2*t_size;                        /* Eigen workspace   */

827:   /***************************************************************************/
828:   /* Compute the minimum eigenvalue of T.                                    */
829:   /***************************************************************************/

831:   vl = 0.0;
832:   vu = 0.0;
833:   il = 1;
834:   iu = 1;

836: #if defined(PETSC_MISSING_LAPACK_STEBZ)
837:   SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"STEBZ - Lapack routine is unavailable.");
838: #else
839:   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));

841:   if ((0 != info) || (1 != e_valus)) {
842:     /*************************************************************************/
843:     /* Calculation of the minimum eigenvalue failed.  Return the             */
844:     /* Steihaug-Toint direction.                                             */
845:     /*************************************************************************/

847:     PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute eigenvalue.\n");
848:     ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
849:     return(0);
850:   }

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

854:   /***************************************************************************/
855:   /* Compute the initial value of lambda to make (T + lamba I) positive      */
856:   /* definite.                                                               */
857:   /***************************************************************************/

859:   pert = cg->init_pert;
860:   if (e_valu[0] < 0.0) cg->lambda = pert - e_valu[0];
861: #endif

863:   while (1) {
864:     for (i = 0; i < t_size; ++i) {
865:       t_diag[i] = cg->diag[i] + cg->lambda;
866:       t_offd[i] = cg->offd[i];
867:     }

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

874:     if (0 == info) break;

876:     pert      += pert;
877:     cg->lambda = cg->lambda * (1.0 + pert) + pert;
878: #endif
879:   }

881:   /***************************************************************************/
882:   /* Compute the initial step and its norm.                                  */
883:   /***************************************************************************/

885:   nrhs = 1;
886:   nldb = t_size;

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

891: #if defined(PETSC_MISSING_LAPACK_PTTRS)
892:   SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"PTTRS - Lapack routine is unavailable.");
893: #else
894:   PetscStackCallBLAS("LAPACKpttrs",LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, t_soln, &nldb, &info));
895: #endif

897:   if (0 != info) {
898:     /*************************************************************************/
899:     /* Calculation of the initial step failed; return the Steihaug-Toint     */
900:     /* direction.                                                            */
901:     /*************************************************************************/

903:     PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute step.\n");
904:     ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
905:     return(0);
906:   }

908:   norm_t = 0.;
909:   for (i = 0; i < t_size; ++i) norm_t += t_soln[i] * t_soln[i];
910:   norm_t = PetscSqrtReal(norm_t);

912:   /***************************************************************************/
913:   /* Determine the case we are in.                                           */
914:   /***************************************************************************/

916:   if (norm_t <= cg->radius) {
917:     /*************************************************************************/
918:     /* The step is within the trust region; check if we are in the hard case */
919:     /* and need to move to the boundary by following a direction of negative */
920:     /* curvature.                                                            */
921:     /*************************************************************************/

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

929: #if defined(PETSC_MISSING_LAPACK_STEIN)
930:       SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"STEIN - Lapack routine is unavailable.");
931: #else
932:       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));
933: #endif

935:       if (0 != info) {
936:         /*********************************************************************/
937:         /* Calculation of the minimum eigenvalue failed.  Return the         */
938:         /* Steihaug-Toint direction.                                         */
939:         /*********************************************************************/

941:         PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute eigenvector.\n");
942:         ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
943:         return(0);
944:       }

946:       coef1 = 0.0;
947:       coef2 = 0.0;
948:       coef3 = -cg->radius * cg->radius;
949:       for (i = 0; i < t_size; ++i) {
950:         coef1 += e_vect[i] * e_vect[i];
951:         coef2 += e_vect[i] * t_soln[i];
952:         coef3 += t_soln[i] * t_soln[i];
953:       }

955:       coef3 = PetscSqrtReal(coef2 * coef2 - coef1 * coef3);
956:       root1 = (-coef2 + coef3) / coef1;
957:       root2 = (-coef2 - coef3) / coef1;

959:       /***********************************************************************/
960:       /* Compute objective value for (t_soln + root1 * e_vect)               */
961:       /***********************************************************************/

963:       for (i = 0; i < t_size; ++i) {
964:         e_rwrk[i] = t_soln[i] + root1 * e_vect[i];
965:       }

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

977:       /***********************************************************************/
978:       /* Compute objective value for (t_soln + root2 * e_vect)               */
979:       /***********************************************************************/

981:       for (i = 0; i < t_size; ++i) {
982:         e_rwrk[i] = t_soln[i] + root2 * e_vect[i];
983:       }

985:       obj2 = e_rwrk[0]*(0.5*(cg->diag[0]*e_rwrk[0]+
986:                              cg->offd[1]*e_rwrk[1])+cg->norm_r[0]);
987:       for (i = 1; i < t_size - 1; ++i) {
988:         obj2 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
989:                                cg->diag[i]*e_rwrk[i]+
990:                                cg->offd[i+1]*e_rwrk[i+1]);
991:       }
992:       obj2 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
993:                              cg->diag[i]*e_rwrk[i]);

995:       /***********************************************************************/
996:       /* Choose the point with the best objective function value.            */
997:       /***********************************************************************/

999:       if (obj1 <= obj2) {
1000:         for (i = 0; i < t_size; ++i) {
1001:           t_soln[i] += root1 * e_vect[i];
1002:         }
1003:       }
1004:       else {
1005:         for (i = 0; i < t_size; ++i) {
1006:           t_soln[i] += root2 * e_vect[i];
1007:         }
1008:       }
1009:     } else {
1010:       /***********************************************************************/
1011:       /* The matrix is positive definite or there was no room to move; the   */
1012:       /* solution is already contained in t_soln.                            */
1013:       /***********************************************************************/
1014:     }
1015:   } else {
1016:     /*************************************************************************/
1017:     /* The step is outside the trust-region.  Compute the correct value for  */
1018:     /* lambda by performing Newton's method.                                 */
1019:     /*************************************************************************/

1021:     for (i = 0; i < max_newton_its; ++i) {
1022:       /***********************************************************************/
1023:       /* Check for convergence.                                              */
1024:       /***********************************************************************/

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

1028:       /***********************************************************************/
1029:       /* Compute the update.                                                 */
1030:       /***********************************************************************/

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

1034: #if defined(PETSC_MISSING_LAPACK_PTTRS)
1035:       SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"PTTRS - Lapack routine is unavailable.");
1036: #else
1037:       PetscStackCallBLAS("LAPACKpttrs",LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, e_rwrk, &nldb, &info));
1038: #endif

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

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

1051:       /***********************************************************************/
1052:       /* Modify lambda.                                                      */
1053:       /***********************************************************************/

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

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

1060:       /***********************************************************************/
1061:       /* Factor T + lambda I                                                 */
1062:       /***********************************************************************/

1064:       for (j = 0; j < t_size; ++j) {
1065:         t_diag[j] = cg->diag[j] + cg->lambda;
1066:         t_offd[j] = cg->offd[j];
1067:       }

1069: #if defined(PETSC_MISSING_LAPACK_PTTRF)
1070:       SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"PTTRF - Lapack routine is unavailable.");
1071: #else
1072:       PetscStackCallBLAS("LAPACKpttrf",LAPACKpttrf_(&t_size, t_diag, t_offd + 1, &info));
1073: #endif

1075:       if (0 != info) {
1076:         /*********************************************************************/
1077:         /* Calculation of factorization failed; return the Steihaug-Toint    */
1078:         /* direction.                                                        */
1079:         /*********************************************************************/

1081:         PetscInfo(ksp, "KSPCGSolve_GLTR: factorization failed.\n");
1082:         ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1083:         return(0);
1084:       }

1086:       /***********************************************************************/
1087:       /* Compute the new step and its norm.                                  */
1088:       /***********************************************************************/

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

1093: #if defined(PETSC_MISSING_LAPACK_PTTRS)
1094:       SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"PTTRS - Lapack routine is unavailable.");
1095: #else
1096:       PetscStackCallBLAS("LAPACKpttrs",LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, t_soln, &nldb, &info));
1097: #endif

1099:       if (0 != info) {
1100:         /*********************************************************************/
1101:         /* Calculation of the step failed; return the Steihaug-Toint         */
1102:         /* direction.                                                        */
1103:         /*********************************************************************/

1105:         PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute step.\n");
1106:         ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1107:         return(0);
1108:       }

1110:       norm_t = 0.;
1111:       for (j = 0; j < t_size; ++j) norm_t += t_soln[j] * t_soln[j];
1112:       norm_t = PetscSqrtReal(norm_t);
1113:     }

1115:     /*************************************************************************/
1116:     /* Check for convergence.                                                */
1117:     /*************************************************************************/

1119:     if (PetscAbsReal(norm_t - cg->radius) > cg->newton_tol * cg->radius) {
1120:       /***********************************************************************/
1121:       /* Newton method failed to converge in iteration limit.                */
1122:       /***********************************************************************/

1124:       PetscInfo(ksp, "KSPCGSolve_GLTR: failed to converge.\n");
1125:       ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1126:       return(0);
1127:     }
1128:   }

1130:   /***************************************************************************/
1131:   /* Recover the norm of the direction and objective function value.         */
1132:   /***************************************************************************/

1134:   cg->norm_d = norm_t;

1136:   cg->o_fcn = t_soln[0]*(0.5*(cg->diag[0]*t_soln[0]+cg->offd[1]*t_soln[1])+cg->norm_r[0]);
1137:   for (i = 1; i < t_size - 1; ++i) {
1138:     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]);
1139:   }
1140:   cg->o_fcn += 0.5*t_soln[i]*(cg->offd[i]*t_soln[i-1]+cg->diag[i]*t_soln[i]);

1142:   /***************************************************************************/
1143:   /* Recover the direction.                                                  */
1144:   /***************************************************************************/

1146:   sigma = -1;

1148:   /***************************************************************************/
1149:   /* Start conjugate gradient method from the beginning                      */
1150:   /***************************************************************************/

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

1155:   /***************************************************************************/
1156:   /* Accumulate Q * s                                                        */
1157:   /***************************************************************************/

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

1162:   /***************************************************************************/
1163:   /* Compute the first direction.                                            */
1164:   /***************************************************************************/

1166:   VecCopy(z, p);                   /* p = z             */
1167:   KSP_MatMult(ksp, Qmat, p, z);    /* z = Q * p         */
1168:   ++ksp->its;

1170:   for (i = 0; i < l_size - 1; ++i) {
1171:     /*************************************************************************/
1172:     /* Update the residual and direction.                                    */
1173:     /*************************************************************************/

1175:     alpha = cg->alpha[i];
1176:     if (alpha >= 0.0) sigma = -sigma;

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

1181:     /*************************************************************************/
1182:     /* Accumulate Q * s                                                      */
1183:     /*************************************************************************/

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

1187:     /*************************************************************************/
1188:     /* Update p.                                                             */
1189:     /*************************************************************************/

1191:     beta = cg->beta[i];
1192:     VecAYPX(p, beta, z);          /* p = z + beta p    */
1193:     KSP_MatMult(ksp, Qmat, p, z);  /* z = Q * p         */
1194:     ++ksp->its;
1195:   }

1197:   /***************************************************************************/
1198:   /* Update the residual and direction.                                      */
1199:   /***************************************************************************/

1201:   alpha = cg->alpha[i];
1202:   if (alpha >= 0.0) sigma = -sigma;

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

1207:   /***************************************************************************/
1208:   /* Accumulate Q * s                                                        */
1209:   /***************************************************************************/

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

1213:   /***************************************************************************/
1214:   /* Set the termination reason.                                             */
1215:   /***************************************************************************/

1217:   ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1218:   return(0);
1219: #endif
1220: }

1222: static PetscErrorCode KSPCGSetUp_GLTR(KSP ksp)
1223: {
1224:   KSPCG_GLTR     *cg = (KSPCG_GLTR*)ksp->data;
1225:   PetscInt       max_its;

1229:   /***************************************************************************/
1230:   /* Determine the total maximum number of iterations.                       */
1231:   /***************************************************************************/

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

1235:   /***************************************************************************/
1236:   /* Set work vectors needed by conjugate gradient method and allocate       */
1237:   /* workspace for Lanczos matrix.                                           */
1238:   /***************************************************************************/

1240:   KSPSetWorkVecs(ksp, 3);

1242:   PetscCalloc5(max_its,&cg->diag,max_its,&cg->offd,max_its,&cg->alpha,max_its,&cg->beta,max_its,&cg->norm_r);
1243:   PetscLogObjectMemory((PetscObject)ksp, 5*max_its*sizeof(PetscReal));
1244:   return(0);
1245: }

1247: static PetscErrorCode KSPCGDestroy_GLTR(KSP ksp)
1248: {
1249:   KSPCG_GLTR     *cg = (KSPCG_GLTR*)ksp->data;

1253:   /***************************************************************************/
1254:   /* Free memory allocated for the data.                                     */
1255:   /***************************************************************************/

1257:   PetscFree5(cg->diag,cg->offd,cg->alpha,cg->beta,cg->norm_r);
1258:   if (cg->alloced) {
1259:     PetscFree2(cg->rwork,cg->iwork);
1260:   }

1262:   /***************************************************************************/
1263:   /* Clear composed functions                                                */
1264:   /***************************************************************************/

1266:   PetscObjectComposeFunction((PetscObject)ksp,"KSPCGSetRadius_C",NULL);
1267:   PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetNormD_C",NULL);
1268:   PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetObjFcn_C",NULL);
1269:   PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGLTRGetMinEig_C",NULL);
1270:   PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGLTRGetLambda_C",NULL);

1272:   /***************************************************************************/
1273:   /* Destroy KSP object.                                                     */
1274:   /***************************************************************************/
1275:   KSPDestroyDefault(ksp);
1276:   return(0);
1277: }

1279: static PetscErrorCode  KSPCGSetRadius_GLTR(KSP ksp, PetscReal radius)
1280: {
1281:   KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;

1284:   cg->radius = radius;
1285:   return(0);
1286: }

1288: static PetscErrorCode  KSPCGGetNormD_GLTR(KSP ksp, PetscReal *norm_d)
1289: {
1290:   KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;

1293:   *norm_d = cg->norm_d;
1294:   return(0);
1295: }

1297: static PetscErrorCode  KSPCGGetObjFcn_GLTR(KSP ksp, PetscReal *o_fcn)
1298: {
1299:   KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;

1302:   *o_fcn = cg->o_fcn;
1303:   return(0);
1304: }

1306: static PetscErrorCode  KSPCGGLTRGetMinEig_GLTR(KSP ksp, PetscReal *e_min)
1307: {
1308:   KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;

1311:   *e_min = cg->e_min;
1312:   return(0);
1313: }

1315: static PetscErrorCode  KSPCGGLTRGetLambda_GLTR(KSP ksp, PetscReal *lambda)
1316: {
1317:   KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;

1320:   *lambda = cg->lambda;
1321:   return(0);
1322: }

1324: static PetscErrorCode KSPCGSetFromOptions_GLTR(PetscOptionItems *PetscOptionsObject,KSP ksp)
1325: {
1327:   KSPCG_GLTR       *cg = (KSPCG_GLTR*)ksp->data;

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

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

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

1336:   PetscOptionsReal("-ksp_cg_gltr_init_pert", "Initial perturbation", "", cg->init_pert, &cg->init_pert, NULL);
1337:   PetscOptionsReal("-ksp_cg_gltr_eigen_tol", "Eigenvalue tolerance", "", cg->eigen_tol, &cg->eigen_tol, NULL);
1338:   PetscOptionsReal("-ksp_cg_gltr_newton_tol", "Newton tolerance", "", cg->newton_tol, &cg->newton_tol, NULL);

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

1343:   PetscOptionsTail();
1344:   return(0);
1345: }

1347: /*MC
1348:      KSPCGGLTR -   Code to run conjugate gradient method subject to a constraint
1349:          on the solution norm. This is used in Trust Region methods for
1350:          nonlinear equations, SNESNEWTONTR

1352:    Options Database Keys:
1353: .      -ksp_cg_radius <r> - Trust Region Radius

1355:    Notes: This is rarely used directly

1357:   Use preconditioned conjugate gradient to compute
1358:   an approximate minimizer of the quadratic function

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

1362:    subject to the trust region constraint

1364:             || s || <= delta,

1366:    where

1368:      delta is the trust region radius,
1369:      g is the gradient vector,
1370:      H is the Hessian approximation,
1371:      M is the positive definite preconditioner matrix.

1373:    KSPConvergedReason may be
1374: $  KSP_CONVERGED_CG_NEG_CURVE if convergence is reached along a negative curvature direction,
1375: $  KSP_CONVERGED_CG_CONSTRAINED if convergence is reached along a constrained step,
1376: $  other KSP converged/diverged reasons

1378:   Notes:
1379:   The preconditioner supplied should be symmetric and positive definite.

1381:    Level: developer

1383: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPCGSetRadius(), KSPCGGetNormD(), KSPCGGetObjFcn(), KSPCGGLTRGetMinEig(), KSPCGGLTRGetLambda()
1384: M*/

1386: PETSC_EXTERN PetscErrorCode KSPCreate_CGGLTR(KSP ksp)
1387: {
1389:   KSPCG_GLTR       *cg;

1392:   PetscNewLog(ksp,&cg);
1393:   cg->radius = 0.0;
1394:   cg->dtype  = GLTR_UNPRECONDITIONED_DIRECTION;

1396:   cg->init_pert  = 1.0e-8;
1397:   cg->eigen_tol  = 1.0e-10;
1398:   cg->newton_tol = 1.0e-6;

1400:   cg->alloced    = 0;
1401:   cg->init_alloc = 1024;

1403:   cg->max_lanczos_its = 20;
1404:   cg->max_newton_its  = 10;

1406:   ksp->data = (void*) cg;
1407:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,3);
1408:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);
1409:   KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,2);

1411:   /***************************************************************************/
1412:   /* Sets the functions that are associated with this data structure         */
1413:   /* (in C++ this is the same as defining virtual functions).                */
1414:   /***************************************************************************/

1416:   ksp->ops->setup          = KSPCGSetUp_GLTR;
1417:   ksp->ops->solve          = KSPCGSolve_GLTR;
1418:   ksp->ops->destroy        = KSPCGDestroy_GLTR;
1419:   ksp->ops->setfromoptions = KSPCGSetFromOptions_GLTR;
1420:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
1421:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
1422:   ksp->ops->view           = 0;

1424:   PetscObjectComposeFunction((PetscObject)ksp,"KSPCGSetRadius_C",KSPCGSetRadius_GLTR);
1425:   PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetNormD_C", KSPCGGetNormD_GLTR);
1426:   PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetObjFcn_C",KSPCGGetObjFcn_GLTR);
1427:   PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGLTRGetMinEig_C",KSPCGGLTRGetMinEig_GLTR);
1428:   PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGLTRGetLambda_C",KSPCGGLTRGetLambda_GLTR);
1429:   return(0);
1430: }