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: PetscFunctionBegin;
17: cgls->nwork_m = 2;
18: if (cgls->vwork_m) PetscCall(VecDestroyVecs(cgls->nwork_m, &cgls->vwork_m));
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: }
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;
34: PetscFunctionBegin;
35: PetscCall(KSPGetOperators(ksp, &A, NULL)); /* Matrix of the system */
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];
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];
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 */
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);
83: if (ksp->its >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
84: PetscFunctionReturn(PETSC_SUCCESS);
85: }
87: static PetscErrorCode KSPDestroy_CGLS(KSP ksp)
88: {
89: KSP_CGLS *cgls = (KSP_CGLS *)ksp->data;
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: }
99: /*MC
100: KSPCGLS - Conjugate Gradient method for Least-Squares problems. Supports non-square (rectangular) matrices.
102: Level: beginner
104: Note:
105: This does not use the preconditioner, so one should probably use `KSPLSQR` instead.
107: .seealso: [](ch_ksp), `KSPLSQR`, `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`,
108: `KSPCGSetType()`, `KSPCGUseSingleReduction()`, `KSPPIPECG`, `KSPGROPPCG`
109: M*/
110: PETSC_EXTERN PetscErrorCode KSPCreate_CGLS(KSP ksp)
111: {
112: KSP_CGLS *cgls;
114: PetscFunctionBegin;
115: PetscCall(PetscNew(&cgls));
116: ksp->data = (void *)cgls;
117: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 3));
118: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
119: ksp->ops->setup = KSPSetUp_CGLS;
120: ksp->ops->solve = KSPSolve_CGLS;
121: ksp->ops->destroy = KSPDestroy_CGLS;
122: ksp->ops->buildsolution = KSPBuildSolutionDefault;
123: ksp->ops->buildresidual = KSPBuildResidualDefault;
124: ksp->ops->setfromoptions = NULL;
125: ksp->ops->view = NULL;
126: PetscFunctionReturn(PETSC_SUCCESS);
127: }