Actual source code: qcg.c

      2: #include <../src/ksp/ksp/impls/qcg/qcgimpl.h>

      4: static PetscErrorCode KSPQCGQuadraticRoots(Vec,Vec,PetscReal,PetscReal*,PetscReal*);

      6: /*@
      7:     KSPQCGSetTrustRegionRadius - Sets the radius of the trust region.

      9:     Logically Collective on ksp

     11:     Input Parameters:
     12: +   ksp   - the iterative context
     13: -   delta - the trust region radius (Infinity is the default)

     15:     Options Database Key:
     16: .   -ksp_qcg_trustregionradius <delta>

     18:     Level: advanced

     20: @*/
     21: PetscErrorCode  KSPQCGSetTrustRegionRadius(KSP ksp,PetscReal delta)
     22: {

     27:   if (delta < 0.0) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Tolerance must be non-negative");
     28:   PetscTryMethod(ksp,"KSPQCGSetTrustRegionRadius_C",(KSP,PetscReal),(ksp,delta));
     29:   return(0);
     30: }

     32: /*@
     33:     KSPQCGGetTrialStepNorm - Gets the norm of a trial step vector.  The WCG step may be
     34:     constrained, so this is not necessarily the length of the ultimate step taken in QCG.

     36:     Not Collective

     38:     Input Parameter:
     39: .   ksp - the iterative context

     41:     Output Parameter:
     42: .   tsnorm - the norm

     44:     Level: advanced
     45: @*/
     46: PetscErrorCode  KSPQCGGetTrialStepNorm(KSP ksp,PetscReal *tsnorm)
     47: {

     52:   PetscUseMethod(ksp,"KSPQCGGetTrialStepNorm_C",(KSP,PetscReal*),(ksp,tsnorm));
     53:   return(0);
     54: }

     56: /*@
     57:     KSPQCGGetQuadratic - Gets the value of the quadratic function, evaluated at the new iterate:

     59:        q(s) = g^T * s + 0.5 * s^T * H * s

     61:     which satisfies the Euclidian Norm trust region constraint

     63:        || D * s || <= delta,

     65:     where

     67:      delta is the trust region radius,
     68:      g is the gradient vector, and
     69:      H is Hessian matrix,
     70:      D is a scaling matrix.

     72:     Collective on ksp

     74:     Input Parameter:
     75: .   ksp - the iterative context

     77:     Output Parameter:
     78: .   quadratic - the quadratic function evaluated at the new iterate

     80:     Level: advanced
     81: @*/
     82: PetscErrorCode  KSPQCGGetQuadratic(KSP ksp,PetscReal *quadratic)
     83: {

     88:   PetscUseMethod(ksp,"KSPQCGGetQuadratic_C",(KSP,PetscReal*),(ksp,quadratic));
     89:   return(0);
     90: }

     93: PetscErrorCode KSPSolve_QCG(KSP ksp)
     94: {
     95: /*
     96:    Correpondence with documentation above:
     97:       B = g = gradient,
     98:       X = s = step
     99:    Note:  This is not coded correctly for complex arithmetic!
    100:  */

    102:   KSP_QCG        *pcgP = (KSP_QCG*)ksp->data;
    103:   Mat            Amat,Pmat;
    104:   Vec            W,WA,WA2,R,P,ASP,BS,X,B;
    105:   PetscScalar    scal,beta,rntrn,step;
    106:   PetscReal      q1,q2,xnorm,step1,step2,rnrm = 0.0,btx,xtax;
    107:   PetscReal      ptasp,rtr,wtasp,bstp;
    108:   PetscReal      dzero = 0.0,bsnrm = 0.0;
    110:   PetscInt       i,maxit;
    111:   PC             pc = ksp->pc;
    112:   PCSide         side;
    113:   PetscBool      diagonalscale;

    116:   PCGetDiagonalScale(ksp->pc,&diagonalscale);
    117:   if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);
    118:   if (ksp->transpose_solve) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Currently does not support transpose solve");

    120:   ksp->its = 0;
    121:   maxit    = ksp->max_it;
    122:   WA       = ksp->work[0];
    123:   R        = ksp->work[1];
    124:   P        = ksp->work[2];
    125:   ASP      = ksp->work[3];
    126:   BS       = ksp->work[4];
    127:   W        = ksp->work[5];
    128:   WA2      = ksp->work[6];
    129:   X        = ksp->vec_sol;
    130:   B        = ksp->vec_rhs;

    132:   if (pcgP->delta <= dzero) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Input error: delta <= 0");
    133:   KSPGetPCSide(ksp,&side);
    134:   if (side != PC_SYMMETRIC) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Requires symmetric preconditioner!");

    136:   /* Initialize variables */
    137:   VecSet(W,0.0);  /* W = 0 */
    138:   VecSet(X,0.0);  /* X = 0 */
    139:   PCGetOperators(pc,&Amat,&Pmat);

    141:   /* Compute:  BS = D^{-1} B */
    142:   PCApplySymmetricLeft(pc,B,BS);

    144:   if (ksp->normtype != KSP_NORM_NONE) {
    145:     VecNorm(BS,NORM_2,&bsnrm);
    146:     KSPCheckNorm(ksp,bsnrm);
    147:   }
    148:   PetscObjectSAWsTakeAccess((PetscObject)ksp);
    149:   ksp->its   = 0;
    150:   ksp->rnorm = bsnrm;
    151:   PetscObjectSAWsGrantAccess((PetscObject)ksp);
    152:   KSPLogResidualHistory(ksp,bsnrm);
    153:   KSPMonitor(ksp,0,bsnrm);
    154:   (*ksp->converged)(ksp,0,bsnrm,&ksp->reason,ksp->cnvP);
    155:   if (ksp->reason) return(0);

    157:   /* Compute the initial scaled direction and scaled residual */
    158:   VecCopy(BS,R);
    159:   VecScale(R,-1.0);
    160:   VecCopy(R,P);
    161:   VecDotRealPart(R,R,&rtr);

    163:   for (i=0; i<=maxit; i++) {
    164:     PetscObjectSAWsTakeAccess((PetscObject)ksp);
    165:     ksp->its++;
    166:     PetscObjectSAWsGrantAccess((PetscObject)ksp);

    168:     /* Compute:  asp = D^{-T}*A*D^{-1}*p  */
    169:     PCApplySymmetricRight(pc,P,WA);
    170:     KSP_MatMult(ksp,Amat,WA,WA2);
    171:     PCApplySymmetricLeft(pc,WA2,ASP);

    173:     /* Check for negative curvature */
    174:     VecDotRealPart(P,ASP,&ptasp);
    175:     if (ptasp <= dzero) {

    177:       /* Scaled negative curvature direction:  Compute a step so that
    178:         ||w + step*p|| = delta and QS(w + step*p) is least */

    180:       if (!i) {
    181:         VecCopy(P,X);
    182:         VecNorm(X,NORM_2,&xnorm);
    183:         KSPCheckNorm(ksp,xnorm);
    184:         scal = pcgP->delta / xnorm;
    185:         VecScale(X,scal);
    186:       } else {
    187:         /* Compute roots of quadratic */
    188:         KSPQCGQuadraticRoots(W,P,pcgP->delta,&step1,&step2);
    189:         VecDotRealPart(W,ASP,&wtasp);
    190:         VecDotRealPart(BS,P,&bstp);
    191:         VecCopy(W,X);
    192:         q1   = step1*(bstp + wtasp + .5*step1*ptasp);
    193:         q2   = step2*(bstp + wtasp + .5*step2*ptasp);
    194:         if (q1 <= q2) {
    195:           VecAXPY(X,step1,P);
    196:         } else {
    197:           VecAXPY(X,step2,P);
    198:         }
    199:       }
    200:       pcgP->ltsnrm = pcgP->delta;                       /* convergence in direction of */
    201:       ksp->reason  = KSP_CONVERGED_CG_NEG_CURVE;  /* negative curvature */
    202:       if (!i) {
    203:         PetscInfo1(ksp,"negative curvature: delta=%g\n",(double)pcgP->delta);
    204:       } else {
    205:         PetscInfo3(ksp,"negative curvature: step1=%g, step2=%g, delta=%g\n",(double)step1,(double)step2,(double)pcgP->delta);
    206:       }

    208:     } else {
    209:       /* Compute step along p */
    210:       step = rtr/ptasp;
    211:       VecCopy(W,X);        /*  x = w  */
    212:       VecAXPY(X,step,P);   /*  x <- step*p + x  */
    213:       VecNorm(X,NORM_2,&pcgP->ltsnrm);
    214:       KSPCheckNorm(ksp,pcgP->ltsnrm);

    216:       if (pcgP->ltsnrm > pcgP->delta) {
    217:         /* Since the trial iterate is outside the trust region,
    218:             evaluate a constrained step along p so that
    219:                     ||w + step*p|| = delta
    220:           The positive step is always better in this case. */
    221:         if (!i) {
    222:           scal = pcgP->delta / pcgP->ltsnrm;
    223:           VecScale(X,scal);
    224:         } else {
    225:           /* Compute roots of quadratic */
    226:           KSPQCGQuadraticRoots(W,P,pcgP->delta,&step1,&step2);
    227:           VecCopy(W,X);
    228:           VecAXPY(X,step1,P);  /*  x <- step1*p + x  */
    229:         }
    230:         pcgP->ltsnrm = pcgP->delta;
    231:         ksp->reason  = KSP_CONVERGED_CG_CONSTRAINED; /* convergence along constrained step */
    232:         if (!i) {
    233:           PetscInfo1(ksp,"constrained step: delta=%g\n",(double)pcgP->delta);
    234:         } else {
    235:           PetscInfo3(ksp,"constrained step: step1=%g, step2=%g, delta=%g\n",(double)step1,(double)step2,(double)pcgP->delta);
    236:         }

    238:       } else {
    239:         /* Evaluate the current step */
    240:         VecCopy(X,W);  /* update interior iterate */
    241:         VecAXPY(R,-step,ASP); /* r <- -step*asp + r */
    242:         if (ksp->normtype != KSP_NORM_NONE) {
    243:           VecNorm(R,NORM_2,&rnrm);
    244:           KSPCheckNorm(ksp,rnrm);
    245:         }
    246:         PetscObjectSAWsTakeAccess((PetscObject)ksp);
    247:         ksp->rnorm = rnrm;
    248:         PetscObjectSAWsGrantAccess((PetscObject)ksp);
    249:         KSPLogResidualHistory(ksp,rnrm);
    250:         KSPMonitor(ksp,i+1,rnrm);
    251:         (*ksp->converged)(ksp,i+1,rnrm,&ksp->reason,ksp->cnvP);
    252:         if (ksp->reason) {                 /* convergence for */
    253:           PetscInfo3(ksp,"truncated step: step=%g, rnrm=%g, delta=%g\n",(double)PetscRealPart(step),(double)rnrm,(double)pcgP->delta);
    254:         }
    255:       }
    256:     }
    257:     if (ksp->reason) break;  /* Convergence has been attained */
    258:     else {                   /* Compute a new AS-orthogonal direction */
    259:       VecDot(R,R,&rntrn);
    260:       beta = rntrn/rtr;
    261:       VecAYPX(P,beta,R);  /*  p <- r + beta*p  */
    262:       rtr  = PetscRealPart(rntrn);
    263:     }
    264:   }
    265:   if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;

    267:   /* Unscale x */
    268:   VecCopy(X,WA2);
    269:   PCApplySymmetricRight(pc,WA2,X);

    271:   KSP_MatMult(ksp,Amat,X,WA);
    272:   VecDotRealPart(B,X,&btx);
    273:   VecDotRealPart(X,WA,&xtax);

    275:   pcgP->quadratic = btx + .5*xtax;
    276:   return(0);
    277: }

    279: PetscErrorCode KSPSetUp_QCG(KSP ksp)
    280: {

    284:   /* Get work vectors from user code */
    285:   KSPSetWorkVecs(ksp,7);
    286:   return(0);
    287: }

    289: PetscErrorCode KSPDestroy_QCG(KSP ksp)
    290: {

    294:   PetscObjectComposeFunction((PetscObject)ksp,"KSPQCGGetQuadratic_C",NULL);
    295:   PetscObjectComposeFunction((PetscObject)ksp,"KSPQCGGetTrialStepNorm_C",NULL);
    296:   PetscObjectComposeFunction((PetscObject)ksp,"KSPQCGSetTrustRegionRadius_C",NULL);
    297:   KSPDestroyDefault(ksp);
    298:   return(0);
    299: }

    301: static PetscErrorCode  KSPQCGSetTrustRegionRadius_QCG(KSP ksp,PetscReal delta)
    302: {
    303:   KSP_QCG *cgP = (KSP_QCG*)ksp->data;

    306:   cgP->delta = delta;
    307:   return(0);
    308: }

    310: static PetscErrorCode  KSPQCGGetTrialStepNorm_QCG(KSP ksp,PetscReal *ltsnrm)
    311: {
    312:   KSP_QCG *cgP = (KSP_QCG*)ksp->data;

    315:   *ltsnrm = cgP->ltsnrm;
    316:   return(0);
    317: }

    319: static PetscErrorCode  KSPQCGGetQuadratic_QCG(KSP ksp,PetscReal *quadratic)
    320: {
    321:   KSP_QCG *cgP = (KSP_QCG*)ksp->data;

    324:   *quadratic = cgP->quadratic;
    325:   return(0);
    326: }

    328: PetscErrorCode KSPSetFromOptions_QCG(PetscOptionItems *PetscOptionsObject,KSP ksp)
    329: {
    331:   PetscReal      delta;
    332:   KSP_QCG        *cgP = (KSP_QCG*)ksp->data;
    333:   PetscBool      flg;

    336:   PetscOptionsHead(PetscOptionsObject,"KSP QCG Options");
    337:   PetscOptionsReal("-ksp_qcg_trustregionradius","Trust Region Radius","KSPQCGSetTrustRegionRadius",cgP->delta,&delta,&flg);
    338:   if (flg) { KSPQCGSetTrustRegionRadius(ksp,delta); }
    339:   PetscOptionsTail();
    340:   return(0);
    341: }

    343: /*MC
    344:      KSPQCG -   Code to run conjugate gradient method subject to a constraint
    345:          on the solution norm. This is used in Trust Region methods for nonlinear equations, SNESNEWTONTR

    347:    Options Database Keys:
    348: .      -ksp_qcg_trustregionradius <r> - Trust Region Radius

    350:    Notes:
    351:     This is rarely used directly

    353:    Level: developer

    355:   Notes:
    356:     Use preconditioned conjugate gradient to compute
    357:       an approximate minimizer of the quadratic function

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

    361:    subject to the Euclidean norm trust region constraint

    363:             || D * s || <= delta,

    365:    where

    367:      delta is the trust region radius,
    368:      g is the gradient vector, and
    369:      H is Hessian matrix,
    370:      D is a scaling matrix.

    372:    KSPConvergedReason may be
    373: $  KSP_CONVERGED_CG_NEG_CURVE if convergence is reached along a negative curvature direction,
    374: $  KSP_CONVERGED_CG_CONSTRAINED if convergence is reached along a constrained step,
    375: $  other KSP converged/diverged reasons

    378:   Notes:
    379:   Currently we allow symmetric preconditioning with the following scaling matrices:
    380:       PCNONE:   D = Identity matrix
    381:       PCJACOBI: D = diag [d_1, d_2, ...., d_n], where d_i = sqrt(H[i,i])
    382:       PCICC:    D = L^T, implemented with forward and backward solves.
    383:                 Here L is an incomplete Cholesky factor of H.

    385:   References:
    386: .  1. - Trond Steihaug, The Conjugate Gradient Method and Trust Regions in Large Scale Optimization,
    387:    SIAM Journal on Numerical Analysis, Vol. 20, No. 3 (Jun., 1983).

    389: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPQCGSetTrustRegionRadius()
    390:            KSPQCGGetTrialStepNorm(), KSPQCGGetQuadratic()
    391: M*/

    393: PETSC_EXTERN PetscErrorCode KSPCreate_QCG(KSP ksp)
    394: {
    396:   KSP_QCG        *cgP;

    399:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_SYMMETRIC,3);
    400:   KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_SYMMETRIC,1);
    401:   PetscNewLog(ksp,&cgP);

    403:   ksp->data                = (void*)cgP;
    404:   ksp->ops->setup          = KSPSetUp_QCG;
    405:   ksp->ops->setfromoptions = KSPSetFromOptions_QCG;
    406:   ksp->ops->solve          = KSPSolve_QCG;
    407:   ksp->ops->destroy        = KSPDestroy_QCG;
    408:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
    409:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
    410:   ksp->ops->view           = NULL;

    412:   PetscObjectComposeFunction((PetscObject)ksp,"KSPQCGGetQuadratic_C",KSPQCGGetQuadratic_QCG);
    413:   PetscObjectComposeFunction((PetscObject)ksp,"KSPQCGGetTrialStepNorm_C",KSPQCGGetTrialStepNorm_QCG);
    414:   PetscObjectComposeFunction((PetscObject)ksp,"KSPQCGSetTrustRegionRadius_C",KSPQCGSetTrustRegionRadius_QCG);
    415:   cgP->delta = PETSC_MAX_REAL; /* default trust region radius is infinite */
    416:   return(0);
    417: }

    419: /* ---------------------------------------------------------- */
    420: /*
    421:   KSPQCGQuadraticRoots - Computes the roots of the quadratic,
    422:          ||s + step*p|| - delta = 0
    423:    such that step1 >= 0 >= step2.
    424:    where
    425:       delta:
    426:         On entry delta must contain scalar delta.
    427:         On exit delta is unchanged.
    428:       step1:
    429:         On entry step1 need not be specified.
    430:         On exit step1 contains the non-negative root.
    431:       step2:
    432:         On entry step2 need not be specified.
    433:         On exit step2 contains the non-positive root.
    434:    C code is translated from the Fortran version of the MINPACK-2 Project,
    435:    Argonne National Laboratory, Brett M. Averick and Richard G. Carter.
    436: */
    437: static PetscErrorCode KSPQCGQuadraticRoots(Vec s,Vec p,PetscReal delta,PetscReal *step1,PetscReal *step2)
    438: {
    439:   PetscReal      dsq,ptp,pts,rad,sts;

    443:   VecDotRealPart(p,s,&pts);
    444:   VecDotRealPart(p,p,&ptp);
    445:   VecDotRealPart(s,s,&sts);
    446:   dsq  = delta*delta;
    447:   rad  = PetscSqrtReal((pts*pts) - ptp*(sts - dsq));
    448:   if (pts > 0.0) {
    449:     *step2 = -(pts + rad)/ptp;
    450:     *step1 = (sts - dsq)/(ptp * *step2);
    451:   } else {
    452:     *step1 = -(pts - rad)/ptp;
    453:     *step2 = (sts - dsq)/(ptp * *step1);
    454:   }
    455:   return(0);
    456: }