Actual source code: pod.c

      1: #include <petsc/private/kspimpl.h>
      2: #include <petsc/private/matimpl.h>
      3: #include <petscblaslapack.h>
      4: static PetscBool  cited      = PETSC_FALSE;
      5: static const char citation[] = "@phdthesis{zampini2010non,\n"
      6:                                "  title={Non-overlapping Domain Decomposition Methods for Cardiac Reaction-Diffusion Models and Applications},\n"
      7:                                "  author={Zampini, S},\n"
      8:                                "  year={2010},\n"
      9:                                "  school={PhD thesis, Universita degli Studi di Milano}\n"
     10:                                "}\n";

     12: typedef struct {
     13:   PetscInt      maxn;  /* maximum number of snapshots */
     14:   PetscInt      n;     /* number of active snapshots */
     15:   PetscInt      curr;  /* current tip of snapshots set */
     16:   Vec          *xsnap; /* snapshots */
     17:   Vec          *bsnap; /* rhs snapshots */
     18:   Vec          *work;  /* parallel work vectors */
     19:   PetscScalar  *dots_iallreduce;
     20:   MPI_Request   req_iallreduce;
     21:   PetscInt      ndots_iallreduce; /* if we have iallreduce we can hide the VecMDot communications */
     22:   PetscReal     tol;              /* relative tolerance to retain eigenvalues */
     23:   PetscBool     Aspd;             /* if true, uses the SPD operator as inner product */
     24:   PetscScalar  *corr;             /* correlation matrix */
     25:   PetscReal    *eigs;             /* eigenvalues */
     26:   PetscScalar  *eigv;             /* eigenvectors */
     27:   PetscBLASInt  nen;              /* dimension of lower dimensional system */
     28:   PetscInt      st;               /* first eigenvector of correlation matrix to be retained */
     29:   PetscBLASInt *iwork;            /* integer work vector */
     30:   PetscScalar  *yhay;             /* Y^H * A * Y */
     31:   PetscScalar  *low;              /* lower dimensional linear system */
     32: #if defined(PETSC_USE_COMPLEX)
     33:   PetscReal *rwork;
     34: #endif
     35:   PetscBLASInt lwork;
     36:   PetscScalar *swork;
     37:   PetscBool    monitor;
     38: } KSPGuessPOD;

     40: static PetscErrorCode KSPGuessReset_POD(KSPGuess guess)
     41: {
     42:   KSPGuessPOD *pod  = (KSPGuessPOD *)guess->data;
     43:   PetscLayout  Alay = NULL, vlay = NULL;
     44:   PetscBool    cong;

     46:   PetscFunctionBegin;
     47:   pod->nen  = 0;
     48:   pod->n    = 0;
     49:   pod->curr = 0;
     50:   /* need to wait for completion of outstanding requests */
     51:   if (pod->ndots_iallreduce) PetscCallMPI(MPI_Wait(&pod->req_iallreduce, MPI_STATUS_IGNORE));
     52:   pod->ndots_iallreduce = 0;
     53:   /* destroy vectors if the size of the linear system has changed */
     54:   if (guess->A) PetscCall(MatGetLayouts(guess->A, &Alay, NULL));
     55:   if (pod->xsnap) PetscCall(VecGetLayout(pod->xsnap[0], &vlay));
     56:   cong = PETSC_FALSE;
     57:   if (vlay && Alay) PetscCall(PetscLayoutCompare(Alay, vlay, &cong));
     58:   if (!cong) {
     59:     PetscCall(VecDestroyVecs(pod->maxn, &pod->xsnap));
     60:     PetscCall(VecDestroyVecs(pod->maxn, &pod->bsnap));
     61:     PetscCall(VecDestroyVecs(1, &pod->work));
     62:   }
     63:   PetscFunctionReturn(PETSC_SUCCESS);
     64: }

     66: static PetscErrorCode KSPGuessSetUp_POD(KSPGuess guess)
     67: {
     68:   KSPGuessPOD *pod = (KSPGuessPOD *)guess->data;

     70:   PetscFunctionBegin;
     71:   if (!pod->corr) {
     72:     PetscScalar  sdummy;
     73:     PetscReal    rdummy = 0;
     74:     PetscBLASInt bN, lierr, idummy;

     76:     PetscCall(PetscCalloc6(pod->maxn * pod->maxn, &pod->corr, pod->maxn, &pod->eigs, pod->maxn * pod->maxn, &pod->eigv, 6 * pod->maxn, &pod->iwork, pod->maxn * pod->maxn, &pod->yhay, pod->maxn * pod->maxn, &pod->low));
     77: #if defined(PETSC_USE_COMPLEX)
     78:     PetscCall(PetscMalloc1(7 * pod->maxn, &pod->rwork));
     79: #endif
     80: #if defined(PETSC_HAVE_MPI_NONBLOCKING_COLLECTIVES)
     81:     PetscCall(PetscMalloc1(3 * pod->maxn, &pod->dots_iallreduce));
     82: #endif
     83:     pod->lwork = -1;
     84:     PetscCall(PetscBLASIntCast(pod->maxn, &bN));
     85: #if !defined(PETSC_USE_COMPLEX)
     86:     PetscCallBLAS("LAPACKsyevx", LAPACKsyevx_("V", "A", "L", &bN, pod->corr, &bN, &rdummy, &rdummy, &idummy, &idummy, &rdummy, &idummy, pod->eigs, pod->eigv, &bN, &sdummy, &pod->lwork, pod->iwork, pod->iwork + 5 * bN, &lierr));
     87: #else
     88:     PetscCallBLAS("LAPACKsyevx", LAPACKsyevx_("V", "A", "L", &bN, pod->corr, &bN, &rdummy, &rdummy, &idummy, &idummy, &rdummy, &idummy, pod->eigs, pod->eigv, &bN, &sdummy, &pod->lwork, pod->rwork, pod->iwork, pod->iwork + 5 * bN, &lierr));
     89: #endif
     90:     PetscCheck(!lierr, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in query to SYEV Lapack routine %d", (int)lierr);
     91:     PetscCall(PetscBLASIntCast((PetscInt)PetscRealPart(sdummy), &pod->lwork));
     92:     PetscCall(PetscMalloc1(pod->lwork + PetscMax(bN * bN, 6 * bN), &pod->swork));
     93:   }
     94:   /* work vectors are sequential, we explicitly use MPI_Allreduce */
     95:   if (!pod->xsnap) {
     96:     Vec *v, vseq;

     98:     PetscCall(KSPCreateVecs(guess->ksp, 1, &v, 0, NULL));
     99:     PetscCall(VecCreateLocalVector(v[0], &vseq));
    100:     PetscCall(VecDestroyVecs(1, &v));
    101:     PetscCall(VecDuplicateVecs(vseq, pod->maxn, &pod->xsnap));
    102:     PetscCall(VecDestroy(&vseq));
    103:   }
    104:   if (!pod->bsnap) {
    105:     Vec *v, vseq;

    107:     PetscCall(KSPCreateVecs(guess->ksp, 0, NULL, 1, &v));
    108:     PetscCall(VecCreateLocalVector(v[0], &vseq));
    109:     PetscCall(VecDestroyVecs(1, &v));
    110:     PetscCall(VecDuplicateVecs(vseq, pod->maxn, &pod->bsnap));
    111:     PetscCall(VecDestroy(&vseq));
    112:   }
    113:   if (!pod->work) PetscCall(KSPCreateVecs(guess->ksp, 1, &pod->work, 0, NULL));
    114:   PetscFunctionReturn(PETSC_SUCCESS);
    115: }

    117: static PetscErrorCode KSPGuessDestroy_POD(KSPGuess guess)
    118: {
    119:   KSPGuessPOD *pod = (KSPGuessPOD *)guess->data;

    121:   PetscFunctionBegin;
    122:   PetscCall(PetscFree6(pod->corr, pod->eigs, pod->eigv, pod->iwork, pod->yhay, pod->low));
    123: #if defined(PETSC_USE_COMPLEX)
    124:   PetscCall(PetscFree(pod->rwork));
    125: #endif
    126:   /* need to wait for completion before destroying dots_iallreduce */
    127:   if (pod->ndots_iallreduce) PetscCallMPI(MPI_Wait(&pod->req_iallreduce, MPI_STATUS_IGNORE));
    128:   PetscCall(PetscFree(pod->dots_iallreduce));
    129:   PetscCall(PetscFree(pod->swork));
    130:   PetscCall(VecDestroyVecs(pod->maxn, &pod->bsnap));
    131:   PetscCall(VecDestroyVecs(pod->maxn, &pod->xsnap));
    132:   PetscCall(VecDestroyVecs(1, &pod->work));
    133:   PetscCall(PetscFree(pod));
    134:   PetscFunctionReturn(PETSC_SUCCESS);
    135: }

    137: static PetscErrorCode KSPGuessUpdate_POD(KSPGuess, Vec, Vec);

    139: static PetscErrorCode KSPGuessFormGuess_POD(KSPGuess guess, Vec b, Vec x)
    140: {
    141:   KSPGuessPOD *pod = (KSPGuessPOD *)guess->data;
    142:   PetscScalar  one = 1, zero = 0;
    143:   PetscBLASInt bN, ione      = 1, bNen, lierr;
    144:   PetscInt     i;

    146:   PetscFunctionBegin;
    147:   PetscCall(PetscCitationsRegister(citation, &cited));
    148:   if (pod->ndots_iallreduce) { /* complete communication and project the linear system */
    149:     PetscCall(KSPGuessUpdate_POD(guess, NULL, NULL));
    150:   }
    151:   if (!pod->nen) PetscFunctionReturn(PETSC_SUCCESS);
    152:   /* b_low = S * V^T * X^T * b */
    153:   PetscCall(VecGetLocalVectorRead(b, pod->bsnap[pod->curr]));
    154:   PetscCall(VecMDot(pod->bsnap[pod->curr], pod->n, pod->xsnap, pod->swork));
    155:   PetscCall(VecRestoreLocalVectorRead(b, pod->bsnap[pod->curr]));
    156:   PetscCallMPI(MPIU_Allreduce(pod->swork, pod->swork + pod->n, (PetscMPIInt)pod->n, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)guess)));
    157:   PetscCall(PetscBLASIntCast(pod->n, &bN));
    158:   PetscCall(PetscBLASIntCast(pod->nen, &bNen));
    159:   PetscCallBLAS("BLASgemv", BLASgemv_("T", &bN, &bNen, &one, pod->eigv + pod->st * pod->n, &bN, pod->swork + pod->n, &ione, &zero, pod->swork, &ione));
    160:   if (pod->monitor) {
    161:     PetscCall(PetscPrintf(PetscObjectComm((PetscObject)guess), "  KSPGuessPOD alphas = "));
    162:     for (i = 0; i < pod->nen; i++) {
    163: #if defined(PETSC_USE_COMPLEX)
    164:       PetscCall(PetscPrintf(PetscObjectComm((PetscObject)guess), "%g + %g i", (double)PetscRealPart(pod->swork[i]), (double)PetscImaginaryPart(pod->swork[i])));
    165: #else
    166:       PetscCall(PetscPrintf(PetscObjectComm((PetscObject)guess), "%g ", (double)pod->swork[i]));
    167: #endif
    168:     }
    169:     PetscCall(PetscPrintf(PetscObjectComm((PetscObject)guess), "\n"));
    170:   }
    171:   /* A_low x_low = b_low */
    172:   if (!pod->Aspd) { /* A is spd -> LOW = Identity */
    173:     KSP       pksp = guess->ksp;
    174:     PetscBool tsolve, symm, set;

    176:     if (pod->monitor) {
    177:       PetscMPIInt rank;
    178:       Mat         L;

    180:       PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)guess), &rank));
    181:       if (rank == 0) {
    182:         PetscCall(PetscPrintf(PetscObjectComm((PetscObject)guess), "  L = "));
    183:         PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, pod->nen, pod->nen, pod->low, &L));
    184:         PetscCall(MatView(L, NULL));
    185:         PetscCall(MatDestroy(&L));
    186:       }
    187:     }
    188:     PetscCall(MatIsSymmetricKnown(guess->A, &set, &symm));
    189:     tsolve = (set && symm) ? PETSC_FALSE : pksp->transpose_solve;
    190:     PetscCallBLAS("LAPACKgetrf", LAPACKgetrf_(&bNen, &bNen, pod->low, &bNen, pod->iwork, &lierr));
    191:     PetscCheck(!lierr, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in GETRF Lapack routine %d", (int)lierr);
    192:     PetscCallBLAS("LAPACKgetrs", LAPACKgetrs_(tsolve ? "T" : "N", &bNen, &ione, pod->low, &bNen, pod->iwork, pod->swork, &bNen, &lierr));
    193:     PetscCheck(!lierr, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in GETRS Lapack routine %d", (int)lierr);
    194:   }
    195:   /* x = X * V * S * x_low */
    196:   PetscCallBLAS("BLASgemv", BLASgemv_("N", &bN, &bNen, &one, pod->eigv + pod->st * pod->n, &bN, pod->swork, &ione, &zero, pod->swork + pod->n, &ione));
    197:   if (pod->monitor) {
    198:     PetscCall(PetscPrintf(PetscObjectComm((PetscObject)guess), "  KSPGuessPOD sol = "));
    199:     for (i = 0; i < pod->nen; i++) {
    200: #if defined(PETSC_USE_COMPLEX)
    201:       PetscCall(PetscPrintf(PetscObjectComm((PetscObject)guess), "%g + %g i", (double)PetscRealPart(pod->swork[i + pod->n]), (double)PetscImaginaryPart(pod->swork[i + pod->n])));
    202: #else
    203:       PetscCall(PetscPrintf(PetscObjectComm((PetscObject)guess), "%g ", (double)pod->swork[i + pod->n]));
    204: #endif
    205:     }
    206:     PetscCall(PetscPrintf(PetscObjectComm((PetscObject)guess), "\n"));
    207:   }
    208:   PetscCall(VecGetLocalVector(x, pod->bsnap[pod->curr]));
    209:   PetscCall(VecSet(pod->bsnap[pod->curr], 0));
    210:   PetscCall(VecMAXPY(pod->bsnap[pod->curr], pod->n, pod->swork + pod->n, pod->xsnap));
    211:   PetscCall(VecRestoreLocalVector(x, pod->bsnap[pod->curr]));
    212:   PetscFunctionReturn(PETSC_SUCCESS);
    213: }

    215: static PetscErrorCode KSPGuessUpdate_POD(KSPGuess guess, Vec b, Vec x)
    216: {
    217:   KSPGuessPOD *pod = (KSPGuessPOD *)guess->data;
    218:   PetscScalar  one = 1, zero = 0;
    219:   PetscReal    toten, parten, reps = 0; /* dlamch? */
    220:   PetscBLASInt bN, lierr, idummy;
    221:   PetscInt     i;

    223:   PetscFunctionBegin;
    224:   if (pod->ndots_iallreduce) goto complete_request;
    225:   pod->n = pod->n < pod->maxn ? pod->n + 1 : pod->maxn;
    226:   PetscCall(VecCopy(x, pod->xsnap[pod->curr]));
    227:   PetscCall(KSP_MatMult(guess->ksp, guess->A, x, pod->work[0]));
    228:   PetscCall(VecCopy(pod->work[0], pod->bsnap[pod->curr]));
    229:   if (pod->Aspd) {
    230:     PetscCall(VecMDot(pod->xsnap[pod->curr], pod->n, pod->bsnap, pod->swork));
    231: #if !defined(PETSC_HAVE_MPI_NONBLOCKING_COLLECTIVES)
    232:     PetscCallMPI(MPIU_Allreduce(pod->swork, pod->swork + 3 * pod->n, pod->n, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)guess)));
    233: #else
    234:     PetscCallMPI(MPI_Iallreduce(pod->swork, pod->dots_iallreduce, (PetscMPIInt)pod->n, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)guess), &pod->req_iallreduce));
    235:     pod->ndots_iallreduce = 1;
    236: #endif
    237:   } else {
    238:     PetscInt  off;
    239:     PetscBool set, herm;

    241: #if defined(PETSC_USE_COMPLEX)
    242:     PetscCall(MatIsHermitianKnown(guess->A, &set, &herm));
    243: #else
    244:     PetscCall(MatIsSymmetricKnown(guess->A, &set, &herm));
    245: #endif
    246:     off = (guess->ksp->transpose_solve && (!set || !herm)) ? 2 * pod->n : pod->n;

    248:     /* TODO: we may want to use a user-defined dot for the correlation matrix */
    249:     PetscCall(VecMDot(pod->xsnap[pod->curr], pod->n, pod->xsnap, pod->swork));
    250:     PetscCall(VecMDot(pod->bsnap[pod->curr], pod->n, pod->xsnap, pod->swork + off));
    251:     if (!set || !herm) {
    252:       off = (off == pod->n) ? 2 * pod->n : pod->n;
    253:       PetscCall(VecMDot(pod->xsnap[pod->curr], pod->n, pod->bsnap, pod->swork + off));
    254: #if !defined(PETSC_HAVE_MPI_NONBLOCKING_COLLECTIVES)
    255:       PetscCallMPI(MPIU_Allreduce(pod->swork, pod->swork + 3 * pod->n, (PetscMPIInt)(3 * pod->n), MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)guess)));
    256: #else
    257:       PetscCallMPI(MPI_Iallreduce(pod->swork, pod->dots_iallreduce, (PetscMPIInt)(3 * pod->n), MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)guess), &pod->req_iallreduce));
    258:       pod->ndots_iallreduce = 3;
    259: #endif
    260:     } else {
    261: #if !defined(PETSC_HAVE_MPI_NONBLOCKING_COLLECTIVES)
    262:       PetscCallMPI(MPIU_Allreduce(pod->swork, pod->swork + 3 * pod->n, (PetscMPIInt)(2 * pod->n), MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)guess)));
    263:       for (i = 0; i < pod->n; i++) pod->swork[5 * pod->n + i] = pod->swork[4 * pod->n + i];
    264: #else
    265:       PetscCallMPI(MPI_Iallreduce(pod->swork, pod->dots_iallreduce, (PetscMPIInt)(2 * pod->n), MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)guess), &pod->req_iallreduce));
    266:       pod->ndots_iallreduce = 2;
    267: #endif
    268:     }
    269:   }
    270:   if (pod->ndots_iallreduce) PetscFunctionReturn(PETSC_SUCCESS);

    272: complete_request:
    273:   if (pod->ndots_iallreduce) {
    274:     PetscCallMPI(MPI_Wait(&pod->req_iallreduce, MPI_STATUS_IGNORE));
    275:     switch (pod->ndots_iallreduce) {
    276:     case 3:
    277:       for (i = 0; i < pod->n; i++) pod->swork[3 * pod->n + i] = pod->dots_iallreduce[i];
    278:       for (i = 0; i < pod->n; i++) pod->swork[4 * pod->n + i] = pod->dots_iallreduce[pod->n + i];
    279:       for (i = 0; i < pod->n; i++) pod->swork[5 * pod->n + i] = pod->dots_iallreduce[2 * pod->n + i];
    280:       break;
    281:     case 2:
    282:       for (i = 0; i < pod->n; i++) pod->swork[3 * pod->n + i] = pod->dots_iallreduce[i];
    283:       for (i = 0; i < pod->n; i++) pod->swork[4 * pod->n + i] = pod->dots_iallreduce[pod->n + i];
    284:       for (i = 0; i < pod->n; i++) pod->swork[5 * pod->n + i] = pod->dots_iallreduce[pod->n + i];
    285:       break;
    286:     case 1:
    287:       for (i = 0; i < pod->n; i++) pod->swork[3 * pod->n + i] = pod->dots_iallreduce[i];
    288:       break;
    289:     default:
    290:       SETERRQ(PetscObjectComm((PetscObject)guess), PETSC_ERR_PLIB, "Invalid number of outstanding dots operations: %" PetscInt_FMT, pod->ndots_iallreduce);
    291:     }
    292:   }
    293:   pod->ndots_iallreduce = 0;

    295:   /* correlation matrix and Y^H A Y (Galerkin) */
    296:   for (i = 0; i < pod->n; i++) {
    297:     pod->corr[pod->curr * pod->maxn + i] = pod->swork[3 * pod->n + i];
    298:     pod->corr[i * pod->maxn + pod->curr] = PetscConj(pod->swork[3 * pod->n + i]);
    299:     if (!pod->Aspd) {
    300:       pod->yhay[pod->curr * pod->maxn + i] = pod->swork[4 * pod->n + i];
    301:       pod->yhay[i * pod->maxn + pod->curr] = PetscConj(pod->swork[5 * pod->n + i]);
    302:     }
    303:   }
    304:   /* syevx changes the input matrix */
    305:   for (i = 0; i < pod->n; i++) {
    306:     PetscInt j;
    307:     for (j = i; j < pod->n; j++) pod->swork[i * pod->n + j] = pod->corr[i * pod->maxn + j];
    308:   }
    309:   PetscCall(PetscBLASIntCast(pod->n, &bN));
    310: #if !defined(PETSC_USE_COMPLEX)
    311:   PetscCallBLAS("LAPACKsyevx", LAPACKsyevx_("V", "A", "L", &bN, pod->swork, &bN, &reps, &reps, &idummy, &idummy, &reps, &idummy, pod->eigs, pod->eigv, &bN, pod->swork + bN * bN, &pod->lwork, pod->iwork, pod->iwork + 5 * bN, &lierr));
    312: #else
    313:   PetscCallBLAS("LAPACKsyevx", LAPACKsyevx_("V", "A", "L", &bN, pod->swork, &bN, &reps, &reps, &idummy, &idummy, &reps, &idummy, pod->eigs, pod->eigv, &bN, pod->swork + bN * bN, &pod->lwork, pod->rwork, pod->iwork, pod->iwork + 5 * bN, &lierr));
    314: #endif
    315:   PetscCheck(lierr >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in SYEV Lapack routine: illegal argument %d", -(int)lierr);
    316:   PetscCheck(!lierr, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in SYEV Lapack routine: %d eigenvectors failed to converge", (int)lierr);

    318:   /* dimension of lower dimensional system */
    319:   pod->st = -1;
    320:   for (i = 0, toten = 0; i < pod->n; i++) {
    321:     pod->eigs[i] = PetscMax(pod->eigs[i], 0.0);
    322:     toten += pod->eigs[i];
    323:     if (!pod->eigs[i]) pod->st = i;
    324:   }
    325:   pod->nen = 0;
    326:   for (i = pod->n - 1, parten = 0; i > pod->st && toten > 0; i--) {
    327:     pod->nen++;
    328:     parten += pod->eigs[i];
    329:     if (parten + toten * pod->tol >= toten) break;
    330:   }
    331:   pod->st = pod->n - pod->nen;

    333:   /* Compute eigv = V * S */
    334:   for (i = pod->st; i < pod->n; i++) {
    335:     const PetscReal v  = 1.0 / PetscSqrtReal(pod->eigs[i]);
    336:     const PetscInt  st = pod->n * i;
    337:     PetscInt        j;

    339:     for (j = 0; j < pod->n; j++) pod->eigv[st + j] *= v;
    340:   }

    342:   /* compute S * V^T * X^T * A * X * V * S if needed */
    343:   if (pod->nen && !pod->Aspd) {
    344:     PetscBLASInt bNen, bMaxN;
    345:     PetscInt     st = pod->st * pod->n;
    346:     PetscCall(PetscBLASIntCast(pod->nen, &bNen));
    347:     PetscCall(PetscBLASIntCast(pod->maxn, &bMaxN));
    348:     PetscCallBLAS("BLASgemm", BLASgemm_("T", "N", &bNen, &bN, &bN, &one, pod->eigv + st, &bN, pod->yhay, &bMaxN, &zero, pod->swork, &bNen));
    349:     PetscCallBLAS("BLASgemm", BLASgemm_("N", "N", &bNen, &bNen, &bN, &one, pod->swork, &bNen, pod->eigv + st, &bN, &zero, pod->low, &bNen));
    350:   }

    352:   if (pod->monitor) {
    353:     PetscMPIInt rank;
    354:     Mat         C;

    356:     PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)guess), &rank));
    357:     if (rank == 0) {
    358:       PetscCall(PetscPrintf(PetscObjectComm((PetscObject)guess), "  C = "));
    359:       PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, pod->n, pod->n, pod->corr, &C));
    360:       PetscCall(MatDenseSetLDA(C, pod->maxn));
    361:       PetscCall(MatView(C, NULL));
    362:       PetscCall(MatDestroy(&C));
    363:       PetscCall(PetscPrintf(PetscObjectComm((PetscObject)guess), "  YHAY = "));
    364:       PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, pod->n, pod->n, pod->yhay, &C));
    365:       PetscCall(MatDenseSetLDA(C, pod->maxn));
    366:       PetscCall(MatView(C, NULL));
    367:       PetscCall(MatDestroy(&C));
    368:     }
    369:     PetscCall(PetscPrintf(PetscObjectComm((PetscObject)guess), "  KSPGuessPOD: basis %" PetscBLASInt_FMT ", energy fractions = ", pod->nen));
    370:     for (i = pod->n - 1; i >= 0; i--) PetscCall(PetscPrintf(PetscObjectComm((PetscObject)guess), "%1.6e (%d) ", (double)(pod->eigs[i] / toten), i >= pod->st ? 1 : 0));
    371:     PetscCall(PetscPrintf(PetscObjectComm((PetscObject)guess), "\n"));
    372:     if (PetscDefined(USE_DEBUG)) {
    373:       for (i = 0; i < pod->n; i++) {
    374:         Vec          v;
    375:         PetscInt     j;
    376:         PetscBLASInt bNen, ione = 1;

    378:         PetscCall(VecDuplicate(pod->xsnap[i], &v));
    379:         PetscCall(VecCopy(pod->xsnap[i], v));
    380:         PetscCall(PetscBLASIntCast(pod->nen, &bNen));
    381:         PetscCallBLAS("BLASgemv", BLASgemv_("T", &bN, &bNen, &one, pod->eigv + pod->st * pod->n, &bN, pod->corr + pod->maxn * i, &ione, &zero, pod->swork, &ione));
    382:         PetscCallBLAS("BLASgemv", BLASgemv_("N", &bN, &bNen, &one, pod->eigv + pod->st * pod->n, &bN, pod->swork, &ione, &zero, pod->swork + pod->n, &ione));
    383:         for (j = 0; j < pod->n; j++) pod->swork[j] = -pod->swork[pod->n + j];
    384:         PetscCall(VecMAXPY(v, pod->n, pod->swork, pod->xsnap));
    385:         PetscCall(VecDot(v, v, pod->swork));
    386:         PetscCallMPI(MPIU_Allreduce(pod->swork, pod->swork + 1, 1, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)guess)));
    387:         PetscCall(PetscPrintf(PetscObjectComm((PetscObject)guess), "  Error projection %" PetscInt_FMT ": %g (expected lower than %g)\n", i, (double)PetscRealPart(pod->swork[1]), (double)(toten - parten)));
    388:         PetscCall(VecDestroy(&v));
    389:       }
    390:     }
    391:   }
    392:   /* new tip */
    393:   pod->curr = (pod->curr + 1) % pod->maxn;
    394:   PetscFunctionReturn(PETSC_SUCCESS);
    395: }

    397: static PetscErrorCode KSPGuessSetFromOptions_POD(KSPGuess guess)
    398: {
    399:   KSPGuessPOD *pod = (KSPGuessPOD *)guess->data;

    401:   PetscFunctionBegin;
    402:   PetscOptionsBegin(PetscObjectComm((PetscObject)guess), ((PetscObject)guess)->prefix, "POD initial guess options", "KSPGuess");
    403:   PetscCall(PetscOptionsInt("-ksp_guess_pod_size", "Number of snapshots", NULL, pod->maxn, &pod->maxn, NULL));
    404:   PetscCall(PetscOptionsBool("-ksp_guess_pod_monitor", "Monitor initial guess generator", NULL, pod->monitor, &pod->monitor, NULL));
    405:   PetscCall(PetscOptionsReal("-ksp_guess_pod_tol", "Tolerance to retain eigenvectors", "KSPGuessSetTolerance", pod->tol, &pod->tol, NULL));
    406:   PetscCall(PetscOptionsBool("-ksp_guess_pod_Ainner", "Use the operator as inner product (must be SPD)", NULL, pod->Aspd, &pod->Aspd, NULL));
    407:   PetscOptionsEnd();
    408:   PetscFunctionReturn(PETSC_SUCCESS);
    409: }

    411: static PetscErrorCode KSPGuessSetTolerance_POD(KSPGuess guess, PetscReal tol)
    412: {
    413:   KSPGuessPOD *pod = (KSPGuessPOD *)guess->data;

    415:   PetscFunctionBegin;
    416:   pod->tol = tol;
    417:   PetscFunctionReturn(PETSC_SUCCESS);
    418: }

    420: static PetscErrorCode KSPGuessView_POD(KSPGuess guess, PetscViewer viewer)
    421: {
    422:   KSPGuessPOD *pod = (KSPGuessPOD *)guess->data;
    423:   PetscBool    isascii;

    425:   PetscFunctionBegin;
    426:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    427:   if (isascii) PetscCall(PetscViewerASCIIPrintf(viewer, "Max size %" PetscInt_FMT ", tolerance %g, Ainner %d\n", pod->maxn, (double)pod->tol, pod->Aspd));
    428:   PetscFunctionReturn(PETSC_SUCCESS);
    429: }

    431: /*MC
    432:     KSPGUESSPOD - Implements a proper orthogonal decomposition based Galerkin scheme for repeated linear system solves.

    434:   Options Database Keys:
    435: +  -ksp_guess_pod_size <size> - Number of snapshots
    436: .  -ksp_guess_pod_monitor <true or false> - Monitor initial guess generator
    437: .  -ksp_guess_pod_tol <tol> - Tolerance to retain eigenvectors
    438: -  -ksp_guess_pod_Ainner <true> - Use the operator as inner product (must be SPD)

    440:   Level: intermediate

    442:   Note:
    443:   The initial guess is obtained by solving a small and dense linear system, obtained by Galerkin projection on a lower dimensional space generated by the previous solutions as presented in {cite}`volkwein2013proper`.

    445: .seealso: [](ch_ksp), `KSPGuess`, `KSPGuessType`, `KSPGuessCreate()`, `KSPSetGuess()`, `KSPGetGuess()`
    446: M*/
    447: PetscErrorCode KSPGuessCreate_POD(KSPGuess guess)
    448: {
    449:   KSPGuessPOD *pod;

    451:   PetscFunctionBegin;
    452:   PetscCall(PetscNew(&pod));
    453:   pod->maxn   = 10;
    454:   pod->tol    = PETSC_MACHINE_EPSILON;
    455:   guess->data = pod;

    457:   guess->ops->setfromoptions = KSPGuessSetFromOptions_POD;
    458:   guess->ops->destroy        = KSPGuessDestroy_POD;
    459:   guess->ops->settolerance   = KSPGuessSetTolerance_POD;
    460:   guess->ops->setup          = KSPGuessSetUp_POD;
    461:   guess->ops->view           = KSPGuessView_POD;
    462:   guess->ops->reset          = KSPGuessReset_POD;
    463:   guess->ops->update         = KSPGuessUpdate_POD;
    464:   guess->ops->formguess      = KSPGuessFormGuess_POD;
    465:   PetscFunctionReturn(PETSC_SUCCESS);
    466: }