Actual source code: fischer.c

petsc-3.14.6 2021-03-30
Report Typos and Errors
  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: }