Actual source code: fischer.c
petsc-3.14.6 2021-03-30
1: #include <petsc/private/kspimpl.h>
3: typedef struct {
4: PetscInt method; /* 1 or 2 */
5: PetscInt curl; /* Current number of basis vectors */
6: PetscInt maxl; /* Maximum number of basis vectors */
7: PetscBool monitor;
8: PetscScalar *alpha; /* */
9: Vec *xtilde; /* Saved x vectors */
10: Vec *btilde; /* Saved b vectors, method 1 */
11: Vec Ax; /* method 2 */
12: Vec guess;
13: } KSPGuessFischer;
15: static PetscErrorCode KSPGuessReset_Fischer(KSPGuess guess)
16: {
17: KSPGuessFischer *itg = (KSPGuessFischer*)guess->data;
18: PetscLayout Alay = NULL,vlay = NULL;
19: PetscBool cong;
20: PetscErrorCode ierr;
23: itg->curl = 0;
24: /* destroy vectors if the size of the linear system has changed */
25: if (guess->A) {
26: MatGetLayouts(guess->A,&Alay,NULL);
27: }
28: if (itg->xtilde) {
29: VecGetLayout(itg->xtilde[0],&vlay);
30: }
31: cong = PETSC_FALSE;
32: if (vlay && Alay) {
33: PetscLayoutCompare(Alay,vlay,&cong);
34: }
35: if (!cong) {
36: VecDestroyVecs(itg->maxl,&itg->btilde);
37: VecDestroyVecs(itg->maxl,&itg->xtilde);
38: VecDestroy(&itg->guess);
39: VecDestroy(&itg->Ax);
40: }
41: return(0);
42: }
44: static PetscErrorCode KSPGuessSetUp_Fischer(KSPGuess guess)
45: {
46: KSPGuessFischer *itg = (KSPGuessFischer*)guess->data;
47: PetscErrorCode ierr;
50: if (!itg->alpha) {
51: PetscMalloc1(itg->maxl,&itg->alpha);
52: PetscLogObjectMemory((PetscObject)guess,itg->maxl*sizeof(PetscScalar));
53: }
54: if (!itg->xtilde) {
55: KSPCreateVecs(guess->ksp,itg->maxl,&itg->xtilde,0,NULL);
56: PetscLogObjectParents(guess,itg->maxl,itg->xtilde);
57: }
58: if (!itg->btilde && itg->method == 1) {
59: KSPCreateVecs(guess->ksp,itg->maxl,&itg->btilde,0,NULL);
60: PetscLogObjectParents(guess,itg->maxl,itg->btilde);
61: }
62: if (!itg->Ax && itg->method == 2) {
63: VecDuplicate(itg->xtilde[0],&itg->Ax);
64: PetscLogObjectParent((PetscObject)guess,(PetscObject)itg->Ax);
65: }
66: if (!itg->guess) {
67: VecDuplicate(itg->xtilde[0],&itg->guess);
68: PetscLogObjectParent((PetscObject)guess,(PetscObject)itg->guess);
69: }
70: return(0);
71: }
73: static PetscErrorCode KSPGuessDestroy_Fischer(KSPGuess guess)
74: {
75: KSPGuessFischer *itg = (KSPGuessFischer*)guess->data;
76: PetscErrorCode ierr;
79: PetscFree(itg->alpha);
80: VecDestroyVecs(itg->maxl,&itg->btilde);
81: VecDestroyVecs(itg->maxl,&itg->xtilde);
82: VecDestroy(&itg->guess);
83: VecDestroy(&itg->Ax);
84: PetscFree(itg);
85: PetscObjectComposeFunction((PetscObject)guess,"KSPGuessFischerSetModel_C",NULL);
86: return(0);
87: }
89: /* Note: do not change the b right hand side as is done in the publication */
90: static PetscErrorCode KSPGuessFormGuess_Fischer_1(KSPGuess guess,Vec b,Vec x)
91: {
92: KSPGuessFischer *itg = (KSPGuessFischer*)guess->data;
93: PetscErrorCode ierr;
94: PetscInt i;
97: VecSet(x,0.0);
98: VecMDot(b,itg->curl,itg->btilde,itg->alpha);
99: if (itg->monitor) {
100: PetscPrintf(((PetscObject)guess)->comm,"KSPFischerGuess alphas = ");
101: for (i=0; i<itg->curl; i++) {
102: PetscPrintf(((PetscObject)guess)->comm,"%g ",(double)PetscAbsScalar(itg->alpha[i]));
103: }
104: PetscPrintf(((PetscObject)guess)->comm,"\n");
105: }
106: VecMAXPY(x,itg->curl,itg->alpha,itg->xtilde);
107: VecCopy(x,itg->guess);
108: return(0);
109: }
111: static PetscErrorCode KSPGuessUpdate_Fischer_1(KSPGuess guess, Vec b, Vec x)
112: {
113: KSPGuessFischer *itg = (KSPGuessFischer*)guess->data;
114: PetscReal norm;
115: PetscErrorCode ierr;
116: int curl = itg->curl,i;
119: if (curl == itg->maxl) {
120: KSP_MatMult(guess->ksp,guess->A,x,itg->btilde[0]);
121: /* VecCopy(b,itg->btilde[0]); */
122: VecNormalize(itg->btilde[0],&norm);
123: VecCopy(x,itg->xtilde[0]);
124: VecScale(itg->xtilde[0],1.0/norm);
125: itg->curl = 1;
126: } else {
127: if (!curl) {
128: VecCopy(x,itg->xtilde[curl]);
129: } else {
130: VecWAXPY(itg->xtilde[curl],-1.0,itg->guess,x);
131: }
132: KSP_MatMult(guess->ksp,guess->A,itg->xtilde[curl],itg->btilde[curl]);
133: VecMDot(itg->btilde[curl],curl,itg->btilde,itg->alpha);
134: for (i=0; i<curl; i++) itg->alpha[i] = -itg->alpha[i];
135: VecMAXPY(itg->btilde[curl],curl,itg->alpha,itg->btilde);
136: VecMAXPY(itg->xtilde[curl],curl,itg->alpha,itg->xtilde);
137: VecNormalize(itg->btilde[curl],&norm);
138: if (norm) {
139: VecScale(itg->xtilde[curl],1.0/norm);
140: itg->curl++;
141: } else {
142: PetscInfo(guess->ksp,"Not increasing dimension of Fischer space because new direction is identical to previous\n");
143: }
144: }
145: return(0);
146: }
148: /*
149: Given a basis generated already this computes a new guess x from the new right hand side b
150: Figures out the components of b in each btilde direction and adds them to x
151: Note: do not change the b right hand side as is done in the publication
152: */
153: static PetscErrorCode KSPGuessFormGuess_Fischer_2(KSPGuess guess, Vec b, Vec x)
154: {
155: KSPGuessFischer *itg = (KSPGuessFischer*)guess->data;
156: PetscErrorCode ierr;
157: PetscInt i;
160: VecSet(x,0.0);
161: VecMDot(b,itg->curl,itg->xtilde,itg->alpha);
162: if (itg->monitor) {
163: PetscPrintf(((PetscObject)guess)->comm,"KSPFischerGuess alphas = ");
164: for (i=0; i<itg->curl; i++) {
165: PetscPrintf(((PetscObject)guess)->comm,"%g ",(double)PetscAbsScalar(itg->alpha[i]));
166: }
167: PetscPrintf(((PetscObject)guess)->comm,"\n");
168: }
169: VecMAXPY(x,itg->curl,itg->alpha,itg->xtilde);
170: VecCopy(x,itg->guess);
171: return(0);
172: }
174: static PetscErrorCode KSPGuessUpdate_Fischer_2(KSPGuess guess, Vec b, Vec x)
175: {
176: KSPGuessFischer *itg = (KSPGuessFischer*)guess->data;
177: PetscScalar norm;
178: PetscErrorCode ierr;
179: int curl = itg->curl,i;
182: if (curl == itg->maxl) {
183: KSP_MatMult(guess->ksp,guess->A,x,itg->Ax); /* norm = sqrt(x'Ax) */
184: VecDot(x,itg->Ax,&norm);
185: VecCopy(x,itg->xtilde[0]);
186: VecScale(itg->xtilde[0],1.0/PetscSqrtScalar(norm));
187: itg->curl = 1;
188: } else {
189: if (!curl) {
190: VecCopy(x,itg->xtilde[curl]);
191: } else {
192: VecWAXPY(itg->xtilde[curl],-1.0,itg->guess,x);
193: }
194: KSP_MatMult(guess->ksp,guess->A,itg->xtilde[curl],itg->Ax);
195: VecMDot(itg->Ax,curl,itg->xtilde,itg->alpha);
196: for (i=0; i<curl; i++) itg->alpha[i] = -itg->alpha[i];
197: VecMAXPY(itg->xtilde[curl],curl,itg->alpha,itg->xtilde);
199: KSP_MatMult(guess->ksp,guess->A,itg->xtilde[curl],itg->Ax); /* norm = sqrt(xtilde[curl]'Axtilde[curl]) */
200: VecDot(itg->xtilde[curl],itg->Ax,&norm);
201: if (PetscAbsScalar(norm) != 0.0) {
202: VecScale(itg->xtilde[curl],1.0/PetscSqrtScalar(norm));
203: itg->curl++;
204: } else {
205: PetscInfo(guess->ksp,"Not increasing dimension of Fischer space because new direction is identical to previous\n");
206: }
207: }
208: return(0);
209: }
211: static PetscErrorCode KSPGuessSetFromOptions_Fischer(KSPGuess guess)
212: {
213: KSPGuessFischer *ITG = (KSPGuessFischer *)guess->data;
214: PetscInt nmax = 2, model[2];
215: PetscBool flg;
216: PetscErrorCode ierr;
219: model[0] = ITG->method;
220: model[1] = ITG->maxl;
221: PetscOptionsBegin(PetscObjectComm((PetscObject)guess),((PetscObject)guess)->prefix,"Fischer guess options","KSPGuess");
222: PetscOptionsIntArray("-ksp_guess_fischer_model","Model type and dimension of basis","KSPGuessFischerSetModel",model,&nmax,&flg);
223: if (flg) {
224: KSPGuessFischerSetModel(guess,model[0],model[1]);
225: }
226: PetscOptionsBool("-ksp_guess_fischer_monitor","Monitor the guess",NULL,ITG->monitor,&ITG->monitor,NULL);
227: PetscOptionsEnd();
228: return(0);
229: }
231: static PetscErrorCode KSPGuessView_Fischer(KSPGuess guess,PetscViewer viewer)
232: {
233: KSPGuessFischer *itg = (KSPGuessFischer*)guess->data;
234: PetscBool isascii;
235: PetscErrorCode ierr;
238: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
239: if (isascii) {
240: PetscViewerASCIIPrintf(viewer,"Model %D, size %D\n",itg->method,itg->maxl);
241: }
242: return(0);
243: }
245: /*@
246: KSPGuessFischerSetModel - Use the Paul Fischer algorithm
248: Logically Collective on guess
250: Input Parameters:
251: + guess - the initial guess context
252: . model - use model 1, model 2 or any other number to turn it off
253: - size - size of subspace used to generate initial guess
255: Options Database:
256: . -ksp_guess_fischer_model <model,size> - uses the Fischer initial guess generator for repeated linear solves
258: Level: advanced
260: .seealso: KSPGuess, KSPGuessCreate(), KSPSetUseFischerGuess(), KSPSetGuess(), KSPGetGuess(), KSP
261: @*/
262: PetscErrorCode KSPGuessFischerSetModel(KSPGuess guess,PetscInt model,PetscInt size)
263: {
270: PetscTryMethod(guess,"KSPGuessFischerSetModel_C",(KSPGuess,PetscInt,PetscInt),(guess,model,size));
271: return(0);
272: }
274: static PetscErrorCode KSPGuessFischerSetModel_Fischer(KSPGuess guess,PetscInt model,PetscInt size)
275: {
276: KSPGuessFischer *itg = (KSPGuessFischer*)guess->data;
280: if (model == 1) {
281: guess->ops->update = KSPGuessUpdate_Fischer_1;
282: guess->ops->formguess = KSPGuessFormGuess_Fischer_1;
283: } else if (model == 2) {
284: guess->ops->update = KSPGuessUpdate_Fischer_2;
285: guess->ops->formguess = KSPGuessFormGuess_Fischer_2;
286: } else {
287: guess->ops->update = NULL;
288: guess->ops->formguess = NULL;
289: itg->method = 0;
290: return(0);
291: }
292: if (size != itg->maxl) {
293: PetscFree(itg->alpha);
294: VecDestroyVecs(itg->maxl,&itg->btilde);
295: VecDestroyVecs(itg->maxl,&itg->xtilde);
296: VecDestroy(&itg->guess);
297: VecDestroy(&itg->Ax);
298: }
299: itg->method = model;
300: itg->maxl = size;
301: return(0);
302: }
304: /*
305: KSPGUESSFISCHER - Implements Paul Fischer's initial guess algorithm Method 1 and 2 for situations where
306: a linear system is solved repeatedly
308: References:
309: . 1. - https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19940020363_1994020363.pdf
311: Notes:
312: the algorithm is different from the paper because we do not CHANGE the right hand side of the new
313: problem and solve the problem with an initial guess of zero, rather we solve the original problem
314: with a nonzero initial guess (this is done so that the linear solver convergence tests are based on
315: the original RHS). We use the xtilde = x - xguess as the new direction so that it is not
316: mostly orthogonal to the previous solutions.
318: 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
319: .vb
320: KSPGetGuess(ksp,&guess);
321: KSPGuessSetType(guess,KSPGUESSFISCHER);
322: KSPGuessFischerSetModel(guess,model,basis);
324: Method 2 is only for positive definite matrices, since it uses the A norm.
326: Developer note: the option -ksp_fischer_guess <int,int> is still available for backward compatibility
328: Level: intermediate
330: @*/
331: PetscErrorCode KSPGuessCreate_Fischer(KSPGuess guess)
332: {
333: KSPGuessFischer *fischer;
334: PetscErrorCode ierr;
337: PetscNewLog(guess,&fischer);
338: fischer->method = 1; /* defaults to method 1 */
339: fischer->maxl = 10;
340: guess->data = fischer;
342: guess->ops->setfromoptions = KSPGuessSetFromOptions_Fischer;
343: guess->ops->destroy = KSPGuessDestroy_Fischer;
344: guess->ops->setup = KSPGuessSetUp_Fischer;
345: guess->ops->view = KSPGuessView_Fischer;
346: guess->ops->reset = KSPGuessReset_Fischer;
347: guess->ops->update = KSPGuessUpdate_Fischer_1;
348: guess->ops->formguess = KSPGuessFormGuess_Fischer_1;
350: PetscObjectComposeFunction((PetscObject)guess,"KSPGuessFischerSetModel_C",KSPGuessFischerSetModel_Fischer);
351: return(0);
352: }