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) { c = a; a = b; b = c; }
7: #include <petsc/private/kspimpl.h>
8: #include <petscdraw.h>
10: typedef struct {
11: PetscInt nwork_n,nwork_m;
12: Vec *vwork_m; /* work vectors of length m, where the system is size m x n */
13: Vec *vwork_n; /* work vectors of length n */
14: Vec se; /* Optional standard error vector */
15: PetscBool se_flg; /* flag for -ksp_lsqr_set_standard_error */
16: PetscBool exact_norm; /* flag for -ksp_lsqr_exact_mat_norm */
17: PetscReal arnorm; /* Good estimate of norm((A*inv(Pmat))'*r), where r = A*x - b, used in specific stopping criterion */
18: PetscReal anorm; /* Poor estimate of norm(A*inv(Pmat),'fro') used in specific stopping criterion */
19: /* Backup previous convergence test */
20: PetscErrorCode (*converged)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*);
21: PetscErrorCode (*convergeddestroy)(void*);
22: void *cnvP;
23: } KSP_LSQR;
25: static PetscErrorCode VecSquare(Vec v)
26: {
27: PetscScalar *x;
28: PetscInt i, n;
30: VecGetLocalSize(v, &n);
31: VecGetArray(v, &x);
32: for (i = 0; i < n; i++) x[i] *= PetscConj(x[i]);
33: VecRestoreArray(v, &x);
34: return 0;
35: }
37: static PetscErrorCode KSPSetUp_LSQR(KSP ksp)
38: {
39: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
40: PetscBool nopreconditioner;
42: PetscObjectTypeCompare((PetscObject)ksp->pc,PCNONE,&nopreconditioner);
44: if (lsqr->vwork_m) {
45: VecDestroyVecs(lsqr->nwork_m,&lsqr->vwork_m);
46: }
48: if (lsqr->vwork_n) {
49: VecDestroyVecs(lsqr->nwork_n,&lsqr->vwork_n);
50: }
52: lsqr->nwork_m = 2;
53: if (nopreconditioner) lsqr->nwork_n = 4;
54: else lsqr->nwork_n = 5;
55: KSPCreateVecs(ksp,lsqr->nwork_n,&lsqr->vwork_n,lsqr->nwork_m,&lsqr->vwork_m);
57: if (lsqr->se_flg && !lsqr->se) {
58: VecDuplicate(lsqr->vwork_n[0],&lsqr->se);
59: VecSet(lsqr->se,PETSC_INFINITY);
60: } else if (!lsqr->se_flg) {
61: VecDestroy(&lsqr->se);
62: }
63: return 0;
64: }
66: static PetscErrorCode KSPSolve_LSQR(KSP ksp)
67: {
68: PetscInt i,size1,size2;
69: PetscScalar rho,rhobar,phi,phibar,theta,c,s,tmp,tau;
70: PetscReal beta,alpha,rnorm;
71: Vec X,B,V,V1,U,U1,TMP,W,W2,Z = NULL;
72: Mat Amat,Pmat;
73: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
74: PetscBool diagonalscale,nopreconditioner;
76: PCGetDiagonalScale(ksp->pc,&diagonalscale);
79: PCGetOperators(ksp->pc,&Amat,&Pmat);
80: PetscObjectTypeCompare((PetscObject)ksp->pc,PCNONE,&nopreconditioner);
82: /* vectors of length m, where system size is mxn */
83: B = ksp->vec_rhs;
84: U = lsqr->vwork_m[0];
85: U1 = lsqr->vwork_m[1];
87: /* vectors of length n */
88: X = ksp->vec_sol;
89: W = lsqr->vwork_n[0];
90: V = lsqr->vwork_n[1];
91: V1 = lsqr->vwork_n[2];
92: W2 = lsqr->vwork_n[3];
93: if (!nopreconditioner) Z = lsqr->vwork_n[4];
95: /* standard error vector */
96: if (lsqr->se) {
97: VecSet(lsqr->se,0.0);
98: }
100: /* Compute initial residual, temporarily use work vector u */
101: if (!ksp->guess_zero) {
102: KSP_MatMult(ksp,Amat,X,U); /* u <- b - Ax */
103: VecAYPX(U,-1.0,B);
104: } else {
105: VecCopy(B,U); /* u <- b (x is 0) */
106: }
108: /* Test for nothing to do */
109: VecNorm(U,NORM_2,&rnorm);
110: KSPCheckNorm(ksp,rnorm);
111: PetscObjectSAWsTakeAccess((PetscObject)ksp);
112: ksp->its = 0;
113: ksp->rnorm = rnorm;
114: PetscObjectSAWsGrantAccess((PetscObject)ksp);
115: KSPLogResidualHistory(ksp,rnorm);
116: KSPMonitor(ksp,0,rnorm);
117: (*ksp->converged)(ksp,0,rnorm,&ksp->reason,ksp->cnvP);
118: if (ksp->reason) return 0;
120: beta = rnorm;
121: VecScale(U,1.0/beta);
122: KSP_MatMultHermitianTranspose(ksp,Amat,U,V);
123: if (nopreconditioner) {
124: VecNorm(V,NORM_2,&alpha);
125: KSPCheckNorm(ksp,rnorm);
126: } else {
127: /* this is an application of the preconditioner for the normal equations; not the operator, see the manual page */
128: PCApply(ksp->pc,V,Z);
129: VecDotRealPart(V,Z,&alpha);
130: if (alpha <= 0.0) {
131: ksp->reason = KSP_DIVERGED_BREAKDOWN;
132: return 0;
133: }
134: alpha = PetscSqrtReal(alpha);
135: VecScale(Z,1.0/alpha);
136: }
137: VecScale(V,1.0/alpha);
139: if (nopreconditioner) {
140: VecCopy(V,W);
141: } else {
142: VecCopy(Z,W);
143: }
145: if (lsqr->exact_norm) {
146: MatNorm(Amat,NORM_FROBENIUS,&lsqr->anorm);
147: } else lsqr->anorm = 0.0;
149: lsqr->arnorm = alpha * beta;
150: phibar = beta;
151: rhobar = alpha;
152: i = 0;
153: do {
154: if (nopreconditioner) {
155: KSP_MatMult(ksp,Amat,V,U1);
156: } else {
157: KSP_MatMult(ksp,Amat,Z,U1);
158: }
159: VecAXPY(U1,-alpha,U);
160: VecNorm(U1,NORM_2,&beta);
161: KSPCheckNorm(ksp,beta);
162: if (beta > 0.0) {
163: VecScale(U1,1.0/beta); /* beta*U1 = Amat*V - alpha*U */
164: if (!lsqr->exact_norm) {
165: lsqr->anorm = PetscSqrtReal(PetscSqr(lsqr->anorm) + PetscSqr(alpha) + PetscSqr(beta));
166: }
167: }
169: KSP_MatMultHermitianTranspose(ksp,Amat,U1,V1);
170: VecAXPY(V1,-beta,V);
171: if (nopreconditioner) {
172: VecNorm(V1,NORM_2,&alpha);
173: KSPCheckNorm(ksp,alpha);
174: } else {
175: PCApply(ksp->pc,V1,Z);
176: VecDotRealPart(V1,Z,&alpha);
177: if (alpha <= 0.0) {
178: ksp->reason = KSP_DIVERGED_BREAKDOWN;
179: break;
180: }
181: alpha = PetscSqrtReal(alpha);
182: VecScale(Z,1.0/alpha);
183: }
184: 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: VecAXPY(X,phi/rho,W); /* x <- x + (phi/rho) w */
196: if (lsqr->se) {
197: VecCopy(W,W2);
198: VecSquare(W2);
199: VecScale(W2,1.0/(rho*rho));
200: VecAXPY(lsqr->se, 1.0, W2); /* lsqr->se <- lsqr->se + (w^2/rho^2) */
201: }
202: if (nopreconditioner) {
203: VecAYPX(W,-theta/rho,V1); /* w <- v - (theta/rho) w */
204: } else {
205: VecAYPX(W,-theta/rho,Z); /* w <- z - (theta/rho) w */
206: }
208: lsqr->arnorm = alpha*PetscAbsScalar(tau);
209: rnorm = PetscRealPart(phibar);
211: PetscObjectSAWsTakeAccess((PetscObject)ksp);
212: ksp->its++;
213: ksp->rnorm = rnorm;
214: PetscObjectSAWsGrantAccess((PetscObject)ksp);
215: KSPLogResidualHistory(ksp,rnorm);
216: KSPMonitor(ksp,i+1,rnorm);
217: (*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: MatGetSize(Amat,&size1,&size2);
230: if (size1 > size2) tmp = size1 - size2;
231: tmp = rnorm / PetscSqrtScalar(tmp);
232: VecSqrtAbs(lsqr->se);
233: VecScale(lsqr->se,tmp);
234: }
235: return 0;
236: }
238: PetscErrorCode KSPDestroy_LSQR(KSP ksp)
239: {
240: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
242: /* Free work vectors */
243: if (lsqr->vwork_n) {
244: VecDestroyVecs(lsqr->nwork_n,&lsqr->vwork_n);
245: }
246: if (lsqr->vwork_m) {
247: VecDestroyVecs(lsqr->nwork_m,&lsqr->vwork_m);
248: }
249: VecDestroy(&lsqr->se);
250: /* Revert convergence test */
251: KSPSetConvergenceTest(ksp,lsqr->converged,lsqr->cnvP,lsqr->convergeddestroy);
252: /* Free the KSP_LSQR context */
253: PetscFree(ksp->data);
254: PetscObjectComposeFunction((PetscObject)ksp,"KSPLSQRMonitorResidual_C",NULL);
255: PetscObjectComposeFunction((PetscObject)ksp,"KSPLSQRMonitorResidualDrawLG_C",NULL);
256: return 0;
257: }
259: /*@
260: KSPLSQRSetComputeStandardErrorVec - Compute vector of standard error estimates during KSPSolve_LSQR().
262: Not Collective
264: Input Parameters:
265: + ksp - iterative context
266: - flg - compute the vector of standard estimates or not
268: Developer notes:
269: Vaclav: I'm not sure whether this vector is useful for anything.
271: Level: intermediate
273: .seealso: KSPSolve(), KSPLSQR, KSPLSQRGetStandardErrorVec()
274: @*/
275: PetscErrorCode KSPLSQRSetComputeStandardErrorVec(KSP ksp, PetscBool flg)
276: {
277: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
279: lsqr->se_flg = flg;
280: return 0;
281: }
283: /*@
284: KSPLSQRSetExactMatNorm - Compute exact matrix norm instead of iteratively refined estimate.
286: Not Collective
288: Input Parameters:
289: + ksp - iterative context
290: - flg - compute exact matrix norm or not
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: Level: intermediate
299: .seealso: KSPSolve(), KSPLSQR, KSPLSQRGetNorms(), KSPLSQRConvergedDefault()
300: @*/
301: PetscErrorCode KSPLSQRSetExactMatNorm(KSP ksp, PetscBool flg)
302: {
303: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
305: lsqr->exact_norm = flg;
306: return 0;
307: }
309: /*@
310: KSPLSQRGetStandardErrorVec - Get vector of standard error estimates.
311: Only available if -ksp_lsqr_set_standard_error was set to true
312: or KSPLSQRSetComputeStandardErrorVec(ksp, PETSC_TRUE) was called.
313: Otherwise returns NULL.
315: Not Collective
317: Input Parameters:
318: . ksp - iterative context
320: Output Parameters:
321: . se - vector of standard estimates
323: Options Database Keys:
324: . -ksp_lsqr_set_standard_error - set standard error estimates of solution
326: Developer notes:
327: Vaclav: I'm not sure whether this vector is useful for anything.
329: Level: intermediate
331: .seealso: KSPSolve(), KSPLSQR, KSPLSQRSetComputeStandardErrorVec()
332: @*/
333: PetscErrorCode KSPLSQRGetStandardErrorVec(KSP ksp,Vec *se)
334: {
335: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
337: *se = lsqr->se;
338: return 0;
339: }
341: /*@
342: KSPLSQRGetNorms - Get norm estimates that LSQR computes internally during KSPSolve().
344: Not Collective
346: Input Parameter:
347: . ksp - iterative context
349: Output Parameters:
350: + arnorm - good estimate of norm((A*inv(Pmat))'*r), where r = A*x - b, used in specific stopping criterion
351: - anorm - poor estimate of norm(A*inv(Pmat),'fro') used in specific stopping criterion
353: Notes:
354: Output parameters are meaningful only after KSPSolve().
355: 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.
356: If -ksp_lsqr_exact_mat_norm is set or KSPLSQRSetExactMatNorm(ksp, PETSC_TRUE) called, then anorm is exact Frobenius norm.
358: Level: intermediate
360: .seealso: KSPSolve(), KSPLSQR, KSPLSQRSetExactMatNorm()
361: @*/
362: PetscErrorCode KSPLSQRGetNorms(KSP ksp,PetscReal *arnorm, PetscReal *anorm)
363: {
364: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
366: if (arnorm) *arnorm = lsqr->arnorm;
367: if (anorm) *anorm = lsqr->anorm;
368: return 0;
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: PetscObjectGetTabLevel((PetscObject) ksp, &tablevel);
381: PetscObjectGetOptionsPrefix((PetscObject) ksp, &prefix);
382: PetscStrncpy(normtype, KSPNormTypes[ksp->normtype], sizeof(normtype));
383: PetscStrtolower(normtype);
384: PetscViewerPushFormat(viewer, format);
385: PetscViewerASCIIAddTab(viewer, tablevel);
386: if (n == 0 && prefix) PetscViewerASCIIPrintf(viewer, " Residual norm, norm of normal equations, and matrix norm for %s solve.\n", prefix);
387: if (!n) {
388: PetscViewerASCIIPrintf(viewer,"%3D KSP resid norm %14.12e\n",n,(double)rnorm);
389: } else {
390: PetscViewerASCIIPrintf(viewer,"%3D 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);
391: }
392: PetscViewerASCIISubtractTab(viewer, tablevel);
393: PetscViewerPopFormat(viewer);
394: return 0;
395: }
397: /*@C
398: KSPLSQRMonitorResidual - Prints the residual norm, as well as the normal equation residual norm, at each iteration of an iterative solver.
400: Collective on ksp
402: Input Parameters:
403: + ksp - iterative context
404: . n - iteration number
405: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
406: - vf - The viewer context
408: Options Database Key:
409: . -ksp_lsqr_monitor - Activates KSPLSQRMonitorResidual()
411: Level: intermediate
413: .seealso: KSPMonitorSet(), KSPMonitorResidual(),KSPMonitorTrueResidualMaxNorm()
414: @*/
415: PetscErrorCode KSPLSQRMonitorResidual(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
416: {
420: PetscTryMethod(ksp, "KSPLSQRMonitorResidual_C", (KSP,PetscInt,PetscReal,PetscViewerAndFormat*), (ksp,n,rnorm,vf));
421: return 0;
422: }
424: PetscErrorCode KSPLSQRMonitorResidualDrawLG_LSQR(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
425: {
426: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
427: PetscViewer viewer = vf->viewer;
428: PetscViewerFormat format = vf->format;
429: PetscDrawLG lg = vf->lg;
430: KSPConvergedReason reason;
431: PetscReal x[2], y[2];
433: PetscViewerPushFormat(viewer, format);
434: if (!n) PetscDrawLGReset(lg);
435: x[0] = (PetscReal) n;
436: if (rnorm > 0.0) y[0] = PetscLog10Real(rnorm);
437: else y[0] = -15.0;
438: x[1] = (PetscReal) n;
439: if (lsqr->arnorm > 0.0) y[1] = PetscLog10Real(lsqr->arnorm);
440: else y[1] = -15.0;
441: PetscDrawLGAddPoint(lg, x, y);
442: KSPGetConvergedReason(ksp, &reason);
443: if (n <= 20 || !(n % 5) || reason) {
444: PetscDrawLGDraw(lg);
445: PetscDrawLGSave(lg);
446: }
447: PetscViewerPopFormat(viewer);
448: return 0;
449: }
451: /*@C
452: KSPLSQRMonitorResidualDrawLG - Plots the true residual norm at each iteration of an iterative solver.
454: Collective on ksp
456: Input Parameters:
457: + ksp - iterative context
458: . n - iteration number
459: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
460: - vf - The viewer context
462: Options Database Key:
463: . -ksp_lsqr_monitor draw::draw_lg - Activates KSPMonitorTrueResidualDrawLG()
465: Level: intermediate
467: .seealso: KSPMonitorSet(), KSPMonitorTrueResidual()
468: @*/
469: PetscErrorCode KSPLSQRMonitorResidualDrawLG(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
470: {
475: PetscTryMethod(ksp, "KSPLSQRMonitorResidualDrawLG_C", (KSP,PetscInt,PetscReal,PetscViewerAndFormat*), (ksp,n,rnorm,vf));
476: return 0;
477: }
479: /*@C
480: KSPLSQRMonitorResidualDrawLGCreate - Creates the plotter for the LSQR residual and normal eqn residual.
482: Collective on ksp
484: Input Parameters:
485: + viewer - The PetscViewer
486: . format - The viewer format
487: - ctx - An optional user context
489: Output Parameter:
490: . vf - The viewer context
492: Level: intermediate
494: .seealso: KSPMonitorSet(), KSPLSQRMonitorResidual()
495: @*/
496: PetscErrorCode KSPLSQRMonitorResidualDrawLGCreate(PetscViewer viewer, PetscViewerFormat format, void *ctx, PetscViewerAndFormat **vf)
497: {
498: const char *names[] = {"residual", "normal eqn residual"};
500: PetscViewerAndFormatCreate(viewer, format, vf);
501: (*vf)->data = ctx;
502: KSPMonitorLGCreate(PetscObjectComm((PetscObject) viewer), NULL, NULL, "Log Residual Norm", 2, names, PETSC_DECIDE, PETSC_DECIDE, 400, 300, &(*vf)->lg);
503: return 0;
504: }
506: PetscErrorCode KSPSetFromOptions_LSQR(PetscOptionItems *PetscOptionsObject,KSP ksp)
507: {
508: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
510: PetscOptionsHead(PetscOptionsObject,"KSP LSQR Options");
511: PetscOptionsBool("-ksp_lsqr_compute_standard_error","Set Standard Error Estimates of Solution","KSPLSQRSetComputeStandardErrorVec",lsqr->se_flg,&lsqr->se_flg,NULL);
512: PetscOptionsBool("-ksp_lsqr_exact_mat_norm","Compute exact matrix norm instead of iteratively refined estimate","KSPLSQRSetExactMatNorm",lsqr->exact_norm,&lsqr->exact_norm,NULL);
513: KSPMonitorSetFromOptions(ksp, "-ksp_lsqr_monitor", "lsqr_residual", NULL);
514: PetscOptionsTail();
515: return 0;
516: }
518: PetscErrorCode KSPView_LSQR(KSP ksp,PetscViewer viewer)
519: {
520: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
521: PetscBool iascii;
523: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
524: if (iascii) {
525: if (lsqr->se) {
526: PetscReal rnorm;
527: VecNorm(lsqr->se,NORM_2,&rnorm);
528: PetscViewerASCIIPrintf(viewer," norm of standard error %g, iterations %d\n",(double)rnorm,ksp->its);
529: } else {
530: PetscViewerASCIIPrintf(viewer," standard error not computed\n");
531: }
532: if (lsqr->exact_norm) {
533: PetscViewerASCIIPrintf(viewer," using exact matrix norm\n");
534: } else {
535: PetscViewerASCIIPrintf(viewer," using inexact matrix norm\n");
536: }
537: }
538: return 0;
539: }
541: /*@C
542: KSPLSQRConvergedDefault - Determines convergence of the LSQR Krylov method.
544: Collective on ksp
546: Input Parameters:
547: + ksp - iterative context
548: . n - iteration number
549: . rnorm - 2-norm residual value (may be estimated)
550: - ctx - convergence context which must be created by KSPConvergedDefaultCreate()
552: reason is set to:
553: + positive - if the iteration has converged;
554: . negative - if residual norm exceeds divergence threshold;
555: - 0 - otherwise.
557: Notes:
558: KSPConvergedDefault() is called first to check for convergence in A*x=b.
559: If that does not determine convergence then checks convergence for the least squares problem, i.e. in min{|b-A*x|}.
560: Possible convergence for the least squares problem (which is based on the residual of the normal equations) are KSP_CONVERGED_RTOL_NORMAL norm and KSP_CONVERGED_ATOL_NORMAL.
561: KSP_CONVERGED_RTOL_NORMAL is returned if ||A'*r|| < rtol * ||A|| * ||r||.
562: Matrix norm ||A|| is iteratively refined estimate, see KSPLSQRGetNorms().
563: This criterion is now largely compatible with that in MATLAB lsqr().
565: Level: intermediate
567: .seealso: KSPLSQR, KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(),
568: KSPConvergedDefaultSetUIRNorm(), KSPConvergedDefaultSetUMIRNorm(), KSPConvergedDefaultCreate(), KSPConvergedDefaultDestroy(), KSPConvergedDefault(), KSPLSQRGetNorms(), KSPLSQRSetExactMatNorm()
569: @*/
570: PetscErrorCode KSPLSQRConvergedDefault(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *ctx)
571: {
572: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
574: /* check for convergence in A*x=b */
575: KSPConvergedDefault(ksp,n,rnorm,reason,ctx);
576: if (!n || *reason) return 0;
578: /* check for convergence in min{|b-A*x|} */
579: if (lsqr->arnorm < ksp->abstol) {
580: PetscInfo(ksp,"LSQR solver has converged. Normal equation residual %14.12e is less than absolute tolerance %14.12e at iteration %D\n",(double)lsqr->arnorm,(double)ksp->abstol,n);
581: *reason = KSP_CONVERGED_ATOL_NORMAL;
582: } else if (lsqr->arnorm < ksp->rtol * lsqr->anorm * rnorm) {
583: 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 %D\n",(double)lsqr->arnorm,(double)ksp->rtol,lsqr->exact_norm?"exact":"approx.",(double)lsqr->anorm,(double)rnorm,n);
584: *reason = KSP_CONVERGED_RTOL_NORMAL;
585: }
586: return 0;
587: }
589: /*MC
590: KSPLSQR - This implements LSQR
592: Options Database Keys:
593: + -ksp_lsqr_set_standard_error - set standard error estimates of solution, see KSPLSQRSetComputeStandardErrorVec() and KSPLSQRGetStandardErrorVec()
594: . -ksp_lsqr_exact_mat_norm - compute exact matrix norm instead of iteratively refined estimate, see KSPLSQRSetExactMatNorm()
595: - -ksp_lsqr_monitor - monitor residual norm, norm of residual of normal equations A'*A x = A' b, and estimate of matrix norm ||A||
597: Level: beginner
599: Notes:
600: Supports non-square (rectangular) matrices.
602: This variant, when applied with no preconditioning is identical to the original algorithm in exact arithematic; however, in practice, with no preconditioning
603: due to inexact arithmetic, it can converge differently. Hence when no preconditioner is used (PCType PCNONE) it automatically reverts to the original algorithm.
605: With the PETSc built-in preconditioners, such as ICC, one should call KSPSetOperators(ksp,A,A'*A)) since the preconditioner needs to work
606: for the normal equations A'*A.
608: Supports only left preconditioning.
610: 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().
612: References:
613: . * - The original unpreconditioned algorithm can be found in Paige and Saunders, ACM Transactions on Mathematical Software, Vol 8, 1982.
615: In exact arithmetic the LSQR method (with no preconditioning) is identical to the KSPCG algorithm applied to the normal equations.
616: The preconditioned variant was implemented by Bas van't Hof and is essentially a left preconditioning for the Normal Equations. It appears the implementation with preconditioner
617: track the true norm of the residual and uses that in the convergence test.
619: Developer Notes:
620: How is this related to the KSPCGNE implementation? One difference is that KSPCGNE applies
621: the preconditioner transpose times the preconditioner, so one does not need to pass A'*A as the third argument to KSPSetOperators().
623: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPSolve(), KSPLSQRConvergedDefault(), KSPLSQRSetComputeStandardErrorVec(), KSPLSQRGetStandardErrorVec(), KSPLSQRSetExactMatNorm()
625: M*/
626: PETSC_EXTERN PetscErrorCode KSPCreate_LSQR(KSP ksp)
627: {
628: KSP_LSQR *lsqr;
629: void *ctx;
631: PetscNewLog(ksp,&lsqr);
632: lsqr->se = NULL;
633: lsqr->se_flg = PETSC_FALSE;
634: lsqr->exact_norm = PETSC_FALSE;
635: lsqr->anorm = -1.0;
636: lsqr->arnorm = -1.0;
637: ksp->data = (void*)lsqr;
638: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,3);
640: ksp->ops->setup = KSPSetUp_LSQR;
641: ksp->ops->solve = KSPSolve_LSQR;
642: ksp->ops->destroy = KSPDestroy_LSQR;
643: ksp->ops->setfromoptions = KSPSetFromOptions_LSQR;
644: ksp->ops->view = KSPView_LSQR;
646: /* Backup current convergence test; remove destroy routine from KSP to prevent destroying the convergence context in KSPSetConvergenceTest() */
647: KSPGetAndClearConvergenceTest(ksp,&lsqr->converged,&lsqr->cnvP,&lsqr->convergeddestroy);
648: /* Override current convergence test */
649: KSPConvergedDefaultCreate(&ctx);
650: KSPSetConvergenceTest(ksp,KSPLSQRConvergedDefault,ctx,KSPConvergedDefaultDestroy);
651: PetscObjectComposeFunction((PetscObject)ksp,"KSPLSQRMonitorResidual_C",KSPLSQRMonitorResidual_LSQR);
652: PetscObjectComposeFunction((PetscObject)ksp,"KSPLSQRMonitorResidualDrawLG_C",KSPLSQRMonitorResidualDrawLG_LSQR);
653: return 0;
654: }