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) {
 30:     MatGetLayouts(guess->A,&Alay,NULL);
 31:   }
 32:   if (itg->xtilde) {
 33:     VecGetLayout(itg->xtilde[0],&vlay);
 34:   }
 35:   cong = PETSC_FALSE;
 36:   if (vlay && Alay) {
 37:     PetscLayoutCompare(Alay,vlay,&cong);
 38:   }
 39:   if (!cong) {
 40:     VecDestroyVecs(itg->maxl,&itg->btilde);
 41:     VecDestroyVecs(itg->maxl,&itg->xtilde);
 42:     VecDestroy(&itg->guess);
 43:     VecDestroy(&itg->Ax);
 44:   }
 45:   if (itg->corr) {
 46:     PetscMemzero(itg->corr,sizeof(*itg->corr)*itg->maxl*itg->maxl);
 47:   }
 48:   itg->last_b = NULL;
 49:   itg->last_b_state = 0;
 50:   if (itg->last_b_coefs) {
 51:     PetscMemzero(itg->last_b_coefs,sizeof(*itg->last_b_coefs)*itg->maxl);
 52:   }
 53:   return 0;
 54: }

 56: static PetscErrorCode KSPGuessSetUp_Fischer(KSPGuess guess)
 57: {
 58:   KSPGuessFischer *itg = (KSPGuessFischer*)guess->data;

 60:   if (!itg->alpha) {
 61:     PetscMalloc1(itg->maxl,&itg->alpha);
 62:     PetscLogObjectMemory((PetscObject)guess,itg->maxl*sizeof(PetscScalar));
 63:   }
 64:   if (!itg->xtilde) {
 65:     KSPCreateVecs(guess->ksp,itg->maxl,&itg->xtilde,0,NULL);
 66:     PetscLogObjectParents(guess,itg->maxl,itg->xtilde);
 67:   }
 68:   if (!itg->btilde && (itg->method == 1 || itg->method == 3)) {
 69:     KSPCreateVecs(guess->ksp,itg->maxl,&itg->btilde,0,NULL);
 70:     PetscLogObjectParents(guess,itg->maxl,itg->btilde);
 71:   }
 72:   if (!itg->Ax && itg->method == 2) {
 73:     VecDuplicate(itg->xtilde[0],&itg->Ax);
 74:     PetscLogObjectParent((PetscObject)guess,(PetscObject)itg->Ax);
 75:   }
 76:   if (!itg->guess && (itg->method == 1 || itg->method == 2)) {
 77:     VecDuplicate(itg->xtilde[0],&itg->guess);
 78:     PetscLogObjectParent((PetscObject)guess,(PetscObject)itg->guess);
 79:   }
 80:   if (!itg->corr && itg->method == 3) {
 81:     PetscCalloc1(itg->maxl*itg->maxl,&itg->corr);
 82:     PetscLogObjectMemory((PetscObject)guess,itg->maxl*itg->maxl*sizeof(PetscScalar));
 83:   }
 84:   if (!itg->last_b_coefs && itg->method == 3) {
 85:     PetscCalloc1(itg->maxl,&itg->last_b_coefs);
 86:     PetscLogObjectMemory((PetscObject)guess,itg->maxl*sizeof(PetscScalar));
 87:   }
 88:   return 0;
 89: }

 91: static PetscErrorCode KSPGuessDestroy_Fischer(KSPGuess guess)
 92: {
 93:   KSPGuessFischer *itg = (KSPGuessFischer*)guess->data;

 95:   PetscFree(itg->alpha);
 96:   VecDestroyVecs(itg->maxl,&itg->btilde);
 97:   VecDestroyVecs(itg->maxl,&itg->xtilde);
 98:   VecDestroy(&itg->guess);
 99:   VecDestroy(&itg->Ax);
100:   PetscFree(itg->corr);
101:   PetscFree(itg->last_b_coefs);
102:   PetscFree(itg);
103:   PetscObjectComposeFunction((PetscObject)guess,"KSPGuessFischerSetModel_C",NULL);
104:   return 0;
105: }

107: /* Note: do not change the b right hand side as is done in the publication */
108: static PetscErrorCode KSPGuessFormGuess_Fischer_1(KSPGuess guess,Vec b,Vec x)
109: {
110:   KSPGuessFischer *itg = (KSPGuessFischer*)guess->data;
111:   PetscInt        i;

113:   VecSet(x,0.0);
114:   VecMDot(b,itg->curl,itg->btilde,itg->alpha);
115:   if (itg->monitor) {
116:     PetscPrintf(((PetscObject)guess)->comm,"KSPFischerGuess alphas =");
117:     for (i=0; i<itg->curl; i++) {
118:       PetscPrintf(((PetscObject)guess)->comm," %g",(double)PetscAbsScalar(itg->alpha[i]));
119:     }
120:     PetscPrintf(((PetscObject)guess)->comm,"\n");
121:   }
122:   VecMAXPY(x,itg->curl,itg->alpha,itg->xtilde);
123:   VecCopy(x,itg->guess);
124:   return 0;
125: }

127: static PetscErrorCode KSPGuessUpdate_Fischer_1(KSPGuess guess, Vec b, Vec x)
128: {
129:   KSPGuessFischer *itg = (KSPGuessFischer*)guess->data;
130:   PetscReal       norm;
131:   int             curl = itg->curl,i;

133:   if (curl == itg->maxl) {
134:     KSP_MatMult(guess->ksp,guess->A,x,itg->btilde[0]);
135:     /* VecCopy(b,itg->btilde[0]); */
136:     VecNormalize(itg->btilde[0],&norm);
137:     VecCopy(x,itg->xtilde[0]);
138:     VecScale(itg->xtilde[0],1.0/norm);
139:     itg->curl = 1;
140:   } else {
141:     if (!curl) {
142:       VecCopy(x,itg->xtilde[curl]);
143:     } else {
144:       VecWAXPY(itg->xtilde[curl],-1.0,itg->guess,x);
145:     }
146:     KSP_MatMult(guess->ksp,guess->A,itg->xtilde[curl],itg->btilde[curl]);
147:     VecMDot(itg->btilde[curl],curl,itg->btilde,itg->alpha);
148:     for (i=0; i<curl; i++) itg->alpha[i] = -itg->alpha[i];
149:     VecMAXPY(itg->btilde[curl],curl,itg->alpha,itg->btilde);
150:     VecMAXPY(itg->xtilde[curl],curl,itg->alpha,itg->xtilde);
151:     VecNormalize(itg->btilde[curl],&norm);
152:     if (norm) {
153:       VecScale(itg->xtilde[curl],1.0/norm);
154:       itg->curl++;
155:     } else {
156:       PetscInfo(guess->ksp,"Not increasing dimension of Fischer space because new direction is identical to previous\n");
157:     }
158:   }
159:   return 0;
160: }

162: /*
163:   Given a basis generated already this computes a new guess x from the new right hand side b
164:   Figures out the components of b in each btilde direction and adds them to x
165:   Note: do not change the b right hand side as is done in the publication
166: */
167: static PetscErrorCode KSPGuessFormGuess_Fischer_2(KSPGuess guess, Vec b, Vec x)
168: {
169:   KSPGuessFischer *itg = (KSPGuessFischer*)guess->data;
170:   PetscInt        i;

172:   VecSet(x,0.0);
173:   VecMDot(b,itg->curl,itg->xtilde,itg->alpha);
174:   if (itg->monitor) {
175:     PetscPrintf(((PetscObject)guess)->comm,"KSPFischerGuess alphas =");
176:     for (i=0; i<itg->curl; i++) {
177:       PetscPrintf(((PetscObject)guess)->comm," %g",(double)PetscAbsScalar(itg->alpha[i]));
178:     }
179:     PetscPrintf(((PetscObject)guess)->comm,"\n");
180:   }
181:   VecMAXPY(x,itg->curl,itg->alpha,itg->xtilde);
182:   VecCopy(x,itg->guess);
183:   return 0;
184: }

186: static PetscErrorCode KSPGuessUpdate_Fischer_2(KSPGuess guess, Vec b, Vec x)
187: {
188:   KSPGuessFischer *itg = (KSPGuessFischer*)guess->data;
189:   PetscScalar     norm;
190:   int             curl = itg->curl,i;

192:   if (curl == itg->maxl) {
193:     KSP_MatMult(guess->ksp,guess->A,x,itg->Ax); /* norm = sqrt(x'Ax) */
194:     VecDot(x,itg->Ax,&norm);
195:     VecCopy(x,itg->xtilde[0]);
196:     VecScale(itg->xtilde[0],1.0/PetscSqrtScalar(norm));
197:     itg->curl = 1;
198:   } else {
199:     if (!curl) {
200:       VecCopy(x,itg->xtilde[curl]);
201:     } else {
202:       VecWAXPY(itg->xtilde[curl],-1.0,itg->guess,x);
203:     }
204:     KSP_MatMult(guess->ksp,guess->A,itg->xtilde[curl],itg->Ax);
205:     VecMDot(itg->Ax,curl,itg->xtilde,itg->alpha);
206:     for (i=0; i<curl; i++) itg->alpha[i] = -itg->alpha[i];
207:     VecMAXPY(itg->xtilde[curl],curl,itg->alpha,itg->xtilde);

209:     KSP_MatMult(guess->ksp,guess->A,itg->xtilde[curl],itg->Ax); /* norm = sqrt(xtilde[curl]'Axtilde[curl]) */
210:     VecDot(itg->xtilde[curl],itg->Ax,&norm);
211:     if (PetscAbsScalar(norm) != 0.0) {
212:       VecScale(itg->xtilde[curl],1.0/PetscSqrtScalar(norm));
213:       itg->curl++;
214:     } else {
215:       PetscInfo(guess->ksp,"Not increasing dimension of Fischer space because new direction is identical to previous\n");
216:     }
217:   }
218:   return 0;
219: }

221: /*
222:   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.
223: */
224: static PetscErrorCode KSPGuessFormGuess_Fischer_3(KSPGuess guess, Vec b, Vec x)
225: {
226:   KSPGuessFischer *itg = (KSPGuessFischer*)guess->data;
227:   PetscInt        i,j,m;
228:   PetscReal       *s_values;
229:   PetscScalar     *corr,*work,*scratch_vec,zero=0.0,one=1.0;
230:   PetscBLASInt    blas_m,blas_info,blas_rank=0,blas_lwork,blas_one = 1;
231: #if defined(PETSC_USE_COMPLEX)
232:   PetscReal       *rwork;
233: #endif

235:   /* project provided b onto space of stored btildes */
236:   VecSet(x,0.0);
237:   m = itg->curl;
238:   itg->last_b = b;
239:   PetscObjectStateGet((PetscObject)b,&itg->last_b_state);
240:   if (m > 0) {
241:     PetscBLASIntCast(m,&blas_m);
242:     blas_lwork = (/* assume a block size of m */blas_m+2)*blas_m;
243: #if defined(PETSC_USE_COMPLEX)
244:     PetscCalloc5(m*m,&corr,m,&s_values,blas_lwork,&work,3*m-2,&rwork,m,&scratch_vec);
245: #else
246:     PetscCalloc4(m*m,&corr,m,&s_values,blas_lwork,&work,m,&scratch_vec);
247: #endif
248:     VecMDot(b,itg->curl,itg->btilde,itg->last_b_coefs);
249:     for (j=0;j<m;++j) {
250:       for (i=0;i<m;++i) {
251:         corr[m*j+i] = itg->corr[(itg->maxl)*j+i];
252:       }
253:     }
254:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
255:     PetscReal max_s_value = 0.0;
256: #if defined(PETSC_USE_COMPLEX)
257:     PetscStackCallBLAS("LAPACKheev", LAPACKheev_("V", "L", &blas_m, corr, &blas_m, s_values, work, &blas_lwork, rwork, &blas_info));
258: #else
259:     PetscStackCallBLAS(
260:       "LAPACKsyev", LAPACKsyev_("V", "L", &blas_m, corr, &blas_m, s_values, work, &blas_lwork, &blas_info));
261: #endif

263:     if (blas_info == 0) {
264:       /* make corr store singular vectors and s_values store singular values */
265:       for (j=0; j<m; ++j) {
266:         if (s_values[j] < 0.0) {
267:           s_values[j] = PetscAbsReal(s_values[j]);
268:           for (i=0; i<m; ++i) {
269:             corr[m*j + i] *= -1.0;
270:           }
271:         }
272:         max_s_value = PetscMax(max_s_value, s_values[j]);
273:       }

275:       /* manually apply the action of the pseudoinverse */
276:       PetscStackCallBLAS("BLASgemv", BLASgemv_("T", &blas_m, &blas_m, &one, corr, &blas_m, itg->last_b_coefs, &blas_one, &zero, scratch_vec, &blas_one));
277:       for (j=0; j<m; ++j) {
278:         if (s_values[j] > itg->tol*max_s_value) {
279:           scratch_vec[j] /= s_values[j];
280:           blas_rank += 1;
281:         } else {
282:           scratch_vec[j] = 0.0;
283:         }
284:       }
285:       PetscStackCallBLAS("BLASgemv", BLASgemv_("N", &blas_m, &blas_m, &one, corr, &blas_m, scratch_vec, &blas_one, &zero, itg->alpha, &blas_one));

287:     } else {
288:       PetscInfo(guess, "Warning eigenvalue solver failed with error code %d - setting initial guess to zero\n", (int)blas_info);
289:       PetscMemzero(itg->alpha,sizeof(*itg->alpha)*itg->maxl);
290:     }
291:     PetscFPTrapPop();

293:     if (itg->monitor && blas_info == 0) {
294:       PetscPrintf(((PetscObject)guess)->comm,"KSPFischerGuess correlation rank = %d\n",(int)blas_rank);
295:       PetscPrintf(((PetscObject)guess)->comm,"KSPFischerGuess singular values = ");
296:       for (i=0; i<itg->curl; i++) {
297:         PetscPrintf(((PetscObject)guess)->comm," %g",(double)s_values[i]);
298:       }
299:       PetscPrintf(((PetscObject)guess)->comm,"\n");

301:       PetscPrintf(((PetscObject)guess)->comm,"KSPFischerGuess alphas =");
302:       for (i=0; i<itg->curl; i++) {
303:         PetscPrintf(((PetscObject)guess)->comm," %g",(double)PetscAbsScalar(itg->alpha[i]));
304:       }
305:       PetscPrintf(((PetscObject)guess)->comm,"\n");
306:     }
307:     /* Form the initial guess by using b's projection coefficients with the xs */
308:     VecMAXPY(x,itg->curl,itg->alpha,itg->xtilde);
309: #if defined(PETSC_USE_COMPLEX)
310:     PetscFree5(corr, s_values, work, rwork, scratch_vec);
311: #else
312:     PetscFree4(corr, s_values, work, scratch_vec);
313: #endif
314:   }
315:   return 0;
316: }

318: static PetscErrorCode KSPGuessUpdate_Fischer_3(KSPGuess guess, Vec b, Vec x)
319: {
320:   KSPGuessFischer  *itg = (KSPGuessFischer*)guess->data;
321:   PetscBool        rotate = itg->curl == itg->maxl ? PETSC_TRUE : PETSC_FALSE;
322:   PetscInt         i,j;
323:   PetscObjectState b_state;
324:   PetscScalar      *last_column;
325:   Vec              oldest;

327:   if (rotate) {
328:     /* we have the maximum number of vectors so rotate: oldest vector is at index 0 */
329:     oldest = itg->xtilde[0];
330:     for (i=1;i<itg->curl;++i) {
331:       itg->xtilde[i-1] = itg->xtilde[i];
332:     }
333:     itg->xtilde[itg->curl-1] = oldest;
334:     VecCopy(x,itg->xtilde[itg->curl-1]);

336:     oldest = itg->btilde[0];
337:     for (i=1;i<itg->curl;++i) {
338:       itg->btilde[i-1] = itg->btilde[i];
339:     }
340:     itg->btilde[itg->curl-1] = oldest;
341:     VecCopy(b,itg->btilde[itg->curl-1]);
342:     /* shift correlation matrix up and left */
343:     for (j=1; j<itg->maxl; ++j) {
344:       for (i=1; i<itg->maxl; ++i) {
345:         itg->corr[(j-1)*itg->maxl+i-1]=itg->corr[j*itg->maxl+i];
346:       }
347:     }
348:   } else {
349:     /* append new vectors */
350:     VecCopy(x,itg->xtilde[itg->curl]);
351:     VecCopy(b,itg->btilde[itg->curl]);
352:     itg->curl++;
353:   }

355:   /*
356:       Populate new column of the correlation matrix and then copy it into the
357:       row. itg->maxl is the allocated length per column: itg->curl is the actual
358:       column length.
359:       If possible reuse the dot products from FormGuess
360:   */
361:   last_column = itg->corr+(itg->curl-1)*itg->maxl;
362:   PetscObjectStateGet((PetscObject)b,&b_state);
363:   if (b_state == itg->last_b_state && b == itg->last_b) {
364:     if (rotate) {
365:       for (i=1; i<itg->maxl; ++i) {
366:         itg->last_b_coefs[i-1] = itg->last_b_coefs[i];
367:       }
368:     }
369:     VecDot(b,b,&itg->last_b_coefs[itg->curl-1]);
370:     PetscArraycpy(last_column,itg->last_b_coefs,itg->curl);
371:   } else {
372:     VecMDot(b,itg->curl,itg->btilde,last_column);
373:   }
374:   for (i=0;i<itg->curl;++i) {
375:     itg->corr[i*itg->maxl+itg->curl-1] = last_column[i];
376:   }
377:   return 0;
378: }

380: static PetscErrorCode KSPGuessSetFromOptions_Fischer(KSPGuess guess)
381: {
382:   KSPGuessFischer *ITG = (KSPGuessFischer *)guess->data;
383:   PetscInt        nmax = 2, model[2];
384:   PetscBool       flg;
385:   PetscErrorCode  ierr;

387:   model[0] = ITG->method;
388:   model[1] = ITG->maxl;
389:   PetscOptionsBegin(PetscObjectComm((PetscObject)guess),((PetscObject)guess)->prefix,"Fischer guess options","KSPGuess");
390:   PetscOptionsIntArray("-ksp_guess_fischer_model","Model type and dimension of basis","KSPGuessFischerSetModel",model,&nmax,&flg);
391:   if (flg) {
392:     KSPGuessFischerSetModel(guess,model[0],model[1]);
393:   }
394:   PetscOptionsReal("-ksp_guess_fischer_tol","Tolerance to determine rank via ratio of singular values","KSPGuessSetTolerance",ITG->tol,&ITG->tol,NULL);
395:   PetscOptionsBool("-ksp_guess_fischer_monitor","Monitor the guess",NULL,ITG->monitor,&ITG->monitor,NULL);
396:   PetscOptionsEnd();
397:   return 0;
398: }

400: static PetscErrorCode KSPGuessSetTolerance_Fischer(KSPGuess guess,PetscReal tol)
401: {
402:   KSPGuessFischer *itg = (KSPGuessFischer *)guess->data;

404:   itg->tol = tol;
405:   return 0;
406: }

408: static PetscErrorCode KSPGuessView_Fischer(KSPGuess guess,PetscViewer viewer)
409: {
410:   KSPGuessFischer *itg = (KSPGuessFischer*)guess->data;
411:   PetscBool       isascii;

413:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
414:   if (isascii) {
415:     PetscViewerASCIIPrintf(viewer,"Model %D, size %D\n",itg->method,itg->maxl);
416:   }
417:   return 0;
418: }

420: /*@
421:    KSPGuessFischerSetModel - Use the Paul Fischer algorithm or its variants

423:    Logically Collective on guess

425:    Input Parameters:
426: +  guess - the initial guess context
427: .  model - use model 1, model 2, model 3, or any other number to turn it off
428: -  size  - size of subspace used to generate initial guess

430:     Options Database:
431: .   -ksp_guess_fischer_model <model,size> - uses the Fischer initial guess generator for repeated linear solves

433:    Level: advanced

435: .seealso: KSPGuess, KSPGuessCreate(), KSPSetUseFischerGuess(), KSPSetGuess(), KSPGetGuess(), KSP
436: @*/
437: PetscErrorCode  KSPGuessFischerSetModel(KSPGuess guess,PetscInt model,PetscInt size)
438: {
441:   PetscTryMethod(guess,"KSPGuessFischerSetModel_C",(KSPGuess,PetscInt,PetscInt),(guess,model,size));
442:   return 0;
443: }

445: static PetscErrorCode KSPGuessFischerSetModel_Fischer(KSPGuess guess,PetscInt model,PetscInt size)
446: {
447:   KSPGuessFischer *itg = (KSPGuessFischer*)guess->data;

449:   if (model == 1) {
450:     guess->ops->update    = KSPGuessUpdate_Fischer_1;
451:     guess->ops->formguess = KSPGuessFormGuess_Fischer_1;
452:   } else if (model == 2) {
453:     guess->ops->update    = KSPGuessUpdate_Fischer_2;
454:     guess->ops->formguess = KSPGuessFormGuess_Fischer_2;
455:   } else if (model == 3) {
456:     guess->ops->update    = KSPGuessUpdate_Fischer_3;
457:     guess->ops->formguess = KSPGuessFormGuess_Fischer_3;
458:   } else {
459:     guess->ops->update    = NULL;
460:     guess->ops->formguess = NULL;
461:     itg->method           = 0;
462:     return 0;
463:   }
464:   if (size != itg->maxl) {
465:     PetscFree(itg->alpha);
466:     VecDestroyVecs(itg->maxl,&itg->btilde);
467:     VecDestroyVecs(itg->maxl,&itg->xtilde);
468:     VecDestroy(&itg->guess);
469:     VecDestroy(&itg->Ax);
470:   }
471:   itg->method = model;
472:   itg->maxl   = size;
473:   return 0;
474: }

476: /*
477:     KSPGUESSFISCHER - Implements Paul Fischer's two initial guess algorithms and a nonorthogonalizing variant for situations where
478:     a linear system is solved repeatedly

480:   References:
481: . * - https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19940020363_1994020363.pdf

483:    Notes:
484:     the algorithm is different from Fischer's paper because we do not CHANGE the right hand side of the new
485:     problem and solve the problem with an initial guess of zero, rather we solve the original problem
486:     with a nonzero initial guess (this is done so that the linear solver convergence tests are based on
487:     the original RHS). We use the xtilde = x - xguess as the new direction so that it is not
488:     mostly orthogonal to the previous solutions.

490:     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
491: .vb
492:     KSPGetGuess(ksp,&guess);
493:     KSPGuessSetType(guess,KSPGUESSFISCHER);
494:     KSPGuessFischerSetModel(guess,model,basis);
495:     KSPGuessSetTolerance(guess,PETSC_MACHINE_EPSILON);

497:     The default tolerance (which is only used in Method 3) is 32*PETSC_MACHINE_EPSILON. This value was chosen
498:     empirically by trying a range of tolerances and picking the one that lowered the solver iteration count the most
499:     with five vectors.

501:     Method 2 is only for positive definite matrices, since it uses the A norm.

503:     Method 3 is not in the original paper. It is the same as the first two methods except that it
504:     does not orthogonalize the input vectors or use A at all. This choice is faster but provides a
505:     less effective initial guess for large (about 10) numbers of stored vectors.

507:     Developer note:
508:       The option -ksp_fischer_guess <int,int> is still available for backward compatibility

510:     Level: intermediate

512: @*/
513: PetscErrorCode KSPGuessCreate_Fischer(KSPGuess guess)
514: {
515:   KSPGuessFischer *fischer;

517:   PetscNewLog(guess,&fischer);
518:   fischer->method = 1;  /* defaults to method 1 */
519:   fischer->maxl   = 10;
520:   fischer->tol    = 32.0*PETSC_MACHINE_EPSILON;
521:   guess->data     = fischer;

523:   guess->ops->setfromoptions = KSPGuessSetFromOptions_Fischer;
524:   guess->ops->destroy        = KSPGuessDestroy_Fischer;
525:   guess->ops->settolerance   = KSPGuessSetTolerance_Fischer;
526:   guess->ops->setup          = KSPGuessSetUp_Fischer;
527:   guess->ops->view           = KSPGuessView_Fischer;
528:   guess->ops->reset          = KSPGuessReset_Fischer;
529:   guess->ops->update         = KSPGuessUpdate_Fischer_1;
530:   guess->ops->formguess      = KSPGuessFormGuess_Fischer_1;

532:   PetscObjectComposeFunction((PetscObject)guess,"KSPGuessFischerSetModel_C",KSPGuessFischerSetModel_Fischer);
533:   return 0;
534: }