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: PetscFunctionBegin;
36: PetscCall(VecGetLocalSize(v, &n));
37: PetscCall(VecGetArray(v, &x));
38: for (i = 0; i < n; i++) x[i] *= PetscConj(x[i]);
39: PetscCall(VecRestoreArray(v, &x));
40: PetscFunctionReturn(PETSC_SUCCESS);
41: }
43: static PetscErrorCode KSPSetUp_LSQR(KSP ksp)
44: {
45: KSP_LSQR *lsqr = (KSP_LSQR *)ksp->data;
46: PetscBool nopreconditioner;
48: PetscFunctionBegin;
49: PetscCall(PetscObjectTypeCompare((PetscObject)ksp->pc, PCNONE, &nopreconditioner));
51: if (lsqr->vwork_m) PetscCall(VecDestroyVecs(lsqr->nwork_m, &lsqr->vwork_m));
53: if (lsqr->vwork_n) PetscCall(VecDestroyVecs(lsqr->nwork_n, &lsqr->vwork_n));
55: lsqr->nwork_m = 2;
56: if (nopreconditioner) lsqr->nwork_n = 4;
57: else lsqr->nwork_n = 5;
58: PetscCall(KSPCreateVecs(ksp, lsqr->nwork_n, &lsqr->vwork_n, lsqr->nwork_m, &lsqr->vwork_m));
60: if (lsqr->se_flg && !lsqr->se) {
61: PetscCall(VecDuplicate(lsqr->vwork_n[0], &lsqr->se));
62: PetscCall(VecSet(lsqr->se, PETSC_INFINITY));
63: } else if (!lsqr->se_flg) {
64: PetscCall(VecDestroy(&lsqr->se));
65: }
66: PetscFunctionReturn(PETSC_SUCCESS);
67: }
69: static PetscErrorCode KSPSolve_LSQR(KSP ksp)
70: {
71: PetscInt i, size1, size2;
72: PetscScalar rho, rhobar, phi, phibar, theta, c, s, tmp, tau;
73: PetscReal beta, alpha, rnorm;
74: Vec X, B, V, V1, U, U1, TMP, W, W2, Z = NULL;
75: Mat Amat, Pmat;
76: KSP_LSQR *lsqr = (KSP_LSQR *)ksp->data;
77: PetscBool diagonalscale, nopreconditioner;
79: PetscFunctionBegin;
80: PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
81: PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
83: PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));
84: PetscCall(PetscObjectTypeCompare((PetscObject)ksp->pc, PCNONE, &nopreconditioner));
86: /* vectors of length m, where system size is mxn */
87: B = ksp->vec_rhs;
88: U = lsqr->vwork_m[0];
89: U1 = lsqr->vwork_m[1];
91: /* vectors of length n */
92: X = ksp->vec_sol;
93: W = lsqr->vwork_n[0];
94: V = lsqr->vwork_n[1];
95: V1 = lsqr->vwork_n[2];
96: W2 = lsqr->vwork_n[3];
97: if (!nopreconditioner) Z = lsqr->vwork_n[4];
99: /* standard error vector */
100: if (lsqr->se) PetscCall(VecSet(lsqr->se, 0.0));
102: /* Compute initial residual, temporarily use work vector u */
103: if (!ksp->guess_zero) {
104: PetscCall(KSP_MatMult(ksp, Amat, X, U)); /* u <- b - Ax */
105: PetscCall(VecAYPX(U, -1.0, B));
106: } else {
107: PetscCall(VecCopy(B, U)); /* u <- b (x is 0) */
108: }
110: /* Test for nothing to do */
111: PetscCall(VecNorm(U, NORM_2, &rnorm));
112: KSPCheckNorm(ksp, rnorm);
113: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
114: ksp->its = 0;
115: ksp->rnorm = rnorm;
116: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
117: PetscCall(KSPLogResidualHistory(ksp, rnorm));
118: PetscCall(KSPMonitor(ksp, 0, rnorm));
119: PetscCall((*ksp->converged)(ksp, 0, rnorm, &ksp->reason, ksp->cnvP));
120: if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);
122: beta = rnorm;
123: PetscCall(VecScale(U, 1.0 / beta));
124: PetscCall(KSP_MatMultHermitianTranspose(ksp, Amat, U, V));
125: if (nopreconditioner) {
126: PetscCall(VecNorm(V, NORM_2, &alpha));
127: KSPCheckNorm(ksp, rnorm);
128: } else {
129: /* this is an application of the preconditioner for the normal equations; not the operator, see the manual page */
130: PetscCall(PCApply(ksp->pc, V, Z));
131: PetscCall(VecDotRealPart(V, Z, &alpha));
132: if (alpha <= 0.0) {
133: ksp->reason = KSP_DIVERGED_BREAKDOWN;
134: PetscFunctionReturn(PETSC_SUCCESS);
135: }
136: alpha = PetscSqrtReal(alpha);
137: PetscCall(VecScale(Z, 1.0 / alpha));
138: }
139: PetscCall(VecScale(V, 1.0 / alpha));
141: if (nopreconditioner) {
142: PetscCall(VecCopy(V, W));
143: } else {
144: PetscCall(VecCopy(Z, W));
145: }
147: if (lsqr->exact_norm) {
148: PetscCall(MatNorm(Amat, NORM_FROBENIUS, &lsqr->anorm));
149: } else lsqr->anorm = 0.0;
151: lsqr->arnorm = alpha * beta;
152: phibar = beta;
153: rhobar = alpha;
154: i = 0;
155: do {
156: if (nopreconditioner) {
157: PetscCall(KSP_MatMult(ksp, Amat, V, U1));
158: } else {
159: PetscCall(KSP_MatMult(ksp, Amat, Z, U1));
160: }
161: PetscCall(VecAXPY(U1, -alpha, U));
162: PetscCall(VecNorm(U1, NORM_2, &beta));
163: KSPCheckNorm(ksp, beta);
164: if (beta > 0.0) {
165: PetscCall(VecScale(U1, 1.0 / beta)); /* beta*U1 = Amat*V - alpha*U */
166: if (!lsqr->exact_norm) lsqr->anorm = PetscSqrtReal(PetscSqr(lsqr->anorm) + PetscSqr(alpha) + PetscSqr(beta));
167: }
169: PetscCall(KSP_MatMultHermitianTranspose(ksp, Amat, U1, V1));
170: PetscCall(VecAXPY(V1, -beta, V));
171: if (nopreconditioner) {
172: PetscCall(VecNorm(V1, NORM_2, &alpha));
173: KSPCheckNorm(ksp, alpha);
174: } else {
175: PetscCall(PCApply(ksp->pc, V1, Z));
176: PetscCall(VecDotRealPart(V1, Z, &alpha));
177: if (alpha <= 0.0) {
178: ksp->reason = KSP_DIVERGED_BREAKDOWN;
179: break;
180: }
181: alpha = PetscSqrtReal(alpha);
182: PetscCall(VecScale(Z, 1.0 / alpha));
183: }
184: PetscCall(VecScale(V1, 1.0 / alpha)); /* alpha*V1 = Amat^T*U1 - beta*V */
185: rho = PetscSqrtScalar(rhobar * rhobar + beta * beta);
186: c = rhobar / rho;
187: s = beta / rho;
188: theta = s * alpha;
189: rhobar = -c * alpha;
190: phi = c * phibar;
191: phibar = s * phibar;
192: tau = s * phi;
194: PetscCall(VecAXPY(X, phi / rho, W)); /* x <- x + (phi/rho) w */
196: if (lsqr->se) {
197: PetscCall(VecCopy(W, W2));
198: PetscCall(VecSquare(W2));
199: PetscCall(VecScale(W2, 1.0 / (rho * rho)));
200: PetscCall(VecAXPY(lsqr->se, 1.0, W2)); /* lsqr->se <- lsqr->se + (w^2/rho^2) */
201: }
202: if (nopreconditioner) {
203: PetscCall(VecAYPX(W, -theta / rho, V1)); /* w <- v - (theta/rho) w */
204: } else {
205: PetscCall(VecAYPX(W, -theta / rho, Z)); /* w <- z - (theta/rho) w */
206: }
208: lsqr->arnorm = alpha * PetscAbsScalar(tau);
209: rnorm = PetscRealPart(phibar);
211: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
212: ksp->its++;
213: ksp->rnorm = rnorm;
214: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
215: PetscCall(KSPLogResidualHistory(ksp, rnorm));
216: PetscCall(KSPMonitor(ksp, i + 1, rnorm));
217: PetscCall((*ksp->converged)(ksp, i + 1, rnorm, &ksp->reason, ksp->cnvP));
218: if (ksp->reason) break;
219: SWAP(U1, U, TMP);
220: SWAP(V1, V, TMP);
222: i++;
223: } while (i < ksp->max_it);
224: if (i >= ksp->max_it && !ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
226: /* Finish off the standard error estimates */
227: if (lsqr->se) {
228: tmp = 1.0;
229: PetscCall(MatGetSize(Amat, &size1, &size2));
230: if (size1 > size2) tmp = size1 - size2;
231: tmp = rnorm / PetscSqrtScalar(tmp);
232: PetscCall(VecSqrtAbs(lsqr->se));
233: PetscCall(VecScale(lsqr->se, tmp));
234: }
235: PetscFunctionReturn(PETSC_SUCCESS);
236: }
238: PetscErrorCode KSPDestroy_LSQR(KSP ksp)
239: {
240: KSP_LSQR *lsqr = (KSP_LSQR *)ksp->data;
242: PetscFunctionBegin;
243: /* Free work vectors */
244: if (lsqr->vwork_n) PetscCall(VecDestroyVecs(lsqr->nwork_n, &lsqr->vwork_n));
245: if (lsqr->vwork_m) PetscCall(VecDestroyVecs(lsqr->nwork_m, &lsqr->vwork_m));
246: PetscCall(VecDestroy(&lsqr->se));
247: /* Revert convergence test */
248: PetscCall(KSPSetConvergenceTest(ksp, lsqr->converged, lsqr->cnvP, lsqr->convergeddestroy));
249: /* Free the KSP_LSQR context */
250: PetscCall(PetscFree(ksp->data));
251: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPLSQRMonitorResidual_C", NULL));
252: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPLSQRMonitorResidualDrawLG_C", NULL));
253: PetscFunctionReturn(PETSC_SUCCESS);
254: }
256: /*@
257: KSPLSQRSetComputeStandardErrorVec - Compute a vector of standard error estimates during `KSPSolve()` for `KSPLSQR`.
259: Logically Collective
261: Input Parameters:
262: + ksp - iterative context
263: - flg - compute the vector of standard estimates or not
265: Level: intermediate
267: Developer Note:
268: Vaclav: I'm not sure whether this vector is useful for anything.
270: .seealso: [](ch_ksp), `KSPSolve()`, `KSPLSQR`, `KSPLSQRGetStandardErrorVec()`
271: @*/
272: PetscErrorCode KSPLSQRSetComputeStandardErrorVec(KSP ksp, PetscBool flg)
273: {
274: KSP_LSQR *lsqr = (KSP_LSQR *)ksp->data;
276: PetscFunctionBegin;
277: lsqr->se_flg = flg;
278: PetscFunctionReturn(PETSC_SUCCESS);
279: }
281: /*@
282: KSPLSQRSetExactMatNorm - Compute exact matrix norm instead of iteratively refined estimate.
284: Not Collective
286: Input Parameters:
287: + ksp - iterative context
288: - flg - compute exact matrix norm or not
290: Level: intermediate
292: Notes:
293: By default, flg = `PETSC_FALSE`. This is usually preferred to avoid possibly expensive computation of the norm.
294: For flg = `PETSC_TRUE`, we call `MatNorm`(Amat,`NORM_FROBENIUS`,&lsqr->anorm) which will work only for some types of explicitly assembled matrices.
295: This can affect convergence rate as `KSPLSQRConvergedDefault()` assumes different value of ||A|| used in normal equation stopping criterion.
297: .seealso: [](ch_ksp), `KSPSolve()`, `KSPLSQR`, `KSPLSQRGetNorms()`, `KSPLSQRConvergedDefault()`
298: @*/
299: PetscErrorCode KSPLSQRSetExactMatNorm(KSP ksp, PetscBool flg)
300: {
301: KSP_LSQR *lsqr = (KSP_LSQR *)ksp->data;
303: PetscFunctionBegin;
304: lsqr->exact_norm = flg;
305: PetscFunctionReturn(PETSC_SUCCESS);
306: }
308: /*@
309: KSPLSQRGetStandardErrorVec - Get vector of standard error estimates.
310: Only available if -ksp_lsqr_set_standard_error was set to true
311: or `KSPLSQRSetComputeStandardErrorVec`(ksp, `PETSC_TRUE`) was called.
312: Otherwise returns NULL.
314: Not Collective
316: Input Parameter:
317: . ksp - iterative context
319: Output Parameter:
320: . se - vector of standard estimates
322: Level: intermediate
324: Developer Note:
325: Vaclav: I'm not sure whether this vector is useful for anything.
327: .seealso: [](ch_ksp), `KSPSolve()`, `KSPLSQR`, `KSPLSQRSetComputeStandardErrorVec()`
328: @*/
329: PetscErrorCode KSPLSQRGetStandardErrorVec(KSP ksp, Vec *se)
330: {
331: KSP_LSQR *lsqr = (KSP_LSQR *)ksp->data;
333: PetscFunctionBegin;
334: *se = lsqr->se;
335: PetscFunctionReturn(PETSC_SUCCESS);
336: }
338: /*@
339: KSPLSQRGetNorms - Get the norm estimates that `KSPLSQR` computes internally during `KSPSolve()`.
341: Not Collective
343: Input Parameter:
344: . ksp - iterative context
346: Output Parameters:
347: + arnorm - good estimate of norm((A*inv(Pmat))'*r), where r = A*x - b, used in specific stopping criterion
348: - anorm - poor estimate of norm(A*inv(Pmat),'fro') used in specific stopping criterion
350: Notes:
351: Output parameters are meaningful only after `KSPSolve()`.
353: 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.
355: If -ksp_lsqr_exact_mat_norm is set or `KSPLSQRSetExactMatNorm`(ksp, `PETSC_TRUE`) called, then anorm is the exact Frobenius norm.
357: Level: intermediate
359: .seealso: [](ch_ksp), `KSPSolve()`, `KSPLSQR`, `KSPLSQRSetExactMatNorm()`
360: @*/
361: PetscErrorCode KSPLSQRGetNorms(KSP ksp, PetscReal *arnorm, PetscReal *anorm)
362: {
363: KSP_LSQR *lsqr = (KSP_LSQR *)ksp->data;
365: PetscFunctionBegin;
366: if (arnorm) *arnorm = lsqr->arnorm;
367: if (anorm) *anorm = lsqr->anorm;
368: PetscFunctionReturn(PETSC_SUCCESS);
369: }
371: PetscErrorCode KSPLSQRMonitorResidual_LSQR(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
372: {
373: KSP_LSQR *lsqr = (KSP_LSQR *)ksp->data;
374: PetscViewer viewer = vf->viewer;
375: PetscViewerFormat format = vf->format;
376: char normtype[256];
377: PetscInt tablevel;
378: const char *prefix;
380: PetscFunctionBegin;
381: PetscCall(PetscObjectGetTabLevel((PetscObject)ksp, &tablevel));
382: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)ksp, &prefix));
383: PetscCall(PetscStrncpy(normtype, KSPNormTypes[ksp->normtype], sizeof(normtype)));
384: PetscCall(PetscStrtolower(normtype));
385: PetscCall(PetscViewerPushFormat(viewer, format));
386: PetscCall(PetscViewerASCIIAddTab(viewer, tablevel));
387: if (n == 0 && prefix) PetscCall(PetscViewerASCIIPrintf(viewer, " Residual norm, norm of normal equations, and matrix norm for %s solve.\n", prefix));
388: if (!n) {
389: PetscCall(PetscViewerASCIIPrintf(viewer, "%3" PetscInt_FMT " KSP resid norm %14.12e\n", n, (double)rnorm));
390: } else {
391: PetscCall(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));
392: }
393: PetscCall(PetscViewerASCIISubtractTab(viewer, tablevel));
394: PetscCall(PetscViewerPopFormat(viewer));
395: PetscFunctionReturn(PETSC_SUCCESS);
396: }
398: /*@C
399: KSPLSQRMonitorResidual - Prints the residual norm, as well as the normal equation residual norm, at each iteration of an iterative solver for the `KSPLSQR` solver
401: Collective
403: Input Parameters:
404: + ksp - iterative context
405: . n - iteration number
406: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
407: - vf - The viewer context
409: Options Database Key:
410: . -ksp_lsqr_monitor - Activates `KSPLSQRMonitorResidual()`
412: Level: intermediate
414: .seealso: [](ch_ksp), `KSPLSQR`, `KSPMonitorSet()`, `KSPMonitorResidual()`, `KSPMonitorTrueResidualMaxNorm()`, `KSPLSQRMonitorResidualDrawLG()`
415: @*/
416: PetscErrorCode KSPLSQRMonitorResidual(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
417: {
418: PetscFunctionBegin;
422: PetscTryMethod(ksp, "KSPLSQRMonitorResidual_C", (KSP, PetscInt, PetscReal, PetscViewerAndFormat *), (ksp, n, rnorm, vf));
423: PetscFunctionReturn(PETSC_SUCCESS);
424: }
426: PetscErrorCode KSPLSQRMonitorResidualDrawLG_LSQR(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
427: {
428: KSP_LSQR *lsqr = (KSP_LSQR *)ksp->data;
429: PetscViewer viewer = vf->viewer;
430: PetscViewerFormat format = vf->format;
431: PetscDrawLG lg = vf->lg;
432: KSPConvergedReason reason;
433: PetscReal x[2], y[2];
435: PetscFunctionBegin;
436: PetscCall(PetscViewerPushFormat(viewer, format));
437: if (!n) PetscCall(PetscDrawLGReset(lg));
438: x[0] = (PetscReal)n;
439: if (rnorm > 0.0) y[0] = PetscLog10Real(rnorm);
440: else y[0] = -15.0;
441: x[1] = (PetscReal)n;
442: if (lsqr->arnorm > 0.0) y[1] = PetscLog10Real(lsqr->arnorm);
443: else y[1] = -15.0;
444: PetscCall(PetscDrawLGAddPoint(lg, x, y));
445: PetscCall(KSPGetConvergedReason(ksp, &reason));
446: if (n <= 20 || !(n % 5) || reason) {
447: PetscCall(PetscDrawLGDraw(lg));
448: PetscCall(PetscDrawLGSave(lg));
449: }
450: PetscCall(PetscViewerPopFormat(viewer));
451: PetscFunctionReturn(PETSC_SUCCESS);
452: }
454: /*@C
455: KSPLSQRMonitorResidualDrawLG - Plots the true residual norm at each iteration of an iterative solver for the `KSPLSQR` solver
457: Collective
459: Input Parameters:
460: + ksp - iterative context
461: . n - iteration number
462: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
463: - vf - The viewer context
465: Options Database Key:
466: . -ksp_lsqr_monitor draw::draw_lg - Activates `KSPMonitorTrueResidualDrawLG()`
468: Level: intermediate
470: .seealso: [](ch_ksp), `KSPLSQR`, `KSPMonitorSet()`, `KSPMonitorTrueResidual()`, `KSPLSQRMonitorResidual()`, `KSPLSQRMonitorResidualDrawLGCreate()`
471: @*/
472: PetscErrorCode KSPLSQRMonitorResidualDrawLG(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
473: {
474: PetscFunctionBegin;
479: PetscTryMethod(ksp, "KSPLSQRMonitorResidualDrawLG_C", (KSP, PetscInt, PetscReal, PetscViewerAndFormat *), (ksp, n, rnorm, vf));
480: PetscFunctionReturn(PETSC_SUCCESS);
481: }
483: /*@C
484: KSPLSQRMonitorResidualDrawLGCreate - Creates the plotter for the `KSPLSQR` residual and normal equation residual norm
486: Collective
488: Input Parameters:
489: + viewer - The PetscViewer
490: . format - The viewer format
491: - ctx - An optional user context
493: Output Parameter:
494: . vf - The viewer context
496: Level: intermediate
498: .seealso: [](ch_ksp), `KSPLSQR`, `KSPMonitorSet()`, `KSPLSQRMonitorResidual()`, `KSPLSQRMonitorResidualDrawLG()`
499: @*/
500: PetscErrorCode KSPLSQRMonitorResidualDrawLGCreate(PetscViewer viewer, PetscViewerFormat format, void *ctx, PetscViewerAndFormat **vf)
501: {
502: const char *names[] = {"residual", "normal eqn residual"};
504: PetscFunctionBegin;
505: PetscCall(PetscViewerAndFormatCreate(viewer, format, vf));
506: (*vf)->data = ctx;
507: PetscCall(KSPMonitorLGCreate(PetscObjectComm((PetscObject)viewer), NULL, NULL, "Log Residual Norm", 2, names, PETSC_DECIDE, PETSC_DECIDE, 400, 300, &(*vf)->lg));
508: PetscFunctionReturn(PETSC_SUCCESS);
509: }
511: PetscErrorCode KSPSetFromOptions_LSQR(KSP ksp, PetscOptionItems *PetscOptionsObject)
512: {
513: KSP_LSQR *lsqr = (KSP_LSQR *)ksp->data;
515: PetscFunctionBegin;
516: PetscOptionsHeadBegin(PetscOptionsObject, "KSP LSQR Options");
517: PetscCall(PetscOptionsBool("-ksp_lsqr_compute_standard_error", "Set Standard Error Estimates of Solution", "KSPLSQRSetComputeStandardErrorVec", lsqr->se_flg, &lsqr->se_flg, NULL));
518: PetscCall(PetscOptionsBool("-ksp_lsqr_exact_mat_norm", "Compute exact matrix norm instead of iteratively refined estimate", "KSPLSQRSetExactMatNorm", lsqr->exact_norm, &lsqr->exact_norm, NULL));
519: PetscCall(KSPMonitorSetFromOptions(ksp, "-ksp_lsqr_monitor", "lsqr_residual", NULL));
520: PetscOptionsHeadEnd();
521: PetscFunctionReturn(PETSC_SUCCESS);
522: }
524: PetscErrorCode KSPView_LSQR(KSP ksp, PetscViewer viewer)
525: {
526: KSP_LSQR *lsqr = (KSP_LSQR *)ksp->data;
527: PetscBool iascii;
529: PetscFunctionBegin;
530: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
531: if (iascii) {
532: if (lsqr->se) {
533: PetscReal rnorm;
534: PetscCall(VecNorm(lsqr->se, NORM_2, &rnorm));
535: PetscCall(PetscViewerASCIIPrintf(viewer, " norm of standard error %g, iterations %" PetscInt_FMT "\n", (double)rnorm, ksp->its));
536: } else {
537: PetscCall(PetscViewerASCIIPrintf(viewer, " standard error not computed\n"));
538: }
539: if (lsqr->exact_norm) {
540: PetscCall(PetscViewerASCIIPrintf(viewer, " using exact matrix norm\n"));
541: } else {
542: PetscCall(PetscViewerASCIIPrintf(viewer, " using inexact matrix norm\n"));
543: }
544: }
545: PetscFunctionReturn(PETSC_SUCCESS);
546: }
548: /*@C
549: KSPLSQRConvergedDefault - Determines convergence of the `KSPLSQR` Krylov method.
551: Collective
553: Input Parameters:
554: + ksp - iterative context
555: . n - iteration number
556: . rnorm - 2-norm residual value (may be estimated)
557: - ctx - convergence context which must be created by `KSPConvergedDefaultCreate()`
559: reason is set to:
560: + positive - if the iteration has converged;
561: . negative - if residual norm exceeds divergence threshold;
562: - 0 - otherwise.
564: Notes:
565: `KSPConvergedDefault()` is called first to check for convergence in A*x=b.
566: If that does not determine convergence then checks convergence for the least squares problem, i.e. in min{|b-A*x|}.
567: Possible convergence for the least squares problem (which is based on the residual of the normal equations) are `KSP_CONVERGED_RTOL_NORMAL` norm
568: and `KSP_CONVERGED_ATOL_NORMAL`.
570: `KSP_CONVERGED_RTOL_NORMAL` is returned if ||A'*r|| < rtol * ||A|| * ||r||.
571: Matrix norm ||A|| is iteratively refined estimate, see `KSPLSQRGetNorms()`.
572: This criterion is now largely compatible with that in MATLAB `lsqr()`.
574: Level: intermediate
576: .seealso: [](ch_ksp), `KSPLSQR`, `KSPSetConvergenceTest()`, `KSPSetTolerances()`, `KSPConvergedSkip()`, `KSPConvergedReason`, `KSPGetConvergedReason()`,
577: `KSPConvergedDefaultSetUIRNorm()`, `KSPConvergedDefaultSetUMIRNorm()`, `KSPConvergedDefaultCreate()`, `KSPConvergedDefaultDestroy()`, `KSPConvergedDefault()`, `KSPLSQRGetNorms()`, `KSPLSQRSetExactMatNorm()`
578: @*/
579: PetscErrorCode KSPLSQRConvergedDefault(KSP ksp, PetscInt n, PetscReal rnorm, KSPConvergedReason *reason, void *ctx)
580: {
581: KSP_LSQR *lsqr = (KSP_LSQR *)ksp->data;
583: PetscFunctionBegin;
584: /* check for convergence in A*x=b */
585: PetscCall(KSPConvergedDefault(ksp, n, rnorm, reason, ctx));
586: if (!n || *reason) PetscFunctionReturn(PETSC_SUCCESS);
588: /* check for convergence in min{|b-A*x|} */
589: if (lsqr->arnorm < ksp->abstol) {
590: PetscCall(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));
591: *reason = KSP_CONVERGED_ATOL_NORMAL;
592: } else if (lsqr->arnorm < ksp->rtol * lsqr->anorm * rnorm) {
593: 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,
594: (double)ksp->rtol, lsqr->exact_norm ? "exact" : "approx.", (double)lsqr->anorm, (double)rnorm, n));
595: *reason = KSP_CONVERGED_RTOL_NORMAL;
596: }
597: PetscFunctionReturn(PETSC_SUCCESS);
598: }
600: /*MC
601: KSPLSQR - Implements LSQR
603: Options Database Keys:
604: + -ksp_lsqr_set_standard_error - set standard error estimates of solution, see `KSPLSQRSetComputeStandardErrorVec()` and `KSPLSQRGetStandardErrorVec()`
605: . -ksp_lsqr_exact_mat_norm - compute exact matrix norm instead of iteratively refined estimate, see `KSPLSQRSetExactMatNorm()`
606: - -ksp_lsqr_monitor - monitor residual norm, norm of residual of normal equations A'*A x = A' b, and estimate of matrix norm ||A||
608: Level: beginner
610: Notes:
611: Supports non-square (rectangular) matrices.
613: This variant, when applied with no preconditioning is identical to the original algorithm in exact arithmetic; however, in practice, with no preconditioning
614: due to inexact arithmetic, it can converge differently. Hence when no preconditioner is used (`PCType` `PCNONE`) it automatically reverts to the original algorithm.
616: With the PETSc built-in preconditioners, such as `PCICC`, one should call `KSPSetOperators`(ksp,A,A'*A)) since the preconditioner needs to work
617: for the normal equations A'*A.
619: Supports only left preconditioning.
621: 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()`.
623: In exact arithmetic the LSQR method (with no preconditioning) is identical to the `KSPCG` algorithm applied to the normal equations.
624: The preconditioned variant was implemented by Bas van't Hof and is essentially a left preconditioning for the Normal Equations.
625: It appears the implementation with preconditioning tracks the true norm of the residual and uses that in the convergence test.
627: Developer Note:
628: How is this related to the `KSPCGNE` implementation? One difference is that `KSPCGNE` applies
629: the preconditioner transpose times the preconditioner, so one does not need to pass A'*A as the third argument to `KSPSetOperators()`.
631: Reference:
632: . * - The original unpreconditioned algorithm can be found in Paige and Saunders, ACM Transactions on Mathematical Software, Vol 8, 1982.
634: .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPSolve()`, `KSPLSQRConvergedDefault()`, `KSPLSQRSetComputeStandardErrorVec()`, `KSPLSQRGetStandardErrorVec()`, `KSPLSQRSetExactMatNorm()`, `KSPLSQRMonitorResidualDrawLGCreate()`, `KSPLSQRMonitorResidualDrawLG()`, `KSPLSQRMonitorResidual()`
635: M*/
636: PETSC_EXTERN PetscErrorCode KSPCreate_LSQR(KSP ksp)
637: {
638: KSP_LSQR *lsqr;
639: void *ctx;
641: PetscFunctionBegin;
642: PetscCall(PetscNew(&lsqr));
643: lsqr->se = NULL;
644: lsqr->se_flg = PETSC_FALSE;
645: lsqr->exact_norm = PETSC_FALSE;
646: lsqr->anorm = -1.0;
647: lsqr->arnorm = -1.0;
648: ksp->data = (void *)lsqr;
649: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 3));
651: ksp->ops->setup = KSPSetUp_LSQR;
652: ksp->ops->solve = KSPSolve_LSQR;
653: ksp->ops->destroy = KSPDestroy_LSQR;
654: ksp->ops->setfromoptions = KSPSetFromOptions_LSQR;
655: ksp->ops->view = KSPView_LSQR;
657: /* Backup current convergence test; remove destroy routine from KSP to prevent destroying the convergence context in KSPSetConvergenceTest() */
658: PetscCall(KSPGetAndClearConvergenceTest(ksp, &lsqr->converged, &lsqr->cnvP, &lsqr->convergeddestroy));
659: /* Override current convergence test */
660: PetscCall(KSPConvergedDefaultCreate(&ctx));
661: PetscCall(KSPSetConvergenceTest(ksp, KSPLSQRConvergedDefault, ctx, KSPConvergedDefaultDestroy));
662: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPLSQRMonitorResidual_C", KSPLSQRMonitorResidual_LSQR));
663: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPLSQRMonitorResidualDrawLG_C", KSPLSQRMonitorResidualDrawLG_LSQR));
664: PetscFunctionReturn(PETSC_SUCCESS);
665: }