Actual source code: lsqr.c


  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 <petscdraw.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:   PetscBool exact_norm; /* flag for -ksp_lsqr_exact_mat_norm */
 17:   PetscReal arnorm;     /* Good estimate of norm((A*inv(Pmat))'*r), where r = A*x - b, used in specific stopping criterion */
 18:   PetscReal anorm;      /* Poor estimate of norm(A*inv(Pmat),'fro') used in specific stopping criterion */
 19:   /* Backup previous convergence test */
 20:   PetscErrorCode        (*converged)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*);
 21:   PetscErrorCode        (*convergeddestroy)(void*);
 22:   void                  *cnvP;
 23: } KSP_LSQR;

 25: static PetscErrorCode  VecSquare(Vec v)
 26: {
 27:   PetscScalar    *x;
 28:   PetscInt       i, n;

 30:   VecGetLocalSize(v, &n);
 31:   VecGetArray(v, &x);
 32:   for (i = 0; i < n; i++) x[i] *= PetscConj(x[i]);
 33:   VecRestoreArray(v, &x);
 34:   return 0;
 35: }

 37: static PetscErrorCode KSPSetUp_LSQR(KSP ksp)
 38: {
 39:   KSP_LSQR       *lsqr = (KSP_LSQR*)ksp->data;
 40:   PetscBool      nopreconditioner;

 42:   PetscObjectTypeCompare((PetscObject)ksp->pc,PCNONE,&nopreconditioner);

 44:   if (lsqr->vwork_m) {
 45:     VecDestroyVecs(lsqr->nwork_m,&lsqr->vwork_m);
 46:   }

 48:   if (lsqr->vwork_n) {
 49:     VecDestroyVecs(lsqr->nwork_n,&lsqr->vwork_n);
 50:   }

 52:   lsqr->nwork_m = 2;
 53:   if (nopreconditioner) lsqr->nwork_n = 4;
 54:   else lsqr->nwork_n = 5;
 55:   KSPCreateVecs(ksp,lsqr->nwork_n,&lsqr->vwork_n,lsqr->nwork_m,&lsqr->vwork_m);

 57:   if (lsqr->se_flg && !lsqr->se) {
 58:     VecDuplicate(lsqr->vwork_n[0],&lsqr->se);
 59:     VecSet(lsqr->se,PETSC_INFINITY);
 60:   } else if (!lsqr->se_flg) {
 61:     VecDestroy(&lsqr->se);
 62:   }
 63:   return 0;
 64: }

 66: static PetscErrorCode KSPSolve_LSQR(KSP ksp)
 67: {
 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;

 76:   PCGetDiagonalScale(ksp->pc,&diagonalscale);

 79:   PCGetOperators(ksp->pc,&Amat,&Pmat);
 80:   PetscObjectTypeCompare((PetscObject)ksp->pc,PCNONE,&nopreconditioner);

 82:   /* vectors of length m, where system size is mxn */
 83:   B  = ksp->vec_rhs;
 84:   U  = lsqr->vwork_m[0];
 85:   U1 = lsqr->vwork_m[1];

 87:   /* vectors of length n */
 88:   X  = ksp->vec_sol;
 89:   W  = lsqr->vwork_n[0];
 90:   V  = lsqr->vwork_n[1];
 91:   V1 = lsqr->vwork_n[2];
 92:   W2 = lsqr->vwork_n[3];
 93:   if (!nopreconditioner) Z = lsqr->vwork_n[4];

 95:   /* standard error vector */
 96:   if (lsqr->se) {
 97:     VecSet(lsqr->se,0.0);
 98:   }

100:   /* Compute initial residual, temporarily use work vector u */
101:   if (!ksp->guess_zero) {
102:     KSP_MatMult(ksp,Amat,X,U);       /*   u <- b - Ax     */
103:     VecAYPX(U,-1.0,B);
104:   } else {
105:     VecCopy(B,U);            /*   u <- b (x is 0) */
106:   }

108:   /* Test for nothing to do */
109:   VecNorm(U,NORM_2,&rnorm);
110:   KSPCheckNorm(ksp,rnorm);
111:   PetscObjectSAWsTakeAccess((PetscObject)ksp);
112:   ksp->its   = 0;
113:   ksp->rnorm = rnorm;
114:   PetscObjectSAWsGrantAccess((PetscObject)ksp);
115:   KSPLogResidualHistory(ksp,rnorm);
116:   KSPMonitor(ksp,0,rnorm);
117:   (*ksp->converged)(ksp,0,rnorm,&ksp->reason,ksp->cnvP);
118:   if (ksp->reason) return 0;

120:   beta = rnorm;
121:   VecScale(U,1.0/beta);
122:   KSP_MatMultHermitianTranspose(ksp,Amat,U,V);
123:   if (nopreconditioner) {
124:     VecNorm(V,NORM_2,&alpha);
125:     KSPCheckNorm(ksp,rnorm);
126:   } else {
127:     /* this is an application of the preconditioner for the normal equations; not the operator, see the manual page */
128:     PCApply(ksp->pc,V,Z);
129:     VecDotRealPart(V,Z,&alpha);
130:     if (alpha <= 0.0) {
131:       ksp->reason = KSP_DIVERGED_BREAKDOWN;
132:       return 0;
133:     }
134:     alpha = PetscSqrtReal(alpha);
135:     VecScale(Z,1.0/alpha);
136:   }
137:   VecScale(V,1.0/alpha);

139:   if (nopreconditioner) {
140:     VecCopy(V,W);
141:   } else {
142:     VecCopy(Z,W);
143:   }

145:   if (lsqr->exact_norm) {
146:     MatNorm(Amat,NORM_FROBENIUS,&lsqr->anorm);
147:   } else lsqr->anorm = 0.0;

149:   lsqr->arnorm = alpha * beta;
150:   phibar       = beta;
151:   rhobar       = alpha;
152:   i            = 0;
153:   do {
154:     if (nopreconditioner) {
155:       KSP_MatMult(ksp,Amat,V,U1);
156:     } else {
157:       KSP_MatMult(ksp,Amat,Z,U1);
158:     }
159:     VecAXPY(U1,-alpha,U);
160:     VecNorm(U1,NORM_2,&beta);
161:     KSPCheckNorm(ksp,beta);
162:     if (beta > 0.0) {
163:       VecScale(U1,1.0/beta); /* beta*U1 = Amat*V - alpha*U */
164:       if (!lsqr->exact_norm) {
165:         lsqr->anorm = PetscSqrtReal(PetscSqr(lsqr->anorm) + PetscSqr(alpha) + PetscSqr(beta));
166:       }
167:     }

169:     KSP_MatMultHermitianTranspose(ksp,Amat,U1,V1);
170:     VecAXPY(V1,-beta,V);
171:     if (nopreconditioner) {
172:       VecNorm(V1,NORM_2,&alpha);
173:       KSPCheckNorm(ksp,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 (lsqr->se) {
197:       VecCopy(W,W2);
198:       VecSquare(W2);
199:       VecScale(W2,1.0/(rho*rho));
200:       VecAXPY(lsqr->se, 1.0, W2); /* lsqr->se <- lsqr->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 (lsqr->se) {
228:     tmp  = 1.0;
229:     MatGetSize(Amat,&size1,&size2);
230:     if (size1 > size2) tmp = size1 - size2;
231:     tmp  = rnorm / PetscSqrtScalar(tmp);
232:     VecSqrtAbs(lsqr->se);
233:     VecScale(lsqr->se,tmp);
234:   }
235:   return 0;
236: }

238: PetscErrorCode KSPDestroy_LSQR(KSP ksp)
239: {
240:   KSP_LSQR       *lsqr = (KSP_LSQR*)ksp->data;

242:   /* Free work vectors */
243:   if (lsqr->vwork_n) {
244:     VecDestroyVecs(lsqr->nwork_n,&lsqr->vwork_n);
245:   }
246:   if (lsqr->vwork_m) {
247:     VecDestroyVecs(lsqr->nwork_m,&lsqr->vwork_m);
248:   }
249:   VecDestroy(&lsqr->se);
250:   /* Revert convergence test */
251:   KSPSetConvergenceTest(ksp,lsqr->converged,lsqr->cnvP,lsqr->convergeddestroy);
252:   /* Free the KSP_LSQR context */
253:   PetscFree(ksp->data);
254:   PetscObjectComposeFunction((PetscObject)ksp,"KSPLSQRMonitorResidual_C",NULL);
255:   PetscObjectComposeFunction((PetscObject)ksp,"KSPLSQRMonitorResidualDrawLG_C",NULL);
256:   return 0;
257: }

259: /*@
260:    KSPLSQRSetComputeStandardErrorVec - Compute vector of standard error estimates during KSPSolve_LSQR().

262:    Not Collective

264:    Input Parameters:
265: +  ksp   - iterative context
266: -  flg   - compute the vector of standard estimates or not

268:    Developer notes:
269:    Vaclav: I'm not sure whether this vector is useful for anything.

271:    Level: intermediate

273: .seealso: KSPSolve(), KSPLSQR, KSPLSQRGetStandardErrorVec()
274: @*/
275: PetscErrorCode  KSPLSQRSetComputeStandardErrorVec(KSP ksp, PetscBool flg)
276: {
277:   KSP_LSQR       *lsqr = (KSP_LSQR*)ksp->data;

279:   lsqr->se_flg = flg;
280:   return 0;
281: }

283: /*@
284:    KSPLSQRSetExactMatNorm - Compute exact matrix norm instead of iteratively refined estimate.

286:    Not Collective

288:    Input Parameters:
289: +  ksp   - iterative context
290: -  flg   - compute exact matrix norm or not

292:    Notes:
293:    By default, flg=PETSC_FALSE. This is usually preferred to avoid possibly expensive computation of the norm.
294:    For flg=PETSC_TRUE, we call MatNorm(Amat,NORM_FROBENIUS,&lsqr->anorm) which will work only for some types of explicitly assembled matrices.
295:    This can affect convergence rate as KSPLSQRConvergedDefault() assumes different value of ||A|| used in normal equation stopping criterion.

297:    Level: intermediate

299: .seealso: KSPSolve(), KSPLSQR, KSPLSQRGetNorms(), KSPLSQRConvergedDefault()
300: @*/
301: PetscErrorCode  KSPLSQRSetExactMatNorm(KSP ksp, PetscBool flg)
302: {
303:   KSP_LSQR       *lsqr = (KSP_LSQR*)ksp->data;

305:   lsqr->exact_norm = flg;
306:   return 0;
307: }

309: /*@
310:    KSPLSQRGetStandardErrorVec - Get vector of standard error estimates.
311:    Only available if -ksp_lsqr_set_standard_error was set to true
312:    or KSPLSQRSetComputeStandardErrorVec(ksp, PETSC_TRUE) was called.
313:    Otherwise returns NULL.

315:    Not Collective

317:    Input Parameters:
318: .  ksp   - iterative context

320:    Output Parameters:
321: .  se - vector of standard estimates

323:    Options Database Keys:
324: .   -ksp_lsqr_set_standard_error  - set standard error estimates of solution

326:    Developer notes:
327:    Vaclav: I'm not sure whether this vector is useful for anything.

329:    Level: intermediate

331: .seealso: KSPSolve(), KSPLSQR, KSPLSQRSetComputeStandardErrorVec()
332: @*/
333: PetscErrorCode  KSPLSQRGetStandardErrorVec(KSP ksp,Vec *se)
334: {
335:   KSP_LSQR *lsqr = (KSP_LSQR*)ksp->data;

337:   *se = lsqr->se;
338:   return 0;
339: }

341: /*@
342:    KSPLSQRGetNorms - Get norm estimates that LSQR computes internally during KSPSolve().

344:    Not Collective

346:    Input Parameter:
347: .  ksp   - iterative context

349:    Output Parameters:
350: +  arnorm - good estimate of norm((A*inv(Pmat))'*r), where r = A*x - b, used in specific stopping criterion
351: -  anorm - poor estimate of norm(A*inv(Pmat),'fro') used in specific stopping criterion

353:    Notes:
354:    Output parameters are meaningful only after KSPSolve().
355:    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.
356:    If -ksp_lsqr_exact_mat_norm is set or KSPLSQRSetExactMatNorm(ksp, PETSC_TRUE) called, then anorm is exact Frobenius norm.

358:    Level: intermediate

360: .seealso: KSPSolve(), KSPLSQR, KSPLSQRSetExactMatNorm()
361: @*/
362: PetscErrorCode  KSPLSQRGetNorms(KSP ksp,PetscReal *arnorm, PetscReal *anorm)
363: {
364:   KSP_LSQR       *lsqr = (KSP_LSQR*)ksp->data;

366:   if (arnorm)   *arnorm = lsqr->arnorm;
367:   if (anorm)    *anorm = lsqr->anorm;
368:   return 0;
369: }

371: PetscErrorCode KSPLSQRMonitorResidual_LSQR(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
372: {
373:   KSP_LSQR          *lsqr  = (KSP_LSQR*)ksp->data;
374:   PetscViewer       viewer = vf->viewer;
375:   PetscViewerFormat format = vf->format;
376:   char              normtype[256];
377:   PetscInt          tablevel;
378:   const char        *prefix;

380:   PetscObjectGetTabLevel((PetscObject) ksp, &tablevel);
381:   PetscObjectGetOptionsPrefix((PetscObject) ksp, &prefix);
382:   PetscStrncpy(normtype, KSPNormTypes[ksp->normtype], sizeof(normtype));
383:   PetscStrtolower(normtype);
384:   PetscViewerPushFormat(viewer, format);
385:   PetscViewerASCIIAddTab(viewer, tablevel);
386:   if (n == 0 && prefix) PetscViewerASCIIPrintf(viewer, "  Residual norm, norm of normal equations, and matrix norm for %s solve.\n", prefix);
387:   if (!n) {
388:     PetscViewerASCIIPrintf(viewer,"%3D KSP resid norm %14.12e\n",n,(double)rnorm);
389:   } else {
390:     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);
391:   }
392:   PetscViewerASCIISubtractTab(viewer, tablevel);
393:   PetscViewerPopFormat(viewer);
394:   return 0;
395: }

397: /*@C
398:   KSPLSQRMonitorResidual - Prints the residual norm, as well as the normal equation residual norm, at each iteration of an iterative solver.

400:   Collective on ksp

402:   Input Parameters:
403: + ksp   - iterative context
404: . n     - iteration number
405: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
406: - vf    - The viewer context

408:   Options Database Key:
409: . -ksp_lsqr_monitor - Activates KSPLSQRMonitorResidual()

411:   Level: intermediate

413: .seealso: KSPMonitorSet(), KSPMonitorResidual(),KSPMonitorTrueResidualMaxNorm()
414: @*/
415: PetscErrorCode KSPLSQRMonitorResidual(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
416: {
420:   PetscTryMethod(ksp, "KSPLSQRMonitorResidual_C", (KSP,PetscInt,PetscReal,PetscViewerAndFormat*), (ksp,n,rnorm,vf));
421:   return 0;
422: }

424: PetscErrorCode KSPLSQRMonitorResidualDrawLG_LSQR(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
425: {
426:   KSP_LSQR           *lsqr  = (KSP_LSQR*)ksp->data;
427:   PetscViewer        viewer = vf->viewer;
428:   PetscViewerFormat  format = vf->format;
429:   PetscDrawLG        lg     = vf->lg;
430:   KSPConvergedReason reason;
431:   PetscReal          x[2], y[2];

433:   PetscViewerPushFormat(viewer, format);
434:   if (!n) PetscDrawLGReset(lg);
435:   x[0] = (PetscReal) n;
436:   if (rnorm > 0.0) y[0] = PetscLog10Real(rnorm);
437:   else y[0] = -15.0;
438:   x[1] = (PetscReal) n;
439:   if (lsqr->arnorm > 0.0) y[1] = PetscLog10Real(lsqr->arnorm);
440:   else y[1] = -15.0;
441:   PetscDrawLGAddPoint(lg, x, y);
442:   KSPGetConvergedReason(ksp, &reason);
443:   if (n <= 20 || !(n % 5) || reason) {
444:     PetscDrawLGDraw(lg);
445:     PetscDrawLGSave(lg);
446:   }
447:   PetscViewerPopFormat(viewer);
448:   return 0;
449: }

451: /*@C
452:   KSPLSQRMonitorResidualDrawLG - Plots the true residual norm at each iteration of an iterative solver.

454:   Collective on ksp

456:   Input Parameters:
457: + ksp   - iterative context
458: . n     - iteration number
459: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
460: - vf    - The viewer context

462:   Options Database Key:
463: . -ksp_lsqr_monitor draw::draw_lg - Activates KSPMonitorTrueResidualDrawLG()

465:   Level: intermediate

467: .seealso: KSPMonitorSet(), KSPMonitorTrueResidual()
468: @*/
469: PetscErrorCode KSPLSQRMonitorResidualDrawLG(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
470: {
475:   PetscTryMethod(ksp, "KSPLSQRMonitorResidualDrawLG_C", (KSP,PetscInt,PetscReal,PetscViewerAndFormat*), (ksp,n,rnorm,vf));
476:   return 0;
477: }

479: /*@C
480:   KSPLSQRMonitorResidualDrawLGCreate - Creates the plotter for the LSQR residual and normal eqn residual.

482:   Collective on ksp

484:   Input Parameters:
485: + viewer - The PetscViewer
486: . format - The viewer format
487: - ctx    - An optional user context

489:   Output Parameter:
490: . vf    - The viewer context

492:   Level: intermediate

494: .seealso: KSPMonitorSet(), KSPLSQRMonitorResidual()
495: @*/
496: PetscErrorCode KSPLSQRMonitorResidualDrawLGCreate(PetscViewer viewer, PetscViewerFormat format, void *ctx, PetscViewerAndFormat **vf)
497: {
498:   const char    *names[] = {"residual", "normal eqn residual"};

500:   PetscViewerAndFormatCreate(viewer, format, vf);
501:   (*vf)->data = ctx;
502:   KSPMonitorLGCreate(PetscObjectComm((PetscObject) viewer), NULL, NULL, "Log Residual Norm", 2, names, PETSC_DECIDE, PETSC_DECIDE, 400, 300, &(*vf)->lg);
503:   return 0;
504: }

506: PetscErrorCode KSPSetFromOptions_LSQR(PetscOptionItems *PetscOptionsObject,KSP ksp)
507: {
508:   KSP_LSQR       *lsqr = (KSP_LSQR*)ksp->data;

510:   PetscOptionsHead(PetscOptionsObject,"KSP LSQR Options");
511:   PetscOptionsBool("-ksp_lsqr_compute_standard_error","Set Standard Error Estimates of Solution","KSPLSQRSetComputeStandardErrorVec",lsqr->se_flg,&lsqr->se_flg,NULL);
512:   PetscOptionsBool("-ksp_lsqr_exact_mat_norm","Compute exact matrix norm instead of iteratively refined estimate","KSPLSQRSetExactMatNorm",lsqr->exact_norm,&lsqr->exact_norm,NULL);
513:   KSPMonitorSetFromOptions(ksp, "-ksp_lsqr_monitor", "lsqr_residual", NULL);
514:   PetscOptionsTail();
515:   return 0;
516: }

518: PetscErrorCode KSPView_LSQR(KSP ksp,PetscViewer viewer)
519: {
520:   KSP_LSQR       *lsqr = (KSP_LSQR*)ksp->data;
521:   PetscBool      iascii;

523:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
524:   if (iascii) {
525:     if (lsqr->se) {
526:       PetscReal rnorm;
527:       VecNorm(lsqr->se,NORM_2,&rnorm);
528:       PetscViewerASCIIPrintf(viewer,"  norm of standard error %g, iterations %d\n",(double)rnorm,ksp->its);
529:     } else {
530:       PetscViewerASCIIPrintf(viewer,"  standard error not computed\n");
531:     }
532:     if (lsqr->exact_norm) {
533:       PetscViewerASCIIPrintf(viewer,"  using exact matrix norm\n");
534:     } else {
535:       PetscViewerASCIIPrintf(viewer,"  using inexact matrix norm\n");
536:     }
537:   }
538:   return 0;
539: }

541: /*@C
542:    KSPLSQRConvergedDefault - Determines convergence of the LSQR Krylov method.

544:    Collective on ksp

546:    Input Parameters:
547: +  ksp   - iterative context
548: .  n     - iteration number
549: .  rnorm - 2-norm residual value (may be estimated)
550: -  ctx - convergence context which must be created by KSPConvergedDefaultCreate()

552:    reason is set to:
553: +   positive - if the iteration has converged;
554: .   negative - if residual norm exceeds divergence threshold;
555: -   0 - otherwise.

557:    Notes:
558:    KSPConvergedDefault() is called first to check for convergence in A*x=b.
559:    If that does not determine convergence then checks convergence for the least squares problem, i.e. in min{|b-A*x|}.
560:    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.
561:    KSP_CONVERGED_RTOL_NORMAL is returned if ||A'*r|| < rtol * ||A|| * ||r||.
562:    Matrix norm ||A|| is iteratively refined estimate, see KSPLSQRGetNorms().
563:    This criterion is now largely compatible with that in MATLAB lsqr().

565:    Level: intermediate

567: .seealso: KSPLSQR, KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(),
568:           KSPConvergedDefaultSetUIRNorm(), KSPConvergedDefaultSetUMIRNorm(), KSPConvergedDefaultCreate(), KSPConvergedDefaultDestroy(), KSPConvergedDefault(), KSPLSQRGetNorms(), KSPLSQRSetExactMatNorm()
569: @*/
570: PetscErrorCode  KSPLSQRConvergedDefault(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *ctx)
571: {
572:   KSP_LSQR       *lsqr = (KSP_LSQR*)ksp->data;

574:   /* check for convergence in A*x=b */
575:   KSPConvergedDefault(ksp,n,rnorm,reason,ctx);
576:   if (!n || *reason) return 0;

578:   /* check for convergence in min{|b-A*x|} */
579:   if (lsqr->arnorm < ksp->abstol) {
580:     PetscInfo(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);
581:     *reason = KSP_CONVERGED_ATOL_NORMAL;
582:   } else if (lsqr->arnorm < ksp->rtol * lsqr->anorm * rnorm) {
583:     PetscInfo(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);
584:     *reason = KSP_CONVERGED_RTOL_NORMAL;
585:   }
586:   return 0;
587: }

589: /*MC
590:      KSPLSQR - This implements LSQR

592:    Options Database Keys:
593: +   -ksp_lsqr_set_standard_error  - set standard error estimates of solution, see KSPLSQRSetComputeStandardErrorVec() and KSPLSQRGetStandardErrorVec()
594: .   -ksp_lsqr_exact_mat_norm - compute exact matrix norm instead of iteratively refined estimate, see KSPLSQRSetExactMatNorm()
595: -   -ksp_lsqr_monitor - monitor residual norm, norm of residual of normal equations A'*A x = A' b, and estimate of matrix norm ||A||

597:    Level: beginner

599:    Notes:
600:      Supports non-square (rectangular) matrices.

602:      This variant, when applied with no preconditioning is identical to the original algorithm in exact arithematic; however, in practice, with no preconditioning
603:      due to inexact arithmetic, it can converge differently. Hence when no preconditioner is used (PCType PCNONE) it automatically reverts to the original algorithm.

605:      With the PETSc built-in preconditioners, such as ICC, one should call KSPSetOperators(ksp,A,A'*A)) since the preconditioner needs to work
606:      for the normal equations A'*A.

608:      Supports only left preconditioning.

610:      For least squares problems with nonzero residual A*x - b, there are additional convergence tests for the residual of the normal equations, A'*(b - Ax), see KSPLSQRConvergedDefault().

612:    References:
613: .  * - The original unpreconditioned algorithm can be found in Paige and Saunders, ACM Transactions on Mathematical Software, Vol 8, 1982.

615:      In exact arithmetic the LSQR method (with no preconditioning) is identical to the KSPCG algorithm applied to the normal equations.
616:      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
617:      track the true norm of the residual and uses that in the convergence test.

619:    Developer Notes:
620:     How is this related to the KSPCGNE implementation? One difference is that KSPCGNE applies
621:             the preconditioner transpose times the preconditioner,  so one does not need to pass A'*A as the third argument to KSPSetOperators().

623: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPSolve(), KSPLSQRConvergedDefault(), KSPLSQRSetComputeStandardErrorVec(), KSPLSQRGetStandardErrorVec(), KSPLSQRSetExactMatNorm()

625: M*/
626: PETSC_EXTERN PetscErrorCode KSPCreate_LSQR(KSP ksp)
627: {
628:   KSP_LSQR       *lsqr;
629:   void           *ctx;

631:   PetscNewLog(ksp,&lsqr);
632:   lsqr->se     = NULL;
633:   lsqr->se_flg = PETSC_FALSE;
634:   lsqr->exact_norm = PETSC_FALSE;
635:   lsqr->anorm  = -1.0;
636:   lsqr->arnorm = -1.0;
637:   ksp->data    = (void*)lsqr;
638:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,3);

640:   ksp->ops->setup          = KSPSetUp_LSQR;
641:   ksp->ops->solve          = KSPSolve_LSQR;
642:   ksp->ops->destroy        = KSPDestroy_LSQR;
643:   ksp->ops->setfromoptions = KSPSetFromOptions_LSQR;
644:   ksp->ops->view           = KSPView_LSQR;

646:   /* Backup current convergence test; remove destroy routine from KSP to prevent destroying the convergence context in KSPSetConvergenceTest() */
647:   KSPGetAndClearConvergenceTest(ksp,&lsqr->converged,&lsqr->cnvP,&lsqr->convergeddestroy);
648:   /* Override current convergence test */
649:   KSPConvergedDefaultCreate(&ctx);
650:   KSPSetConvergenceTest(ksp,KSPLSQRConvergedDefault,ctx,KSPConvergedDefaultDestroy);
651:   PetscObjectComposeFunction((PetscObject)ksp,"KSPLSQRMonitorResidual_C",KSPLSQRMonitorResidual_LSQR);
652:   PetscObjectComposeFunction((PetscObject)ksp,"KSPLSQRMonitorResidualDrawLG_C",KSPLSQRMonitorResidualDrawLG_LSQR);
653:   return 0;
654: }