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