Actual source code: lsqr.c
petsc-3.11.4 2019-09-28
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: } 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);
43: if (lsqr->vwork_m) {
44: VecDestroyVecs(lsqr->nwork_m,&lsqr->vwork_m);
45: }
47: if (lsqr->vwork_n) {
48: VecDestroyVecs(lsqr->nwork_n,&lsqr->vwork_n);
49: }
51: lsqr->nwork_m = 2;
52: if (nopreconditioner) lsqr->nwork_n = 4;
53: else lsqr->nwork_n = 5;
54: KSPCreateVecs(ksp,lsqr->nwork_n,&lsqr->vwork_n,lsqr->nwork_m,&lsqr->vwork_m);
56: if (lsqr->se_flg && !lsqr->se) {
57: VecDuplicate(lsqr->vwork_n[0],&lsqr->se);
58: VecSet(lsqr->se,PETSC_INFINITY);
59: } else if (!lsqr->se_flg) {
60: VecDestroy(&lsqr->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,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: /* vectors of length m, where system size is mxn */
84: B = ksp->vec_rhs;
85: U = lsqr->vwork_m[0];
86: U1 = lsqr->vwork_m[1];
88: /* vectors of length n */
89: X = ksp->vec_sol;
90: W = lsqr->vwork_n[0];
91: V = lsqr->vwork_n[1];
92: V1 = lsqr->vwork_n[2];
93: W2 = lsqr->vwork_n[3];
94: if (!nopreconditioner) Z = lsqr->vwork_n[4];
96: /* standard error vector */
97: if (lsqr->se) {
98: VecSet(lsqr->se,0.0);
99: }
101: /* Compute initial residual, temporarily use work vector u */
102: if (!ksp->guess_zero) {
103: KSP_MatMult(ksp,Amat,X,U); /* u <- b - Ax */
104: VecAYPX(U,-1.0,B);
105: } else {
106: VecCopy(B,U); /* u <- b (x is 0) */
107: }
109: /* Test for nothing to do */
110: VecNorm(U,NORM_2,&rnorm);
111: KSPCheckNorm(ksp,rnorm);
112: PetscObjectSAWsTakeAccess((PetscObject)ksp);
113: ksp->its = 0;
114: ksp->rnorm = rnorm;
115: PetscObjectSAWsGrantAccess((PetscObject)ksp);
116: KSPLogResidualHistory(ksp,rnorm);
117: KSPMonitor(ksp,0,rnorm);
118: (*ksp->converged)(ksp,0,rnorm,&ksp->reason,ksp->cnvP);
119: if (ksp->reason) return(0);
121: beta = rnorm;
122: VecScale(U,1.0/beta);
123: KSP_MatMultTranspose(ksp,Amat,U,V);
124: if (nopreconditioner) {
125: VecNorm(V,NORM_2,&alpha);
126: KSPCheckNorm(ksp,rnorm);
127: } else {
128: /* this is an application of the preconditioner for the normal equations; not the operator, see the manual page */
129: PCApply(ksp->pc,V,Z);
130: VecDotRealPart(V,Z,&alpha);
131: if (alpha <= 0.0) {
132: ksp->reason = KSP_DIVERGED_BREAKDOWN;
133: return(0);
134: }
135: alpha = PetscSqrtReal(alpha);
136: VecScale(Z,1.0/alpha);
137: }
138: VecScale(V,1.0/alpha);
140: if (nopreconditioner) {
141: VecCopy(V,W);
142: } else {
143: VecCopy(Z,W);
144: }
146: if (lsqr->exact_norm) {
147: MatNorm(Amat,NORM_FROBENIUS,&lsqr->anorm);
148: } else lsqr->anorm = 0.0;
150: lsqr->arnorm = alpha * beta;
151: phibar = beta;
152: rhobar = alpha;
153: i = 0;
154: do {
155: if (nopreconditioner) {
156: KSP_MatMult(ksp,Amat,V,U1);
157: } else {
158: KSP_MatMult(ksp,Amat,Z,U1);
159: }
160: VecAXPY(U1,-alpha,U);
161: VecNorm(U1,NORM_2,&beta);
162: KSPCheckNorm(ksp,beta);
163: if (beta > 0.0) {
164: VecScale(U1,1.0/beta); /* beta*U1 = Amat*V - alpha*U */
165: if (!lsqr->exact_norm) {
166: lsqr->anorm = PetscSqrtScalar(PetscSqr(lsqr->anorm) + PetscSqr(alpha) + PetscSqr(beta));
167: }
168: }
170: KSP_MatMultTranspose(ksp,Amat,U1,V1);
171: VecAXPY(V1,-beta,V);
172: if (nopreconditioner) {
173: VecNorm(V1,NORM_2,&alpha);
174: KSPCheckNorm(ksp,alpha);
175: } else {
176: PCApply(ksp->pc,V1,Z);
177: VecDotRealPart(V1,Z,&alpha);
178: if (alpha <= 0.0) {
179: ksp->reason = KSP_DIVERGED_BREAKDOWN;
180: break;
181: }
182: alpha = PetscSqrtReal(alpha);
183: VecScale(Z,1.0/alpha);
184: }
185: VecScale(V1,1.0/alpha); /* alpha*V1 = Amat^T*U1 - beta*V */
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 (lsqr->se) {
198: VecCopy(W,W2);
199: VecSquare(W2);
200: VecScale(W2,1.0/(rho*rho));
201: VecAXPY(lsqr->se, 1.0, W2); /* lsqr->se <- lsqr->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: PetscObjectSAWsTakeAccess((PetscObject)ksp);
213: ksp->its++;
214: ksp->rnorm = rnorm;
215: PetscObjectSAWsGrantAccess((PetscObject)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) ksp->reason = KSP_DIVERGED_ITS;
227: /* Finish off the standard error estimates */
228: if (lsqr->se) {
229: tmp = 1.0;
230: MatGetSize(Amat,&size1,&size2);
231: if (size1 > size2) tmp = size1 - size2;
232: tmp = rnorm / PetscSqrtScalar(tmp);
233: VecSqrtAbs(lsqr->se);
234: VecScale(lsqr->se,tmp);
235: }
236: return(0);
237: }
240: PetscErrorCode KSPDestroy_LSQR(KSP ksp)
241: {
242: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
246: /* Free work vectors */
247: if (lsqr->vwork_n) {
248: VecDestroyVecs(lsqr->nwork_n,&lsqr->vwork_n);
249: }
250: if (lsqr->vwork_m) {
251: VecDestroyVecs(lsqr->nwork_m,&lsqr->vwork_m);
252: }
253: VecDestroy(&lsqr->se);
254: PetscFree(ksp->data);
255: return(0);
256: }
258: /*@
259: KSPLSQRSetComputeStandardErrorVec - Compute vector of standard error estimates during KSPSolve_LSQR().
261: Not Collective
263: Input Parameters:
264: + ksp - iterative context
265: - flg - compute the vector of standard estimates or not
267: Developer notes:
268: Vaclav: I'm not sure whether this vector is useful for anything.
270: Level: intermediate
272: .keywords: KSP, KSPLSQR
274: .seealso: KSPSolve(), KSPLSQR, KSPLSQRGetStandardErrorVec()
275: @*/
276: PetscErrorCode KSPLSQRSetComputeStandardErrorVec(KSP ksp, PetscBool flg)
277: {
278: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
281: lsqr->se_flg = flg;
282: return(0);
283: }
285: /*@
286: KSPLSQRSetExactMatNorm - Compute exact matrix norm instead of iteratively refined estimate.
288: Not Collective
290: Input Parameters:
291: + ksp - iterative context
292: - flg - compute exact matrix norm or not
294: Notes:
295: By default, flg=PETSC_FALSE. This is usually preferred to avoid possibly expensive computation of the norm.
296: For flg=PETSC_TRUE, we call MatNorm(Amat,NORM_FROBENIUS,&lsqr->anorm) which will work only for some types of explicitly assembled matrices.
297: This can affect convergence rate as KSPLSQRConvergedDefault() assumes different value of ||A|| used in normal equation stopping criterion.
299: Level: intermediate
301: .keywords: KSP, KSPLSQR
303: .seealso: KSPSolve(), KSPLSQR, KSPLSQRGetNorms(), KSPLSQRConvergedDefault()
304: @*/
305: PetscErrorCode KSPLSQRSetExactMatNorm(KSP ksp, PetscBool flg)
306: {
307: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
310: lsqr->exact_norm = flg;
311: return(0);
312: }
314: /*@
315: KSPLSQRGetStandardErrorVec - Get vector of standard error estimates.
316: Only available if -ksp_lsqr_set_standard_error was set to true
317: or KSPLSQRSetComputeStandardErrorVec(ksp, PETSC_TRUE) was called.
318: Otherwise returns NULL.
320: Not Collective
322: Input Parameters:
323: . ksp - iterative context
325: Output Parameters:
326: . se - vector of standard estimates
328: Options Database Keys:
329: . -ksp_lsqr_set_standard_error - set standard error estimates of solution
331: Developer notes:
332: Vaclav: I'm not sure whether this vector is useful for anything.
334: Level: intermediate
336: .keywords: KSP, KSPLSQR
338: .seealso: KSPSolve(), KSPLSQR, KSPLSQRSetComputeStandardErrorVec()
339: @*/
340: PetscErrorCode KSPLSQRGetStandardErrorVec(KSP ksp,Vec *se)
341: {
342: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
345: *se = lsqr->se;
346: return(0);
347: }
349: /*@
350: KSPLSQRGetNorms - Get norm estimates that LSQR computes internally during KSPSolve().
352: Not Collective
354: Input Parameters:
355: . ksp - iterative context
357: Output Parameters:
358: + arnorm - good estimate of norm((A*inv(Pmat))'*r), where r = A*x - b, used in specific stopping criterion
359: - anorm - poor estimate of norm(A*inv(Pmat),'fro') used in specific stopping criterion
361: Notes:
362: Output parameters are meaningful only after KSPSolve().
363: 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.
364: If -ksp_lsqr_exact_mat_norm is set or KSPLSQRSetExactMatNorm(ksp, PETSC_TRUE) called, then anorm is exact Frobenius norm.
366: Level: intermediate
368: .keywords: KSP, KSPLSQR
370: .seealso: KSPSolve(), KSPLSQR, KSPLSQRSetExactMatNorm()
371: @*/
372: PetscErrorCode KSPLSQRGetNorms(KSP ksp,PetscReal *arnorm, PetscReal *anorm)
373: {
374: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
377: if (arnorm) *arnorm = lsqr->arnorm;
378: if (anorm) *anorm = lsqr->anorm;
379: return(0);
380: }
382: /*@C
383: KSPLSQRMonitorDefault - Print the residual norm at each iteration of the LSQR method,
384: norm of the residual of the normal equations A'*A x = A' b, and estimate of matrix norm ||A||.
386: Collective on KSP
388: Input Parameters:
389: + ksp - iterative context
390: . n - iteration number
391: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
392: - dummy - viewer and format context
394: Level: intermediate
396: .keywords: KSP, KSPLSQR, default, monitor, residual
398: .seealso: KSPLSQR, KSPMonitorSet(), KSPMonitorTrueResidualNorm(), KSPMonitorLGResidualNormCreate(), KSPMonitorDefault()
399: @*/
400: PetscErrorCode KSPLSQRMonitorDefault(KSP ksp,PetscInt n,PetscReal rnorm,PetscViewerAndFormat *dummy)
401: {
403: PetscViewer viewer = dummy->viewer;
404: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
407: PetscViewerPushFormat(viewer,dummy->format);
408: PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
409: if (!n && ((PetscObject)ksp)->prefix) {
410: PetscViewerASCIIPrintf(viewer," Residual norm, norm of normal equations, and matrix norm for %s solve.\n",((PetscObject)ksp)->prefix);
411: }
413: if (!n) {
414: PetscViewerASCIIPrintf(viewer,"%3D KSP resid norm %14.12e\n",n,(double)rnorm);
415: } else {
416: 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);
417: }
419: PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
420: PetscViewerPopFormat(viewer);
421: return(0);
422: }
424: PetscErrorCode KSPSetFromOptions_LSQR(PetscOptionItems *PetscOptionsObject,KSP ksp)
425: {
427: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
430: PetscOptionsHead(PetscOptionsObject,"KSP LSQR Options");
431: PetscOptionsBool("-ksp_lsqr_compute_standard_error","Set Standard Error Estimates of Solution","KSPLSQRSetComputeStandardErrorVec",lsqr->se_flg,&lsqr->se_flg,NULL);
432: PetscOptionsBool("-ksp_lsqr_exact_mat_norm","Compute exact matrix norm instead of iteratively refined estimate","KSPLSQRSetExactMatNorm",lsqr->exact_norm,&lsqr->exact_norm,NULL);
433: KSPMonitorSetFromOptions(ksp,"-ksp_lsqr_monitor","Monitor residual norm and norm of residual of normal equations","KSPMonitorSet",KSPLSQRMonitorDefault);
434: PetscOptionsTail();
435: return(0);
436: }
438: PetscErrorCode KSPView_LSQR(KSP ksp,PetscViewer viewer)
439: {
440: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
442: PetscBool iascii;
445: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
446: if (iascii) {
447: if (lsqr->se) {
448: PetscReal rnorm;
449: VecNorm(lsqr->se,NORM_2,&rnorm);
450: PetscViewerASCIIPrintf(viewer," norm of standard error %g, iterations %d\n",(double)rnorm,ksp->its);
451: } else {
452: PetscViewerASCIIPrintf(viewer," standard error not computed\n");
453: }
454: if (lsqr->exact_norm) {
455: PetscViewerASCIIPrintf(viewer," using exact matrix norm\n");
456: } else {
457: PetscViewerASCIIPrintf(viewer," using inexact matrix norm\n");
458: }
459: }
460: return(0);
461: }
463: /*@C
464: KSPLSQRConvergedDefault - Determines convergence of the LSQR Krylov method.
466: Collective on KSP
468: Input Parameters:
469: + ksp - iterative context
470: . n - iteration number
471: . rnorm - 2-norm residual value (may be estimated)
472: - ctx - convergence context which must be created by KSPConvergedDefaultCreate()
474: reason is set to:
475: + positive - if the iteration has converged;
476: . negative - if residual norm exceeds divergence threshold;
477: - 0 - otherwise.
479: Notes:
480: KSPConvergedDefault() is called first to check for convergence in A*x=b.
481: If that does not determine convergence then checks convergence for the least squares problem, i.e. in min{|b-A*x|}.
482: 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.
483: KSP_CONVERGED_RTOL_NORMAL is returned if ||A'*r|| < rtol * ||A|| * ||r||.
484: Matrix norm ||A|| is iteratively refined estimate, see KSPLSQRGetNorms().
485: This criterion is now largely compatible with that in MATLAB lsqr().
487: Level: intermediate
489: .keywords: KSP, KSPLSQR, default, convergence, residual
491: .seealso: KSPLSQR, KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(),
492: KSPConvergedDefaultSetUIRNorm(), KSPConvergedDefaultSetUMIRNorm(), KSPConvergedDefaultCreate(), KSPConvergedDefaultDestroy(), KSPConvergedDefault(), KSPLSQRGetNorms(), KSPLSQRSetExactMatNorm()
493: @*/
494: PetscErrorCode KSPLSQRConvergedDefault(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *ctx)
495: {
497: KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;
500: /* check for convergence in A*x=b */
501: KSPConvergedDefault(ksp,n,rnorm,reason,ctx);
502: if (!n || *reason) return(0);
504: /* check for convergence in min{|b-A*x|} */
505: if (lsqr->arnorm < ksp->abstol) {
506: 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);
507: *reason = KSP_CONVERGED_ATOL_NORMAL;
508: } else if (lsqr->arnorm < ksp->rtol * lsqr->anorm * rnorm) {
509: 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);
510: *reason = KSP_CONVERGED_RTOL_NORMAL;
511: }
512: return(0);
513: }
515: /*MC
516: KSPLSQR - This implements LSQR
518: Options Database Keys:
519: + -ksp_lsqr_set_standard_error - set standard error estimates of solution, see KSPLSQRSetComputeStandardErrorVec() and KSPLSQRGetStandardErrorVec()
520: . -ksp_lsqr_exact_mat_norm - compute exact matrix norm instead of iteratively refined estimate, see KSPLSQRSetExactMatNorm()
521: - -ksp_lsqr_monitor - monitor residual norm, norm of residual of normal equations A'*A x = A' b, and estimate of matrix norm ||A||
523: Level: beginner
525: Notes:
526: Supports non-square (rectangular) matrices.
528: This varient, when applied with no preconditioning is identical to the original algorithm in exact arithematic; however, in practice, with no preconditioning
529: due to inexact arithematic, it can converge differently. Hence when no preconditioner is used (PCType PCNONE) it automatically reverts to the original algorithm.
531: With the PETSc built-in preconditioners, such as ICC, one should call KSPSetOperators(ksp,A,A'*A)) since the preconditioner needs to work
532: for the normal equations A'*A.
534: Supports only left preconditioning.
536: 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().
538: References:
539: . 1. - The original unpreconditioned algorithm can be found in Paige and Saunders, ACM Transactions on Mathematical Software, Vol 8, 1982.
541: In exact arithmetic the LSQR method (with no preconditioning) is identical to the KSPCG algorithm applied to the normal equations.
542: 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
543: track the true norm of the residual and uses that in the convergence test.
545: Developer Notes:
546: How is this related to the KSPCGNE implementation? One difference is that KSPCGNE applies
547: the preconditioner transpose times the preconditioner, so one does not need to pass A'*A as the third argument to KSPSetOperators().
551: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPSolve(), KSPLSQRConvergedDefault(), KSPLSQRSetComputeStandardErrorVec(), KSPLSQRGetStandardErrorVec(), KSPLSQRSetExactMatNorm()
553: M*/
554: PETSC_EXTERN PetscErrorCode KSPCreate_LSQR(KSP ksp)
555: {
556: KSP_LSQR *lsqr;
557: void *ctx;
561: PetscNewLog(ksp,&lsqr);
562: lsqr->se = NULL;
563: lsqr->se_flg = PETSC_FALSE;
564: lsqr->exact_norm = PETSC_FALSE;
565: lsqr->anorm = -1.0;
566: lsqr->arnorm = -1.0;
567: ksp->data = (void*)lsqr;
568: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,3);
570: ksp->ops->setup = KSPSetUp_LSQR;
571: ksp->ops->solve = KSPSolve_LSQR;
572: ksp->ops->destroy = KSPDestroy_LSQR;
573: ksp->ops->setfromoptions = KSPSetFromOptions_LSQR;
574: ksp->ops->view = KSPView_LSQR;
576: KSPConvergedDefaultCreate(&ctx);
577: KSPSetConvergenceTest(ksp,KSPLSQRConvergedDefault,ctx,KSPConvergedDefaultDestroy);
578: return(0);
579: }