Actual source code: lsqr.c
petsc-3.3-p7 2013-05-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>
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) {
40: lsqr->nwork_n = 4;
41: } else {
42: lsqr->nwork_n = 5;
43: }
44: if (lsqr->vwork_n) {
45: VecDestroyVecs(lsqr->nwork_n,&lsqr->vwork_n);
46: }
47: KSPGetVecs(ksp,lsqr->nwork_n,&lsqr->vwork_n,lsqr->nwork_m,&lsqr->vwork_m);
48: if (lsqr->se_flg && !lsqr->se){
49: /* lsqr->se is not set by user, get it from pmat */
50: Vec *se;
51: KSPGetVecs(ksp,1,&se,0,PETSC_NULL);
52: lsqr->se = *se;
53: PetscFree(se);
54: }
55: return(0);
56: }
60: static PetscErrorCode KSPSolve_LSQR(KSP ksp)
61: {
63: PetscInt i,size1,size2;
64: PetscScalar rho,rhobar,phi,phibar,theta,c,s,tmp,tau,alphac;
65: PetscReal beta,alpha,rnorm;
66: Vec X,B,V,V1,U,U1,TMP,W,W2,SE,Z = PETSC_NULL;
67: Mat Amat,Pmat;
68: MatStructure pflag;
69: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
70: PetscBool diagonalscale,nopreconditioner;
71:
73: PCGetDiagonalScale(ksp->pc,&diagonalscale);
74: if (diagonalscale) SETERRQ1(((PetscObject)ksp)->comm,PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);
76: PCGetOperators(ksp->pc,&Amat,&Pmat,&pflag);
77: PetscObjectTypeCompare((PetscObject)ksp->pc,PCNONE,&nopreconditioner);
79: /* nopreconditioner =PETSC_FALSE; */
80: /* Calculate norm of right hand side */
81: VecNorm(ksp->vec_rhs,NORM_2,&lsqr->rhs_norm);
83: /* mark norm of matrix with negative number to indicate it has not yet been computed */
84: lsqr->anorm = -1.0;
86: /* vectors of length m, where system size is mxn */
87: B = ksp->vec_rhs;
88: U = lsqr->vwork_m[0];
89: U1 = lsqr->vwork_m[1];
91: /* vectors of length n */
92: X = ksp->vec_sol;
93: W = lsqr->vwork_n[0];
94: V = lsqr->vwork_n[1];
95: V1 = lsqr->vwork_n[2];
96: W2 = lsqr->vwork_n[3];
97: if (!nopreconditioner) {
98: Z = lsqr->vwork_n[4];
99: }
101: /* standard error vector */
102: SE = lsqr->se;
103: if (SE){
104: VecGetSize(SE,&size1);
105: VecGetSize(X ,&size2);
106: 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);
107: VecSet(SE,0.0);
108: }
110: /* Compute initial residual, temporarily use work vector u */
111: if (!ksp->guess_zero) {
112: KSP_MatMult(ksp,Amat,X,U); /* u <- b - Ax */
113: VecAYPX(U,-1.0,B);
114: } else {
115: VecCopy(B,U); /* u <- b (x is 0) */
116: }
118: /* Test for nothing to do */
119: VecNorm(U,NORM_2,&rnorm);
120: PetscObjectTakeAccess(ksp);
121: ksp->its = 0;
122: ksp->rnorm = rnorm;
123: PetscObjectGrantAccess(ksp);
124: KSPLogResidualHistory(ksp,rnorm);
125: KSPMonitor(ksp,0,rnorm);
126: (*ksp->converged)(ksp,0,rnorm,&ksp->reason,ksp->cnvP);
127: if (ksp->reason) return(0);
129: beta = rnorm;
130: VecScale(U,1.0/beta);
131: KSP_MatMultTranspose(ksp,Amat,U,V);
132: if (nopreconditioner) {
133: VecNorm(V,NORM_2,&alpha);
134: } else {
135: PCApply(ksp->pc,V,Z);
136: VecDot(V,Z,&alphac);
137: if (PetscRealPart(alphac) <= 0.0) {
138: ksp->reason = KSP_DIVERGED_BREAKDOWN;
139: return(0);
140: }
141: alpha = PetscSqrtReal(PetscRealPart(alphac));
142: VecScale(Z,1.0/alpha);
143: }
144: VecScale(V,1.0/alpha);
146: if (nopreconditioner){
147: VecCopy(V,W);
148: } else {
149: VecCopy(Z,W);
150: }
152: lsqr->arnorm = alpha * beta;
153: phibar = beta;
154: rhobar = alpha;
155: tau = -beta;
156: i = 0;
157: do {
158: if (nopreconditioner) {
159: KSP_MatMult(ksp,Amat,V,U1);
160: } else {
161: KSP_MatMult(ksp,Amat,Z,U1);
162: }
163: VecAXPY(U1,-alpha,U);
164: VecNorm(U1,NORM_2,&beta);
165: if (beta == 0.0){
166: ksp->reason = KSP_DIVERGED_BREAKDOWN;
167: break;
168: }
169: VecScale(U1,1.0/beta);
171: KSP_MatMultTranspose(ksp,Amat,U1,V1);
172: VecAXPY(V1,-beta,V);
173: if (nopreconditioner) {
174: VecNorm(V1,NORM_2,&alpha);
175: } else {
176: PCApply(ksp->pc,V1,Z);
177: VecDot(V1,Z,&alphac);
178: if (PetscRealPart(alphac) <= 0.0) {
179: ksp->reason = KSP_DIVERGED_BREAKDOWN;
180: break;
181: }
182: alpha = PetscSqrtReal(PetscRealPart(alphac));
183: VecScale(Z,1.0/alpha);
184: }
185: VecScale(V1,1.0/alpha);
186: rho = PetscSqrtScalar(rhobar*rhobar + beta*beta);
187: c = rhobar / rho;
188: s = beta / rho;
189: theta = s * alpha;
190: rhobar = - c * alpha;
191: phi = c * phibar;
192: phibar = s * phibar;
193: tau = s * phi;
195: VecAXPY(X,phi/rho,W); /* x <- x + (phi/rho) w */
197: if (SE) {
198: VecCopy(W,W2);
199: VecSquare(W2);
200: VecScale(W2,1.0/(rho*rho));
201: VecAXPY(SE, 1.0, W2); /* SE <- SE + (w^2/rho^2) */
202: }
203: if (nopreconditioner) {
204: VecAYPX(W,-theta/rho,V1); /* w <- v - (theta/rho) w */
205: } else {
206: VecAYPX(W,-theta/rho,Z); /* w <- z - (theta/rho) w */
207: }
209: lsqr->arnorm = alpha*PetscAbsScalar(tau);
210: rnorm = PetscRealPart(phibar);
212: PetscObjectTakeAccess(ksp);
213: ksp->its++;
214: ksp->rnorm = rnorm;
215: PetscObjectGrantAccess(ksp);
216: KSPLogResidualHistory(ksp,rnorm);
217: KSPMonitor(ksp,i+1,rnorm);
218: (*ksp->converged)(ksp,i+1,rnorm,&ksp->reason,ksp->cnvP);
219: if (ksp->reason) break;
220: SWAP(U1,U,TMP);
221: SWAP(V1,V,TMP);
223: i++;
224: } while (i<ksp->max_it);
225: if (i >= ksp->max_it && !ksp->reason) {
226: ksp->reason = KSP_DIVERGED_ITS;
227: }
229: /* Finish off the standard error estimates */
230: if (SE) {
231: tmp = 1.0;
232: MatGetSize(Amat,&size1,&size2);
233: if ( size1 > size2 ) tmp = size1 - size2;
234: tmp = rnorm / PetscSqrtScalar(tmp);
235: VecSqrtAbs(SE);
236: VecScale(SE,tmp);
237: }
239: return(0);
240: }
245: PetscErrorCode KSPDestroy_LSQR(KSP ksp)
246: {
247: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
252: /* Free work vectors */
253: if (lsqr->vwork_n) {
254: VecDestroyVecs(lsqr->nwork_n,&lsqr->vwork_n);
255: }
256: if (lsqr->vwork_m) {
257: VecDestroyVecs(lsqr->nwork_m,&lsqr->vwork_m);
258: }
259: if (lsqr->se_flg){
260: VecDestroy(&lsqr->se);
261: }
262: PetscFree(ksp->data);
263: return(0);
264: }
268: PetscErrorCode KSPLSQRSetStandardErrorVec( KSP ksp, Vec se )
269: {
270: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
274: VecDestroy(&lsqr->se);
275: lsqr->se = se;
276: return(0);
277: }
281: PetscErrorCode KSPLSQRGetStandardErrorVec( KSP ksp,Vec *se )
282: {
283: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
286: *se = lsqr->se;
287: return(0);
288: }
292: PetscErrorCode KSPLSQRGetArnorm( KSP ksp,PetscReal *arnorm, PetscReal *rhs_norm , PetscReal *anorm)
293: {
294: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
298: *arnorm = lsqr->arnorm;
299: if (anorm) {
300: if (lsqr->anorm < 0.0) {
301: PC pc;
302: Mat Amat;
303: KSPGetPC(ksp,&pc);
304: PCGetOperators(pc,&Amat,PETSC_NULL,PETSC_NULL);
305: MatNorm(Amat,NORM_FROBENIUS,&lsqr->anorm);
306: }
307: *anorm = lsqr->anorm;
308: }
309: if (rhs_norm) {
310: *rhs_norm = lsqr->rhs_norm;
311: }
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 - unused monitor context
328: Level: intermediate
330: .keywords: KSP, default, monitor, residual
332: .seealso: KSPMonitorSet(), KSPMonitorTrueResidualNorm(), KSPMonitorLGCreate(), KSPMonitorDefault()
333: @*/
334: PetscErrorCode KSPLSQRMonitorDefault(KSP ksp,PetscInt n,PetscReal rnorm,void *dummy)
335: {
336: PetscErrorCode ierr;
337: PetscViewer viewer = dummy ? (PetscViewer) dummy : PETSC_VIEWER_STDOUT_(((PetscObject)ksp)->comm);
338: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
341: PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
342: if (((PetscObject)ksp)->prefix) {
343: PetscViewerASCIIPrintf(viewer," Residual norm and norm of normal equations for %s solve.\n",((PetscObject)ksp)->prefix);
344: }
345: if (!n) {
346: PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %14.12e\n",n,rnorm);
347: } else {
348: PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %14.12e Residual norm normal equations %14.12e\n",n,rnorm,lsqr->arnorm);
349: }
350: PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
351: return(0);
352: }
356: PetscErrorCode KSPSetFromOptions_LSQR(KSP ksp)
357: {
359: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
360: char monfilename[PETSC_MAX_PATH_LEN];
361: PetscViewer monviewer;
362: PetscBool flg;
365: PetscOptionsHead("KSP LSQR Options");
366: PetscOptionsName("-ksp_lsqr_set_standard_error","Set Standard Error Estimates of Solution","KSPLSQRSetStandardErrorVec",&lsqr->se_flg);
367: PetscOptionsString("-ksp_lsqr_monitor","Monitor residual norm and norm of residual of normal equations","KSPMonitorSet","stdout",monfilename,PETSC_MAX_PATH_LEN,&flg);
368: if (flg) {
369: PetscViewerASCIIOpen(((PetscObject)ksp)->comm,monfilename,&monviewer);
370: KSPMonitorSet(ksp,KSPLSQRMonitorDefault,monviewer,(PetscErrorCode (*)(void**))PetscViewerDestroy);
371: }
372: PetscOptionsTail();
373: return(0);
374: }
378: PetscErrorCode KSPView_LSQR(KSP ksp,PetscViewer viewer)
379: {
380: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
382: PetscBool iascii;
385: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
386: if (iascii) {
387: if (lsqr->se) {
388: PetscReal rnorm;
389: KSPLSQRGetStandardErrorVec(ksp,&lsqr->se);
390: VecNorm(lsqr->se,NORM_2,&rnorm);
391: PetscViewerASCIIPrintf(viewer," Norm of Standard Error %G, Iterations %D\n",rnorm,ksp->its);
392: }
393: } else SETERRQ1(((PetscObject)ksp)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for KSP LSQR",((PetscObject)viewer)->type_name);
394: return(0);
395: }
399: /*@C
400: KSPLSQRDefaultConverged - Determines convergence of the LSQR Krylov method. This calls KSPDefaultConverged() and if that does not determine convergence then checks
401: convergence for the least squares problem.
403: Collective on KSP
405: Input Parameters:
406: + ksp - iterative context
407: . n - iteration number
408: . rnorm - 2-norm residual value (may be estimated)
409: - ctx - convergence context which must be created by KSPDefaultConvergedCreate()
411: reason is set to:
412: + positive - if the iteration has converged;
413: . negative - if residual norm exceeds divergence threshold;
414: - 0 - otherwise.
416: Notes:
417: 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.
419: Level: intermediate
421: .keywords: KSP, default, convergence, residual
423: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPSkipConverged(), KSPConvergedReason, KSPGetConvergedReason(),
424: KSPDefaultConvergedSetUIRNorm(), KSPDefaultConvergedSetUMIRNorm(), KSPDefaultConvergedCreate(), KSPDefaultConvergedDestroy(), KSPDefaultConverged()
425: @*/
426: PetscErrorCode KSPLSQRDefaultConverged(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *ctx)
427: {
428: PetscErrorCode ierr;
429: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
432: KSPDefaultConverged(ksp,n,rnorm,reason,ctx);
433: if (!n || *reason) return(0);
434: if (lsqr->arnorm/lsqr->rhs_norm < ksp->rtol) {
435: *reason = KSP_CONVERGED_RTOL_NORMAL;
436: }
437: if (lsqr->arnorm < ksp->abstol) {
438: *reason = KSP_CONVERGED_ATOL_NORMAL;
439: }
440: return(0);
441: }
445: /*MC
446: KSPLSQR - This implements LSQR
448: Options Database Keys:
449: + -ksp_lsqr_set_standard_error - Set Standard Error Estimates of Solution see KSPLSQRSetStandardErrorVec()
450: . -ksp_lsqr_monitor - Monitor residual norm and norm of residual of normal equations
451: - see KSPSolve()
453: Level: beginner
455: Notes:
456: This varient, when applied with no preconditioning is identical to the original algorithm in exact arithematic; however, in practice, with no preconditioning
457: due to inexact arithematic, it can converge differently. Hence when no preconditioner is used (PCType PCNONE) it automatically reverts to the original algorithm.
459: With the PETSc built-in preconditioners, such as ICC, one should call KSPSetOperators(ksp,A,A'*A,...) since the preconditioner needs to work
460: for the normal equations A'*A.
462: Supports only left preconditioning.
464: References:The original unpreconditioned algorithm can be found in Paige and Saunders, ACM Transactions on Mathematical Software, Vol 8, pp 43-71, 1982.
465: In exact arithmetic the LSQR method (with no preconditioning) is identical to the KSPCG algorithm applied to the normal equations.
466: 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
467: track the true norm of the residual and uses that in the convergence test.
469: Developer Notes: How is this related to the KSPCGNE implementation? One difference is that KSPCGNE applies
470: the preconditioner transpose times the preconditioner, so one does not need to pass A'*A as the third argument to KSPSetOperators().
471:
473: 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()
475: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPLSQRDefaultConverged()
477: M*/
478: EXTERN_C_BEGIN
481: PetscErrorCode KSPCreate_LSQR(KSP ksp)
482: {
483: KSP_LSQR *lsqr;
487: PetscNewLog(ksp,KSP_LSQR,&lsqr);
488: lsqr->se = PETSC_NULL;
489: lsqr->se_flg = PETSC_FALSE;
490: lsqr->arnorm = 0.0;
491: ksp->data = (void*)lsqr;
492: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);
494: ksp->ops->setup = KSPSetUp_LSQR;
495: ksp->ops->solve = KSPSolve_LSQR;
496: ksp->ops->destroy = KSPDestroy_LSQR;
497: ksp->ops->buildsolution = KSPDefaultBuildSolution;
498: ksp->ops->buildresidual = KSPDefaultBuildResidual;
499: ksp->ops->setfromoptions = KSPSetFromOptions_LSQR;
500: ksp->ops->view = KSPView_LSQR;
501: ksp->converged = KSPLSQRDefaultConverged;
502: return(0);
503: }
504: EXTERN_C_END
508: PetscErrorCode VecSquare(Vec v)
509: {
511: PetscScalar *x;
512: PetscInt i, n;
515: VecGetLocalSize(v, &n);
516: VecGetArray(v, &x);
517: for(i = 0; i < n; i++) {
518: x[i] *= x[i];
519: }
520: VecRestoreArray(v, &x);
521: return(0);
522: }