Actual source code: lsqr.c
2: /* lourens.vanzanen@shell.com contributed the standard error estimates of the solution, Jul 25, 2006 */
3: /* Bas van't Hof contributed the preconditioned aspects Feb 10, 2010 */
5: #define SWAP(a, b, c) \
6: { \
7: c = a; \
8: a = b; \
9: b = c; \
10: }
12: #include <petsc/private/kspimpl.h>
13: #include <petscdraw.h>
15: typedef struct {
16: PetscInt nwork_n, nwork_m;
17: Vec *vwork_m; /* work vectors of length m, where the system is size m x n */
18: Vec *vwork_n; /* work vectors of length n */
19: Vec se; /* Optional standard error vector */
20: PetscBool se_flg; /* flag for -ksp_lsqr_set_standard_error */
21: PetscBool exact_norm; /* flag for -ksp_lsqr_exact_mat_norm */
22: PetscReal arnorm; /* Good estimate of norm((A*inv(Pmat))'*r), where r = A*x - b, used in specific stopping criterion */
23: PetscReal anorm; /* Poor estimate of norm(A*inv(Pmat),'fro') used in specific stopping criterion */
24: /* Backup previous convergence test */
25: PetscErrorCode (*converged)(KSP, PetscInt, PetscReal, KSPConvergedReason *, void *);
26: PetscErrorCode (*convergeddestroy)(void *);
27: void *cnvP;
28: } KSP_LSQR;
30: static PetscErrorCode VecSquare(Vec v)
31: {
32: PetscScalar *x;
33: PetscInt i, n;
35: VecGetLocalSize(v, &n);
36: VecGetArray(v, &x);
37: for (i = 0; i < n; i++) x[i] *= PetscConj(x[i]);
38: VecRestoreArray(v, &x);
39: return 0;
40: }
42: static PetscErrorCode KSPSetUp_LSQR(KSP ksp)
43: {
44: KSP_LSQR *lsqr = (KSP_LSQR *)ksp->data;
45: PetscBool nopreconditioner;
47: PetscObjectTypeCompare((PetscObject)ksp->pc, PCNONE, &nopreconditioner);
49: if (lsqr->vwork_m) VecDestroyVecs(lsqr->nwork_m, &lsqr->vwork_m);
51: if (lsqr->vwork_n) VecDestroyVecs(lsqr->nwork_n, &lsqr->vwork_n);
53: lsqr->nwork_m = 2;
54: if (nopreconditioner) lsqr->nwork_n = 4;
55: else lsqr->nwork_n = 5;
56: KSPCreateVecs(ksp, lsqr->nwork_n, &lsqr->vwork_n, lsqr->nwork_m, &lsqr->vwork_m);
58: if (lsqr->se_flg && !lsqr->se) {
59: VecDuplicate(lsqr->vwork_n[0], &lsqr->se);
60: VecSet(lsqr->se, PETSC_INFINITY);
61: } else if (!lsqr->se_flg) {
62: VecDestroy(&lsqr->se);
63: }
64: return 0;
65: }
67: static PetscErrorCode KSPSolve_LSQR(KSP ksp)
68: {
69: PetscInt i, size1, size2;
70: PetscScalar rho, rhobar, phi, phibar, theta, c, s, tmp, tau;
71: PetscReal beta, alpha, rnorm;
72: Vec X, B, V, V1, U, U1, TMP, W, W2, Z = NULL;
73: Mat Amat, Pmat;
74: KSP_LSQR *lsqr = (KSP_LSQR *)ksp->data;
75: PetscBool diagonalscale, nopreconditioner;
77: PCGetDiagonalScale(ksp->pc, &diagonalscale);
80: PCGetOperators(ksp->pc, &Amat, &Pmat);
81: PetscObjectTypeCompare((PetscObject)ksp->pc, PCNONE, &nopreconditioner);
83: /* vectors of length m, where system size is mxn */
84: B = ksp->vec_rhs;
85: U = lsqr->vwork_m[0];
86: U1 = lsqr->vwork_m[1];
88: /* vectors of length n */
89: X = ksp->vec_sol;
90: W = lsqr->vwork_n[0];
91: V = lsqr->vwork_n[1];
92: V1 = lsqr->vwork_n[2];
93: W2 = lsqr->vwork_n[3];
94: if (!nopreconditioner) Z = lsqr->vwork_n[4];
96: /* standard error vector */
97: if (lsqr->se) VecSet(lsqr->se, 0.0);
99: /* Compute initial residual, temporarily use work vector u */
100: if (!ksp->guess_zero) {
101: KSP_MatMult(ksp, Amat, X, U); /* u <- b - Ax */
102: VecAYPX(U, -1.0, B);
103: } else {
104: VecCopy(B, U); /* u <- b (x is 0) */
105: }
107: /* Test for nothing to do */
108: VecNorm(U, NORM_2, &rnorm);
109: KSPCheckNorm(ksp, rnorm);
110: PetscObjectSAWsTakeAccess((PetscObject)ksp);
111: ksp->its = 0;
112: ksp->rnorm = rnorm;
113: PetscObjectSAWsGrantAccess((PetscObject)ksp);
114: KSPLogResidualHistory(ksp, rnorm);
115: KSPMonitor(ksp, 0, rnorm);
116: (*ksp->converged)(ksp, 0, rnorm, &ksp->reason, ksp->cnvP);
117: if (ksp->reason) return 0;
119: beta = rnorm;
120: VecScale(U, 1.0 / beta);
121: KSP_MatMultHermitianTranspose(ksp, Amat, U, V);
122: if (nopreconditioner) {
123: VecNorm(V, NORM_2, &alpha);
124: KSPCheckNorm(ksp, rnorm);
125: } else {
126: /* this is an application of the preconditioner for the normal equations; not the operator, see the manual page */
127: PCApply(ksp->pc, V, Z);
128: VecDotRealPart(V, Z, &alpha);
129: if (alpha <= 0.0) {
130: ksp->reason = KSP_DIVERGED_BREAKDOWN;
131: return 0;
132: }
133: alpha = PetscSqrtReal(alpha);
134: VecScale(Z, 1.0 / alpha);
135: }
136: VecScale(V, 1.0 / alpha);
138: if (nopreconditioner) {
139: VecCopy(V, W);
140: } else {
141: VecCopy(Z, W);
142: }
144: if (lsqr->exact_norm) {
145: MatNorm(Amat, NORM_FROBENIUS, &lsqr->anorm);
146: } else lsqr->anorm = 0.0;
148: lsqr->arnorm = alpha * beta;
149: phibar = beta;
150: rhobar = alpha;
151: i = 0;
152: do {
153: if (nopreconditioner) {
154: KSP_MatMult(ksp, Amat, V, U1);
155: } else {
156: KSP_MatMult(ksp, Amat, Z, U1);
157: }
158: VecAXPY(U1, -alpha, U);
159: VecNorm(U1, NORM_2, &beta);
160: KSPCheckNorm(ksp, beta);
161: if (beta > 0.0) {
162: VecScale(U1, 1.0 / beta); /* beta*U1 = Amat*V - alpha*U */
163: if (!lsqr->exact_norm) lsqr->anorm = PetscSqrtReal(PetscSqr(lsqr->anorm) + PetscSqr(alpha) + PetscSqr(beta));
164: }
166: KSP_MatMultHermitianTranspose(ksp, Amat, U1, V1);
167: VecAXPY(V1, -beta, V);
168: if (nopreconditioner) {
169: VecNorm(V1, NORM_2, &alpha);
170: KSPCheckNorm(ksp, alpha);
171: } else {
172: PCApply(ksp->pc, V1, Z);
173: VecDotRealPart(V1, Z, &alpha);
174: if (alpha <= 0.0) {
175: ksp->reason = KSP_DIVERGED_BREAKDOWN;
176: break;
177: }
178: alpha = PetscSqrtReal(alpha);
179: VecScale(Z, 1.0 / alpha);
180: }
181: VecScale(V1, 1.0 / alpha); /* alpha*V1 = Amat^T*U1 - beta*V */
182: rho = PetscSqrtScalar(rhobar * rhobar + beta * beta);
183: c = rhobar / rho;
184: s = beta / rho;
185: theta = s * alpha;
186: rhobar = -c * alpha;
187: phi = c * phibar;
188: phibar = s * phibar;
189: tau = s * phi;
191: VecAXPY(X, phi / rho, W); /* x <- x + (phi/rho) w */
193: if (lsqr->se) {
194: VecCopy(W, W2);
195: VecSquare(W2);
196: VecScale(W2, 1.0 / (rho * rho));
197: VecAXPY(lsqr->se, 1.0, W2); /* lsqr->se <- lsqr->se + (w^2/rho^2) */
198: }
199: if (nopreconditioner) {
200: VecAYPX(W, -theta / rho, V1); /* w <- v - (theta/rho) w */
201: } else {
202: VecAYPX(W, -theta / rho, Z); /* w <- z - (theta/rho) w */
203: }
205: lsqr->arnorm = alpha * PetscAbsScalar(tau);
206: rnorm = PetscRealPart(phibar);
208: PetscObjectSAWsTakeAccess((PetscObject)ksp);
209: ksp->its++;
210: ksp->rnorm = rnorm;
211: PetscObjectSAWsGrantAccess((PetscObject)ksp);
212: KSPLogResidualHistory(ksp, rnorm);
213: KSPMonitor(ksp, i + 1, rnorm);
214: (*ksp->converged)(ksp, i + 1, rnorm, &ksp->reason, ksp->cnvP);
215: if (ksp->reason) break;
216: SWAP(U1, U, TMP);
217: SWAP(V1, V, TMP);
219: i++;
220: } while (i < ksp->max_it);
221: if (i >= ksp->max_it && !ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
223: /* Finish off the standard error estimates */
224: if (lsqr->se) {
225: tmp = 1.0;
226: MatGetSize(Amat, &size1, &size2);
227: if (size1 > size2) tmp = size1 - size2;
228: tmp = rnorm / PetscSqrtScalar(tmp);
229: VecSqrtAbs(lsqr->se);
230: VecScale(lsqr->se, tmp);
231: }
232: return 0;
233: }
235: PetscErrorCode KSPDestroy_LSQR(KSP ksp)
236: {
237: KSP_LSQR *lsqr = (KSP_LSQR *)ksp->data;
239: /* Free work vectors */
240: if (lsqr->vwork_n) VecDestroyVecs(lsqr->nwork_n, &lsqr->vwork_n);
241: if (lsqr->vwork_m) VecDestroyVecs(lsqr->nwork_m, &lsqr->vwork_m);
242: VecDestroy(&lsqr->se);
243: /* Revert convergence test */
244: KSPSetConvergenceTest(ksp, lsqr->converged, lsqr->cnvP, lsqr->convergeddestroy);
245: /* Free the KSP_LSQR context */
246: PetscFree(ksp->data);
247: PetscObjectComposeFunction((PetscObject)ksp, "KSPLSQRMonitorResidual_C", NULL);
248: PetscObjectComposeFunction((PetscObject)ksp, "KSPLSQRMonitorResidualDrawLG_C", NULL);
249: return 0;
250: }
252: /*@
253: KSPLSQRSetComputeStandardErrorVec - Compute a vector of standard error estimates during `KSPSolve()` for `KSPLSQR`.
255: Logically Collective
257: Input Parameters:
258: + ksp - iterative context
259: - flg - compute the vector of standard estimates or not
261: Level: intermediate
263: Developer Note:
264: Vaclav: I'm not sure whether this vector is useful for anything.
266: .seealso: [](chapter_ksp), `KSPSolve()`, `KSPLSQR`, `KSPLSQRGetStandardErrorVec()`
267: @*/
268: PetscErrorCode KSPLSQRSetComputeStandardErrorVec(KSP ksp, PetscBool flg)
269: {
270: KSP_LSQR *lsqr = (KSP_LSQR *)ksp->data;
272: lsqr->se_flg = flg;
273: return 0;
274: }
276: /*@
277: KSPLSQRSetExactMatNorm - Compute exact matrix norm instead of iteratively refined estimate.
279: Not Collective
281: Input Parameters:
282: + ksp - iterative context
283: - flg - compute exact matrix norm or not
285: Level: intermediate
287: Notes:
288: By default, flg = `PETSC_FALSE`. This is usually preferred to avoid possibly expensive computation of the norm.
289: For flg = `PETSC_TRUE`, we call `MatNorm`(Amat,`NORM_FROBENIUS`,&lsqr->anorm) which will work only for some types of explicitly assembled matrices.
290: This can affect convergence rate as `KSPLSQRConvergedDefault()` assumes different value of ||A|| used in normal equation stopping criterion.
292: .seealso: [](chapter_ksp), `KSPSolve()`, `KSPLSQR`, `KSPLSQRGetNorms()`, `KSPLSQRConvergedDefault()`
293: @*/
294: PetscErrorCode KSPLSQRSetExactMatNorm(KSP ksp, PetscBool flg)
295: {
296: KSP_LSQR *lsqr = (KSP_LSQR *)ksp->data;
298: lsqr->exact_norm = flg;
299: return 0;
300: }
302: /*@
303: KSPLSQRGetStandardErrorVec - Get vector of standard error estimates.
304: Only available if -ksp_lsqr_set_standard_error was set to true
305: or `KSPLSQRSetComputeStandardErrorVec`(ksp, `PETSC_TRUE`) was called.
306: Otherwise returns NULL.
308: Not Collective
310: Input Parameters:
311: . ksp - iterative context
313: Output Parameters:
314: . se - vector of standard estimates
316: Level: intermediate
318: Developer Note:
319: Vaclav: I'm not sure whether this vector is useful for anything.
321: .seealso: [](chapter_ksp), `KSPSolve()`, `KSPLSQR`, `KSPLSQRSetComputeStandardErrorVec()`
322: @*/
323: PetscErrorCode KSPLSQRGetStandardErrorVec(KSP ksp, Vec *se)
324: {
325: KSP_LSQR *lsqr = (KSP_LSQR *)ksp->data;
327: *se = lsqr->se;
328: return 0;
329: }
331: /*@
332: KSPLSQRGetNorms - Get the norm estimates that `KSPLSQR` computes internally during `KSPSolve()`.
334: Not Collective
336: Input Parameter:
337: . ksp - iterative context
339: Output Parameters:
340: + arnorm - good estimate of norm((A*inv(Pmat))'*r), where r = A*x - b, used in specific stopping criterion
341: - anorm - poor estimate of norm(A*inv(Pmat),'fro') used in specific stopping criterion
343: Notes:
344: Output parameters are meaningful only after `KSPSolve()`.
346: These are the same quantities as normar and norma in MATLAB's `lsqr()`, whose output lsvec is a vector of normar / norma for all iterations.
348: If -ksp_lsqr_exact_mat_norm is set or `KSPLSQRSetExactMatNorm`(ksp, `PETSC_TRUE`) called, then anorm is the exact Frobenius norm.
350: Level: intermediate
352: .seealso: [](chapter_ksp), `KSPSolve()`, `KSPLSQR`, `KSPLSQRSetExactMatNorm()`
353: @*/
354: PetscErrorCode KSPLSQRGetNorms(KSP ksp, PetscReal *arnorm, PetscReal *anorm)
355: {
356: KSP_LSQR *lsqr = (KSP_LSQR *)ksp->data;
358: if (arnorm) *arnorm = lsqr->arnorm;
359: if (anorm) *anorm = lsqr->anorm;
360: return 0;
361: }
363: PetscErrorCode KSPLSQRMonitorResidual_LSQR(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
364: {
365: KSP_LSQR *lsqr = (KSP_LSQR *)ksp->data;
366: PetscViewer viewer = vf->viewer;
367: PetscViewerFormat format = vf->format;
368: char normtype[256];
369: PetscInt tablevel;
370: const char *prefix;
372: PetscObjectGetTabLevel((PetscObject)ksp, &tablevel);
373: PetscObjectGetOptionsPrefix((PetscObject)ksp, &prefix);
374: PetscStrncpy(normtype, KSPNormTypes[ksp->normtype], sizeof(normtype));
375: PetscStrtolower(normtype);
376: PetscViewerPushFormat(viewer, format);
377: PetscViewerASCIIAddTab(viewer, tablevel);
378: if (n == 0 && prefix) PetscViewerASCIIPrintf(viewer, " Residual norm, norm of normal equations, and matrix norm for %s solve.\n", prefix);
379: if (!n) {
380: PetscViewerASCIIPrintf(viewer, "%3" PetscInt_FMT " KSP resid norm %14.12e\n", n, (double)rnorm);
381: } else {
382: PetscViewerASCIIPrintf(viewer, "%3" PetscInt_FMT " KSP resid norm %14.12e normal eq resid norm %14.12e matrix norm %14.12e\n", n, (double)rnorm, (double)lsqr->arnorm, (double)lsqr->anorm);
383: }
384: PetscViewerASCIISubtractTab(viewer, tablevel);
385: PetscViewerPopFormat(viewer);
386: return 0;
387: }
389: /*@C
390: KSPLSQRMonitorResidual - Prints the residual norm, as well as the normal equation residual norm, at each iteration of an iterative solver for the `KSPLSQR` solver
392: Collective
394: Input Parameters:
395: + ksp - iterative context
396: . n - iteration number
397: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
398: - vf - The viewer context
400: Options Database Key:
401: . -ksp_lsqr_monitor - Activates `KSPLSQRMonitorResidual()`
403: Level: intermediate
405: .seealso: [](chapter_ksp), `KSPLSQR`, `KSPMonitorSet()`, `KSPMonitorResidual()`, `KSPMonitorTrueResidualMaxNorm()`, `KSPLSQRMonitorResidualDrawLG()`
406: @*/
407: PetscErrorCode KSPLSQRMonitorResidual(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
408: {
412: PetscTryMethod(ksp, "KSPLSQRMonitorResidual_C", (KSP, PetscInt, PetscReal, PetscViewerAndFormat *), (ksp, n, rnorm, vf));
413: return 0;
414: }
416: PetscErrorCode KSPLSQRMonitorResidualDrawLG_LSQR(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
417: {
418: KSP_LSQR *lsqr = (KSP_LSQR *)ksp->data;
419: PetscViewer viewer = vf->viewer;
420: PetscViewerFormat format = vf->format;
421: PetscDrawLG lg = vf->lg;
422: KSPConvergedReason reason;
423: PetscReal x[2], y[2];
425: PetscViewerPushFormat(viewer, format);
426: if (!n) PetscDrawLGReset(lg);
427: x[0] = (PetscReal)n;
428: if (rnorm > 0.0) y[0] = PetscLog10Real(rnorm);
429: else y[0] = -15.0;
430: x[1] = (PetscReal)n;
431: if (lsqr->arnorm > 0.0) y[1] = PetscLog10Real(lsqr->arnorm);
432: else y[1] = -15.0;
433: PetscDrawLGAddPoint(lg, x, y);
434: KSPGetConvergedReason(ksp, &reason);
435: if (n <= 20 || !(n % 5) || reason) {
436: PetscDrawLGDraw(lg);
437: PetscDrawLGSave(lg);
438: }
439: PetscViewerPopFormat(viewer);
440: return 0;
441: }
443: /*@C
444: KSPLSQRMonitorResidualDrawLG - Plots the true residual norm at each iteration of an iterative solver for the `KSPLSQR` solver
446: Collective
448: Input Parameters:
449: + ksp - iterative context
450: . n - iteration number
451: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
452: - vf - The viewer context
454: Options Database Key:
455: . -ksp_lsqr_monitor draw::draw_lg - Activates `KSPMonitorTrueResidualDrawLG()`
457: Level: intermediate
459: .seealso: [](chapter_ksp), `KSPLSQR`, `KSPMonitorSet()`, `KSPMonitorTrueResidual()`, `KSPLSQRMonitorResidual()`, `KSPLSQRMonitorResidualDrawLGCreate()`
460: @*/
461: PetscErrorCode KSPLSQRMonitorResidualDrawLG(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
462: {
467: PetscTryMethod(ksp, "KSPLSQRMonitorResidualDrawLG_C", (KSP, PetscInt, PetscReal, PetscViewerAndFormat *), (ksp, n, rnorm, vf));
468: return 0;
469: }
471: /*@C
472: KSPLSQRMonitorResidualDrawLGCreate - Creates the plotter for the `KSPLSQR` residual and normal equation residual norm
474: Collective
476: Input Parameters:
477: + viewer - The PetscViewer
478: . format - The viewer format
479: - ctx - An optional user context
481: Output Parameter:
482: . vf - The viewer context
484: Level: intermediate
486: .seealso: [](chapter_ksp), `KSPLSQR`, `KSPMonitorSet()`, `KSPLSQRMonitorResidual()`, `KSPLSQRMonitorResidualDrawLG()`
487: @*/
488: PetscErrorCode KSPLSQRMonitorResidualDrawLGCreate(PetscViewer viewer, PetscViewerFormat format, void *ctx, PetscViewerAndFormat **vf)
489: {
490: const char *names[] = {"residual", "normal eqn residual"};
492: PetscViewerAndFormatCreate(viewer, format, vf);
493: (*vf)->data = ctx;
494: KSPMonitorLGCreate(PetscObjectComm((PetscObject)viewer), NULL, NULL, "Log Residual Norm", 2, names, PETSC_DECIDE, PETSC_DECIDE, 400, 300, &(*vf)->lg);
495: return 0;
496: }
498: PetscErrorCode KSPSetFromOptions_LSQR(KSP ksp, PetscOptionItems *PetscOptionsObject)
499: {
500: KSP_LSQR *lsqr = (KSP_LSQR *)ksp->data;
502: PetscOptionsHeadBegin(PetscOptionsObject, "KSP LSQR Options");
503: PetscOptionsBool("-ksp_lsqr_compute_standard_error", "Set Standard Error Estimates of Solution", "KSPLSQRSetComputeStandardErrorVec", lsqr->se_flg, &lsqr->se_flg, NULL);
504: PetscOptionsBool("-ksp_lsqr_exact_mat_norm", "Compute exact matrix norm instead of iteratively refined estimate", "KSPLSQRSetExactMatNorm", lsqr->exact_norm, &lsqr->exact_norm, NULL);
505: KSPMonitorSetFromOptions(ksp, "-ksp_lsqr_monitor", "lsqr_residual", NULL);
506: PetscOptionsHeadEnd();
507: return 0;
508: }
510: PetscErrorCode KSPView_LSQR(KSP ksp, PetscViewer viewer)
511: {
512: KSP_LSQR *lsqr = (KSP_LSQR *)ksp->data;
513: PetscBool iascii;
515: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
516: if (iascii) {
517: if (lsqr->se) {
518: PetscReal rnorm;
519: VecNorm(lsqr->se, NORM_2, &rnorm);
520: PetscViewerASCIIPrintf(viewer, " norm of standard error %g, iterations %" PetscInt_FMT "\n", (double)rnorm, ksp->its);
521: } else {
522: PetscViewerASCIIPrintf(viewer, " standard error not computed\n");
523: }
524: if (lsqr->exact_norm) {
525: PetscViewerASCIIPrintf(viewer, " using exact matrix norm\n");
526: } else {
527: PetscViewerASCIIPrintf(viewer, " using inexact matrix norm\n");
528: }
529: }
530: return 0;
531: }
533: /*@C
534: KSPLSQRConvergedDefault - Determines convergence of the `KSPLSQR` Krylov method.
536: Collective
538: Input Parameters:
539: + ksp - iterative context
540: . n - iteration number
541: . rnorm - 2-norm residual value (may be estimated)
542: - ctx - convergence context which must be created by `KSPConvergedDefaultCreate()`
544: reason is set to:
545: + positive - if the iteration has converged;
546: . negative - if residual norm exceeds divergence threshold;
547: - 0 - otherwise.
549: Notes:
550: `KSPConvergedDefault()` is called first to check for convergence in A*x=b.
551: If that does not determine convergence then checks convergence for the least squares problem, i.e. in min{|b-A*x|}.
552: Possible convergence for the least squares problem (which is based on the residual of the normal equations) are `KSP_CONVERGED_RTOL_NORMAL` norm
553: and `KSP_CONVERGED_ATOL_NORMAL`.
555: `KSP_CONVERGED_RTOL_NORMAL` is returned if ||A'*r|| < rtol * ||A|| * ||r||.
556: Matrix norm ||A|| is iteratively refined estimate, see `KSPLSQRGetNorms()`.
557: This criterion is now largely compatible with that in MATLAB `lsqr()`.
559: Level: intermediate
561: .seealso: [](chapter_ksp), `KSPLSQR`, `KSPSetConvergenceTest()`, `KSPSetTolerances()`, `KSPConvergedSkip()`, `KSPConvergedReason`, `KSPGetConvergedReason()`,
562: `KSPConvergedDefaultSetUIRNorm()`, `KSPConvergedDefaultSetUMIRNorm()`, `KSPConvergedDefaultCreate()`, `KSPConvergedDefaultDestroy()`, `KSPConvergedDefault()`, `KSPLSQRGetNorms()`, `KSPLSQRSetExactMatNorm()`
563: @*/
564: PetscErrorCode KSPLSQRConvergedDefault(KSP ksp, PetscInt n, PetscReal rnorm, KSPConvergedReason *reason, void *ctx)
565: {
566: KSP_LSQR *lsqr = (KSP_LSQR *)ksp->data;
568: /* check for convergence in A*x=b */
569: KSPConvergedDefault(ksp, n, rnorm, reason, ctx);
570: if (!n || *reason) return 0;
572: /* check for convergence in min{|b-A*x|} */
573: if (lsqr->arnorm < ksp->abstol) {
574: PetscInfo(ksp, "LSQR solver has converged. Normal equation residual %14.12e is less than absolute tolerance %14.12e at iteration %" PetscInt_FMT "\n", (double)lsqr->arnorm, (double)ksp->abstol, n);
575: *reason = KSP_CONVERGED_ATOL_NORMAL;
576: } else if (lsqr->arnorm < ksp->rtol * lsqr->anorm * rnorm) {
577: PetscCall(PetscInfo(ksp, "LSQR solver has converged. Normal equation residual %14.12e is less than rel. tol. %14.12e times %s Frobenius norm of matrix %14.12e times residual %14.12e at iteration %" PetscInt_FMT "\n", (double)lsqr->arnorm,
578: (double)ksp->rtol, lsqr->exact_norm ? "exact" : "approx.", (double)lsqr->anorm, (double)rnorm, n));
579: *reason = KSP_CONVERGED_RTOL_NORMAL;
580: }
581: return 0;
582: }
584: /*MC
585: KSPLSQR - Implements LSQR
587: Options Database Keys:
588: + -ksp_lsqr_set_standard_error - set standard error estimates of solution, see `KSPLSQRSetComputeStandardErrorVec()` and `KSPLSQRGetStandardErrorVec()`
589: . -ksp_lsqr_exact_mat_norm - compute exact matrix norm instead of iteratively refined estimate, see `KSPLSQRSetExactMatNorm()`
590: - -ksp_lsqr_monitor - monitor residual norm, norm of residual of normal equations A'*A x = A' b, and estimate of matrix norm ||A||
592: Level: beginner
594: Notes:
595: Supports non-square (rectangular) matrices.
597: This variant, when applied with no preconditioning is identical to the original algorithm in exact arithmetic; however, in practice, with no preconditioning
598: due to inexact arithmetic, it can converge differently. Hence when no preconditioner is used (`PCType` `PCNONE`) it automatically reverts to the original algorithm.
600: With the PETSc built-in preconditioners, such as `PCICC`, one should call `KSPSetOperators`(ksp,A,A'*A)) since the preconditioner needs to work
601: for the normal equations A'*A.
603: Supports only left preconditioning.
605: For least squares problems with nonzero residual A*x - b, there are additional convergence tests for the residual of the normal equations, A'*(b - Ax), see `KSPLSQRConvergedDefault()`.
607: In exact arithmetic the LSQR method (with no preconditioning) is identical to the `KSPCG` algorithm applied to the normal equations.
608: The preconditioned variant was implemented by Bas van't Hof and is essentially a left preconditioning for the Normal Equations.
609: It appears the implementation with preconditioning tracks the true norm of the residual and uses that in the convergence test.
611: Developer Note:
612: How is this related to the `KSPCGNE` implementation? One difference is that `KSPCGNE` applies
613: the preconditioner transpose times the preconditioner, so one does not need to pass A'*A as the third argument to `KSPSetOperators()`.
615: Reference:
616: . * - The original unpreconditioned algorithm can be found in Paige and Saunders, ACM Transactions on Mathematical Software, Vol 8, 1982.
618: .seealso: [](chapter_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPSolve()`, `KSPLSQRConvergedDefault()`, `KSPLSQRSetComputeStandardErrorVec()`, `KSPLSQRGetStandardErrorVec()`, `KSPLSQRSetExactMatNorm()`, `KSPLSQRMonitorResidualDrawLGCreate()`, `KSPLSQRMonitorResidualDrawLG()`, `KSPLSQRMonitorResidual()`
619: M*/
620: PETSC_EXTERN PetscErrorCode KSPCreate_LSQR(KSP ksp)
621: {
622: KSP_LSQR *lsqr;
623: void *ctx;
625: PetscNew(&lsqr);
626: lsqr->se = NULL;
627: lsqr->se_flg = PETSC_FALSE;
628: lsqr->exact_norm = PETSC_FALSE;
629: lsqr->anorm = -1.0;
630: lsqr->arnorm = -1.0;
631: ksp->data = (void *)lsqr;
632: KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 3);
634: ksp->ops->setup = KSPSetUp_LSQR;
635: ksp->ops->solve = KSPSolve_LSQR;
636: ksp->ops->destroy = KSPDestroy_LSQR;
637: ksp->ops->setfromoptions = KSPSetFromOptions_LSQR;
638: ksp->ops->view = KSPView_LSQR;
640: /* Backup current convergence test; remove destroy routine from KSP to prevent destroying the convergence context in KSPSetConvergenceTest() */
641: KSPGetAndClearConvergenceTest(ksp, &lsqr->converged, &lsqr->cnvP, &lsqr->convergeddestroy);
642: /* Override current convergence test */
643: KSPConvergedDefaultCreate(&ctx);
644: KSPSetConvergenceTest(ksp, KSPLSQRConvergedDefault, ctx, KSPConvergedDefaultDestroy);
645: PetscObjectComposeFunction((PetscObject)ksp, "KSPLSQRMonitorResidual_C", KSPLSQRMonitorResidual_LSQR);
646: PetscObjectComposeFunction((PetscObject)ksp, "KSPLSQRMonitorResidualDrawLG_C", KSPLSQRMonitorResidualDrawLG_LSQR);
647: return 0;
648: }