Actual source code: fischer.c
1: #include <petsc/private/kspimpl.h>
2: #include <petscblaslapack.h>
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;
21: static PetscErrorCode KSPGuessReset_Fischer(KSPGuess guess)
22: {
23: KSPGuessFischer *itg = (KSPGuessFischer *)guess->data;
24: PetscLayout Alay = NULL, vlay = NULL;
25: PetscBool cong;
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: }
47: static PetscErrorCode KSPGuessSetUp_Fischer(KSPGuess guess)
48: {
49: KSPGuessFischer *itg = (KSPGuessFischer *)guess->data;
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: }
62: static PetscErrorCode KSPGuessDestroy_Fischer(KSPGuess guess)
63: {
64: KSPGuessFischer *itg = (KSPGuessFischer *)guess->data;
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: }
79: /* Note: do not change the b right hand side as is done in the publication */
80: static PetscErrorCode KSPGuessFormGuess_Fischer_1(KSPGuess guess, Vec b, Vec x)
81: {
82: KSPGuessFischer *itg = (KSPGuessFischer *)guess->data;
83: PetscInt i;
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: }
98: static PetscErrorCode KSPGuessUpdate_Fischer_1(KSPGuess guess, Vec b, Vec x)
99: {
100: KSPGuessFischer *itg = (KSPGuessFischer *)guess->data;
101: PetscReal norm;
102: int curl = itg->curl, i;
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: }
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: */
139: static PetscErrorCode KSPGuessFormGuess_Fischer_2(KSPGuess guess, Vec b, Vec x)
140: {
141: KSPGuessFischer *itg = (KSPGuessFischer *)guess->data;
142: PetscInt i;
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: }
157: static PetscErrorCode KSPGuessUpdate_Fischer_2(KSPGuess guess, Vec b, Vec x)
158: {
159: KSPGuessFischer *itg = (KSPGuessFischer *)guess->data;
160: PetscScalar norm;
161: int curl = itg->curl, i;
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));
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: }
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: */
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
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
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: }
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));
255: } else {
256: PetscCall(PetscInfo(guess, "Warning eigenvalue solver failed with error code %d - setting initial guess to zero\n", (int)blas_info));
257: PetscCall(PetscMemzero(itg->alpha, sizeof(*itg->alpha) * itg->maxl));
258: }
259: PetscCall(PetscFPTrapPop());
261: if (itg->monitor && blas_info == 0) {
262: PetscCall(PetscPrintf(((PetscObject)guess)->comm, "KSPFischerGuess correlation rank = %d\n", (int)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"));
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: }
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;
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]));
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: }
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: }
335: static PetscErrorCode KSPGuessSetFromOptions_Fischer(KSPGuess guess)
336: {
337: KSPGuessFischer *ITG = (KSPGuessFischer *)guess->data;
338: PetscInt nmax = 2, model[2];
339: PetscBool flg;
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: }
353: static PetscErrorCode KSPGuessSetTolerance_Fischer(KSPGuess guess, PetscReal tol)
354: {
355: KSPGuessFischer *itg = (KSPGuessFischer *)guess->data;
357: PetscFunctionBegin;
358: itg->tol = tol;
359: PetscFunctionReturn(PETSC_SUCCESS);
360: }
362: static PetscErrorCode KSPGuessView_Fischer(KSPGuess guess, PetscViewer viewer)
363: {
364: KSPGuessFischer *itg = (KSPGuessFischer *)guess->data;
365: PetscBool isascii;
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: }
373: /*@
374: KSPGuessFischerSetModel - Set the Paul Fischer algorithm or its variants to compute the initial guess for a `KSPSolve()`
376: Logically Collective
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
383: Options Database Key:
384: . -ksp_guess_fischer_model <model,size> - uses the Fischer initial guess generator for repeated linear solves
386: Level: advanced
388: .seealso: [](ch_ksp), `KSPGuess`, `KSPGuessCreate()`, `KSPSetUseFischerGuess()`, `KSPSetGuess()`, `KSPGetGuess()`, `KSP`
389: @*/
390: PetscErrorCode KSPGuessFischerSetModel(KSPGuess guess, PetscInt model, PetscInt size)
391: {
392: PetscFunctionBegin;
395: PetscTryMethod(guess, "KSPGuessFischerSetModel_C", (KSPGuess, PetscInt, PetscInt), (guess, model, size));
396: PetscFunctionReturn(PETSC_SUCCESS);
397: }
399: static PetscErrorCode KSPGuessFischerSetModel_Fischer(KSPGuess guess, PetscInt model, PetscInt size)
400: {
401: KSPGuessFischer *itg = (KSPGuessFischer *)guess->data;
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: }
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
435: Level: intermediate
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.
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.
456: Method 2 is only for positive definite matrices, since it uses the energy norm.
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.
462: Developer Note:
463: The option `-ksp_fischer_guess <int,int>` is still available for backward compatibility
465: .seealso: [](ch_ksp), `KSPGuess`, `KSPGuessType`, `KSP`
466: M*/
468: PetscErrorCode KSPGuessCreate_Fischer(KSPGuess guess)
469: {
470: KSPGuessFischer *fischer;
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;
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;
488: PetscCall(PetscObjectComposeFunction((PetscObject)guess, "KSPGuessFischerSetModel_C", KSPGuessFischerSetModel_Fischer));
489: PetscFunctionReturn(PETSC_SUCCESS);
490: }