Actual source code: gltr.c

petsc-master 2019-06-22
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: #if !defined(PETSC_MISSING_LAPACK_STEBZ)
 81:   PetscBLASInt e_valus, e_splts;
 82: #endif
 83:   PetscBool diagonalscale;

 86:   /***************************************************************************/
 87:   /* Check the arguments and parameters.                                     */
 88:   /***************************************************************************/

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

 94:   /***************************************************************************/
 95:   /* Get the workspace vectors and initialize variables                      */
 96:   /***************************************************************************/

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

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

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

113:   /***************************************************************************/
114:   /* Initialize objective function direction, and minimum eigenvalue.        */
115:   /***************************************************************************/

117:   cg->o_fcn = 0.0;

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

122:   cg->e_min  = 0.0;
123:   cg->lambda = 0.0;

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

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

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

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

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

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

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

165:       /***********************************************************************/
166:       /* Compute objective function.                                         */
167:       /***********************************************************************/

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

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

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

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

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

200:       /***********************************************************************/
201:       /* Compute objective function.                                         */
202:       /***********************************************************************/

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

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

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

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

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

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

235:   default:
236:     norm_r = 0.0;
237:     break;
238:   }

240:   KSPLogResidualHistory(ksp, norm_r);
241:   KSPMonitor(ksp, ksp->its, norm_r);
242:   ksp->rnorm = norm_r;

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

246:   /***************************************************************************/
247:   /* Compute the first direction and update the iteration.                   */
248:   /***************************************************************************/

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

254:   /***************************************************************************/
255:   /* Check the matrix for numerical problems.                                */
256:   /***************************************************************************/

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

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

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

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

280:       /***********************************************************************/
281:       /* Compute objective function.                                         */
282:       /***********************************************************************/

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

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

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

307:   default:
308:     VecDot(p, p, &norm_p);
309:     break;
310:   }

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

316:   piv = 1.0;

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

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

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

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

338:       step       = PetscSqrtReal(r2 / norm_p);
339:       cg->norm_d = cg->radius;

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

343:       /***********************************************************************/
344:       /* Update objective function.                                          */
345:       /***********************************************************************/

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

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

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

364:       /***********************************************************************/
365:       /* Compute objective function.                                         */
366:       /***********************************************************************/

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

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

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

390:     alpha             = rz / kappa;
391:     cg->alpha[l_size] = alpha;

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

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

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

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

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

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

420:         /*********************************************************************/
421:         /* Update objective function.                                        */
422:         /*********************************************************************/

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

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

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

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

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

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

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

458:         /*********************************************************************/
459:         /* Update objective function.                                        */
460:         /*********************************************************************/

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

471:     /*************************************************************************/
472:     /* Now we can update the direction and residual.                         */
473:     /*************************************************************************/

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

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

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

490:     /*************************************************************************/
491:     /* Update objective function.                                            */
492:     /*************************************************************************/

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

496:     /*************************************************************************/
497:     /* Check that the preconditioner appears positive semidefinite.          */
498:     /*************************************************************************/

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

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

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

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

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

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

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

532:     default:
533:       norm_r = 0.0;
534:       break;
535:     }

537:     KSPLogResidualHistory(ksp, norm_r);
538:     KSPMonitor(ksp, ksp->its, norm_r);
539:     ksp->rnorm = norm_r;

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

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

551:     /*************************************************************************/
552:     /* We have not converged yet.  Check for breakdown.                      */
553:     /*************************************************************************/

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

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

566:     /*************************************************************************/
567:     /* Check iteration limit.                                                */
568:     /*************************************************************************/

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

576:     /*************************************************************************/
577:     /* Update p and the norms.                                               */
578:     /*************************************************************************/

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

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

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

595:     /*************************************************************************/
596:     /* Compute the new direction and update the iteration.                   */
597:     /*************************************************************************/

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

603:     /*************************************************************************/
604:     /* Update the Lanczos tridiagonal matrix.                            */
605:     /*************************************************************************/

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

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

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

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

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

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

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

636:         /*********************************************************************/
637:         /* Update objective function.                                        */
638:         /*********************************************************************/

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

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

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

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

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

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

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

687:     /*************************************************************************/
688:     /* Update the direction and residual.                                    */
689:     /*************************************************************************/

691:     alpha             = rz / kappa;
692:     cg->alpha[l_size] = alpha;

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

697:     /*************************************************************************/
698:     /* Check that the preconditioner appears positive semidefinite.          */
699:     /*************************************************************************/

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

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

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

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

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

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

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

733:     default:
734:       norm_r = 0.0;
735:       break;
736:     }

738:     KSPLogResidualHistory(ksp, norm_r);
739:     KSPMonitor(ksp, ksp->its, norm_r);
740:     ksp->rnorm = norm_r;

742:     /*************************************************************************/
743:     /* Check for breakdown.                                                  */
744:     /*************************************************************************/

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

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

757:     /*************************************************************************/
758:     /* Update p and the norms.                                               */
759:     /*************************************************************************/

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

764:     /*************************************************************************/
765:     /* Compute the new direction and update the iteration.                   */
766:     /*************************************************************************/

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

772:     /*************************************************************************/
773:     /* Update the Lanczos tridiagonal matrix.                                */
774:     /*************************************************************************/

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

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

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

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

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

810:   /***************************************************************************/
811:   /* Set up the required vectors.                                            */
812:   /***************************************************************************/

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

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

825:   /***************************************************************************/
826:   /* Compute the minimum eigenvalue of T.                                    */
827:   /***************************************************************************/

829:   vl = 0.0;
830:   vu = 0.0;
831:   il = 1;
832:   iu = 1;

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

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

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

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

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

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

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

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

872:     if (0 == info) break;

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

879:   /***************************************************************************/
880:   /* Compute the initial step and its norm.                                  */
881:   /***************************************************************************/

883:   nrhs = 1;
884:   nldb = t_size;

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

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

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

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

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

910:   /***************************************************************************/
911:   /* Determine the case we are in.                                           */
912:   /***************************************************************************/

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

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

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

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

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

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

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

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

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

965:       obj1 = 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:         obj1 += 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:       obj1 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
973:                              cg->diag[i]*e_rwrk[i]);

975:       /***********************************************************************/
976:       /* Compute objective value for (t_soln + root2 * e_vect)               */
977:       /***********************************************************************/

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

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

993:       /***********************************************************************/
994:       /* Choose the point with the best objective function value.            */
995:       /***********************************************************************/

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

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

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

1026:       /***********************************************************************/
1027:       /* Compute the update.                                                 */
1028:       /***********************************************************************/

1030:       PetscArraycpy(e_rwrk, t_soln, t_size);

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

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

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

1049:       /***********************************************************************/
1050:       /* Modify lambda.                                                      */
1051:       /***********************************************************************/

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

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

1058:       /***********************************************************************/
1059:       /* Factor T + lambda I                                                 */
1060:       /***********************************************************************/

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

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

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

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

1084:       /***********************************************************************/
1085:       /* Compute the new step and its norm.                                  */
1086:       /***********************************************************************/

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

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

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

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

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

1113:     /*************************************************************************/
1114:     /* Check for convergence.                                                */
1115:     /*************************************************************************/

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

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

1128:   /***************************************************************************/
1129:   /* Recover the norm of the direction and objective function value.         */
1130:   /***************************************************************************/

1132:   cg->norm_d = norm_t;

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

1140:   /***************************************************************************/
1141:   /* Recover the direction.                                                  */
1142:   /***************************************************************************/

1144:   sigma = -1;

1146:   /***************************************************************************/
1147:   /* Start conjugate gradient method from the beginning                      */
1148:   /***************************************************************************/

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

1153:   /***************************************************************************/
1154:   /* Accumulate Q * s                                                        */
1155:   /***************************************************************************/

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

1160:   /***************************************************************************/
1161:   /* Compute the first direction.                                            */
1162:   /***************************************************************************/

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

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

1173:     alpha = cg->alpha[i];
1174:     if (alpha >= 0.0) sigma = -sigma;

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

1179:     /*************************************************************************/
1180:     /* Accumulate Q * s                                                      */
1181:     /*************************************************************************/

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

1185:     /*************************************************************************/
1186:     /* Update p.                                                             */
1187:     /*************************************************************************/

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

1195:   /***************************************************************************/
1196:   /* Update the residual and direction.                                      */
1197:   /***************************************************************************/

1199:   alpha = cg->alpha[i];
1200:   if (alpha >= 0.0) sigma = -sigma;

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

1205:   /***************************************************************************/
1206:   /* Accumulate Q * s                                                        */
1207:   /***************************************************************************/

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

1211:   /***************************************************************************/
1212:   /* Set the termination reason.                                             */
1213:   /***************************************************************************/

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

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

1227:   /***************************************************************************/
1228:   /* Determine the total maximum number of iterations.                       */
1229:   /***************************************************************************/

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

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

1238:   KSPSetWorkVecs(ksp, 3);
1239:   if (cg->diag) {
1240:     PetscArrayzero(cg->diag, max_its);
1241:     PetscArrayzero(cg->offd, max_its);
1242:     PetscArrayzero(cg->alpha, max_its);
1243:     PetscArrayzero(cg->beta, max_its);
1244:     PetscArrayzero(cg->norm_r, max_its);
1245:   } else {
1246:     PetscCalloc5(max_its,&cg->diag,max_its,&cg->offd,max_its,&cg->alpha,max_its,&cg->beta,max_its,&cg->norm_r);
1247:     PetscLogObjectMemory((PetscObject)ksp, 5*max_its*sizeof(PetscReal));
1248:   }
1249:   return(0);
1250: }

1252: static PetscErrorCode KSPCGDestroy_GLTR(KSP ksp)
1253: {
1254:   KSPCG_GLTR     *cg = (KSPCG_GLTR*)ksp->data;

1258:   /***************************************************************************/
1259:   /* Free memory allocated for the data.                                     */
1260:   /***************************************************************************/

1262:   PetscFree5(cg->diag,cg->offd,cg->alpha,cg->beta,cg->norm_r);
1263:   if (cg->alloced) {
1264:     PetscFree2(cg->rwork,cg->iwork);
1265:   }

1267:   /***************************************************************************/
1268:   /* Clear composed functions                                                */
1269:   /***************************************************************************/

1271:   PetscObjectComposeFunction((PetscObject)ksp,"KSPCGSetRadius_C",NULL);
1272:   PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetNormD_C",NULL);
1273:   PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetObjFcn_C",NULL);
1274:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetMinEig_C",NULL);
1275:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetLambda_C",NULL);

1277:   /***************************************************************************/
1278:   /* Destroy KSP object.                                                     */
1279:   /***************************************************************************/
1280:   KSPDestroyDefault(ksp);
1281:   return(0);
1282: }

1284: static PetscErrorCode  KSPCGSetRadius_GLTR(KSP ksp, PetscReal radius)
1285: {
1286:   KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;

1289:   cg->radius = radius;
1290:   return(0);
1291: }

1293: static PetscErrorCode  KSPCGGetNormD_GLTR(KSP ksp, PetscReal *norm_d)
1294: {
1295:   KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;

1298:   *norm_d = cg->norm_d;
1299:   return(0);
1300: }

1302: static PetscErrorCode  KSPCGGetObjFcn_GLTR(KSP ksp, PetscReal *o_fcn)
1303: {
1304:   KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;

1307:   *o_fcn = cg->o_fcn;
1308:   return(0);
1309: }

1311: static PetscErrorCode  KSPGLTRGetMinEig_GLTR(KSP ksp, PetscReal *e_min)
1312: {
1313:   KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;

1316:   *e_min = cg->e_min;
1317:   return(0);
1318: }

1320: static PetscErrorCode  KSPGLTRGetLambda_GLTR(KSP ksp, PetscReal *lambda)
1321: {
1322:   KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;

1325:   *lambda = cg->lambda;
1326:   return(0);
1327: }

1329: static PetscErrorCode KSPCGSetFromOptions_GLTR(PetscOptionItems *PetscOptionsObject,KSP ksp)
1330: {
1332:   KSPCG_GLTR       *cg = (KSPCG_GLTR*)ksp->data;

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

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

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

1341:   PetscOptionsReal("-ksp_cg_gltr_init_pert", "Initial perturbation", "", cg->init_pert, &cg->init_pert, NULL);
1342:   PetscOptionsReal("-ksp_cg_gltr_eigen_tol", "Eigenvalue tolerance", "", cg->eigen_tol, &cg->eigen_tol, NULL);
1343:   PetscOptionsReal("-ksp_cg_gltr_newton_tol", "Newton tolerance", "", cg->newton_tol, &cg->newton_tol, NULL);

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

1348:   PetscOptionsTail();
1349:   return(0);
1350: }

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

1357:    Options Database Keys:
1358: .      -ksp_cg_radius <r> - Trust Region Radius

1360:    Notes:
1361:     This is rarely used directly

1363:   Use preconditioned conjugate gradient to compute
1364:   an approximate minimizer of the quadratic function

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

1368:    subject to the trust region constraint

1370:             || s || <= delta,

1372:    where

1374:      delta is the trust region radius,
1375:      g is the gradient vector,
1376:      H is the Hessian approximation,
1377:      M is the positive definite preconditioner matrix.

1379:    KSPConvergedReason may be
1380: $  KSP_CONVERGED_CG_NEG_CURVE if convergence is reached along a negative curvature direction,
1381: $  KSP_CONVERGED_CG_CONSTRAINED if convergence is reached along a constrained step,
1382: $  other KSP converged/diverged reasons

1384:   Notes:
1385:   The preconditioner supplied should be symmetric and positive definite.

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

1391:    Level: developer

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

1396: PETSC_EXTERN PetscErrorCode KSPCreate_GLTR(KSP ksp)
1397: {
1399:   KSPCG_GLTR       *cg;

1402:   PetscNewLog(ksp,&cg);
1403:   cg->radius = 0.0;
1404:   cg->dtype  = GLTR_UNPRECONDITIONED_DIRECTION;

1406:   cg->init_pert  = 1.0e-8;
1407:   cg->eigen_tol  = 1.0e-10;
1408:   cg->newton_tol = 1.0e-6;

1410:   cg->alloced    = 0;
1411:   cg->init_alloc = 1024;

1413:   cg->max_lanczos_its = 20;
1414:   cg->max_newton_its  = 10;

1416:   ksp->data = (void*) cg;
1417:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,3);
1418:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);
1419:   KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,2);

1421:   /***************************************************************************/
1422:   /* Sets the functions that are associated with this data structure         */
1423:   /* (in C++ this is the same as defining virtual functions).                */
1424:   /***************************************************************************/

1426:   ksp->ops->setup          = KSPCGSetUp_GLTR;
1427:   ksp->ops->solve          = KSPCGSolve_GLTR;
1428:   ksp->ops->destroy        = KSPCGDestroy_GLTR;
1429:   ksp->ops->setfromoptions = KSPCGSetFromOptions_GLTR;
1430:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
1431:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
1432:   ksp->ops->view           = 0;

1434:   PetscObjectComposeFunction((PetscObject)ksp,"KSPCGSetRadius_C",KSPCGSetRadius_GLTR);
1435:   PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetNormD_C", KSPCGGetNormD_GLTR);
1436:   PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetObjFcn_C",KSPCGGetObjFcn_GLTR);
1437:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetMinEig_C",KSPGLTRGetMinEig_GLTR);
1438:   PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetLambda_C",KSPGLTRGetLambda_GLTR);
1439:   return(0);
1440: }