Actual source code: lsqr.c
petsc-master 2020-07-04
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: PetscBool exact_norm; /* flag for -ksp_lsqr_exact_mat_norm */
16: PetscReal arnorm; /* Good estimate of norm((A*inv(Pmat))'*r), where r = A*x - b, used in specific stopping criterion */
17: PetscReal anorm; /* Poor estimate of norm(A*inv(Pmat),'fro') used in specific stopping criterion */
18: /* Backup previous convergence test */
19: PetscErrorCode (*converged)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*);
20: PetscErrorCode (*convergeddestroy)(void*);
21: void *cnvP;
22: } KSP_LSQR;
24: static PetscErrorCode VecSquare(Vec v)
25: {
27: PetscScalar *x;
28: PetscInt i, n;
31: VecGetLocalSize(v, &n);
32: VecGetArray(v, &x);
33: for (i = 0; i < n; i++) x[i] *= PetscConj(x[i]);
34: VecRestoreArray(v, &x);
35: return(0);
36: }
38: static PetscErrorCode KSPSetUp_LSQR(KSP ksp)
39: {
41: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
42: PetscBool nopreconditioner;
45: PetscObjectTypeCompare((PetscObject)ksp->pc,PCNONE,&nopreconditioner);
47: if (lsqr->vwork_m) {
48: VecDestroyVecs(lsqr->nwork_m,&lsqr->vwork_m);
49: }
51: if (lsqr->vwork_n) {
52: VecDestroyVecs(lsqr->nwork_n,&lsqr->vwork_n);
53: }
55: lsqr->nwork_m = 2;
56: if (nopreconditioner) lsqr->nwork_n = 4;
57: else lsqr->nwork_n = 5;
58: KSPCreateVecs(ksp,lsqr->nwork_n,&lsqr->vwork_n,lsqr->nwork_m,&lsqr->vwork_m);
60: if (lsqr->se_flg && !lsqr->se) {
61: VecDuplicate(lsqr->vwork_n[0],&lsqr->se);
62: VecSet(lsqr->se,PETSC_INFINITY);
63: } else if (!lsqr->se_flg) {
64: VecDestroy(&lsqr->se);
65: }
66: return(0);
67: }
69: static PetscErrorCode KSPSolve_LSQR(KSP ksp)
70: {
72: PetscInt i,size1,size2;
73: PetscScalar rho,rhobar,phi,phibar,theta,c,s,tmp,tau;
74: PetscReal beta,alpha,rnorm;
75: Vec X,B,V,V1,U,U1,TMP,W,W2,Z = NULL;
76: Mat Amat,Pmat;
77: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
78: PetscBool diagonalscale,nopreconditioner;
81: PCGetDiagonalScale(ksp->pc,&diagonalscale);
82: if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);
84: PCGetOperators(ksp->pc,&Amat,&Pmat);
85: PetscObjectTypeCompare((PetscObject)ksp->pc,PCNONE,&nopreconditioner);
87: /* vectors of length m, where system size is mxn */
88: B = ksp->vec_rhs;
89: U = lsqr->vwork_m[0];
90: U1 = lsqr->vwork_m[1];
92: /* vectors of length n */
93: X = ksp->vec_sol;
94: W = lsqr->vwork_n[0];
95: V = lsqr->vwork_n[1];
96: V1 = lsqr->vwork_n[2];
97: W2 = lsqr->vwork_n[3];
98: if (!nopreconditioner) Z = lsqr->vwork_n[4];
100: /* standard error vector */
101: if (lsqr->se) {
102: VecSet(lsqr->se,0.0);
103: }
105: /* Compute initial residual, temporarily use work vector u */
106: if (!ksp->guess_zero) {
107: KSP_MatMult(ksp,Amat,X,U); /* u <- b - Ax */
108: VecAYPX(U,-1.0,B);
109: } else {
110: VecCopy(B,U); /* u <- b (x is 0) */
111: }
113: /* Test for nothing to do */
114: VecNorm(U,NORM_2,&rnorm);
115: KSPCheckNorm(ksp,rnorm);
116: PetscObjectSAWsTakeAccess((PetscObject)ksp);
117: ksp->its = 0;
118: ksp->rnorm = rnorm;
119: PetscObjectSAWsGrantAccess((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: KSPCheckNorm(ksp,rnorm);
131: } else {
132: /* this is an application of the preconditioner for the normal equations; not the operator, see the manual page */
133: PCApply(ksp->pc,V,Z);
134: VecDotRealPart(V,Z,&alpha);
135: if (alpha <= 0.0) {
136: ksp->reason = KSP_DIVERGED_BREAKDOWN;
137: return(0);
138: }
139: alpha = PetscSqrtReal(alpha);
140: VecScale(Z,1.0/alpha);
141: }
142: VecScale(V,1.0/alpha);
144: if (nopreconditioner) {
145: VecCopy(V,W);
146: } else {
147: VecCopy(Z,W);
148: }
150: if (lsqr->exact_norm) {
151: MatNorm(Amat,NORM_FROBENIUS,&lsqr->anorm);
152: } else lsqr->anorm = 0.0;
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: KSPCheckNorm(ksp,beta);
167: if (beta > 0.0) {
168: VecScale(U1,1.0/beta); /* beta*U1 = Amat*V - alpha*U */
169: if (!lsqr->exact_norm) {
170: lsqr->anorm = PetscSqrtScalar(PetscSqr(lsqr->anorm) + PetscSqr(alpha) + PetscSqr(beta));
171: }
172: }
174: KSP_MatMultTranspose(ksp,Amat,U1,V1);
175: VecAXPY(V1,-beta,V);
176: if (nopreconditioner) {
177: VecNorm(V1,NORM_2,&alpha);
178: KSPCheckNorm(ksp,alpha);
179: } else {
180: PCApply(ksp->pc,V1,Z);
181: VecDotRealPart(V1,Z,&alpha);
182: if (alpha <= 0.0) {
183: ksp->reason = KSP_DIVERGED_BREAKDOWN;
184: break;
185: }
186: alpha = PetscSqrtReal(alpha);
187: VecScale(Z,1.0/alpha);
188: }
189: VecScale(V1,1.0/alpha); /* alpha*V1 = Amat^T*U1 - beta*V */
190: rho = PetscSqrtScalar(rhobar*rhobar + beta*beta);
191: c = rhobar / rho;
192: s = beta / rho;
193: theta = s * alpha;
194: rhobar = -c * alpha;
195: phi = c * phibar;
196: phibar = s * phibar;
197: tau = s * phi;
199: VecAXPY(X,phi/rho,W); /* x <- x + (phi/rho) w */
201: if (lsqr->se) {
202: VecCopy(W,W2);
203: VecSquare(W2);
204: VecScale(W2,1.0/(rho*rho));
205: VecAXPY(lsqr->se, 1.0, W2); /* lsqr->se <- lsqr->se + (w^2/rho^2) */
206: }
207: if (nopreconditioner) {
208: VecAYPX(W,-theta/rho,V1); /* w <- v - (theta/rho) w */
209: } else {
210: VecAYPX(W,-theta/rho,Z); /* w <- z - (theta/rho) w */
211: }
213: lsqr->arnorm = alpha*PetscAbsScalar(tau);
214: rnorm = PetscRealPart(phibar);
216: PetscObjectSAWsTakeAccess((PetscObject)ksp);
217: ksp->its++;
218: ksp->rnorm = rnorm;
219: PetscObjectSAWsGrantAccess((PetscObject)ksp);
220: KSPLogResidualHistory(ksp,rnorm);
221: KSPMonitor(ksp,i+1,rnorm);
222: (*ksp->converged)(ksp,i+1,rnorm,&ksp->reason,ksp->cnvP);
223: if (ksp->reason) break;
224: SWAP(U1,U,TMP);
225: SWAP(V1,V,TMP);
227: i++;
228: } while (i<ksp->max_it);
229: if (i >= ksp->max_it && !ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
231: /* Finish off the standard error estimates */
232: if (lsqr->se) {
233: tmp = 1.0;
234: MatGetSize(Amat,&size1,&size2);
235: if (size1 > size2) tmp = size1 - size2;
236: tmp = rnorm / PetscSqrtScalar(tmp);
237: VecSqrtAbs(lsqr->se);
238: VecScale(lsqr->se,tmp);
239: }
240: return(0);
241: }
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: return(0);
263: }
265: /*@
266: KSPLSQRSetComputeStandardErrorVec - Compute vector of standard error estimates during KSPSolve_LSQR().
268: Not Collective
270: Input Parameters:
271: + ksp - iterative context
272: - flg - compute the vector of standard estimates or not
274: Developer notes:
275: Vaclav: I'm not sure whether this vector is useful for anything.
277: Level: intermediate
279: .seealso: KSPSolve(), KSPLSQR, KSPLSQRGetStandardErrorVec()
280: @*/
281: PetscErrorCode KSPLSQRSetComputeStandardErrorVec(KSP ksp, PetscBool flg)
282: {
283: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
286: lsqr->se_flg = flg;
287: return(0);
288: }
290: /*@
291: KSPLSQRSetExactMatNorm - Compute exact matrix norm instead of iteratively refined estimate.
293: Not Collective
295: Input Parameters:
296: + ksp - iterative context
297: - flg - compute exact matrix norm or not
299: Notes:
300: By default, flg=PETSC_FALSE. This is usually preferred to avoid possibly expensive computation of the norm.
301: For flg=PETSC_TRUE, we call MatNorm(Amat,NORM_FROBENIUS,&lsqr->anorm) which will work only for some types of explicitly assembled matrices.
302: This can affect convergence rate as KSPLSQRConvergedDefault() assumes different value of ||A|| used in normal equation stopping criterion.
304: Level: intermediate
306: .seealso: KSPSolve(), KSPLSQR, KSPLSQRGetNorms(), KSPLSQRConvergedDefault()
307: @*/
308: PetscErrorCode KSPLSQRSetExactMatNorm(KSP ksp, PetscBool flg)
309: {
310: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
313: lsqr->exact_norm = flg;
314: return(0);
315: }
317: /*@
318: KSPLSQRGetStandardErrorVec - Get vector of standard error estimates.
319: Only available if -ksp_lsqr_set_standard_error was set to true
320: or KSPLSQRSetComputeStandardErrorVec(ksp, PETSC_TRUE) was called.
321: Otherwise returns NULL.
323: Not Collective
325: Input Parameters:
326: . ksp - iterative context
328: Output Parameters:
329: . se - vector of standard estimates
331: Options Database Keys:
332: . -ksp_lsqr_set_standard_error - set standard error estimates of solution
334: Developer notes:
335: Vaclav: I'm not sure whether this vector is useful for anything.
337: Level: intermediate
339: .seealso: KSPSolve(), KSPLSQR, KSPLSQRSetComputeStandardErrorVec()
340: @*/
341: PetscErrorCode KSPLSQRGetStandardErrorVec(KSP ksp,Vec *se)
342: {
343: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
346: *se = lsqr->se;
347: return(0);
348: }
350: /*@
351: KSPLSQRGetNorms - Get norm estimates that LSQR computes internally during KSPSolve().
353: Not Collective
355: Input Parameters:
356: . ksp - iterative context
358: Output Parameters:
359: + arnorm - good estimate of norm((A*inv(Pmat))'*r), where r = A*x - b, used in specific stopping criterion
360: - anorm - poor estimate of norm(A*inv(Pmat),'fro') used in specific stopping criterion
362: Notes:
363: Output parameters are meaningful only after KSPSolve().
364: 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.
365: If -ksp_lsqr_exact_mat_norm is set or KSPLSQRSetExactMatNorm(ksp, PETSC_TRUE) called, then anorm is exact Frobenius norm.
367: Level: intermediate
369: .seealso: KSPSolve(), KSPLSQR, KSPLSQRSetExactMatNorm()
370: @*/
371: PetscErrorCode KSPLSQRGetNorms(KSP ksp,PetscReal *arnorm, PetscReal *anorm)
372: {
373: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
376: if (arnorm) *arnorm = lsqr->arnorm;
377: if (anorm) *anorm = lsqr->anorm;
378: return(0);
379: }
381: /*@C
382: KSPLSQRMonitorDefault - Print the residual norm at each iteration of the LSQR method,
383: norm of the residual of the normal equations A'*A x = A' b, and estimate of matrix norm ||A||.
385: Collective on ksp
387: Input Parameters:
388: + ksp - iterative context
389: . n - iteration number
390: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
391: - dummy - viewer and format context
393: Level: intermediate
395: .seealso: KSPLSQR, KSPMonitorSet(), KSPMonitorTrueResidualNorm(), KSPMonitorLGResidualNormCreate(), KSPMonitorDefault()
396: @*/
397: PetscErrorCode KSPLSQRMonitorDefault(KSP ksp,PetscInt n,PetscReal rnorm,PetscViewerAndFormat *dummy)
398: {
400: PetscViewer viewer = dummy->viewer;
401: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
404: PetscViewerPushFormat(viewer,dummy->format);
405: PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
406: if (!n && ((PetscObject)ksp)->prefix) {
407: PetscViewerASCIIPrintf(viewer," Residual norm, norm of normal equations, and matrix norm for %s solve.\n",((PetscObject)ksp)->prefix);
408: }
410: if (!n) {
411: PetscViewerASCIIPrintf(viewer,"%3D KSP resid norm %14.12e\n",n,(double)rnorm);
412: } else {
413: 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);
414: }
416: PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
417: PetscViewerPopFormat(viewer);
418: return(0);
419: }
421: PetscErrorCode KSPSetFromOptions_LSQR(PetscOptionItems *PetscOptionsObject,KSP ksp)
422: {
424: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
427: PetscOptionsHead(PetscOptionsObject,"KSP LSQR Options");
428: PetscOptionsBool("-ksp_lsqr_compute_standard_error","Set Standard Error Estimates of Solution","KSPLSQRSetComputeStandardErrorVec",lsqr->se_flg,&lsqr->se_flg,NULL);
429: PetscOptionsBool("-ksp_lsqr_exact_mat_norm","Compute exact matrix norm instead of iteratively refined estimate","KSPLSQRSetExactMatNorm",lsqr->exact_norm,&lsqr->exact_norm,NULL);
430: KSPMonitorSetFromOptions(ksp,"-ksp_lsqr_monitor","Monitor residual norm and norm of residual of normal equations","KSPMonitorSet",KSPLSQRMonitorDefault);
431: PetscOptionsTail();
432: return(0);
433: }
435: PetscErrorCode KSPView_LSQR(KSP ksp,PetscViewer viewer)
436: {
437: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
439: PetscBool iascii;
442: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
443: if (iascii) {
444: if (lsqr->se) {
445: PetscReal rnorm;
446: VecNorm(lsqr->se,NORM_2,&rnorm);
447: PetscViewerASCIIPrintf(viewer," norm of standard error %g, iterations %d\n",(double)rnorm,ksp->its);
448: } else {
449: PetscViewerASCIIPrintf(viewer," standard error not computed\n");
450: }
451: if (lsqr->exact_norm) {
452: PetscViewerASCIIPrintf(viewer," using exact matrix norm\n");
453: } else {
454: PetscViewerASCIIPrintf(viewer," using inexact matrix norm\n");
455: }
456: }
457: return(0);
458: }
460: /*@C
461: KSPLSQRConvergedDefault - Determines convergence of the LSQR Krylov method.
463: Collective on ksp
465: Input Parameters:
466: + ksp - iterative context
467: . n - iteration number
468: . rnorm - 2-norm residual value (may be estimated)
469: - ctx - convergence context which must be created by KSPConvergedDefaultCreate()
471: reason is set to:
472: + positive - if the iteration has converged;
473: . negative - if residual norm exceeds divergence threshold;
474: - 0 - otherwise.
476: Notes:
477: KSPConvergedDefault() is called first to check for convergence in A*x=b.
478: If that does not determine convergence then checks convergence for the least squares problem, i.e. in min{|b-A*x|}.
479: 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.
480: KSP_CONVERGED_RTOL_NORMAL is returned if ||A'*r|| < rtol * ||A|| * ||r||.
481: Matrix norm ||A|| is iteratively refined estimate, see KSPLSQRGetNorms().
482: This criterion is now largely compatible with that in MATLAB lsqr().
484: Level: intermediate
486: .seealso: KSPLSQR, KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(),
487: KSPConvergedDefaultSetUIRNorm(), KSPConvergedDefaultSetUMIRNorm(), KSPConvergedDefaultCreate(), KSPConvergedDefaultDestroy(), KSPConvergedDefault(), KSPLSQRGetNorms(), KSPLSQRSetExactMatNorm()
488: @*/
489: PetscErrorCode KSPLSQRConvergedDefault(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *ctx)
490: {
492: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
495: /* check for convergence in A*x=b */
496: KSPConvergedDefault(ksp,n,rnorm,reason,ctx);
497: if (!n || *reason) return(0);
499: /* check for convergence in min{|b-A*x|} */
500: if (lsqr->arnorm < ksp->abstol) {
501: 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);
502: *reason = KSP_CONVERGED_ATOL_NORMAL;
503: } else if (lsqr->arnorm < ksp->rtol * lsqr->anorm * rnorm) {
504: 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);
505: *reason = KSP_CONVERGED_RTOL_NORMAL;
506: }
507: return(0);
508: }
510: /*MC
511: KSPLSQR - This implements LSQR
513: Options Database Keys:
514: + -ksp_lsqr_set_standard_error - set standard error estimates of solution, see KSPLSQRSetComputeStandardErrorVec() and KSPLSQRGetStandardErrorVec()
515: . -ksp_lsqr_exact_mat_norm - compute exact matrix norm instead of iteratively refined estimate, see KSPLSQRSetExactMatNorm()
516: - -ksp_lsqr_monitor - monitor residual norm, norm of residual of normal equations A'*A x = A' b, and estimate of matrix norm ||A||
518: Level: beginner
520: Notes:
521: Supports non-square (rectangular) matrices.
523: This varient, when applied with no preconditioning is identical to the original algorithm in exact arithematic; however, in practice, with no preconditioning
524: due to inexact arithematic, it can converge differently. Hence when no preconditioner is used (PCType PCNONE) it automatically reverts to the original algorithm.
526: With the PETSc built-in preconditioners, such as ICC, one should call KSPSetOperators(ksp,A,A'*A)) since the preconditioner needs to work
527: for the normal equations A'*A.
529: Supports only left preconditioning.
531: For least squares problems wit nonzero residual A*x - b, there are additional convergence tests for the residual of the normal equations, A'*(b - Ax), see KSPLSQRConvergedDefault().
533: References:
534: . 1. - The original unpreconditioned algorithm can be found in Paige and Saunders, ACM Transactions on Mathematical Software, Vol 8, 1982.
536: In exact arithmetic the LSQR method (with no preconditioning) is identical to the KSPCG algorithm applied to the normal equations.
537: 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
538: track the true norm of the residual and uses that in the convergence test.
540: Developer Notes:
541: How is this related to the KSPCGNE implementation? One difference is that KSPCGNE applies
542: the preconditioner transpose times the preconditioner, so one does not need to pass A'*A as the third argument to KSPSetOperators().
546: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPSolve(), KSPLSQRConvergedDefault(), KSPLSQRSetComputeStandardErrorVec(), KSPLSQRGetStandardErrorVec(), KSPLSQRSetExactMatNorm()
548: M*/
549: PETSC_EXTERN PetscErrorCode KSPCreate_LSQR(KSP ksp)
550: {
551: KSP_LSQR *lsqr;
552: void *ctx;
556: PetscNewLog(ksp,&lsqr);
557: lsqr->se = NULL;
558: lsqr->se_flg = PETSC_FALSE;
559: lsqr->exact_norm = PETSC_FALSE;
560: lsqr->anorm = -1.0;
561: lsqr->arnorm = -1.0;
562: ksp->data = (void*)lsqr;
563: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,3);
565: ksp->ops->setup = KSPSetUp_LSQR;
566: ksp->ops->solve = KSPSolve_LSQR;
567: ksp->ops->destroy = KSPDestroy_LSQR;
568: ksp->ops->setfromoptions = KSPSetFromOptions_LSQR;
569: ksp->ops->view = KSPView_LSQR;
571: /* Backup current convergence test; remove destroy routine from KSP to prevent destroying the convergence context in KSPSetConvergenceTest() */
572: KSPGetAndClearConvergenceTest(ksp,&lsqr->converged,&lsqr->cnvP,&lsqr->convergeddestroy);
573: /* Override current convergence test */
574: KSPConvergedDefaultCreate(&ctx);
575: KSPSetConvergenceTest(ksp,KSPLSQRConvergedDefault,ctx,KSPConvergedDefaultDestroy);
576: return(0);
577: }