Actual source code: hpddm.cxx

      1: #include <petsc/private/petschpddm.h>

      3: /* static array length */
      4: #define ALEN(a) (sizeof(a)/sizeof((a)[0]))

      6: const char *const KSPHPDDMTypes[]           = { KSPGMRES, "bgmres", KSPCG, "bcg", "gcrodr", "bgcrodr", "bfbcg", KSPPREONLY };
      7: static const char *HPDDMOrthogonalization[] = { "cgs", "mgs" };
      8: static const char *HPDDMQR[]                = { "cholqr", "cgs", "mgs" };
      9: static const char *HPDDMVariant[]           = { "left", "right", "flexible" };
     10: static const char *HPDDMRecycleTarget[]     = { "SM", "LM", "SR", "LR", "SI", "LI" };
     11: static const char *HPDDMRecycleStrategy[]   = { "A", "B" };

     13: PetscBool HPDDMCite = PETSC_FALSE;
     14: const char HPDDMCitation[] = "@article{jolivet2020petsc,\n"
     15: "  Author = {Jolivet, Pierre and Roman, Jose E. and Zampini, Stefano},\n"
     16: "  Title = {{KSPHPDDM} and {PCHPDDM}: Extending {PETSc} with Robust Overlapping {Schwarz} Preconditioners and Advanced {Krylov} Methods},\n"
     17: "  Year = {2021},\n"
     18: "  Publisher = {Elsevier},\n"
     19: "  Journal = {Computer \\& Mathematics with Applications},\n"
     20: "  Volume = {84},\n"
     21: "  Pages = {277--295},\n"
     22: "  Url = {https://github.com/prj-/jolivet2020petsc}\n"
     23: "}\n";

     25: #if defined(PETSC_HAVE_SLEPC) && defined(PETSC_USE_SHARED_LIBRARIES)
     26: static PetscBool loadedDL = PETSC_FALSE;
     27: #endif

     29: static PetscErrorCode KSPSetFromOptions_HPDDM(PetscOptionItems *PetscOptionsObject, KSP ksp)
     30: {
     31:   KSP_HPDDM      *data = (KSP_HPDDM*)ksp->data;
     32:   PetscInt       i, j;
     33:   PetscMPIInt    size;

     37:   PetscOptionsHead(PetscOptionsObject, "KSPHPDDM options, cf. https://github.com/hpddm/hpddm");
     38:   i = (data->cntl[0] == static_cast<char>(PETSC_DECIDE) ? HPDDM_KRYLOV_METHOD_GMRES : data->cntl[0]);
     39:   PetscOptionsEList("-ksp_hpddm_type", "Type of Krylov method", "KSPHPDDM", KSPHPDDMTypes, ALEN(KSPHPDDMTypes), KSPHPDDMTypes[HPDDM_KRYLOV_METHOD_GMRES], &i, NULL);
     40:   if (i == ALEN(KSPHPDDMTypes) - 1)
     41:     i = HPDDM_KRYLOV_METHOD_NONE; /* need to shift the value since HPDDM_KRYLOV_METHOD_RICHARDSON is not registered in PETSc */
     42:   data->cntl[0] = i;
     43:   if (data->cntl[0] != HPDDM_KRYLOV_METHOD_NONE) {
     44:     if (data->cntl[0] != HPDDM_KRYLOV_METHOD_BCG && data->cntl[0] != HPDDM_KRYLOV_METHOD_BFBCG) {
     45:       i = (data->cntl[1] == static_cast<char>(PETSC_DECIDE) ? HPDDM_VARIANT_LEFT : data->cntl[1]);
     46:       if (ksp->pc_side_set == PC_SIDE_DEFAULT) {
     47:         PetscOptionsEList("-ksp_hpddm_variant", "Left, right, or variable preconditioning", "KSPHPDDM", HPDDMVariant, ALEN(HPDDMVariant), HPDDMVariant[HPDDM_VARIANT_LEFT], &i, NULL);
     48:       } else if (ksp->pc_side_set == PC_RIGHT) i = HPDDM_VARIANT_RIGHT;
     49:       else if (ksp->pc_side_set == PC_SYMMETRIC) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Symmetric preconditioning not implemented");
     50:       data->cntl[1] = i;
     51:       if (i > 0) {
     52:         KSPSetPCSide(ksp, PC_RIGHT);
     53:       }
     54:     }
     55:     if (data->cntl[0] == HPDDM_KRYLOV_METHOD_BGMRES || data->cntl[0] == HPDDM_KRYLOV_METHOD_BGCRODR || data->cntl[0] == HPDDM_KRYLOV_METHOD_BFBCG) {
     56:       data->rcntl[0] = (std::abs(data->rcntl[0] - static_cast<PetscReal>(PETSC_DECIDE)) < PETSC_SMALL ? -1.0 : data->rcntl[0]);
     57:       PetscOptionsReal("-ksp_hpddm_deflation_tol", "Tolerance when deflating right-hand sides inside block methods", "KSPHPDDM", data->rcntl[0], data->rcntl, NULL);
     58:       i = (data->scntl[data->cntl[0] != HPDDM_KRYLOV_METHOD_BFBCG] == static_cast<unsigned short>(PETSC_DECIDE) ? 1 : PetscMax(1, data->scntl[data->cntl[0] != HPDDM_KRYLOV_METHOD_BFBCG]));
     59:       PetscOptionsRangeInt("-ksp_hpddm_enlarge_krylov_subspace", "Split the initial right-hand side into multiple vectors", "KSPHPDDM", i, &i, NULL, 1, std::numeric_limits<unsigned short>::max() - 1);
     60:       data->scntl[data->cntl[0] != HPDDM_KRYLOV_METHOD_BFBCG] = i;
     61:     } else data->scntl[data->cntl[0] != HPDDM_KRYLOV_METHOD_BCG] = 0;
     62:     if (data->cntl[0] == HPDDM_KRYLOV_METHOD_GMRES || data->cntl[0] == HPDDM_KRYLOV_METHOD_BGMRES || data->cntl[0] == HPDDM_KRYLOV_METHOD_GCRODR || data->cntl[0] == HPDDM_KRYLOV_METHOD_BGCRODR) {
     63:       i = (data->cntl[2] == static_cast<char>(PETSC_DECIDE) ? HPDDM_ORTHOGONALIZATION_CGS : data->cntl[2] & 3);
     64:       PetscOptionsEList("-ksp_hpddm_orthogonalization", "Classical (faster) or Modified (more robust) Gram--Schmidt process", "KSPHPDDM", HPDDMOrthogonalization, ALEN(HPDDMOrthogonalization), HPDDMOrthogonalization[HPDDM_ORTHOGONALIZATION_CGS], &i, NULL);
     65:       j = (data->cntl[2] == static_cast<char>(PETSC_DECIDE) ? HPDDM_QR_CHOLQR : ((data->cntl[2] >> 2) & 7));
     66:       PetscOptionsEList("-ksp_hpddm_qr", "Distributed QR factorizations computed with Cholesky QR, Classical or Modified Gram--Schmidt process", "KSPHPDDM", HPDDMQR, ALEN(HPDDMQR), HPDDMQR[HPDDM_QR_CHOLQR], &j, NULL);
     67:       data->cntl[2] = static_cast<char>(i) + (static_cast<char>(j) << 2);
     68:       i = (data->scntl[0] == static_cast<unsigned short>(PETSC_DECIDE) ? PetscMin(30, ksp->max_it) : data->scntl[0]);
     69:       PetscOptionsRangeInt("-ksp_gmres_restart", "Maximum number of Arnoldi vectors generated per cycle", "KSPHPDDM", i, &i, NULL, PetscMin(1, ksp->max_it), PetscMin(ksp->max_it, std::numeric_limits<unsigned short>::max() - 1));
     70:       data->scntl[0] = i;
     71:     }
     72:     if (data->cntl[0] == HPDDM_KRYLOV_METHOD_BCG || data->cntl[0] == HPDDM_KRYLOV_METHOD_BFBCG) {
     73:       j = (data->cntl[1] == static_cast<char>(PETSC_DECIDE) ? HPDDM_QR_CHOLQR : data->cntl[1]);
     74:       PetscOptionsEList("-ksp_hpddm_qr", "Distributed QR factorizations computed with Cholesky QR, Classical or Modified Gram--Schmidt process", "KSPHPDDM", HPDDMQR, ALEN(HPDDMQR), HPDDMQR[HPDDM_QR_CHOLQR], &j, NULL);
     75:       data->cntl[1] = j;
     76:     }
     77:     if (data->cntl[0] == HPDDM_KRYLOV_METHOD_GCRODR || data->cntl[0] == HPDDM_KRYLOV_METHOD_BGCRODR) {
     78:       i = (data->icntl[0] == static_cast<int>(PETSC_DECIDE) ? PetscMin(20, data->scntl[0] - 1) : data->icntl[0]);
     79:       PetscOptionsRangeInt("-ksp_hpddm_recycle", "Number of harmonic Ritz vectors to compute", "KSPHPDDM", i, &i, NULL, 1, data->scntl[0] - 1);
     80:       data->icntl[0] = i;
     81:       if (!PetscDefined(HAVE_SLEPC) || !PetscDefined(USE_SHARED_LIBRARIES) || data->cntl[0] == HPDDM_KRYLOV_METHOD_GCRODR) {
     82:         i = (data->cntl[3] == static_cast<char>(PETSC_DECIDE) ? HPDDM_RECYCLE_TARGET_SM : data->cntl[3]);
     83:         PetscOptionsEList("-ksp_hpddm_recycle_target", "Criterion to select harmonic Ritz vectors", "KSPHPDDM", HPDDMRecycleTarget, ALEN(HPDDMRecycleTarget), HPDDMRecycleTarget[HPDDM_RECYCLE_TARGET_SM], &i, NULL);
     84:         data->cntl[3] = i;
     85:       } else {
     86:         MPI_Comm_size(PetscObjectComm((PetscObject)ksp), &size);
     87:         i = (data->cntl[3] == static_cast<char>(PETSC_DECIDE) ? 1 : data->cntl[3]);
     88:         PetscOptionsRangeInt("-ksp_hpddm_recycle_redistribute", "Number of processes used to solve eigenvalue problems when recycling in BGCRODR", "KSPHPDDM", i, &i, NULL, 1, PetscMin(size, 192));
     89:         data->cntl[3] = i;
     90:       }
     91:       i = (data->cntl[4] == static_cast<char>(PETSC_DECIDE) ? HPDDM_RECYCLE_STRATEGY_A : data->cntl[4]);
     92:       PetscOptionsEList("-ksp_hpddm_recycle_strategy", "Generalized eigenvalue problem to solve for recycling", "KSPHPDDM", HPDDMRecycleStrategy, ALEN(HPDDMRecycleStrategy), HPDDMRecycleStrategy[HPDDM_RECYCLE_STRATEGY_A], &i, NULL);
     93:       data->cntl[4] = i;
     94:     }
     95:   } else {
     96:     data->cntl[0] = HPDDM_KRYLOV_METHOD_NONE;
     97:     data->scntl[1] = 1;
     98:   }
     99:   if (ksp->nmax > std::numeric_limits<int>::max() || ksp->nmax < std::numeric_limits<int>::min()) SETERRQ1(PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "KSPMatSolve() block size %D not representable by an integer, which is not handled by KSPHPDDM", ksp->nmax);
    100:   else data->icntl[1] = static_cast<int>(ksp->nmax);
    101:   PetscOptionsTail();
    102:   return(0);
    103: }

    105: static PetscErrorCode KSPView_HPDDM(KSP ksp, PetscViewer viewer)
    106: {
    107:   KSP_HPDDM            *data = (KSP_HPDDM*)ksp->data;
    108:   HPDDM::PETScOperator *op = data->op;
    109:   const PetscScalar    *array = op ? op->storage() : NULL;
    110:   PetscBool            ascii;
    111:   PetscErrorCode       ierr;

    114:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &ascii);
    115:   if (op && ascii) {
    116:     PetscViewerASCIIPrintf(viewer, "HPDDM type: %s\n", KSPHPDDMTypes[std::min(static_cast<PetscInt>(data->cntl[0]), static_cast<PetscInt>(ALEN(KSPHPDDMTypes) - 1))]);
    117:     if (data->cntl[0] == HPDDM_KRYLOV_METHOD_BGMRES || data->cntl[0] == HPDDM_KRYLOV_METHOD_BGCRODR || data->cntl[0] == HPDDM_KRYLOV_METHOD_BFBCG) {
    118:       if (std::abs(data->rcntl[0] - static_cast<PetscReal>(PETSC_DECIDE)) < PETSC_SMALL) {
    119:         PetscViewerASCIIPrintf(viewer, "no deflation at restarts\n", PetscBools[array ? PETSC_TRUE : PETSC_FALSE]);
    120:       } else {
    121:         PetscViewerASCIIPrintf(viewer, "deflation tolerance: %g\n", data->rcntl[0]);
    122:       }
    123:     }
    124:     if (data->cntl[0] == HPDDM_KRYLOV_METHOD_GCRODR || data->cntl[0] == HPDDM_KRYLOV_METHOD_BGCRODR) {
    125:       PetscViewerASCIIPrintf(viewer, "deflation subspace attached? %s\n", PetscBools[array ? PETSC_TRUE : PETSC_FALSE]);
    126:       if (!PetscDefined(HAVE_SLEPC) || !PetscDefined(USE_SHARED_LIBRARIES) || data->cntl[0] == HPDDM_KRYLOV_METHOD_GCRODR) {
    127:         PetscViewerASCIIPrintf(viewer, "deflation target: %s\n", HPDDMRecycleTarget[static_cast<PetscInt>(data->cntl[3])]);
    128:       } else {
    129:         PetscViewerASCIIPrintf(viewer, "redistribution size: %D\n", static_cast<PetscInt>(data->cntl[3]));
    130:       }
    131:     }
    132:     if (data->icntl[1] != static_cast<int>(PETSC_DECIDE)) {
    133:       PetscViewerASCIIPrintf(viewer, "  block size is %d\n", data->icntl[1]);
    134:     }
    135:   }
    136:   return(0);
    137: }

    139: static PetscErrorCode KSPSetUp_HPDDM(KSP ksp)
    140: {
    141:   KSP_HPDDM      *data = (KSP_HPDDM*)ksp->data;
    142:   Mat            A;
    143:   PetscInt       n, bs;
    144:   PetscBool      match;

    148:   KSPGetOperators(ksp, &A, NULL);
    149:   MatGetLocalSize(A, &n, NULL);
    150:   MatGetBlockSize(A, &bs);
    151:   PetscObjectTypeCompareAny((PetscObject)A, &match, MATSEQKAIJ, MATMPIKAIJ, "");
    152:   if (match) n /= bs;
    153:   data->op = new HPDDM::PETScOperator(ksp, n);
    154:   if (PetscUnlikely(!ksp->setfromoptionscalled || data->cntl[0] == static_cast<char>(PETSC_DECIDE))) { /* what follows is basically a copy/paste of KSPSetFromOptions_HPDDM, with no call to PetscOptions() */
    155:     PetscInfo(ksp, "KSPSetFromOptions() not called or uninitialized internal structure, hardwiring default KSPHPDDM options\n");
    156:     if (data->cntl[0] == static_cast<char>(PETSC_DECIDE))
    157:       data->cntl[0] = 0; /* GMRES by default */
    158:     if (data->cntl[0] != HPDDM_KRYLOV_METHOD_NONE) { /* following options do not matter with PREONLY */
    159:       if (data->cntl[0] != HPDDM_KRYLOV_METHOD_BCG && data->cntl[0] != HPDDM_KRYLOV_METHOD_BFBCG) {
    160:         data->cntl[1] = HPDDM_VARIANT_LEFT; /* left preconditioning by default */
    161:         if (ksp->pc_side_set == PC_RIGHT) data->cntl[1] = HPDDM_VARIANT_RIGHT;
    162:         else if (ksp->pc_side_set == PC_SYMMETRIC) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Symmetric preconditioning not implemented");
    163:         if (data->cntl[1] > 0) {
    164:           KSPSetPCSide(ksp, PC_RIGHT);
    165:         }
    166:       }
    167:       if (data->cntl[0] == HPDDM_KRYLOV_METHOD_BGMRES || data->cntl[0] == HPDDM_KRYLOV_METHOD_BGCRODR || data->cntl[0] == HPDDM_KRYLOV_METHOD_BFBCG) {
    168:         data->rcntl[0] = -1.0; /* no deflation by default */
    169:         data->scntl[data->cntl[0] != HPDDM_KRYLOV_METHOD_BFBCG] = 1; /* Krylov subspace not enlarged by default */
    170:       } else data->scntl[data->cntl[0] != HPDDM_KRYLOV_METHOD_BCG] = 0;
    171:       if (data->cntl[0] == HPDDM_KRYLOV_METHOD_GMRES || data->cntl[0] == HPDDM_KRYLOV_METHOD_BGMRES || data->cntl[0] == HPDDM_KRYLOV_METHOD_GCRODR || data->cntl[0] == HPDDM_KRYLOV_METHOD_BGCRODR) {
    172:         data->cntl[2] = static_cast<char>(HPDDM_ORTHOGONALIZATION_CGS) + (static_cast<char>(HPDDM_QR_CHOLQR) << 2); /* CGS and CholQR by default */
    173:         data->scntl[0] = PetscMin(30, ksp->max_it); /* restart parameter of 30 by default */
    174:       }
    175:       if (data->cntl[0] == HPDDM_KRYLOV_METHOD_BCG || data->cntl[0] == HPDDM_KRYLOV_METHOD_BFBCG) {
    176:         data->cntl[1] = HPDDM_QR_CHOLQR; /* CholQR by default */
    177:       }
    178:       if (data->cntl[0] == HPDDM_KRYLOV_METHOD_GCRODR || data->cntl[0] == HPDDM_KRYLOV_METHOD_BGCRODR) {
    179:         data->icntl[0] = PetscMin(20, data->scntl[0] - 1); /* recycled subspace of size 20 by default */
    180:         if (!PetscDefined(HAVE_SLEPC) || !PetscDefined(USE_SHARED_LIBRARIES) || data->cntl[0] == HPDDM_KRYLOV_METHOD_GCRODR) {
    181:           data->cntl[3] = HPDDM_RECYCLE_TARGET_SM; /* default recycling target */
    182:         } else {
    183:           data->cntl[3] = 1; /* redistribution parameter of 1 by default */
    184:         }
    185:         data->cntl[4] = HPDDM_RECYCLE_STRATEGY_A; /* default recycling strategy */
    186:       }
    187:     } else data->scntl[1] = 1;
    188:   }
    189:   if (ksp->nmax > std::numeric_limits<int>::max() || ksp->nmax < std::numeric_limits<int>::min()) SETERRQ1(PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "KSPMatSolve() block size %D not representable by an integer, which is not handled by KSPHPDDM", ksp->nmax);
    190:   else data->icntl[1] = static_cast<int>(ksp->nmax);
    191:   return(0);
    192: }

    194: PETSC_STATIC_INLINE PetscErrorCode KSPHPDDMReset_Private(KSP ksp)
    195: {
    196:   KSP_HPDDM *data = (KSP_HPDDM*)ksp->data;

    199:   /* cast PETSC_DECIDE into the appropriate types to avoid compiler warnings */
    200:   std::fill_n(data->rcntl, ALEN(data->rcntl), static_cast<PetscReal>(PETSC_DECIDE));
    201:   std::fill_n(data->icntl, ALEN(data->icntl), static_cast<int>(PETSC_DECIDE));
    202:   std::fill_n(data->scntl, ALEN(data->scntl), static_cast<unsigned short>(PETSC_DECIDE));
    203:   std::fill_n(data->cntl , ALEN(data->cntl) , static_cast<char>(PETSC_DECIDE));
    204:   return(0);
    205: }

    207: static PetscErrorCode KSPReset_HPDDM(KSP ksp)
    208: {
    209:   KSP_HPDDM      *data = (KSP_HPDDM*)ksp->data;

    213:   if (data->op) {
    214:     delete data->op;
    215:     data->op = NULL;
    216:   }
    217:   KSPHPDDMReset_Private(ksp);
    218:   return(0);
    219: }

    221: static PetscErrorCode KSPDestroy_HPDDM(KSP ksp)
    222: {

    226:   KSPReset_HPDDM(ksp);
    227:   KSPDestroyDefault(ksp);
    228:   PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMSetDeflationSpace_C", NULL);
    229:   PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMGetDeflationSpace_C", NULL);
    230:   PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMSetType_C", NULL);
    231:   PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMGetType_C", NULL);
    232:   return(0);
    233: }

    235: PETSC_STATIC_INLINE PetscErrorCode KSPSolve_HPDDM_Private(KSP ksp, const PetscScalar *b, PetscScalar *x, PetscInt n)
    236: {
    237:   KSP_HPDDM              *data = (KSP_HPDDM*)ksp->data;
    238:   KSPConvergedDefaultCtx *ctx = (KSPConvergedDefaultCtx*)ksp->cnvP;
    239:   PetscBool              scale;
    240:   PetscErrorCode         ierr;

    243:   PCGetDiagonalScale(ksp->pc, &scale);
    244:   if (scale) SETERRQ1(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
    245:   if (n > 1) {
    246:     if (ksp->converged == KSPConvergedDefault) {
    247:       if (ctx->mininitialrtol) SETERRQ1(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support KSPConvergedDefaultSetUMIRNorm()", ((PetscObject)ksp)->type_name);
    248:       if (!ctx->initialrtol) {
    249:         PetscInfo(ksp, "Forcing KSPConvergedDefaultSetUIRNorm() since KSPConvergedDefault() cannot handle multiple norms\n");
    250:         ctx->initialrtol = PETSC_TRUE;
    251:       }
    252:     } else {
    253:       PetscInfo(ksp, "Using a special \"converged\" callback, be careful, it is used in KSPHPDDM to track blocks of residuals\n");
    254:     }
    255:   }
    256:   /* initial guess is always nonzero with recycling methods if there is a deflation subspace available */
    257:   if ((data->cntl[0] == HPDDM_KRYLOV_METHOD_GCRODR || data->cntl[0] == HPDDM_KRYLOV_METHOD_BGCRODR) && data->op->storage()) ksp->guess_zero = PETSC_FALSE;
    258:   ksp->its = 0;
    259:   ksp->reason = KSP_CONVERGED_ITERATING;
    260:   static_cast<PetscErrorCode>(HPDDM::IterativeMethod::solve(*data->op, b, x, n, PetscObjectComm((PetscObject)ksp)));
    261:   if (!ksp->reason) { /* KSPConvergedDefault() is still returning 0 (= KSP_CONVERGED_ITERATING) */
    262:     if (ksp->its >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
    263:     else ksp->reason = KSP_CONVERGED_RTOL; /* early exit by HPDDM, which only happens on breakdowns or convergence */
    264:   }
    265:   ksp->its = PetscMin(ksp->its, ksp->max_it);
    266:   return(0);
    267: }

    269: static PetscErrorCode KSPSolve_HPDDM(KSP ksp)
    270: {
    271:   KSP_HPDDM         *data = (KSP_HPDDM*)ksp->data;
    272:   Mat               A, B;
    273:   PetscScalar       *x, *bt = NULL, **ptr;
    274:   const PetscScalar *b;
    275:   PetscInt          i, j, n;
    276:   PetscBool         flg;
    277:   PetscErrorCode    ierr;

    280:   PetscCitationsRegister(HPDDMCitation, &HPDDMCite);
    281:   KSPGetOperators(ksp, &A, NULL);
    282:   PetscObjectTypeCompareAny((PetscObject)A, &flg, MATSEQKAIJ, MATMPIKAIJ, "");
    283:   VecGetArrayWrite(ksp->vec_sol, &x);
    284:   VecGetArrayRead(ksp->vec_rhs, &b);
    285:   if (!flg) {
    286:     KSPSolve_HPDDM_Private(ksp, b, x, 1);
    287:   } else {
    288:     MatKAIJGetScaledIdentity(A, &flg);
    289:     MatKAIJGetAIJ(A, &B);
    290:     MatGetBlockSize(A, &n);
    291:     MatGetLocalSize(B, &i, NULL);
    292:     j = data->op->getDof();
    293:     if (!flg) i *= n; /* S and T are not scaled identities, cannot use block methods */
    294:     if (i != j) { /* switching between block and standard methods */
    295:       delete data->op;
    296:       data->op = new HPDDM::PETScOperator(ksp, i);
    297:     }
    298:     if (flg && n > 1) {
    299:       PetscMalloc1(i * n, &bt);
    300:       /* from row- to column-major to be consistent with HPDDM */
    301:       HPDDM::Wrapper<PetscScalar>::omatcopy<'T'>(i, n, b, n, bt, i);
    302:       ptr = const_cast<PetscScalar**>(&b);
    303:       std::swap(*ptr, bt);
    304:       HPDDM::Wrapper<PetscScalar>::imatcopy<'T'>(i, n, x, n, i);
    305:     }
    306:     KSPSolve_HPDDM_Private(ksp, b, x, flg ? n : 1);
    307:     if (flg && n > 1) {
    308:       std::swap(*ptr, bt);
    309:       PetscFree(bt);
    310:       /* from column- to row-major to be consistent with MatKAIJ format */
    311:       HPDDM::Wrapper<PetscScalar>::imatcopy<'T'>(n, i, x, i, n);
    312:     }
    313:   }
    314:   VecRestoreArrayRead(ksp->vec_rhs, &b);
    315:   VecRestoreArrayWrite(ksp->vec_sol, &x);
    316:   return(0);
    317: }

    319: /*@
    320:      KSPHPDDMSetDeflationSpace - Sets the deflation space used by Krylov methods with recycling. This space is viewed as a set of vectors stored in a MATDENSE (column major).

    322:    Input Parameters:
    323: +     ksp - iterative context
    324: -     U - deflation space to be used during KSPSolve()

    326:    Level: intermediate

    328: .seealso:  KSPCreate(), KSPType (for list of available types), KSPHPDDMGetDeflationSpace()
    329: @*/
    330: PetscErrorCode KSPHPDDMSetDeflationSpace(KSP ksp, Mat U)
    331: {

    338:   PetscUseMethod(ksp, "KSPHPDDMSetDeflationSpace_C", (KSP, Mat), (ksp, U));
    339:   return(0);
    340: }

    342: /*@
    343:      KSPHPDDMGetDeflationSpace - Gets the deflation space computed by Krylov methods with recycling or NULL if KSPSolve() has not been called yet. This space is viewed as a set of vectors stored in a MATDENSE (column major). It is the responsibility of the user to free the returned Mat.

    345:    Input Parameter:
    346: .     ksp - iterative context

    348:    Output Parameter:
    349: .     U - deflation space generated during KSPSolve()

    351:    Level: intermediate

    353: .seealso:  KSPCreate(), KSPType (for list of available types), KSPHPDDMSetDeflationSpace()
    354: @*/
    355: PetscErrorCode KSPHPDDMGetDeflationSpace(KSP ksp, Mat *U)
    356: {

    361:   PetscUseMethod(ksp, "KSPHPDDMGetDeflationSpace_C", (KSP, Mat*), (ksp, U));
    362:   return(0);
    363: }

    365: static PetscErrorCode KSPHPDDMSetDeflationSpace_HPDDM(KSP ksp, Mat U)
    366: {
    367:   KSP_HPDDM            *data = (KSP_HPDDM*)ksp->data;
    368:   HPDDM::PETScOperator *op = data->op;
    369:   Mat                  A;
    370:   const PetscScalar    *array;
    371:   PetscScalar          *copy;
    372:   PetscInt             m1, M1, m2, M2, n2, N2, ldu;
    373:   PetscBool            match;
    374:   PetscErrorCode       ierr;

    377:   if (!op) {
    378:     KSPSetUp(ksp);
    379:     op = data->op;
    380:   }
    381:   KSPGetOperators(ksp, &A, NULL);
    382:   MatGetLocalSize(A, &m1, NULL);
    383:   MatGetLocalSize(U, &m2, &n2);
    384:   MatGetSize(A, &M1, NULL);
    385:   MatGetSize(U, &M2, &N2);
    386:   if (m1 != m2 || M1 != M2) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Cannot use a deflation space with (m2,M2) = (%D,%D) for a linear system with (m1,M1) = (%D,%D)", m2, M2, m1, M1);
    387:   PetscObjectTypeCompareAny((PetscObject)U, &match, MATSEQDENSE, MATMPIDENSE, "");
    388:   if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Provided deflation space not stored in a dense Mat");
    389:   MatDenseGetArrayRead(U, &array);
    390:   copy = op->allocate(m2, 1, N2);
    391:   if (!copy) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_POINTER, "Memory allocation error");
    392:   MatDenseGetLDA(U, &ldu);
    393:   HPDDM::Wrapper<PetscScalar>::omatcopy<'N'>(N2, m2, array, ldu, copy, m2);
    394:   MatDenseRestoreArrayRead(U, &array);
    395:   return(0);
    396: }

    398: static PetscErrorCode KSPHPDDMGetDeflationSpace_HPDDM(KSP ksp, Mat *U)
    399: {
    400:   KSP_HPDDM            *data = (KSP_HPDDM*)ksp->data;
    401:   HPDDM::PETScOperator *op = data->op;
    402:   Mat                  A;
    403:   const PetscScalar    *array;
    404:   PetscScalar          *copy;
    405:   PetscInt             m1, M1, N2;
    406:   PetscErrorCode       ierr;

    409:   if (!op) {
    410:     KSPSetUp(ksp);
    411:     op = data->op;
    412:   }
    413:   array = op->storage();
    414:   N2 = op->k().first * op->k().second;
    415:   if (!array) *U = NULL;
    416:   else {
    417:     KSPGetOperators(ksp, &A, NULL);
    418:     MatGetLocalSize(A, &m1, NULL);
    419:     MatGetSize(A, &M1, NULL);
    420:     MatCreateDense(PetscObjectComm((PetscObject)ksp), m1, PETSC_DECIDE, M1, N2, NULL, U);
    421:     MatDenseGetArrayWrite(*U, &copy);
    422:     PetscArraycpy(copy, array, m1 * N2);
    423:     MatDenseRestoreArrayWrite(*U, &copy);
    424:   }
    425:   return(0);
    426: }

    428: static PetscErrorCode KSPMatSolve_HPDDM(KSP ksp, Mat B, Mat X)
    429: {
    430:   KSP_HPDDM            *data = (KSP_HPDDM*)ksp->data;
    431:   HPDDM::PETScOperator *op = data->op;
    432:   Mat                  A;
    433:   const PetscScalar    *b;
    434:   PetscScalar          *x;
    435:   PetscInt             n, lda;
    436:   PetscErrorCode       ierr;

    439:   if (!op) {
    440:     KSPSetUp(ksp);
    441:     op = data->op;
    442:   }
    443:   KSPGetOperators(ksp, &A, NULL);
    444:   MatGetLocalSize(B, &n, NULL);
    445:   MatDenseGetLDA(B, &lda);
    446:   if (n != lda) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Unhandled leading dimension lda = %D with n = %D", lda, n);
    447:   MatGetLocalSize(A, &n, NULL);
    448:   MatDenseGetLDA(X, &lda);
    449:   if (n != lda) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Unhandled leading dimension lda = %D with n = %D", lda, n);
    450:   MatDenseGetArrayRead(B, &b);
    451:   MatDenseGetArrayWrite(X, &x);
    452:   MatGetSize(X, NULL, &n);
    453:   KSPSolve_HPDDM_Private(ksp, b, x, n);
    454:   MatDenseRestoreArrayWrite(X, &x);
    455:   MatDenseRestoreArrayRead(B, &b);
    456:   return(0);
    457: }

    459: /*@
    460:      KSPHPDDMSetType - Sets the type of Krylov method used in KSPHPDDM.

    462:    Input Parameters:
    463: +     ksp - iterative context
    464: -     type - any of gmres, bgmres, cg, bcg, gcrodr, bgcrodr, bfbcg, or preonly

    466:    Level: intermediate

    468:    Notes:
    469:      Unlike KSPReset(), this function does not destroy any deflation space attached to the KSP.
    470:      As an example, in the following sequence: KSPHPDDMSetType(ksp, KSPGCRODR); KSPSolve(ksp, b, x); KSPHPDDMSetType(ksp, KSPGMRES); KSPHPDDMSetType(ksp, KSPGCRODR); KSPSolve(ksp, b, x); the recycled space is reused in the second KSPSolve().

    472: .seealso:  KSPCreate(), KSPType (for list of available types), KSPHPDDMType, KSPHPDDMGetType()
    473: @*/
    474: PetscErrorCode KSPHPDDMSetType(KSP ksp, KSPHPDDMType type)
    475: {

    480:   PetscUseMethod(ksp, "KSPHPDDMSetType_C", (KSP, KSPHPDDMType), (ksp, type));
    481:   return(0);
    482: }

    484: /*@
    485:      KSPHPDDMGetType - Gets the type of Krylov method used in KSPHPDDM.

    487:    Input Parameter:
    488: .     ksp - iterative context

    490:    Output Parameter:
    491: .     type - any of gmres, bgmres, cg, bcg, gcrodr, bgcrodr, bfbcg, or preonly

    493:    Level: intermediate

    495: .seealso:  KSPCreate(), KSPType (for list of available types), KSPHPDDMType, KSPHPDDMSetType()
    496: @*/
    497: PetscErrorCode KSPHPDDMGetType(KSP ksp, KSPHPDDMType *type)
    498: {

    504:   PetscUseMethod(ksp, "KSPHPDDMGetType_C", (KSP, KSPHPDDMType*), (ksp, type));
    505:   return(0);
    506: }

    508: static PetscErrorCode KSPHPDDMSetType_HPDDM(KSP ksp, KSPHPDDMType type)
    509: {
    510:   KSP_HPDDM      *data = (KSP_HPDDM*)ksp->data;
    511:   PetscInt       i;
    512:   PetscBool      flg = PETSC_FALSE;

    516:   for (i = 0; i < static_cast<PetscInt>(ALEN(KSPHPDDMTypes)); ++i) {
    517:     PetscStrcmp(KSPHPDDMTypes[type], KSPHPDDMTypes[i], &flg);
    518:     if (flg) break;
    519:   }
    520:   if (i == ALEN(KSPHPDDMTypes)) SETERRQ1(PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown KSPHPDDMType %s", type);
    521:   if (data->cntl[0] != static_cast<char>(PETSC_DECIDE) && data->cntl[0] != i) {
    522:     KSPHPDDMReset_Private(ksp);
    523:   }
    524:   data->cntl[0] = i;
    525:   return(0);
    526: }

    528: static PetscErrorCode KSPHPDDMGetType_HPDDM(KSP ksp, KSPHPDDMType *type)
    529: {
    530:   KSP_HPDDM *data = (KSP_HPDDM*)ksp->data;

    533:   if (data->cntl[0] == static_cast<char>(PETSC_DECIDE)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ORDER, "KSPHPDDMType not set yet");
    534:   /* need to shift by -1 for HPDDM_KRYLOV_METHOD_NONE */
    535:   *type = static_cast<KSPHPDDMType>(PetscMin(data->cntl[0], static_cast<char>(ALEN(KSPHPDDMTypes) - 1)));
    536:   return(0);
    537: }

    539: /*MC
    540:      KSPHPDDM - Interface with the HPDDM library.

    542:    This KSP may be used to further select methods that are currently not implemented natively in PETSc, e.g., GCRODR [2006], a recycled Krylov method which is similar to KSPLGMRES, see [2016] for a comparison. ex75.c shows how to reproduce the results from the aforementioned paper [2006]. A chronological bibliography of relevant publications linked with KSP available in HPDDM through KSPHPDDM, and not available directly in PETSc, may be found below. The interface is explained in details in [2021].

    544:    Options Database Keys:
    545: +   -ksp_gmres_restart <restart, default=30> - see KSPGMRES
    546: .   -ksp_hpddm_type <type, default=gmres> - any of gmres, bgmres, cg, bcg, gcrodr, bgcrodr, bfbcg, or preonly, see KSPHPDDMType
    547: .   -ksp_hpddm_deflation_tol <eps, default=\-1.0> - tolerance when deflating right-hand sides inside block methods (no deflation by default, only relevant with block methods)
    548: .   -ksp_hpddm_enlarge_krylov_subspace <p, default=1> - split the initial right-hand side into multiple vectors (only relevant with nonblock methods)
    549: .   -ksp_hpddm_orthogonalization <type, default=cgs> - any of cgs or mgs, see KSPGMRES
    550: .   -ksp_hpddm_qr <type, default=cholqr> - distributed QR factorizations with any of cholqr, cgs, or mgs (only relevant with block methods)
    551: .   -ksp_hpddm_variant <type, default=left> - any of left, right, or flexible (this option is superseded by KSPSetPCSide())
    552: .   -ksp_hpddm_recycle <n, default=0> - number of harmonic Ritz vectors to compute (only relevant with GCRODR or BGCRODR)
    553: .   -ksp_hpddm_recycle_target <type, default=SM> - criterion to select harmonic Ritz vectors using either SM, LM, SR, LR, SI, or LI (only relevant with GCRODR or BGCRODR). For BGCRODR, if PETSc is compiled with SLEPc, this option is not relevant, since SLEPc is used instead. Options are set with the prefix -ksp_hpddm_recycle_eps_
    554: .   -ksp_hpddm_recycle_strategy <type, default=A> - generalized eigenvalue problem A or B to solve for recycling (only relevant with flexible GCRODR or BGCRODR)
    555: -   -ksp_hpddm_recycle_symmetric <true, default=false> - symmetric generalized eigenproblems in BGCRODR, useful to switch to distributed solvers like EPSELEMENTAL or EPSSCALAPACK (only relevant when PETSc is compiled with SLEPc)

    557:    References:
    558: +   1980 - The Block Conjugate Gradient Algorithm and Related Methods. O'Leary. Linear Algebra and its Applications.
    559: .   2006 - Recycling Krylov Subspaces for Sequences of Linear Systems. Parks, de Sturler, Mackey, Johnson, and Maiti. SIAM Journal on Scientific Computing
    560: .   2013 - A Modified Block Flexible GMRES Method with Deflation at Each Iteration for the Solution of Non-Hermitian Linear Systems with Multiple Right-Hand Sides. Calandra, Gratton, Lago, Vasseur, and Carvalho. SIAM Journal on Scientific Computing.
    561: .   2016 - Block Iterative Methods and Recycling for Improved Scalability of Linear Solvers. Jolivet and Tournier. SC16.
    562: .   2017 - A breakdown-free block conjugate gradient method. Ji and Li. BIT Numerical Mathematics.
    563: -   2021 - KSPHPDDM and PCHPDDM: extending PETSc with advanced Krylov methods and robust multilevel overlapping Schwarz preconditioners. Jolivet, Roman, and Zampini. Computer & Mathematics with Applications.

    565:    Level: intermediate

    567: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPGMRES, KSPCG, KSPLGMRES, KSPDGMRES
    568: M*/
    569: PETSC_EXTERN PetscErrorCode KSPCreate_HPDDM(KSP ksp)
    570: {
    571:   KSP_HPDDM      *data;
    572:   PetscInt       i;
    573:   const char     *common[] = { KSPGMRES, KSPCG, KSPPREONLY };
    574:   PetscBool      flg = PETSC_FALSE;

    578:   PetscNewLog(ksp, &data);
    579:   ksp->data = (void*)data;
    580:   KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 2);
    581:   KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 1);
    582:   ksp->ops->solve          = KSPSolve_HPDDM;
    583:   ksp->ops->matsolve       = KSPMatSolve_HPDDM;
    584:   ksp->ops->setup          = KSPSetUp_HPDDM;
    585:   ksp->ops->setfromoptions = KSPSetFromOptions_HPDDM;
    586:   ksp->ops->destroy        = KSPDestroy_HPDDM;
    587:   ksp->ops->view           = KSPView_HPDDM;
    588:   ksp->ops->reset          = KSPReset_HPDDM;
    589:   KSPHPDDMReset_Private(ksp);
    590:   for (i = 0; i < static_cast<PetscInt>(ALEN(common)); ++i) {
    591:     PetscStrcmp(((PetscObject)ksp)->type_name, common[i], &flg);
    592:     if (flg) break;
    593:   }
    594:   if (!i) data->cntl[0] = HPDDM_KRYLOV_METHOD_GMRES;
    595:   else if (i == 1) data->cntl[0] = HPDDM_KRYLOV_METHOD_CG;
    596:   else if (i == 2) data->cntl[0] = HPDDM_KRYLOV_METHOD_NONE;
    597:   if (data->cntl[0] != static_cast<char>(PETSC_DECIDE)) {
    598:     PetscInfo1(ksp, "Using the previously set KSPType %s\n", common[i]);
    599:   }
    600:   PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMSetDeflationSpace_C", KSPHPDDMSetDeflationSpace_HPDDM);
    601:   PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMGetDeflationSpace_C", KSPHPDDMGetDeflationSpace_HPDDM);
    602:   PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMSetType_C", KSPHPDDMSetType_HPDDM);
    603:   PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMGetType_C", KSPHPDDMGetType_HPDDM);
    604: #if defined(PETSC_HAVE_SLEPC) && defined(PETSC_USE_SHARED_LIBRARIES)
    605:   if (!loadedDL) {
    606:     HPDDMLoadDL_Private(&loadedDL);
    607:   }
    608: #endif
    609:   return(0);
    610: }