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) \
  6:   { \
  7:     c = a; \
  8:     a = b; \
  9:     b = c; \
 10:   }

 12: #include <petsc/private/kspimpl.h>
 13: #include <petscdraw.h>

 15: typedef struct {
 16:   PetscInt  nwork_n, nwork_m;
 17:   Vec      *vwork_m;    /* work vectors of length m, where the system is size m x n */
 18:   Vec      *vwork_n;    /* work vectors of length n */
 19:   Vec       se;         /* Optional standard error vector */
 20:   PetscBool se_flg;     /* flag for -ksp_lsqr_set_standard_error */
 21:   PetscBool exact_norm; /* flag for -ksp_lsqr_exact_mat_norm */
 22:   PetscReal arnorm;     /* Good estimate of norm((A*inv(Pmat))'*r), where r = A*x - b, used in specific stopping criterion */
 23:   PetscReal anorm;      /* Poor estimate of norm(A*inv(Pmat),'fro') used in specific stopping criterion */
 24:   /* Backup previous convergence test */
 25:   PetscErrorCode (*converged)(KSP, PetscInt, PetscReal, KSPConvergedReason *, void *);
 26:   PetscErrorCode (*convergeddestroy)(void *);
 27:   void *cnvP;
 28: } KSP_LSQR;

 30: static PetscErrorCode VecSquare(Vec v)
 31: {
 32:   PetscScalar *x;
 33:   PetscInt     i, n;

 35:   VecGetLocalSize(v, &n);
 36:   VecGetArray(v, &x);
 37:   for (i = 0; i < n; i++) x[i] *= PetscConj(x[i]);
 38:   VecRestoreArray(v, &x);
 39:   return 0;
 40: }

 42: static PetscErrorCode KSPSetUp_LSQR(KSP ksp)
 43: {
 44:   KSP_LSQR *lsqr = (KSP_LSQR *)ksp->data;
 45:   PetscBool nopreconditioner;

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

 49:   if (lsqr->vwork_m) VecDestroyVecs(lsqr->nwork_m, &lsqr->vwork_m);

 51:   if (lsqr->vwork_n) VecDestroyVecs(lsqr->nwork_n, &lsqr->vwork_n);

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

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

 67: static PetscErrorCode KSPSolve_LSQR(KSP ksp)
 68: {
 69:   PetscInt    i, size1, size2;
 70:   PetscScalar rho, rhobar, phi, phibar, theta, c, s, tmp, tau;
 71:   PetscReal   beta, alpha, rnorm;
 72:   Vec         X, B, V, V1, U, U1, TMP, W, W2, Z = NULL;
 73:   Mat         Amat, Pmat;
 74:   KSP_LSQR   *lsqr = (KSP_LSQR *)ksp->data;
 75:   PetscBool   diagonalscale, nopreconditioner;

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

 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) VecSet(lsqr->se, 0.0);

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

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

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

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

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

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

166:     KSP_MatMultHermitianTranspose(ksp, Amat, U1, V1);
167:     VecAXPY(V1, -beta, V);
168:     if (nopreconditioner) {
169:       VecNorm(V1, NORM_2, &alpha);
170:       KSPCheckNorm(ksp, alpha);
171:     } else {
172:       PCApply(ksp->pc, V1, Z);
173:       VecDotRealPart(V1, Z, &alpha);
174:       if (alpha <= 0.0) {
175:         ksp->reason = KSP_DIVERGED_BREAKDOWN;
176:         break;
177:       }
178:       alpha = PetscSqrtReal(alpha);
179:       VecScale(Z, 1.0 / alpha);
180:     }
181:     VecScale(V1, 1.0 / alpha); /* alpha*V1 = Amat^T*U1 - beta*V */
182:     rho    = PetscSqrtScalar(rhobar * rhobar + beta * beta);
183:     c      = rhobar / rho;
184:     s      = beta / rho;
185:     theta  = s * alpha;
186:     rhobar = -c * alpha;
187:     phi    = c * phibar;
188:     phibar = s * phibar;
189:     tau    = s * phi;

191:     VecAXPY(X, phi / rho, W); /*    x <- x + (phi/rho) w   */

193:     if (lsqr->se) {
194:       VecCopy(W, W2);
195:       VecSquare(W2);
196:       VecScale(W2, 1.0 / (rho * rho));
197:       VecAXPY(lsqr->se, 1.0, W2); /* lsqr->se <- lsqr->se + (w^2/rho^2) */
198:     }
199:     if (nopreconditioner) {
200:       VecAYPX(W, -theta / rho, V1); /* w <- v - (theta/rho) w */
201:     } else {
202:       VecAYPX(W, -theta / rho, Z); /* w <- z - (theta/rho) w */
203:     }

205:     lsqr->arnorm = alpha * PetscAbsScalar(tau);
206:     rnorm        = PetscRealPart(phibar);

208:     PetscObjectSAWsTakeAccess((PetscObject)ksp);
209:     ksp->its++;
210:     ksp->rnorm = rnorm;
211:     PetscObjectSAWsGrantAccess((PetscObject)ksp);
212:     KSPLogResidualHistory(ksp, rnorm);
213:     KSPMonitor(ksp, i + 1, rnorm);
214:     (*ksp->converged)(ksp, i + 1, rnorm, &ksp->reason, ksp->cnvP);
215:     if (ksp->reason) break;
216:     SWAP(U1, U, TMP);
217:     SWAP(V1, V, TMP);

219:     i++;
220:   } while (i < ksp->max_it);
221:   if (i >= ksp->max_it && !ksp->reason) ksp->reason = KSP_DIVERGED_ITS;

223:   /* Finish off the standard error estimates */
224:   if (lsqr->se) {
225:     tmp = 1.0;
226:     MatGetSize(Amat, &size1, &size2);
227:     if (size1 > size2) tmp = size1 - size2;
228:     tmp = rnorm / PetscSqrtScalar(tmp);
229:     VecSqrtAbs(lsqr->se);
230:     VecScale(lsqr->se, tmp);
231:   }
232:   return 0;
233: }

235: PetscErrorCode KSPDestroy_LSQR(KSP ksp)
236: {
237:   KSP_LSQR *lsqr = (KSP_LSQR *)ksp->data;

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

252: /*@
253:    KSPLSQRSetComputeStandardErrorVec - Compute a vector of standard error estimates during `KSPSolve()` for  `KSPLSQR`.

255:    Logically Collective

257:    Input Parameters:
258: +  ksp   - iterative context
259: -  flg   - compute the vector of standard estimates or not

261:    Level: intermediate

263:    Developer Note:
264:    Vaclav: I'm not sure whether this vector is useful for anything.

266: .seealso: [](chapter_ksp), `KSPSolve()`, `KSPLSQR`, `KSPLSQRGetStandardErrorVec()`
267: @*/
268: PetscErrorCode KSPLSQRSetComputeStandardErrorVec(KSP ksp, PetscBool flg)
269: {
270:   KSP_LSQR *lsqr = (KSP_LSQR *)ksp->data;

272:   lsqr->se_flg = flg;
273:   return 0;
274: }

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

279:    Not Collective

281:    Input Parameters:
282: +  ksp   - iterative context
283: -  flg   - compute exact matrix norm or not

285:    Level: intermediate

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

292: .seealso: [](chapter_ksp), `KSPSolve()`, `KSPLSQR`, `KSPLSQRGetNorms()`, `KSPLSQRConvergedDefault()`
293: @*/
294: PetscErrorCode KSPLSQRSetExactMatNorm(KSP ksp, PetscBool flg)
295: {
296:   KSP_LSQR *lsqr = (KSP_LSQR *)ksp->data;

298:   lsqr->exact_norm = flg;
299:   return 0;
300: }

302: /*@
303:    KSPLSQRGetStandardErrorVec - Get vector of standard error estimates.
304:    Only available if -ksp_lsqr_set_standard_error was set to true
305:    or `KSPLSQRSetComputeStandardErrorVec`(ksp, `PETSC_TRUE`) was called.
306:    Otherwise returns NULL.

308:    Not Collective

310:    Input Parameters:
311: .  ksp   - iterative context

313:    Output Parameters:
314: .  se - vector of standard estimates

316:    Level: intermediate

318:    Developer Note:
319:    Vaclav: I'm not sure whether this vector is useful for anything.

321: .seealso: [](chapter_ksp), `KSPSolve()`, `KSPLSQR`, `KSPLSQRSetComputeStandardErrorVec()`
322: @*/
323: PetscErrorCode KSPLSQRGetStandardErrorVec(KSP ksp, Vec *se)
324: {
325:   KSP_LSQR *lsqr = (KSP_LSQR *)ksp->data;

327:   *se = lsqr->se;
328:   return 0;
329: }

331: /*@
332:    KSPLSQRGetNorms - Get the norm estimates that `KSPLSQR` computes internally during `KSPSolve()`.

334:    Not Collective

336:    Input Parameter:
337: .  ksp   - iterative context

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

343:    Notes:
344:    Output parameters are meaningful only after `KSPSolve()`.

346:    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.

348:    If -ksp_lsqr_exact_mat_norm is set or `KSPLSQRSetExactMatNorm`(ksp, `PETSC_TRUE`) called, then anorm is the exact Frobenius norm.

350:    Level: intermediate

352: .seealso: [](chapter_ksp), `KSPSolve()`, `KSPLSQR`, `KSPLSQRSetExactMatNorm()`
353: @*/
354: PetscErrorCode KSPLSQRGetNorms(KSP ksp, PetscReal *arnorm, PetscReal *anorm)
355: {
356:   KSP_LSQR *lsqr = (KSP_LSQR *)ksp->data;

358:   if (arnorm) *arnorm = lsqr->arnorm;
359:   if (anorm) *anorm = lsqr->anorm;
360:   return 0;
361: }

363: PetscErrorCode KSPLSQRMonitorResidual_LSQR(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
364: {
365:   KSP_LSQR         *lsqr   = (KSP_LSQR *)ksp->data;
366:   PetscViewer       viewer = vf->viewer;
367:   PetscViewerFormat format = vf->format;
368:   char              normtype[256];
369:   PetscInt          tablevel;
370:   const char       *prefix;

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

389: /*@C
390:   KSPLSQRMonitorResidual - Prints the residual norm, as well as the normal equation residual norm, at each iteration of an iterative solver for the `KSPLSQR` solver

392:   Collective

394:   Input Parameters:
395: + ksp   - iterative context
396: . n     - iteration number
397: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
398: - vf    - The viewer context

400:   Options Database Key:
401: . -ksp_lsqr_monitor - Activates `KSPLSQRMonitorResidual()`

403:   Level: intermediate

405: .seealso: [](chapter_ksp), `KSPLSQR`, `KSPMonitorSet()`, `KSPMonitorResidual()`, `KSPMonitorTrueResidualMaxNorm()`, `KSPLSQRMonitorResidualDrawLG()`
406: @*/
407: PetscErrorCode KSPLSQRMonitorResidual(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
408: {
412:   PetscTryMethod(ksp, "KSPLSQRMonitorResidual_C", (KSP, PetscInt, PetscReal, PetscViewerAndFormat *), (ksp, n, rnorm, vf));
413:   return 0;
414: }

416: PetscErrorCode KSPLSQRMonitorResidualDrawLG_LSQR(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
417: {
418:   KSP_LSQR          *lsqr   = (KSP_LSQR *)ksp->data;
419:   PetscViewer        viewer = vf->viewer;
420:   PetscViewerFormat  format = vf->format;
421:   PetscDrawLG        lg     = vf->lg;
422:   KSPConvergedReason reason;
423:   PetscReal          x[2], y[2];

425:   PetscViewerPushFormat(viewer, format);
426:   if (!n) PetscDrawLGReset(lg);
427:   x[0] = (PetscReal)n;
428:   if (rnorm > 0.0) y[0] = PetscLog10Real(rnorm);
429:   else y[0] = -15.0;
430:   x[1] = (PetscReal)n;
431:   if (lsqr->arnorm > 0.0) y[1] = PetscLog10Real(lsqr->arnorm);
432:   else y[1] = -15.0;
433:   PetscDrawLGAddPoint(lg, x, y);
434:   KSPGetConvergedReason(ksp, &reason);
435:   if (n <= 20 || !(n % 5) || reason) {
436:     PetscDrawLGDraw(lg);
437:     PetscDrawLGSave(lg);
438:   }
439:   PetscViewerPopFormat(viewer);
440:   return 0;
441: }

443: /*@C
444:   KSPLSQRMonitorResidualDrawLG - Plots the true residual norm at each iteration of an iterative solver for the `KSPLSQR` solver

446:   Collective

448:   Input Parameters:
449: + ksp   - iterative context
450: . n     - iteration number
451: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
452: - vf    - The viewer context

454:   Options Database Key:
455: . -ksp_lsqr_monitor draw::draw_lg - Activates `KSPMonitorTrueResidualDrawLG()`

457:   Level: intermediate

459: .seealso: [](chapter_ksp), `KSPLSQR`, `KSPMonitorSet()`, `KSPMonitorTrueResidual()`, `KSPLSQRMonitorResidual()`, `KSPLSQRMonitorResidualDrawLGCreate()`
460: @*/
461: PetscErrorCode KSPLSQRMonitorResidualDrawLG(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
462: {
467:   PetscTryMethod(ksp, "KSPLSQRMonitorResidualDrawLG_C", (KSP, PetscInt, PetscReal, PetscViewerAndFormat *), (ksp, n, rnorm, vf));
468:   return 0;
469: }

471: /*@C
472:   KSPLSQRMonitorResidualDrawLGCreate - Creates the plotter for the `KSPLSQR` residual and normal equation residual norm

474:   Collective

476:   Input Parameters:
477: + viewer - The PetscViewer
478: . format - The viewer format
479: - ctx    - An optional user context

481:   Output Parameter:
482: . vf    - The viewer context

484:   Level: intermediate

486: .seealso: [](chapter_ksp), `KSPLSQR`, `KSPMonitorSet()`, `KSPLSQRMonitorResidual()`, `KSPLSQRMonitorResidualDrawLG()`
487: @*/
488: PetscErrorCode KSPLSQRMonitorResidualDrawLGCreate(PetscViewer viewer, PetscViewerFormat format, void *ctx, PetscViewerAndFormat **vf)
489: {
490:   const char *names[] = {"residual", "normal eqn residual"};

492:   PetscViewerAndFormatCreate(viewer, format, vf);
493:   (*vf)->data = ctx;
494:   KSPMonitorLGCreate(PetscObjectComm((PetscObject)viewer), NULL, NULL, "Log Residual Norm", 2, names, PETSC_DECIDE, PETSC_DECIDE, 400, 300, &(*vf)->lg);
495:   return 0;
496: }

498: PetscErrorCode KSPSetFromOptions_LSQR(KSP ksp, PetscOptionItems *PetscOptionsObject)
499: {
500:   KSP_LSQR *lsqr = (KSP_LSQR *)ksp->data;

502:   PetscOptionsHeadBegin(PetscOptionsObject, "KSP LSQR Options");
503:   PetscOptionsBool("-ksp_lsqr_compute_standard_error", "Set Standard Error Estimates of Solution", "KSPLSQRSetComputeStandardErrorVec", lsqr->se_flg, &lsqr->se_flg, NULL);
504:   PetscOptionsBool("-ksp_lsqr_exact_mat_norm", "Compute exact matrix norm instead of iteratively refined estimate", "KSPLSQRSetExactMatNorm", lsqr->exact_norm, &lsqr->exact_norm, NULL);
505:   KSPMonitorSetFromOptions(ksp, "-ksp_lsqr_monitor", "lsqr_residual", NULL);
506:   PetscOptionsHeadEnd();
507:   return 0;
508: }

510: PetscErrorCode KSPView_LSQR(KSP ksp, PetscViewer viewer)
511: {
512:   KSP_LSQR *lsqr = (KSP_LSQR *)ksp->data;
513:   PetscBool iascii;

515:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
516:   if (iascii) {
517:     if (lsqr->se) {
518:       PetscReal rnorm;
519:       VecNorm(lsqr->se, NORM_2, &rnorm);
520:       PetscViewerASCIIPrintf(viewer, "  norm of standard error %g, iterations %" PetscInt_FMT "\n", (double)rnorm, ksp->its);
521:     } else {
522:       PetscViewerASCIIPrintf(viewer, "  standard error not computed\n");
523:     }
524:     if (lsqr->exact_norm) {
525:       PetscViewerASCIIPrintf(viewer, "  using exact matrix norm\n");
526:     } else {
527:       PetscViewerASCIIPrintf(viewer, "  using inexact matrix norm\n");
528:     }
529:   }
530:   return 0;
531: }

533: /*@C
534:    KSPLSQRConvergedDefault - Determines convergence of the `KSPLSQR` Krylov method.

536:    Collective

538:    Input Parameters:
539: +  ksp   - iterative context
540: .  n     - iteration number
541: .  rnorm - 2-norm residual value (may be estimated)
542: -  ctx - convergence context which must be created by `KSPConvergedDefaultCreate()`

544:    reason is set to:
545: +   positive - if the iteration has converged;
546: .   negative - if residual norm exceeds divergence threshold;
547: -   0 - otherwise.

549:    Notes:
550:    `KSPConvergedDefault()` is called first to check for convergence in A*x=b.
551:    If that does not determine convergence then checks convergence for the least squares problem, i.e. in min{|b-A*x|}.
552:    Possible convergence for the least squares problem (which is based on the residual of the normal equations) are `KSP_CONVERGED_RTOL_NORMAL` norm
553:    and `KSP_CONVERGED_ATOL_NORMAL`.

555:    `KSP_CONVERGED_RTOL_NORMAL` is returned if ||A'*r|| < rtol * ||A|| * ||r||.
556:    Matrix norm ||A|| is iteratively refined estimate, see `KSPLSQRGetNorms()`.
557:    This criterion is now largely compatible with that in MATLAB `lsqr()`.

559:    Level: intermediate

561: .seealso: [](chapter_ksp), `KSPLSQR`, `KSPSetConvergenceTest()`, `KSPSetTolerances()`, `KSPConvergedSkip()`, `KSPConvergedReason`, `KSPGetConvergedReason()`,
562:           `KSPConvergedDefaultSetUIRNorm()`, `KSPConvergedDefaultSetUMIRNorm()`, `KSPConvergedDefaultCreate()`, `KSPConvergedDefaultDestroy()`, `KSPConvergedDefault()`, `KSPLSQRGetNorms()`, `KSPLSQRSetExactMatNorm()`
563: @*/
564: PetscErrorCode KSPLSQRConvergedDefault(KSP ksp, PetscInt n, PetscReal rnorm, KSPConvergedReason *reason, void *ctx)
565: {
566:   KSP_LSQR *lsqr = (KSP_LSQR *)ksp->data;

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

572:   /* check for convergence in min{|b-A*x|} */
573:   if (lsqr->arnorm < ksp->abstol) {
574:     PetscInfo(ksp, "LSQR solver has converged. Normal equation residual %14.12e is less than absolute tolerance %14.12e at iteration %" PetscInt_FMT "\n", (double)lsqr->arnorm, (double)ksp->abstol, n);
575:     *reason = KSP_CONVERGED_ATOL_NORMAL;
576:   } else if (lsqr->arnorm < ksp->rtol * lsqr->anorm * rnorm) {
577:     PetscCall(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 %" PetscInt_FMT "\n", (double)lsqr->arnorm,
578:                         (double)ksp->rtol, lsqr->exact_norm ? "exact" : "approx.", (double)lsqr->anorm, (double)rnorm, n));
579:     *reason = KSP_CONVERGED_RTOL_NORMAL;
580:   }
581:   return 0;
582: }

584: /*MC
585:      KSPLSQR - Implements LSQR

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

592:    Level: beginner

594:    Notes:
595:      Supports non-square (rectangular) matrices.

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

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

603:      Supports only left preconditioning.

605:      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()`.

607:      In exact arithmetic the LSQR method (with no preconditioning) is identical to the `KSPCG` algorithm applied to the normal equations.
608:      The preconditioned variant was implemented by Bas van't Hof and is essentially a left preconditioning for the Normal Equations.
609:      It appears the implementation with preconditioning tracks the true norm of the residual and uses that in the convergence test.

611:    Developer Note:
612:     How is this related to the `KSPCGNE` implementation? One difference is that `KSPCGNE` applies
613:     the preconditioner transpose times the preconditioner,  so one does not need to pass A'*A as the third argument to `KSPSetOperators()`.

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

618: .seealso: [](chapter_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPSolve()`, `KSPLSQRConvergedDefault()`, `KSPLSQRSetComputeStandardErrorVec()`, `KSPLSQRGetStandardErrorVec()`, `KSPLSQRSetExactMatNorm()`, `KSPLSQRMonitorResidualDrawLGCreate()`, `KSPLSQRMonitorResidualDrawLG()`, `KSPLSQRMonitorResidual()`
619: M*/
620: PETSC_EXTERN PetscErrorCode KSPCreate_LSQR(KSP ksp)
621: {
622:   KSP_LSQR *lsqr;
623:   void     *ctx;

625:   PetscNew(&lsqr);
626:   lsqr->se         = NULL;
627:   lsqr->se_flg     = PETSC_FALSE;
628:   lsqr->exact_norm = PETSC_FALSE;
629:   lsqr->anorm      = -1.0;
630:   lsqr->arnorm     = -1.0;
631:   ksp->data        = (void *)lsqr;
632:   KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 3);

634:   ksp->ops->setup          = KSPSetUp_LSQR;
635:   ksp->ops->solve          = KSPSolve_LSQR;
636:   ksp->ops->destroy        = KSPDestroy_LSQR;
637:   ksp->ops->setfromoptions = KSPSetFromOptions_LSQR;
638:   ksp->ops->view           = KSPView_LSQR;

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