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:   PetscFunctionBegin;
     28:   itg->curl = 0;
     29:   /* destroy vectors if the size of the linear system has changed */
     30:   if (guess->A) PetscCall(MatGetLayouts(guess->A, &Alay, NULL));
     31:   if (itg->xtilde) PetscCall(VecGetLayout(itg->xtilde[0], &vlay));
     32:   cong = PETSC_FALSE;
     33:   if (vlay && Alay) PetscCall(PetscLayoutCompare(Alay, vlay, &cong));
     34:   if (!cong) {
     35:     PetscCall(VecDestroyVecs(itg->maxl, &itg->btilde));
     36:     PetscCall(VecDestroyVecs(itg->maxl, &itg->xtilde));
     37:     PetscCall(VecDestroy(&itg->guess));
     38:     PetscCall(VecDestroy(&itg->Ax));
     39:   }
     40:   if (itg->corr) PetscCall(PetscMemzero(itg->corr, sizeof(*itg->corr) * itg->maxl * itg->maxl));
     41:   itg->last_b       = NULL;
     42:   itg->last_b_state = 0;
     43:   if (itg->last_b_coefs) PetscCall(PetscMemzero(itg->last_b_coefs, sizeof(*itg->last_b_coefs) * itg->maxl));
     44:   PetscFunctionReturn(PETSC_SUCCESS);
     45: }

     47: static PetscErrorCode KSPGuessSetUp_Fischer(KSPGuess guess)
     48: {
     49:   KSPGuessFischer *itg = (KSPGuessFischer *)guess->data;

     51:   PetscFunctionBegin;
     52:   if (!itg->alpha) PetscCall(PetscMalloc1(itg->maxl, &itg->alpha));
     53:   if (!itg->xtilde) PetscCall(KSPCreateVecs(guess->ksp, itg->maxl, &itg->xtilde, 0, NULL));
     54:   if (!itg->btilde && (itg->method == 1 || itg->method == 3)) PetscCall(KSPCreateVecs(guess->ksp, itg->maxl, &itg->btilde, 0, NULL));
     55:   if (!itg->Ax && itg->method == 2) PetscCall(VecDuplicate(itg->xtilde[0], &itg->Ax));
     56:   if (!itg->guess && (itg->method == 1 || itg->method == 2)) PetscCall(VecDuplicate(itg->xtilde[0], &itg->guess));
     57:   if (!itg->corr && itg->method == 3) PetscCall(PetscCalloc1(itg->maxl * itg->maxl, &itg->corr));
     58:   if (!itg->last_b_coefs && itg->method == 3) PetscCall(PetscCalloc1(itg->maxl, &itg->last_b_coefs));
     59:   PetscFunctionReturn(PETSC_SUCCESS);
     60: }

     62: static PetscErrorCode KSPGuessDestroy_Fischer(KSPGuess guess)
     63: {
     64:   KSPGuessFischer *itg = (KSPGuessFischer *)guess->data;

     66:   PetscFunctionBegin;
     67:   PetscCall(PetscFree(itg->alpha));
     68:   PetscCall(VecDestroyVecs(itg->maxl, &itg->btilde));
     69:   PetscCall(VecDestroyVecs(itg->maxl, &itg->xtilde));
     70:   PetscCall(VecDestroy(&itg->guess));
     71:   PetscCall(VecDestroy(&itg->Ax));
     72:   PetscCall(PetscFree(itg->corr));
     73:   PetscCall(PetscFree(itg->last_b_coefs));
     74:   PetscCall(PetscFree(itg));
     75:   PetscCall(PetscObjectComposeFunction((PetscObject)guess, "KSPGuessFischerSetModel_C", NULL));
     76:   PetscFunctionReturn(PETSC_SUCCESS);
     77: }

     79: /* Note: do not change the b right-hand side as is done in the publication */
     80: static PetscErrorCode KSPGuessFormGuess_Fischer_1(KSPGuess guess, Vec b, Vec x)
     81: {
     82:   KSPGuessFischer *itg = (KSPGuessFischer *)guess->data;
     83:   PetscInt         i;

     85:   PetscFunctionBegin;
     86:   PetscCall(VecSet(x, 0.0));
     87:   PetscCall(VecMDot(b, itg->curl, itg->btilde, itg->alpha));
     88:   if (itg->monitor) {
     89:     PetscCall(PetscPrintf(((PetscObject)guess)->comm, "KSPFischerGuess alphas ="));
     90:     for (i = 0; i < itg->curl; i++) PetscCall(PetscPrintf(((PetscObject)guess)->comm, " %g", (double)PetscAbsScalar(itg->alpha[i])));
     91:     PetscCall(PetscPrintf(((PetscObject)guess)->comm, "\n"));
     92:   }
     93:   PetscCall(VecMAXPY(x, itg->curl, itg->alpha, itg->xtilde));
     94:   PetscCall(VecCopy(x, itg->guess));
     95:   PetscFunctionReturn(PETSC_SUCCESS);
     96: }

     98: static PetscErrorCode KSPGuessUpdate_Fischer_1(KSPGuess guess, Vec b, Vec x)
     99: {
    100:   KSPGuessFischer *itg = (KSPGuessFischer *)guess->data;
    101:   PetscReal        norm;
    102:   PetscInt         curl = itg->curl, i;

    104:   PetscFunctionBegin;
    105:   if (curl == itg->maxl) {
    106:     PetscCall(KSP_MatMult(guess->ksp, guess->A, x, itg->btilde[0]));
    107:     /* PetscCall(VecCopy(b,itg->btilde[0])); */
    108:     PetscCall(VecNormalize(itg->btilde[0], &norm));
    109:     PetscCall(VecCopy(x, itg->xtilde[0]));
    110:     PetscCall(VecScale(itg->xtilde[0], 1.0 / norm));
    111:     itg->curl = 1;
    112:   } else {
    113:     if (!curl) {
    114:       PetscCall(VecCopy(x, itg->xtilde[curl]));
    115:     } else {
    116:       PetscCall(VecWAXPY(itg->xtilde[curl], -1.0, itg->guess, x));
    117:     }
    118:     PetscCall(KSP_MatMult(guess->ksp, guess->A, itg->xtilde[curl], itg->btilde[curl]));
    119:     PetscCall(VecMDot(itg->btilde[curl], curl, itg->btilde, itg->alpha));
    120:     for (i = 0; i < curl; i++) itg->alpha[i] = -itg->alpha[i];
    121:     PetscCall(VecMAXPY(itg->btilde[curl], curl, itg->alpha, itg->btilde));
    122:     PetscCall(VecMAXPY(itg->xtilde[curl], curl, itg->alpha, itg->xtilde));
    123:     PetscCall(VecNormalize(itg->btilde[curl], &norm));
    124:     if (norm) {
    125:       PetscCall(VecScale(itg->xtilde[curl], 1.0 / norm));
    126:       itg->curl++;
    127:     } else {
    128:       PetscCall(PetscInfo(guess->ksp, "Not increasing dimension of Fischer space because new direction is identical to previous\n"));
    129:     }
    130:   }
    131:   PetscFunctionReturn(PETSC_SUCCESS);
    132: }

    134: /*
    135:   Given a basis generated already this computes a new guess x from the new right-hand side b
    136:   Figures out the components of b in each btilde direction and adds them to x
    137:   Note: do not change the b right-hand side as is done in the publication
    138: */
    139: static PetscErrorCode KSPGuessFormGuess_Fischer_2(KSPGuess guess, Vec b, Vec x)
    140: {
    141:   KSPGuessFischer *itg = (KSPGuessFischer *)guess->data;
    142:   PetscInt         i;

    144:   PetscFunctionBegin;
    145:   PetscCall(VecSet(x, 0.0));
    146:   PetscCall(VecMDot(b, itg->curl, itg->xtilde, itg->alpha));
    147:   if (itg->monitor) {
    148:     PetscCall(PetscPrintf(((PetscObject)guess)->comm, "KSPFischerGuess alphas ="));
    149:     for (i = 0; i < itg->curl; i++) PetscCall(PetscPrintf(((PetscObject)guess)->comm, " %g", (double)PetscAbsScalar(itg->alpha[i])));
    150:     PetscCall(PetscPrintf(((PetscObject)guess)->comm, "\n"));
    151:   }
    152:   PetscCall(VecMAXPY(x, itg->curl, itg->alpha, itg->xtilde));
    153:   PetscCall(VecCopy(x, itg->guess));
    154:   PetscFunctionReturn(PETSC_SUCCESS);
    155: }

    157: static PetscErrorCode KSPGuessUpdate_Fischer_2(KSPGuess guess, Vec b, Vec x)
    158: {
    159:   KSPGuessFischer *itg = (KSPGuessFischer *)guess->data;
    160:   PetscScalar      norm;
    161:   PetscInt         curl = itg->curl, i;

    163:   PetscFunctionBegin;
    164:   if (curl == itg->maxl) {
    165:     PetscCall(KSP_MatMult(guess->ksp, guess->A, x, itg->Ax)); /* norm = sqrt(x'Ax) */
    166:     PetscCall(VecDot(x, itg->Ax, &norm));
    167:     PetscCall(VecCopy(x, itg->xtilde[0]));
    168:     PetscCall(VecScale(itg->xtilde[0], 1.0 / PetscSqrtScalar(norm)));
    169:     itg->curl = 1;
    170:   } else {
    171:     if (!curl) {
    172:       PetscCall(VecCopy(x, itg->xtilde[curl]));
    173:     } else {
    174:       PetscCall(VecWAXPY(itg->xtilde[curl], -1.0, itg->guess, x));
    175:     }
    176:     PetscCall(KSP_MatMult(guess->ksp, guess->A, itg->xtilde[curl], itg->Ax));
    177:     PetscCall(VecMDot(itg->Ax, curl, itg->xtilde, itg->alpha));
    178:     for (i = 0; i < curl; i++) itg->alpha[i] = -itg->alpha[i];
    179:     PetscCall(VecMAXPY(itg->xtilde[curl], curl, itg->alpha, itg->xtilde));

    181:     PetscCall(KSP_MatMult(guess->ksp, guess->A, itg->xtilde[curl], itg->Ax)); /* norm = sqrt(xtilde[curl]'Axtilde[curl]) */
    182:     PetscCall(VecDot(itg->xtilde[curl], itg->Ax, &norm));
    183:     if (PetscAbsScalar(norm) != 0.0) {
    184:       PetscCall(VecScale(itg->xtilde[curl], 1.0 / PetscSqrtScalar(norm)));
    185:       itg->curl++;
    186:     } else {
    187:       PetscCall(PetscInfo(guess->ksp, "Not increasing dimension of Fischer space because new direction is identical to previous\n"));
    188:     }
    189:   }
    190:   PetscFunctionReturn(PETSC_SUCCESS);
    191: }

    193: /*
    194:   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.
    195: */
    196: static PetscErrorCode KSPGuessFormGuess_Fischer_3(KSPGuess guess, Vec b, Vec x)
    197: {
    198:   KSPGuessFischer *itg = (KSPGuessFischer *)guess->data;
    199:   PetscInt         i, j, m;
    200:   PetscReal       *s_values;
    201:   PetscScalar     *corr, *work, *scratch_vec, zero = 0.0, one = 1.0;
    202:   PetscBLASInt     blas_m, blas_info, blas_rank = 0, blas_lwork, blas_one = 1;
    203: #if defined(PETSC_USE_COMPLEX)
    204:   PetscReal *rwork;
    205: #endif

    207:   /* project provided b onto space of stored btildes */
    208:   PetscFunctionBegin;
    209:   PetscCall(VecSet(x, 0.0));
    210:   m           = itg->curl;
    211:   itg->last_b = b;
    212:   PetscCall(PetscObjectStateGet((PetscObject)b, &itg->last_b_state));
    213:   if (m > 0) {
    214:     PetscCall(PetscBLASIntCast(m, &blas_m));
    215:     blas_lwork = (/* assume a block size of m */ blas_m + 2) * blas_m;
    216: #if defined(PETSC_USE_COMPLEX)
    217:     PetscCall(PetscCalloc5(m * m, &corr, m, &s_values, blas_lwork, &work, 3 * m - 2, &rwork, m, &scratch_vec));
    218: #else
    219:     PetscCall(PetscCalloc4(m * m, &corr, m, &s_values, blas_lwork, &work, m, &scratch_vec));
    220: #endif
    221:     PetscCall(VecMDot(b, itg->curl, itg->btilde, itg->last_b_coefs));
    222:     for (j = 0; j < m; ++j) {
    223:       for (i = 0; i < m; ++i) corr[m * j + i] = itg->corr[(itg->maxl) * j + i];
    224:     }
    225:     PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
    226:     PetscReal max_s_value = 0.0;
    227: #if defined(PETSC_USE_COMPLEX)
    228:     PetscCallBLAS("LAPACKheev", LAPACKheev_("V", "L", &blas_m, corr, &blas_m, s_values, work, &blas_lwork, rwork, &blas_info));
    229: #else
    230:     PetscCallBLAS("LAPACKsyev", LAPACKsyev_("V", "L", &blas_m, corr, &blas_m, s_values, work, &blas_lwork, &blas_info));
    231: #endif

    233:     if (blas_info == 0) {
    234:       /* make corr store singular vectors and s_values store singular values */
    235:       for (j = 0; j < m; ++j) {
    236:         if (s_values[j] < 0.0) {
    237:           s_values[j] = PetscAbsReal(s_values[j]);
    238:           for (i = 0; i < m; ++i) corr[m * j + i] *= -1.0;
    239:         }
    240:         max_s_value = PetscMax(max_s_value, s_values[j]);
    241:       }

    243:       /* manually apply the action of the pseudoinverse */
    244:       PetscCallBLAS("BLASgemv", BLASgemv_("T", &blas_m, &blas_m, &one, corr, &blas_m, itg->last_b_coefs, &blas_one, &zero, scratch_vec, &blas_one));
    245:       for (j = 0; j < m; ++j) {
    246:         if (s_values[j] > itg->tol * max_s_value) {
    247:           scratch_vec[j] /= s_values[j];
    248:           blas_rank += 1;
    249:         } else {
    250:           scratch_vec[j] = 0.0;
    251:         }
    252:       }
    253:       PetscCallBLAS("BLASgemv", BLASgemv_("N", &blas_m, &blas_m, &one, corr, &blas_m, scratch_vec, &blas_one, &zero, itg->alpha, &blas_one));

    255:     } else {
    256:       PetscCall(PetscInfo(guess, "Warning eigenvalue solver failed with error code %d - setting initial guess to zero\n", (int)blas_info));
    257:       PetscCall(PetscMemzero(itg->alpha, sizeof(*itg->alpha) * itg->maxl));
    258:     }
    259:     PetscCall(PetscFPTrapPop());

    261:     if (itg->monitor && blas_info == 0) {
    262:       PetscCall(PetscPrintf(((PetscObject)guess)->comm, "KSPFischerGuess correlation rank = %d\n", (int)blas_rank));
    263:       PetscCall(PetscPrintf(((PetscObject)guess)->comm, "KSPFischerGuess singular values = "));
    264:       for (i = 0; i < itg->curl; i++) PetscCall(PetscPrintf(((PetscObject)guess)->comm, " %g", (double)s_values[i]));
    265:       PetscCall(PetscPrintf(((PetscObject)guess)->comm, "\n"));

    267:       PetscCall(PetscPrintf(((PetscObject)guess)->comm, "KSPFischerGuess alphas ="));
    268:       for (i = 0; i < itg->curl; i++) PetscCall(PetscPrintf(((PetscObject)guess)->comm, " %g", (double)PetscAbsScalar(itg->alpha[i])));
    269:       PetscCall(PetscPrintf(((PetscObject)guess)->comm, "\n"));
    270:     }
    271:     /* Form the initial guess by using b's projection coefficients with the xs */
    272:     PetscCall(VecMAXPY(x, itg->curl, itg->alpha, itg->xtilde));
    273: #if defined(PETSC_USE_COMPLEX)
    274:     PetscCall(PetscFree5(corr, s_values, work, rwork, scratch_vec));
    275: #else
    276:     PetscCall(PetscFree4(corr, s_values, work, scratch_vec));
    277: #endif
    278:   }
    279:   PetscFunctionReturn(PETSC_SUCCESS);
    280: }

    282: static PetscErrorCode KSPGuessUpdate_Fischer_3(KSPGuess guess, Vec b, Vec x)
    283: {
    284:   KSPGuessFischer *itg    = (KSPGuessFischer *)guess->data;
    285:   PetscBool        rotate = itg->curl == itg->maxl ? PETSC_TRUE : PETSC_FALSE;
    286:   PetscInt         i, j;
    287:   PetscObjectState b_state;
    288:   PetscScalar     *last_column;
    289:   Vec              oldest;

    291:   PetscFunctionBegin;
    292:   if (rotate) {
    293:     /* we have the maximum number of vectors so rotate: oldest vector is at index 0 */
    294:     oldest = itg->xtilde[0];
    295:     for (i = 1; i < itg->curl; ++i) itg->xtilde[i - 1] = itg->xtilde[i];
    296:     itg->xtilde[itg->curl - 1] = oldest;
    297:     PetscCall(VecCopy(x, itg->xtilde[itg->curl - 1]));

    299:     oldest = itg->btilde[0];
    300:     for (i = 1; i < itg->curl; ++i) itg->btilde[i - 1] = itg->btilde[i];
    301:     itg->btilde[itg->curl - 1] = oldest;
    302:     PetscCall(VecCopy(b, itg->btilde[itg->curl - 1]));
    303:     /* shift correlation matrix up and left */
    304:     for (j = 1; j < itg->maxl; ++j) {
    305:       for (i = 1; i < itg->maxl; ++i) itg->corr[(j - 1) * itg->maxl + i - 1] = itg->corr[j * itg->maxl + i];
    306:     }
    307:   } else {
    308:     /* append new vectors */
    309:     PetscCall(VecCopy(x, itg->xtilde[itg->curl]));
    310:     PetscCall(VecCopy(b, itg->btilde[itg->curl]));
    311:     itg->curl++;
    312:   }

    314:   /*
    315:       Populate new column of the correlation matrix and then copy it into the
    316:       row. itg->maxl is the allocated length per column: itg->curl is the actual
    317:       column length.
    318:       If possible reuse the dot products from FormGuess
    319:   */
    320:   last_column = itg->corr + (itg->curl - 1) * itg->maxl;
    321:   PetscCall(PetscObjectStateGet((PetscObject)b, &b_state));
    322:   if (b_state == itg->last_b_state && b == itg->last_b) {
    323:     if (rotate) {
    324:       for (i = 1; i < itg->maxl; ++i) itg->last_b_coefs[i - 1] = itg->last_b_coefs[i];
    325:     }
    326:     PetscCall(VecDot(b, b, &itg->last_b_coefs[itg->curl - 1]));
    327:     PetscCall(PetscArraycpy(last_column, itg->last_b_coefs, itg->curl));
    328:   } else {
    329:     PetscCall(VecMDot(b, itg->curl, itg->btilde, last_column));
    330:   }
    331:   for (i = 0; i < itg->curl; ++i) itg->corr[i * itg->maxl + itg->curl - 1] = last_column[i];
    332:   PetscFunctionReturn(PETSC_SUCCESS);
    333: }

    335: static PetscErrorCode KSPGuessSetFromOptions_Fischer(KSPGuess guess)
    336: {
    337:   KSPGuessFischer *ITG  = (KSPGuessFischer *)guess->data;
    338:   PetscInt         nmax = 2, model[2];
    339:   PetscBool        flg;

    341:   PetscFunctionBegin;
    342:   model[0] = ITG->method;
    343:   model[1] = ITG->maxl;
    344:   PetscOptionsBegin(PetscObjectComm((PetscObject)guess), ((PetscObject)guess)->prefix, "Fischer guess options", "KSPGuess");
    345:   PetscCall(PetscOptionsIntArray("-ksp_guess_fischer_model", "Model type and dimension of basis", "KSPGuessFischerSetModel", model, &nmax, &flg));
    346:   if (flg) PetscCall(KSPGuessFischerSetModel(guess, model[0], model[1]));
    347:   PetscCall(PetscOptionsReal("-ksp_guess_fischer_tol", "Tolerance to determine rank via ratio of singular values", "KSPGuessSetTolerance", ITG->tol, &ITG->tol, NULL));
    348:   PetscCall(PetscOptionsBool("-ksp_guess_fischer_monitor", "Monitor the guess", NULL, ITG->monitor, &ITG->monitor, NULL));
    349:   PetscOptionsEnd();
    350:   PetscFunctionReturn(PETSC_SUCCESS);
    351: }

    353: static PetscErrorCode KSPGuessSetTolerance_Fischer(KSPGuess guess, PetscReal tol)
    354: {
    355:   KSPGuessFischer *itg = (KSPGuessFischer *)guess->data;

    357:   PetscFunctionBegin;
    358:   itg->tol = tol;
    359:   PetscFunctionReturn(PETSC_SUCCESS);
    360: }

    362: static PetscErrorCode KSPGuessView_Fischer(KSPGuess guess, PetscViewer viewer)
    363: {
    364:   KSPGuessFischer *itg = (KSPGuessFischer *)guess->data;
    365:   PetscBool        isascii;

    367:   PetscFunctionBegin;
    368:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    369:   if (isascii) PetscCall(PetscViewerASCIIPrintf(viewer, "Model %" PetscInt_FMT ", size %" PetscInt_FMT "\n", itg->method, itg->maxl));
    370:   PetscFunctionReturn(PETSC_SUCCESS);
    371: }

    373: /*@
    374:   KSPGuessFischerSetModel - Set the Paul Fischer algorithm or its variants to compute the initial guess for a `KSPSolve()`

    376:   Logically Collective

    378:   Input Parameters:
    379: + guess - the initial guess context
    380: . model - use model 1, model 2, model 3, or any other number to turn it off
    381: - size  - size of subspace used to generate initial guess

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

    386:   Level: advanced

    388: .seealso: [](ch_ksp), `KSPGuess`, `KSPGuessCreate()`, `KSPSetUseFischerGuess()`, `KSPSetGuess()`, `KSPGetGuess()`, `KSP`
    389: @*/
    390: PetscErrorCode KSPGuessFischerSetModel(KSPGuess guess, PetscInt model, PetscInt size)
    391: {
    392:   PetscFunctionBegin;
    395:   PetscTryMethod(guess, "KSPGuessFischerSetModel_C", (KSPGuess, PetscInt, PetscInt), (guess, model, size));
    396:   PetscFunctionReturn(PETSC_SUCCESS);
    397: }

    399: static PetscErrorCode KSPGuessFischerSetModel_Fischer(KSPGuess guess, PetscInt model, PetscInt size)
    400: {
    401:   KSPGuessFischer *itg = (KSPGuessFischer *)guess->data;

    403:   PetscFunctionBegin;
    404:   if (model == 1) {
    405:     guess->ops->update    = KSPGuessUpdate_Fischer_1;
    406:     guess->ops->formguess = KSPGuessFormGuess_Fischer_1;
    407:   } else if (model == 2) {
    408:     guess->ops->update    = KSPGuessUpdate_Fischer_2;
    409:     guess->ops->formguess = KSPGuessFormGuess_Fischer_2;
    410:   } else if (model == 3) {
    411:     guess->ops->update    = KSPGuessUpdate_Fischer_3;
    412:     guess->ops->formguess = KSPGuessFormGuess_Fischer_3;
    413:   } else {
    414:     guess->ops->update    = NULL;
    415:     guess->ops->formguess = NULL;
    416:     itg->method           = 0;
    417:     PetscFunctionReturn(PETSC_SUCCESS);
    418:   }
    419:   if (size != itg->maxl) {
    420:     PetscCall(PetscFree(itg->alpha));
    421:     PetscCall(VecDestroyVecs(itg->maxl, &itg->btilde));
    422:     PetscCall(VecDestroyVecs(itg->maxl, &itg->xtilde));
    423:     PetscCall(VecDestroy(&itg->guess));
    424:     PetscCall(VecDestroy(&itg->Ax));
    425:   }
    426:   itg->method = model;
    427:   itg->maxl   = size;
    428:   PetscFunctionReturn(PETSC_SUCCESS);
    429: }

    431: /*MC
    432:     KSPGUESSFISCHER - Implements Paul Fischer's initial guess algorithms {cite}`fischer1998projection`
    433:     and a non-orthogonalizing variant for situations where a linear system is solved repeatedly

    435:     Level: intermediate

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

    444:     These are not intended to be used directly, they are called by `KSPSolve()` automatically with the command
    445:     line options `-ksp_guess_type fischer` `-ksp_guess_fischer_model <int,int>` or programmatically with
    446: .vb
    447:     KSPGetGuess(ksp,&guess);
    448:     KSPGuessSetType(guess,KSPGUESSFISCHER);
    449:     KSPGuessFischerSetModel(guess,model,basis);
    450:     KSPGuessSetTolerance(guess,PETSC_MACHINE_EPSILON);
    451: .ve
    452:     The default tolerance (which is only used in Method 3) is 32*`PETSC_MACHINE_EPSILON`. This value was chosen
    453:     empirically by trying a range of tolerances and picking the one that lowered the solver iteration count the most
    454:     with five vectors.

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

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

    462:     Developer Note:
    463:     The option `-ksp_fischer_guess <int,int>` is still available for backward compatibility

    465: .seealso: [](ch_ksp), `KSPGuess`, `KSPGuessType`, `KSP`
    466: M*/

    468: PetscErrorCode KSPGuessCreate_Fischer(KSPGuess guess)
    469: {
    470:   KSPGuessFischer *fischer;

    472:   PetscFunctionBegin;
    473:   PetscCall(PetscNew(&fischer));
    474:   fischer->method = 1; /* defaults to method 1 */
    475:   fischer->maxl   = 10;
    476:   fischer->tol    = 32.0 * PETSC_MACHINE_EPSILON;
    477:   guess->data     = fischer;

    479:   guess->ops->setfromoptions = KSPGuessSetFromOptions_Fischer;
    480:   guess->ops->destroy        = KSPGuessDestroy_Fischer;
    481:   guess->ops->settolerance   = KSPGuessSetTolerance_Fischer;
    482:   guess->ops->setup          = KSPGuessSetUp_Fischer;
    483:   guess->ops->view           = KSPGuessView_Fischer;
    484:   guess->ops->reset          = KSPGuessReset_Fischer;
    485:   guess->ops->update         = KSPGuessUpdate_Fischer_1;
    486:   guess->ops->formguess      = KSPGuessFormGuess_Fischer_1;

    488:   PetscCall(PetscObjectComposeFunction((PetscObject)guess, "KSPGuessFischerSetModel_C", KSPGuessFischerSetModel_Fischer));
    489:   PetscFunctionReturn(PETSC_SUCCESS);
    490: }