1 #include <../src/ksp/ksp/impls/lcd/lcdimpl.h>
2
KSPSetUp_LCD(KSP ksp)3 static PetscErrorCode KSPSetUp_LCD(KSP ksp)
4 {
5 KSP_LCD *lcd = (KSP_LCD *)ksp->data;
6 PetscInt restart = lcd->restart;
7
8 PetscFunctionBegin;
9 /* get work vectors needed by LCD */
10 PetscCall(KSPSetWorkVecs(ksp, 2));
11
12 PetscCall(VecDuplicateVecs(ksp->work[0], restart + 1, &lcd->P));
13 PetscCall(VecDuplicateVecs(ksp->work[0], restart + 1, &lcd->Q));
14 PetscFunctionReturn(PETSC_SUCCESS);
15 }
16
17 /* KSPSolve_LCD - This routine actually applies the left conjugate
18 direction method
19
20 Input Parameter:
21 . ksp - the Krylov space object that was set to use LCD, by, for
22 example, KSPCreate(MPI_Comm,KSP *ksp); KSPSetType(ksp,KSPLCD);
23
24 Output Parameter:
25 . its - number of iterations used
26
27 */
KSPSolve_LCD(KSP ksp)28 static PetscErrorCode KSPSolve_LCD(KSP ksp)
29 {
30 PetscInt it, j, max_k;
31 PetscScalar alfa, beta, num, den, mone;
32 PetscReal rnorm = 0.0;
33 Vec X, B, R, Z;
34 KSP_LCD *lcd;
35 Mat Amat, Pmat;
36 PetscBool diagonalscale;
37
38 PetscFunctionBegin;
39 PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
40 PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
41
42 lcd = (KSP_LCD *)ksp->data;
43 X = ksp->vec_sol;
44 B = ksp->vec_rhs;
45 R = ksp->work[0];
46 Z = ksp->work[1];
47 max_k = lcd->restart;
48 mone = -1;
49
50 PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));
51
52 ksp->its = 0;
53 if (!ksp->guess_zero) {
54 PetscCall(KSP_MatMult(ksp, Amat, X, Z)); /* z <- b - Ax */
55 PetscCall(VecAYPX(Z, mone, B));
56 } else {
57 PetscCall(VecCopy(B, Z)); /* z <- b (x is 0) */
58 }
59
60 PetscCall(KSP_PCApply(ksp, Z, R)); /* r <- M^-1z */
61 if (ksp->normtype != KSP_NORM_NONE) {
62 PetscCall(VecNorm(R, NORM_2, &rnorm));
63 KSPCheckNorm(ksp, rnorm);
64 }
65 PetscCall(KSPLogResidualHistory(ksp, rnorm));
66 PetscCall(KSPMonitor(ksp, 0, rnorm));
67 ksp->rnorm = rnorm;
68
69 /* test for convergence */
70 PetscCall((*ksp->converged)(ksp, 0, rnorm, &ksp->reason, ksp->cnvP));
71 if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);
72
73 PetscCall(VecCopy(R, lcd->P[0]));
74
75 while (!ksp->reason && ksp->its < ksp->max_it) {
76 it = 0;
77 PetscCall(KSP_MatMult(ksp, Amat, lcd->P[it], Z));
78 PetscCall(KSP_PCApply(ksp, Z, lcd->Q[it]));
79
80 while (!ksp->reason && it < max_k && ksp->its < ksp->max_it) {
81 ksp->its++;
82 PetscCall(VecDot(lcd->P[it], R, &num));
83 PetscCall(VecDot(lcd->P[it], lcd->Q[it], &den));
84 KSPCheckDot(ksp, den);
85 alfa = num / den;
86 PetscCall(VecAXPY(X, alfa, lcd->P[it]));
87 PetscCall(VecAXPY(R, -alfa, lcd->Q[it]));
88 if (ksp->normtype != KSP_NORM_NONE) {
89 PetscCall(VecNorm(R, NORM_2, &rnorm));
90 KSPCheckNorm(ksp, rnorm);
91 }
92
93 ksp->rnorm = rnorm;
94 PetscCall(KSPLogResidualHistory(ksp, rnorm));
95 PetscCall(KSPMonitor(ksp, ksp->its, rnorm));
96 PetscCall((*ksp->converged)(ksp, ksp->its, rnorm, &ksp->reason, ksp->cnvP));
97
98 if (ksp->reason) break;
99
100 PetscCall(VecCopy(R, lcd->P[it + 1]));
101 PetscCall(KSP_MatMult(ksp, Amat, lcd->P[it + 1], Z));
102 PetscCall(KSP_PCApply(ksp, Z, lcd->Q[it + 1]));
103
104 for (j = 0; j <= it; j++) {
105 PetscCall(VecDot(lcd->P[j], lcd->Q[it + 1], &num));
106 KSPCheckDot(ksp, num);
107 PetscCall(VecDot(lcd->P[j], lcd->Q[j], &den));
108 beta = -num / den;
109 PetscCall(VecAXPY(lcd->P[it + 1], beta, lcd->P[j]));
110 PetscCall(VecAXPY(lcd->Q[it + 1], beta, lcd->Q[j]));
111 }
112 it++;
113 }
114 PetscCall(VecCopy(lcd->P[it], lcd->P[0]));
115 }
116 if (ksp->its >= ksp->max_it && !ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
117 PetscCall(VecCopy(X, ksp->vec_sol));
118 PetscFunctionReturn(PETSC_SUCCESS);
119 }
120 /*
121 KSPDestroy_LCD - Frees all memory space used by the Krylov method
122
123 */
KSPReset_LCD(KSP ksp)124 static PetscErrorCode KSPReset_LCD(KSP ksp)
125 {
126 KSP_LCD *lcd = (KSP_LCD *)ksp->data;
127
128 PetscFunctionBegin;
129 if (lcd->P) PetscCall(VecDestroyVecs(lcd->restart + 1, &lcd->P));
130 if (lcd->Q) PetscCall(VecDestroyVecs(lcd->restart + 1, &lcd->Q));
131 PetscFunctionReturn(PETSC_SUCCESS);
132 }
133
KSPDestroy_LCD(KSP ksp)134 static PetscErrorCode KSPDestroy_LCD(KSP ksp)
135 {
136 PetscFunctionBegin;
137 PetscCall(KSPReset_LCD(ksp));
138 PetscCall(PetscFree(ksp->data));
139 PetscFunctionReturn(PETSC_SUCCESS);
140 }
141
142 /*
143 KSPView_LCD - Prints information about the current Krylov method being used
144
145 Currently this only prints information to a file (or stdout) about the
146 symmetry of the problem. If your Krylov method has special options or
147 flags that information should be printed here.
148
149 */
KSPView_LCD(KSP ksp,PetscViewer viewer)150 static PetscErrorCode KSPView_LCD(KSP ksp, PetscViewer viewer)
151 {
152 KSP_LCD *lcd = (KSP_LCD *)ksp->data;
153 PetscBool isascii;
154
155 PetscFunctionBegin;
156 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
157 if (isascii) {
158 PetscCall(PetscViewerASCIIPrintf(viewer, " restart=%" PetscInt_FMT "\n", lcd->restart));
159 PetscCall(PetscViewerASCIIPrintf(viewer, " happy breakdown tolerance=%g\n", (double)lcd->haptol));
160 }
161 PetscFunctionReturn(PETSC_SUCCESS);
162 }
163
164 /*
165 KSPSetFromOptions_LCD - Checks the options database for options related to the
166 LCD method.
167 */
KSPSetFromOptions_LCD(KSP ksp,PetscOptionItems PetscOptionsObject)168 static PetscErrorCode KSPSetFromOptions_LCD(KSP ksp, PetscOptionItems PetscOptionsObject)
169 {
170 PetscBool flg;
171 KSP_LCD *lcd = (KSP_LCD *)ksp->data;
172
173 PetscFunctionBegin;
174 PetscOptionsHeadBegin(PetscOptionsObject, "KSP LCD options");
175 PetscCall(PetscOptionsBoundedInt("-ksp_lcd_restart", "Number of vectors conjugate", "KSPLCDSetRestart", lcd->restart, &lcd->restart, &flg, 1));
176 PetscCall(PetscOptionsBoundedReal("-ksp_lcd_haptol", "Tolerance for exact convergence (happy ending)", "KSPLCDSetHapTol", lcd->haptol, &lcd->haptol, &flg, 0.0));
177 PetscFunctionReturn(PETSC_SUCCESS);
178 }
179
180 /*MC
181 KSPLCD - Implements the LCD (left conjugate direction) method
182
183 Options Database Keys:
184 + -ksp_lcd_restart - number of vectors conjugate
185 - -ksp_lcd_haptol - tolerance for exact convergence (happy ending)
186
187 Level: beginner
188
189 Notes:
190 Support only for left preconditioning
191
192 See {cite}`yuan2004semi`, {cite}`dai2004study`, {cite}`catabriga2004evaluating`, and {cite}`catabriga2006performance`
193
194 Contributed by:
195 Lucia Catabriga <luciac@ices.utexas.edu>
196
197 .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPCG`,
198 `KSPCGSetType()`, `KSPLCDSetRestart()`, `KSPLCDSetHapTol()`
199 M*/
200
KSPCreate_LCD(KSP ksp)201 PETSC_EXTERN PetscErrorCode KSPCreate_LCD(KSP ksp)
202 {
203 KSP_LCD *lcd;
204
205 PetscFunctionBegin;
206 PetscCall(PetscNew(&lcd));
207 ksp->data = (void *)lcd;
208 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
209 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
210 lcd->restart = 30;
211 lcd->haptol = 1.0e-30;
212
213 /*
214 Sets the functions that are associated with this data structure
215 (in C++ this is the same as defining virtual functions)
216 */
217 ksp->ops->setup = KSPSetUp_LCD;
218 ksp->ops->solve = KSPSolve_LCD;
219 ksp->ops->reset = KSPReset_LCD;
220 ksp->ops->destroy = KSPDestroy_LCD;
221 ksp->ops->view = KSPView_LCD;
222 ksp->ops->setfromoptions = KSPSetFromOptions_LCD;
223 ksp->ops->buildsolution = KSPBuildSolutionDefault;
224 ksp->ops->buildresidual = KSPBuildResidualDefault;
225 PetscFunctionReturn(PETSC_SUCCESS);
226 }
227