Actual source code: fischer.c
petsc-3.9.4 2018-09-11
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 KSP
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: .keywords: set, options, prefix, database
262: .seealso: KSPGuess, KSPGuessCreate(), KSPSetUseFischerGuess(), KSPSetGuess(), KSPGetGuess(), KSP
263: @*/
264: PetscErrorCode KSPGuessFischerSetModel(KSPGuess guess,PetscInt model,PetscInt size)
265: {
272: PetscTryMethod(guess,"KSPGuessFischerSetModel_C",(KSPGuess,PetscInt,PetscInt),(guess,model,size));
273: return(0);
274: }
276: static PetscErrorCode KSPGuessFischerSetModel_Fischer(KSPGuess guess,PetscInt model,PetscInt size)
277: {
278: KSPGuessFischer *itg = (KSPGuessFischer*)guess->data;
282: if (model == 1) {
283: guess->ops->update = KSPGuessUpdate_Fischer_1;
284: guess->ops->formguess = KSPGuessFormGuess_Fischer_1;
285: } else if (model == 2) {
286: guess->ops->update = KSPGuessUpdate_Fischer_2;
287: guess->ops->formguess = KSPGuessFormGuess_Fischer_2;
288: } else {
289: guess->ops->update = NULL;
290: guess->ops->formguess = NULL;
291: itg->method = 0;
292: return(0);
293: }
294: if (size != itg->maxl) {
295: PetscFree(itg->alpha);
296: VecDestroyVecs(itg->maxl,&itg->btilde);
297: VecDestroyVecs(itg->maxl,&itg->xtilde);
298: VecDestroy(&itg->guess);
299: VecDestroy(&itg->Ax);
300: }
301: itg->method = model;
302: itg->maxl = size;
303: return(0);
304: }
306: /*
307: KSPGUESSFISCHER - Implements Paul Fischer's initial guess algorithm Method 1 and 2 for situations where
308: a linear system is solved repeatedly
310: References:
311: . 1. - http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19940020363_1994020363.pdf
313: Notes: the algorithm is different from the paper because we do not CHANGE the right hand side of the new
314: problem and solve the problem with an initial guess of zero, rather we solve the original new problem
315: with a nonzero initial guess (this is done so that the linear solver convergence tests are based on
316: the original RHS.) But we use the xtilde = x - xguess as the new direction so that it is not
317: mostly orthogonal to the previous solutions.
319: 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
320: .vb
321: KSPGetGuess(ksp,&guess);
322: KSPGuessSetType(guess,KSPGUESSFISCHER);
323: KSPGuessFischerSetModel(guess,model,basis);
325: Method 2 is only for positive definite matrices, since it uses the A norm.
327: Developer note: the option -ksp_fischer_guess <int,int> is still available for backward compatibility
329: Level: intermediate
331: @*/
332: PetscErrorCode KSPGuessCreate_Fischer(KSPGuess guess)
333: {
334: KSPGuessFischer *fischer;
335: PetscErrorCode ierr;
338: PetscNewLog(guess,&fischer);
339: fischer->method = 1; /* defaults to method 1 */
340: fischer->maxl = 10;
341: guess->data = fischer;
343: guess->ops->setfromoptions = KSPGuessSetFromOptions_Fischer;
344: guess->ops->destroy = KSPGuessDestroy_Fischer;
345: guess->ops->setup = KSPGuessSetUp_Fischer;
346: guess->ops->view = KSPGuessView_Fischer;
347: guess->ops->reset = KSPGuessReset_Fischer;
348: guess->ops->update = KSPGuessUpdate_Fischer_1;
349: guess->ops->formguess = KSPGuessFormGuess_Fischer_1;
351: PetscObjectComposeFunction((PetscObject)guess,"KSPGuessFischerSetModel_C",KSPGuessFischerSetModel_Fischer);
352: return(0);
353: }