Actual source code: qcg.c

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

  3: /*
  4:   KSPQCGQuadraticRoots - Computes the roots of the quadratic,
  5:          ||s + step*p|| - delta = 0
  6:    such that step1 >= 0 >= step2.
  7:    where
  8:       delta:
  9:         On entry delta must contain scalar delta.
 10:         On exit delta is unchanged.
 11:       step1:
 12:         On entry step1 need not be specified.
 13:         On exit step1 contains the non-negative root.
 14:       step2:
 15:         On entry step2 need not be specified.
 16:         On exit step2 contains the non-positive root.
 17:    C code is translated from the Fortran version of the MINPACK-2 Project,
 18:    Argonne National Laboratory, Brett M. Averick and Richard G. Carter.
 19: */
 20: static PetscErrorCode KSPQCGQuadraticRoots(Vec s, Vec p, PetscReal delta, PetscReal *step1, PetscReal *step2)
 21: {
 22:   PetscReal dsq, ptp, pts, rad, sts;

 24:   PetscFunctionBegin;
 25:   PetscCall(VecDotRealPart(p, s, &pts));
 26:   PetscCall(VecDotRealPart(p, p, &ptp));
 27:   PetscCall(VecDotRealPart(s, s, &sts));
 28:   dsq = delta * delta;
 29:   rad = PetscSqrtReal((pts * pts) - ptp * (sts - dsq));
 30:   if (pts > 0.0) {
 31:     *step2 = -(pts + rad) / ptp;
 32:     *step1 = (sts - dsq) / (ptp * *step2);
 33:   } else {
 34:     *step1 = -(pts - rad) / ptp;
 35:     *step2 = (sts - dsq) / (ptp * *step1);
 36:   }
 37:   PetscFunctionReturn(PETSC_SUCCESS);
 38: }

 40: /*@
 41:   KSPQCGSetTrustRegionRadius - Sets the radius of the trust region for `KSPQCG`

 43:   Logically Collective

 45:   Input Parameters:
 46: + ksp   - the iterative context
 47: - delta - the trust region radius (Infinity is the default)

 49:   Options Database Key:
 50: . -ksp_qcg_trustregionradius <delta> - trust region radius

 52:   Level: advanced

 54: .seealso: `KSPQCG`
 55: @*/
 56: PetscErrorCode KSPQCGSetTrustRegionRadius(KSP ksp, PetscReal delta)
 57: {
 58:   PetscFunctionBegin;
 60:   PetscCheck(delta >= 0.0, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "Tolerance must be non-negative");
 61:   PetscTryMethod(ksp, "KSPQCGSetTrustRegionRadius_C", (KSP, PetscReal), (ksp, delta));
 62:   PetscFunctionReturn(PETSC_SUCCESS);
 63: }

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

 69:   Not Collective

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

 74:   Output Parameter:
 75: . tsnorm - the norm

 77:   Level: advanced

 79: .seealso: `KSPQCG`
 80: @*/
 81: PetscErrorCode KSPQCGGetTrialStepNorm(KSP ksp, PetscReal *tsnorm)
 82: {
 83:   PetscFunctionBegin;
 85:   PetscUseMethod(ksp, "KSPQCGGetTrialStepNorm_C", (KSP, PetscReal *), (ksp, tsnorm));
 86:   PetscFunctionReturn(PETSC_SUCCESS);
 87: }

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

 92:   Collective

 94:   Input Parameter:
 95: . ksp - the iterative context

 97:   Output Parameter:
 98: . quadratic - the quadratic function evaluated at the new iterate

100:   Level: advanced

102:   Note:
103:   $ q(s) = g^T * s + 0.5 * s^T * H * s $ which satisfies the Euclidean Norm trust region constraint
104: .vb
105:   || D * s || <= delta,
106: .ve
107:   where
108: .vb
109:   delta is the trust region radius,
110:   g is the gradient vector, and
111:   H is Hessian matrix,
112:   D is a scaling matrix.
113: .ve

115: .seealso: [](ch_ksp), `KSPQCG`
116: @*/
117: PetscErrorCode KSPQCGGetQuadratic(KSP ksp, PetscReal *quadratic)
118: {
119:   PetscFunctionBegin;
121:   PetscUseMethod(ksp, "KSPQCGGetQuadratic_C", (KSP, PetscReal *), (ksp, quadratic));
122:   PetscFunctionReturn(PETSC_SUCCESS);
123: }

125: static PetscErrorCode KSPSolve_QCG(KSP ksp)
126: {
127:   /*
128:    Correpondence with documentation above:
129:       B = g = gradient,
130:       X = s = step
131:    Note:  This is not coded correctly for complex arithmetic!
132:  */

134:   KSP_QCG    *pcgP = (KSP_QCG *)ksp->data;
135:   Mat         Amat, Pmat;
136:   Vec         W, WA, WA2, R, P, ASP, BS, X, B;
137:   PetscScalar scal, beta, rntrn, step;
138:   PetscReal   q1, q2, xnorm, step1, step2, rnrm = 0.0, btx, xtax;
139:   PetscReal   ptasp, rtr, wtasp, bstp;
140:   PetscReal   dzero = 0.0, bsnrm = 0.0;
141:   PetscInt    i, maxit;
142:   PC          pc = ksp->pc;
143:   PetscBool   diagonalscale;

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

150:   ksp->its = 0;
151:   maxit    = ksp->max_it;
152:   WA       = ksp->work[0];
153:   R        = ksp->work[1];
154:   P        = ksp->work[2];
155:   ASP      = ksp->work[3];
156:   BS       = ksp->work[4];
157:   W        = ksp->work[5];
158:   WA2      = ksp->work[6];
159:   X        = ksp->vec_sol;
160:   B        = ksp->vec_rhs;

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

164:   /* Initialize variables */
165:   PetscCall(VecSet(W, 0.0)); /* W = 0 */
166:   PetscCall(VecSet(X, 0.0)); /* X = 0 */
167:   PetscCall(PCGetOperators(pc, &Amat, &Pmat));

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

172:   if (ksp->normtype != KSP_NORM_NONE) {
173:     PetscCall(VecNorm(BS, NORM_2, &bsnrm));
174:     KSPCheckNorm(ksp, bsnrm);
175:   }
176:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
177:   ksp->its   = 0;
178:   ksp->rnorm = bsnrm;
179:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
180:   PetscCall(KSPLogResidualHistory(ksp, bsnrm));
181:   PetscCall(KSPMonitor(ksp, 0, bsnrm));
182:   PetscCall((*ksp->converged)(ksp, 0, bsnrm, &ksp->reason, ksp->cnvP));
183:   if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);

185:   /* Compute the initial scaled direction and scaled residual */
186:   PetscCall(VecCopy(BS, R));
187:   PetscCall(VecScale(R, -1.0));
188:   PetscCall(VecCopy(R, P));
189:   PetscCall(VecDotRealPart(R, R, &rtr));

191:   for (i = 0; i <= maxit; i++) {
192:     PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
193:     ksp->its++;
194:     PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));

196:     /* Compute:  asp = D^{-T}*A*D^{-1}*p  */
197:     PetscCall(PCApplySymmetricRight(pc, P, WA));
198:     PetscCall(KSP_MatMult(ksp, Amat, WA, WA2));
199:     PetscCall(PCApplySymmetricLeft(pc, WA2, ASP));

201:     /* Check for negative curvature */
202:     PetscCall(VecDotRealPart(P, ASP, &ptasp));
203:     if (ptasp <= dzero) {
204:       /* Scaled negative curvature direction:  Compute a step so that
205:         ||w + step*p|| = delta and QS(w + step*p) is least */

207:       if (!i) {
208:         PetscCall(VecCopy(P, X));
209:         PetscCall(VecNorm(X, NORM_2, &xnorm));
210:         KSPCheckNorm(ksp, xnorm);
211:         scal = pcgP->delta / xnorm;
212:         PetscCall(VecScale(X, scal));
213:       } else {
214:         /* Compute roots of quadratic */
215:         PetscCall(KSPQCGQuadraticRoots(W, P, pcgP->delta, &step1, &step2));
216:         PetscCall(VecDotRealPart(W, ASP, &wtasp));
217:         PetscCall(VecDotRealPart(BS, P, &bstp));
218:         PetscCall(VecCopy(W, X));
219:         q1 = step1 * (bstp + wtasp + .5 * step1 * ptasp);
220:         q2 = step2 * (bstp + wtasp + .5 * step2 * ptasp);
221:         if (q1 <= q2) {
222:           PetscCall(VecAXPY(X, step1, P));
223:         } else {
224:           PetscCall(VecAXPY(X, step2, P));
225:         }
226:       }
227:       pcgP->ltsnrm = pcgP->delta; /* convergence in direction of */
228:       ksp->reason  = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
229:       if (!i) {
230:         PetscCall(PetscInfo(ksp, "negative curvature: delta=%g\n", (double)pcgP->delta));
231:       } else {
232:         PetscCall(PetscInfo(ksp, "negative curvature: step1=%g, step2=%g, delta=%g\n", (double)step1, (double)step2, (double)pcgP->delta));
233:       }

235:     } else {
236:       /* Compute step along p */
237:       step = rtr / ptasp;
238:       PetscCall(VecCopy(W, X));       /*  x = w  */
239:       PetscCall(VecAXPY(X, step, P)); /*  x <- step*p + x  */
240:       PetscCall(VecNorm(X, NORM_2, &pcgP->ltsnrm));
241:       KSPCheckNorm(ksp, pcgP->ltsnrm);

243:       if (pcgP->ltsnrm > pcgP->delta) {
244:         /* Since the trial iterate is outside the trust region,
245:             evaluate a constrained step along p so that
246:                     ||w + step*p|| = delta
247:           The positive step is always better in this case. */
248:         if (!i) {
249:           scal = pcgP->delta / pcgP->ltsnrm;
250:           PetscCall(VecScale(X, scal));
251:         } else {
252:           /* Compute roots of quadratic */
253:           PetscCall(KSPQCGQuadraticRoots(W, P, pcgP->delta, &step1, &step2));
254:           PetscCall(VecCopy(W, X));
255:           PetscCall(VecAXPY(X, step1, P)); /*  x <- step1*p + x  */
256:         }
257:         pcgP->ltsnrm = pcgP->delta;
258:         ksp->reason  = KSP_CONVERGED_STEP_LENGTH; /* convergence along constrained step */
259:         if (!i) {
260:           PetscCall(PetscInfo(ksp, "constrained step: delta=%g\n", (double)pcgP->delta));
261:         } else {
262:           PetscCall(PetscInfo(ksp, "constrained step: step1=%g, step2=%g, delta=%g\n", (double)step1, (double)step2, (double)pcgP->delta));
263:         }

265:       } else {
266:         /* Evaluate the current step */
267:         PetscCall(VecCopy(X, W));          /* update interior iterate */
268:         PetscCall(VecAXPY(R, -step, ASP)); /* r <- -step*asp + r */
269:         if (ksp->normtype != KSP_NORM_NONE) {
270:           PetscCall(VecNorm(R, NORM_2, &rnrm));
271:           KSPCheckNorm(ksp, rnrm);
272:         }
273:         PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
274:         ksp->rnorm = rnrm;
275:         PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
276:         PetscCall(KSPLogResidualHistory(ksp, rnrm));
277:         PetscCall(KSPMonitor(ksp, i + 1, rnrm));
278:         PetscCall((*ksp->converged)(ksp, i + 1, rnrm, &ksp->reason, ksp->cnvP));
279:         if (ksp->reason) { /* convergence for */
280:           PetscCall(PetscInfo(ksp, "truncated step: step=%g, rnrm=%g, delta=%g\n", (double)PetscRealPart(step), (double)rnrm, (double)pcgP->delta));
281:         }
282:       }
283:     }
284:     if (ksp->reason) break; /* Convergence has been attained */
285:     else { /* Compute a new AS-orthogonal direction */ PetscCall(VecDot(R, R, &rntrn));
286:       beta = rntrn / rtr;
287:       PetscCall(VecAYPX(P, beta, R)); /*  p <- r + beta*p  */
288:       rtr = PetscRealPart(rntrn);
289:     }
290:   }
291:   if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;

293:   /* Unscale x */
294:   PetscCall(VecCopy(X, WA2));
295:   PetscCall(PCApplySymmetricRight(pc, WA2, X));

297:   PetscCall(KSP_MatMult(ksp, Amat, X, WA));
298:   PetscCall(VecDotRealPart(B, X, &btx));
299:   PetscCall(VecDotRealPart(X, WA, &xtax));

301:   pcgP->quadratic = btx + .5 * xtax;
302:   PetscFunctionReturn(PETSC_SUCCESS);
303: }

305: static PetscErrorCode KSPSetUp_QCG(KSP ksp)
306: {
307:   PetscFunctionBegin;
308:   /* Get work vectors from user code */
309:   PetscCall(KSPSetWorkVecs(ksp, 7));
310:   PetscFunctionReturn(PETSC_SUCCESS);
311: }

313: static PetscErrorCode KSPDestroy_QCG(KSP ksp)
314: {
315:   PetscFunctionBegin;
316:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGGetQuadratic_C", NULL));
317:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGGetTrialStepNorm_C", NULL));
318:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGSetTrustRegionRadius_C", NULL));
319:   PetscCall(KSPDestroyDefault(ksp));
320:   PetscFunctionReturn(PETSC_SUCCESS);
321: }

323: static PetscErrorCode KSPQCGSetTrustRegionRadius_QCG(KSP ksp, PetscReal delta)
324: {
325:   KSP_QCG *cgP = (KSP_QCG *)ksp->data;

327:   PetscFunctionBegin;
328:   cgP->delta = delta;
329:   PetscFunctionReturn(PETSC_SUCCESS);
330: }

332: static PetscErrorCode KSPQCGGetTrialStepNorm_QCG(KSP ksp, PetscReal *ltsnrm)
333: {
334:   KSP_QCG *cgP = (KSP_QCG *)ksp->data;

336:   PetscFunctionBegin;
337:   *ltsnrm = cgP->ltsnrm;
338:   PetscFunctionReturn(PETSC_SUCCESS);
339: }

341: static PetscErrorCode KSPQCGGetQuadratic_QCG(KSP ksp, PetscReal *quadratic)
342: {
343:   KSP_QCG *cgP = (KSP_QCG *)ksp->data;

345:   PetscFunctionBegin;
346:   *quadratic = cgP->quadratic;
347:   PetscFunctionReturn(PETSC_SUCCESS);
348: }

350: static PetscErrorCode KSPSetFromOptions_QCG(KSP ksp, PetscOptionItems *PetscOptionsObject)
351: {
352:   PetscReal delta;
353:   KSP_QCG  *cgP = (KSP_QCG *)ksp->data;
354:   PetscBool flg;

356:   PetscFunctionBegin;
357:   PetscOptionsHeadBegin(PetscOptionsObject, "KSP QCG Options");
358:   PetscCall(PetscOptionsReal("-ksp_qcg_trustregionradius", "Trust Region Radius", "KSPQCGSetTrustRegionRadius", cgP->delta, &delta, &flg));
359:   if (flg) PetscCall(KSPQCGSetTrustRegionRadius(ksp, delta));
360:   PetscOptionsHeadEnd();
361:   PetscFunctionReturn(PETSC_SUCCESS);
362: }

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

367:    Options Database Key:
368: .      -ksp_qcg_trustregionradius <r> - Trust Region Radius

370:    Level: developer

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

375:     Uses preconditioned conjugate gradient to compute
376:       an approximate minimizer of the quadratic function $ q(s) = g^T * s + .5 * s^T * H * s $   subject to the Euclidean norm trust region constraint
377: .vb
378:             || D * s || <= delta,
379: .ve
380:    where
381: .vb
382:      delta is the trust region radius,
383:      g is the gradient vector, and
384:      H is Hessian matrix,
385:      D is a scaling matrix.
386: .ve

388:    `KSPConvergedReason` may be
389: .vb
390:    KSP_CONVERGED_NEG_CURVE if convergence is reached along a negative curvature direction,
391:    KSP_CONVERGED_STEP_LENGTH if convergence is reached along a constrained step,
392: .ve
393:    and other `KSP` converged/diverged reasons

395:   Notes:
396:   Currently we allow symmetric preconditioning with the following scaling matrices:
397: .vb
398:       `PCNONE`:   D = Identity matrix
399:       `PCJACOBI`: D = diag [d_1, d_2, ...., d_n], where d_i = sqrt(H[i,i])
400:       `PCICC`:    D = L^T, implemented with forward and backward solves. Here L is an incomplete Cholesky factor of H.
401: .ve

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

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

411: PETSC_EXTERN PetscErrorCode KSPCreate_QCG(KSP ksp)
412: {
413:   KSP_QCG *cgP;

415:   PetscFunctionBegin;
416:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_SYMMETRIC, 3));
417:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_SYMMETRIC, 1));
418:   PetscCall(KSPSetConvergedNegativeCurvature(ksp, PETSC_TRUE));
419:   PetscCall(PetscNew(&cgP));

421:   ksp->data                = (void *)cgP;
422:   ksp->ops->setup          = KSPSetUp_QCG;
423:   ksp->ops->setfromoptions = KSPSetFromOptions_QCG;
424:   ksp->ops->solve          = KSPSolve_QCG;
425:   ksp->ops->destroy        = KSPDestroy_QCG;
426:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
427:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
428:   ksp->ops->view           = NULL;

430:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGGetQuadratic_C", KSPQCGGetQuadratic_QCG));
431:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGGetTrialStepNorm_C", KSPQCGGetTrialStepNorm_QCG));
432:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPQCGSetTrustRegionRadius_C", KSPQCGSetTrustRegionRadius_QCG));
433:   cgP->delta = PETSC_MAX_REAL; /* default trust region radius is infinite */
434:   PetscFunctionReturn(PETSC_SUCCESS);
435: }