1 #include <petsc/private/kspimpl.h> /*I "petscksp.h" I*/
2 #include <petscblaslapack.h>
3
4 typedef struct {
5 PetscInt method; /* 1, 2 or 3 */
6 PetscInt curl; /* Current number of basis vectors */
7 PetscInt maxl; /* Maximum number of basis vectors */
8 PetscBool monitor;
9 PetscScalar *alpha; /* */
10 Vec *xtilde; /* Saved x vectors */
11 Vec *btilde; /* Saved b vectors, methods 1 and 3 */
12 Vec Ax; /* method 2 */
13 Vec guess;
14 PetscScalar *corr; /* correlation matrix in column-major format, method 3 */
15 PetscReal tol; /* tolerance for determining rank, method 3 */
16 Vec last_b; /* last b provided to FormGuess (not owned by this object), method 3 */
17 PetscObjectState last_b_state; /* state of last_b as of the last call to FormGuess, method 3 */
18 PetscScalar *last_b_coefs; /* dot products of last_b and btilde, method 3 */
19 } KSPGuessFischer;
20
KSPGuessReset_Fischer(KSPGuess guess)21 static PetscErrorCode KSPGuessReset_Fischer(KSPGuess guess)
22 {
23 KSPGuessFischer *itg = (KSPGuessFischer *)guess->data;
24 PetscLayout Alay = NULL, vlay = NULL;
25 PetscBool cong;
26
27 PetscFunctionBegin;
28 itg->curl = 0;
29 /* destroy vectors if the size of the linear system has changed */
30 if (guess->A) PetscCall(MatGetLayouts(guess->A, &Alay, NULL));
31 if (itg->xtilde) PetscCall(VecGetLayout(itg->xtilde[0], &vlay));
32 cong = PETSC_FALSE;
33 if (vlay && Alay) PetscCall(PetscLayoutCompare(Alay, vlay, &cong));
34 if (!cong) {
35 PetscCall(VecDestroyVecs(itg->maxl, &itg->btilde));
36 PetscCall(VecDestroyVecs(itg->maxl, &itg->xtilde));
37 PetscCall(VecDestroy(&itg->guess));
38 PetscCall(VecDestroy(&itg->Ax));
39 }
40 if (itg->corr) PetscCall(PetscMemzero(itg->corr, sizeof(*itg->corr) * itg->maxl * itg->maxl));
41 itg->last_b = NULL;
42 itg->last_b_state = 0;
43 if (itg->last_b_coefs) PetscCall(PetscMemzero(itg->last_b_coefs, sizeof(*itg->last_b_coefs) * itg->maxl));
44 PetscFunctionReturn(PETSC_SUCCESS);
45 }
46
KSPGuessSetUp_Fischer(KSPGuess guess)47 static PetscErrorCode KSPGuessSetUp_Fischer(KSPGuess guess)
48 {
49 KSPGuessFischer *itg = (KSPGuessFischer *)guess->data;
50
51 PetscFunctionBegin;
52 if (!itg->alpha) PetscCall(PetscMalloc1(itg->maxl, &itg->alpha));
53 if (!itg->xtilde) PetscCall(KSPCreateVecs(guess->ksp, itg->maxl, &itg->xtilde, 0, NULL));
54 if (!itg->btilde && (itg->method == 1 || itg->method == 3)) PetscCall(KSPCreateVecs(guess->ksp, itg->maxl, &itg->btilde, 0, NULL));
55 if (!itg->Ax && itg->method == 2) PetscCall(VecDuplicate(itg->xtilde[0], &itg->Ax));
56 if (!itg->guess && (itg->method == 1 || itg->method == 2)) PetscCall(VecDuplicate(itg->xtilde[0], &itg->guess));
57 if (!itg->corr && itg->method == 3) PetscCall(PetscCalloc1(itg->maxl * itg->maxl, &itg->corr));
58 if (!itg->last_b_coefs && itg->method == 3) PetscCall(PetscCalloc1(itg->maxl, &itg->last_b_coefs));
59 PetscFunctionReturn(PETSC_SUCCESS);
60 }
61
KSPGuessDestroy_Fischer(KSPGuess guess)62 static PetscErrorCode KSPGuessDestroy_Fischer(KSPGuess guess)
63 {
64 KSPGuessFischer *itg = (KSPGuessFischer *)guess->data;
65
66 PetscFunctionBegin;
67 PetscCall(PetscFree(itg->alpha));
68 PetscCall(VecDestroyVecs(itg->maxl, &itg->btilde));
69 PetscCall(VecDestroyVecs(itg->maxl, &itg->xtilde));
70 PetscCall(VecDestroy(&itg->guess));
71 PetscCall(VecDestroy(&itg->Ax));
72 PetscCall(PetscFree(itg->corr));
73 PetscCall(PetscFree(itg->last_b_coefs));
74 PetscCall(PetscFree(itg));
75 PetscCall(PetscObjectComposeFunction((PetscObject)guess, "KSPGuessFischerSetModel_C", NULL));
76 PetscFunctionReturn(PETSC_SUCCESS);
77 }
78
79 /* Note: do not change the b right-hand side as is done in the publication */
KSPGuessFormGuess_Fischer_1(KSPGuess guess,Vec b,Vec x)80 static PetscErrorCode KSPGuessFormGuess_Fischer_1(KSPGuess guess, Vec b, Vec x)
81 {
82 KSPGuessFischer *itg = (KSPGuessFischer *)guess->data;
83 PetscInt i;
84
85 PetscFunctionBegin;
86 PetscCall(VecSet(x, 0.0));
87 PetscCall(VecMDot(b, itg->curl, itg->btilde, itg->alpha));
88 if (itg->monitor) {
89 PetscCall(PetscPrintf(((PetscObject)guess)->comm, "KSPFischerGuess alphas ="));
90 for (i = 0; i < itg->curl; i++) PetscCall(PetscPrintf(((PetscObject)guess)->comm, " %g", (double)PetscAbsScalar(itg->alpha[i])));
91 PetscCall(PetscPrintf(((PetscObject)guess)->comm, "\n"));
92 }
93 PetscCall(VecMAXPY(x, itg->curl, itg->alpha, itg->xtilde));
94 PetscCall(VecCopy(x, itg->guess));
95 PetscFunctionReturn(PETSC_SUCCESS);
96 }
97
KSPGuessUpdate_Fischer_1(KSPGuess guess,Vec b,Vec x)98 static PetscErrorCode KSPGuessUpdate_Fischer_1(KSPGuess guess, Vec b, Vec x)
99 {
100 KSPGuessFischer *itg = (KSPGuessFischer *)guess->data;
101 PetscReal norm;
102 PetscInt curl = itg->curl, i;
103
104 PetscFunctionBegin;
105 if (curl == itg->maxl) {
106 PetscCall(KSP_MatMult(guess->ksp, guess->A, x, itg->btilde[0]));
107 /* PetscCall(VecCopy(b,itg->btilde[0])); */
108 PetscCall(VecNormalize(itg->btilde[0], &norm));
109 PetscCall(VecCopy(x, itg->xtilde[0]));
110 PetscCall(VecScale(itg->xtilde[0], 1.0 / norm));
111 itg->curl = 1;
112 } else {
113 if (!curl) {
114 PetscCall(VecCopy(x, itg->xtilde[curl]));
115 } else {
116 PetscCall(VecWAXPY(itg->xtilde[curl], -1.0, itg->guess, x));
117 }
118 PetscCall(KSP_MatMult(guess->ksp, guess->A, itg->xtilde[curl], itg->btilde[curl]));
119 PetscCall(VecMDot(itg->btilde[curl], curl, itg->btilde, itg->alpha));
120 for (i = 0; i < curl; i++) itg->alpha[i] = -itg->alpha[i];
121 PetscCall(VecMAXPY(itg->btilde[curl], curl, itg->alpha, itg->btilde));
122 PetscCall(VecMAXPY(itg->xtilde[curl], curl, itg->alpha, itg->xtilde));
123 PetscCall(VecNormalize(itg->btilde[curl], &norm));
124 if (norm) {
125 PetscCall(VecScale(itg->xtilde[curl], 1.0 / norm));
126 itg->curl++;
127 } else {
128 PetscCall(PetscInfo(guess->ksp, "Not increasing dimension of Fischer space because new direction is identical to previous\n"));
129 }
130 }
131 PetscFunctionReturn(PETSC_SUCCESS);
132 }
133
134 /*
135 Given a basis generated already this computes a new guess x from the new right-hand side b
136 Figures out the components of b in each btilde direction and adds them to x
137 Note: do not change the b right-hand side as is done in the publication
138 */
KSPGuessFormGuess_Fischer_2(KSPGuess guess,Vec b,Vec x)139 static PetscErrorCode KSPGuessFormGuess_Fischer_2(KSPGuess guess, Vec b, Vec x)
140 {
141 KSPGuessFischer *itg = (KSPGuessFischer *)guess->data;
142 PetscInt i;
143
144 PetscFunctionBegin;
145 PetscCall(VecSet(x, 0.0));
146 PetscCall(VecMDot(b, itg->curl, itg->xtilde, itg->alpha));
147 if (itg->monitor) {
148 PetscCall(PetscPrintf(((PetscObject)guess)->comm, "KSPFischerGuess alphas ="));
149 for (i = 0; i < itg->curl; i++) PetscCall(PetscPrintf(((PetscObject)guess)->comm, " %g", (double)PetscAbsScalar(itg->alpha[i])));
150 PetscCall(PetscPrintf(((PetscObject)guess)->comm, "\n"));
151 }
152 PetscCall(VecMAXPY(x, itg->curl, itg->alpha, itg->xtilde));
153 PetscCall(VecCopy(x, itg->guess));
154 PetscFunctionReturn(PETSC_SUCCESS);
155 }
156
KSPGuessUpdate_Fischer_2(KSPGuess guess,Vec b,Vec x)157 static PetscErrorCode KSPGuessUpdate_Fischer_2(KSPGuess guess, Vec b, Vec x)
158 {
159 KSPGuessFischer *itg = (KSPGuessFischer *)guess->data;
160 PetscScalar norm;
161 PetscInt curl = itg->curl, i;
162
163 PetscFunctionBegin;
164 if (curl == itg->maxl) {
165 PetscCall(KSP_MatMult(guess->ksp, guess->A, x, itg->Ax)); /* norm = sqrt(x'Ax) */
166 PetscCall(VecDot(x, itg->Ax, &norm));
167 PetscCall(VecCopy(x, itg->xtilde[0]));
168 PetscCall(VecScale(itg->xtilde[0], 1.0 / PetscSqrtScalar(norm)));
169 itg->curl = 1;
170 } else {
171 if (!curl) {
172 PetscCall(VecCopy(x, itg->xtilde[curl]));
173 } else {
174 PetscCall(VecWAXPY(itg->xtilde[curl], -1.0, itg->guess, x));
175 }
176 PetscCall(KSP_MatMult(guess->ksp, guess->A, itg->xtilde[curl], itg->Ax));
177 PetscCall(VecMDot(itg->Ax, curl, itg->xtilde, itg->alpha));
178 for (i = 0; i < curl; i++) itg->alpha[i] = -itg->alpha[i];
179 PetscCall(VecMAXPY(itg->xtilde[curl], curl, itg->alpha, itg->xtilde));
180
181 PetscCall(KSP_MatMult(guess->ksp, guess->A, itg->xtilde[curl], itg->Ax)); /* norm = sqrt(xtilde[curl]'Axtilde[curl]) */
182 PetscCall(VecDot(itg->xtilde[curl], itg->Ax, &norm));
183 if (PetscAbsScalar(norm) != 0.0) {
184 PetscCall(VecScale(itg->xtilde[curl], 1.0 / PetscSqrtScalar(norm)));
185 itg->curl++;
186 } else {
187 PetscCall(PetscInfo(guess->ksp, "Not increasing dimension of Fischer space because new direction is identical to previous\n"));
188 }
189 }
190 PetscFunctionReturn(PETSC_SUCCESS);
191 }
192
193 /*
194 Rather than the standard algorithm implemented in 2, we treat the provided x and b vectors to be spanning sets (not necessarily linearly independent) and use them to compute a windowed correlation matrix. Since the correlation matrix may be singular we solve it with the pseudoinverse, provided by SYEV/HEEV.
195 */
KSPGuessFormGuess_Fischer_3(KSPGuess guess,Vec b,Vec x)196 static PetscErrorCode KSPGuessFormGuess_Fischer_3(KSPGuess guess, Vec b, Vec x)
197 {
198 KSPGuessFischer *itg = (KSPGuessFischer *)guess->data;
199 PetscInt i, j, m;
200 PetscReal *s_values;
201 PetscScalar *corr, *work, *scratch_vec, zero = 0.0, one = 1.0;
202 PetscBLASInt blas_m, blas_info, blas_rank = 0, blas_lwork, blas_one = 1;
203 #if defined(PETSC_USE_COMPLEX)
204 PetscReal *rwork;
205 #endif
206
207 /* project provided b onto space of stored btildes */
208 PetscFunctionBegin;
209 PetscCall(VecSet(x, 0.0));
210 m = itg->curl;
211 itg->last_b = b;
212 PetscCall(PetscObjectStateGet((PetscObject)b, &itg->last_b_state));
213 if (m > 0) {
214 PetscCall(PetscBLASIntCast(m, &blas_m));
215 blas_lwork = (/* assume a block size of m */ blas_m + 2) * blas_m;
216 #if defined(PETSC_USE_COMPLEX)
217 PetscCall(PetscCalloc5(m * m, &corr, m, &s_values, blas_lwork, &work, 3 * m - 2, &rwork, m, &scratch_vec));
218 #else
219 PetscCall(PetscCalloc4(m * m, &corr, m, &s_values, blas_lwork, &work, m, &scratch_vec));
220 #endif
221 PetscCall(VecMDot(b, itg->curl, itg->btilde, itg->last_b_coefs));
222 for (j = 0; j < m; ++j) {
223 for (i = 0; i < m; ++i) corr[m * j + i] = itg->corr[(itg->maxl) * j + i];
224 }
225 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
226 PetscReal max_s_value = 0.0;
227 #if defined(PETSC_USE_COMPLEX)
228 PetscCallBLAS("LAPACKheev", LAPACKheev_("V", "L", &blas_m, corr, &blas_m, s_values, work, &blas_lwork, rwork, &blas_info));
229 #else
230 PetscCallBLAS("LAPACKsyev", LAPACKsyev_("V", "L", &blas_m, corr, &blas_m, s_values, work, &blas_lwork, &blas_info));
231 #endif
232
233 if (blas_info == 0) {
234 /* make corr store singular vectors and s_values store singular values */
235 for (j = 0; j < m; ++j) {
236 if (s_values[j] < 0.0) {
237 s_values[j] = PetscAbsReal(s_values[j]);
238 for (i = 0; i < m; ++i) corr[m * j + i] *= -1.0;
239 }
240 max_s_value = PetscMax(max_s_value, s_values[j]);
241 }
242
243 /* manually apply the action of the pseudoinverse */
244 PetscCallBLAS("BLASgemv", BLASgemv_("T", &blas_m, &blas_m, &one, corr, &blas_m, itg->last_b_coefs, &blas_one, &zero, scratch_vec, &blas_one));
245 for (j = 0; j < m; ++j) {
246 if (s_values[j] > itg->tol * max_s_value) {
247 scratch_vec[j] /= s_values[j];
248 blas_rank += 1;
249 } else {
250 scratch_vec[j] = 0.0;
251 }
252 }
253 PetscCallBLAS("BLASgemv", BLASgemv_("N", &blas_m, &blas_m, &one, corr, &blas_m, scratch_vec, &blas_one, &zero, itg->alpha, &blas_one));
254
255 } else {
256 PetscCall(PetscInfo(guess, "Warning eigenvalue solver failed with error code %" PetscBLASInt_FMT " - setting initial guess to zero\n", blas_info));
257 PetscCall(PetscMemzero(itg->alpha, sizeof(*itg->alpha) * itg->maxl));
258 }
259 PetscCall(PetscFPTrapPop());
260
261 if (itg->monitor && blas_info == 0) {
262 PetscCall(PetscPrintf(((PetscObject)guess)->comm, "KSPFischerGuess correlation rank = %" PetscBLASInt_FMT "\n", blas_rank));
263 PetscCall(PetscPrintf(((PetscObject)guess)->comm, "KSPFischerGuess singular values = "));
264 for (i = 0; i < itg->curl; i++) PetscCall(PetscPrintf(((PetscObject)guess)->comm, " %g", (double)s_values[i]));
265 PetscCall(PetscPrintf(((PetscObject)guess)->comm, "\n"));
266
267 PetscCall(PetscPrintf(((PetscObject)guess)->comm, "KSPFischerGuess alphas ="));
268 for (i = 0; i < itg->curl; i++) PetscCall(PetscPrintf(((PetscObject)guess)->comm, " %g", (double)PetscAbsScalar(itg->alpha[i])));
269 PetscCall(PetscPrintf(((PetscObject)guess)->comm, "\n"));
270 }
271 /* Form the initial guess by using b's projection coefficients with the xs */
272 PetscCall(VecMAXPY(x, itg->curl, itg->alpha, itg->xtilde));
273 #if defined(PETSC_USE_COMPLEX)
274 PetscCall(PetscFree5(corr, s_values, work, rwork, scratch_vec));
275 #else
276 PetscCall(PetscFree4(corr, s_values, work, scratch_vec));
277 #endif
278 }
279 PetscFunctionReturn(PETSC_SUCCESS);
280 }
281
KSPGuessUpdate_Fischer_3(KSPGuess guess,Vec b,Vec x)282 static PetscErrorCode KSPGuessUpdate_Fischer_3(KSPGuess guess, Vec b, Vec x)
283 {
284 KSPGuessFischer *itg = (KSPGuessFischer *)guess->data;
285 PetscBool rotate = itg->curl == itg->maxl ? PETSC_TRUE : PETSC_FALSE;
286 PetscInt i, j;
287 PetscObjectState b_state;
288 PetscScalar *last_column;
289 Vec oldest;
290
291 PetscFunctionBegin;
292 if (rotate) {
293 /* we have the maximum number of vectors so rotate: oldest vector is at index 0 */
294 oldest = itg->xtilde[0];
295 for (i = 1; i < itg->curl; ++i) itg->xtilde[i - 1] = itg->xtilde[i];
296 itg->xtilde[itg->curl - 1] = oldest;
297 PetscCall(VecCopy(x, itg->xtilde[itg->curl - 1]));
298
299 oldest = itg->btilde[0];
300 for (i = 1; i < itg->curl; ++i) itg->btilde[i - 1] = itg->btilde[i];
301 itg->btilde[itg->curl - 1] = oldest;
302 PetscCall(VecCopy(b, itg->btilde[itg->curl - 1]));
303 /* shift correlation matrix up and left */
304 for (j = 1; j < itg->maxl; ++j) {
305 for (i = 1; i < itg->maxl; ++i) itg->corr[(j - 1) * itg->maxl + i - 1] = itg->corr[j * itg->maxl + i];
306 }
307 } else {
308 /* append new vectors */
309 PetscCall(VecCopy(x, itg->xtilde[itg->curl]));
310 PetscCall(VecCopy(b, itg->btilde[itg->curl]));
311 itg->curl++;
312 }
313
314 /*
315 Populate new column of the correlation matrix and then copy it into the
316 row. itg->maxl is the allocated length per column: itg->curl is the actual
317 column length.
318 If possible reuse the dot products from FormGuess
319 */
320 last_column = itg->corr + (itg->curl - 1) * itg->maxl;
321 PetscCall(PetscObjectStateGet((PetscObject)b, &b_state));
322 if (b_state == itg->last_b_state && b == itg->last_b) {
323 if (rotate) {
324 for (i = 1; i < itg->maxl; ++i) itg->last_b_coefs[i - 1] = itg->last_b_coefs[i];
325 }
326 PetscCall(VecDot(b, b, &itg->last_b_coefs[itg->curl - 1]));
327 PetscCall(PetscArraycpy(last_column, itg->last_b_coefs, itg->curl));
328 } else {
329 PetscCall(VecMDot(b, itg->curl, itg->btilde, last_column));
330 }
331 for (i = 0; i < itg->curl; ++i) itg->corr[i * itg->maxl + itg->curl - 1] = last_column[i];
332 PetscFunctionReturn(PETSC_SUCCESS);
333 }
334
KSPGuessSetFromOptions_Fischer(KSPGuess guess)335 static PetscErrorCode KSPGuessSetFromOptions_Fischer(KSPGuess guess)
336 {
337 KSPGuessFischer *ITG = (KSPGuessFischer *)guess->data;
338 PetscInt nmax = 2, model[2];
339 PetscBool flg;
340
341 PetscFunctionBegin;
342 model[0] = ITG->method;
343 model[1] = ITG->maxl;
344 PetscOptionsBegin(PetscObjectComm((PetscObject)guess), ((PetscObject)guess)->prefix, "Fischer guess options", "KSPGuess");
345 PetscCall(PetscOptionsIntArray("-ksp_guess_fischer_model", "Model type and dimension of basis", "KSPGuessFischerSetModel", model, &nmax, &flg));
346 if (flg) PetscCall(KSPGuessFischerSetModel(guess, model[0], model[1]));
347 PetscCall(PetscOptionsReal("-ksp_guess_fischer_tol", "Tolerance to determine rank via ratio of singular values", "KSPGuessSetTolerance", ITG->tol, &ITG->tol, NULL));
348 PetscCall(PetscOptionsBool("-ksp_guess_fischer_monitor", "Monitor the guess", NULL, ITG->monitor, &ITG->monitor, NULL));
349 PetscOptionsEnd();
350 PetscFunctionReturn(PETSC_SUCCESS);
351 }
352
KSPGuessSetTolerance_Fischer(KSPGuess guess,PetscReal tol)353 static PetscErrorCode KSPGuessSetTolerance_Fischer(KSPGuess guess, PetscReal tol)
354 {
355 KSPGuessFischer *itg = (KSPGuessFischer *)guess->data;
356
357 PetscFunctionBegin;
358 itg->tol = tol;
359 PetscFunctionReturn(PETSC_SUCCESS);
360 }
361
KSPGuessView_Fischer(KSPGuess guess,PetscViewer viewer)362 static PetscErrorCode KSPGuessView_Fischer(KSPGuess guess, PetscViewer viewer)
363 {
364 KSPGuessFischer *itg = (KSPGuessFischer *)guess->data;
365 PetscBool isascii;
366
367 PetscFunctionBegin;
368 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
369 if (isascii) PetscCall(PetscViewerASCIIPrintf(viewer, "Model %" PetscInt_FMT ", size %" PetscInt_FMT "\n", itg->method, itg->maxl));
370 PetscFunctionReturn(PETSC_SUCCESS);
371 }
372
373 /*@
374 KSPGuessFischerSetModel - Set the Paul Fischer algorithm or its variants to compute the initial guess for a `KSPSolve()`
375
376 Logically Collective
377
378 Input Parameters:
379 + guess - the initial guess context
380 . model - use model 1, model 2, model 3, or any other number to turn it off
381 - size - size of subspace used to generate initial guess
382
383 Options Database Key:
384 . -ksp_guess_fischer_model <model,size> - uses the Fischer initial guess generator for repeated linear solves
385
386 Level: advanced
387
388 .seealso: [](ch_ksp), `KSPGuess`, `KSPGuessCreate()`, `KSPSetUseFischerGuess()`, `KSPSetGuess()`, `KSPGetGuess()`, `KSP`
389 @*/
KSPGuessFischerSetModel(KSPGuess guess,PetscInt model,PetscInt size)390 PetscErrorCode KSPGuessFischerSetModel(KSPGuess guess, PetscInt model, PetscInt size)
391 {
392 PetscFunctionBegin;
393 PetscValidHeaderSpecific(guess, KSPGUESS_CLASSID, 1);
394 PetscValidLogicalCollectiveInt(guess, model, 2);
395 PetscTryMethod(guess, "KSPGuessFischerSetModel_C", (KSPGuess, PetscInt, PetscInt), (guess, model, size));
396 PetscFunctionReturn(PETSC_SUCCESS);
397 }
398
KSPGuessFischerSetModel_Fischer(KSPGuess guess,PetscInt model,PetscInt size)399 static PetscErrorCode KSPGuessFischerSetModel_Fischer(KSPGuess guess, PetscInt model, PetscInt size)
400 {
401 KSPGuessFischer *itg = (KSPGuessFischer *)guess->data;
402
403 PetscFunctionBegin;
404 if (model == 1) {
405 guess->ops->update = KSPGuessUpdate_Fischer_1;
406 guess->ops->formguess = KSPGuessFormGuess_Fischer_1;
407 } else if (model == 2) {
408 guess->ops->update = KSPGuessUpdate_Fischer_2;
409 guess->ops->formguess = KSPGuessFormGuess_Fischer_2;
410 } else if (model == 3) {
411 guess->ops->update = KSPGuessUpdate_Fischer_3;
412 guess->ops->formguess = KSPGuessFormGuess_Fischer_3;
413 } else {
414 guess->ops->update = NULL;
415 guess->ops->formguess = NULL;
416 itg->method = 0;
417 PetscFunctionReturn(PETSC_SUCCESS);
418 }
419 if (size != itg->maxl) {
420 PetscCall(PetscFree(itg->alpha));
421 PetscCall(VecDestroyVecs(itg->maxl, &itg->btilde));
422 PetscCall(VecDestroyVecs(itg->maxl, &itg->xtilde));
423 PetscCall(VecDestroy(&itg->guess));
424 PetscCall(VecDestroy(&itg->Ax));
425 }
426 itg->method = model;
427 itg->maxl = size;
428 PetscFunctionReturn(PETSC_SUCCESS);
429 }
430
431 /*MC
432 KSPGUESSFISCHER - Implements Paul Fischer's initial guess algorithms {cite}`fischer1998projection`
433 and a non-orthogonalizing variant for situations where a linear system is solved repeatedly
434
435 Level: intermediate
436
437 Notes:
438 The algorithm is different from Fischer's paper because we do not CHANGE the right-hand side of the new
439 problem and solve the problem with an initial guess of zero, rather we solve the original problem
440 with a nonzero initial guess (this is done so that the linear solver convergence tests are based on
441 the original RHS). We use the $xtilde = x - xguess$ as the new direction so that it is not
442 mostly orthogonal to the previous solutions.
443
444 These are not intended to be used directly, they are called by `KSPSolve()` automatically with the command
445 line options `-ksp_guess_type fischer` `-ksp_guess_fischer_model <int,int>` or programmatically with
446 .vb
447 KSPGetGuess(ksp,&guess);
448 KSPGuessSetType(guess,KSPGUESSFISCHER);
449 KSPGuessFischerSetModel(guess,model,basis);
450 KSPGuessSetTolerance(guess,PETSC_MACHINE_EPSILON);
451 .ve
452 The default tolerance (which is only used in Method 3) is 32*`PETSC_MACHINE_EPSILON`. This value was chosen
453 empirically by trying a range of tolerances and picking the one that lowered the solver iteration count the most
454 with five vectors.
455
456 Method 2 is only for positive definite matrices, since it uses the energy norm.
457
458 Method 3 is not in the original paper. It is the same as the first two methods except that it
459 does not orthogonalize the input vectors or use A at all. This choice is faster but provides a
460 less effective initial guess for large (about 10) numbers of stored vectors.
461
462 Developer Note:
463 The option `-ksp_fischer_guess <int,int>` is still available for backward compatibility
464
465 .seealso: [](ch_ksp), `KSPGuess`, `KSPGuessType`, `KSP`
466 M*/
467
KSPGuessCreate_Fischer(KSPGuess guess)468 PetscErrorCode KSPGuessCreate_Fischer(KSPGuess guess)
469 {
470 KSPGuessFischer *fischer;
471
472 PetscFunctionBegin;
473 PetscCall(PetscNew(&fischer));
474 fischer->method = 1; /* defaults to method 1 */
475 fischer->maxl = 10;
476 fischer->tol = 32.0 * PETSC_MACHINE_EPSILON;
477 guess->data = fischer;
478
479 guess->ops->setfromoptions = KSPGuessSetFromOptions_Fischer;
480 guess->ops->destroy = KSPGuessDestroy_Fischer;
481 guess->ops->settolerance = KSPGuessSetTolerance_Fischer;
482 guess->ops->setup = KSPGuessSetUp_Fischer;
483 guess->ops->view = KSPGuessView_Fischer;
484 guess->ops->reset = KSPGuessReset_Fischer;
485 guess->ops->update = KSPGuessUpdate_Fischer_1;
486 guess->ops->formguess = KSPGuessFormGuess_Fischer_1;
487
488 PetscCall(PetscObjectComposeFunction((PetscObject)guess, "KSPGuessFischerSetModel_C", KSPGuessFischerSetModel_Fischer));
489 PetscFunctionReturn(PETSC_SUCCESS);
490 }
491