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