Actual source code: qcg.c

      1: #define PETSCKSP_DLL

     3:  #include src/ksp/ksp/kspimpl.h
     4:  #include src/ksp/ksp/impls/qcg/qcg.h

      6: static PetscErrorCode QuadraticRoots_Private(Vec,Vec,PetscReal*,PetscReal*,PetscReal*);

     10: /*@
     11:     KSPQCGSetTrustRegionRadius - Sets the radius of the trust region.

     13:     Collective on KSP

     15:     Input Parameters:
     16: +   ksp   - the iterative context
     17: -   delta - the trust region radius (Infinity is the default)

     19:     Options Database Key:
     20: .   -ksp_qcg_trustregionradius <delta>

     22:     Level: advanced

     24: .keywords: KSP, QCG, set, trust region radius
     25: @*/
     26: PetscErrorCode PETSCKSP_DLLEXPORT KSPQCGSetTrustRegionRadius(KSP ksp,PetscReal delta)
     27: {
     28:   PetscErrorCode ierr,(*f)(KSP,PetscReal);

     32:   if (delta < 0.0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Tolerance must be non-negative");
     33:   PetscObjectQueryFunction((PetscObject)ksp,"KSPQCGSetTrustRegionRadius_C",(void (**)(void))&f);
     34:   if (f) {
     35:     (*f)(ksp,delta);
     36:   }

     38:   return(0);
     39: }

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

     47:     Collective on KSP

     49:     Input Parameter:
     50: .   ksp - the iterative context

     52:     Output Parameter:
     53: .   tsnorm - the norm

     55:     Level: advanced
     56: @*/
     57: PetscErrorCode PETSCKSP_DLLEXPORT KSPQCGGetTrialStepNorm(KSP ksp,PetscReal *tsnorm)
     58: {
     59:   PetscErrorCode ierr,(*f)(KSP,PetscReal*);

     63:   PetscObjectQueryFunction((PetscObject)ksp,"KSPQCGGetTrialStepNorm_C",(void (**)(void))&f);
     64:   if (f) {
     65:     (*f)(ksp,tsnorm);
     66:   }
     67:   return(0);
     68: }

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

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

     77:     which satisfies the Euclidian Norm trust region constraint

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

     81:     where

     83:      delta is the trust region radius,
     84:      g is the gradient vector, and
     85:      H is Hessian matrix,
     86:      D is a scaling matrix.

     88:     Collective on KSP

     90:     Input Parameter:
     91: .   ksp - the iterative context

     93:     Output Parameter:
     94: .   quadratic - the quadratic function evaluated at the new iterate

     96:     Level: advanced
     97: @*/
     98: PetscErrorCode PETSCKSP_DLLEXPORT KSPQCGGetQuadratic(KSP ksp,PetscReal *quadratic)
     99: {
    100:   PetscErrorCode ierr,(*f)(KSP,PetscReal*);

    104:   PetscObjectQueryFunction((PetscObject)ksp,"KSPQCGGetQuadratic_C",(void (**)(void))&f);
    105:   if (f) {
    106:     (*f)(ksp,quadratic);
    107:   }
    108:   return(0);
    109: }

    113: /*
    114:   KSPSolve_QCG - Use preconditioned conjugate gradient to compute
    115:   an approximate minimizer of the quadratic function

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

    119:    subject to the Euclidean norm trust region constraint

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

    123:    where

    125:      delta is the trust region radius,
    126:      g is the gradient vector, and
    127:      H is Hessian matrix,
    128:      D is a scaling matrix.

    130:    KSPConvergedReason may be
    131: $  KSP_CONVERGED_QCG_NEG_CURVE if convergence is reached along a negative curvature direction,
    132: $  KSP_CONVERGED_QCG_CONSTRAINED if convergence is reached along a constrained step,
    133: $  other KSP converged/diverged reasons

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

    143:  We should perhaps rewrite using PCApplyBAorAB().
    144:  */
    145: PetscErrorCode KSPSolve_QCG(KSP ksp)
    146: {
    147: /*
    148:    Correpondence with documentation above:
    149:       B = g = gradient,
    150:       X = s = step
    151:    Note:  This is not coded correctly for complex arithmetic!
    152:  */

    154:   KSP_QCG        *pcgP = (KSP_QCG*)ksp->data;
    155:   MatStructure   pflag;
    156:   Mat            Amat,Pmat;
    157:   Vec            W,WA,WA2,R,P,ASP,BS,X,B;
    158:   PetscScalar    zero = 0.0,negone = -1.0,scal,nstep,btx,xtax,beta,rntrn,step;
    159:   PetscReal      ptasp,q1,q2,wtasp,bstp,rtr,xnorm,step1,step2,rnrm,p5 = 0.5;
    160:   PetscReal      dzero = 0.0,bsnrm;
    162:   PetscInt       i,maxit;
    163:   PC             pc = ksp->pc;
    164:   PCSide         side;
    165: #if defined(PETSC_USE_COMPLEX)
    166:   PetscScalar    cstep1,cstep2,cbstp,crtr,cwtasp,cptasp;
    167: #endif
    168:   PetscTruth     diagonalscale;

    171:   PCDiagonalScale(ksp->pc,&diagonalscale);
    172:   if (diagonalscale) SETERRQ1(PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",ksp->type_name);
    173:   if (ksp->transpose_solve) {
    174:     SETERRQ(PETSC_ERR_SUP,"Currently does not support transpose solve");
    175:   }

    177:   ksp->its = 0;
    178:   maxit    = ksp->max_it;
    179:   WA       = ksp->work[0];
    180:   R        = ksp->work[1];
    181:   P        = ksp->work[2];
    182:   ASP      = ksp->work[3];
    183:   BS       = ksp->work[4];
    184:   W        = ksp->work[5];
    185:   WA2      = ksp->work[6];
    186:   X        = ksp->vec_sol;
    187:   B        = ksp->vec_rhs;

    189:   if (pcgP->delta <= dzero) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Input error: delta <= 0");
    190:   KSPGetPreconditionerSide(ksp,&side);
    191:   if (side != PC_SYMMETRIC) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Requires symmetric preconditioner!");

    193:   /* Initialize variables */
    194:   VecSet(W,zero);        /* W = 0 */
    195:   VecSet(X,zero);        /* X = 0 */
    196:   PCGetOperators(pc,&Amat,&Pmat,&pflag);

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

    201:   VecNorm(BS,NORM_2,&bsnrm);
    202:   PetscObjectTakeAccess(ksp);
    203:   ksp->its    = 0;
    204:   ksp->rnorm  = bsnrm;
    205:   PetscObjectGrantAccess(ksp);
    206:   KSPLogResidualHistory(ksp,bsnrm);
    207:   KSPMonitor(ksp,0,bsnrm);
    208:   (*ksp->converged)(ksp,0,bsnrm,&ksp->reason,ksp->cnvP);
    209:   if (ksp->reason) return(0);

    211:   /* Compute the initial scaled direction and scaled residual */
    212:   VecCopy(BS,R);
    213:   VecScale(R,negone);
    214:   VecCopy(R,P);
    215: #if defined(PETSC_USE_COMPLEX)
    216:   VecDot(R,R,&crtr); rtr = PetscRealPart(crtr);
    217: #else
    218:   VecDot(R,R,&rtr);
    219: #endif

    221:   for (i=0; i<=maxit; i++) {
    222:     PetscObjectTakeAccess(ksp);
    223:     ksp->its++;
    224:     PetscObjectGrantAccess(ksp);

    226:     /* Compute:  asp = D^{-T}*A*D^{-1}*p  */
    227:     PCApplySymmetricRight(pc,P,WA);
    228:     MatMult(Amat,WA,WA2);
    229:     PCApplySymmetricLeft(pc,WA2,ASP);

    231:     /* Check for negative curvature */
    232: #if defined(PETSC_USE_COMPLEX)
    233:     VecDot(P,ASP,&cptasp);
    234:     ptasp = PetscRealPart(cptasp);
    235: #else
    236:     VecDot(P,ASP,&ptasp);        /* ptasp = p^T asp */
    237: #endif
    238:     if (ptasp <= dzero) {

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

    243:        if (!i) {
    244:          VecCopy(P,X);
    245:          VecNorm(X,NORM_2,&xnorm);
    246:          scal = pcgP->delta / xnorm;
    247:          VecScale(X,scal);
    248:        } else {
    249:          /* Compute roots of quadratic */
    250:          QuadraticRoots_Private(W,P,&pcgP->delta,&step1,&step2);
    251: #if defined(PETSC_USE_COMPLEX)
    252:          VecDot(W,ASP,&cwtasp); wtasp = PetscRealPart(cwtasp);
    253:          VecDot(BS,P,&cbstp);   bstp  = PetscRealPart(cbstp);
    254: #else
    255:          VecDot(W,ASP,&wtasp);
    256:          VecDot(BS,P,&bstp);
    257: #endif
    258:          VecCopy(W,X);
    259:          q1 = step1*(bstp + wtasp + p5*step1*ptasp);
    260:          q2 = step2*(bstp + wtasp + p5*step2*ptasp);
    261: #if defined(PETSC_USE_COMPLEX)
    262:          if (q1 <= q2) {
    263:            cstep1 = step1; VecAXPY(X,cstep1,P);
    264:          } else {
    265:            cstep2 = step2; VecAXPY(X,cstep2,P);
    266:          }
    267: #else
    268:          if (q1 <= q2) {VecAXPY(X,step1,P);}
    269:          else          {VecAXPY(X,step2,P);}
    270: #endif
    271:        }
    272:        pcgP->ltsnrm = pcgP->delta;                       /* convergence in direction of */
    273:        ksp->reason  = KSP_CONVERGED_QCG_NEG_CURVE;  /* negative curvature */
    274:        if (!i) {
    275:          PetscLogInfo((ksp,"KSPSolve_QCG: negative curvature: delta=%g\n",pcgP->delta));
    276:        } else {
    277:          PetscLogInfo((ksp,"KSPSolve_QCG: negative curvature: step1=%g, step2=%g, delta=%g\n",step1,step2,pcgP->delta));
    278:        }
    279:
    280:     } else {
    281:
    282:        /* Compute step along p */

    284:        step = rtr/ptasp;
    285:        VecCopy(W,X);           /*  x = w  */
    286:        VecAXPY(X,step,P);   /*  x <- step*p + x  */
    287:        VecNorm(X,NORM_2,&pcgP->ltsnrm);

    289:        if (pcgP->ltsnrm > pcgP->delta) {

    291:          /* Since the trial iterate is outside the trust region,
    292:              evaluate a constrained step along p so that
    293:                       ||w + step*p|| = delta
    294:             The positive step is always better in this case. */

    296:          if (!i) {
    297:            scal = pcgP->delta / pcgP->ltsnrm;
    298:            VecScale(X,scal);
    299:          } else {
    300:            /* Compute roots of quadratic */
    301:            QuadraticRoots_Private(W,P,&pcgP->delta,&step1,&step2);
    302:            VecCopy(W,X);
    303: #if defined(PETSC_USE_COMPLEX)
    304:            cstep1 = step1; VecAXPY(X,cstep1,P);
    305: #else
    306:            VecAXPY(X,step1,P);  /*  x <- step1*p + x  */
    307: #endif
    308:          }
    309:          pcgP->ltsnrm = pcgP->delta;
    310:          ksp->reason  = KSP_CONVERGED_QCG_CONSTRAINED;        /* convergence along constrained step */
    311:          if (!i) {
    312:            PetscLogInfo((ksp,"KSPSolve_QCG: constrained step: delta=%g\n",pcgP->delta));
    313:          } else {
    314:            PetscLogInfo((ksp,"KSPSolve_QCG: constrained step: step1=%g, step2=%g, delta=%g\n",step1,step2,pcgP->delta));
    315:          }

    317:        } else {

    319:          /* Evaluate the current step */

    321:          VecCopy(X,W);        /* update interior iterate */
    322:          nstep = -step;
    323:          VecAXPY(R,nstep,ASP); /* r <- -step*asp + r */
    324:          VecNorm(R,NORM_2,&rnrm);

    326:          PetscObjectTakeAccess(ksp);
    327:          ksp->rnorm                                    = rnrm;
    328:          PetscObjectGrantAccess(ksp);
    329:          KSPLogResidualHistory(ksp,rnrm);
    330:          KSPMonitor(ksp,i+1,rnrm);
    331:          (*ksp->converged)(ksp,i+1,rnrm,&ksp->reason,ksp->cnvP);
    332:          if (ksp->reason) {                 /* convergence for */
    333: #if defined(PETSC_USE_COMPLEX)
    334:            PetscLogInfo((ksp,"KSPSolve_QCG: truncated step: step=%g, rnrm=%g, delta=%g\n",PetscRealPart(step),rnrm,pcgP->delta));
    335: #else
    336:            PetscLogInfo((ksp,"KSPSolve_QCG: truncated step: step=%g, rnrm=%g, delta=%g\n",step,rnrm,pcgP->delta));
    337: #endif
    338:          }
    339:       }
    340:     }
    341:     if (ksp->reason) break;        /* Convergence has been attained */
    342:     else {                /* Compute a new AS-orthogonal direction */
    343:       VecDot(R,R,&rntrn);
    344:       beta = rntrn/rtr;
    345:       VecAYPX(P,beta,R);        /*  p <- r + beta*p  */
    346: #if defined(PETSC_USE_COMPLEX)
    347:       rtr = PetscRealPart(rntrn);
    348: #else
    349:       rtr = rntrn;
    350: #endif
    351:     }
    352:   }
    353:   if (!ksp->reason) {
    354:     ksp->reason = KSP_DIVERGED_ITS;
    355:   }

    357:   /* Unscale x */
    358:   VecCopy(X,WA2);
    359:   PCApplySymmetricRight(pc,WA2,X);

    361:   MatMult(Amat,X,WA);
    362:   VecDot(B,X,&btx);
    363:   VecDot(X,WA,&xtax);
    364: #if defined(PETSC_USE_COMPLEX)
    365:   pcgP->quadratic = PetscRealPart(btx) + p5* PetscRealPart(xtax);
    366: #else
    367:   pcgP->quadratic = btx + p5*xtax;              /* Compute q(x) */
    368: #endif
    369:   return(0);
    370: }

    374: PetscErrorCode KSPSetUp_QCG(KSP ksp)
    375: {

    379:   /* Check user parameters and functions */
    380:   if (ksp->pc_side == PC_RIGHT) {
    381:     SETERRQ(PETSC_ERR_SUP,"no right preconditioning for QCG");
    382:   } else if (ksp->pc_side == PC_LEFT) {
    383:     SETERRQ(PETSC_ERR_SUP,"no left preconditioning for QCG");
    384:   }

    386:   /* Get work vectors from user code */
    387:   KSPDefaultGetWork(ksp,7);
    388:   return(0);
    389: }

    393: PetscErrorCode KSPDestroy_QCG(KSP ksp)
    394: {
    395:   KSP_QCG        *cgP = (KSP_QCG*)ksp->data;

    399:   KSPDefaultFreeWork(ksp);
    400:
    401:   /* Free the context variable */
    402:   PetscFree(cgP);
    403:   return(0);
    404: }

    409: PetscErrorCode PETSCKSP_DLLEXPORT KSPQCGSetTrustRegionRadius_QCG(KSP ksp,PetscReal delta)
    410: {
    411:   KSP_QCG *cgP = (KSP_QCG*)ksp->data;

    414:   cgP->delta = delta;
    415:   return(0);
    416: }

    422: PetscErrorCode PETSCKSP_DLLEXPORT KSPQCGGetTrialStepNorm_QCG(KSP ksp,PetscReal *ltsnrm)
    423: {
    424:   KSP_QCG *cgP = (KSP_QCG*)ksp->data;

    427:   *ltsnrm = cgP->ltsnrm;
    428:   return(0);
    429: }

    435: PetscErrorCode PETSCKSP_DLLEXPORT KSPQCGGetQuadratic_QCG(KSP ksp,PetscReal *quadratic)
    436: {
    437:   KSP_QCG *cgP = (KSP_QCG*)ksp->data;

    440:   *quadratic = cgP->quadratic;
    441:   return(0);
    442: }

    447: PetscErrorCode KSPSetFromOptions_QCG(KSP ksp)
    448: {
    450:   PetscReal      delta;
    451:   KSP_QCG        *cgP = (KSP_QCG*)ksp->data;
    452:   PetscTruth     flg;

    455:   PetscOptionsHead("KSP QCG Options");
    456:   PetscOptionsReal("-ksp_qcg_trustregionradius","Trust Region Radius","KSPQCGSetTrustRegionRadius",cgP->delta,&delta,&flg);
    457:   if (flg) { KSPQCGSetTrustRegionRadius(ksp,delta); }
    458:   PetscOptionsTail();
    459:   return(0);
    460: }

    462: /*MC
    463:      KSPQCG -   Code to run conjugate gradient method subject to a constraint
    464:          on the solution norm. This is used in Trust Region methods for nonlinear equations, SNESTR

    466:    Options Database Keys:
    467: .      -ksp_qcg_trustregionradius <r> - Trust Region Radius

    469:    Notes: This is rarely used directly

    471:    Level: developer

    473: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPQCGSetTrustRegionRadius()
    474:            KSPQCGGetTrialStepNorm(), KSPQCGGetQuadratic()
    475: M*/

    480: PetscErrorCode PETSCKSP_DLLEXPORT KSPCreate_QCG(KSP ksp)
    481: {
    483:   KSP_QCG        *cgP;

    486:   PetscMalloc(sizeof(KSP_QCG),&cgP);
    487:   PetscMemzero(cgP,sizeof(KSP_QCG));
    488:   PetscLogObjectMemory(ksp,sizeof(KSP_QCG));
    489:   ksp->data                      = (void*)cgP;
    490:   ksp->pc_side                   = PC_SYMMETRIC;
    491:   ksp->ops->setup                = KSPSetUp_QCG;
    492:   ksp->ops->setfromoptions       = KSPSetFromOptions_QCG;
    493:   ksp->ops->solve                = KSPSolve_QCG;
    494:   ksp->ops->destroy              = KSPDestroy_QCG;
    495:   ksp->ops->buildsolution        = KSPDefaultBuildSolution;
    496:   ksp->ops->buildresidual        = KSPDefaultBuildResidual;
    497:   ksp->ops->setfromoptions       = 0;
    498:   ksp->ops->view                 = 0;

    500:   PetscObjectComposeFunctionDynamic((PetscObject)ksp,"KSPQCGGetQuadratic_C",
    501:                                     "KSPQCGGetQuadratic_QCG",
    502:                                      KSPQCGGetQuadratic_QCG);
    503:   PetscObjectComposeFunctionDynamic((PetscObject)ksp,"KSPQCGGetTrialStepNorm_C",
    504:                                     "KSPQCGGetTrialStepNorm_QCG",
    505:                                      KSPQCGGetTrialStepNorm_QCG);
    506:   PetscObjectComposeFunctionDynamic((PetscObject)ksp,"KSPQCGSetTrustRegionRadius_C",
    507:                                     "KSPQCGSetTrustRegionRadius_QCG",
    508:                                      KSPQCGSetTrustRegionRadius_QCG);
    509:   cgP->delta = PETSC_MAX; /* default trust region radius is infinite */
    510:   return(0);
    511: }

    514: /* ---------------------------------------------------------- */
    517: /*
    518:   QuadraticRoots_Private - Computes the roots of the quadratic,
    519:          ||s + step*p|| - delta = 0
    520:    such that step1 >= 0 >= step2.
    521:    where
    522:       delta:
    523:         On entry delta must contain scalar delta.
    524:         On exit delta is unchanged.
    525:       step1:
    526:         On entry step1 need not be specified.
    527:         On exit step1 contains the non-negative root.
    528:       step2:
    529:         On entry step2 need not be specified.
    530:         On exit step2 contains the non-positive root.
    531:    C code is translated from the Fortran version of the MINPACK-2 Project,
    532:    Argonne National Laboratory, Brett M. Averick and Richard G. Carter.
    533: */
    534: static PetscErrorCode QuadraticRoots_Private(Vec s,Vec p,PetscReal *delta,PetscReal *step1,PetscReal *step2)
    535: {
    536:   PetscReal      zero = 0.0,dsq,ptp,pts,rad,sts;
    538: #if defined(PETSC_USE_COMPLEX)
    539:   PetscScalar    cptp,cpts,csts;
    540: #endif

    543: #if defined(PETSC_USE_COMPLEX)
    544:   VecDot(p,s,&cpts); pts = PetscRealPart(cpts);
    545:   VecDot(p,p,&cptp); ptp = PetscRealPart(cptp);
    546:   VecDot(s,s,&csts); sts = PetscRealPart(csts);
    547: #else
    548:   VecDot(p,s,&pts);
    549:   VecDot(p,p,&ptp);
    550:   VecDot(s,s,&sts);
    551: #endif
    552:   dsq  = (*delta)*(*delta);
    553:   rad  = sqrt((pts*pts) - ptp*(sts - dsq));
    554:   if (pts > zero) {
    555:     *step2 = -(pts + rad)/ptp;
    556:     *step1 = (sts - dsq)/(ptp * *step2);
    557:   } else {
    558:     *step1 = -(pts - rad)/ptp;
    559:     *step2 = (sts - dsq)/(ptp * *step1);
    560:   }
    561:   return(0);
    562: }