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