Actual source code: lsqr.c
petsc-3.7.5 2017-01-01
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;
23: static PetscErrorCode VecSquare(Vec v)
24: {
26: PetscScalar *x;
27: 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: }
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);
47: /* nopreconditioner =PETSC_FALSE; */
49: lsqr->nwork_m = 2;
50: if (lsqr->vwork_m) {
51: VecDestroyVecs(lsqr->nwork_m,&lsqr->vwork_m);
52: }
53: if (nopreconditioner) lsqr->nwork_n = 4;
54: else lsqr->nwork_n = 5;
56: if (lsqr->vwork_n) {
57: VecDestroyVecs(lsqr->nwork_n,&lsqr->vwork_n);
58: }
59: KSPCreateVecs(ksp,lsqr->nwork_n,&lsqr->vwork_n,lsqr->nwork_m,&lsqr->vwork_m);
60: if (lsqr->se_flg && !lsqr->se) {
61: /* lsqr->se is not set by user, get it from pmat */
62: Vec *se;
63: KSPCreateVecs(ksp,1,&se,0,NULL);
64: lsqr->se = *se;
65: PetscFree(se);
66: }
67: return(0);
68: }
72: static PetscErrorCode KSPSolve_LSQR(KSP ksp)
73: {
75: PetscInt i,size1,size2;
76: PetscScalar rho,rhobar,phi,phibar,theta,c,s,tmp,tau;
77: PetscReal beta,alpha,rnorm;
78: Vec X,B,V,V1,U,U1,TMP,W,W2,SE,Z = NULL;
79: Mat Amat,Pmat;
80: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
81: PetscBool diagonalscale,nopreconditioner;
84: PCGetDiagonalScale(ksp->pc,&diagonalscale);
85: if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);
87: PCGetOperators(ksp->pc,&Amat,&Pmat);
88: PetscObjectTypeCompare((PetscObject)ksp->pc,PCNONE,&nopreconditioner);
90: /* nopreconditioner =PETSC_FALSE; */
91: /* Calculate norm of right hand side */
92: VecNorm(ksp->vec_rhs,NORM_2,&lsqr->rhs_norm);
94: /* mark norm of matrix with negative number to indicate it has not yet been computed */
95: lsqr->anorm = -1.0;
97: /* vectors of length m, where system size is mxn */
98: B = ksp->vec_rhs;
99: U = lsqr->vwork_m[0];
100: U1 = lsqr->vwork_m[1];
102: /* vectors of length n */
103: X = ksp->vec_sol;
104: W = lsqr->vwork_n[0];
105: V = lsqr->vwork_n[1];
106: V1 = lsqr->vwork_n[2];
107: W2 = lsqr->vwork_n[3];
108: if (!nopreconditioner) Z = lsqr->vwork_n[4];
110: /* standard error vector */
111: SE = lsqr->se;
112: if (SE) {
113: VecGetSize(SE,&size1);
114: VecGetSize(X,&size2);
115: 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);
116: VecSet(SE,0.0);
117: }
119: /* Compute initial residual, temporarily use work vector u */
120: if (!ksp->guess_zero) {
121: KSP_MatMult(ksp,Amat,X,U); /* u <- b - Ax */
122: VecAYPX(U,-1.0,B);
123: } else {
124: VecCopy(B,U); /* u <- b (x is 0) */
125: }
127: /* Test for nothing to do */
128: VecNorm(U,NORM_2,&rnorm);
129: PetscObjectSAWsTakeAccess((PetscObject)ksp);
130: ksp->its = 0;
131: ksp->rnorm = rnorm;
132: PetscObjectSAWsGrantAccess((PetscObject)ksp);
133: KSPLogResidualHistory(ksp,rnorm);
134: KSPMonitor(ksp,0,rnorm);
135: (*ksp->converged)(ksp,0,rnorm,&ksp->reason,ksp->cnvP);
136: if (ksp->reason) return(0);
138: beta = rnorm;
139: VecScale(U,1.0/beta);
140: KSP_MatMultTranspose(ksp,Amat,U,V);
141: if (nopreconditioner) {
142: VecNorm(V,NORM_2,&alpha);
143: } else {
144: PCApply(ksp->pc,V,Z);
145: VecDotRealPart(V,Z,&alpha);
146: if (alpha <= 0.0) {
147: ksp->reason = KSP_DIVERGED_BREAKDOWN;
148: return(0);
149: }
150: alpha = PetscSqrtReal(alpha);
151: VecScale(Z,1.0/alpha);
152: }
153: VecScale(V,1.0/alpha);
155: if (nopreconditioner) {
156: VecCopy(V,W);
157: } else {
158: VecCopy(Z,W);
159: }
161: lsqr->arnorm = alpha * beta;
162: phibar = beta;
163: rhobar = alpha;
164: i = 0;
165: do {
166: if (nopreconditioner) {
167: KSP_MatMult(ksp,Amat,V,U1);
168: } else {
169: KSP_MatMult(ksp,Amat,Z,U1);
170: }
171: VecAXPY(U1,-alpha,U);
172: VecNorm(U1,NORM_2,&beta);
173: if (beta > 0.0) {
174: VecScale(U1,1.0/beta); /* beta*U1 = Amat*V - alpha*U */
175: }
177: KSP_MatMultTranspose(ksp,Amat,U1,V1);
178: VecAXPY(V1,-beta,V);
179: if (nopreconditioner) {
180: VecNorm(V1,NORM_2,&alpha);
181: } else {
182: PCApply(ksp->pc,V1,Z);
183: VecDotRealPart(V1,Z,&alpha);
184: if (alpha <= 0.0) {
185: ksp->reason = KSP_DIVERGED_BREAKDOWN;
186: break;
187: }
188: alpha = PetscSqrtReal(alpha);
189: VecScale(Z,1.0/alpha);
190: }
191: VecScale(V1,1.0/alpha); /* alpha*V1 = Amat^T*U1 - beta*V */
192: rho = PetscSqrtScalar(rhobar*rhobar + beta*beta);
193: c = rhobar / rho;
194: s = beta / rho;
195: theta = s * alpha;
196: rhobar = -c * alpha;
197: phi = c * phibar;
198: phibar = s * phibar;
199: tau = s * phi;
201: VecAXPY(X,phi/rho,W); /* x <- x + (phi/rho) w */
203: if (SE) {
204: VecCopy(W,W2);
205: VecSquare(W2);
206: VecScale(W2,1.0/(rho*rho));
207: VecAXPY(SE, 1.0, W2); /* SE <- SE + (w^2/rho^2) */
208: }
209: if (nopreconditioner) {
210: VecAYPX(W,-theta/rho,V1); /* w <- v - (theta/rho) w */
211: } else {
212: VecAYPX(W,-theta/rho,Z); /* w <- z - (theta/rho) w */
213: }
215: lsqr->arnorm = alpha*PetscAbsScalar(tau);
216: rnorm = PetscRealPart(phibar);
218: PetscObjectSAWsTakeAccess((PetscObject)ksp);
219: ksp->its++;
220: ksp->rnorm = rnorm;
221: PetscObjectSAWsGrantAccess((PetscObject)ksp);
222: KSPLogResidualHistory(ksp,rnorm);
223: KSPMonitor(ksp,i+1,rnorm);
224: (*ksp->converged)(ksp,i+1,rnorm,&ksp->reason,ksp->cnvP);
225: if (ksp->reason) break;
226: SWAP(U1,U,TMP);
227: SWAP(V1,V,TMP);
229: i++;
230: } while (i<ksp->max_it);
231: if (i >= ksp->max_it && !ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
233: /* Finish off the standard error estimates */
234: if (SE) {
235: tmp = 1.0;
236: MatGetSize(Amat,&size1,&size2);
237: if (size1 > size2) tmp = size1 - size2;
238: tmp = rnorm / PetscSqrtScalar(tmp);
239: VecSqrtAbs(SE);
240: VecScale(SE,tmp);
241: }
242: return(0);
243: }
248: PetscErrorCode KSPDestroy_LSQR(KSP ksp)
249: {
250: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
254: /* Free work vectors */
255: if (lsqr->vwork_n) {
256: VecDestroyVecs(lsqr->nwork_n,&lsqr->vwork_n);
257: }
258: if (lsqr->vwork_m) {
259: VecDestroyVecs(lsqr->nwork_m,&lsqr->vwork_m);
260: }
261: if (lsqr->se_flg) {
262: VecDestroy(&lsqr->se);
263: }
264: PetscFree(ksp->data);
265: return(0);
266: }
270: PetscErrorCode KSPLSQRSetStandardErrorVec(KSP ksp, Vec se)
271: {
272: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
276: VecDestroy(&lsqr->se);
277: lsqr->se = se;
278: return(0);
279: }
283: PetscErrorCode KSPLSQRGetStandardErrorVec(KSP ksp,Vec *se)
284: {
285: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
288: *se = lsqr->se;
289: return(0);
290: }
294: PetscErrorCode KSPLSQRGetArnorm(KSP ksp,PetscReal *arnorm, PetscReal *rhs_norm, PetscReal *anorm)
295: {
296: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
300: *arnorm = lsqr->arnorm;
301: if (anorm) {
302: if (lsqr->anorm < 0.0) {
303: PC pc;
304: Mat Amat;
305: KSPGetPC(ksp,&pc);
306: PCGetOperators(pc,&Amat,NULL);
307: MatNorm(Amat,NORM_FROBENIUS,&lsqr->anorm);
308: }
309: *anorm = lsqr->anorm;
310: }
311: if (rhs_norm) *rhs_norm = lsqr->rhs_norm;
312: return(0);
313: }
317: /*@C
318: 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
320: Collective on KSP
322: Input Parameters:
323: + ksp - iterative context
324: . n - iteration number
325: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
326: - dummy - viewer and format context
328: Level: intermediate
330: .keywords: KSP, default, monitor, residual
332: .seealso: KSPMonitorSet(), KSPMonitorTrueResidualNorm(), KSPMonitorLGResidualNormCreate(), KSPMonitorDefault()
333: @*/
334: PetscErrorCode KSPLSQRMonitorDefault(KSP ksp,PetscInt n,PetscReal rnorm,PetscViewerAndFormat *dummy)
335: {
337: PetscViewer viewer = dummy->viewer;
338: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
341: PetscViewerPushFormat(viewer,dummy->format);
342: PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
343: if (((PetscObject)ksp)->prefix) {
344: PetscViewerASCIIPrintf(viewer," Residual norm and norm of normal equations for %s solve.\n",((PetscObject)ksp)->prefix);
345: }
346: if (!n) {
347: PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %14.12e\n",n,rnorm);
348: } else {
349: PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %14.12e Residual norm normal equations %14.12e\n",n,rnorm,lsqr->arnorm);
350: }
351: PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
352: PetscViewerPopFormat(viewer);
353: return(0);
354: }
358: PetscErrorCode KSPSetFromOptions_LSQR(PetscOptionItems *PetscOptionsObject,KSP ksp)
359: {
361: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
364: PetscOptionsHead(PetscOptionsObject,"KSP LSQR Options");
365: PetscOptionsName("-ksp_lsqr_set_standard_error","Set Standard Error Estimates of Solution","KSPLSQRSetStandardErrorVec",&lsqr->se_flg);
366: KSPMonitorSetFromOptions(ksp,"-ksp_lsqr_monitor","Monitor residual norm and norm of residual of normal equations","KSPMonitorSet",KSPLSQRMonitorDefault);
367: PetscOptionsTail();
368: return(0);
369: }
373: PetscErrorCode KSPView_LSQR(KSP ksp,PetscViewer viewer)
374: {
375: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
377: PetscBool iascii;
380: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
381: if (iascii) {
382: if (lsqr->se) {
383: PetscReal rnorm;
384: KSPLSQRGetStandardErrorVec(ksp,&lsqr->se);
385: VecNorm(lsqr->se,NORM_2,&rnorm);
386: PetscViewerASCIIPrintf(viewer," Norm of Standard Error %g, Iterations %D\n",(double)rnorm,ksp->its);
387: }
388: }
389: return(0);
390: }
394: /*@C
395: KSPLSQRDefaultConverged - Determines convergence of the LSQR Krylov method. This calls KSPConvergedDefault() and if that does not determine convergence then checks
396: convergence for the least squares problem.
398: Collective on KSP
400: Input Parameters:
401: + ksp - iterative context
402: . n - iteration number
403: . rnorm - 2-norm residual value (may be estimated)
404: - ctx - convergence context which must be created by KSPConvergedDefaultCreate()
406: reason is set to:
407: + positive - if the iteration has converged;
408: . negative - if residual norm exceeds divergence threshold;
409: - 0 - otherwise.
411: Notes:
412: 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.
414: Level: intermediate
416: .keywords: KSP, default, convergence, residual
418: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(),
419: KSPConvergedDefaultSetUIRNorm(), KSPConvergedDefaultSetUMIRNorm(), KSPConvergedDefaultCreate(), KSPConvergedDefaultDestroy(), KSPConvergedDefault()
420: @*/
421: PetscErrorCode KSPLSQRDefaultConverged(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *ctx)
422: {
424: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
427: KSPConvergedDefault(ksp,n,rnorm,reason,ctx);
428: if (!n || *reason) return(0);
429: if (lsqr->arnorm/lsqr->rhs_norm < ksp->rtol) *reason = KSP_CONVERGED_RTOL_NORMAL;
430: if (lsqr->arnorm < ksp->abstol) *reason = KSP_CONVERGED_ATOL_NORMAL;
431: return(0);
432: }
434: /*MC
435: KSPLSQR - This implements LSQR
437: Options Database Keys:
438: + -ksp_lsqr_set_standard_error - Set Standard Error Estimates of Solution see KSPLSQRSetStandardErrorVec()
439: . -ksp_lsqr_monitor - Monitor residual norm and norm of residual of normal equations
440: - see KSPSolve()
442: Level: beginner
444: Notes:
445: Supports non-square (rectangular) matrices.
447: This varient, when applied with no preconditioning is identical to the original algorithm in exact arithematic; however, in practice, with no preconditioning
448: due to inexact arithematic, it can converge differently. Hence when no preconditioner is used (PCType PCNONE) it automatically reverts to the original algorithm.
450: With the PETSc built-in preconditioners, such as ICC, one should call KSPSetOperators(ksp,A,A'*A)) since the preconditioner needs to work
451: for the normal equations A'*A.
453: Supports only left preconditioning.
455: References:
456: . 1. - The original unpreconditioned algorithm can be found in Paige and Saunders, ACM Transactions on Mathematical Software, Vol 8, 1982.
458: In exact arithmetic the LSQR method (with no preconditioning) is identical to the KSPCG algorithm applied to the normal equations.
459: 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
460: track the true norm of the residual and uses that in the convergence test.
462: Developer Notes: How is this related to the KSPCGNE implementation? One difference is that KSPCGNE applies
463: the preconditioner transpose times the preconditioner, so one does not need to pass A'*A as the third argument to KSPSetOperators().
466: 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()
468: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPLSQRDefaultConverged()
470: M*/
473: PETSC_EXTERN PetscErrorCode KSPCreate_LSQR(KSP ksp)
474: {
475: KSP_LSQR *lsqr;
479: PetscNewLog(ksp,&lsqr);
480: lsqr->se = NULL;
481: lsqr->se_flg = PETSC_FALSE;
482: lsqr->arnorm = 0.0;
483: ksp->data = (void*)lsqr;
484: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,3);
486: ksp->ops->setup = KSPSetUp_LSQR;
487: ksp->ops->solve = KSPSolve_LSQR;
488: ksp->ops->destroy = KSPDestroy_LSQR;
489: ksp->ops->buildsolution = KSPBuildSolutionDefault;
490: ksp->ops->buildresidual = KSPBuildResidualDefault;
491: ksp->ops->setfromoptions = KSPSetFromOptions_LSQR;
492: ksp->ops->view = KSPView_LSQR;
493: ksp->converged = KSPLSQRDefaultConverged;
494: return(0);
495: }