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: {
28: PetscScalar *x;
29: PetscInt i, n;
32: VecGetLocalSize(v, &n);
33: VecGetArray(v, &x);
34: for (i = 0; i < n; i++) x[i] *= PetscConj(x[i]);
35: VecRestoreArray(v, &x);
36: return(0);
37: }
39: static PetscErrorCode KSPSetUp_LSQR(KSP ksp)
40: {
42: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
43: PetscBool nopreconditioner;
46: PetscObjectTypeCompare((PetscObject)ksp->pc,PCNONE,&nopreconditioner);
48: if (lsqr->vwork_m) {
49: VecDestroyVecs(lsqr->nwork_m,&lsqr->vwork_m);
50: }
52: if (lsqr->vwork_n) {
53: VecDestroyVecs(lsqr->nwork_n,&lsqr->vwork_n);
54: }
56: lsqr->nwork_m = 2;
57: if (nopreconditioner) lsqr->nwork_n = 4;
58: else lsqr->nwork_n = 5;
59: KSPCreateVecs(ksp,lsqr->nwork_n,&lsqr->vwork_n,lsqr->nwork_m,&lsqr->vwork_m);
61: if (lsqr->se_flg && !lsqr->se) {
62: VecDuplicate(lsqr->vwork_n[0],&lsqr->se);
63: VecSet(lsqr->se,PETSC_INFINITY);
64: } else if (!lsqr->se_flg) {
65: VecDestroy(&lsqr->se);
66: }
67: return(0);
68: }
70: static PetscErrorCode KSPSolve_LSQR(KSP ksp)
71: {
73: PetscInt i,size1,size2;
74: PetscScalar rho,rhobar,phi,phibar,theta,c,s,tmp,tau;
75: PetscReal beta,alpha,rnorm;
76: Vec X,B,V,V1,U,U1,TMP,W,W2,Z = NULL;
77: Mat Amat,Pmat;
78: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
79: PetscBool diagonalscale,nopreconditioner;
82: PCGetDiagonalScale(ksp->pc,&diagonalscale);
83: if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);
85: PCGetOperators(ksp->pc,&Amat,&Pmat);
86: PetscObjectTypeCompare((PetscObject)ksp->pc,PCNONE,&nopreconditioner);
88: /* vectors of length m, where system size is mxn */
89: B = ksp->vec_rhs;
90: U = lsqr->vwork_m[0];
91: U1 = lsqr->vwork_m[1];
93: /* vectors of length n */
94: X = ksp->vec_sol;
95: W = lsqr->vwork_n[0];
96: V = lsqr->vwork_n[1];
97: V1 = lsqr->vwork_n[2];
98: W2 = lsqr->vwork_n[3];
99: if (!nopreconditioner) Z = lsqr->vwork_n[4];
101: /* standard error vector */
102: if (lsqr->se) {
103: VecSet(lsqr->se,0.0);
104: }
106: /* Compute initial residual, temporarily use work vector u */
107: if (!ksp->guess_zero) {
108: KSP_MatMult(ksp,Amat,X,U); /* u <- b - Ax */
109: VecAYPX(U,-1.0,B);
110: } else {
111: VecCopy(B,U); /* u <- b (x is 0) */
112: }
114: /* Test for nothing to do */
115: VecNorm(U,NORM_2,&rnorm);
116: KSPCheckNorm(ksp,rnorm);
117: PetscObjectSAWsTakeAccess((PetscObject)ksp);
118: ksp->its = 0;
119: ksp->rnorm = rnorm;
120: PetscObjectSAWsGrantAccess((PetscObject)ksp);
121: KSPLogResidualHistory(ksp,rnorm);
122: KSPMonitor(ksp,0,rnorm);
123: (*ksp->converged)(ksp,0,rnorm,&ksp->reason,ksp->cnvP);
124: if (ksp->reason) return(0);
126: beta = rnorm;
127: VecScale(U,1.0/beta);
128: KSP_MatMultHermitianTranspose(ksp,Amat,U,V);
129: if (nopreconditioner) {
130: VecNorm(V,NORM_2,&alpha);
131: KSPCheckNorm(ksp,rnorm);
132: } else {
133: /* this is an application of the preconditioner for the normal equations; not the operator, see the manual page */
134: PCApply(ksp->pc,V,Z);
135: VecDotRealPart(V,Z,&alpha);
136: if (alpha <= 0.0) {
137: ksp->reason = KSP_DIVERGED_BREAKDOWN;
138: return(0);
139: }
140: alpha = PetscSqrtReal(alpha);
141: VecScale(Z,1.0/alpha);
142: }
143: VecScale(V,1.0/alpha);
145: if (nopreconditioner) {
146: VecCopy(V,W);
147: } else {
148: VecCopy(Z,W);
149: }
151: if (lsqr->exact_norm) {
152: MatNorm(Amat,NORM_FROBENIUS,&lsqr->anorm);
153: } else lsqr->anorm = 0.0;
155: lsqr->arnorm = alpha * beta;
156: phibar = beta;
157: rhobar = alpha;
158: i = 0;
159: do {
160: if (nopreconditioner) {
161: KSP_MatMult(ksp,Amat,V,U1);
162: } else {
163: KSP_MatMult(ksp,Amat,Z,U1);
164: }
165: VecAXPY(U1,-alpha,U);
166: VecNorm(U1,NORM_2,&beta);
167: KSPCheckNorm(ksp,beta);
168: if (beta > 0.0) {
169: VecScale(U1,1.0/beta); /* beta*U1 = Amat*V - alpha*U */
170: if (!lsqr->exact_norm) {
171: lsqr->anorm = PetscSqrtReal(PetscSqr(lsqr->anorm) + PetscSqr(alpha) + PetscSqr(beta));
172: }
173: }
175: KSP_MatMultHermitianTranspose(ksp,Amat,U1,V1);
176: VecAXPY(V1,-beta,V);
177: if (nopreconditioner) {
178: VecNorm(V1,NORM_2,&alpha);
179: KSPCheckNorm(ksp,alpha);
180: } else {
181: PCApply(ksp->pc,V1,Z);
182: VecDotRealPart(V1,Z,&alpha);
183: if (alpha <= 0.0) {
184: ksp->reason = KSP_DIVERGED_BREAKDOWN;
185: break;
186: }
187: alpha = PetscSqrtReal(alpha);
188: VecScale(Z,1.0/alpha);
189: }
190: VecScale(V1,1.0/alpha); /* alpha*V1 = Amat^T*U1 - beta*V */
191: rho = PetscSqrtScalar(rhobar*rhobar + beta*beta);
192: c = rhobar / rho;
193: s = beta / rho;
194: theta = s * alpha;
195: rhobar = -c * alpha;
196: phi = c * phibar;
197: phibar = s * phibar;
198: tau = s * phi;
200: VecAXPY(X,phi/rho,W); /* x <- x + (phi/rho) w */
202: if (lsqr->se) {
203: VecCopy(W,W2);
204: VecSquare(W2);
205: VecScale(W2,1.0/(rho*rho));
206: VecAXPY(lsqr->se, 1.0, W2); /* lsqr->se <- lsqr->se + (w^2/rho^2) */
207: }
208: if (nopreconditioner) {
209: VecAYPX(W,-theta/rho,V1); /* w <- v - (theta/rho) w */
210: } else {
211: VecAYPX(W,-theta/rho,Z); /* w <- z - (theta/rho) w */
212: }
214: lsqr->arnorm = alpha*PetscAbsScalar(tau);
215: rnorm = PetscRealPart(phibar);
217: PetscObjectSAWsTakeAccess((PetscObject)ksp);
218: ksp->its++;
219: ksp->rnorm = rnorm;
220: PetscObjectSAWsGrantAccess((PetscObject)ksp);
221: KSPLogResidualHistory(ksp,rnorm);
222: KSPMonitor(ksp,i+1,rnorm);
223: (*ksp->converged)(ksp,i+1,rnorm,&ksp->reason,ksp->cnvP);
224: if (ksp->reason) break;
225: SWAP(U1,U,TMP);
226: SWAP(V1,V,TMP);
228: i++;
229: } while (i<ksp->max_it);
230: if (i >= ksp->max_it && !ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
232: /* Finish off the standard error estimates */
233: if (lsqr->se) {
234: tmp = 1.0;
235: MatGetSize(Amat,&size1,&size2);
236: if (size1 > size2) tmp = size1 - size2;
237: tmp = rnorm / PetscSqrtScalar(tmp);
238: VecSqrtAbs(lsqr->se);
239: VecScale(lsqr->se,tmp);
240: }
241: return(0);
242: }
244: PetscErrorCode KSPDestroy_LSQR(KSP ksp)
245: {
246: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
250: /* Free work vectors */
251: if (lsqr->vwork_n) {
252: VecDestroyVecs(lsqr->nwork_n,&lsqr->vwork_n);
253: }
254: if (lsqr->vwork_m) {
255: VecDestroyVecs(lsqr->nwork_m,&lsqr->vwork_m);
256: }
257: VecDestroy(&lsqr->se);
258: /* Revert convergence test */
259: KSPSetConvergenceTest(ksp,lsqr->converged,lsqr->cnvP,lsqr->convergeddestroy);
260: /* Free the KSP_LSQR context */
261: PetscFree(ksp->data);
262: PetscObjectComposeFunction((PetscObject)ksp,"KSPLSQRMonitorResidual_C",NULL);
263: PetscObjectComposeFunction((PetscObject)ksp,"KSPLSQRMonitorResidualDrawLG_C",NULL);
264: return(0);
265: }
267: /*@
268: KSPLSQRSetComputeStandardErrorVec - Compute vector of standard error estimates during KSPSolve_LSQR().
270: Not Collective
272: Input Parameters:
273: + ksp - iterative context
274: - flg - compute the vector of standard estimates or not
276: Developer notes:
277: Vaclav: I'm not sure whether this vector is useful for anything.
279: Level: intermediate
281: .seealso: KSPSolve(), KSPLSQR, KSPLSQRGetStandardErrorVec()
282: @*/
283: PetscErrorCode KSPLSQRSetComputeStandardErrorVec(KSP ksp, PetscBool flg)
284: {
285: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
288: lsqr->se_flg = flg;
289: return(0);
290: }
292: /*@
293: KSPLSQRSetExactMatNorm - Compute exact matrix norm instead of iteratively refined estimate.
295: Not Collective
297: Input Parameters:
298: + ksp - iterative context
299: - flg - compute exact matrix norm or not
301: Notes:
302: By default, flg=PETSC_FALSE. This is usually preferred to avoid possibly expensive computation of the norm.
303: For flg=PETSC_TRUE, we call MatNorm(Amat,NORM_FROBENIUS,&lsqr->anorm) which will work only for some types of explicitly assembled matrices.
304: This can affect convergence rate as KSPLSQRConvergedDefault() assumes different value of ||A|| used in normal equation stopping criterion.
306: Level: intermediate
308: .seealso: KSPSolve(), KSPLSQR, KSPLSQRGetNorms(), KSPLSQRConvergedDefault()
309: @*/
310: PetscErrorCode KSPLSQRSetExactMatNorm(KSP ksp, PetscBool flg)
311: {
312: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
315: lsqr->exact_norm = flg;
316: return(0);
317: }
319: /*@
320: KSPLSQRGetStandardErrorVec - Get vector of standard error estimates.
321: Only available if -ksp_lsqr_set_standard_error was set to true
322: or KSPLSQRSetComputeStandardErrorVec(ksp, PETSC_TRUE) was called.
323: Otherwise returns NULL.
325: Not Collective
327: Input Parameters:
328: . ksp - iterative context
330: Output Parameters:
331: . se - vector of standard estimates
333: Options Database Keys:
334: . -ksp_lsqr_set_standard_error - set standard error estimates of solution
336: Developer notes:
337: Vaclav: I'm not sure whether this vector is useful for anything.
339: Level: intermediate
341: .seealso: KSPSolve(), KSPLSQR, KSPLSQRSetComputeStandardErrorVec()
342: @*/
343: PetscErrorCode KSPLSQRGetStandardErrorVec(KSP ksp,Vec *se)
344: {
345: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
348: *se = lsqr->se;
349: return(0);
350: }
352: /*@
353: KSPLSQRGetNorms - Get norm estimates that LSQR computes internally during KSPSolve().
355: Not Collective
357: Input Parameter:
358: . ksp - iterative context
360: Output Parameters:
361: + arnorm - good estimate of norm((A*inv(Pmat))'*r), where r = A*x - b, used in specific stopping criterion
362: - anorm - poor estimate of norm(A*inv(Pmat),'fro') used in specific stopping criterion
364: Notes:
365: Output parameters are meaningful only after KSPSolve().
366: 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.
367: If -ksp_lsqr_exact_mat_norm is set or KSPLSQRSetExactMatNorm(ksp, PETSC_TRUE) called, then anorm is exact Frobenius norm.
369: Level: intermediate
371: .seealso: KSPSolve(), KSPLSQR, KSPLSQRSetExactMatNorm()
372: @*/
373: PetscErrorCode KSPLSQRGetNorms(KSP ksp,PetscReal *arnorm, PetscReal *anorm)
374: {
375: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
378: if (arnorm) *arnorm = lsqr->arnorm;
379: if (anorm) *anorm = lsqr->anorm;
380: return(0);
381: }
383: PetscErrorCode KSPLSQRMonitorResidual_LSQR(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
384: {
385: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
386: PetscViewer viewer = vf->viewer;
387: PetscViewerFormat format = vf->format;
388: char normtype[256];
389: PetscInt tablevel;
390: const char *prefix;
391: PetscErrorCode ierr;
394: PetscObjectGetTabLevel((PetscObject) ksp, &tablevel);
395: PetscObjectGetOptionsPrefix((PetscObject) ksp, &prefix);
396: PetscStrncpy(normtype, KSPNormTypes[ksp->normtype], sizeof(normtype));
397: PetscStrtolower(normtype);
398: PetscViewerPushFormat(viewer, format);
399: PetscViewerASCIIAddTab(viewer, tablevel);
400: if (n == 0 && prefix) {PetscViewerASCIIPrintf(viewer, " Residual norm, norm of normal equations, and matrix norm for %s solve.\n", prefix);}
401: if (!n) {
402: PetscViewerASCIIPrintf(viewer,"%3D KSP resid norm %14.12e\n",n,(double)rnorm);
403: } else {
404: 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);
405: }
406: PetscViewerASCIISubtractTab(viewer, tablevel);
407: PetscViewerPopFormat(viewer);
408: return(0);
409: }
411: /*@C
412: KSPLSQRMonitorResidual - Prints the residual norm, as well as the normal equation residual norm, at each iteration of an iterative solver.
414: Collective on ksp
416: Input Parameters:
417: + ksp - iterative context
418: . n - iteration number
419: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
420: - vf - The viewer context
422: Options Database Key:
423: . -ksp_lsqr_monitor - Activates KSPLSQRMonitorResidual()
425: Level: intermediate
427: .seealso: KSPMonitorSet(), KSPMonitorResidual(),KSPMonitorTrueResidualMaxNorm()
428: @*/
429: PetscErrorCode KSPLSQRMonitorResidual(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
430: {
437: PetscTryMethod(ksp, "KSPLSQRMonitorResidual_C", (KSP,PetscInt,PetscReal,PetscViewerAndFormat*), (ksp,n,rnorm,vf));
438: return(0);
439: }
441: PetscErrorCode KSPLSQRMonitorResidualDrawLG_LSQR(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
442: {
443: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
444: PetscViewer viewer = vf->viewer;
445: PetscViewerFormat format = vf->format;
446: PetscDrawLG lg = vf->lg;
447: KSPConvergedReason reason;
448: PetscReal x[2], y[2];
449: PetscErrorCode ierr;
452: PetscViewerPushFormat(viewer, format);
453: if (!n) {PetscDrawLGReset(lg);}
454: x[0] = (PetscReal) n;
455: if (rnorm > 0.0) y[0] = PetscLog10Real(rnorm);
456: else y[0] = -15.0;
457: x[1] = (PetscReal) n;
458: if (lsqr->arnorm > 0.0) y[1] = PetscLog10Real(lsqr->arnorm);
459: else y[1] = -15.0;
460: PetscDrawLGAddPoint(lg, x, y);
461: KSPGetConvergedReason(ksp, &reason);
462: if (n <= 20 || !(n % 5) || reason) {
463: PetscDrawLGDraw(lg);
464: PetscDrawLGSave(lg);
465: }
466: PetscViewerPopFormat(viewer);
467: return(0);
468: }
470: /*@C
471: KSPLSQRMonitorResidualDrawLG - Plots the true residual norm at each iteration of an iterative solver.
473: Collective on ksp
475: Input Parameters:
476: + ksp - iterative context
477: . n - iteration number
478: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
479: - vf - The viewer context
481: Options Database Key:
482: . -ksp_lsqr_monitor draw::draw_lg - Activates KSPMonitorTrueResidualDrawLG()
484: Level: intermediate
486: .seealso: KSPMonitorSet(), KSPMonitorTrueResidual()
487: @*/
488: PetscErrorCode KSPLSQRMonitorResidualDrawLG(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
489: {
497: PetscTryMethod(ksp, "KSPLSQRMonitorResidualDrawLG_C", (KSP,PetscInt,PetscReal,PetscViewerAndFormat*), (ksp,n,rnorm,vf));
498: return(0);
499: }
501: /*@C
502: KSPLSQRMonitorResidualDrawLGCreate - Creates the plotter for the LSQR residual and normal eqn residual.
504: Collective on ksp
506: Input Parameters:
507: + viewer - The PetscViewer
508: . format - The viewer format
509: - ctx - An optional user context
511: Output Parameter:
512: . vf - The viewer context
514: Level: intermediate
516: .seealso: KSPMonitorSet(), KSPLSQRMonitorResidual()
517: @*/
518: PetscErrorCode KSPLSQRMonitorResidualDrawLGCreate(PetscViewer viewer, PetscViewerFormat format, void *ctx, PetscViewerAndFormat **vf)
519: {
520: const char *names[] = {"residual", "normal eqn residual"};
524: PetscViewerAndFormatCreate(viewer, format, vf);
525: (*vf)->data = ctx;
526: KSPMonitorLGCreate(PetscObjectComm((PetscObject) viewer), NULL, NULL, "Log Residual Norm", 2, names, PETSC_DECIDE, PETSC_DECIDE, 400, 300, &(*vf)->lg);
527: return(0);
528: }
530: PetscErrorCode KSPSetFromOptions_LSQR(PetscOptionItems *PetscOptionsObject,KSP ksp)
531: {
533: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
536: PetscOptionsHead(PetscOptionsObject,"KSP LSQR Options");
537: PetscOptionsBool("-ksp_lsqr_compute_standard_error","Set Standard Error Estimates of Solution","KSPLSQRSetComputeStandardErrorVec",lsqr->se_flg,&lsqr->se_flg,NULL);
538: PetscOptionsBool("-ksp_lsqr_exact_mat_norm","Compute exact matrix norm instead of iteratively refined estimate","KSPLSQRSetExactMatNorm",lsqr->exact_norm,&lsqr->exact_norm,NULL);
539: KSPMonitorSetFromOptions(ksp, "-ksp_lsqr_monitor", "lsqr_residual", NULL);
540: PetscOptionsTail();
541: return(0);
542: }
544: PetscErrorCode KSPView_LSQR(KSP ksp,PetscViewer viewer)
545: {
546: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
548: PetscBool iascii;
551: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
552: if (iascii) {
553: if (lsqr->se) {
554: PetscReal rnorm;
555: VecNorm(lsqr->se,NORM_2,&rnorm);
556: PetscViewerASCIIPrintf(viewer," norm of standard error %g, iterations %d\n",(double)rnorm,ksp->its);
557: } else {
558: PetscViewerASCIIPrintf(viewer," standard error not computed\n");
559: }
560: if (lsqr->exact_norm) {
561: PetscViewerASCIIPrintf(viewer," using exact matrix norm\n");
562: } else {
563: PetscViewerASCIIPrintf(viewer," using inexact matrix norm\n");
564: }
565: }
566: return(0);
567: }
569: /*@C
570: KSPLSQRConvergedDefault - Determines convergence of the LSQR Krylov method.
572: Collective on ksp
574: Input Parameters:
575: + ksp - iterative context
576: . n - iteration number
577: . rnorm - 2-norm residual value (may be estimated)
578: - ctx - convergence context which must be created by KSPConvergedDefaultCreate()
580: reason is set to:
581: + positive - if the iteration has converged;
582: . negative - if residual norm exceeds divergence threshold;
583: - 0 - otherwise.
585: Notes:
586: KSPConvergedDefault() is called first to check for convergence in A*x=b.
587: If that does not determine convergence then checks convergence for the least squares problem, i.e. in min{|b-A*x|}.
588: 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.
589: KSP_CONVERGED_RTOL_NORMAL is returned if ||A'*r|| < rtol * ||A|| * ||r||.
590: Matrix norm ||A|| is iteratively refined estimate, see KSPLSQRGetNorms().
591: This criterion is now largely compatible with that in MATLAB lsqr().
593: Level: intermediate
595: .seealso: KSPLSQR, KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(),
596: KSPConvergedDefaultSetUIRNorm(), KSPConvergedDefaultSetUMIRNorm(), KSPConvergedDefaultCreate(), KSPConvergedDefaultDestroy(), KSPConvergedDefault(), KSPLSQRGetNorms(), KSPLSQRSetExactMatNorm()
597: @*/
598: PetscErrorCode KSPLSQRConvergedDefault(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *ctx)
599: {
601: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
604: /* check for convergence in A*x=b */
605: KSPConvergedDefault(ksp,n,rnorm,reason,ctx);
606: if (!n || *reason) return(0);
608: /* check for convergence in min{|b-A*x|} */
609: if (lsqr->arnorm < ksp->abstol) {
610: PetscInfo3(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);
611: *reason = KSP_CONVERGED_ATOL_NORMAL;
612: } else if (lsqr->arnorm < ksp->rtol * lsqr->anorm * rnorm) {
613: PetscInfo6(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);
614: *reason = KSP_CONVERGED_RTOL_NORMAL;
615: }
616: return(0);
617: }
619: /*MC
620: KSPLSQR - This implements LSQR
622: Options Database Keys:
623: + -ksp_lsqr_set_standard_error - set standard error estimates of solution, see KSPLSQRSetComputeStandardErrorVec() and KSPLSQRGetStandardErrorVec()
624: . -ksp_lsqr_exact_mat_norm - compute exact matrix norm instead of iteratively refined estimate, see KSPLSQRSetExactMatNorm()
625: - -ksp_lsqr_monitor - monitor residual norm, norm of residual of normal equations A'*A x = A' b, and estimate of matrix norm ||A||
627: Level: beginner
629: Notes:
630: Supports non-square (rectangular) matrices.
632: This variant, when applied with no preconditioning is identical to the original algorithm in exact arithematic; however, in practice, with no preconditioning
633: due to inexact arithmetic, it can converge differently. Hence when no preconditioner is used (PCType PCNONE) it automatically reverts to the original algorithm.
635: With the PETSc built-in preconditioners, such as ICC, one should call KSPSetOperators(ksp,A,A'*A)) since the preconditioner needs to work
636: for the normal equations A'*A.
638: Supports only left preconditioning.
640: 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().
642: References:
643: . 1. - The original unpreconditioned algorithm can be found in Paige and Saunders, ACM Transactions on Mathematical Software, Vol 8, 1982.
645: In exact arithmetic the LSQR method (with no preconditioning) is identical to the KSPCG algorithm applied to the normal equations.
646: 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
647: track the true norm of the residual and uses that in the convergence test.
649: Developer Notes:
650: How is this related to the KSPCGNE implementation? One difference is that KSPCGNE applies
651: the preconditioner transpose times the preconditioner, so one does not need to pass A'*A as the third argument to KSPSetOperators().
653: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPSolve(), KSPLSQRConvergedDefault(), KSPLSQRSetComputeStandardErrorVec(), KSPLSQRGetStandardErrorVec(), KSPLSQRSetExactMatNorm()
655: M*/
656: PETSC_EXTERN PetscErrorCode KSPCreate_LSQR(KSP ksp)
657: {
658: KSP_LSQR *lsqr;
659: void *ctx;
663: PetscNewLog(ksp,&lsqr);
664: lsqr->se = NULL;
665: lsqr->se_flg = PETSC_FALSE;
666: lsqr->exact_norm = PETSC_FALSE;
667: lsqr->anorm = -1.0;
668: lsqr->arnorm = -1.0;
669: ksp->data = (void*)lsqr;
670: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,3);
672: ksp->ops->setup = KSPSetUp_LSQR;
673: ksp->ops->solve = KSPSolve_LSQR;
674: ksp->ops->destroy = KSPDestroy_LSQR;
675: ksp->ops->setfromoptions = KSPSetFromOptions_LSQR;
676: ksp->ops->view = KSPView_LSQR;
678: /* Backup current convergence test; remove destroy routine from KSP to prevent destroying the convergence context in KSPSetConvergenceTest() */
679: KSPGetAndClearConvergenceTest(ksp,&lsqr->converged,&lsqr->cnvP,&lsqr->convergeddestroy);
680: /* Override current convergence test */
681: KSPConvergedDefaultCreate(&ctx);
682: KSPSetConvergenceTest(ksp,KSPLSQRConvergedDefault,ctx,KSPConvergedDefaultDestroy);
683: PetscObjectComposeFunction((PetscObject)ksp,"KSPLSQRMonitorResidual_C",KSPLSQRMonitorResidual_LSQR);
684: PetscObjectComposeFunction((PetscObject)ksp,"KSPLSQRMonitorResidualDrawLG_C",KSPLSQRMonitorResidualDrawLG_LSQR);
685: return(0);
686: }