Actual source code: lsqr.c
petsc-3.4.5 2014-06-29
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 <../src/ksp/ksp/impls/lsqr/lsqr.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: PetscReal arnorm; /* Norm of the vector A.r */
17: PetscReal anorm; /* Frobenius norm of the matrix A */
18: PetscReal rhs_norm; /* Norm of the right hand side */
19: } KSP_LSQR;
21: extern PetscErrorCode VecSquare(Vec);
25: static PetscErrorCode KSPSetUp_LSQR(KSP ksp)
26: {
28: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
29: PetscBool nopreconditioner;
32: PetscObjectTypeCompare((PetscObject)ksp->pc,PCNONE,&nopreconditioner);
33: /* nopreconditioner =PETSC_FALSE; */
35: lsqr->nwork_m = 2;
36: if (lsqr->vwork_m) {
37: VecDestroyVecs(lsqr->nwork_m,&lsqr->vwork_m);
38: }
39: if (nopreconditioner) lsqr->nwork_n = 4;
40: else lsqr->nwork_n = 5;
42: if (lsqr->vwork_n) {
43: VecDestroyVecs(lsqr->nwork_n,&lsqr->vwork_n);
44: }
45: KSPGetVecs(ksp,lsqr->nwork_n,&lsqr->vwork_n,lsqr->nwork_m,&lsqr->vwork_m);
46: if (lsqr->se_flg && !lsqr->se) {
47: /* lsqr->se is not set by user, get it from pmat */
48: Vec *se;
49: KSPGetVecs(ksp,1,&se,0,NULL);
50: lsqr->se = *se;
51: PetscFree(se);
52: }
53: return(0);
54: }
58: static PetscErrorCode KSPSolve_LSQR(KSP ksp)
59: {
61: PetscInt i,size1,size2;
62: PetscScalar rho,rhobar,phi,phibar,theta,c,s,tmp,tau;
63: PetscReal beta,alpha,rnorm;
64: Vec X,B,V,V1,U,U1,TMP,W,W2,SE,Z = NULL;
65: Mat Amat,Pmat;
66: MatStructure pflag;
67: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
68: PetscBool diagonalscale,nopreconditioner;
71: PCGetDiagonalScale(ksp->pc,&diagonalscale);
72: if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);
74: PCGetOperators(ksp->pc,&Amat,&Pmat,&pflag);
75: PetscObjectTypeCompare((PetscObject)ksp->pc,PCNONE,&nopreconditioner);
77: /* nopreconditioner =PETSC_FALSE; */
78: /* Calculate norm of right hand side */
79: VecNorm(ksp->vec_rhs,NORM_2,&lsqr->rhs_norm);
81: /* mark norm of matrix with negative number to indicate it has not yet been computed */
82: lsqr->anorm = -1.0;
84: /* vectors of length m, where system size is mxn */
85: B = ksp->vec_rhs;
86: U = lsqr->vwork_m[0];
87: U1 = lsqr->vwork_m[1];
89: /* vectors of length n */
90: X = ksp->vec_sol;
91: W = lsqr->vwork_n[0];
92: V = lsqr->vwork_n[1];
93: V1 = lsqr->vwork_n[2];
94: W2 = lsqr->vwork_n[3];
95: if (!nopreconditioner) Z = lsqr->vwork_n[4];
97: /* standard error vector */
98: SE = lsqr->se;
99: if (SE) {
100: VecGetSize(SE,&size1);
101: VecGetSize(X,&size2);
102: if (size1 != size2) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Standard error vector (size %d) does not match solution vector (size %d)",size1,size2);
103: VecSet(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: PetscObjectAMSTakeAccess((PetscObject)ksp);
117: ksp->its = 0;
118: ksp->rnorm = rnorm;
119: PetscObjectAMSGrantAccess((PetscObject)ksp);
120: KSPLogResidualHistory(ksp,rnorm);
121: KSPMonitor(ksp,0,rnorm);
122: (*ksp->converged)(ksp,0,rnorm,&ksp->reason,ksp->cnvP);
123: if (ksp->reason) return(0);
125: beta = rnorm;
126: VecScale(U,1.0/beta);
127: KSP_MatMultTranspose(ksp,Amat,U,V);
128: if (nopreconditioner) {
129: VecNorm(V,NORM_2,&alpha);
130: } else {
131: PCApply(ksp->pc,V,Z);
132: VecDotRealPart(V,Z,&alpha);
133: if (alpha <= 0.0) {
134: ksp->reason = KSP_DIVERGED_BREAKDOWN;
135: return(0);
136: }
137: alpha = PetscSqrtReal(alpha);
138: VecScale(Z,1.0/alpha);
139: }
140: VecScale(V,1.0/alpha);
142: if (nopreconditioner) {
143: VecCopy(V,W);
144: } else {
145: VecCopy(Z,W);
146: }
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: if (beta == 0.0) {
161: ksp->reason = KSP_DIVERGED_BREAKDOWN;
162: break;
163: }
164: VecScale(U1,1.0/beta);
166: KSP_MatMultTranspose(ksp,Amat,U1,V1);
167: VecAXPY(V1,-beta,V);
168: if (nopreconditioner) {
169: VecNorm(V1,NORM_2,&alpha);
170: } else {
171: PCApply(ksp->pc,V1,Z);
172: VecDotRealPart(V1,Z,&alpha);
173: if (alpha <= 0.0) {
174: ksp->reason = KSP_DIVERGED_BREAKDOWN;
175: break;
176: }
177: alpha = PetscSqrtReal(alpha);
178: VecScale(Z,1.0/alpha);
179: }
180: VecScale(V1,1.0/alpha);
181: rho = PetscSqrtScalar(rhobar*rhobar + beta*beta);
182: c = rhobar / rho;
183: s = beta / rho;
184: theta = s * alpha;
185: rhobar = -c * alpha;
186: phi = c * phibar;
187: phibar = s * phibar;
188: tau = s * phi;
190: VecAXPY(X,phi/rho,W); /* x <- x + (phi/rho) w */
192: if (SE) {
193: VecCopy(W,W2);
194: VecSquare(W2);
195: VecScale(W2,1.0/(rho*rho));
196: VecAXPY(SE, 1.0, W2); /* SE <- SE + (w^2/rho^2) */
197: }
198: if (nopreconditioner) {
199: VecAYPX(W,-theta/rho,V1); /* w <- v - (theta/rho) w */
200: } else {
201: VecAYPX(W,-theta/rho,Z); /* w <- z - (theta/rho) w */
202: }
204: lsqr->arnorm = alpha*PetscAbsScalar(tau);
205: rnorm = PetscRealPart(phibar);
207: PetscObjectAMSTakeAccess((PetscObject)ksp);
208: ksp->its++;
209: ksp->rnorm = rnorm;
210: PetscObjectAMSGrantAccess((PetscObject)ksp);
211: KSPLogResidualHistory(ksp,rnorm);
212: KSPMonitor(ksp,i+1,rnorm);
213: (*ksp->converged)(ksp,i+1,rnorm,&ksp->reason,ksp->cnvP);
214: if (ksp->reason) break;
215: SWAP(U1,U,TMP);
216: SWAP(V1,V,TMP);
218: i++;
219: } while (i<ksp->max_it);
220: if (i >= ksp->max_it && !ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
222: /* Finish off the standard error estimates */
223: if (SE) {
224: tmp = 1.0;
225: MatGetSize(Amat,&size1,&size2);
226: if (size1 > size2) tmp = size1 - size2;
227: tmp = rnorm / PetscSqrtScalar(tmp);
228: VecSqrtAbs(SE);
229: VecScale(SE,tmp);
230: }
231: return(0);
232: }
237: PetscErrorCode KSPDestroy_LSQR(KSP ksp)
238: {
239: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
243: /* Free work vectors */
244: if (lsqr->vwork_n) {
245: VecDestroyVecs(lsqr->nwork_n,&lsqr->vwork_n);
246: }
247: if (lsqr->vwork_m) {
248: VecDestroyVecs(lsqr->nwork_m,&lsqr->vwork_m);
249: }
250: if (lsqr->se_flg) {
251: VecDestroy(&lsqr->se);
252: }
253: PetscFree(ksp->data);
254: return(0);
255: }
259: PetscErrorCode KSPLSQRSetStandardErrorVec(KSP ksp, Vec se)
260: {
261: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
265: VecDestroy(&lsqr->se);
266: lsqr->se = se;
267: return(0);
268: }
272: PetscErrorCode KSPLSQRGetStandardErrorVec(KSP ksp,Vec *se)
273: {
274: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
277: *se = lsqr->se;
278: return(0);
279: }
283: PetscErrorCode KSPLSQRGetArnorm(KSP ksp,PetscReal *arnorm, PetscReal *rhs_norm, PetscReal *anorm)
284: {
285: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
289: *arnorm = lsqr->arnorm;
290: if (anorm) {
291: if (lsqr->anorm < 0.0) {
292: PC pc;
293: Mat Amat;
294: KSPGetPC(ksp,&pc);
295: PCGetOperators(pc,&Amat,NULL,NULL);
296: MatNorm(Amat,NORM_FROBENIUS,&lsqr->anorm);
297: }
298: *anorm = lsqr->anorm;
299: }
300: if (rhs_norm) *rhs_norm = lsqr->rhs_norm;
301: return(0);
302: }
306: /*@C
307: KSPLSQRMonitorDefault - Print the residual norm at each iteration of the LSQR method and the norm of the residual of the normal equations A'*A x = A' b
309: Collective on KSP
311: Input Parameters:
312: + ksp - iterative context
313: . n - iteration number
314: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
315: - dummy - unused monitor context
317: Level: intermediate
319: .keywords: KSP, default, monitor, residual
321: .seealso: KSPMonitorSet(), KSPMonitorTrueResidualNorm(), KSPMonitorLGResidualNormCreate(), KSPMonitorDefault()
322: @*/
323: PetscErrorCode KSPLSQRMonitorDefault(KSP ksp,PetscInt n,PetscReal rnorm,void *dummy)
324: {
326: PetscViewer viewer = dummy ? (PetscViewer) dummy : PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ksp));
327: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
330: PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
331: if (((PetscObject)ksp)->prefix) {
332: PetscViewerASCIIPrintf(viewer," Residual norm and norm of normal equations for %s solve.\n",((PetscObject)ksp)->prefix);
333: }
334: if (!n) {
335: PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %14.12e\n",n,rnorm);
336: } else {
337: PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %14.12e Residual norm normal equations %14.12e\n",n,rnorm,lsqr->arnorm);
338: }
339: PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
340: return(0);
341: }
345: PetscErrorCode KSPSetFromOptions_LSQR(KSP ksp)
346: {
348: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
349: char monfilename[PETSC_MAX_PATH_LEN];
350: PetscViewer monviewer;
351: PetscBool flg;
354: PetscOptionsHead("KSP LSQR Options");
355: PetscOptionsName("-ksp_lsqr_set_standard_error","Set Standard Error Estimates of Solution","KSPLSQRSetStandardErrorVec",&lsqr->se_flg);
356: PetscOptionsString("-ksp_lsqr_monitor","Monitor residual norm and norm of residual of normal equations","KSPMonitorSet","stdout",monfilename,PETSC_MAX_PATH_LEN,&flg);
357: if (flg) {
358: PetscViewerASCIIOpen(PetscObjectComm((PetscObject)ksp),monfilename,&monviewer);
359: KSPMonitorSet(ksp,KSPLSQRMonitorDefault,monviewer,(PetscErrorCode (*)(void**))PetscViewerDestroy);
360: }
361: PetscOptionsTail();
362: return(0);
363: }
367: PetscErrorCode KSPView_LSQR(KSP ksp,PetscViewer viewer)
368: {
369: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
371: PetscBool iascii;
374: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
375: if (iascii) {
376: if (lsqr->se) {
377: PetscReal rnorm;
378: KSPLSQRGetStandardErrorVec(ksp,&lsqr->se);
379: VecNorm(lsqr->se,NORM_2,&rnorm);
380: PetscViewerASCIIPrintf(viewer," Norm of Standard Error %G, Iterations %D\n",rnorm,ksp->its);
381: }
382: }
383: return(0);
384: }
388: /*@C
389: KSPLSQRDefaultConverged - Determines convergence of the LSQR Krylov method. This calls KSPDefaultConverged() and if that does not determine convergence then checks
390: convergence for the least squares problem.
392: Collective on KSP
394: Input Parameters:
395: + ksp - iterative context
396: . n - iteration number
397: . rnorm - 2-norm residual value (may be estimated)
398: - ctx - convergence context which must be created by KSPDefaultConvergedCreate()
400: reason is set to:
401: + positive - if the iteration has converged;
402: . negative - if residual norm exceeds divergence threshold;
403: - 0 - otherwise.
405: Notes:
406: 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.
408: Level: intermediate
410: .keywords: KSP, default, convergence, residual
412: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPSkipConverged(), KSPConvergedReason, KSPGetConvergedReason(),
413: KSPDefaultConvergedSetUIRNorm(), KSPDefaultConvergedSetUMIRNorm(), KSPDefaultConvergedCreate(), KSPDefaultConvergedDestroy(), KSPDefaultConverged()
414: @*/
415: PetscErrorCode KSPLSQRDefaultConverged(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *ctx)
416: {
418: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
421: KSPDefaultConverged(ksp,n,rnorm,reason,ctx);
422: if (!n || *reason) return(0);
423: if (lsqr->arnorm/lsqr->rhs_norm < ksp->rtol) *reason = KSP_CONVERGED_RTOL_NORMAL;
424: if (lsqr->arnorm < ksp->abstol) *reason = KSP_CONVERGED_ATOL_NORMAL;
425: return(0);
426: }
430: /*MC
431: KSPLSQR - This implements LSQR
433: Options Database Keys:
434: + -ksp_lsqr_set_standard_error - Set Standard Error Estimates of Solution see KSPLSQRSetStandardErrorVec()
435: . -ksp_lsqr_monitor - Monitor residual norm and norm of residual of normal equations
436: - see KSPSolve()
438: Level: beginner
440: Notes:
441: This varient, when applied with no preconditioning is identical to the original algorithm in exact arithematic; however, in practice, with no preconditioning
442: due to inexact arithematic, it can converge differently. Hence when no preconditioner is used (PCType PCNONE) it automatically reverts to the original algorithm.
444: With the PETSc built-in preconditioners, such as ICC, one should call KSPSetOperators(ksp,A,A'*A,...) since the preconditioner needs to work
445: for the normal equations A'*A.
447: Supports only left preconditioning.
449: References:The original unpreconditioned algorithm can be found in Paige and Saunders, ACM Transactions on Mathematical Software, Vol 8, pp 43-71, 1982.
450: In exact arithmetic the LSQR method (with no preconditioning) is identical to the KSPCG algorithm applied to the normal equations.
451: The preconditioned varient was implemented by Bas van't Hof and is essentially a left preconditioning for the Normal Equations. It appears the implementation with preconditioner
452: track the true norm of the residual and uses that in the convergence test.
454: Developer Notes: How is this related to the KSPCGNE implementation? One difference is that KSPCGNE applies
455: the preconditioner transpose times the preconditioner, so one does not need to pass A'*A as the third argument to KSPSetOperators().
458: For least squares problems without a zero to A*x = b, there are additional convergence tests for the residual of the normal equations, A'*(b - Ax), see KSPLSQRDefaultConverged()
460: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPLSQRDefaultConverged()
462: M*/
465: PETSC_EXTERN PetscErrorCode KSPCreate_LSQR(KSP ksp)
466: {
467: KSP_LSQR *lsqr;
471: PetscNewLog(ksp,KSP_LSQR,&lsqr);
472: lsqr->se = NULL;
473: lsqr->se_flg = PETSC_FALSE;
474: lsqr->arnorm = 0.0;
475: ksp->data = (void*)lsqr;
476: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);
478: ksp->ops->setup = KSPSetUp_LSQR;
479: ksp->ops->solve = KSPSolve_LSQR;
480: ksp->ops->destroy = KSPDestroy_LSQR;
481: ksp->ops->buildsolution = KSPBuildSolutionDefault;
482: ksp->ops->buildresidual = KSPBuildResidualDefault;
483: ksp->ops->setfromoptions = KSPSetFromOptions_LSQR;
484: ksp->ops->view = KSPView_LSQR;
485: ksp->converged = KSPLSQRDefaultConverged;
486: return(0);
487: }
491: PetscErrorCode VecSquare(Vec v)
492: {
494: PetscScalar *x;
495: PetscInt i, n;
498: VecGetLocalSize(v, &n);
499: VecGetArray(v, &x);
500: for (i = 0; i < n; i++) x[i] *= x[i];
501: VecRestoreArray(v, &x);
502: return(0);
503: }