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:   PetscFunctionBegin;
 36:   PetscCall(VecGetLocalSize(v, &n));
 37:   PetscCall(VecGetArray(v, &x));
 38:   for (i = 0; i < n; i++) x[i] *= PetscConj(x[i]);
 39:   PetscCall(VecRestoreArray(v, &x));
 40:   PetscFunctionReturn(PETSC_SUCCESS);
 41: }

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

 48:   PetscFunctionBegin;
 49:   PetscCall(PetscObjectTypeCompare((PetscObject)ksp->pc, PCNONE, &nopreconditioner));

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

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

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

 60:   if (lsqr->se_flg && !lsqr->se) {
 61:     PetscCall(VecDuplicate(lsqr->vwork_n[0], &lsqr->se));
 62:     PetscCall(VecSet(lsqr->se, PETSC_INFINITY));
 63:   } else if (!lsqr->se_flg) {
 64:     PetscCall(VecDestroy(&lsqr->se));
 65:   }
 66:   PetscFunctionReturn(PETSC_SUCCESS);
 67: }

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

 79:   PetscFunctionBegin;
 80:   PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
 81:   PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);

 83:   PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));
 84:   PetscCall(PetscObjectTypeCompare((PetscObject)ksp->pc, PCNONE, &nopreconditioner));

 86:   /* vectors of length m, where system size is mxn */
 87:   B  = ksp->vec_rhs;
 88:   U  = lsqr->vwork_m[0];
 89:   U1 = lsqr->vwork_m[1];

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

 99:   /* standard error vector */
100:   if (lsqr->se) PetscCall(VecSet(lsqr->se, 0.0));

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

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

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

141:   if (nopreconditioner) {
142:     PetscCall(VecCopy(V, W));
143:   } else {
144:     PetscCall(VecCopy(Z, W));
145:   }

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

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

169:     PetscCall(KSP_MatMultHermitianTranspose(ksp, Amat, U1, V1));
170:     PetscCall(VecAXPY(V1, -beta, V));
171:     if (nopreconditioner) {
172:       PetscCall(VecNorm(V1, NORM_2, &alpha));
173:       KSPCheckNorm(ksp, alpha);
174:     } else {
175:       PetscCall(PCApply(ksp->pc, V1, Z));
176:       PetscCall(VecDotRealPart(V1, Z, &alpha));
177:       if (alpha <= 0.0) {
178:         ksp->reason = KSP_DIVERGED_BREAKDOWN;
179:         break;
180:       }
181:       alpha = PetscSqrtReal(alpha);
182:       PetscCall(VecScale(Z, 1.0 / alpha));
183:     }
184:     PetscCall(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:     PetscCall(VecAXPY(X, phi / rho, W)); /*    x <- x + (phi/rho) w   */

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

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

211:     PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
212:     ksp->its++;
213:     ksp->rnorm = rnorm;
214:     PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
215:     PetscCall(KSPLogResidualHistory(ksp, rnorm));
216:     PetscCall(KSPMonitor(ksp, i + 1, rnorm));
217:     PetscCall((*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:     PetscCall(MatGetSize(Amat, &size1, &size2));
230:     if (size1 > size2) tmp = size1 - size2;
231:     tmp = rnorm / PetscSqrtScalar(tmp);
232:     PetscCall(VecSqrtAbs(lsqr->se));
233:     PetscCall(VecScale(lsqr->se, tmp));
234:   }
235:   PetscFunctionReturn(PETSC_SUCCESS);
236: }

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

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

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

259:    Logically Collective

261:    Input Parameters:
262: +  ksp   - iterative context
263: -  flg   - compute the vector of standard estimates or not

265:    Level: intermediate

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

270: .seealso: [](ch_ksp), `KSPSolve()`, `KSPLSQR`, `KSPLSQRGetStandardErrorVec()`
271: @*/
272: PetscErrorCode KSPLSQRSetComputeStandardErrorVec(KSP ksp, PetscBool flg)
273: {
274:   KSP_LSQR *lsqr = (KSP_LSQR *)ksp->data;

276:   PetscFunctionBegin;
277:   lsqr->se_flg = flg;
278:   PetscFunctionReturn(PETSC_SUCCESS);
279: }

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

284:    Not Collective

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

290:    Level: intermediate

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: .seealso: [](ch_ksp), `KSPSolve()`, `KSPLSQR`, `KSPLSQRGetNorms()`, `KSPLSQRConvergedDefault()`
298: @*/
299: PetscErrorCode KSPLSQRSetExactMatNorm(KSP ksp, PetscBool flg)
300: {
301:   KSP_LSQR *lsqr = (KSP_LSQR *)ksp->data;

303:   PetscFunctionBegin;
304:   lsqr->exact_norm = flg;
305:   PetscFunctionReturn(PETSC_SUCCESS);
306: }

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

314:    Not Collective

316:    Input Parameter:
317: .  ksp   - iterative context

319:    Output Parameter:
320: .  se - vector of standard estimates

322:    Level: intermediate

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

327: .seealso: [](ch_ksp), `KSPSolve()`, `KSPLSQR`, `KSPLSQRSetComputeStandardErrorVec()`
328: @*/
329: PetscErrorCode KSPLSQRGetStandardErrorVec(KSP ksp, Vec *se)
330: {
331:   KSP_LSQR *lsqr = (KSP_LSQR *)ksp->data;

333:   PetscFunctionBegin;
334:   *se = lsqr->se;
335:   PetscFunctionReturn(PETSC_SUCCESS);
336: }

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

341:    Not Collective

343:    Input Parameter:
344: .  ksp   - iterative context

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

350:    Notes:
351:    Output parameters are meaningful only after `KSPSolve()`.

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

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

357:    Level: intermediate

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

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

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

401:   Collective

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

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

412:   Level: intermediate

414: .seealso: [](ch_ksp), `KSPLSQR`, `KSPMonitorSet()`, `KSPMonitorResidual()`, `KSPMonitorTrueResidualMaxNorm()`, `KSPLSQRMonitorResidualDrawLG()`
415: @*/
416: PetscErrorCode KSPLSQRMonitorResidual(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
417: {
418:   PetscFunctionBegin;
422:   PetscTryMethod(ksp, "KSPLSQRMonitorResidual_C", (KSP, PetscInt, PetscReal, PetscViewerAndFormat *), (ksp, n, rnorm, vf));
423:   PetscFunctionReturn(PETSC_SUCCESS);
424: }

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

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

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

457:   Collective

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

465:   Options Database Key:
466: . -ksp_lsqr_monitor draw::draw_lg - Activates `KSPMonitorTrueResidualDrawLG()`

468:   Level: intermediate

470: .seealso: [](ch_ksp), `KSPLSQR`, `KSPMonitorSet()`, `KSPMonitorTrueResidual()`, `KSPLSQRMonitorResidual()`, `KSPLSQRMonitorResidualDrawLGCreate()`
471: @*/
472: PetscErrorCode KSPLSQRMonitorResidualDrawLG(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
473: {
474:   PetscFunctionBegin;
479:   PetscTryMethod(ksp, "KSPLSQRMonitorResidualDrawLG_C", (KSP, PetscInt, PetscReal, PetscViewerAndFormat *), (ksp, n, rnorm, vf));
480:   PetscFunctionReturn(PETSC_SUCCESS);
481: }

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

486:   Collective

488:   Input Parameters:
489: + viewer - The PetscViewer
490: . format - The viewer format
491: - ctx    - An optional user context

493:   Output Parameter:
494: . vf    - The viewer context

496:   Level: intermediate

498: .seealso: [](ch_ksp), `KSPLSQR`, `KSPMonitorSet()`, `KSPLSQRMonitorResidual()`, `KSPLSQRMonitorResidualDrawLG()`
499: @*/
500: PetscErrorCode KSPLSQRMonitorResidualDrawLGCreate(PetscViewer viewer, PetscViewerFormat format, void *ctx, PetscViewerAndFormat **vf)
501: {
502:   const char *names[] = {"residual", "normal eqn residual"};

504:   PetscFunctionBegin;
505:   PetscCall(PetscViewerAndFormatCreate(viewer, format, vf));
506:   (*vf)->data = ctx;
507:   PetscCall(KSPMonitorLGCreate(PetscObjectComm((PetscObject)viewer), NULL, NULL, "Log Residual Norm", 2, names, PETSC_DECIDE, PETSC_DECIDE, 400, 300, &(*vf)->lg));
508:   PetscFunctionReturn(PETSC_SUCCESS);
509: }

511: PetscErrorCode KSPSetFromOptions_LSQR(KSP ksp, PetscOptionItems *PetscOptionsObject)
512: {
513:   KSP_LSQR *lsqr = (KSP_LSQR *)ksp->data;

515:   PetscFunctionBegin;
516:   PetscOptionsHeadBegin(PetscOptionsObject, "KSP LSQR Options");
517:   PetscCall(PetscOptionsBool("-ksp_lsqr_compute_standard_error", "Set Standard Error Estimates of Solution", "KSPLSQRSetComputeStandardErrorVec", lsqr->se_flg, &lsqr->se_flg, NULL));
518:   PetscCall(PetscOptionsBool("-ksp_lsqr_exact_mat_norm", "Compute exact matrix norm instead of iteratively refined estimate", "KSPLSQRSetExactMatNorm", lsqr->exact_norm, &lsqr->exact_norm, NULL));
519:   PetscCall(KSPMonitorSetFromOptions(ksp, "-ksp_lsqr_monitor", "lsqr_residual", NULL));
520:   PetscOptionsHeadEnd();
521:   PetscFunctionReturn(PETSC_SUCCESS);
522: }

524: PetscErrorCode KSPView_LSQR(KSP ksp, PetscViewer viewer)
525: {
526:   KSP_LSQR *lsqr = (KSP_LSQR *)ksp->data;
527:   PetscBool iascii;

529:   PetscFunctionBegin;
530:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
531:   if (iascii) {
532:     if (lsqr->se) {
533:       PetscReal rnorm;
534:       PetscCall(VecNorm(lsqr->se, NORM_2, &rnorm));
535:       PetscCall(PetscViewerASCIIPrintf(viewer, "  norm of standard error %g, iterations %" PetscInt_FMT "\n", (double)rnorm, ksp->its));
536:     } else {
537:       PetscCall(PetscViewerASCIIPrintf(viewer, "  standard error not computed\n"));
538:     }
539:     if (lsqr->exact_norm) {
540:       PetscCall(PetscViewerASCIIPrintf(viewer, "  using exact matrix norm\n"));
541:     } else {
542:       PetscCall(PetscViewerASCIIPrintf(viewer, "  using inexact matrix norm\n"));
543:     }
544:   }
545:   PetscFunctionReturn(PETSC_SUCCESS);
546: }

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

551:    Collective

553:    Input Parameters:
554: +  ksp   - iterative context
555: .  n     - iteration number
556: .  rnorm - 2-norm residual value (may be estimated)
557: -  ctx - convergence context which must be created by `KSPConvergedDefaultCreate()`

559:    reason is set to:
560: +   positive - if the iteration has converged;
561: .   negative - if residual norm exceeds divergence threshold;
562: -   0 - otherwise.

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

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

574:    Level: intermediate

576: .seealso: [](ch_ksp), `KSPLSQR`, `KSPSetConvergenceTest()`, `KSPSetTolerances()`, `KSPConvergedSkip()`, `KSPConvergedReason`, `KSPGetConvergedReason()`,
577:           `KSPConvergedDefaultSetUIRNorm()`, `KSPConvergedDefaultSetUMIRNorm()`, `KSPConvergedDefaultCreate()`, `KSPConvergedDefaultDestroy()`, `KSPConvergedDefault()`, `KSPLSQRGetNorms()`, `KSPLSQRSetExactMatNorm()`
578: @*/
579: PetscErrorCode KSPLSQRConvergedDefault(KSP ksp, PetscInt n, PetscReal rnorm, KSPConvergedReason *reason, void *ctx)
580: {
581:   KSP_LSQR *lsqr = (KSP_LSQR *)ksp->data;

583:   PetscFunctionBegin;
584:   /* check for convergence in A*x=b */
585:   PetscCall(KSPConvergedDefault(ksp, n, rnorm, reason, ctx));
586:   if (!n || *reason) PetscFunctionReturn(PETSC_SUCCESS);

588:   /* check for convergence in min{|b-A*x|} */
589:   if (lsqr->arnorm < ksp->abstol) {
590:     PetscCall(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));
591:     *reason = KSP_CONVERGED_ATOL_NORMAL;
592:   } else if (lsqr->arnorm < ksp->rtol * lsqr->anorm * rnorm) {
593:     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,
594:                         (double)ksp->rtol, lsqr->exact_norm ? "exact" : "approx.", (double)lsqr->anorm, (double)rnorm, n));
595:     *reason = KSP_CONVERGED_RTOL_NORMAL;
596:   }
597:   PetscFunctionReturn(PETSC_SUCCESS);
598: }

600: /*MC
601:      KSPLSQR - Implements LSQR

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

608:    Level: beginner

610:    Notes:
611:      Supports non-square (rectangular) matrices.

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

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

619:      Supports only left preconditioning.

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

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

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

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

634: .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPSolve()`, `KSPLSQRConvergedDefault()`, `KSPLSQRSetComputeStandardErrorVec()`, `KSPLSQRGetStandardErrorVec()`, `KSPLSQRSetExactMatNorm()`, `KSPLSQRMonitorResidualDrawLGCreate()`, `KSPLSQRMonitorResidualDrawLG()`, `KSPLSQRMonitorResidual()`
635: M*/
636: PETSC_EXTERN PetscErrorCode KSPCreate_LSQR(KSP ksp)
637: {
638:   KSP_LSQR *lsqr;
639:   void     *ctx;

641:   PetscFunctionBegin;
642:   PetscCall(PetscNew(&lsqr));
643:   lsqr->se         = NULL;
644:   lsqr->se_flg     = PETSC_FALSE;
645:   lsqr->exact_norm = PETSC_FALSE;
646:   lsqr->anorm      = -1.0;
647:   lsqr->arnorm     = -1.0;
648:   ksp->data        = (void *)lsqr;
649:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 3));

651:   ksp->ops->setup          = KSPSetUp_LSQR;
652:   ksp->ops->solve          = KSPSolve_LSQR;
653:   ksp->ops->destroy        = KSPDestroy_LSQR;
654:   ksp->ops->setfromoptions = KSPSetFromOptions_LSQR;
655:   ksp->ops->view           = KSPView_LSQR;

657:   /* Backup current convergence test; remove destroy routine from KSP to prevent destroying the convergence context in KSPSetConvergenceTest() */
658:   PetscCall(KSPGetAndClearConvergenceTest(ksp, &lsqr->converged, &lsqr->cnvP, &lsqr->convergeddestroy));
659:   /* Override current convergence test */
660:   PetscCall(KSPConvergedDefaultCreate(&ctx));
661:   PetscCall(KSPSetConvergenceTest(ksp, KSPLSQRConvergedDefault, ctx, KSPConvergedDefaultDestroy));
662:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPLSQRMonitorResidual_C", KSPLSQRMonitorResidual_LSQR));
663:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPLSQRMonitorResidualDrawLG_C", KSPLSQRMonitorResidualDrawLG_LSQR));
664:   PetscFunctionReturn(PETSC_SUCCESS);
665: }