1 /*
2 cgimpl.h defines the simple data structured used to store information
3 related to the type of matrix (e.g. complex symmetric) being solved and
4 data used during the optional Lanczos process used to compute eigenvalues
5 */
6 #include <../src/ksp/ksp/impls/cg/cgimpl.h> /*I "petscksp.h" I*/
7 extern PetscErrorCode KSPComputeExtremeSingularValues_CG(KSP, PetscReal *, PetscReal *);
8 extern PetscErrorCode KSPComputeEigenvalues_CG(KSP, PetscInt, PetscReal *, PetscReal *, PetscInt *);
9
KSPCGSetType_CGNE(KSP ksp,KSPCGType type)10 static PetscErrorCode KSPCGSetType_CGNE(KSP ksp, KSPCGType type)
11 {
12 KSP_CG *cg = (KSP_CG *)ksp->data;
13
14 PetscFunctionBegin;
15 cg->type = type;
16 PetscFunctionReturn(PETSC_SUCCESS);
17 }
18
KSPSetUp_CGNE(KSP ksp)19 static PetscErrorCode KSPSetUp_CGNE(KSP ksp)
20 {
21 KSP_CG *cgP = (KSP_CG *)ksp->data;
22 PetscInt maxit = ksp->max_it;
23
24 PetscFunctionBegin;
25 /* get work vectors needed by CGNE */
26 PetscCall(KSPSetWorkVecs(ksp, 4));
27
28 /*
29 If user requested computations of eigenvalues then allocate work space needed
30 */
31 if (ksp->calc_sings) {
32 /* get space to store tridiagonal matrix for Lanczos */
33 PetscCall(PetscMalloc4(maxit, &cgP->e, maxit, &cgP->d, maxit, &cgP->ee, maxit, &cgP->dd));
34
35 ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_CG;
36 ksp->ops->computeeigenvalues = KSPComputeEigenvalues_CG;
37 }
38 PetscFunctionReturn(PETSC_SUCCESS);
39 }
40
KSPSolve_CGNE(KSP ksp)41 static PetscErrorCode KSPSolve_CGNE(KSP ksp)
42 {
43 PetscInt i, stored_max_it, eigs;
44 PetscScalar dpi, a = 1.0, beta, betaold = 1.0, b = 0, *e = NULL, *d = NULL;
45 PetscReal dp = 0.0;
46 Vec X, B, Z, R, P, T;
47 KSP_CG *cg;
48 Mat Amat, Pmat;
49 PetscBool diagonalscale, transpose_pc;
50
51 PetscFunctionBegin;
52 PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
53 PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
54 PetscCall(PCApplyTransposeExists(ksp->pc, &transpose_pc));
55
56 cg = (KSP_CG *)ksp->data;
57 eigs = ksp->calc_sings;
58 stored_max_it = ksp->max_it;
59 X = ksp->vec_sol;
60 B = ksp->vec_rhs;
61 R = ksp->work[0];
62 Z = ksp->work[1];
63 P = ksp->work[2];
64 T = ksp->work[3];
65
66 #define VecXDot(x, y, a) (cg->type == KSP_CG_HERMITIAN ? VecDot(x, y, a) : VecTDot(x, y, a))
67
68 if (eigs) {
69 e = cg->e;
70 d = cg->d;
71 e[0] = 0.0;
72 }
73 PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));
74
75 ksp->its = 0;
76 PetscCall(KSP_MatMultTranspose(ksp, Amat, B, T));
77 if (!ksp->guess_zero) {
78 PetscCall(KSP_MatMult(ksp, Amat, X, P));
79 PetscCall(KSP_MatMultTranspose(ksp, Amat, P, R));
80 PetscCall(VecAYPX(R, -1.0, T));
81 } else {
82 PetscCall(VecCopy(T, R)); /* r <- b (x is 0) */
83 }
84 if (transpose_pc) {
85 PetscCall(KSP_PCApplyTranspose(ksp, R, T));
86 } else {
87 PetscCall(KSP_PCApply(ksp, R, T));
88 }
89 PetscCall(KSP_PCApply(ksp, T, Z));
90
91 if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
92 PetscCall(VecNorm(Z, NORM_2, &dp)); /* dp <- z'*z */
93 } else if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
94 PetscCall(VecNorm(R, NORM_2, &dp)); /* dp <- r'*r */
95 } else if (ksp->normtype == KSP_NORM_NATURAL) {
96 PetscCall(VecXDot(Z, R, &beta));
97 KSPCheckDot(ksp, beta);
98 dp = PetscSqrtReal(PetscAbsScalar(beta));
99 } else dp = 0.0;
100 PetscCall(KSPLogResidualHistory(ksp, dp));
101 PetscCall(KSPMonitor(ksp, 0, dp));
102 ksp->rnorm = dp;
103 PetscCall((*ksp->converged)(ksp, 0, dp, &ksp->reason, ksp->cnvP)); /* test for convergence */
104 if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);
105
106 i = 0;
107 do {
108 ksp->its = i + 1;
109 PetscCall(VecXDot(Z, R, &beta)); /* beta <- r'z */
110 KSPCheckDot(ksp, beta);
111 if (beta == 0.0) {
112 ksp->reason = KSP_CONVERGED_ATOL;
113 PetscCall(PetscInfo(ksp, "converged due to beta = 0\n"));
114 break;
115 #if !defined(PETSC_USE_COMPLEX)
116 } else if (beta < 0.0) {
117 ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
118 PetscCall(PetscInfo(ksp, "diverging due to indefinite preconditioner\n"));
119 break;
120 #endif
121 }
122 if (!i) {
123 PetscCall(VecCopy(Z, P)); /* p <- z */
124 b = 0.0;
125 } else {
126 b = beta / betaold;
127 if (eigs) {
128 PetscCheck(ksp->max_it == stored_max_it, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Can not change maxit AND calculate eigenvalues");
129 e[i] = PetscSqrtReal(PetscAbsScalar(b)) / a;
130 }
131 PetscCall(VecAYPX(P, b, Z)); /* p <- z + b* p */
132 }
133 betaold = beta;
134 PetscCall(KSP_MatMult(ksp, Amat, P, T));
135 PetscCall(KSP_MatMultTranspose(ksp, Amat, T, Z));
136 PetscCall(VecXDot(P, Z, &dpi)); /* dpi <- z'p */
137 KSPCheckDot(ksp, dpi);
138 a = beta / dpi; /* a = beta/p'z */
139 if (eigs) d[i] = PetscSqrtReal(PetscAbsScalar(b)) * e[i] + 1.0 / a;
140 PetscCall(VecAXPY(X, a, P)); /* x <- x + ap */
141 PetscCall(VecAXPY(R, -a, Z)); /* r <- r - az */
142 if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
143 if (transpose_pc) {
144 PetscCall(KSP_PCApplyTranspose(ksp, R, T));
145 } else {
146 PetscCall(KSP_PCApply(ksp, R, T));
147 }
148 PetscCall(KSP_PCApply(ksp, T, Z));
149 PetscCall(VecNorm(Z, NORM_2, &dp)); /* dp <- z'*z */
150 } else if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
151 PetscCall(VecNorm(R, NORM_2, &dp));
152 } else if (ksp->normtype == KSP_NORM_NATURAL) {
153 dp = PetscSqrtReal(PetscAbsScalar(beta));
154 } else dp = 0.0;
155 ksp->rnorm = dp;
156 PetscCall(KSPLogResidualHistory(ksp, dp));
157 PetscCall(KSPMonitor(ksp, i + 1, dp));
158 PetscCall((*ksp->converged)(ksp, i + 1, dp, &ksp->reason, ksp->cnvP));
159 if (ksp->reason) break;
160 if (ksp->normtype != KSP_NORM_PRECONDITIONED) {
161 if (transpose_pc) {
162 PetscCall(KSP_PCApplyTranspose(ksp, R, T));
163 } else {
164 PetscCall(KSP_PCApply(ksp, R, T));
165 }
166 PetscCall(KSP_PCApply(ksp, T, Z));
167 }
168 i++;
169 } while (i < ksp->max_it);
170 if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
171 PetscFunctionReturn(PETSC_SUCCESS);
172 }
173
174 /*MC
175 KSPCGNE - Applies the preconditioned conjugate gradient method to the normal equations
176 without explicitly forming $A^T A$.
177
178 Options Database Key:
179 . -ksp_cg_type <Hermitian or symmetric - (for complex matrices only) indicates the matrix is Hermitian or symmetric
180
181 Level: beginner
182
183 Notes:
184 Eigenvalue computation routines including `KSPSetComputeEigenvalues()` and `KSPComputeEigenvalues()` will return information about the
185 spectrum of $A^T*A$, rather than $A$.
186
187 `KSPCGNE` is a general-purpose non-symmetric method. It works well when the singular values are much better behaved than
188 eigenvalues. A unitary matrix is a classic example where `KSPCGNE` converges in one iteration, but `KSPGMRES` and `KSPCGS` need N
189 iterations, see {cite}`nachtigal90`. If you intend to solve least squares problems, use `KSPLSQR`.
190
191 This is NOT a different algorithm than used with `KSPCG`, it merely uses that algorithm with the
192 matrix defined by $A^T A$ and preconditioner defined by $B^T B$ where $B$ is the preconditioner for $A$.
193
194 See `PETSCREGRESSORLINEAR` for the PETSc toolkit for solving linear regression problems including least squares.
195
196 This method requires that one be able to apply the transpose of the preconditioner and operator
197 as well as the operator and preconditioner. If the transpose of the preconditioner is not available then
198 the preconditioner is used in its place so one ends up preconditioning $A^T A$ with $B B$. Seems odd?
199
200 This only supports left preconditioning.
201
202 Developer Note:
203 This object is subclassed off of `KSPCG`, see the source code in src/ksp/ksp/impls/cg for comments on the structure of the code
204
205 .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPCG`, `KSPLSQR`, `KSPCGLS`,
206 `KSPCGSetType()`, `KSPBICG`, `KSPSetComputeEigenvalues()`, `KSPComputeEigenvalues()`, `PETSCREGRESSORLINEAR`
207 M*/
208
KSPCreate_CGNE(KSP ksp)209 PETSC_EXTERN PetscErrorCode KSPCreate_CGNE(KSP ksp)
210 {
211 KSP_CG *cg;
212
213 PetscFunctionBegin;
214 PetscCall(PetscNew(&cg));
215 cg->type = !PetscDefined(USE_COMPLEX) ? KSP_CG_SYMMETRIC : KSP_CG_HERMITIAN;
216 ksp->data = (void *)cg;
217 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
218 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 2));
219 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 2));
220 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
221
222 /*
223 Sets the functions that are associated with this data structure
224 (in C++ this is the same as defining virtual functions)
225 */
226 ksp->ops->setup = KSPSetUp_CGNE;
227 ksp->ops->solve = KSPSolve_CGNE;
228 ksp->ops->destroy = KSPDestroy_CG;
229 ksp->ops->view = KSPView_CG;
230 ksp->ops->setfromoptions = KSPSetFromOptions_CG;
231 ksp->ops->buildsolution = KSPBuildSolutionDefault;
232 ksp->ops->buildresidual = KSPBuildResidualDefault;
233
234 /*
235 Attach the function KSPCGSetType_CGNE() to this object. The routine
236 KSPCGSetType() checks for this attached function and calls it if it finds
237 it. (Sort of like a dynamic member function that can be added at run time
238 */
239 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetType_C", KSPCGSetType_CGNE));
240 PetscFunctionReturn(PETSC_SUCCESS);
241 }
242