Actual source code: lsqr.c
petsc-3.9.4 2018-09-11
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>
9: typedef struct {
10: PetscInt nwork_n,nwork_m;
11: Vec *vwork_m; /* work vectors of length m, where the system is size m x n */
12: Vec *vwork_n; /* work vectors of length n */
13: Vec se; /* Optional standard error vector */
14: PetscBool se_flg; /* flag for -ksp_lsqr_set_standard_error */
15: PetscReal arnorm; /* Norm of the vector A.r */
16: PetscReal anorm; /* Frobenius norm of the matrix A */
17: PetscReal rhs_norm; /* Norm of the right hand side */
18: } KSP_LSQR;
20: static PetscErrorCode VecSquare(Vec v)
21: {
23: PetscScalar *x;
24: PetscInt i, n;
27: VecGetLocalSize(v, &n);
28: VecGetArray(v, &x);
29: for (i = 0; i < n; i++) x[i] *= PetscConj(x[i]);
30: VecRestoreArray(v, &x);
31: return(0);
32: }
34: static PetscErrorCode KSPSetUp_LSQR(KSP ksp)
35: {
37: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
38: PetscBool nopreconditioner;
41: PetscObjectTypeCompare((PetscObject)ksp->pc,PCNONE,&nopreconditioner);
42: /* nopreconditioner =PETSC_FALSE; */
44: lsqr->nwork_m = 2;
45: if (lsqr->vwork_m) {
46: VecDestroyVecs(lsqr->nwork_m,&lsqr->vwork_m);
47: }
48: if (nopreconditioner) lsqr->nwork_n = 4;
49: else lsqr->nwork_n = 5;
51: if (lsqr->vwork_n) {
52: VecDestroyVecs(lsqr->nwork_n,&lsqr->vwork_n);
53: }
54: KSPCreateVecs(ksp,lsqr->nwork_n,&lsqr->vwork_n,lsqr->nwork_m,&lsqr->vwork_m);
55: if (lsqr->se_flg && !lsqr->se) {
56: /* lsqr->se is not set by user, get it from pmat */
57: Vec *se;
58: KSPCreateVecs(ksp,1,&se,0,NULL);
59: lsqr->se = *se;
60: PetscFree(se);
61: }
62: return(0);
63: }
65: static PetscErrorCode KSPSolve_LSQR(KSP ksp)
66: {
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,SE,Z = NULL;
72: Mat Amat,Pmat;
73: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
74: PetscBool diagonalscale,nopreconditioner;
77: PCGetDiagonalScale(ksp->pc,&diagonalscale);
78: if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);
80: PCGetOperators(ksp->pc,&Amat,&Pmat);
81: PetscObjectTypeCompare((PetscObject)ksp->pc,PCNONE,&nopreconditioner);
83: /* nopreconditioner =PETSC_FALSE; */
84: /* Calculate norm of right hand side */
85: VecNorm(ksp->vec_rhs,NORM_2,&lsqr->rhs_norm);
87: /* mark norm of matrix with negative number to indicate it has not yet been computed */
88: lsqr->anorm = -1.0;
90: /* vectors of length m, where system size is mxn */
91: B = ksp->vec_rhs;
92: U = lsqr->vwork_m[0];
93: U1 = lsqr->vwork_m[1];
95: /* vectors of length n */
96: X = ksp->vec_sol;
97: W = lsqr->vwork_n[0];
98: V = lsqr->vwork_n[1];
99: V1 = lsqr->vwork_n[2];
100: W2 = lsqr->vwork_n[3];
101: if (!nopreconditioner) Z = lsqr->vwork_n[4];
103: /* standard error vector */
104: SE = lsqr->se;
105: if (SE) {
106: VecGetSize(SE,&size1);
107: VecGetSize(X,&size2);
108: 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);
109: VecSet(SE,0.0);
110: }
112: /* Compute initial residual, temporarily use work vector u */
113: if (!ksp->guess_zero) {
114: KSP_MatMult(ksp,Amat,X,U); /* u <- b - Ax */
115: VecAYPX(U,-1.0,B);
116: } else {
117: VecCopy(B,U); /* u <- b (x is 0) */
118: }
120: /* Test for nothing to do */
121: VecNorm(U,NORM_2,&rnorm);
122: PetscObjectSAWsTakeAccess((PetscObject)ksp);
123: ksp->its = 0;
124: ksp->rnorm = rnorm;
125: PetscObjectSAWsGrantAccess((PetscObject)ksp);
126: KSPLogResidualHistory(ksp,rnorm);
127: KSPMonitor(ksp,0,rnorm);
128: (*ksp->converged)(ksp,0,rnorm,&ksp->reason,ksp->cnvP);
129: if (ksp->reason) return(0);
131: beta = rnorm;
132: VecScale(U,1.0/beta);
133: KSP_MatMultTranspose(ksp,Amat,U,V);
134: if (nopreconditioner) {
135: VecNorm(V,NORM_2,&alpha);
136: } else {
137: PCApply(ksp->pc,V,Z);
138: VecDotRealPart(V,Z,&alpha);
139: if (alpha <= 0.0) {
140: ksp->reason = KSP_DIVERGED_BREAKDOWN;
141: return(0);
142: }
143: alpha = PetscSqrtReal(alpha);
144: VecScale(Z,1.0/alpha);
145: }
146: VecScale(V,1.0/alpha);
148: if (nopreconditioner) {
149: VecCopy(V,W);
150: } else {
151: VecCopy(Z,W);
152: }
154: lsqr->arnorm = alpha * beta;
155: phibar = beta;
156: rhobar = alpha;
157: i = 0;
158: do {
159: if (nopreconditioner) {
160: KSP_MatMult(ksp,Amat,V,U1);
161: } else {
162: KSP_MatMult(ksp,Amat,Z,U1);
163: }
164: VecAXPY(U1,-alpha,U);
165: VecNorm(U1,NORM_2,&beta);
166: if (beta > 0.0) {
167: VecScale(U1,1.0/beta); /* beta*U1 = Amat*V - alpha*U */
168: }
170: KSP_MatMultTranspose(ksp,Amat,U1,V1);
171: VecAXPY(V1,-beta,V);
172: if (nopreconditioner) {
173: VecNorm(V1,NORM_2,&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 (SE) {
197: VecCopy(W,W2);
198: VecSquare(W2);
199: VecScale(W2,1.0/(rho*rho));
200: VecAXPY(SE, 1.0, W2); /* SE <- 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 (SE) {
228: tmp = 1.0;
229: MatGetSize(Amat,&size1,&size2);
230: if (size1 > size2) tmp = size1 - size2;
231: tmp = rnorm / PetscSqrtScalar(tmp);
232: VecSqrtAbs(SE);
233: VecScale(SE,tmp);
234: }
235: return(0);
236: }
239: PetscErrorCode KSPDestroy_LSQR(KSP ksp)
240: {
241: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
245: /* Free work vectors */
246: if (lsqr->vwork_n) {
247: VecDestroyVecs(lsqr->nwork_n,&lsqr->vwork_n);
248: }
249: if (lsqr->vwork_m) {
250: VecDestroyVecs(lsqr->nwork_m,&lsqr->vwork_m);
251: }
252: if (lsqr->se_flg) {
253: VecDestroy(&lsqr->se);
254: }
255: PetscFree(ksp->data);
256: return(0);
257: }
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: }
270: PetscErrorCode KSPLSQRGetStandardErrorVec(KSP ksp,Vec *se)
271: {
272: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
275: *se = lsqr->se;
276: return(0);
277: }
279: PetscErrorCode KSPLSQRGetArnorm(KSP ksp,PetscReal *arnorm, PetscReal *rhs_norm, PetscReal *anorm)
280: {
281: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
285: *arnorm = lsqr->arnorm;
286: if (anorm) {
287: if (lsqr->anorm < 0.0) {
288: PC pc;
289: Mat Amat;
290: KSPGetPC(ksp,&pc);
291: PCGetOperators(pc,&Amat,NULL);
292: MatNorm(Amat,NORM_FROBENIUS,&lsqr->anorm);
293: }
294: *anorm = lsqr->anorm;
295: }
296: if (rhs_norm) *rhs_norm = lsqr->rhs_norm;
297: return(0);
298: }
300: /*@C
301: 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
303: Collective on KSP
305: Input Parameters:
306: + ksp - iterative context
307: . n - iteration number
308: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
309: - dummy - viewer and format context
311: Level: intermediate
313: .keywords: KSP, default, monitor, residual
315: .seealso: KSPMonitorSet(), KSPMonitorTrueResidualNorm(), KSPMonitorLGResidualNormCreate(), KSPMonitorDefault()
316: @*/
317: PetscErrorCode KSPLSQRMonitorDefault(KSP ksp,PetscInt n,PetscReal rnorm,PetscViewerAndFormat *dummy)
318: {
320: PetscViewer viewer = dummy->viewer;
321: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
324: PetscViewerPushFormat(viewer,dummy->format);
325: PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
326: if (((PetscObject)ksp)->prefix) {
327: PetscViewerASCIIPrintf(viewer," Residual norm and norm of normal equations for %s solve.\n",((PetscObject)ksp)->prefix);
328: }
329: if (!n) {
330: PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %14.12e\n",n,rnorm);
331: } else {
332: PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %14.12e Residual norm normal equations %14.12e\n",n,rnorm,lsqr->arnorm);
333: }
334: PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
335: PetscViewerPopFormat(viewer);
336: return(0);
337: }
339: PetscErrorCode KSPSetFromOptions_LSQR(PetscOptionItems *PetscOptionsObject,KSP ksp)
340: {
342: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
345: PetscOptionsHead(PetscOptionsObject,"KSP LSQR Options");
346: PetscOptionsName("-ksp_lsqr_set_standard_error","Set Standard Error Estimates of Solution","KSPLSQRSetStandardErrorVec",&lsqr->se_flg);
347: KSPMonitorSetFromOptions(ksp,"-ksp_lsqr_monitor","Monitor residual norm and norm of residual of normal equations","KSPMonitorSet",KSPLSQRMonitorDefault);
348: PetscOptionsTail();
349: return(0);
350: }
352: PetscErrorCode KSPView_LSQR(KSP ksp,PetscViewer viewer)
353: {
354: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
356: PetscBool iascii;
359: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
360: if (iascii) {
361: if (lsqr->se) {
362: PetscReal rnorm;
363: KSPLSQRGetStandardErrorVec(ksp,&lsqr->se);
364: VecNorm(lsqr->se,NORM_2,&rnorm);
365: PetscViewerASCIIPrintf(viewer," Norm of Standard Error %g, Iterations %D\n",(double)rnorm,ksp->its);
366: }
367: }
368: return(0);
369: }
371: /*@C
372: KSPLSQRDefaultConverged - Determines convergence of the LSQR Krylov method. This calls KSPConvergedDefault() and if that does not determine convergence then checks
373: convergence for the least squares problem.
375: Collective on KSP
377: Input Parameters:
378: + ksp - iterative context
379: . n - iteration number
380: . rnorm - 2-norm residual value (may be estimated)
381: - ctx - convergence context which must be created by KSPConvergedDefaultCreate()
383: reason is set to:
384: + positive - if the iteration has converged;
385: . negative - if residual norm exceeds divergence threshold;
386: - 0 - otherwise.
388: Notes:
389: 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.
391: Level: intermediate
393: .keywords: KSP, default, convergence, residual
395: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(),
396: KSPConvergedDefaultSetUIRNorm(), KSPConvergedDefaultSetUMIRNorm(), KSPConvergedDefaultCreate(), KSPConvergedDefaultDestroy(), KSPConvergedDefault()
397: @*/
398: PetscErrorCode KSPLSQRDefaultConverged(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *ctx)
399: {
401: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
404: KSPConvergedDefault(ksp,n,rnorm,reason,ctx);
405: if (!n || *reason) return(0);
406: if (lsqr->arnorm/lsqr->rhs_norm < ksp->rtol) *reason = KSP_CONVERGED_RTOL_NORMAL;
407: if (lsqr->arnorm < ksp->abstol) *reason = KSP_CONVERGED_ATOL_NORMAL;
408: return(0);
409: }
411: /*MC
412: KSPLSQR - This implements LSQR
414: Options Database Keys:
415: + -ksp_lsqr_set_standard_error - Set Standard Error Estimates of Solution see KSPLSQRSetStandardErrorVec()
416: . -ksp_lsqr_monitor - Monitor residual norm and norm of residual of normal equations
417: - see KSPSolve()
419: Level: beginner
421: Notes:
422: Supports non-square (rectangular) matrices.
424: This varient, when applied with no preconditioning is identical to the original algorithm in exact arithematic; however, in practice, with no preconditioning
425: due to inexact arithematic, it can converge differently. Hence when no preconditioner is used (PCType PCNONE) it automatically reverts to the original algorithm.
427: With the PETSc built-in preconditioners, such as ICC, one should call KSPSetOperators(ksp,A,A'*A)) since the preconditioner needs to work
428: for the normal equations A'*A.
430: Supports only left preconditioning.
432: References:
433: . 1. - The original unpreconditioned algorithm can be found in Paige and Saunders, ACM Transactions on Mathematical Software, Vol 8, 1982.
435: In exact arithmetic the LSQR method (with no preconditioning) is identical to the KSPCG algorithm applied to the normal equations.
436: 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
437: track the true norm of the residual and uses that in the convergence test.
439: Developer Notes: How is this related to the KSPCGNE implementation? One difference is that KSPCGNE applies
440: the preconditioner transpose times the preconditioner, so one does not need to pass A'*A as the third argument to KSPSetOperators().
443: 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()
445: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPLSQRDefaultConverged()
447: M*/
448: PETSC_EXTERN PetscErrorCode KSPCreate_LSQR(KSP ksp)
449: {
450: KSP_LSQR *lsqr;
454: PetscNewLog(ksp,&lsqr);
455: lsqr->se = NULL;
456: lsqr->se_flg = PETSC_FALSE;
457: lsqr->arnorm = 0.0;
458: ksp->data = (void*)lsqr;
459: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,3);
461: ksp->ops->setup = KSPSetUp_LSQR;
462: ksp->ops->solve = KSPSolve_LSQR;
463: ksp->ops->destroy = KSPDestroy_LSQR;
464: ksp->ops->buildsolution = KSPBuildSolutionDefault;
465: ksp->ops->buildresidual = KSPBuildResidualDefault;
466: ksp->ops->setfromoptions = KSPSetFromOptions_LSQR;
467: ksp->ops->view = KSPView_LSQR;
468: ksp->converged = KSPLSQRDefaultConverged;
469: return(0);
470: }