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