Actual source code: lsqr.c

petsc-3.14.6 2021-03-30
Report Typos and Errors

  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: }