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 for `KSPQCG`

  9:     Logically Collective

 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> - trust region radius

 18:     Level: advanced

 20: @*/
 21: PetscErrorCode KSPQCGSetTrustRegionRadius(KSP ksp, PetscReal delta)
 22: {
 23:   PetscFunctionBegin;
 25:   PetscCheck(delta >= 0.0, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "Tolerance must be non-negative");
 26:   PetscTryMethod(ksp, "KSPQCGSetTrustRegionRadius_C", (KSP, PetscReal), (ksp, delta));
 27:   PetscFunctionReturn(PETSC_SUCCESS);
 28: }

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

 34:     Not Collective

 36:     Input Parameter:
 37: .   ksp - the iterative context

 39:     Output Parameter:
 40: .   tsnorm - the norm

 42:     Level: advanced
 43: @*/
 44: PetscErrorCode KSPQCGGetTrialStepNorm(KSP ksp, PetscReal *tsnorm)
 45: {
 46:   PetscFunctionBegin;
 48:   PetscUseMethod(ksp, "KSPQCGGetTrialStepNorm_C", (KSP, PetscReal *), (ksp, tsnorm));
 49:   PetscFunctionReturn(PETSC_SUCCESS);
 50: }

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

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

 57:     which satisfies the Euclidean Norm trust region constraint

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

 61:     where

 63:      delta is the trust region radius,
 64:      g is the gradient vector, and
 65:      H is Hessian matrix,
 66:      D is a scaling matrix.

 68:     Collective

 70:     Input Parameter:
 71: .   ksp - the iterative context

 73:     Output Parameter:
 74: .   quadratic - the quadratic function evaluated at the new iterate

 76:     Level: advanced

 78: .seealso: [](ch_ksp), `KSPQCG`
 79: @*/
 80: PetscErrorCode KSPQCGGetQuadratic(KSP ksp, PetscReal *quadratic)
 81: {
 82:   PetscFunctionBegin;
 84:   PetscUseMethod(ksp, "KSPQCGGetQuadratic_C", (KSP, PetscReal *), (ksp, quadratic));
 85:   PetscFunctionReturn(PETSC_SUCCESS);
 86: }

 88: PetscErrorCode KSPSolve_QCG(KSP ksp)
 89: {
 90:   /*
 91:    Correpondence with documentation above:
 92:       B = g = gradient,
 93:       X = s = step
 94:    Note:  This is not coded correctly for complex arithmetic!
 95:  */

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

108:   PetscFunctionBegin;
109:   PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
110:   PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
111:   PetscCheck(!ksp->transpose_solve, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Currently does not support transpose solve");

113:   ksp->its = 0;
114:   maxit    = ksp->max_it;
115:   WA       = ksp->work[0];
116:   R        = ksp->work[1];
117:   P        = ksp->work[2];
118:   ASP      = ksp->work[3];
119:   BS       = ksp->work[4];
120:   W        = ksp->work[5];
121:   WA2      = ksp->work[6];
122:   X        = ksp->vec_sol;
123:   B        = ksp->vec_rhs;

125:   PetscCheck(pcgP->delta > dzero, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "Input error: delta <= 0");

127:   /* Initialize variables */
128:   PetscCall(VecSet(W, 0.0)); /* W = 0 */
129:   PetscCall(VecSet(X, 0.0)); /* X = 0 */
130:   PetscCall(PCGetOperators(pc, &Amat, &Pmat));

132:   /* Compute:  BS = D^{-1} B */
133:   PetscCall(PCApplySymmetricLeft(pc, B, BS));

135:   if (ksp->normtype != KSP_NORM_NONE) {
136:     PetscCall(VecNorm(BS, NORM_2, &bsnrm));
137:     KSPCheckNorm(ksp, bsnrm);
138:   }
139:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
140:   ksp->its   = 0;
141:   ksp->rnorm = bsnrm;
142:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
143:   PetscCall(KSPLogResidualHistory(ksp, bsnrm));
144:   PetscCall(KSPMonitor(ksp, 0, bsnrm));
145:   PetscCall((*ksp->converged)(ksp, 0, bsnrm, &ksp->reason, ksp->cnvP));
146:   if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);

148:   /* Compute the initial scaled direction and scaled residual */
149:   PetscCall(VecCopy(BS, R));
150:   PetscCall(VecScale(R, -1.0));
151:   PetscCall(VecCopy(R, P));
152:   PetscCall(VecDotRealPart(R, R, &rtr));

154:   for (i = 0; i <= maxit; i++) {
155:     PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
156:     ksp->its++;
157:     PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));

159:     /* Compute:  asp = D^{-T}*A*D^{-1}*p  */
160:     PetscCall(PCApplySymmetricRight(pc, P, WA));
161:     PetscCall(KSP_MatMult(ksp, Amat, WA, WA2));
162:     PetscCall(PCApplySymmetricLeft(pc, WA2, ASP));

164:     /* Check for negative curvature */
165:     PetscCall(VecDotRealPart(P, ASP, &ptasp));
166:     if (ptasp <= dzero) {
167:       /* Scaled negative curvature direction:  Compute a step so that
168:         ||w + step*p|| = delta and QS(w + step*p) is least */

170:       if (!i) {
171:         PetscCall(VecCopy(P, X));
172:         PetscCall(VecNorm(X, NORM_2, &xnorm));
173:         KSPCheckNorm(ksp, xnorm);
174:         scal = pcgP->delta / xnorm;
175:         PetscCall(VecScale(X, scal));
176:       } else {
177:         /* Compute roots of quadratic */
178:         PetscCall(KSPQCGQuadraticRoots(W, P, pcgP->delta, &step1, &step2));
179:         PetscCall(VecDotRealPart(W, ASP, &wtasp));
180:         PetscCall(VecDotRealPart(BS, P, &bstp));
181:         PetscCall(VecCopy(W, X));
182:         q1 = step1 * (bstp + wtasp + .5 * step1 * ptasp);
183:         q2 = step2 * (bstp + wtasp + .5 * step2 * ptasp);
184:         if (q1 <= q2) {
185:           PetscCall(VecAXPY(X, step1, P));
186:         } else {
187:           PetscCall(VecAXPY(X, step2, P));
188:         }
189:       }
190:       pcgP->ltsnrm = pcgP->delta; /* convergence in direction of */
191:       ksp->reason  = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
192:       if (!i) {
193:         PetscCall(PetscInfo(ksp, "negative curvature: delta=%g\n", (double)pcgP->delta));
194:       } else {
195:         PetscCall(PetscInfo(ksp, "negative curvature: step1=%g, step2=%g, delta=%g\n", (double)step1, (double)step2, (double)pcgP->delta));
196:       }

198:     } else {
199:       /* Compute step along p */
200:       step = rtr / ptasp;
201:       PetscCall(VecCopy(W, X));       /*  x = w  */
202:       PetscCall(VecAXPY(X, step, P)); /*  x <- step*p + x  */
203:       PetscCall(VecNorm(X, NORM_2, &pcgP->ltsnrm));
204:       KSPCheckNorm(ksp, pcgP->ltsnrm);

206:       if (pcgP->ltsnrm > pcgP->delta) {
207:         /* Since the trial iterate is outside the trust region,
208:             evaluate a constrained step along p so that
209:                     ||w + step*p|| = delta
210:           The positive step is always better in this case. */
211:         if (!i) {
212:           scal = pcgP->delta / pcgP->ltsnrm;
213:           PetscCall(VecScale(X, scal));
214:         } else {
215:           /* Compute roots of quadratic */
216:           PetscCall(KSPQCGQuadraticRoots(W, P, pcgP->delta, &step1, &step2));
217:           PetscCall(VecCopy(W, X));
218:           PetscCall(VecAXPY(X, step1, P)); /*  x <- step1*p + x  */
219:         }
220:         pcgP->ltsnrm = pcgP->delta;
221:         ksp->reason  = KSP_CONVERGED_STEP_LENGTH; /* convergence along constrained step */
222:         if (!i) {
223:           PetscCall(PetscInfo(ksp, "constrained step: delta=%g\n", (double)pcgP->delta));
224:         } else {
225:           PetscCall(PetscInfo(ksp, "constrained step: step1=%g, step2=%g, delta=%g\n", (double)step1, (double)step2, (double)pcgP->delta));
226:         }

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

256:   /* Unscale x */
257:   PetscCall(VecCopy(X, WA2));
258:   PetscCall(PCApplySymmetricRight(pc, WA2, X));

260:   PetscCall(KSP_MatMult(ksp, Amat, X, WA));
261:   PetscCall(VecDotRealPart(B, X, &btx));
262:   PetscCall(VecDotRealPart(X, WA, &xtax));

264:   pcgP->quadratic = btx + .5 * xtax;
265:   PetscFunctionReturn(PETSC_SUCCESS);
266: }

268: PetscErrorCode KSPSetUp_QCG(KSP ksp)
269: {
270:   PetscFunctionBegin;
271:   /* Get work vectors from user code */
272:   PetscCall(KSPSetWorkVecs(ksp, 7));
273:   PetscFunctionReturn(PETSC_SUCCESS);
274: }

276: PetscErrorCode KSPDestroy_QCG(KSP ksp)
277: {
278:   PetscFunctionBegin;
279:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGGetQuadratic_C", NULL));
280:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGGetTrialStepNorm_C", NULL));
281:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGSetTrustRegionRadius_C", NULL));
282:   PetscCall(KSPDestroyDefault(ksp));
283:   PetscFunctionReturn(PETSC_SUCCESS);
284: }

286: static PetscErrorCode KSPQCGSetTrustRegionRadius_QCG(KSP ksp, PetscReal delta)
287: {
288:   KSP_QCG *cgP = (KSP_QCG *)ksp->data;

290:   PetscFunctionBegin;
291:   cgP->delta = delta;
292:   PetscFunctionReturn(PETSC_SUCCESS);
293: }

295: static PetscErrorCode KSPQCGGetTrialStepNorm_QCG(KSP ksp, PetscReal *ltsnrm)
296: {
297:   KSP_QCG *cgP = (KSP_QCG *)ksp->data;

299:   PetscFunctionBegin;
300:   *ltsnrm = cgP->ltsnrm;
301:   PetscFunctionReturn(PETSC_SUCCESS);
302: }

304: static PetscErrorCode KSPQCGGetQuadratic_QCG(KSP ksp, PetscReal *quadratic)
305: {
306:   KSP_QCG *cgP = (KSP_QCG *)ksp->data;

308:   PetscFunctionBegin;
309:   *quadratic = cgP->quadratic;
310:   PetscFunctionReturn(PETSC_SUCCESS);
311: }

313: PetscErrorCode KSPSetFromOptions_QCG(KSP ksp, PetscOptionItems *PetscOptionsObject)
314: {
315:   PetscReal delta;
316:   KSP_QCG  *cgP = (KSP_QCG *)ksp->data;
317:   PetscBool flg;

319:   PetscFunctionBegin;
320:   PetscOptionsHeadBegin(PetscOptionsObject, "KSP QCG Options");
321:   PetscCall(PetscOptionsReal("-ksp_qcg_trustregionradius", "Trust Region Radius", "KSPQCGSetTrustRegionRadius", cgP->delta, &delta, &flg));
322:   if (flg) PetscCall(KSPQCGSetTrustRegionRadius(ksp, delta));
323:   PetscOptionsHeadEnd();
324:   PetscFunctionReturn(PETSC_SUCCESS);
325: }

327: /*MC
328:      KSPQCG -   Code to run conjugate gradient method subject to a constraint on the solution norm.

330:    Options Database Key:
331: .      -ksp_qcg_trustregionradius <r> - Trust Region Radius

333:    Level: developer

335:    Notes:
336:     This is rarely used directly, ir is used in Trust Region methods for nonlinear equations, `SNESNEWTONTR`

338:     Uses preconditioned conjugate gradient to compute
339:       an approximate minimizer of the quadratic function

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

343:    subject to the Euclidean norm trust region constraint

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

347:    where

349:      delta is the trust region radius,
350:      g is the gradient vector, and
351:      H is Hessian matrix,
352:      D is a scaling matrix.

354:    `KSPConvergedReason` may be
355: .vb
356:    KSP_CONVERGED_NEG_CURVE if convergence is reached along a negative curvature direction,
357:    KSP_CONVERGED_STEP_LENGTH if convergence is reached along a constrained step,
358: .ve
359:    and other `KSP` converged/diverged reasons

361:   Notes:
362:   Currently we allow symmetric preconditioning with the following scaling matrices:
363: .vb
364:       `PCNONE`:   D = Identity matrix
365:       `PCJACOBI`: D = diag [d_1, d_2, ...., d_n], where d_i = sqrt(H[i,i])
366:       `PCICC`:    D = L^T, implemented with forward and backward solves. Here L is an incomplete Cholesky factor of H.
367: .ve

369:   References:
370: . * - Trond Steihaug, The Conjugate Gradient Method and Trust Regions in Large Scale Optimization,
371:    SIAM Journal on Numerical Analysis, Vol. 20, No. 3 (Jun., 1983).

373: .seealso: [](ch_ksp), 'KSPNASH`, `KSPGLTR`, `KSPSTCG`, `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPQCGSetTrustRegionRadius()`
374:           `KSPQCGGetTrialStepNorm()`, `KSPQCGGetQuadratic()`
375: M*/

377: PETSC_EXTERN PetscErrorCode KSPCreate_QCG(KSP ksp)
378: {
379:   KSP_QCG *cgP;

381:   PetscFunctionBegin;
382:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_SYMMETRIC, 3));
383:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_SYMMETRIC, 1));
384:   PetscCall(KSPSetConvergedNegativeCurvature(ksp, PETSC_TRUE));
385:   PetscCall(PetscNew(&cgP));

387:   ksp->data                = (void *)cgP;
388:   ksp->ops->setup          = KSPSetUp_QCG;
389:   ksp->ops->setfromoptions = KSPSetFromOptions_QCG;
390:   ksp->ops->solve          = KSPSolve_QCG;
391:   ksp->ops->destroy        = KSPDestroy_QCG;
392:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
393:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
394:   ksp->ops->view           = NULL;

396:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGGetQuadratic_C", KSPQCGGetQuadratic_QCG));
397:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGGetTrialStepNorm_C", KSPQCGGetTrialStepNorm_QCG));
398:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGSetTrustRegionRadius_C", KSPQCGSetTrustRegionRadius_QCG));
399:   cgP->delta = PETSC_MAX_REAL; /* default trust region radius is infinite */
400:   PetscFunctionReturn(PETSC_SUCCESS);
401: }

403: /* ---------------------------------------------------------- */
404: /*
405:   KSPQCGQuadraticRoots - Computes the roots of the quadratic,
406:          ||s + step*p|| - delta = 0
407:    such that step1 >= 0 >= step2.
408:    where
409:       delta:
410:         On entry delta must contain scalar delta.
411:         On exit delta is unchanged.
412:       step1:
413:         On entry step1 need not be specified.
414:         On exit step1 contains the non-negative root.
415:       step2:
416:         On entry step2 need not be specified.
417:         On exit step2 contains the non-positive root.
418:    C code is translated from the Fortran version of the MINPACK-2 Project,
419:    Argonne National Laboratory, Brett M. Averick and Richard G. Carter.
420: */
421: static PetscErrorCode KSPQCGQuadraticRoots(Vec s, Vec p, PetscReal delta, PetscReal *step1, PetscReal *step2)
422: {
423:   PetscReal dsq, ptp, pts, rad, sts;

425:   PetscFunctionBegin;
426:   PetscCall(VecDotRealPart(p, s, &pts));
427:   PetscCall(VecDotRealPart(p, p, &ptp));
428:   PetscCall(VecDotRealPart(s, s, &sts));
429:   dsq = delta * delta;
430:   rad = PetscSqrtReal((pts * pts) - ptp * (sts - dsq));
431:   if (pts > 0.0) {
432:     *step2 = -(pts + rad) / ptp;
433:     *step1 = (sts - dsq) / (ptp * *step2);
434:   } else {
435:     *step1 = -(pts - rad) / ptp;
436:     *step2 = (sts - dsq) / (ptp * *step1);
437:   }
438:   PetscFunctionReturn(PETSC_SUCCESS);
439: }