Actual source code: cgls.c
1: /*
2: This file implements CGLS, the Conjugate Gradient method for Least-Squares problems.
3: */
4: #include <petsc/private/kspimpl.h>
6: typedef struct {
7: PetscInt nwork_n, nwork_m;
8: Vec *vwork_m; /* work vectors of length m, where the system is size m x n */
9: Vec *vwork_n; /* work vectors of length n */
10: } KSP_CGLS;
12: static PetscErrorCode KSPSetUp_CGLS(KSP ksp)
13: {
14: KSP_CGLS *cgls = (KSP_CGLS *)ksp->data;
16: cgls->nwork_m = 2;
17: if (cgls->vwork_m) VecDestroyVecs(cgls->nwork_m, &cgls->vwork_m);
19: cgls->nwork_n = 2;
20: if (cgls->vwork_n) VecDestroyVecs(cgls->nwork_n, &cgls->vwork_n);
21: KSPCreateVecs(ksp, cgls->nwork_n, &cgls->vwork_n, cgls->nwork_m, &cgls->vwork_m);
22: return 0;
23: }
25: static PetscErrorCode KSPSolve_CGLS(KSP ksp)
26: {
27: KSP_CGLS *cgls = (KSP_CGLS *)ksp->data;
28: Mat A;
29: Vec x, b, r, p, q, ss;
30: PetscScalar beta;
31: PetscReal alpha, gamma, oldgamma;
33: KSPGetOperators(ksp, &A, NULL); /* Matrix of the system */
35: /* vectors of length n, where system size is mxn */
36: x = ksp->vec_sol; /* Solution vector */
37: p = cgls->vwork_n[0];
38: ss = cgls->vwork_n[1];
40: /* vectors of length m, where system size is mxn */
41: b = ksp->vec_rhs; /* Right-hand side vector */
42: r = cgls->vwork_m[0];
43: q = cgls->vwork_m[1];
45: /* Minimization with the CGLS method */
46: ksp->its = 0;
47: ksp->rnorm = 0;
48: MatMult(A, x, r);
49: VecAYPX(r, -1, b); /* r_0 = b - A * x_0 */
50: KSP_MatMultHermitianTranspose(ksp, A, r, p); /* p_0 = A^T * r_0 */
51: VecCopy(p, ss); /* s_0 = p_0 */
52: VecNorm(ss, NORM_2, &gamma);
53: KSPCheckNorm(ksp, gamma);
54: if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = gamma;
55: (*ksp->converged)(ksp, ksp->its, ksp->rnorm, &ksp->reason, ksp->cnvP);
56: if (ksp->reason) return 0;
57: gamma = gamma * gamma; /* gamma = norm2(s)^2 */
59: do {
60: MatMult(A, p, q); /* q = A * p */
61: VecNorm(q, NORM_2, &alpha);
62: KSPCheckNorm(ksp, alpha);
63: alpha = alpha * alpha; /* alpha = norm2(q)^2 */
64: alpha = gamma / alpha; /* alpha = gamma / alpha */
65: VecAXPY(x, alpha, p); /* x += alpha * p */
66: VecAXPY(r, -alpha, q); /* r -= alpha * q */
67: KSP_MatMultHermitianTranspose(ksp, A, r, ss); /* ss = A^T * r */
68: oldgamma = gamma; /* oldgamma = gamma */
69: VecNorm(ss, NORM_2, &gamma);
70: KSPCheckNorm(ksp, gamma);
71: ksp->its++;
72: if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = gamma;
73: KSPMonitor(ksp, ksp->its, ksp->rnorm);
74: (*ksp->converged)(ksp, ksp->its, ksp->rnorm, &ksp->reason, ksp->cnvP);
75: if (ksp->reason) return 0;
76: gamma = gamma * gamma; /* gamma = norm2(s)^2 */
77: beta = gamma / oldgamma; /* beta = gamma / oldgamma */
78: VecAYPX(p, beta, ss); /* p = s + beta * p */
79: } while (ksp->its < ksp->max_it);
81: if (ksp->its >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
82: return 0;
83: }
85: static PetscErrorCode KSPDestroy_CGLS(KSP ksp)
86: {
87: KSP_CGLS *cgls = (KSP_CGLS *)ksp->data;
89: /* Free work vectors */
90: if (cgls->vwork_n) VecDestroyVecs(cgls->nwork_n, &cgls->vwork_n);
91: if (cgls->vwork_m) VecDestroyVecs(cgls->nwork_m, &cgls->vwork_m);
92: PetscFree(ksp->data);
93: return 0;
94: }
96: /*MC
97: KSPCGLS - Conjugate Gradient method for Least-Squares problems. Supports non-square (rectangular) matrices.
99: Level: beginner
101: Note:
102: This does not use the preconditioner, so one should probably use `KSPLSQR` instead.
104: .seealso: [](chapter_ksp), `KSPLSQR`, `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`,
105: `KSPCGSetType()`, `KSPCGUseSingleReduction()`, `KSPPIPECG`, `KSPGROPPCG`
106: M*/
107: PETSC_EXTERN PetscErrorCode KSPCreate_CGLS(KSP ksp)
108: {
109: KSP_CGLS *cgls;
111: PetscNew(&cgls);
112: ksp->data = (void *)cgls;
113: KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 3);
114: KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1);
115: ksp->ops->setup = KSPSetUp_CGLS;
116: ksp->ops->solve = KSPSolve_CGLS;
117: ksp->ops->destroy = KSPDestroy_CGLS;
118: ksp->ops->buildsolution = KSPBuildSolutionDefault;
119: ksp->ops->buildresidual = KSPBuildResidualDefault;
120: ksp->ops->setfromoptions = NULL;
121: ksp->ops->view = NULL;
122: return 0;
123: }