1 #include <petsc/private/vecimpl.h> 2 #include <petsc/private/matimpl.h> 3 #include <petsc/private/petschpddm.h> /*I "petscpc.h" I*/ 4 #include <petsc/private/pcimpl.h> 5 #include <petsc/private/dmimpl.h> /* this must be included after petschpddm.h so that DM_MAX_WORK_VECTORS is not defined */ 6 /* otherwise, it is assumed that one is compiling libhpddm_petsc => circular dependency */ 7 #if PetscDefined(USE_FORTRAN_BINDINGS) 8 #include <petsc/private/fortranimpl.h> 9 #endif 10 11 static PetscErrorCode (*loadedSym)(HPDDM::Schwarz<PetscScalar> *const, IS, Mat, Mat, Mat, std::vector<Vec>, PC_HPDDM_Level **const) = nullptr; 12 13 static PetscBool PCHPDDMPackageInitialized = PETSC_FALSE; 14 15 PetscLogEvent PC_HPDDM_Strc; 16 PetscLogEvent PC_HPDDM_PtAP; 17 PetscLogEvent PC_HPDDM_PtBP; 18 PetscLogEvent PC_HPDDM_Next; 19 PetscLogEvent PC_HPDDM_SetUp[PETSC_PCHPDDM_MAXLEVELS]; 20 PetscLogEvent PC_HPDDM_Solve[PETSC_PCHPDDM_MAXLEVELS]; 21 22 const char *const PCHPDDMCoarseCorrectionTypes[] = {"DEFLATED", "ADDITIVE", "BALANCED", "NONE", "PCHPDDMCoarseCorrectionType", "PC_HPDDM_COARSE_CORRECTION_", nullptr}; 23 const char *const PCHPDDMSchurPreTypes[] = {"LEAST_SQUARES", "GENEO", "PCHPDDMSchurPreType", "PC_HPDDM_SCHUR_PRE", nullptr}; 24 25 static PetscErrorCode PCReset_HPDDM(PC pc) 26 { 27 PC_HPDDM *data = (PC_HPDDM *)pc->data; 28 PetscInt i; 29 30 PetscFunctionBegin; 31 if (data->levels) { 32 for (i = 0; i < PETSC_PCHPDDM_MAXLEVELS && data->levels[i]; ++i) { 33 PetscCall(KSPDestroy(&data->levels[i]->ksp)); 34 PetscCall(PCDestroy(&data->levels[i]->pc)); 35 PetscCall(PetscFree(data->levels[i])); 36 } 37 PetscCall(PetscFree(data->levels)); 38 } 39 40 PetscCall(ISDestroy(&data->is)); 41 PetscCall(MatDestroy(&data->aux)); 42 PetscCall(MatDestroy(&data->B)); 43 PetscCall(VecDestroy(&data->normal)); 44 data->correction = PC_HPDDM_COARSE_CORRECTION_DEFLATED; 45 data->Neumann = PETSC_BOOL3_UNKNOWN; 46 data->deflation = PETSC_FALSE; 47 data->setup = nullptr; 48 data->setup_ctx = nullptr; 49 PetscFunctionReturn(PETSC_SUCCESS); 50 } 51 52 static PetscErrorCode PCDestroy_HPDDM(PC pc) 53 { 54 PC_HPDDM *data = (PC_HPDDM *)pc->data; 55 56 PetscFunctionBegin; 57 PetscCall(PCReset_HPDDM(pc)); 58 PetscCall(PetscFree(data)); 59 PetscCall(PetscObjectChangeTypeName((PetscObject)pc, nullptr)); 60 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetAuxiliaryMat_C", nullptr)); 61 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMHasNeumannMat_C", nullptr)); 62 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetRHSMat_C", nullptr)); 63 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetCoarseCorrectionType_C", nullptr)); 64 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMGetCoarseCorrectionType_C", nullptr)); 65 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetSTShareSubKSP_C", nullptr)); 66 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMGetSTShareSubKSP_C", nullptr)); 67 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetDeflationMat_C", nullptr)); 68 PetscCall(PetscObjectCompose((PetscObject)pc, "_PCHPDDM_Schur", nullptr)); 69 PetscFunctionReturn(PETSC_SUCCESS); 70 } 71 72 static inline PetscErrorCode PCHPDDMSetAuxiliaryMat_Private(PC pc, IS is, Mat A, PetscBool deflation) 73 { 74 PC_HPDDM *data = (PC_HPDDM *)pc->data; 75 76 PetscFunctionBegin; 77 if (is) { 78 PetscCall(PetscObjectReference((PetscObject)is)); 79 if (data->is) { /* new overlap definition resets the PC */ 80 PetscCall(PCReset_HPDDM(pc)); 81 pc->setfromoptionscalled = 0; 82 pc->setupcalled = 0; 83 } 84 PetscCall(ISDestroy(&data->is)); 85 data->is = is; 86 } 87 if (A) { 88 PetscCall(PetscObjectReference((PetscObject)A)); 89 PetscCall(MatDestroy(&data->aux)); 90 data->aux = A; 91 } 92 data->deflation = deflation; 93 PetscFunctionReturn(PETSC_SUCCESS); 94 } 95 96 static inline PetscErrorCode PCHPDDMSetAuxiliaryMatNormal_Private(PC pc, Mat A, Mat N, Mat *B, const char *pcpre, Vec *diagonal = nullptr) 97 { 98 PC_HPDDM *data = (PC_HPDDM *)pc->data; 99 Mat *splitting, *sub, aux; 100 Vec d; 101 IS is, cols[2], rows; 102 PetscReal norm; 103 PetscBool flg; 104 char type[256] = {}; /* same size as in src/ksp/pc/interface/pcset.c */ 105 106 PetscFunctionBegin; 107 PetscCall(MatConvert(N, MATAIJ, MAT_INITIAL_MATRIX, B)); 108 PetscCall(ISCreateStride(PETSC_COMM_SELF, A->cmap->n, A->cmap->rstart, 1, cols)); 109 PetscCall(ISCreateStride(PETSC_COMM_SELF, A->rmap->N, 0, 1, &rows)); 110 PetscCall(MatSetOption(A, MAT_SUBMAT_SINGLEIS, PETSC_TRUE)); 111 PetscCall(MatIncreaseOverlap(*B, 1, cols, 1)); 112 PetscCall(MatCreateSubMatrices(A, 1, &rows, cols, MAT_INITIAL_MATRIX, &splitting)); 113 PetscCall(ISCreateStride(PETSC_COMM_SELF, A->cmap->n, A->cmap->rstart, 1, &is)); 114 PetscCall(ISEmbed(*cols, is, PETSC_TRUE, cols + 1)); 115 PetscCall(ISDestroy(&is)); 116 PetscCall(MatCreateSubMatrices(*splitting, 1, &rows, cols + 1, MAT_INITIAL_MATRIX, &sub)); 117 PetscCall(ISDestroy(cols + 1)); 118 PetscCall(MatFindZeroRows(*sub, &is)); 119 PetscCall(MatDestroySubMatrices(1, &sub)); 120 PetscCall(ISDestroy(&rows)); 121 PetscCall(MatSetOption(*splitting, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE)); 122 PetscCall(MatZeroRowsIS(*splitting, is, 0.0, nullptr, nullptr)); 123 PetscCall(ISDestroy(&is)); 124 PetscCall(PetscOptionsGetString(nullptr, pcpre, "-pc_hpddm_levels_1_sub_pc_type", type, sizeof(type), nullptr)); 125 PetscCall(PetscStrcmp(type, PCQR, &flg)); 126 if (!flg) { 127 Mat conjugate = *splitting; 128 if (PetscDefined(USE_COMPLEX)) { 129 PetscCall(MatDuplicate(*splitting, MAT_COPY_VALUES, &conjugate)); 130 PetscCall(MatConjugate(conjugate)); 131 } 132 PetscCall(MatTransposeMatMult(conjugate, *splitting, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &aux)); 133 if (PetscDefined(USE_COMPLEX)) PetscCall(MatDestroy(&conjugate)); 134 PetscCall(MatNorm(aux, NORM_FROBENIUS, &norm)); 135 PetscCall(MatSetOption(aux, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE)); 136 if (diagonal) { 137 PetscCall(VecNorm(*diagonal, NORM_INFINITY, &norm)); 138 if (norm > PETSC_SMALL) { 139 VecScatter scatter; 140 PetscInt n; 141 PetscCall(ISGetLocalSize(*cols, &n)); 142 PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)pc), n, PETSC_DECIDE, &d)); 143 PetscCall(VecScatterCreate(*diagonal, *cols, d, nullptr, &scatter)); 144 PetscCall(VecScatterBegin(scatter, *diagonal, d, INSERT_VALUES, SCATTER_FORWARD)); 145 PetscCall(VecScatterEnd(scatter, *diagonal, d, INSERT_VALUES, SCATTER_FORWARD)); 146 PetscCall(VecScatterDestroy(&scatter)); 147 PetscCall(MatScale(aux, -1.0)); 148 PetscCall(MatDiagonalSet(aux, d, ADD_VALUES)); 149 PetscCall(VecDestroy(&d)); 150 } else PetscCall(VecDestroy(diagonal)); 151 } 152 if (!diagonal) PetscCall(MatShift(aux, PETSC_SMALL * norm)); 153 } else { 154 PetscBool flg; 155 if (diagonal) { 156 PetscCall(VecNorm(*diagonal, NORM_INFINITY, &norm)); 157 PetscCheck(norm < PETSC_SMALL, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Nonzero diagonal A11 block"); 158 PetscCall(VecDestroy(diagonal)); 159 } 160 PetscCall(PetscObjectTypeCompare((PetscObject)N, MATNORMAL, &flg)); 161 if (flg) PetscCall(MatCreateNormal(*splitting, &aux)); 162 else PetscCall(MatCreateNormalHermitian(*splitting, &aux)); 163 } 164 PetscCall(MatDestroySubMatrices(1, &splitting)); 165 PetscCall(PCHPDDMSetAuxiliaryMat(pc, *cols, aux, nullptr, nullptr)); 166 data->Neumann = PETSC_BOOL3_TRUE; 167 PetscCall(ISDestroy(cols)); 168 PetscCall(MatDestroy(&aux)); 169 PetscFunctionReturn(PETSC_SUCCESS); 170 } 171 172 static PetscErrorCode PCHPDDMSetAuxiliaryMat_HPDDM(PC pc, IS is, Mat A, PetscErrorCode (*setup)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void *), void *setup_ctx) 173 { 174 PC_HPDDM *data = (PC_HPDDM *)pc->data; 175 176 PetscFunctionBegin; 177 PetscCall(PCHPDDMSetAuxiliaryMat_Private(pc, is, A, PETSC_FALSE)); 178 if (setup) { 179 data->setup = setup; 180 data->setup_ctx = setup_ctx; 181 } 182 PetscFunctionReturn(PETSC_SUCCESS); 183 } 184 185 typedef struct { 186 KSP ksp; 187 PetscInt its; 188 } PC_KSP; 189 190 static inline PetscErrorCode PCSetUp_KSP_Private(PC pc) 191 { 192 PC_KSP *data = (PC_KSP *)pc->data; 193 const std::string prefix = ((PetscObject)data->ksp)->prefix; 194 std::string sub; 195 196 PetscFunctionBegin; 197 PetscCheck(prefix.size() >= 9, PETSC_COMM_SELF, PETSC_ERR_PLIB, "The prefix of this PCKSP should be of length at least 9 to hold the suffix pc_hpddm_"); 198 sub = prefix.substr(0, prefix.size() - 9); 199 PetscCall(PCSetType(pc, PCHPDDM)); 200 PetscCall(PCSetOptionsPrefix(pc, sub.c_str())); 201 PetscFunctionReturn(PETSC_SUCCESS); 202 } 203 204 static PetscErrorCode PCSetUp_KSP(PC pc) 205 { 206 PetscFunctionBegin; 207 PetscCall(PCSetUp_KSP_Private(pc)); 208 PetscCall(PCSetFromOptions(pc)); 209 PetscCall(PCSetUp(pc)); 210 PetscFunctionReturn(PETSC_SUCCESS); 211 } 212 213 /*@C 214 PCHPDDMSetAuxiliaryMat - Sets the auxiliary matrix used by `PCHPDDM` for the concurrent GenEO problems at the finest level. 215 216 Input Parameters: 217 + pc - preconditioner context 218 . is - index set of the local auxiliary, e.g., Neumann, matrix 219 . A - auxiliary sequential matrix 220 . setup - function for generating the auxiliary matrix entries, may be `NULL` 221 - ctx - context for `setup`, may be `NULL` 222 223 Calling sequence of `setup`: 224 + J - matrix whose values are to be set 225 . t - time 226 . X - linearization point 227 . X_t - time-derivative of the linearization point 228 . s - step 229 . ovl - index set of the local auxiliary, e.g., Neumann, matrix 230 - ctx - context for `setup`, may be `NULL` 231 232 Level: intermediate 233 234 Note: 235 As an example, in a finite element context with nonoverlapping subdomains plus (overlapping) ghost elements, this could be the unassembled (Neumann) 236 local overlapping operator. As opposed to the assembled (Dirichlet) local overlapping operator obtained by summing neighborhood contributions 237 at the interface of ghost elements. 238 239 Fortran Notes: 240 Only `PETSC_NULL_FUNCTION` is supported for `setup` and `ctx` is never accessed 241 242 .seealso: [](ch_ksp), `PCHPDDM`, `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCHPDDMSetRHSMat()`, `MATIS` 243 @*/ 244 PetscErrorCode PCHPDDMSetAuxiliaryMat(PC pc, IS is, Mat A, PetscErrorCode (*setup)(Mat J, PetscReal t, Vec X, Vec X_t, PetscReal s, IS ovl, void *ctx), void *ctx) 245 { 246 PetscFunctionBegin; 247 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 248 if (is) PetscValidHeaderSpecific(is, IS_CLASSID, 2); 249 if (A) PetscValidHeaderSpecific(A, MAT_CLASSID, 3); 250 if (pc->ops->setup == PCSetUp_KSP) PetscCall(PCSetUp_KSP_Private(pc)); 251 PetscTryMethod(pc, "PCHPDDMSetAuxiliaryMat_C", (PC, IS, Mat, PetscErrorCode(*)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void *), void *), (pc, is, A, setup, ctx)); 252 PetscFunctionReturn(PETSC_SUCCESS); 253 } 254 255 static PetscErrorCode PCHPDDMHasNeumannMat_HPDDM(PC pc, PetscBool has) 256 { 257 PC_HPDDM *data = (PC_HPDDM *)pc->data; 258 259 PetscFunctionBegin; 260 data->Neumann = PetscBoolToBool3(has); 261 PetscFunctionReturn(PETSC_SUCCESS); 262 } 263 264 /*@ 265 PCHPDDMHasNeumannMat - Informs `PCHPDDM` that the `Mat` passed to `PCHPDDMSetAuxiliaryMat()` is the local Neumann matrix. 266 267 Input Parameters: 268 + pc - preconditioner context 269 - has - Boolean value 270 271 Level: intermediate 272 273 Notes: 274 This may be used to bypass a call to `MatCreateSubMatrices()` and to `MatConvert()` for `MATSBAIJ` matrices. 275 276 If a function is composed with DMCreateNeumannOverlap_C implementation is available in the `DM` attached to the Pmat, or the Amat, or the `PC`, the flag is internally set to `PETSC_TRUE`. Its default value is otherwise `PETSC_FALSE`. 277 278 .seealso: [](ch_ksp), `PCHPDDM`, `PCHPDDMSetAuxiliaryMat()` 279 @*/ 280 PetscErrorCode PCHPDDMHasNeumannMat(PC pc, PetscBool has) 281 { 282 PetscFunctionBegin; 283 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 284 PetscTryMethod(pc, "PCHPDDMHasNeumannMat_C", (PC, PetscBool), (pc, has)); 285 PetscFunctionReturn(PETSC_SUCCESS); 286 } 287 288 static PetscErrorCode PCHPDDMSetRHSMat_HPDDM(PC pc, Mat B) 289 { 290 PC_HPDDM *data = (PC_HPDDM *)pc->data; 291 292 PetscFunctionBegin; 293 PetscCall(PetscObjectReference((PetscObject)B)); 294 PetscCall(MatDestroy(&data->B)); 295 data->B = B; 296 PetscFunctionReturn(PETSC_SUCCESS); 297 } 298 299 /*@ 300 PCHPDDMSetRHSMat - Sets the right-hand side matrix used by `PCHPDDM` for the concurrent GenEO problems at the finest level. 301 302 Input Parameters: 303 + pc - preconditioner context 304 - B - right-hand side sequential matrix 305 306 Level: advanced 307 308 Note: 309 Must be used in conjunction with `PCHPDDMSetAuxiliaryMat`(N), so that Nv = lambda Bv is solved using `EPSSetOperators`(N, B). 310 It is assumed that N and `B` are provided using the same numbering. This provides a means to try more advanced methods such as GenEO-II or H-GenEO. 311 312 .seealso: [](ch_ksp), `PCHPDDMSetAuxiliaryMat()`, `PCHPDDM` 313 @*/ 314 PetscErrorCode PCHPDDMSetRHSMat(PC pc, Mat B) 315 { 316 PetscFunctionBegin; 317 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 318 if (B) { 319 PetscValidHeaderSpecific(B, MAT_CLASSID, 2); 320 PetscTryMethod(pc, "PCHPDDMSetRHSMat_C", (PC, Mat), (pc, B)); 321 } 322 PetscFunctionReturn(PETSC_SUCCESS); 323 } 324 325 static PetscErrorCode PCSetFromOptions_HPDDM(PC pc, PetscOptionItems *PetscOptionsObject) 326 { 327 PC_HPDDM *data = (PC_HPDDM *)pc->data; 328 PC_HPDDM_Level **levels = data->levels; 329 char prefix[256]; 330 int i = 1; 331 PetscMPIInt size, previous; 332 PetscInt n, overlap = 1; 333 PCHPDDMCoarseCorrectionType type; 334 PetscBool flg = PETSC_TRUE, set; 335 336 PetscFunctionBegin; 337 if (!data->levels) { 338 PetscCall(PetscCalloc1(PETSC_PCHPDDM_MAXLEVELS, &levels)); 339 data->levels = levels; 340 } 341 PetscOptionsHeadBegin(PetscOptionsObject, "PCHPDDM options"); 342 PetscCall(PetscOptionsBoundedInt("-pc_hpddm_harmonic_overlap", "Overlap prior to computing local harmonic extensions", "PCHPDDM", overlap, &overlap, &set, 1)); 343 if (!set) overlap = -1; 344 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size)); 345 previous = size; 346 while (i < PETSC_PCHPDDM_MAXLEVELS) { 347 PetscInt p = 1; 348 349 if (!data->levels[i - 1]) PetscCall(PetscNew(data->levels + i - 1)); 350 data->levels[i - 1]->parent = data; 351 /* if the previous level has a single process, it is not possible to coarsen further */ 352 if (previous == 1 || !flg) break; 353 data->levels[i - 1]->nu = 0; 354 data->levels[i - 1]->threshold = -1.0; 355 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_%d_eps_nev", i)); 356 PetscCall(PetscOptionsBoundedInt(prefix, "Local number of deflation vectors computed by SLEPc", "EPSSetDimensions", data->levels[i - 1]->nu, &data->levels[i - 1]->nu, nullptr, 0)); 357 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_%d_eps_threshold", i)); 358 PetscCall(PetscOptionsReal(prefix, "Local threshold for selecting deflation vectors returned by SLEPc", "PCHPDDM", data->levels[i - 1]->threshold, &data->levels[i - 1]->threshold, nullptr)); 359 if (i == 1) { 360 PetscCheck(overlap == -1 || PetscAbsReal(data->levels[i - 1]->threshold + static_cast<PetscReal>(1.0)) < PETSC_MACHINE_EPSILON, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot supply both -pc_hpddm_levels_1_eps_threshold and -pc_hpddm_harmonic_overlap"); 361 if (overlap != -1) { 362 PetscInt nsv = 0; 363 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_%d_svd_nsv", i)); 364 PetscCall(PetscOptionsBoundedInt(prefix, "Local number of deflation vectors computed by SLEPc", "SVDSetDimensions", nsv, &nsv, nullptr, 0)); 365 PetscCheck(bool(data->levels[0]->nu) != bool(nsv), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot supply %s -pc_hpddm_levels_1_eps_nev %s -pc_hpddm_levels_1_svd_nsv", nsv ? "both" : "neither", nsv ? "and" : "nor"); 366 if (nsv) { 367 data->levels[0]->nu = nsv; 368 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_%d_svd_relative_threshold", i)); 369 } else PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_%d_eps_relative_threshold", i)); 370 PetscCall(PetscOptionsReal(prefix, "Local relative threshold for selecting deflation vectors returned by SLEPc", "PCHPDDM", data->levels[0]->threshold, &data->levels[0]->threshold, nullptr)); 371 } 372 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_1_st_share_sub_ksp")); 373 PetscCall(PetscOptionsBool(prefix, "Shared KSP between SLEPc ST and the fine-level subdomain solver", "PCHPDDMSetSTShareSubKSP", PETSC_FALSE, &data->share, nullptr)); 374 } 375 /* if there is no prescribed coarsening, just break out of the loop */ 376 if (data->levels[i - 1]->threshold <= PetscReal() && data->levels[i - 1]->nu <= 0 && !(data->deflation && i == 1)) break; 377 else { 378 ++i; 379 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_%d_eps_nev", i)); 380 PetscCall(PetscOptionsHasName(PetscOptionsObject->options, PetscOptionsObject->prefix, prefix, &flg)); 381 if (!flg) { 382 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_%d_eps_threshold", i)); 383 PetscCall(PetscOptionsHasName(PetscOptionsObject->options, PetscOptionsObject->prefix, prefix, &flg)); 384 } 385 if (flg) { 386 /* if there are coarsening options for the next level, then register it */ 387 /* otherwise, don't to avoid having both options levels_N_p and coarse_p */ 388 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_%d_p", i)); 389 PetscCall(PetscOptionsRangeInt(prefix, "Number of processes used to assemble the coarse operator at this level", "PCHPDDM", p, &p, &flg, 1, PetscMax(1, previous / 2))); 390 previous = p; 391 } 392 } 393 } 394 data->N = i; 395 n = 1; 396 if (i > 1) { 397 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_coarse_p")); 398 PetscCall(PetscOptionsRangeInt(prefix, "Number of processes used to assemble the coarsest operator", "PCHPDDM", n, &n, nullptr, 1, PetscMax(1, previous / 2))); 399 #if PetscDefined(HAVE_MUMPS) 400 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "pc_hpddm_coarse_")); 401 PetscCall(PetscOptionsHasName(nullptr, prefix, "-mat_mumps_use_omp_threads", &flg)); 402 if (flg) { 403 char type[64]; /* same size as in src/ksp/pc/impls/factor/factimpl.c */ 404 PetscCall(PetscStrncpy(type, n > 1 && PetscDefined(HAVE_MUMPS) ? MATSOLVERMUMPS : MATSOLVERPETSC, sizeof(type))); /* default solver for a MatMPIAIJ or a MatSeqAIJ */ 405 PetscCall(PetscOptionsGetString(nullptr, prefix, "-pc_factor_mat_solver_type", type, sizeof(type), nullptr)); 406 PetscCall(PetscStrcmp(type, MATSOLVERMUMPS, &flg)); 407 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "-%smat_mumps_use_omp_threads and -%spc_factor_mat_solver_type != %s", prefix, prefix, MATSOLVERMUMPS); 408 size = n; 409 n = -1; 410 PetscCall(PetscOptionsGetInt(nullptr, prefix, "-mat_mumps_use_omp_threads", &n, nullptr)); 411 PetscCheck(n >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Need to specify a positive integer for -%smat_mumps_use_omp_threads", prefix); 412 PetscCheck(n * size <= previous, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "%d MPI process%s x %d OpenMP thread%s greater than %d available MPI process%s for the coarsest operator", (int)size, size > 1 ? "es" : "", (int)n, n > 1 ? "s" : "", (int)previous, previous > 1 ? "es" : ""); 413 } 414 #endif 415 PetscCall(PetscOptionsEnum("-pc_hpddm_coarse_correction", "Type of coarse correction applied each iteration", "PCHPDDMSetCoarseCorrectionType", PCHPDDMCoarseCorrectionTypes, (PetscEnum)data->correction, (PetscEnum *)&type, &flg)); 416 if (flg) PetscCall(PCHPDDMSetCoarseCorrectionType(pc, type)); 417 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_has_neumann")); 418 PetscCall(PetscOptionsBool(prefix, "Is the auxiliary Mat the local Neumann matrix?", "PCHPDDMHasNeumannMat", PetscBool3ToBool(data->Neumann), &flg, &set)); 419 if (set) data->Neumann = PetscBoolToBool3(flg); 420 data->log_separate = PETSC_FALSE; 421 if (PetscDefined(USE_LOG)) { 422 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_log_separate")); 423 PetscCall(PetscOptionsBool(prefix, "Log events level by level instead of inside PCSetUp()/KSPSolve()", nullptr, data->log_separate, &data->log_separate, nullptr)); 424 } 425 } 426 PetscOptionsHeadEnd(); 427 while (i < PETSC_PCHPDDM_MAXLEVELS && data->levels[i]) PetscCall(PetscFree(data->levels[i++])); 428 PetscFunctionReturn(PETSC_SUCCESS); 429 } 430 431 static PetscErrorCode PCApply_HPDDM(PC pc, Vec x, Vec y) 432 { 433 PC_HPDDM *data = (PC_HPDDM *)pc->data; 434 435 PetscFunctionBegin; 436 PetscCall(PetscCitationsRegister(HPDDMCitation, &HPDDMCite)); 437 PetscCheck(data->levels[0]->ksp, PETSC_COMM_SELF, PETSC_ERR_PLIB, "No KSP attached to PCHPDDM"); 438 if (data->log_separate) PetscCall(PetscLogEventBegin(PC_HPDDM_Solve[0], data->levels[0]->ksp, nullptr, nullptr, nullptr)); /* coarser-level events are directly triggered in HPDDM */ 439 PetscCall(KSPSolve(data->levels[0]->ksp, x, y)); 440 if (data->log_separate) PetscCall(PetscLogEventEnd(PC_HPDDM_Solve[0], data->levels[0]->ksp, nullptr, nullptr, nullptr)); 441 PetscFunctionReturn(PETSC_SUCCESS); 442 } 443 444 static PetscErrorCode PCMatApply_HPDDM(PC pc, Mat X, Mat Y) 445 { 446 PC_HPDDM *data = (PC_HPDDM *)pc->data; 447 448 PetscFunctionBegin; 449 PetscCall(PetscCitationsRegister(HPDDMCitation, &HPDDMCite)); 450 PetscCheck(data->levels[0]->ksp, PETSC_COMM_SELF, PETSC_ERR_PLIB, "No KSP attached to PCHPDDM"); 451 PetscCall(KSPMatSolve(data->levels[0]->ksp, X, Y)); 452 PetscFunctionReturn(PETSC_SUCCESS); 453 } 454 455 // PetscClangLinter pragma disable: -fdoc-internal-linkage 456 /*@C 457 PCHPDDMGetComplexities - Computes the grid and operator complexities. 458 459 Input Parameter: 460 . pc - preconditioner context 461 462 Output Parameters: 463 + gc - grid complexity = sum_i(m_i) / m_1 464 - oc - operator complexity = sum_i(nnz_i) / nnz_1 465 466 Level: advanced 467 468 .seealso: [](ch_ksp), `PCMGGetGridComplexity()`, `PCHPDDM`, `PCHYPRE`, `PCGAMG` 469 @*/ 470 static PetscErrorCode PCHPDDMGetComplexities(PC pc, PetscReal *gc, PetscReal *oc) 471 { 472 PC_HPDDM *data = (PC_HPDDM *)pc->data; 473 MatInfo info; 474 PetscInt n, m; 475 PetscLogDouble accumulate[2]{}, nnz1 = 1.0, m1 = 1.0; 476 477 PetscFunctionBegin; 478 for (n = 0, *gc = 0, *oc = 0; n < data->N; ++n) { 479 if (data->levels[n]->ksp) { 480 Mat P, A = nullptr; 481 PetscBool flg = PETSC_FALSE; 482 483 PetscCall(KSPGetOperators(data->levels[n]->ksp, nullptr, &P)); 484 PetscCall(MatGetSize(P, &m, nullptr)); 485 accumulate[0] += m; 486 if (n == 0) { 487 PetscCall(PetscObjectTypeCompareAny((PetscObject)P, &flg, MATNORMAL, MATNORMALHERMITIAN, "")); 488 if (flg) { 489 PetscCall(MatConvert(P, MATAIJ, MAT_INITIAL_MATRIX, &A)); 490 P = A; 491 } else { 492 PetscCall(PetscObjectTypeCompare((PetscObject)P, MATSCHURCOMPLEMENT, &flg)); 493 PetscCall(PetscObjectReference((PetscObject)P)); 494 } 495 } 496 if (!A && flg) accumulate[1] += m * m; /* assumption that a MATSCHURCOMPLEMENT is dense if stored explicitly */ 497 else if (P->ops->getinfo) { 498 PetscCall(MatGetInfo(P, MAT_GLOBAL_SUM, &info)); 499 accumulate[1] += info.nz_used; 500 } 501 if (n == 0) { 502 m1 = m; 503 if (!A && flg) nnz1 = m * m; 504 else if (P->ops->getinfo) nnz1 = info.nz_used; 505 PetscCall(MatDestroy(&P)); 506 } 507 } 508 } 509 *gc = static_cast<PetscReal>(accumulate[0] / m1); 510 *oc = static_cast<PetscReal>(accumulate[1] / nnz1); 511 PetscFunctionReturn(PETSC_SUCCESS); 512 } 513 514 static PetscErrorCode PCView_HPDDM(PC pc, PetscViewer viewer) 515 { 516 PC_HPDDM *data = (PC_HPDDM *)pc->data; 517 PetscViewer subviewer; 518 PetscViewerFormat format; 519 PetscSubcomm subcomm; 520 PetscReal oc, gc; 521 PetscInt i, tabs; 522 PetscMPIInt size, color, rank; 523 PetscBool flg; 524 const char *name; 525 526 PetscFunctionBegin; 527 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &flg)); 528 if (flg) { 529 PetscCall(PetscViewerASCIIPrintf(viewer, "level%s: %" PetscInt_FMT "\n", data->N > 1 ? "s" : "", data->N)); 530 PetscCall(PCHPDDMGetComplexities(pc, &gc, &oc)); 531 if (data->N > 1) { 532 if (!data->deflation) { 533 PetscCall(PetscViewerASCIIPrintf(viewer, "Neumann matrix attached? %s\n", PetscBools[PetscBool3ToBool(data->Neumann)])); 534 PetscCall(PetscViewerASCIIPrintf(viewer, "shared subdomain KSP between SLEPc and PETSc? %s\n", PetscBools[data->share])); 535 } else PetscCall(PetscViewerASCIIPrintf(viewer, "user-supplied deflation matrix\n")); 536 PetscCall(PetscViewerASCIIPrintf(viewer, "coarse correction: %s\n", PCHPDDMCoarseCorrectionTypes[data->correction])); 537 PetscCall(PetscViewerASCIIPrintf(viewer, "on process #0, value%s (+ threshold%s if available) for selecting deflation vectors:", data->N > 2 ? "s" : "", data->N > 2 ? "s" : "")); 538 PetscCall(PetscViewerASCIIGetTab(viewer, &tabs)); 539 PetscCall(PetscViewerASCIISetTab(viewer, 0)); 540 for (i = 1; i < data->N; ++i) { 541 PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, data->levels[i - 1]->nu)); 542 if (data->levels[i - 1]->threshold > static_cast<PetscReal>(-0.1)) PetscCall(PetscViewerASCIIPrintf(viewer, " (%g)", (double)data->levels[i - 1]->threshold)); 543 } 544 PetscCall(PetscViewerASCIIPrintf(viewer, "\n")); 545 PetscCall(PetscViewerASCIISetTab(viewer, tabs)); 546 } 547 PetscCall(PetscViewerASCIIPrintf(viewer, "grid and operator complexities: %g %g\n", (double)gc, (double)oc)); 548 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size)); 549 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank)); 550 if (data->levels[0]->ksp) { 551 PetscCall(KSPView(data->levels[0]->ksp, viewer)); 552 if (data->levels[0]->pc) PetscCall(PCView(data->levels[0]->pc, viewer)); 553 for (i = 1; i < data->N; ++i) { 554 if (data->levels[i]->ksp) color = 1; 555 else color = 0; 556 PetscCall(PetscSubcommCreate(PetscObjectComm((PetscObject)pc), &subcomm)); 557 PetscCall(PetscSubcommSetNumber(subcomm, PetscMin(size, 2))); 558 PetscCall(PetscSubcommSetTypeGeneral(subcomm, color, rank)); 559 PetscCall(PetscViewerASCIIPushTab(viewer)); 560 PetscCall(PetscViewerGetSubViewer(viewer, PetscSubcommChild(subcomm), &subviewer)); 561 if (color == 1) { 562 PetscCall(KSPView(data->levels[i]->ksp, subviewer)); 563 if (data->levels[i]->pc) PetscCall(PCView(data->levels[i]->pc, subviewer)); 564 PetscCall(PetscViewerFlush(subviewer)); 565 } 566 PetscCall(PetscViewerRestoreSubViewer(viewer, PetscSubcommChild(subcomm), &subviewer)); 567 PetscCall(PetscViewerASCIIPopTab(viewer)); 568 PetscCall(PetscSubcommDestroy(&subcomm)); 569 } 570 } 571 PetscCall(PetscViewerGetFormat(viewer, &format)); 572 if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 573 PetscCall(PetscViewerFileGetName(viewer, &name)); 574 if (name) { 575 IS is; 576 PetscInt sizes[4] = {pc->mat->rmap->n, pc->mat->cmap->n, pc->mat->rmap->N, pc->mat->cmap->N}; 577 char *tmp; 578 std::string prefix, suffix; 579 size_t pos; 580 581 PetscCall(PetscStrstr(name, ".", &tmp)); 582 if (tmp) { 583 pos = std::distance(const_cast<char *>(name), tmp); 584 prefix = std::string(name, pos); 585 suffix = std::string(name + pos + 1); 586 } else prefix = name; 587 if (data->aux) { 588 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_SELF, std::string(prefix + "_aux_" + std::to_string(rank) + "_" + std::to_string(size) + (tmp ? ("." + suffix) : "")).c_str(), FILE_MODE_WRITE, &subviewer)); 589 PetscCall(MatView(data->aux, subviewer)); 590 PetscCall(PetscViewerDestroy(&subviewer)); 591 } 592 if (data->is) { 593 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_SELF, std::string(prefix + "_is_" + std::to_string(rank) + "_" + std::to_string(size) + (tmp ? ("." + suffix) : "")).c_str(), FILE_MODE_WRITE, &subviewer)); 594 PetscCall(ISView(data->is, subviewer)); 595 PetscCall(PetscViewerDestroy(&subviewer)); 596 } 597 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, PETSC_STATIC_ARRAY_LENGTH(sizes), sizes, PETSC_USE_POINTER, &is)); 598 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_SELF, std::string(prefix + "_sizes_" + std::to_string(rank) + "_" + std::to_string(size) + (tmp ? ("." + suffix) : "")).c_str(), FILE_MODE_WRITE, &subviewer)); 599 PetscCall(ISView(is, subviewer)); 600 PetscCall(PetscViewerDestroy(&subviewer)); 601 PetscCall(ISDestroy(&is)); 602 } 603 } 604 } 605 PetscFunctionReturn(PETSC_SUCCESS); 606 } 607 608 static PetscErrorCode PCPreSolve_HPDDM(PC pc, KSP ksp, Vec, Vec) 609 { 610 PC_HPDDM *data = (PC_HPDDM *)pc->data; 611 PetscBool flg; 612 Mat A; 613 614 PetscFunctionBegin; 615 if (ksp) { 616 PetscCall(PetscObjectTypeCompare((PetscObject)ksp, KSPLSQR, &flg)); 617 if (flg && !data->normal) { 618 PetscCall(KSPGetOperators(ksp, &A, nullptr)); 619 PetscCall(MatCreateVecs(A, nullptr, &data->normal)); /* temporary Vec used in PCApply_HPDDMShell() for coarse grid corrections */ 620 } 621 } 622 PetscFunctionReturn(PETSC_SUCCESS); 623 } 624 625 static PetscErrorCode PCSetUp_HPDDMShell(PC pc) 626 { 627 PC_HPDDM_Level *ctx; 628 Mat A, P; 629 Vec x; 630 const char *pcpre; 631 632 PetscFunctionBegin; 633 PetscCall(PCShellGetContext(pc, &ctx)); 634 PetscCall(KSPGetOptionsPrefix(ctx->ksp, &pcpre)); 635 PetscCall(KSPGetOperators(ctx->ksp, &A, &P)); 636 /* smoother */ 637 PetscCall(PCSetOptionsPrefix(ctx->pc, pcpre)); 638 PetscCall(PCSetOperators(ctx->pc, A, P)); 639 if (!ctx->v[0]) { 640 PetscCall(VecDuplicateVecs(ctx->D, 1, &ctx->v[0])); 641 if (!std::is_same<PetscScalar, PetscReal>::value) PetscCall(VecDestroy(&ctx->D)); 642 PetscCall(MatCreateVecs(A, &x, nullptr)); 643 PetscCall(VecDuplicateVecs(x, 2, &ctx->v[1])); 644 PetscCall(VecDestroy(&x)); 645 } 646 std::fill_n(ctx->V, 3, nullptr); 647 PetscFunctionReturn(PETSC_SUCCESS); 648 } 649 650 template <class Type, typename std::enable_if<std::is_same<Type, Vec>::value>::type * = nullptr> 651 static inline PetscErrorCode PCHPDDMDeflate_Private(PC pc, Type x, Type y) 652 { 653 PC_HPDDM_Level *ctx; 654 PetscScalar *out; 655 656 PetscFunctionBegin; 657 PetscCall(PCShellGetContext(pc, &ctx)); 658 /* going from PETSc to HPDDM numbering */ 659 PetscCall(VecScatterBegin(ctx->scatter, x, ctx->v[0][0], INSERT_VALUES, SCATTER_FORWARD)); 660 PetscCall(VecScatterEnd(ctx->scatter, x, ctx->v[0][0], INSERT_VALUES, SCATTER_FORWARD)); 661 PetscCall(VecGetArrayWrite(ctx->v[0][0], &out)); 662 ctx->P->deflation<false>(nullptr, out, 1); /* y = Q x */ 663 PetscCall(VecRestoreArrayWrite(ctx->v[0][0], &out)); 664 /* going from HPDDM to PETSc numbering */ 665 PetscCall(VecScatterBegin(ctx->scatter, ctx->v[0][0], y, INSERT_VALUES, SCATTER_REVERSE)); 666 PetscCall(VecScatterEnd(ctx->scatter, ctx->v[0][0], y, INSERT_VALUES, SCATTER_REVERSE)); 667 PetscFunctionReturn(PETSC_SUCCESS); 668 } 669 670 template <class Type, typename std::enable_if<std::is_same<Type, Mat>::value>::type * = nullptr> 671 static inline PetscErrorCode PCHPDDMDeflate_Private(PC pc, Type X, Type Y) 672 { 673 PC_HPDDM_Level *ctx; 674 Vec vX, vY, vC; 675 PetscScalar *out; 676 PetscInt i, N; 677 678 PetscFunctionBegin; 679 PetscCall(PCShellGetContext(pc, &ctx)); 680 PetscCall(MatGetSize(X, nullptr, &N)); 681 /* going from PETSc to HPDDM numbering */ 682 for (i = 0; i < N; ++i) { 683 PetscCall(MatDenseGetColumnVecRead(X, i, &vX)); 684 PetscCall(MatDenseGetColumnVecWrite(ctx->V[0], i, &vC)); 685 PetscCall(VecScatterBegin(ctx->scatter, vX, vC, INSERT_VALUES, SCATTER_FORWARD)); 686 PetscCall(VecScatterEnd(ctx->scatter, vX, vC, INSERT_VALUES, SCATTER_FORWARD)); 687 PetscCall(MatDenseRestoreColumnVecWrite(ctx->V[0], i, &vC)); 688 PetscCall(MatDenseRestoreColumnVecRead(X, i, &vX)); 689 } 690 PetscCall(MatDenseGetArrayWrite(ctx->V[0], &out)); 691 ctx->P->deflation<false>(nullptr, out, N); /* Y = Q X */ 692 PetscCall(MatDenseRestoreArrayWrite(ctx->V[0], &out)); 693 /* going from HPDDM to PETSc numbering */ 694 for (i = 0; i < N; ++i) { 695 PetscCall(MatDenseGetColumnVecRead(ctx->V[0], i, &vC)); 696 PetscCall(MatDenseGetColumnVecWrite(Y, i, &vY)); 697 PetscCall(VecScatterBegin(ctx->scatter, vC, vY, INSERT_VALUES, SCATTER_REVERSE)); 698 PetscCall(VecScatterEnd(ctx->scatter, vC, vY, INSERT_VALUES, SCATTER_REVERSE)); 699 PetscCall(MatDenseRestoreColumnVecWrite(Y, i, &vY)); 700 PetscCall(MatDenseRestoreColumnVecRead(ctx->V[0], i, &vC)); 701 } 702 PetscFunctionReturn(PETSC_SUCCESS); 703 } 704 705 /* 706 PCApply_HPDDMShell - Applies a (2) deflated, (1) additive, (3) balanced, or (4) no coarse correction. In what follows, E = Z Pmat Z^T and Q = Z^T E^-1 Z. 707 708 .vb 709 (1) y = Pmat^-1 x + Q x, 710 (2) y = Pmat^-1 (I - Amat Q) x + Q x (default), 711 (3) y = (I - Q Amat^T) Pmat^-1 (I - Amat Q) x + Q x, 712 (4) y = Pmat^-1 x . 713 .ve 714 715 Input Parameters: 716 + pc - preconditioner context 717 - x - input vector 718 719 Output Parameter: 720 . y - output vector 721 722 Notes: 723 The options of Pmat^1 = pc(Pmat) are prefixed by -pc_hpddm_levels_1_pc_. Z is a tall-and-skiny matrix assembled by HPDDM. The number of processes on which (Z Pmat Z^T) is aggregated is set via -pc_hpddm_coarse_p. 724 The options of (Z Pmat Z^T)^-1 = ksp(Z Pmat Z^T) are prefixed by -pc_hpddm_coarse_ (`KSPPREONLY` and `PCCHOLESKY` by default), unless a multilevel correction is turned on, in which case, this function is called recursively at each level except the coarsest one. 725 (1) and (2) visit the "next" level (in terms of coarsening) once per application, while (3) visits it twice, so it is asymptotically twice costlier. (2) is not symmetric even if both Amat and Pmat are symmetric. 726 727 Level: advanced 728 729 Developer Note: 730 Since this is not an actual manual page the material below should be moved to an appropriate manual page with the appropriate context, i.e. explaining when it is used and how 731 to trigger it. Likely the manual page is `PCHPDDM` 732 733 .seealso: [](ch_ksp), `PCHPDDM`, `PCHPDDMCoarseCorrectionType` 734 */ 735 static PetscErrorCode PCApply_HPDDMShell(PC pc, Vec x, Vec y) 736 { 737 PC_HPDDM_Level *ctx; 738 Mat A; 739 740 PetscFunctionBegin; 741 PetscCall(PCShellGetContext(pc, &ctx)); 742 PetscCheck(ctx->P, PETSC_COMM_SELF, PETSC_ERR_PLIB, "PCSHELL from PCHPDDM called with no HPDDM object"); 743 PetscCall(KSPGetOperators(ctx->ksp, &A, nullptr)); 744 if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_NONE) PetscCall(PCApply(ctx->pc, x, y)); /* y = M^-1 x */ 745 else { 746 PetscCall(PCHPDDMDeflate_Private(pc, x, y)); /* y = Q x */ 747 if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_DEFLATED || ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_BALANCED) { 748 if (!ctx->parent->normal || ctx != ctx->parent->levels[0]) PetscCall(MatMult(A, y, ctx->v[1][0])); /* y = A Q x */ 749 else { /* KSPLSQR and finest level */ PetscCall(MatMult(A, y, ctx->parent->normal)); /* y = A Q x */ 750 PetscCall(MatMultHermitianTranspose(A, ctx->parent->normal, ctx->v[1][0])); /* y = A^T A Q x */ 751 } 752 PetscCall(VecWAXPY(ctx->v[1][1], -1.0, ctx->v[1][0], x)); /* y = (I - A Q) x */ 753 PetscCall(PCApply(ctx->pc, ctx->v[1][1], ctx->v[1][0])); /* y = M^-1 (I - A Q) x */ 754 if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_BALANCED) { 755 if (!ctx->parent->normal || ctx != ctx->parent->levels[0]) PetscCall(MatMultHermitianTranspose(A, ctx->v[1][0], ctx->v[1][1])); /* z = A^T y */ 756 else { 757 PetscCall(MatMult(A, ctx->v[1][0], ctx->parent->normal)); 758 PetscCall(MatMultHermitianTranspose(A, ctx->parent->normal, ctx->v[1][1])); /* z = A^T A y */ 759 } 760 PetscCall(PCHPDDMDeflate_Private(pc, ctx->v[1][1], ctx->v[1][1])); 761 PetscCall(VecAXPBYPCZ(y, -1.0, 1.0, 1.0, ctx->v[1][1], ctx->v[1][0])); /* y = (I - Q A^T) y + Q x */ 762 } else PetscCall(VecAXPY(y, 1.0, ctx->v[1][0])); /* y = Q M^-1 (I - A Q) x + Q x */ 763 } else { 764 PetscCheck(ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_ADDITIVE, PetscObjectComm((PetscObject)pc), PETSC_ERR_PLIB, "PCSHELL from PCHPDDM called with an unknown PCHPDDMCoarseCorrectionType %d", ctx->parent->correction); 765 PetscCall(PCApply(ctx->pc, x, ctx->v[1][0])); 766 PetscCall(VecAXPY(y, 1.0, ctx->v[1][0])); /* y = M^-1 x + Q x */ 767 } 768 } 769 PetscFunctionReturn(PETSC_SUCCESS); 770 } 771 772 /* 773 PCMatApply_HPDDMShell - Variant of PCApply_HPDDMShell() for blocks of vectors. 774 775 Input Parameters: 776 + pc - preconditioner context 777 - X - block of input vectors 778 779 Output Parameter: 780 . Y - block of output vectors 781 782 Level: advanced 783 784 .seealso: [](ch_ksp), `PCHPDDM`, `PCApply_HPDDMShell()`, `PCHPDDMCoarseCorrectionType` 785 */ 786 static PetscErrorCode PCMatApply_HPDDMShell(PC pc, Mat X, Mat Y) 787 { 788 PC_HPDDM_Level *ctx; 789 Mat A, *ptr; 790 PetscContainer container = nullptr; 791 PetscScalar *array; 792 PetscInt m, M, N, prev = 0; 793 PetscBool reset = PETSC_FALSE; 794 795 PetscFunctionBegin; 796 PetscCall(PCShellGetContext(pc, &ctx)); 797 PetscCheck(ctx->P, PETSC_COMM_SELF, PETSC_ERR_PLIB, "PCSHELL from PCHPDDM called with no HPDDM object"); 798 PetscCall(MatGetSize(X, nullptr, &N)); 799 PetscCall(KSPGetOperators(ctx->ksp, &A, nullptr)); 800 if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_NONE) PetscCall(PCMatApply(ctx->pc, X, Y)); 801 else { 802 PetscCall(PetscObjectQuery((PetscObject)A, "_HPDDM_MatProduct", (PetscObject *)&container)); 803 if (container) { /* MatProduct container already attached */ 804 PetscCall(PetscContainerGetPointer(container, (void **)&ptr)); 805 if (ptr[1] != ctx->V[2]) /* Mat has changed or may have been set first in KSPHPDDM */ 806 for (m = 0; m < 2; ++m) { 807 PetscCall(MatDestroy(ctx->V + m + 1)); 808 ctx->V[m + 1] = ptr[m]; 809 PetscCall(PetscObjectReference((PetscObject)ctx->V[m + 1])); 810 } 811 } 812 if (ctx->V[1]) PetscCall(MatGetSize(ctx->V[1], nullptr, &prev)); 813 if (N != prev || !ctx->V[0]) { 814 PetscCall(MatDestroy(ctx->V)); 815 PetscCall(VecGetLocalSize(ctx->v[0][0], &m)); 816 PetscCall(MatCreateDense(PetscObjectComm((PetscObject)pc), m, PETSC_DECIDE, PETSC_DECIDE, N, nullptr, ctx->V)); 817 if (N != prev) { 818 PetscCall(MatDestroy(ctx->V + 1)); 819 PetscCall(MatDestroy(ctx->V + 2)); 820 PetscCall(MatGetLocalSize(X, &m, nullptr)); 821 PetscCall(MatGetSize(X, &M, nullptr)); 822 if (ctx->parent->correction != PC_HPDDM_COARSE_CORRECTION_BALANCED) PetscCall(MatDenseGetArrayWrite(ctx->V[0], &array)); 823 else array = nullptr; 824 PetscCall(MatCreateDense(PetscObjectComm((PetscObject)pc), m, PETSC_DECIDE, M, N, array, ctx->V + 1)); 825 if (ctx->parent->correction != PC_HPDDM_COARSE_CORRECTION_BALANCED) PetscCall(MatDenseRestoreArrayWrite(ctx->V[0], &array)); 826 PetscCall(MatCreateDense(PetscObjectComm((PetscObject)pc), m, PETSC_DECIDE, M, N, nullptr, ctx->V + 2)); 827 PetscCall(MatProductCreateWithMat(A, Y, nullptr, ctx->V[1])); 828 PetscCall(MatProductSetType(ctx->V[1], MATPRODUCT_AB)); 829 PetscCall(MatProductSetFromOptions(ctx->V[1])); 830 PetscCall(MatProductSymbolic(ctx->V[1])); 831 if (!container) { /* no MatProduct container attached, create one to be queried in KSPHPDDM or at the next call to PCMatApply() */ 832 PetscCall(PetscContainerCreate(PetscObjectComm((PetscObject)A), &container)); 833 PetscCall(PetscObjectCompose((PetscObject)A, "_HPDDM_MatProduct", (PetscObject)container)); 834 } 835 PetscCall(PetscContainerSetPointer(container, ctx->V + 1)); /* need to compose B and D from MatProductCreateWithMat(A, B, NULL, D), which are stored in the contiguous array ctx->V */ 836 } 837 if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_BALANCED) { 838 PetscCall(MatProductCreateWithMat(A, ctx->V[1], nullptr, ctx->V[2])); 839 PetscCall(MatProductSetType(ctx->V[2], MATPRODUCT_AtB)); 840 PetscCall(MatProductSetFromOptions(ctx->V[2])); 841 PetscCall(MatProductSymbolic(ctx->V[2])); 842 } 843 ctx->P->start(N); 844 } 845 if (N == prev || container) { /* when MatProduct container is attached, always need to MatProductReplaceMats() since KSPHPDDM may have replaced the Mat as well */ 846 PetscCall(MatProductReplaceMats(nullptr, Y, nullptr, ctx->V[1])); 847 if (container && ctx->parent->correction != PC_HPDDM_COARSE_CORRECTION_BALANCED) { 848 PetscCall(MatDenseGetArrayWrite(ctx->V[0], &array)); 849 PetscCall(MatDensePlaceArray(ctx->V[1], array)); 850 PetscCall(MatDenseRestoreArrayWrite(ctx->V[0], &array)); 851 reset = PETSC_TRUE; 852 } 853 } 854 PetscCall(PCHPDDMDeflate_Private(pc, X, Y)); 855 if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_DEFLATED || ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_BALANCED) { 856 PetscCall(MatProductNumeric(ctx->V[1])); 857 PetscCall(MatCopy(ctx->V[1], ctx->V[2], SAME_NONZERO_PATTERN)); 858 PetscCall(MatAXPY(ctx->V[2], -1.0, X, SAME_NONZERO_PATTERN)); 859 PetscCall(PCMatApply(ctx->pc, ctx->V[2], ctx->V[1])); 860 if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_BALANCED) { 861 PetscCall(MatProductNumeric(ctx->V[2])); 862 PetscCall(PCHPDDMDeflate_Private(pc, ctx->V[2], ctx->V[2])); 863 PetscCall(MatAXPY(ctx->V[1], -1.0, ctx->V[2], SAME_NONZERO_PATTERN)); 864 } 865 PetscCall(MatAXPY(Y, -1.0, ctx->V[1], SAME_NONZERO_PATTERN)); 866 } else { 867 PetscCheck(ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_ADDITIVE, PetscObjectComm((PetscObject)pc), PETSC_ERR_PLIB, "PCSHELL from PCHPDDM called with an unknown PCHPDDMCoarseCorrectionType %d", ctx->parent->correction); 868 PetscCall(PCMatApply(ctx->pc, X, ctx->V[1])); 869 PetscCall(MatAXPY(Y, 1.0, ctx->V[1], SAME_NONZERO_PATTERN)); 870 } 871 if (reset) PetscCall(MatDenseResetArray(ctx->V[1])); 872 } 873 PetscFunctionReturn(PETSC_SUCCESS); 874 } 875 876 static PetscErrorCode PCDestroy_HPDDMShell(PC pc) 877 { 878 PC_HPDDM_Level *ctx; 879 PetscContainer container; 880 881 PetscFunctionBegin; 882 PetscCall(PCShellGetContext(pc, &ctx)); 883 PetscCall(HPDDM::Schwarz<PetscScalar>::destroy(ctx, PETSC_TRUE)); 884 PetscCall(VecDestroyVecs(1, &ctx->v[0])); 885 PetscCall(VecDestroyVecs(2, &ctx->v[1])); 886 PetscCall(PetscObjectQuery((PetscObject)(ctx->pc)->mat, "_HPDDM_MatProduct", (PetscObject *)&container)); 887 PetscCall(PetscContainerDestroy(&container)); 888 PetscCall(PetscObjectCompose((PetscObject)(ctx->pc)->mat, "_HPDDM_MatProduct", nullptr)); 889 PetscCall(MatDestroy(ctx->V)); 890 PetscCall(MatDestroy(ctx->V + 1)); 891 PetscCall(MatDestroy(ctx->V + 2)); 892 PetscCall(VecDestroy(&ctx->D)); 893 PetscCall(VecScatterDestroy(&ctx->scatter)); 894 PetscCall(PCDestroy(&ctx->pc)); 895 PetscFunctionReturn(PETSC_SUCCESS); 896 } 897 898 template <class Type, bool T = false, typename std::enable_if<std::is_same<Type, Vec>::value>::type * = nullptr> 899 static inline PetscErrorCode PCApply_Schur_Private(std::tuple<KSP, IS, Vec[2]> *p, PC factor, Type x, Type y) 900 { 901 PetscFunctionBegin; 902 PetscCall(VecISCopy(std::get<2>(*p)[0], std::get<1>(*p), SCATTER_FORWARD, x)); 903 if (!T) PetscCall(PCApply(factor, std::get<2>(*p)[0], std::get<2>(*p)[1])); 904 else PetscCall(PCApplyTranspose(factor, std::get<2>(*p)[0], std::get<2>(*p)[1])); 905 PetscCall(VecISCopy(std::get<2>(*p)[1], std::get<1>(*p), SCATTER_REVERSE, y)); 906 PetscFunctionReturn(PETSC_SUCCESS); 907 } 908 909 template <class Type, bool = false, typename std::enable_if<std::is_same<Type, Mat>::value>::type * = nullptr> 910 static inline PetscErrorCode PCApply_Schur_Private(std::tuple<KSP, IS, Vec[2]> *p, PC factor, Type X, Type Y) 911 { 912 Mat B[2]; 913 Vec x, y; 914 915 PetscFunctionBegin; 916 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, factor->mat->rmap->n, X->cmap->n, nullptr, B)); 917 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, factor->mat->rmap->n, X->cmap->n, nullptr, B + 1)); 918 for (PetscInt i = 0; i < X->cmap->n; ++i) { 919 PetscCall(MatDenseGetColumnVecRead(X, i, &x)); 920 PetscCall(MatDenseGetColumnVecWrite(B[0], i, &y)); 921 PetscCall(VecISCopy(y, std::get<1>(*p), SCATTER_FORWARD, x)); 922 PetscCall(MatDenseRestoreColumnVecWrite(B[0], i, &y)); 923 PetscCall(MatDenseRestoreColumnVecRead(X, i, &x)); 924 } 925 PetscCall(PCMatApply(factor, B[0], B[1])); 926 PetscCall(MatDestroy(B)); 927 for (PetscInt i = 0; i < X->cmap->n; ++i) { 928 PetscCall(MatDenseGetColumnVecRead(B[1], i, &x)); 929 PetscCall(MatDenseGetColumnVecWrite(Y, i, &y)); 930 PetscCall(VecISCopy(x, std::get<1>(*p), SCATTER_REVERSE, y)); 931 PetscCall(MatDenseRestoreColumnVecWrite(Y, i, &y)); 932 PetscCall(MatDenseRestoreColumnVecRead(B[1], i, &x)); 933 } 934 PetscCall(MatDestroy(B + 1)); 935 PetscFunctionReturn(PETSC_SUCCESS); 936 } 937 938 template <class Type = Vec, bool T = false> 939 static PetscErrorCode PCApply_Schur(PC pc, Type x, Type y) 940 { 941 PC factor; 942 Mat A; 943 MatSolverType type; 944 PetscBool flg; 945 std::tuple<KSP, IS, Vec[2]> *p; 946 947 PetscFunctionBegin; 948 PetscCall(PCShellGetContext(pc, &p)); 949 PetscCall(KSPGetPC(std::get<0>(*p), &factor)); 950 PetscCall(PCFactorGetMatSolverType(factor, &type)); 951 PetscCall(PCFactorGetMatrix(factor, &A)); 952 PetscCall(PetscStrcmp(type, MATSOLVERMUMPS, &flg)); 953 if (flg) { 954 PetscCheck(PetscDefined(HAVE_MUMPS), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent MatSolverType"); 955 #if PetscDefined(HAVE_MUMPS) 956 PetscCall(MatMumpsSetIcntl(A, 26, 0)); 957 #endif 958 } else { 959 PetscCall(PetscStrcmp(type, MATSOLVERMKL_PARDISO, &flg)); 960 PetscCheck(flg && PetscDefined(HAVE_MKL_PARDISO), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent MatSolverType"); 961 flg = PETSC_FALSE; 962 #if PetscDefined(HAVE_MKL_PARDISO) 963 PetscCall(MatMkl_PardisoSetCntl(A, 70, 1)); 964 #endif 965 } 966 PetscCall(PCApply_Schur_Private<Type, T>(p, factor, x, y)); 967 if (flg) { 968 #if PetscDefined(HAVE_MUMPS) 969 PetscCall(MatMumpsSetIcntl(A, 26, -1)); 970 #endif 971 } else { 972 #if PetscDefined(HAVE_MKL_PARDISO) 973 PetscCall(MatMkl_PardisoSetCntl(A, 70, 0)); 974 #endif 975 } 976 PetscFunctionReturn(PETSC_SUCCESS); 977 } 978 979 static PetscErrorCode PCDestroy_Schur(PC pc) 980 { 981 std::tuple<KSP, IS, Vec[2]> *p; 982 983 PetscFunctionBegin; 984 PetscCall(PCShellGetContext(pc, &p)); 985 PetscCall(ISDestroy(&std::get<1>(*p))); 986 PetscCall(VecDestroy(std::get<2>(*p))); 987 PetscCall(VecDestroy(std::get<2>(*p) + 1)); 988 PetscCall(PetscFree(p)); 989 PetscFunctionReturn(PETSC_SUCCESS); 990 } 991 992 static PetscErrorCode PCHPDDMSolve_Private(const PC_HPDDM_Level *ctx, PetscScalar *rhs, const unsigned short &mu) 993 { 994 Mat B, X; 995 PetscInt n, N, j = 0; 996 997 PetscFunctionBegin; 998 PetscCall(KSPGetOperators(ctx->ksp, &B, nullptr)); 999 PetscCall(MatGetLocalSize(B, &n, nullptr)); 1000 PetscCall(MatGetSize(B, &N, nullptr)); 1001 if (ctx->parent->log_separate) { 1002 j = std::distance(ctx->parent->levels, std::find(ctx->parent->levels, ctx->parent->levels + ctx->parent->N, ctx)); 1003 PetscCall(PetscLogEventBegin(PC_HPDDM_Solve[j], ctx->ksp, nullptr, nullptr, nullptr)); 1004 } 1005 if (mu == 1) { 1006 if (!ctx->ksp->vec_rhs) { 1007 PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)ctx->ksp), 1, n, N, nullptr, &ctx->ksp->vec_rhs)); 1008 PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)ctx->ksp), n, N, &ctx->ksp->vec_sol)); 1009 } 1010 PetscCall(VecPlaceArray(ctx->ksp->vec_rhs, rhs)); 1011 PetscCall(KSPSolve(ctx->ksp, nullptr, nullptr)); 1012 PetscCall(VecCopy(ctx->ksp->vec_sol, ctx->ksp->vec_rhs)); 1013 PetscCall(VecResetArray(ctx->ksp->vec_rhs)); 1014 } else { 1015 PetscCall(MatCreateDense(PetscObjectComm((PetscObject)ctx->ksp), n, PETSC_DECIDE, N, mu, rhs, &B)); 1016 PetscCall(MatCreateDense(PetscObjectComm((PetscObject)ctx->ksp), n, PETSC_DECIDE, N, mu, nullptr, &X)); 1017 PetscCall(KSPMatSolve(ctx->ksp, B, X)); 1018 PetscCall(MatCopy(X, B, SAME_NONZERO_PATTERN)); 1019 PetscCall(MatDestroy(&X)); 1020 PetscCall(MatDestroy(&B)); 1021 } 1022 if (ctx->parent->log_separate) PetscCall(PetscLogEventEnd(PC_HPDDM_Solve[j], ctx->ksp, nullptr, nullptr, nullptr)); 1023 PetscFunctionReturn(PETSC_SUCCESS); 1024 } 1025 1026 static PetscErrorCode PCHPDDMSetUpNeumannOverlap_Private(PC pc) 1027 { 1028 PC_HPDDM *data = (PC_HPDDM *)pc->data; 1029 1030 PetscFunctionBegin; 1031 if (data->setup) { 1032 Mat P; 1033 Vec x, xt = nullptr; 1034 PetscReal t = 0.0, s = 0.0; 1035 1036 PetscCall(PCGetOperators(pc, nullptr, &P)); 1037 PetscCall(PetscObjectQuery((PetscObject)P, "__SNES_latest_X", (PetscObject *)&x)); 1038 PetscCallBack("PCHPDDM Neumann callback", (*data->setup)(data->aux, t, x, xt, s, data->is, data->setup_ctx)); 1039 } 1040 PetscFunctionReturn(PETSC_SUCCESS); 1041 } 1042 1043 static PetscErrorCode PCHPDDMCreateSubMatrices_Private(Mat mat, PetscInt n, const IS *, const IS *, MatReuse scall, Mat *submat[]) 1044 { 1045 Mat A; 1046 PetscBool flg; 1047 1048 PetscFunctionBegin; 1049 PetscCheck(n == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "MatCreateSubMatrices() called to extract %" PetscInt_FMT " submatrices, which is different than 1", n); 1050 /* previously composed Mat */ 1051 PetscCall(PetscObjectQuery((PetscObject)mat, "_PCHPDDM_SubMatrices", (PetscObject *)&A)); 1052 PetscCheck(A, PETSC_COMM_SELF, PETSC_ERR_PLIB, "SubMatrices not found in Mat"); 1053 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATSCHURCOMPLEMENT, &flg)); /* MATSCHURCOMPLEMENT has neither a MatDuplicate() nor a MatCopy() implementation */ 1054 if (scall == MAT_INITIAL_MATRIX) { 1055 PetscCall(PetscCalloc1(2, submat)); /* allocate an extra Mat to avoid errors in MatDestroySubMatrices_Dummy() */ 1056 if (!flg) PetscCall(MatDuplicate(A, MAT_COPY_VALUES, *submat)); 1057 } else if (!flg) PetscCall(MatCopy(A, (*submat)[0], SAME_NONZERO_PATTERN)); 1058 if (flg) { 1059 PetscCall(MatDestroy(*submat)); /* previously created Mat has to be destroyed */ 1060 (*submat)[0] = A; 1061 PetscCall(PetscObjectReference((PetscObject)A)); 1062 } 1063 PetscFunctionReturn(PETSC_SUCCESS); 1064 } 1065 1066 static PetscErrorCode PCHPDDMCommunicationAvoidingPCASM_Private(PC pc, Mat C, PetscBool sorted) 1067 { 1068 void (*op)(void); 1069 1070 PetscFunctionBegin; 1071 /* previously-composed Mat */ 1072 PetscCall(PetscObjectCompose((PetscObject)pc->pmat, "_PCHPDDM_SubMatrices", (PetscObject)C)); 1073 PetscCall(MatGetOperation(pc->pmat, MATOP_CREATE_SUBMATRICES, &op)); 1074 /* trick suggested by Barry https://lists.mcs.anl.gov/pipermail/petsc-dev/2020-January/025491.html */ 1075 PetscCall(MatSetOperation(pc->pmat, MATOP_CREATE_SUBMATRICES, (void (*)(void))PCHPDDMCreateSubMatrices_Private)); 1076 if (sorted) PetscCall(PCASMSetSortIndices(pc, PETSC_FALSE)); /* everything is already sorted */ 1077 PetscCall(PCSetFromOptions(pc)); /* otherwise -pc_hpddm_levels_1_pc_asm_sub_mat_type is not used */ 1078 PetscCall(PCSetUp(pc)); 1079 /* reset MatCreateSubMatrices() */ 1080 PetscCall(MatSetOperation(pc->pmat, MATOP_CREATE_SUBMATRICES, op)); 1081 PetscCall(PetscObjectCompose((PetscObject)pc->pmat, "_PCHPDDM_SubMatrices", nullptr)); 1082 PetscFunctionReturn(PETSC_SUCCESS); 1083 } 1084 1085 static PetscErrorCode PCHPDDMPermute_Private(IS is, IS in_is, IS *out_is, Mat in_C, Mat *out_C, IS *p) 1086 { 1087 IS perm; 1088 const PetscInt *ptr; 1089 PetscInt *concatenate, size, n, bs; 1090 std::map<PetscInt, PetscInt> order; 1091 PetscBool sorted; 1092 1093 PetscFunctionBegin; 1094 PetscCall(ISSorted(is, &sorted)); 1095 if (!sorted) { 1096 PetscCall(ISGetLocalSize(is, &size)); 1097 PetscCall(ISGetIndices(is, &ptr)); 1098 PetscCall(ISGetBlockSize(is, &bs)); 1099 /* MatCreateSubMatrices(), called by PCASM, follows the global numbering of Pmat */ 1100 for (n = 0; n < size; n += bs) order.insert(std::make_pair(ptr[n] / bs, n / bs)); 1101 PetscCall(ISRestoreIndices(is, &ptr)); 1102 size /= bs; 1103 if (out_C) { 1104 PetscCall(PetscMalloc1(size, &concatenate)); 1105 for (const std::pair<const PetscInt, PetscInt> &i : order) *concatenate++ = i.second; 1106 concatenate -= size; 1107 PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)in_C), bs, size, concatenate, PETSC_OWN_POINTER, &perm)); 1108 PetscCall(ISSetPermutation(perm)); 1109 /* permute user-provided Mat so that it matches with MatCreateSubMatrices() numbering */ 1110 PetscCall(MatPermute(in_C, perm, perm, out_C)); 1111 if (p) *p = perm; 1112 else PetscCall(ISDestroy(&perm)); /* no need to save the permutation */ 1113 } 1114 if (out_is) { 1115 PetscCall(PetscMalloc1(size, &concatenate)); 1116 for (const std::pair<const PetscInt, PetscInt> &i : order) *concatenate++ = i.first; 1117 concatenate -= size; 1118 /* permute user-provided IS so that it matches with MatCreateSubMatrices() numbering */ 1119 PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)in_is), bs, size, concatenate, PETSC_OWN_POINTER, out_is)); 1120 } 1121 } else { /* input IS is sorted, nothing to permute, simply duplicate inputs when needed */ 1122 if (out_C) PetscCall(MatDuplicate(in_C, MAT_COPY_VALUES, out_C)); 1123 if (out_is) PetscCall(ISDuplicate(in_is, out_is)); 1124 } 1125 PetscFunctionReturn(PETSC_SUCCESS); 1126 } 1127 1128 static PetscErrorCode PCHPDDMCheckSymmetry_Private(PC pc, Mat A01, Mat A10) 1129 { 1130 Mat T, U = nullptr, B = nullptr; 1131 IS z; 1132 PetscBool flg[2]; 1133 1134 PetscFunctionBegin; 1135 PetscCall(PetscObjectTypeCompare((PetscObject)A10, MATTRANSPOSEVIRTUAL, flg)); 1136 if (flg[0]) PetscCall(MatTransposeGetMat(A10, &U)); 1137 else { 1138 PetscCall(PetscObjectTypeCompare((PetscObject)A10, MATHERMITIANTRANSPOSEVIRTUAL, flg + 1)); 1139 if (flg[1]) PetscCall(MatHermitianTransposeGetMat(A10, &U)); 1140 } 1141 if (U) PetscCall(MatDuplicate(U, MAT_COPY_VALUES, &T)); 1142 else PetscCall(MatHermitianTranspose(A10, MAT_INITIAL_MATRIX, &T)); 1143 PetscCall(PetscObjectTypeCompare((PetscObject)A01, MATTRANSPOSEVIRTUAL, flg)); 1144 if (flg[0]) { 1145 PetscCall(MatTransposeGetMat(A01, &A01)); 1146 PetscCall(MatTranspose(A01, MAT_INITIAL_MATRIX, &B)); 1147 A01 = B; 1148 } else { 1149 PetscCall(PetscObjectTypeCompare((PetscObject)A01, MATHERMITIANTRANSPOSEVIRTUAL, flg)); 1150 if (flg[0]) { 1151 PetscCall(MatHermitianTransposeGetMat(A01, &A01)); 1152 PetscCall(MatHermitianTranspose(A01, MAT_INITIAL_MATRIX, &B)); 1153 A01 = B; 1154 } 1155 } 1156 PetscCall(PetscLayoutCompare(T->rmap, A01->rmap, flg)); 1157 if (flg[0]) { 1158 PetscCall(PetscLayoutCompare(T->cmap, A01->cmap, flg)); 1159 if (flg[0]) { 1160 PetscCall(MatFindZeroRows(A01, &z)); /* for essential boundary conditions, some implementations will */ 1161 if (z) { /* zero rows in [P00 A01] except for the diagonal of P00 */ 1162 PetscCall(MatSetOption(T, MAT_NO_OFF_PROC_ZERO_ROWS, PETSC_TRUE)); 1163 PetscCall(MatZeroRowsIS(T, z, 0.0, nullptr, nullptr)); /* corresponding zero rows from A01 */ 1164 PetscCall(ISDestroy(&z)); 1165 } 1166 PetscCall(MatMultEqual(A01, T, 20, flg)); 1167 PetscCheck(flg[0], PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "A01 != A10^T"); 1168 } else PetscCall(PetscInfo(pc, "A01 and A10^T have non-congruent column layouts, cannot test for equality\n")); 1169 } 1170 PetscCall(MatDestroy(&B)); 1171 PetscCall(MatDestroy(&T)); 1172 PetscFunctionReturn(PETSC_SUCCESS); 1173 } 1174 1175 static PetscErrorCode PCHPDDMDestroySubMatrices_Private(PetscBool flg, PetscBool algebraic, Mat *sub) 1176 { 1177 IS is; 1178 1179 PetscFunctionBegin; 1180 if (!flg) { 1181 if (algebraic) { 1182 PetscCall(PetscObjectQuery((PetscObject)sub[0], "_PCHPDDM_Embed", (PetscObject *)&is)); 1183 PetscCall(ISDestroy(&is)); 1184 PetscCall(PetscObjectCompose((PetscObject)sub[0], "_PCHPDDM_Embed", nullptr)); 1185 PetscCall(PetscObjectCompose((PetscObject)sub[0], "_PCHPDDM_Compact", nullptr)); 1186 } 1187 PetscCall(MatDestroySubMatrices(algebraic ? 2 : 1, &sub)); 1188 } 1189 PetscFunctionReturn(PETSC_SUCCESS); 1190 } 1191 1192 static PetscErrorCode PCHPDDMAlgebraicAuxiliaryMat_Private(Mat P, IS *is, Mat *sub[], PetscBool block) 1193 { 1194 IS icol[3], irow[2]; 1195 Mat *M, Q; 1196 PetscReal *ptr; 1197 PetscInt *idx, p = 0, n, bs = PetscAbs(P->cmap->bs); 1198 PetscBool flg; 1199 1200 PetscFunctionBegin; 1201 PetscCall(ISCreateStride(PETSC_COMM_SELF, P->cmap->N, 0, 1, icol + 2)); 1202 PetscCall(ISSetBlockSize(icol[2], bs)); 1203 PetscCall(ISSetIdentity(icol[2])); 1204 PetscCall(PetscObjectTypeCompare((PetscObject)P, MATMPISBAIJ, &flg)); 1205 if (flg) { 1206 /* MatCreateSubMatrices() does not handle MATMPISBAIJ properly when iscol != isrow, so convert first to MATMPIBAIJ */ 1207 PetscCall(MatConvert(P, MATMPIBAIJ, MAT_INITIAL_MATRIX, &Q)); 1208 std::swap(P, Q); 1209 } 1210 PetscCall(MatCreateSubMatrices(P, 1, is, icol + 2, MAT_INITIAL_MATRIX, &M)); 1211 if (flg) { 1212 std::swap(P, Q); 1213 PetscCall(MatDestroy(&Q)); 1214 } 1215 PetscCall(ISDestroy(icol + 2)); 1216 PetscCall(ISCreateStride(PETSC_COMM_SELF, M[0]->rmap->N, 0, 1, irow)); 1217 PetscCall(ISSetBlockSize(irow[0], bs)); 1218 PetscCall(ISSetIdentity(irow[0])); 1219 if (!block) { 1220 PetscCall(PetscMalloc2(P->cmap->N, &ptr, P->cmap->N / bs, &idx)); 1221 PetscCall(MatGetColumnNorms(M[0], NORM_INFINITY, ptr)); 1222 /* check for nonzero columns so that M[0] may be expressed in compact form */ 1223 for (n = 0; n < P->cmap->N; n += bs) { 1224 if (std::find_if(ptr + n, ptr + n + bs, [](PetscReal v) { return v > PETSC_MACHINE_EPSILON; }) != ptr + n + bs) idx[p++] = n / bs; 1225 } 1226 PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, p, idx, PETSC_USE_POINTER, icol + 1)); 1227 PetscCall(ISSetInfo(icol[1], IS_SORTED, IS_GLOBAL, PETSC_TRUE, PETSC_TRUE)); 1228 PetscCall(ISEmbed(*is, icol[1], PETSC_FALSE, icol + 2)); 1229 irow[1] = irow[0]; 1230 /* first Mat will be used in PCASM (if it is used as a PC on this level) and as the left-hand side of GenEO */ 1231 icol[0] = is[0]; 1232 PetscCall(MatCreateSubMatrices(M[0], 2, irow, icol, MAT_INITIAL_MATRIX, sub)); 1233 PetscCall(ISDestroy(icol + 1)); 1234 PetscCall(PetscFree2(ptr, idx)); 1235 /* IS used to go back and forth between the augmented and the original local linear system, see eq. (3.4) of [2022b] */ 1236 PetscCall(PetscObjectCompose((PetscObject)(*sub)[0], "_PCHPDDM_Embed", (PetscObject)icol[2])); 1237 /* Mat used in eq. (3.1) of [2022b] */ 1238 PetscCall(PetscObjectCompose((PetscObject)(*sub)[0], "_PCHPDDM_Compact", (PetscObject)(*sub)[1])); 1239 } else { 1240 Mat aux; 1241 PetscCall(MatSetOption(M[0], MAT_SUBMAT_SINGLEIS, PETSC_TRUE)); 1242 /* diagonal block of the overlapping rows */ 1243 PetscCall(MatCreateSubMatrices(M[0], 1, irow, is, MAT_INITIAL_MATRIX, sub)); 1244 PetscCall(MatDuplicate((*sub)[0], MAT_COPY_VALUES, &aux)); 1245 PetscCall(MatSetOption(aux, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE)); 1246 if (bs == 1) { /* scalar case */ 1247 Vec sum[2]; 1248 PetscCall(MatCreateVecs(aux, sum, sum + 1)); 1249 PetscCall(MatGetRowSum(M[0], sum[0])); 1250 PetscCall(MatGetRowSum(aux, sum[1])); 1251 /* off-diagonal block row sum (full rows - diagonal block rows) */ 1252 PetscCall(VecAXPY(sum[0], -1.0, sum[1])); 1253 /* subdomain matrix plus off-diagonal block row sum */ 1254 PetscCall(MatDiagonalSet(aux, sum[0], ADD_VALUES)); 1255 PetscCall(VecDestroy(sum)); 1256 PetscCall(VecDestroy(sum + 1)); 1257 } else { /* vectorial case */ 1258 /* TODO: missing MatGetValuesBlocked(), so the code below is */ 1259 /* an extension of the scalar case for when bs > 1, but it could */ 1260 /* be more efficient by avoiding all these MatMatMult() */ 1261 Mat sum[2], ones; 1262 PetscScalar *ptr; 1263 PetscCall(PetscCalloc1(M[0]->cmap->n * bs, &ptr)); 1264 PetscCall(MatCreateDense(PETSC_COMM_SELF, M[0]->cmap->n, bs, M[0]->cmap->n, bs, ptr, &ones)); 1265 for (n = 0; n < M[0]->cmap->n; n += bs) { 1266 for (p = 0; p < bs; ++p) ptr[n + p * (M[0]->cmap->n + 1)] = 1.0; 1267 } 1268 PetscCall(MatMatMult(M[0], ones, MAT_INITIAL_MATRIX, PETSC_DEFAULT, sum)); 1269 PetscCall(MatDestroy(&ones)); 1270 PetscCall(MatCreateDense(PETSC_COMM_SELF, aux->cmap->n, bs, aux->cmap->n, bs, ptr, &ones)); 1271 PetscCall(MatDenseSetLDA(ones, M[0]->cmap->n)); 1272 PetscCall(MatMatMult(aux, ones, MAT_INITIAL_MATRIX, PETSC_DEFAULT, sum + 1)); 1273 PetscCall(MatDestroy(&ones)); 1274 PetscCall(PetscFree(ptr)); 1275 /* off-diagonal block row sum (full rows - diagonal block rows) */ 1276 PetscCall(MatAXPY(sum[0], -1.0, sum[1], SAME_NONZERO_PATTERN)); 1277 PetscCall(MatDestroy(sum + 1)); 1278 /* re-order values to be consistent with MatSetValuesBlocked() */ 1279 /* equivalent to MatTranspose() which does not truly handle */ 1280 /* MAT_INPLACE_MATRIX in the rectangular case, as it calls PetscMalloc() */ 1281 PetscCall(MatDenseGetArrayWrite(sum[0], &ptr)); 1282 HPDDM::Wrapper<PetscScalar>::imatcopy<'T'>(bs, sum[0]->rmap->n, ptr, sum[0]->rmap->n, bs); 1283 /* subdomain matrix plus off-diagonal block row sum */ 1284 for (n = 0; n < aux->cmap->n / bs; ++n) PetscCall(MatSetValuesBlocked(aux, 1, &n, 1, &n, ptr + n * bs * bs, ADD_VALUES)); 1285 PetscCall(MatAssemblyBegin(aux, MAT_FINAL_ASSEMBLY)); 1286 PetscCall(MatAssemblyEnd(aux, MAT_FINAL_ASSEMBLY)); 1287 PetscCall(MatDenseRestoreArrayWrite(sum[0], &ptr)); 1288 PetscCall(MatDestroy(sum)); 1289 } 1290 PetscCall(MatSetOption(aux, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE)); 1291 /* left-hand side of GenEO, with the same sparsity pattern as PCASM subdomain solvers */ 1292 PetscCall(PetscObjectCompose((PetscObject)(*sub)[0], "_PCHPDDM_Neumann_Mat", (PetscObject)aux)); 1293 } 1294 PetscCall(ISDestroy(irow)); 1295 PetscCall(MatDestroySubMatrices(1, &M)); 1296 PetscFunctionReturn(PETSC_SUCCESS); 1297 } 1298 1299 static PetscErrorCode PCApply_Nest(PC pc, Vec x, Vec y) 1300 { 1301 Mat A; 1302 MatSolverType type; 1303 IS is[2]; 1304 PetscBool flg; 1305 std::pair<PC, Vec[2]> *p; 1306 1307 PetscFunctionBegin; 1308 PetscCall(PCShellGetContext(pc, &p)); 1309 PetscCall(PCGetOperators(p->first, &A, nullptr)); 1310 PetscCall(MatNestGetISs(A, is, nullptr)); 1311 PetscCall(PCFactorGetMatSolverType(p->first, &type)); 1312 PetscCall(PCFactorGetMatrix(p->first, &A)); 1313 PetscCall(PetscStrcmp(type, MATSOLVERMUMPS, &flg)); 1314 if (flg && A->schur) { 1315 #if PetscDefined(HAVE_MUMPS) 1316 PetscCall(MatMumpsSetIcntl(A, 26, 1)); /* reduction/condensation phase followed by Schur complement solve */ 1317 #endif 1318 } 1319 PetscCall(VecISCopy(p->second[0], is[1], SCATTER_FORWARD, x)); /* assign the RHS associated to the Schur complement */ 1320 PetscCall(PCApply(p->first, p->second[0], p->second[1])); 1321 PetscCall(VecISCopy(p->second[1], is[1], SCATTER_REVERSE, y)); /* retrieve the partial solution associated to the Schur complement */ 1322 if (flg) { 1323 #if PetscDefined(HAVE_MUMPS) 1324 PetscCall(MatMumpsSetIcntl(A, 26, -1)); /* default ICNTL(26) value in PETSc */ 1325 #endif 1326 } 1327 PetscFunctionReturn(PETSC_SUCCESS); 1328 } 1329 1330 static PetscErrorCode PCView_Nest(PC pc, PetscViewer viewer) 1331 { 1332 std::pair<PC, Vec[2]> *p; 1333 1334 PetscFunctionBegin; 1335 PetscCall(PCShellGetContext(pc, &p)); 1336 PetscCall(PCView(p->first, viewer)); 1337 PetscFunctionReturn(PETSC_SUCCESS); 1338 } 1339 1340 static PetscErrorCode PCDestroy_Nest(PC pc) 1341 { 1342 std::pair<PC, Vec[2]> *p; 1343 1344 PetscFunctionBegin; 1345 PetscCall(PCShellGetContext(pc, &p)); 1346 PetscCall(VecDestroy(p->second)); 1347 PetscCall(VecDestroy(p->second + 1)); 1348 PetscCall(PCDestroy(&p->first)); 1349 PetscCall(PetscFree(p)); 1350 PetscFunctionReturn(PETSC_SUCCESS); 1351 } 1352 1353 template <bool T = false> 1354 static PetscErrorCode MatMult_Schur(Mat A, Vec x, Vec y) 1355 { 1356 std::tuple<Mat, VecScatter, Vec[2]> *ctx; 1357 1358 PetscFunctionBegin; 1359 PetscCall(MatShellGetContext(A, &ctx)); 1360 PetscCall(VecScatterBegin(std::get<1>(*ctx), x, std::get<2>(*ctx)[0], INSERT_VALUES, SCATTER_FORWARD)); /* local Vec with overlap */ 1361 PetscCall(VecScatterEnd(std::get<1>(*ctx), x, std::get<2>(*ctx)[0], INSERT_VALUES, SCATTER_FORWARD)); 1362 if (!T) PetscCall(MatMult(std::get<0>(*ctx), std::get<2>(*ctx)[0], std::get<2>(*ctx)[1])); /* local Schur complement */ 1363 else PetscCall(MatMultTranspose(std::get<0>(*ctx), std::get<2>(*ctx)[0], std::get<2>(*ctx)[1])); 1364 PetscCall(VecSet(y, 0.0)); 1365 PetscCall(VecScatterBegin(std::get<1>(*ctx), std::get<2>(*ctx)[1], y, ADD_VALUES, SCATTER_REVERSE)); /* global Vec with summed up contributions on the overlap */ 1366 PetscCall(VecScatterEnd(std::get<1>(*ctx), std::get<2>(*ctx)[1], y, ADD_VALUES, SCATTER_REVERSE)); 1367 PetscFunctionReturn(PETSC_SUCCESS); 1368 } 1369 1370 static PetscErrorCode MatDestroy_Schur(Mat A) 1371 { 1372 std::tuple<Mat, VecScatter, Vec[2]> *ctx; 1373 1374 PetscFunctionBegin; 1375 PetscCall(MatShellGetContext(A, &ctx)); 1376 PetscCall(VecDestroy(std::get<2>(*ctx))); 1377 PetscCall(VecDestroy(std::get<2>(*ctx) + 1)); 1378 PetscCall(PetscFree(ctx)); 1379 PetscFunctionReturn(PETSC_SUCCESS); 1380 } 1381 1382 static PetscErrorCode MatMult_SchurCorrection(Mat A, Vec x, Vec y) 1383 { 1384 PC pc; 1385 std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *ctx; 1386 1387 PetscFunctionBegin; 1388 PetscCall(MatShellGetContext(A, &ctx)); 1389 pc = ((PC_HPDDM *)std::get<0>(*ctx)[0]->data)->levels[0]->ksp->pc; 1390 if (std::get<2>(*ctx) == PC_LEFT || std::get<2>(*ctx) == PC_SIDE_DEFAULT) { /* Q_0 is the coarse correction associated to the A00 block from PCFIELDSPLIT */ 1391 PetscCall(MatMult(std::get<1>(*ctx)[0], x, std::get<3>(*ctx)[1])); /* A_01 x */ 1392 PetscCall(PCHPDDMDeflate_Private(pc, std::get<3>(*ctx)[1], std::get<3>(*ctx)[1])); /* Q_0 A_01 x */ 1393 PetscCall(MatMult(std::get<1>(*ctx)[1], std::get<3>(*ctx)[1], std::get<3>(*ctx)[0])); /* A_10 Q_0 A_01 x */ 1394 PetscCall(PCApply(std::get<0>(*ctx)[1], std::get<3>(*ctx)[0], y)); /* y = M_S^-1 A_10 Q_0 A_01 x */ 1395 } else { 1396 PetscCall(PCApply(std::get<0>(*ctx)[1], x, std::get<3>(*ctx)[0])); /* M_S^-1 x */ 1397 PetscCall(MatMult(std::get<1>(*ctx)[0], std::get<3>(*ctx)[0], std::get<3>(*ctx)[1])); /* A_01 M_S^-1 x */ 1398 PetscCall(PCHPDDMDeflate_Private(pc, std::get<3>(*ctx)[1], std::get<3>(*ctx)[1])); /* Q_0 A_01 M_S^-1 x */ 1399 PetscCall(MatMult(std::get<1>(*ctx)[1], std::get<3>(*ctx)[1], y)); /* y = A_10 Q_0 A_01 M_S^-1 x */ 1400 } 1401 PetscCall(VecAXPY(y, -1.0, x)); /* y -= x, preconditioned eq. (24) of https://hal.science/hal-02343808v6/document (with a sign flip) */ 1402 PetscFunctionReturn(PETSC_SUCCESS); 1403 } 1404 1405 static PetscErrorCode MatView_SchurCorrection(Mat A, PetscViewer viewer) 1406 { 1407 PetscBool ascii; 1408 std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *ctx; 1409 1410 PetscFunctionBegin; 1411 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &ascii)); 1412 if (ascii) { 1413 PetscCall(MatShellGetContext(A, &ctx)); 1414 PetscCall(PetscViewerASCIIPrintf(viewer, "action of %s\n", std::get<2>(*ctx) == PC_LEFT || std::get<2>(*ctx) == PC_SIDE_DEFAULT ? "(I - M_S^-1 A_10 Q_0 A_01)" : "(I - A_10 Q_0 A_01 M_S^-1)")); 1415 PetscCall(PCView(std::get<0>(*ctx)[1], viewer)); /* no need to PCView(Q_0) since it will be done by PCFIELDSPLIT */ 1416 } 1417 PetscFunctionReturn(PETSC_SUCCESS); 1418 } 1419 1420 static PetscErrorCode MatDestroy_SchurCorrection(Mat A) 1421 { 1422 std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *ctx; 1423 1424 PetscFunctionBegin; 1425 PetscCall(MatShellGetContext(A, &ctx)); 1426 PetscCall(VecDestroy(std::get<3>(*ctx))); 1427 PetscCall(VecDestroy(std::get<3>(*ctx) + 1)); 1428 PetscCall(VecDestroy(std::get<3>(*ctx) + 2)); 1429 PetscCall(PCDestroy(std::get<0>(*ctx) + 1)); 1430 PetscCall(PetscFree(ctx)); 1431 PetscFunctionReturn(PETSC_SUCCESS); 1432 } 1433 1434 static PetscErrorCode KSPPreSolve_SchurCorrection(KSP, Vec b, Vec, void *context) 1435 { 1436 std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *ctx = reinterpret_cast<std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *>(context); 1437 1438 PetscFunctionBegin; 1439 if (std::get<2>(*ctx) == PC_LEFT || std::get<2>(*ctx) == PC_SIDE_DEFAULT) { 1440 PetscCall(PCApply(std::get<0>(*ctx)[1], b, std::get<3>(*ctx)[2])); 1441 std::swap(*b, *std::get<3>(*ctx)[2]); /* replace b by M^-1 b, but need to keep a copy of the original RHS, so swap it with the work Vec */ 1442 } 1443 PetscFunctionReturn(PETSC_SUCCESS); 1444 } 1445 1446 static PetscErrorCode KSPPostSolve_SchurCorrection(KSP, Vec b, Vec x, void *context) 1447 { 1448 std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *ctx = reinterpret_cast<std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *>(context); 1449 1450 PetscFunctionBegin; 1451 if (std::get<2>(*ctx) == PC_LEFT || std::get<2>(*ctx) == PC_SIDE_DEFAULT) std::swap(*b, *std::get<3>(*ctx)[2]); /* put back the original RHS where it belongs */ 1452 else { 1453 PetscCall(PCApply(std::get<0>(*ctx)[1], x, std::get<3>(*ctx)[2])); 1454 PetscCall(VecCopy(std::get<3>(*ctx)[2], x)); /* replace x by M^-1 x */ 1455 } 1456 PetscFunctionReturn(PETSC_SUCCESS); 1457 } 1458 1459 static PetscErrorCode MatMult_Harmonic(Mat, Vec, Vec); 1460 static PetscErrorCode MatMultTranspose_Harmonic(Mat, Vec, Vec); 1461 static PetscErrorCode MatProduct_AB_Harmonic(Mat, Mat, Mat, void *); 1462 static PetscErrorCode MatDestroy_Harmonic(Mat); 1463 1464 static PetscErrorCode PCSetUp_HPDDM(PC pc) 1465 { 1466 PC_HPDDM *data = (PC_HPDDM *)pc->data; 1467 PC inner; 1468 KSP *ksp; 1469 Mat *sub, A, P, N, C = nullptr, uaux = nullptr, weighted, subA[2], S; 1470 Vec xin, v; 1471 std::vector<Vec> initial; 1472 IS is[1], loc, uis = data->is, unsorted = nullptr; 1473 ISLocalToGlobalMapping l2g; 1474 char prefix[256]; 1475 const char *pcpre; 1476 const PetscScalar *const *ev; 1477 PetscInt n, requested = data->N, reused = 0, overlap = -1; 1478 MatStructure structure = UNKNOWN_NONZERO_PATTERN; 1479 PetscBool subdomains = PETSC_FALSE, flg = PETSC_FALSE, ismatis, swap = PETSC_FALSE, algebraic = PETSC_FALSE, block = PETSC_FALSE; 1480 DM dm; 1481 std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *ctx = nullptr; 1482 #if PetscDefined(USE_DEBUG) 1483 IS dis = nullptr; 1484 Mat daux = nullptr; 1485 #endif 1486 1487 PetscFunctionBegin; 1488 PetscCheck(data->levels && data->levels[0], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Not a single level allocated"); 1489 PetscCall(PCGetOptionsPrefix(pc, &pcpre)); 1490 PetscCall(PCGetOperators(pc, &A, &P)); 1491 if (!data->levels[0]->ksp) { 1492 PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &data->levels[0]->ksp)); 1493 PetscCall(KSPSetNestLevel(data->levels[0]->ksp, pc->kspnestlevel)); 1494 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_%s_", pcpre ? pcpre : "", data->N > 1 ? "levels_1" : "coarse")); 1495 PetscCall(KSPSetOptionsPrefix(data->levels[0]->ksp, prefix)); 1496 PetscCall(KSPSetType(data->levels[0]->ksp, KSPPREONLY)); 1497 } else if (data->levels[0]->ksp->pc && data->levels[0]->ksp->pc->setupcalled == 1 && data->levels[0]->ksp->pc->reusepreconditioner) { 1498 /* if the fine-level PCSHELL exists, its setup has succeeded, and one wants to reuse it, */ 1499 /* then just propagate the appropriate flag to the coarser levels */ 1500 for (n = 0; n < PETSC_PCHPDDM_MAXLEVELS && data->levels[n]; ++n) { 1501 /* the following KSP and PC may be NULL for some processes, hence the check */ 1502 if (data->levels[n]->ksp) PetscCall(KSPSetReusePreconditioner(data->levels[n]->ksp, PETSC_TRUE)); 1503 if (data->levels[n]->pc) PetscCall(PCSetReusePreconditioner(data->levels[n]->pc, PETSC_TRUE)); 1504 } 1505 /* early bail out because there is nothing to do */ 1506 PetscFunctionReturn(PETSC_SUCCESS); 1507 } else { 1508 /* reset coarser levels */ 1509 for (n = 1; n < PETSC_PCHPDDM_MAXLEVELS && data->levels[n]; ++n) { 1510 if (data->levels[n]->ksp && data->levels[n]->ksp->pc && data->levels[n]->ksp->pc->setupcalled == 1 && data->levels[n]->ksp->pc->reusepreconditioner && n < data->N) { 1511 reused = data->N - n; 1512 break; 1513 } 1514 PetscCall(KSPDestroy(&data->levels[n]->ksp)); 1515 PetscCall(PCDestroy(&data->levels[n]->pc)); 1516 } 1517 /* check if some coarser levels are being reused */ 1518 PetscCall(MPIU_Allreduce(MPI_IN_PLACE, &reused, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)pc))); 1519 const int *addr = data->levels[0]->P ? data->levels[0]->P->getAddrLocal() : &HPDDM::i__0; 1520 1521 if (addr != &HPDDM::i__0 && reused != data->N - 1) { 1522 /* reuse previously computed eigenvectors */ 1523 ev = data->levels[0]->P->getVectors(); 1524 if (ev) { 1525 initial.reserve(*addr); 1526 PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, data->levels[0]->P->getDof(), ev[0], &xin)); 1527 for (n = 0; n < *addr; ++n) { 1528 PetscCall(VecDuplicate(xin, &v)); 1529 PetscCall(VecPlaceArray(xin, ev[n])); 1530 PetscCall(VecCopy(xin, v)); 1531 initial.emplace_back(v); 1532 PetscCall(VecResetArray(xin)); 1533 } 1534 PetscCall(VecDestroy(&xin)); 1535 } 1536 } 1537 } 1538 data->N -= reused; 1539 PetscCall(KSPSetOperators(data->levels[0]->ksp, A, P)); 1540 1541 PetscCall(PetscObjectTypeCompare((PetscObject)P, MATIS, &ismatis)); 1542 if (!data->is && !ismatis) { 1543 PetscErrorCode (*create)(DM, IS *, Mat *, PetscErrorCode (**)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void *), void **) = nullptr; 1544 PetscErrorCode (*usetup)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void *) = nullptr; 1545 void *uctx = nullptr; 1546 1547 /* first see if we can get the data from the DM */ 1548 PetscCall(MatGetDM(P, &dm)); 1549 if (!dm) PetscCall(MatGetDM(A, &dm)); 1550 if (!dm) PetscCall(PCGetDM(pc, &dm)); 1551 if (dm) { /* this is the hook for DMPLEX for which the auxiliary Mat is the local Neumann matrix */ 1552 PetscCall(PetscObjectQueryFunction((PetscObject)dm, "DMCreateNeumannOverlap_C", &create)); 1553 if (create) { 1554 PetscCall((*create)(dm, &uis, &uaux, &usetup, &uctx)); 1555 if (data->Neumann == PETSC_BOOL3_UNKNOWN) data->Neumann = PETSC_BOOL3_TRUE; /* set the value only if it was not already provided by the user */ 1556 } 1557 } 1558 if (!create) { 1559 if (!uis) { 1560 PetscCall(PetscObjectQuery((PetscObject)pc, "_PCHPDDM_Neumann_IS", (PetscObject *)&uis)); 1561 PetscCall(PetscObjectReference((PetscObject)uis)); 1562 } 1563 if (!uaux) { 1564 PetscCall(PetscObjectQuery((PetscObject)pc, "_PCHPDDM_Neumann_Mat", (PetscObject *)&uaux)); 1565 PetscCall(PetscObjectReference((PetscObject)uaux)); 1566 } 1567 /* look inside the Pmat instead of the PC, needed for MatSchurComplementComputeExplicitOperator() */ 1568 if (!uis) { 1569 PetscCall(PetscObjectQuery((PetscObject)P, "_PCHPDDM_Neumann_IS", (PetscObject *)&uis)); 1570 PetscCall(PetscObjectReference((PetscObject)uis)); 1571 } 1572 if (!uaux) { 1573 PetscCall(PetscObjectQuery((PetscObject)P, "_PCHPDDM_Neumann_Mat", (PetscObject *)&uaux)); 1574 PetscCall(PetscObjectReference((PetscObject)uaux)); 1575 } 1576 } 1577 PetscCall(PCHPDDMSetAuxiliaryMat(pc, uis, uaux, usetup, uctx)); 1578 PetscCall(MatDestroy(&uaux)); 1579 PetscCall(ISDestroy(&uis)); 1580 } 1581 1582 if (!ismatis) { 1583 PetscCall(PCHPDDMSetUpNeumannOverlap_Private(pc)); 1584 PetscCall(PetscOptionsGetBool(nullptr, pcpre, "-pc_hpddm_block_splitting", &block, nullptr)); 1585 PetscCall(PetscOptionsGetInt(nullptr, pcpre, "-pc_hpddm_harmonic_overlap", &overlap, nullptr)); 1586 PetscCall(PetscObjectTypeCompare((PetscObject)P, MATSCHURCOMPLEMENT, &flg)); 1587 if (data->is || (data->N > 1 && flg)) { 1588 if (block || overlap != -1) { 1589 PetscCall(ISDestroy(&data->is)); 1590 PetscCall(MatDestroy(&data->aux)); 1591 } else if (data->N > 1 && flg) { 1592 PCHPDDMSchurPreType type = PC_HPDDM_SCHUR_PRE_GENEO; 1593 1594 PetscCall(PetscOptionsGetEnum(nullptr, pcpre, "-pc_hpddm_schur_precondition", PCHPDDMSchurPreTypes, (PetscEnum *)&type, &flg)); 1595 if (type == PC_HPDDM_SCHUR_PRE_LEAST_SQUARES) { 1596 PetscCall(ISDestroy(&data->is)); /* destroy any previously user-set objects since they will be set automatically */ 1597 PetscCall(MatDestroy(&data->aux)); 1598 } else if (type == PC_HPDDM_SCHUR_PRE_GENEO) { 1599 PetscContainer container = nullptr; 1600 1601 PetscCall(PetscObjectQuery((PetscObject)pc, "_PCHPDDM_Schur", (PetscObject *)&container)); 1602 if (!container) { /* first call to PCSetUp() on the PC associated to the Schur complement */ 1603 PC_HPDDM *data_00; 1604 KSP ksp, inner_ksp; 1605 PC pc_00; 1606 char *prefix; 1607 PetscReal norm; 1608 1609 PetscCall(MatSchurComplementGetKSP(P, &ksp)); 1610 PetscCall(KSPGetPC(ksp, &pc_00)); 1611 PetscCall(PetscObjectTypeCompare((PetscObject)pc_00, PCHPDDM, &flg)); 1612 PetscCheck(flg, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-%spc_hpddm_schur_precondition %s and -%spc_type %s (!= %s)", pcpre ? pcpre : "", PCHPDDMSchurPreTypes[type], ((PetscObject)pc_00)->prefix ? ((PetscObject)pc_00)->prefix : "", 1613 ((PetscObject)pc_00)->type_name, PCHPDDM); 1614 data_00 = (PC_HPDDM *)pc_00->data; 1615 PetscCheck(data_00->N == 2, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-%spc_hpddm_schur_precondition %s and %" PetscInt_FMT " level%s instead of 2 for the A00 block -%s", pcpre ? pcpre : "", PCHPDDMSchurPreTypes[type], 1616 data_00->N, data_00->N > 1 ? "s" : "", ((PetscObject)pc_00)->prefix); 1617 PetscCall(PetscObjectTypeCompare((PetscObject)data_00->levels[0]->pc, PCASM, &flg)); 1618 PetscCheck(flg, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-%spc_hpddm_schur_precondition %s and -%spc_type %s (!= %s)", pcpre ? pcpre : "", PCHPDDMSchurPreTypes[type], ((PetscObject)data_00->levels[0]->pc)->prefix, 1619 ((PetscObject)data_00->levels[0]->pc)->type_name, PCASM); 1620 PetscCheck(data->Neumann == PETSC_BOOL3_TRUE, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-%spc_hpddm_schur_precondition %s and -%spc_hpddm_has_neumann != true", pcpre ? pcpre : "", PCHPDDMSchurPreTypes[type], pcpre ? pcpre : ""); 1621 if (PetscDefined(USE_DEBUG) || !data->is) { 1622 Mat A01, A10, B = nullptr, C = nullptr, *sub; 1623 1624 PetscCall(MatSchurComplementGetSubMatrices(P, &A, nullptr, &A01, &A10, nullptr)); 1625 PetscCall(PetscObjectTypeCompare((PetscObject)A10, MATTRANSPOSEVIRTUAL, &flg)); 1626 if (flg) { 1627 PetscCall(MatTransposeGetMat(A10, &C)); 1628 PetscCall(MatTranspose(C, MAT_INITIAL_MATRIX, &B)); 1629 } else { 1630 PetscCall(PetscObjectTypeCompare((PetscObject)A10, MATHERMITIANTRANSPOSEVIRTUAL, &flg)); 1631 if (flg) { 1632 PetscCall(MatHermitianTransposeGetMat(A10, &C)); 1633 PetscCall(MatHermitianTranspose(C, MAT_INITIAL_MATRIX, &B)); 1634 } 1635 } 1636 if (!B) { 1637 B = A10; 1638 PetscCall(PetscObjectReference((PetscObject)B)); 1639 } else if (!data->is) { 1640 PetscCall(PetscObjectTypeCompareAny((PetscObject)A01, &flg, MATTRANSPOSEVIRTUAL, MATHERMITIANTRANSPOSEVIRTUAL, "")); 1641 if (!flg) C = A01; 1642 } 1643 PetscCall(ISCreateStride(PETSC_COMM_SELF, B->rmap->N, 0, 1, &uis)); 1644 PetscCall(ISSetIdentity(uis)); 1645 if (!data->is) { 1646 if (C) PetscCall(PetscObjectReference((PetscObject)C)); 1647 else PetscCall(MatTranspose(B, MAT_INITIAL_MATRIX, &C)); 1648 PetscCall(ISDuplicate(data_00->is, is)); 1649 PetscCall(MatIncreaseOverlap(A, 1, is, 1)); 1650 PetscCall(MatSetOption(C, MAT_SUBMAT_SINGLEIS, PETSC_TRUE)); 1651 PetscCall(MatCreateSubMatrices(C, 1, is, &uis, MAT_INITIAL_MATRIX, &sub)); 1652 PetscCall(MatDestroy(&C)); 1653 PetscCall(MatTranspose(sub[0], MAT_INITIAL_MATRIX, &C)); 1654 PetscCall(MatDestroySubMatrices(1, &sub)); 1655 PetscCall(MatFindNonzeroRows(C, &data->is)); 1656 PetscCall(MatDestroy(&C)); 1657 PetscCall(ISDestroy(is)); 1658 } 1659 if (PetscDefined(USE_DEBUG)) { 1660 PetscCall(PCHPDDMCheckSymmetry_Private(pc, A01, A10)); 1661 PetscCall(MatCreateSubMatrices(B, 1, &uis, &data_00->is, MAT_INITIAL_MATRIX, &sub)); /* expensive check since all processes fetch all rows (but only some columns) of the constraint matrix */ 1662 PetscCall(ISDestroy(&uis)); 1663 PetscCall(ISDuplicate(data->is, &uis)); 1664 PetscCall(ISSort(uis)); 1665 PetscCall(ISComplement(uis, 0, B->rmap->N, is)); 1666 PetscCall(MatDuplicate(sub[0], MAT_COPY_VALUES, &C)); 1667 PetscCall(MatZeroRowsIS(C, is[0], 0.0, nullptr, nullptr)); 1668 PetscCall(ISDestroy(is)); 1669 PetscCall(MatMultEqual(sub[0], C, 20, &flg)); 1670 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "The image of A_10 (R_i^p)^T from the local primal (e.g., velocity) space to the full dual (e.g., pressure) space is not restricted to the local dual space: A_10 (R_i^p)^T != R_i^d (R_i^d)^T A_10 (R_i^p)^T"); /* cf. eq. (9) of https://hal.science/hal-02343808v6/document */ 1671 PetscCall(MatDestroy(&C)); 1672 PetscCall(MatDestroySubMatrices(1, &sub)); 1673 } 1674 PetscCall(ISDestroy(&uis)); 1675 PetscCall(MatDestroy(&B)); 1676 } 1677 if (data->aux) PetscCall(MatNorm(data->aux, NORM_FROBENIUS, &norm)); 1678 else norm = 0.0; 1679 PetscCall(MPIU_Allreduce(MPI_IN_PLACE, &norm, 1, MPIU_REAL, MPI_MAX, PetscObjectComm((PetscObject)P))); 1680 if (norm < PETSC_MACHINE_EPSILON * static_cast<PetscReal>(10.0)) { /* if A11 is near zero, e.g., Stokes equation, build a diagonal auxiliary (Neumann) Mat which is just a small diagonal weighted by the inverse of the multiplicity */ 1681 VecScatter scatter; 1682 Vec x; 1683 const PetscScalar *read; 1684 PetscScalar *write; 1685 1686 PetscCall(MatDestroy(&data->aux)); 1687 PetscCall(ISGetLocalSize(data->is, &n)); 1688 PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)P), n, PETSC_DECIDE, &x)); 1689 PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)P), n, PETSC_DECIDE, &v)); 1690 PetscCall(VecScatterCreate(x, data->is, v, nullptr, &scatter)); 1691 PetscCall(VecSet(v, 1.0)); 1692 PetscCall(VecSet(x, 1.0)); 1693 PetscCall(VecScatterBegin(scatter, v, x, ADD_VALUES, SCATTER_REVERSE)); 1694 PetscCall(VecScatterEnd(scatter, v, x, ADD_VALUES, SCATTER_REVERSE)); /* v has the multiplicity of all unknowns on the overlap */ 1695 PetscCall(VecScatterDestroy(&scatter)); 1696 PetscCall(VecDestroy(&v)); 1697 PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &v)); 1698 PetscCall(VecGetArrayRead(x, &read)); 1699 PetscCall(VecGetArrayWrite(v, &write)); 1700 PetscCallCXX(std::transform(read, read + n, write, [](const PetscScalar &m) { return PETSC_SMALL / (static_cast<PetscReal>(1000.0) * m); })); 1701 PetscCall(VecRestoreArrayRead(x, &read)); 1702 PetscCall(VecRestoreArrayWrite(v, &write)); 1703 PetscCall(VecDestroy(&x)); 1704 PetscCall(MatCreateDiagonal(v, &data->aux)); 1705 PetscCall(VecDestroy(&v)); 1706 } 1707 uis = data->is; 1708 uaux = data->aux; 1709 PetscCall(PetscObjectReference((PetscObject)uis)); 1710 PetscCall(PetscObjectReference((PetscObject)uaux)); 1711 PetscCall(PetscStrallocpy(pcpre, &prefix)); 1712 PetscCall(PCSetOptionsPrefix(pc, nullptr)); 1713 PetscCall(PCSetType(pc, PCKSP)); /* replace the PC associated to the Schur complement by PCKSP */ 1714 PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &inner_ksp)); /* new KSP that will be attached to the previously set PC */ 1715 pc->ops->setup = PCSetUp_KSP; 1716 PetscCall(PetscObjectGetTabLevel((PetscObject)pc, &n)); 1717 PetscCall(PetscObjectSetTabLevel((PetscObject)inner_ksp, n + 2)); 1718 PetscCall(KSPSetOperators(inner_ksp, pc->mat, pc->pmat)); 1719 PetscCall(KSPSetOptionsPrefix(inner_ksp, std::string(std::string(prefix) + "pc_hpddm_").c_str())); 1720 PetscCall(KSPSetSkipPCSetFromOptions(inner_ksp, PETSC_TRUE)); 1721 PetscCall(KSPSetFromOptions(inner_ksp)); 1722 PetscCall(KSPGetPC(inner_ksp, &inner)); 1723 PetscCall(PCSetOptionsPrefix(inner, nullptr)); 1724 PetscCall(PCSetType(inner, PCNONE)); /* no preconditioner since the action of M^-1 A or A M^-1 will be computed by the Amat */ 1725 PetscCall(PCKSPSetKSP(pc, inner_ksp)); 1726 PetscCall(PetscContainerCreate(PetscObjectComm((PetscObject)pc), &container)); 1727 PetscCall(PetscNew(&ctx)); /* context to pass data around for the inner-most PC, which will be a proper PCHPDDM */ 1728 PetscCall(PetscContainerSetPointer(container, ctx)); 1729 std::get<0>(*ctx)[0] = pc_00; /* for coarse correction on the primal (e.g., velocity) space */ 1730 PetscCall(PCCreate(PetscObjectComm((PetscObject)pc), &std::get<0>(*ctx)[1])); 1731 PetscCall(PCSetOptionsPrefix(std::get<0>(*ctx)[1], prefix)); 1732 PetscCall(PetscFree(prefix)); 1733 PetscCall(PCSetOperators(std::get<0>(*ctx)[1], pc->mat, pc->pmat)); 1734 PetscCall(PCSetType(std::get<0>(*ctx)[1], PCHPDDM)); 1735 PetscCall(PCHPDDMSetAuxiliaryMat(std::get<0>(*ctx)[1], uis, uaux, nullptr, nullptr)); /* transfer ownership of the auxiliary inputs from the inner (PCKSP) to the inner-most (PCHPDDM) PC */ 1736 PetscCall(PCSetFromOptions(std::get<0>(*ctx)[1])); 1737 PetscCall(PetscObjectDereference((PetscObject)uis)); 1738 PetscCall(PetscObjectDereference((PetscObject)uaux)); 1739 PetscCall(MatCreateShell(PetscObjectComm((PetscObject)pc), inner->mat->rmap->n, inner->mat->cmap->n, inner->mat->rmap->N, inner->mat->cmap->N, ctx, &S)); /* MatShell computing the action of M^-1 A or A M^-1 */ 1740 PetscCall(MatShellSetOperation(S, MATOP_MULT, (void (*)(void))MatMult_SchurCorrection)); 1741 PetscCall(MatShellSetOperation(S, MATOP_VIEW, (void (*)(void))MatView_SchurCorrection)); 1742 PetscCall(MatShellSetOperation(S, MATOP_DESTROY, (void (*)(void))MatDestroy_SchurCorrection)); 1743 PetscCall(KSPGetPCSide(inner_ksp, &(std::get<2>(*ctx)))); 1744 if (std::get<2>(*ctx) == PC_LEFT || std::get<2>(*ctx) == PC_SIDE_DEFAULT) { 1745 PetscCall(KSPSetPreSolve(inner_ksp, KSPPreSolve_SchurCorrection, ctx)); 1746 } else { /* no support for PC_SYMMETRIC */ 1747 PetscCheck(std::get<2>(*ctx) == PC_RIGHT, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "PCSide %s (!= %s or %s or %s)", PCSides[std::get<2>(*ctx)], PCSides[PC_SIDE_DEFAULT], PCSides[PC_LEFT], PCSides[PC_RIGHT]); 1748 } 1749 PetscCall(KSPSetPostSolve(inner_ksp, KSPPostSolve_SchurCorrection, ctx)); 1750 PetscCall(PetscObjectCompose((PetscObject)(std::get<0>(*ctx)[1]), "_PCHPDDM_Schur", (PetscObject)container)); 1751 PetscCall(PetscObjectDereference((PetscObject)container)); 1752 PetscCall(PCSetUp(std::get<0>(*ctx)[1])); 1753 PetscCall(KSPSetOperators(inner_ksp, S, S)); 1754 PetscCall(MatCreateVecs(std::get<1>(*ctx)[0], std::get<3>(*ctx), std::get<3>(*ctx) + 1)); 1755 PetscCall(VecDuplicate(std::get<3>(*ctx)[0], std::get<3>(*ctx) + 2)); 1756 PetscCall(PetscObjectDereference((PetscObject)inner_ksp)); 1757 PetscCall(PetscObjectDereference((PetscObject)S)); 1758 PetscFunctionReturn(PETSC_SUCCESS); 1759 } else { /* second call to PCSetUp() on the PC associated to the Schur complement, retrieve previously set context */ 1760 PetscCall(PetscContainerGetPointer(container, (void **)&ctx)); 1761 } 1762 } 1763 } 1764 } 1765 if (!data->is && data->N > 1) { 1766 char type[256] = {}; /* same size as in src/ksp/pc/interface/pcset.c */ 1767 PetscCall(PetscObjectTypeCompareAny((PetscObject)P, &flg, MATNORMAL, MATNORMALHERMITIAN, "")); 1768 if (flg || (A->rmap->N != A->cmap->N && P->rmap->N == P->cmap->N && P->rmap->N == A->cmap->N)) { 1769 Mat B; 1770 PetscCall(PCHPDDMSetAuxiliaryMatNormal_Private(pc, A, P, &B, pcpre)); 1771 if (data->correction == PC_HPDDM_COARSE_CORRECTION_DEFLATED) data->correction = PC_HPDDM_COARSE_CORRECTION_BALANCED; 1772 PetscCall(MatDestroy(&B)); 1773 } else { 1774 PetscCall(PetscObjectTypeCompare((PetscObject)P, MATSCHURCOMPLEMENT, &flg)); 1775 if (flg) { 1776 Mat A00, P00, A01, A10, A11, B, N; 1777 PCHPDDMSchurPreType type = PC_HPDDM_SCHUR_PRE_LEAST_SQUARES; 1778 1779 PetscCall(MatSchurComplementGetSubMatrices(P, &A00, &P00, &A01, &A10, &A11)); 1780 if (PetscDefined(USE_DEBUG)) PetscCall(PCHPDDMCheckSymmetry_Private(pc, A01, A10)); 1781 PetscCall(PetscOptionsGetEnum(nullptr, pcpre, "-pc_hpddm_schur_precondition", PCHPDDMSchurPreTypes, (PetscEnum *)&type, &flg)); 1782 if (type == PC_HPDDM_SCHUR_PRE_LEAST_SQUARES) { 1783 Vec diagonal = nullptr; 1784 const PetscScalar *array; 1785 MatSchurComplementAinvType type; 1786 1787 if (A11) { 1788 PetscCall(MatCreateVecs(A11, &diagonal, nullptr)); 1789 PetscCall(MatGetDiagonal(A11, diagonal)); 1790 } 1791 PetscCall(MatCreateVecs(P00, &v, nullptr)); 1792 PetscCall(MatSchurComplementGetAinvType(P, &type)); 1793 PetscCheck(type == MAT_SCHUR_COMPLEMENT_AINV_DIAG || type == MAT_SCHUR_COMPLEMENT_AINV_LUMP, PetscObjectComm((PetscObject)P), PETSC_ERR_SUP, "-%smat_schur_complement_ainv_type %s", ((PetscObject)P)->prefix ? ((PetscObject)P)->prefix : "", MatSchurComplementAinvTypes[type]); 1794 if (type == MAT_SCHUR_COMPLEMENT_AINV_LUMP) { 1795 PetscCall(MatGetRowSum(P00, v)); 1796 if (A00 == P00) PetscCall(PetscObjectReference((PetscObject)A00)); 1797 PetscCall(MatDestroy(&P00)); 1798 PetscCall(VecGetArrayRead(v, &array)); 1799 PetscCall(MatCreateAIJ(PetscObjectComm((PetscObject)A00), A00->rmap->n, A00->cmap->n, A00->rmap->N, A00->cmap->N, 1, nullptr, 0, nullptr, &P00)); 1800 PetscCall(MatSetOption(P00, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE)); 1801 for (n = A00->rmap->rstart; n < A00->rmap->rend; ++n) PetscCall(MatSetValue(P00, n, n, array[n - A00->rmap->rstart], INSERT_VALUES)); 1802 PetscCall(MatAssemblyBegin(P00, MAT_FINAL_ASSEMBLY)); 1803 PetscCall(MatAssemblyEnd(P00, MAT_FINAL_ASSEMBLY)); 1804 PetscCall(VecRestoreArrayRead(v, &array)); 1805 PetscCall(MatSchurComplementUpdateSubMatrices(P, A00, P00, A01, A10, A11)); /* replace P00 by diag(sum of each row of P00) */ 1806 PetscCall(MatDestroy(&P00)); 1807 } else PetscCall(MatGetDiagonal(P00, v)); 1808 PetscCall(VecReciprocal(v)); /* inv(diag(P00)) */ 1809 PetscCall(VecSqrtAbs(v)); /* sqrt(inv(diag(P00))) */ 1810 PetscCall(MatDuplicate(A01, MAT_COPY_VALUES, &B)); 1811 PetscCall(MatDiagonalScale(B, v, nullptr)); 1812 PetscCall(VecDestroy(&v)); 1813 PetscCall(MatCreateNormalHermitian(B, &N)); 1814 PetscCall(PCHPDDMSetAuxiliaryMatNormal_Private(pc, B, N, &P, pcpre, &diagonal)); 1815 PetscCall(PetscObjectTypeCompare((PetscObject)data->aux, MATSEQAIJ, &flg)); 1816 if (!flg) { 1817 PetscCall(MatDestroy(&P)); 1818 P = N; 1819 PetscCall(PetscObjectReference((PetscObject)P)); 1820 } else PetscCall(MatScale(P, -1.0)); 1821 if (diagonal) { 1822 PetscCall(MatDiagonalSet(P, diagonal, ADD_VALUES)); 1823 PetscCall(PCSetOperators(pc, P, P)); /* replace P by diag(P11) - A01^T inv(diag(P00)) A01 */ 1824 PetscCall(VecDestroy(&diagonal)); 1825 } else { 1826 PetscCall(MatScale(N, -1.0)); 1827 PetscCall(PCSetOperators(pc, N, P)); /* replace P by - A01^T inv(diag(P00)) A01 */ 1828 } 1829 PetscCall(MatDestroy(&N)); 1830 PetscCall(MatDestroy(&P)); 1831 PetscCall(MatDestroy(&B)); 1832 } else 1833 PetscCheck(type != PC_HPDDM_SCHUR_PRE_GENEO, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-%spc_hpddm_schur_precondition %s without a prior call to PCHPDDMSetAuxiliaryMat() on the A11 block%s%s", pcpre ? pcpre : "", PCHPDDMSchurPreTypes[type], pcpre ? " -" : "", pcpre ? pcpre : ""); 1834 PetscFunctionReturn(PETSC_SUCCESS); 1835 } else { 1836 PetscCall(PetscOptionsGetString(nullptr, pcpre, "-pc_hpddm_levels_1_st_pc_type", type, sizeof(type), nullptr)); 1837 PetscCall(PetscStrcmp(type, PCMAT, &algebraic)); 1838 PetscCheck(!algebraic || !block, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-pc_hpddm_levels_1_st_pc_type mat and -pc_hpddm_block_splitting"); 1839 if (overlap != -1) { 1840 PetscCheck(!block && !algebraic, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-pc_hpddm_%s and -pc_hpddm_harmonic_overlap", block ? "block_splitting" : "levels_1_st_pc_type mat"); 1841 PetscCheck(overlap >= 1, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_WRONG, "-pc_hpddm_harmonic_overlap %" PetscInt_FMT " < 1", overlap); 1842 } 1843 if (block || overlap != -1) algebraic = PETSC_TRUE; 1844 if (algebraic) { 1845 PetscCall(ISCreateStride(PETSC_COMM_SELF, P->rmap->n, P->rmap->rstart, 1, &data->is)); 1846 PetscCall(MatIncreaseOverlap(P, 1, &data->is, 1)); 1847 PetscCall(ISSort(data->is)); 1848 } else 1849 PetscCall(PetscInfo(pc, "Cannot assemble a fully-algebraic coarse operator with an assembled Pmat and -%spc_hpddm_levels_1_st_pc_type != mat and -%spc_hpddm_block_splitting != true and -%spc_hpddm_harmonic_overlap < 1\n", pcpre ? pcpre : "", pcpre ? pcpre : "", pcpre ? pcpre : "")); 1850 } 1851 } 1852 } 1853 } 1854 #if PetscDefined(USE_DEBUG) 1855 if (data->is) PetscCall(ISDuplicate(data->is, &dis)); 1856 if (data->aux) PetscCall(MatDuplicate(data->aux, MAT_COPY_VALUES, &daux)); 1857 #endif 1858 if (data->is || (ismatis && data->N > 1)) { 1859 if (ismatis) { 1860 std::initializer_list<std::string> list = {MATSEQBAIJ, MATSEQSBAIJ}; 1861 PetscCall(MatISGetLocalMat(P, &N)); 1862 std::initializer_list<std::string>::const_iterator it = std::find(list.begin(), list.end(), ((PetscObject)N)->type_name); 1863 PetscCall(MatISRestoreLocalMat(P, &N)); 1864 switch (std::distance(list.begin(), it)) { 1865 case 0: 1866 PetscCall(MatConvert(P, MATMPIBAIJ, MAT_INITIAL_MATRIX, &C)); 1867 break; 1868 case 1: 1869 /* MatCreateSubMatrices() does not work with MATSBAIJ and unsorted ISes, so convert to MPIBAIJ */ 1870 PetscCall(MatConvert(P, MATMPIBAIJ, MAT_INITIAL_MATRIX, &C)); 1871 PetscCall(MatSetOption(C, MAT_SYMMETRIC, PETSC_TRUE)); 1872 break; 1873 default: 1874 PetscCall(MatConvert(P, MATMPIAIJ, MAT_INITIAL_MATRIX, &C)); 1875 } 1876 PetscCall(MatISGetLocalToGlobalMapping(P, &l2g, nullptr)); 1877 PetscCall(PetscObjectReference((PetscObject)P)); 1878 PetscCall(KSPSetOperators(data->levels[0]->ksp, A, C)); 1879 std::swap(C, P); 1880 PetscCall(ISLocalToGlobalMappingGetSize(l2g, &n)); 1881 PetscCall(ISCreateStride(PETSC_COMM_SELF, n, 0, 1, &loc)); 1882 PetscCall(ISLocalToGlobalMappingApplyIS(l2g, loc, &is[0])); 1883 PetscCall(ISDestroy(&loc)); 1884 /* the auxiliary Mat is _not_ the local Neumann matrix */ 1885 /* it is the local Neumann matrix augmented (with zeros) through MatIncreaseOverlap() */ 1886 data->Neumann = PETSC_BOOL3_FALSE; 1887 structure = SAME_NONZERO_PATTERN; 1888 } else { 1889 is[0] = data->is; 1890 if (algebraic || ctx) subdomains = PETSC_TRUE; 1891 PetscCall(PetscOptionsGetBool(nullptr, pcpre, "-pc_hpddm_define_subdomains", &subdomains, nullptr)); 1892 if (ctx) PetscCheck(subdomains, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-%spc_hpddm_schur_precondition geneo and -%spc_hpddm_define_subdomains false", pcpre, pcpre); 1893 if (PetscBool3ToBool(data->Neumann)) { 1894 PetscCheck(!block, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-pc_hpddm_block_splitting and -pc_hpddm_has_neumann"); 1895 PetscCheck(overlap == -1, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-pc_hpddm_harmonic_overlap %" PetscInt_FMT " and -pc_hpddm_has_neumann", overlap); 1896 PetscCheck(!algebraic, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-pc_hpddm_levels_1_st_pc_type mat and -pc_hpddm_has_neumann"); 1897 } 1898 if (PetscBool3ToBool(data->Neumann) || block) structure = SAME_NONZERO_PATTERN; 1899 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)data->is), P->rmap->n, P->rmap->rstart, 1, &loc)); 1900 } 1901 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_levels_1_", pcpre ? pcpre : "")); 1902 PetscCall(PetscOptionsGetEnum(nullptr, prefix, "-st_matstructure", MatStructures, (PetscEnum *)&structure, &flg)); /* if not user-provided, force its value when possible */ 1903 if (!flg && structure == SAME_NONZERO_PATTERN) { /* cannot call STSetMatStructure() yet, insert the appropriate option in the database, parsed by STSetFromOptions() */ 1904 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-%spc_hpddm_levels_1_st_matstructure", pcpre ? pcpre : "")); 1905 PetscCall(PetscOptionsSetValue(nullptr, prefix, MatStructures[structure])); 1906 } 1907 flg = PETSC_FALSE; 1908 if (data->share) { 1909 data->share = PETSC_FALSE; /* will be reset to PETSC_TRUE if none of the conditions below are true */ 1910 if (!subdomains) PetscCall(PetscInfo(pc, "Cannot share subdomain KSP between SLEPc and PETSc since -%spc_hpddm_define_subdomains is not true\n", pcpre ? pcpre : "")); 1911 else if (data->deflation) PetscCall(PetscInfo(pc, "Nothing to share since PCHPDDMSetDeflationMat() has been called\n")); 1912 else if (ismatis) PetscCall(PetscInfo(pc, "Cannot share subdomain KSP between SLEPc and PETSc with a Pmat of type MATIS\n")); 1913 else if (!algebraic && structure != SAME_NONZERO_PATTERN) 1914 PetscCall(PetscInfo(pc, "Cannot share subdomain KSP between SLEPc and PETSc since -%spc_hpddm_levels_1_st_matstructure %s (!= %s)\n", pcpre ? pcpre : "", MatStructures[structure], MatStructures[SAME_NONZERO_PATTERN])); 1915 else data->share = PETSC_TRUE; 1916 } 1917 if (!ismatis) { 1918 if (data->share || (!PetscBool3ToBool(data->Neumann) && subdomains)) PetscCall(ISDuplicate(is[0], &unsorted)); 1919 else unsorted = is[0]; 1920 } 1921 if (data->N > 1 && (data->aux || ismatis || algebraic)) { 1922 PetscCheck(loadedSym, PETSC_COMM_SELF, PETSC_ERR_PLIB, "HPDDM library not loaded, cannot use more than one level"); 1923 PetscCall(MatSetOption(P, MAT_SUBMAT_SINGLEIS, PETSC_TRUE)); 1924 if (ismatis) { 1925 /* needed by HPDDM (currently) so that the partition of unity is 0 on subdomain interfaces */ 1926 PetscCall(MatIncreaseOverlap(P, 1, is, 1)); 1927 PetscCall(ISDestroy(&data->is)); 1928 data->is = is[0]; 1929 } else { 1930 if (PetscDefined(USE_DEBUG)) { 1931 PetscBool equal; 1932 IS intersect; 1933 1934 PetscCall(ISIntersect(data->is, loc, &intersect)); 1935 PetscCall(ISEqualUnsorted(loc, intersect, &equal)); 1936 PetscCall(ISDestroy(&intersect)); 1937 PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "IS of the auxiliary Mat does not include all local rows of A"); 1938 } 1939 if (overlap == -1) PetscCall(PetscObjectComposeFunction((PetscObject)pc->pmat, "PCHPDDMAlgebraicAuxiliaryMat_Private_C", PCHPDDMAlgebraicAuxiliaryMat_Private)); 1940 if (!PetscBool3ToBool(data->Neumann) && (!algebraic || overlap != -1)) { 1941 PetscCall(PetscObjectTypeCompare((PetscObject)P, MATMPISBAIJ, &flg)); 1942 if (flg) { 1943 /* maybe better to ISSort(is[0]), MatCreateSubMatrices(), and then MatPermute() */ 1944 /* but there is no MatPermute_SeqSBAIJ(), so as before, just use MATMPIBAIJ */ 1945 PetscCall(MatConvert(P, MATMPIBAIJ, MAT_INITIAL_MATRIX, &uaux)); 1946 flg = PETSC_FALSE; 1947 } 1948 } 1949 } 1950 if (algebraic && overlap == -1) { 1951 PetscUseMethod(pc->pmat, "PCHPDDMAlgebraicAuxiliaryMat_Private_C", (Mat, IS *, Mat *[], PetscBool), (P, is, &sub, block)); 1952 if (block) { 1953 PetscCall(PetscObjectQuery((PetscObject)sub[0], "_PCHPDDM_Neumann_Mat", (PetscObject *)&data->aux)); 1954 PetscCall(PetscObjectCompose((PetscObject)sub[0], "_PCHPDDM_Neumann_Mat", nullptr)); 1955 } 1956 } else if (!uaux || overlap != -1) { 1957 if (!ctx) { 1958 if (PetscBool3ToBool(data->Neumann)) sub = &data->aux; 1959 else { 1960 if (overlap != -1) { 1961 Harmonic h; 1962 Mat A0, *a; /* with an SVD: [ A_00 A_01 ] */ 1963 IS ov[2], rows, cols, stride; /* [ A_10 A_11 A_12 ] */ 1964 const PetscInt *i[2], bs = PetscAbs(P->cmap->bs); /* with a GEVP: [ A_00 A_01 ] */ 1965 PetscInt n[2]; /* [ A_10 A_11 A_12 ] */ 1966 std::vector<PetscInt> v[2]; /* [ A_21 A_22 ] */ 1967 PetscBool flg; 1968 1969 PetscCall(ISDuplicate(data->is, ov)); 1970 if (overlap > 1) PetscCall(MatIncreaseOverlap(P, 1, ov, overlap - 1)); 1971 PetscCall(ISDuplicate(ov[0], ov + 1)); 1972 PetscCall(MatIncreaseOverlap(P, 1, ov + 1, 1)); 1973 PetscCall(PetscNew(&h)); 1974 h->ksp = nullptr; 1975 PetscCall(PetscCalloc1(2, &h->A)); 1976 PetscCall(PetscOptionsHasName(nullptr, prefix, "-svd_nsv", &flg)); 1977 if (!flg) PetscCall(PetscOptionsHasName(nullptr, prefix, "-svd_relative_threshold", &flg)); 1978 PetscCall(ISSort(ov[0])); 1979 if (!flg) PetscCall(ISSort(ov[1])); 1980 PetscCall(PetscMalloc1(!flg ? 5 : 3, &h->is)); 1981 PetscCall(MatCreateSubMatrices(uaux ? uaux : P, 1, ov + !flg, ov + 1, MAT_INITIAL_MATRIX, &a)); /* submatrix from above, either square (!flg) or rectangular (flg) */ 1982 for (PetscInt j = 0; j < 2; ++j) { 1983 PetscCall(ISGetIndices(ov[j], i + j)); 1984 PetscCall(ISGetLocalSize(ov[j], n + j)); 1985 } 1986 v[1].reserve((n[1] - n[0]) / bs); 1987 for (PetscInt j = 0; j < n[1]; j += bs) { /* indices of the (2,2) block */ 1988 PetscInt location; 1989 PetscCall(ISLocate(ov[0], i[1][j], &location)); 1990 if (location < 0) v[1].emplace_back(j / bs); 1991 } 1992 if (!flg) { 1993 h->A[1] = a[0]; 1994 PetscCall(PetscObjectReference((PetscObject)h->A[1])); 1995 v[0].reserve((n[0] - P->rmap->n) / bs); 1996 for (PetscInt j = 0; j < n[1]; j += bs) { /* row indices of the (1,2) block */ 1997 PetscInt location; 1998 PetscCall(ISLocate(loc, i[1][j], &location)); 1999 if (location < 0) { 2000 PetscCall(ISLocate(ov[0], i[1][j], &location)); 2001 if (location >= 0) v[0].emplace_back(j / bs); 2002 } 2003 } 2004 PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, v[0].size(), v[0].data(), PETSC_USE_POINTER, &rows)); 2005 PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, v[1].size(), v[1].data(), PETSC_COPY_VALUES, h->is + 4)); 2006 PetscCall(MatCreateSubMatrix(a[0], rows, h->is[4], MAT_INITIAL_MATRIX, h->A)); /* A_12 submatrix from above */ 2007 PetscCall(ISDestroy(&rows)); 2008 if (uaux) PetscCall(MatConvert(a[0], MATSEQSBAIJ, MAT_INPLACE_MATRIX, a)); /* initial Pmat was MATSBAIJ, convert back to the same format since the rectangular A_12 submatrix has been created */ 2009 PetscCall(ISEmbed(ov[0], ov[1], PETSC_TRUE, &rows)); 2010 PetscCall(MatCreateSubMatrix(a[0], rows, cols = rows, MAT_INITIAL_MATRIX, &A0)); /* [ A_00 A_01 ; A_10 A_11 ] submatrix from above */ 2011 PetscCall(ISDestroy(&rows)); 2012 v[0].clear(); 2013 PetscCall(ISEmbed(loc, ov[1], PETSC_TRUE, h->is + 3)); 2014 PetscCall(ISEmbed(data->is, ov[1], PETSC_TRUE, h->is + 2)); 2015 } 2016 v[0].reserve((n[0] - P->rmap->n) / bs); 2017 for (PetscInt j = 0; j < n[0]; j += bs) { 2018 PetscInt location; 2019 PetscCall(ISLocate(loc, i[0][j], &location)); 2020 if (location < 0) v[0].emplace_back(j / bs); 2021 } 2022 PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, v[0].size(), v[0].data(), PETSC_USE_POINTER, &rows)); 2023 for (PetscInt j = 0; j < 2; ++j) PetscCall(ISRestoreIndices(ov[j], i + j)); 2024 if (flg) { 2025 IS is; 2026 PetscCall(ISCreateStride(PETSC_COMM_SELF, a[0]->rmap->n, 0, 1, &is)); 2027 PetscCall(ISEmbed(ov[0], ov[1], PETSC_TRUE, &cols)); 2028 PetscCall(MatCreateSubMatrix(a[0], is, cols, MAT_INITIAL_MATRIX, &A0)); /* [ A_00 A_01 ; A_10 A_11 ] submatrix from above */ 2029 PetscCall(ISDestroy(&cols)); 2030 PetscCall(ISDestroy(&is)); 2031 if (uaux) PetscCall(MatConvert(A0, MATSEQSBAIJ, MAT_INPLACE_MATRIX, &A0)); /* initial Pmat was MATSBAIJ, convert back to the same format since this submatrix is square */ 2032 PetscCall(ISEmbed(loc, data->is, PETSC_TRUE, h->is + 2)); 2033 PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, v[1].size(), v[1].data(), PETSC_USE_POINTER, &cols)); 2034 PetscCall(MatCreateSubMatrix(a[0], rows, cols, MAT_INITIAL_MATRIX, h->A)); /* A_12 submatrix from above */ 2035 PetscCall(ISDestroy(&cols)); 2036 } 2037 PetscCall(ISCreateStride(PETSC_COMM_SELF, A0->rmap->n, 0, 1, &stride)); 2038 PetscCall(ISEmbed(rows, stride, PETSC_TRUE, h->is)); 2039 PetscCall(ISDestroy(&stride)); 2040 PetscCall(ISDestroy(&rows)); 2041 PetscCall(ISEmbed(loc, ov[0], PETSC_TRUE, h->is + 1)); 2042 if (subdomains) { 2043 if (!data->levels[0]->pc) { 2044 PetscCall(PCCreate(PetscObjectComm((PetscObject)pc), &data->levels[0]->pc)); 2045 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_levels_1_", pcpre ? pcpre : "")); 2046 PetscCall(PCSetOptionsPrefix(data->levels[0]->pc, prefix)); 2047 PetscCall(PCSetOperators(data->levels[0]->pc, A, P)); 2048 } 2049 PetscCall(PCSetType(data->levels[0]->pc, PCASM)); 2050 if (!data->levels[0]->pc->setupcalled) PetscCall(PCASMSetLocalSubdomains(data->levels[0]->pc, 1, ov + !flg, &loc)); 2051 PetscCall(PCHPDDMCommunicationAvoidingPCASM_Private(data->levels[0]->pc, flg ? A0 : a[0], PETSC_TRUE)); 2052 if (!flg) ++overlap; 2053 if (data->share) { 2054 PetscInt n = -1; 2055 PetscTryMethod(data->levels[0]->pc, "PCASMGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (data->levels[0]->pc, &n, nullptr, &ksp)); 2056 PetscCheck(n == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of subdomain solver %" PetscInt_FMT " != 1", n); 2057 if (flg) { 2058 h->ksp = ksp[0]; 2059 PetscCall(PetscObjectReference((PetscObject)h->ksp)); 2060 } 2061 } 2062 } 2063 if (!h->ksp) { 2064 PetscBool share = data->share; 2065 PetscCall(KSPCreate(PETSC_COMM_SELF, &h->ksp)); 2066 PetscCall(KSPSetType(h->ksp, KSPPREONLY)); 2067 PetscCall(KSPSetOperators(h->ksp, A0, A0)); 2068 do { 2069 if (!data->share) { 2070 share = PETSC_FALSE; 2071 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_levels_1_%s", pcpre ? pcpre : "", flg ? "svd_" : "eps_")); 2072 PetscCall(KSPSetOptionsPrefix(h->ksp, prefix)); 2073 PetscCall(KSPSetFromOptions(h->ksp)); 2074 } else { 2075 MatSolverType type; 2076 PetscCall(KSPGetPC(ksp[0], &pc)); 2077 PetscCall(PetscObjectTypeCompareAny((PetscObject)pc, &data->share, PCLU, PCCHOLESKY, "")); 2078 if (data->share) { 2079 PetscCall(PCFactorGetMatSolverType(pc, &type)); 2080 if (!type) { 2081 if (PetscDefined(HAVE_MUMPS)) PetscCall(PCFactorSetMatSolverType(pc, MATSOLVERMUMPS)); 2082 else if (PetscDefined(HAVE_MKL_PARDISO)) PetscCall(PCFactorSetMatSolverType(pc, MATSOLVERMKL_PARDISO)); 2083 else data->share = PETSC_FALSE; 2084 if (data->share) PetscCall(PCSetFromOptions(pc)); 2085 } else { 2086 PetscCall(PetscStrcmp(type, MATSOLVERMUMPS, &data->share)); 2087 if (!data->share) PetscCall(PetscStrcmp(type, MATSOLVERMKL_PARDISO, &data->share)); 2088 } 2089 if (data->share) { 2090 std::tuple<KSP, IS, Vec[2]> *p; 2091 PetscCall(PCFactorGetMatrix(pc, &A)); 2092 PetscCall(MatFactorSetSchurIS(A, h->is[4])); 2093 PetscCall(KSPSetUp(ksp[0])); 2094 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_levels_1_eps_shell_", pcpre ? pcpre : "")); 2095 PetscCall(KSPSetOptionsPrefix(h->ksp, prefix)); 2096 PetscCall(KSPSetFromOptions(h->ksp)); 2097 PetscCall(KSPGetPC(h->ksp, &pc)); 2098 PetscCall(PCSetType(pc, PCSHELL)); 2099 PetscCall(PetscNew(&p)); 2100 std::get<0>(*p) = ksp[0]; 2101 PetscCall(ISEmbed(ov[0], ov[1], PETSC_TRUE, &std::get<1>(*p))); 2102 PetscCall(MatCreateVecs(A, std::get<2>(*p), std::get<2>(*p) + 1)); 2103 PetscCall(PCShellSetContext(pc, p)); 2104 PetscCall(PCShellSetApply(pc, PCApply_Schur)); 2105 PetscCall(PCShellSetApplyTranspose(pc, PCApply_Schur<Vec, true>)); 2106 PetscCall(PCShellSetMatApply(pc, PCApply_Schur<Mat>)); 2107 PetscCall(PCShellSetDestroy(pc, PCDestroy_Schur)); 2108 } 2109 } 2110 if (!data->share) PetscCall(PetscInfo(pc, "Cannot share subdomain KSP between SLEPc and PETSc since neither MUMPS nor MKL PARDISO is used\n")); 2111 } 2112 } while (!share != !data->share); /* if data->share is initially PETSC_TRUE, but then reset to PETSC_FALSE, then go back to the beginning of the do loop */ 2113 } 2114 PetscCall(ISDestroy(ov)); 2115 PetscCall(ISDestroy(ov + 1)); 2116 if (overlap == 1 && subdomains && flg) { 2117 *subA = A0; 2118 sub = subA; 2119 if (uaux) PetscCall(MatDestroy(&uaux)); 2120 } else PetscCall(MatDestroy(&A0)); 2121 PetscCall(MatCreateShell(PETSC_COMM_SELF, P->rmap->n, n[1] - n[0], P->rmap->n, n[1] - n[0], h, &data->aux)); 2122 PetscCall(MatCreateVecs(h->ksp->pc->pmat, &h->v, nullptr)); 2123 PetscCall(MatShellSetOperation(data->aux, MATOP_MULT, (void (*)(void))MatMult_Harmonic)); 2124 PetscCall(MatShellSetOperation(data->aux, MATOP_MULT_TRANSPOSE, (void (*)(void))MatMultTranspose_Harmonic)); 2125 PetscCall(MatShellSetMatProductOperation(data->aux, MATPRODUCT_AB, nullptr, MatProduct_AB_Harmonic, nullptr, MATDENSE, MATDENSE)); 2126 PetscCall(MatShellSetOperation(data->aux, MATOP_DESTROY, (void (*)(void))MatDestroy_Harmonic)); 2127 PetscCall(MatDestroySubMatrices(1, &a)); 2128 } 2129 if (overlap != 1 || !subdomains) PetscCall(MatCreateSubMatrices(uaux ? uaux : P, 1, is, is, MAT_INITIAL_MATRIX, &sub)); 2130 if (uaux) { 2131 PetscCall(MatDestroy(&uaux)); 2132 PetscCall(MatConvert(sub[0], MATSEQSBAIJ, MAT_INPLACE_MATRIX, sub)); 2133 } 2134 } 2135 } 2136 } else { 2137 PetscCall(MatCreateSubMatrices(uaux, 1, is, is, MAT_INITIAL_MATRIX, &sub)); 2138 PetscCall(MatDestroy(&uaux)); 2139 PetscCall(MatConvert(sub[0], MATSEQSBAIJ, MAT_INPLACE_MATRIX, sub)); 2140 } 2141 /* Vec holding the partition of unity */ 2142 if (!data->levels[0]->D) { 2143 PetscCall(ISGetLocalSize(data->is, &n)); 2144 PetscCall(VecCreateMPI(PETSC_COMM_SELF, n, PETSC_DETERMINE, &data->levels[0]->D)); 2145 } 2146 if (data->share && overlap == -1) { 2147 Mat D; 2148 IS perm = nullptr; 2149 PetscInt size = -1; 2150 if (!data->levels[0]->pc) { 2151 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_levels_1_", pcpre ? pcpre : "")); 2152 PetscCall(PCCreate(PetscObjectComm((PetscObject)pc), &data->levels[0]->pc)); 2153 PetscCall(PCSetOptionsPrefix(data->levels[0]->pc, prefix)); 2154 PetscCall(PCSetOperators(data->levels[0]->pc, A, P)); 2155 } 2156 PetscCall(PCSetType(data->levels[0]->pc, PCASM)); 2157 if (!ctx) { 2158 if (!data->levels[0]->pc->setupcalled) { 2159 IS sorted; /* PCASM will sort the input IS, duplicate it to return an unmodified (PCHPDDM) input IS */ 2160 PetscCall(ISDuplicate(is[0], &sorted)); 2161 PetscCall(PCASMSetLocalSubdomains(data->levels[0]->pc, 1, &sorted, &loc)); 2162 PetscCall(PetscObjectDereference((PetscObject)sorted)); 2163 } 2164 PetscCall(PCSetFromOptions(data->levels[0]->pc)); 2165 if (block) { 2166 PetscCall(PCHPDDMPermute_Private(unsorted, data->is, &uis, sub[0], &C, &perm)); 2167 PetscCall(PCHPDDMCommunicationAvoidingPCASM_Private(data->levels[0]->pc, C, algebraic)); 2168 } else PetscCall(PCSetUp(data->levels[0]->pc)); 2169 PetscTryMethod(data->levels[0]->pc, "PCASMGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (data->levels[0]->pc, &size, nullptr, &ksp)); 2170 if (size != 1) { 2171 data->share = PETSC_FALSE; 2172 PetscCheck(size == -1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of subdomain solver %" PetscInt_FMT " != 1", size); 2173 PetscCall(PetscInfo(pc, "Cannot share subdomain KSP between SLEPc and PETSc since PCASMGetSubKSP() not found in fine-level PC\n")); 2174 PetscCall(ISDestroy(&unsorted)); 2175 unsorted = is[0]; 2176 } else { 2177 if (!block && !ctx) PetscCall(PCHPDDMPermute_Private(unsorted, data->is, &uis, PetscBool3ToBool(data->Neumann) ? sub[0] : data->aux, &C, &perm)); 2178 if (!PetscBool3ToBool(data->Neumann) && !block) { 2179 PetscCall(MatPermute(sub[0], perm, perm, &D)); /* permute since PCASM will call ISSort() */ 2180 PetscCall(MatHeaderReplace(sub[0], &D)); 2181 } 2182 if (data->B) { /* see PCHPDDMSetRHSMat() */ 2183 PetscCall(MatPermute(data->B, perm, perm, &D)); 2184 PetscCall(MatHeaderReplace(data->B, &D)); 2185 } 2186 PetscCall(ISDestroy(&perm)); 2187 const char *matpre; 2188 PetscBool cmp[4]; 2189 PetscCall(KSPGetOperators(ksp[0], subA, subA + 1)); 2190 PetscCall(MatDuplicate(subA[1], MAT_SHARE_NONZERO_PATTERN, &D)); 2191 PetscCall(MatGetOptionsPrefix(subA[1], &matpre)); 2192 PetscCall(MatSetOptionsPrefix(D, matpre)); 2193 PetscCall(PetscObjectTypeCompare((PetscObject)D, MATNORMAL, cmp)); 2194 PetscCall(PetscObjectTypeCompare((PetscObject)C, MATNORMAL, cmp + 1)); 2195 if (!cmp[0]) PetscCall(PetscObjectTypeCompare((PetscObject)D, MATNORMALHERMITIAN, cmp + 2)); 2196 else cmp[2] = PETSC_FALSE; 2197 if (!cmp[1]) PetscCall(PetscObjectTypeCompare((PetscObject)C, MATNORMALHERMITIAN, cmp + 3)); 2198 else cmp[3] = PETSC_FALSE; 2199 PetscCheck(cmp[0] == cmp[1] && cmp[2] == cmp[3], PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "-%spc_hpddm_levels_1_pc_asm_sub_mat_type %s and auxiliary Mat of type %s", pcpre ? pcpre : "", ((PetscObject)D)->type_name, ((PetscObject)C)->type_name); 2200 if (!cmp[0] && !cmp[2]) { 2201 if (!block) PetscCall(MatAXPY(D, 1.0, C, SUBSET_NONZERO_PATTERN)); 2202 else { 2203 PetscCall(MatMissingDiagonal(D, cmp, nullptr)); 2204 if (cmp[0]) structure = DIFFERENT_NONZERO_PATTERN; /* data->aux has no missing diagonal entry */ 2205 PetscCall(MatAXPY(D, 1.0, data->aux, structure)); 2206 } 2207 } else { 2208 Mat mat[2]; 2209 if (cmp[0]) { 2210 PetscCall(MatNormalGetMat(D, mat)); 2211 PetscCall(MatNormalGetMat(C, mat + 1)); 2212 } else { 2213 PetscCall(MatNormalHermitianGetMat(D, mat)); 2214 PetscCall(MatNormalHermitianGetMat(C, mat + 1)); 2215 } 2216 PetscCall(MatAXPY(mat[0], 1.0, mat[1], SUBSET_NONZERO_PATTERN)); 2217 } 2218 PetscCall(MatPropagateSymmetryOptions(C, D)); 2219 PetscCall(MatDestroy(&C)); 2220 C = D; 2221 /* swap pointers so that variables stay consistent throughout PCSetUp() */ 2222 std::swap(C, data->aux); 2223 std::swap(uis, data->is); 2224 swap = PETSC_TRUE; 2225 } 2226 } 2227 } 2228 if (ctx) { 2229 PC_HPDDM *data_00 = (PC_HPDDM *)std::get<0>(*ctx)[0]->data; 2230 PC s; 2231 Mat A00, P00, A01 = nullptr, A10, A11, N, b[4]; 2232 IS sorted, is[2]; 2233 MatSolverType type; 2234 std::pair<PC, Vec[2]> *p; 2235 2236 PetscCall(PCHPDDMPermute_Private(unsorted, data->is, &uis, data->aux, &C, nullptr)); /* permute since PCASM works with a sorted IS */ 2237 std::swap(C, data->aux); 2238 std::swap(uis, data->is); 2239 swap = PETSC_TRUE; 2240 PetscCall(PCASMSetType(data->levels[0]->pc, PC_ASM_NONE)); /* "Neumann--Neumann" preconditioning with overlap and a Boolean partition of unity */ 2241 PetscCall(PCASMSetLocalSubdomains(data->levels[0]->pc, 1, &data->is, &loc)); 2242 PetscCall(PCSetFromOptions(data->levels[0]->pc)); /* action of eq. (15) of https://hal.science/hal-02343808v6/document (with a sign flip) */ 2243 PetscCall(MatSchurComplementGetSubMatrices(P, &A00, &P00, std::get<1>(*ctx), &A10, &A11)); 2244 std::get<1>(*ctx)[1] = A10; 2245 PetscCall(PetscObjectTypeCompare((PetscObject)A10, MATTRANSPOSEVIRTUAL, &flg)); 2246 if (flg) PetscCall(MatTransposeGetMat(A10, &A01)); 2247 else { 2248 PetscBool flg; 2249 2250 PetscCall(PetscObjectTypeCompare((PetscObject)A10, MATHERMITIANTRANSPOSEVIRTUAL, &flg)); 2251 if (flg) PetscCall(MatHermitianTransposeGetMat(A10, &A01)); 2252 } 2253 PetscCall(ISDuplicate(data_00->is, &sorted)); /* during setup of the PC associated to the A00 block, this IS has already been sorted, but it's put back to its original state at the end of PCSetUp_HPDDM(), which may be unsorted */ 2254 PetscCall(ISSort(sorted)); /* this is to avoid changing users inputs, but it requires a new call to ISSort() here */ 2255 if (!A01) { 2256 PetscCall(MatSetOption(A10, MAT_SUBMAT_SINGLEIS, PETSC_TRUE)); 2257 PetscCall(MatCreateSubMatrices(A10, 1, &data->is, &sorted, MAT_INITIAL_MATRIX, &sub)); 2258 b[2] = sub[0]; 2259 PetscCall(PetscObjectReference((PetscObject)sub[0])); 2260 PetscCall(MatDestroySubMatrices(1, &sub)); 2261 PetscCall(PetscObjectTypeCompare((PetscObject)std::get<1>(*ctx)[0], MATTRANSPOSEVIRTUAL, &flg)); 2262 A10 = nullptr; 2263 if (flg) PetscCall(MatTransposeGetMat(std::get<1>(*ctx)[0], &A10)); 2264 else { 2265 PetscBool flg; 2266 2267 PetscCall(PetscObjectTypeCompare((PetscObject)(PetscObject)std::get<1>(*ctx)[0], MATHERMITIANTRANSPOSEVIRTUAL, &flg)); 2268 if (flg) PetscCall(MatHermitianTransposeGetMat(std::get<1>(*ctx)[0], &A10)); 2269 } 2270 if (!A10) { 2271 PetscCall(MatCreateSubMatrices(std::get<1>(*ctx)[0], 1, &sorted, &data->is, MAT_INITIAL_MATRIX, &sub)); 2272 b[1] = sub[0]; 2273 PetscCall(PetscObjectReference((PetscObject)sub[0])); 2274 } 2275 } else { 2276 PetscCall(MatSetOption(A01, MAT_SUBMAT_SINGLEIS, PETSC_TRUE)); 2277 PetscCall(MatCreateSubMatrices(A01, 1, &sorted, &data->is, MAT_INITIAL_MATRIX, &sub)); 2278 if (flg) PetscCall(MatTranspose(*sub, MAT_INITIAL_MATRIX, b + 2)); 2279 else PetscCall(MatHermitianTranspose(*sub, MAT_INITIAL_MATRIX, b + 2)); 2280 } 2281 PetscCall(MatDestroySubMatrices(1, &sub)); 2282 PetscCall(ISDestroy(&sorted)); 2283 n = -1; 2284 PetscTryMethod(data_00->levels[0]->pc, "PCASMGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (data_00->levels[0]->pc, &n, nullptr, &ksp)); 2285 PetscCheck(n == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of subdomain solver %" PetscInt_FMT " != 1", n); 2286 PetscCall(KSPGetOperators(ksp[0], subA, subA + 1)); 2287 PetscCall(ISGetLocalSize(data_00->is, &n)); 2288 PetscCheck(n == subA[0]->rmap->n && n == subA[0]->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "-%spc_hpddm_schur_precondition geneo and -%spc_hpddm_define_subdomains false", pcpre ? pcpre : "", ((PetscObject)pc)->prefix); 2289 if (A01 || A10) { 2290 if (flg) PetscCall(MatTranspose(b[2], MAT_INITIAL_MATRIX, b + 1)); 2291 else PetscCall(MatHermitianTranspose(b[2], MAT_INITIAL_MATRIX, b + 1)); 2292 } 2293 PetscCall(MatCreateSchurComplement(subA[0], subA[1], b[1], b[2], data->aux, &S)); 2294 PetscCall(MatSchurComplementSetKSP(S, ksp[0])); 2295 PetscCall(PCHPDDMCommunicationAvoidingPCASM_Private(data->levels[0]->pc, S, PETSC_TRUE)); /* the subdomain Mat is already known and the input IS of PCASMSetLocalSubdomains() is already sorted */ 2296 PetscTryMethod(data->levels[0]->pc, "PCASMGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (data->levels[0]->pc, &n, nullptr, &ksp)); 2297 PetscCheck(n == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of subdomain solver %" PetscInt_FMT " != 1", n); 2298 PetscCall(KSPGetPC(ksp[0], &inner)); 2299 PetscCall(PCSetType(inner, PCSHELL)); /* compute the action of the inverse of the local Schur complement with a PCSHELL */ 2300 b[0] = subA[0]; 2301 b[3] = data->aux; 2302 PetscCall(MatCreateNest(PETSC_COMM_SELF, 2, nullptr, 2, nullptr, b, &N)); /* instead of computing inv(A11 - A10 inv(A00) A01), compute inv([A00, A01; A10, A11]) followed by a partial solution associated to the A11 block */ 2303 PetscCall(PetscObjectDereference((PetscObject)b[1])); 2304 PetscCall(PetscObjectDereference((PetscObject)b[2])); 2305 PetscCall(PCCreate(PETSC_COMM_SELF, &s)); 2306 PetscCall(PCSetOptionsPrefix(s, ((PetscObject)inner)->prefix)); 2307 PetscCall(PCSetOptionsPrefix(inner, nullptr)); 2308 PetscCall(KSPSetSkipPCSetFromOptions(ksp[0], PETSC_TRUE)); 2309 PetscCall(PCSetType(s, PCLU)); 2310 if (PetscDefined(HAVE_MUMPS)) { /* only MATSOLVERMUMPS handles MATNEST, so for the others, e.g., MATSOLVERPETSC or MATSOLVERMKL_PARDISO, convert to plain MATAIJ */ 2311 PetscCall(PCFactorSetMatSolverType(s, MATSOLVERMUMPS)); 2312 } 2313 PetscCall(PCSetFromOptions(s)); 2314 PetscCall(PCFactorGetMatSolverType(s, &type)); 2315 PetscCall(PetscStrcmp(type, MATSOLVERMUMPS, &flg)); 2316 if (flg) { 2317 PetscCall(PCSetOperators(s, N, N)); 2318 PetscCall(PCFactorGetMatrix(s, b)); 2319 PetscCall(MatSetOptionsPrefix(*b, ((PetscObject)s)->prefix)); 2320 n = -1; 2321 PetscCall(PetscOptionsGetInt(nullptr, ((PetscObject)s)->prefix, "-mat_mumps_icntl_26", &n, nullptr)); 2322 if (n == 1) { 2323 PetscCall(MatNestGetISs(N, is, nullptr)); /* allocates a square MatDense of size is[1]->map->n, so one */ 2324 PetscCall(MatFactorSetSchurIS(*b, is[1])); /* needs to be able to deactivate this path when dealing */ 2325 } /* with a large constraint space in order to avoid OOM */ 2326 } else { 2327 PetscCall(MatConvert(N, MATAIJ, MAT_INITIAL_MATRIX, b)); 2328 PetscCall(PCSetOperators(s, N, *b)); 2329 PetscCall(PetscObjectDereference((PetscObject)*b)); 2330 PetscCall(PCFactorGetMatrix(s, b)); /* MATSOLVERMKL_PARDISO cannot compute in PETSc (yet) a partial solution associated to the A11 block, only partial solution associated to the A00 block or full solution */ 2331 } 2332 PetscCall(PetscNew(&p)); 2333 p->first = s; 2334 PetscCall(MatCreateVecs(*b, p->second, p->second + 1)); 2335 PetscCall(PCShellSetContext(inner, p)); 2336 PetscCall(PCShellSetApply(inner, PCApply_Nest)); 2337 PetscCall(PCShellSetView(inner, PCView_Nest)); 2338 PetscCall(PCShellSetDestroy(inner, PCDestroy_Nest)); 2339 PetscCall(PetscObjectDereference((PetscObject)N)); 2340 } 2341 if (!data->levels[0]->scatter) { 2342 PetscCall(MatCreateVecs(P, &xin, nullptr)); 2343 if (ismatis) PetscCall(MatDestroy(&P)); 2344 PetscCall(VecScatterCreate(xin, data->is, data->levels[0]->D, nullptr, &data->levels[0]->scatter)); 2345 PetscCall(VecDestroy(&xin)); 2346 } 2347 if (data->levels[0]->P) { 2348 /* if the pattern is the same and PCSetUp() has previously succeeded, reuse HPDDM buffers and connectivity */ 2349 PetscCall(HPDDM::Schwarz<PetscScalar>::destroy(data->levels[0], pc->setupcalled < 1 || pc->flag == DIFFERENT_NONZERO_PATTERN ? PETSC_TRUE : PETSC_FALSE)); 2350 } 2351 if (!data->levels[0]->P) data->levels[0]->P = new HPDDM::Schwarz<PetscScalar>(); 2352 if (data->log_separate) PetscCall(PetscLogEventBegin(PC_HPDDM_SetUp[0], data->levels[0]->ksp, nullptr, nullptr, nullptr)); 2353 else PetscCall(PetscLogEventBegin(PC_HPDDM_Strc, data->levels[0]->ksp, nullptr, nullptr, nullptr)); 2354 /* HPDDM internal data structure */ 2355 PetscCall(data->levels[0]->P->structure(loc, data->is, !ctx ? sub[0] : nullptr, ismatis ? C : data->aux, data->levels)); 2356 if (!data->log_separate) PetscCall(PetscLogEventEnd(PC_HPDDM_Strc, data->levels[0]->ksp, nullptr, nullptr, nullptr)); 2357 /* matrix pencil of the generalized eigenvalue problem on the overlap (GenEO) */ 2358 if (!ctx) { 2359 if (data->deflation || overlap != -1) weighted = data->aux; 2360 else if (!data->B) { 2361 PetscBool cmp; 2362 PetscCall(MatDuplicate(sub[0], MAT_COPY_VALUES, &weighted)); 2363 PetscCall(PetscObjectTypeCompareAny((PetscObject)weighted, &cmp, MATNORMAL, MATNORMALHERMITIAN, "")); 2364 if (cmp) flg = PETSC_FALSE; 2365 PetscCall(MatDiagonalScale(weighted, data->levels[0]->D, data->levels[0]->D)); 2366 /* neither MatDuplicate() nor MatDiagonalScale() handles the symmetry options, so propagate the options explicitly */ 2367 /* only useful for -mat_type baij -pc_hpddm_levels_1_st_pc_type cholesky (no problem with MATAIJ or MATSBAIJ) */ 2368 PetscCall(MatPropagateSymmetryOptions(sub[0], weighted)); 2369 } else weighted = data->B; 2370 } else weighted = nullptr; 2371 /* SLEPc is used inside the loaded symbol */ 2372 PetscCall((*loadedSym)(data->levels[0]->P, data->is, ismatis ? C : (algebraic && !block && overlap == -1 ? sub[0] : (!ctx ? data->aux : S)), weighted, data->B, initial, data->levels)); 2373 if (!ctx && data->share && overlap == -1) { 2374 Mat st[2]; 2375 PetscCall(KSPGetOperators(ksp[0], st, st + 1)); 2376 PetscCall(MatCopy(subA[0], st[0], structure)); 2377 if (subA[1] != subA[0] || st[1] != st[0]) PetscCall(MatCopy(subA[1], st[1], SAME_NONZERO_PATTERN)); 2378 } 2379 if (data->log_separate) PetscCall(PetscLogEventEnd(PC_HPDDM_SetUp[0], data->levels[0]->ksp, nullptr, nullptr, nullptr)); 2380 if (ismatis) PetscCall(MatISGetLocalMat(C, &N)); 2381 else N = data->aux; 2382 if (!ctx) P = sub[0]; 2383 else P = S; 2384 /* going through the grid hierarchy */ 2385 for (n = 1; n < data->N; ++n) { 2386 if (data->log_separate) PetscCall(PetscLogEventBegin(PC_HPDDM_SetUp[n], data->levels[n]->ksp, nullptr, nullptr, nullptr)); 2387 /* method composed in the loaded symbol since there, SLEPc is used as well */ 2388 PetscTryMethod(data->levels[0]->ksp, "PCHPDDMSetUp_Private_C", (Mat *, Mat *, PetscInt, PetscInt *const, PC_HPDDM_Level **const), (&P, &N, n, &data->N, data->levels)); 2389 if (data->log_separate) PetscCall(PetscLogEventEnd(PC_HPDDM_SetUp[n], data->levels[n]->ksp, nullptr, nullptr, nullptr)); 2390 } 2391 /* reset to NULL to avoid any faulty use */ 2392 PetscCall(PetscObjectComposeFunction((PetscObject)data->levels[0]->ksp, "PCHPDDMSetUp_Private_C", nullptr)); 2393 if (!ismatis) PetscCall(PetscObjectComposeFunction((PetscObject)pc->pmat, "PCHPDDMAlgebraicAuxiliaryMat_C", nullptr)); 2394 else PetscCall(PetscObjectDereference((PetscObject)C)); /* matching PetscObjectReference() above */ 2395 for (n = 0; n < data->N - 1; ++n) 2396 if (data->levels[n]->P) { 2397 /* HPDDM internal work buffers */ 2398 data->levels[n]->P->setBuffer(); 2399 data->levels[n]->P->super::start(); 2400 } 2401 if (ismatis || !subdomains) PetscCall(PCHPDDMDestroySubMatrices_Private(PetscBool3ToBool(data->Neumann), PetscBool(algebraic && !block && overlap == -1), sub)); 2402 if (ismatis) data->is = nullptr; 2403 for (n = 0; n < data->N - 1 + (reused > 0); ++n) { 2404 if (data->levels[n]->P) { 2405 PC spc; 2406 2407 /* force the PC to be PCSHELL to do the coarse grid corrections */ 2408 PetscCall(KSPSetSkipPCSetFromOptions(data->levels[n]->ksp, PETSC_TRUE)); 2409 PetscCall(KSPGetPC(data->levels[n]->ksp, &spc)); 2410 PetscCall(PCSetType(spc, PCSHELL)); 2411 PetscCall(PCShellSetContext(spc, data->levels[n])); 2412 PetscCall(PCShellSetSetUp(spc, PCSetUp_HPDDMShell)); 2413 PetscCall(PCShellSetApply(spc, PCApply_HPDDMShell)); 2414 PetscCall(PCShellSetMatApply(spc, PCMatApply_HPDDMShell)); 2415 if (ctx && n == 0) { 2416 Mat Amat, Pmat; 2417 PetscInt m, M; 2418 std::tuple<Mat, VecScatter, Vec[2]> *ctx; 2419 2420 PetscCall(KSPGetOperators(data->levels[n]->ksp, nullptr, &Pmat)); 2421 PetscCall(MatGetLocalSize(Pmat, &m, nullptr)); 2422 PetscCall(MatGetSize(Pmat, &M, nullptr)); 2423 PetscCall(PetscNew(&ctx)); 2424 std::get<0>(*ctx) = S; 2425 std::get<1>(*ctx) = data->levels[n]->scatter; 2426 PetscCall(MatCreateShell(PetscObjectComm((PetscObject)data->levels[n]->ksp), m, m, M, M, ctx, &Amat)); 2427 PetscCall(MatShellSetOperation(Amat, MATOP_MULT, (void (*)(void))MatMult_Schur<false>)); 2428 PetscCall(MatShellSetOperation(Amat, MATOP_MULT_TRANSPOSE, (void (*)(void))MatMult_Schur<true>)); 2429 PetscCall(MatShellSetOperation(Amat, MATOP_DESTROY, (void (*)(void))MatDestroy_Schur)); 2430 PetscCall(MatCreateVecs(S, std::get<2>(*ctx), std::get<2>(*ctx) + 1)); 2431 PetscCall(KSPSetOperators(data->levels[n]->ksp, Amat, Pmat)); 2432 PetscCall(PetscObjectDereference((PetscObject)Amat)); 2433 } 2434 PetscCall(PCShellSetDestroy(spc, PCDestroy_HPDDMShell)); 2435 if (!data->levels[n]->pc) PetscCall(PCCreate(PetscObjectComm((PetscObject)data->levels[n]->ksp), &data->levels[n]->pc)); 2436 if (n < reused) { 2437 PetscCall(PCSetReusePreconditioner(spc, PETSC_TRUE)); 2438 PetscCall(PCSetReusePreconditioner(data->levels[n]->pc, PETSC_TRUE)); 2439 } 2440 PetscCall(PCSetUp(spc)); 2441 } 2442 } 2443 if (ctx) PetscCall(MatDestroy(&S)); 2444 if (overlap == -1) PetscCall(PetscObjectComposeFunction((PetscObject)pc->pmat, "PCHPDDMAlgebraicAuxiliaryMat_Private_C", nullptr)); 2445 } else flg = reused ? PETSC_FALSE : PETSC_TRUE; 2446 if (!ismatis && subdomains) { 2447 if (flg) PetscCall(KSPGetPC(data->levels[0]->ksp, &inner)); 2448 else inner = data->levels[0]->pc; 2449 if (inner) { 2450 if (!inner->setupcalled) PetscCall(PCSetType(inner, PCASM)); 2451 PetscCall(PCSetFromOptions(inner)); 2452 PetscCall(PetscStrcmp(((PetscObject)inner)->type_name, PCASM, &flg)); 2453 if (flg) { 2454 if (!inner->setupcalled) { /* evaluates to PETSC_FALSE when -pc_hpddm_block_splitting */ 2455 IS sorted; /* PCASM will sort the input IS, duplicate it to return an unmodified (PCHPDDM) input IS */ 2456 PetscCall(ISDuplicate(is[0], &sorted)); 2457 PetscCall(PCASMSetLocalSubdomains(inner, 1, &sorted, &loc)); 2458 PetscCall(PetscObjectDereference((PetscObject)sorted)); 2459 } 2460 if (!PetscBool3ToBool(data->Neumann) && data->N > 1) { /* subdomain matrices are already created for the eigenproblem, reuse them for the fine-level PC */ 2461 PetscCall(PCHPDDMPermute_Private(*is, nullptr, nullptr, sub[0], &P, nullptr)); 2462 PetscCall(PCHPDDMCommunicationAvoidingPCASM_Private(inner, P, algebraic)); 2463 PetscCall(PetscObjectDereference((PetscObject)P)); 2464 } 2465 } 2466 } 2467 if (data->N > 1) { 2468 if (overlap != 1) PetscCall(PCHPDDMDestroySubMatrices_Private(PetscBool3ToBool(data->Neumann), PetscBool(algebraic && !block && overlap == -1), sub)); 2469 if (overlap == 1) PetscCall(MatDestroy(subA)); 2470 } 2471 } 2472 PetscCall(ISDestroy(&loc)); 2473 } else data->N = 1 + reused; /* enforce this value to 1 + reused if there is no way to build another level */ 2474 if (requested != data->N + reused) { 2475 PetscCall(PetscInfo(pc, "%" PetscInt_FMT " levels requested, only %" PetscInt_FMT " built + %" PetscInt_FMT " reused. Options for level(s) > %" PetscInt_FMT ", including -%spc_hpddm_coarse_ will not be taken into account\n", requested, data->N, reused, 2476 data->N, pcpre ? pcpre : "")); 2477 PetscCall(PetscInfo(pc, "It is best to tune parameters, e.g., a higher value for -%spc_hpddm_levels_%" PetscInt_FMT "_eps_threshold so that at least one local deflation vector will be selected\n", pcpre ? pcpre : "", data->N)); 2478 /* cannot use PCDestroy_HPDDMShell() because PCSHELL not set for unassembled levels */ 2479 for (n = data->N - 1; n < requested - 1; ++n) { 2480 if (data->levels[n]->P) { 2481 PetscCall(HPDDM::Schwarz<PetscScalar>::destroy(data->levels[n], PETSC_TRUE)); 2482 PetscCall(VecDestroyVecs(1, &data->levels[n]->v[0])); 2483 PetscCall(VecDestroyVecs(2, &data->levels[n]->v[1])); 2484 PetscCall(MatDestroy(data->levels[n]->V)); 2485 PetscCall(MatDestroy(data->levels[n]->V + 1)); 2486 PetscCall(MatDestroy(data->levels[n]->V + 2)); 2487 PetscCall(VecDestroy(&data->levels[n]->D)); 2488 PetscCall(VecScatterDestroy(&data->levels[n]->scatter)); 2489 } 2490 } 2491 if (reused) { 2492 for (n = reused; n < PETSC_PCHPDDM_MAXLEVELS && data->levels[n]; ++n) { 2493 PetscCall(KSPDestroy(&data->levels[n]->ksp)); 2494 PetscCall(PCDestroy(&data->levels[n]->pc)); 2495 } 2496 } 2497 PetscCheck(!PetscDefined(USE_DEBUG), PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "%" PetscInt_FMT " levels requested, only %" PetscInt_FMT " built + %" PetscInt_FMT " reused. Options for level(s) > %" PetscInt_FMT ", including -%spc_hpddm_coarse_ will not be taken into account. It is best to tune parameters, e.g., a higher value for -%spc_hpddm_levels_%" PetscInt_FMT "_eps_threshold so that at least one local deflation vector will be selected. If you don't want this to error out, compile --with-debugging=0", requested, 2498 data->N, reused, data->N, pcpre ? pcpre : "", pcpre ? pcpre : "", data->N); 2499 } 2500 /* these solvers are created after PCSetFromOptions() is called */ 2501 if (pc->setfromoptionscalled) { 2502 for (n = 0; n < data->N; ++n) { 2503 if (data->levels[n]->ksp) PetscCall(KSPSetFromOptions(data->levels[n]->ksp)); 2504 if (data->levels[n]->pc) PetscCall(PCSetFromOptions(data->levels[n]->pc)); 2505 } 2506 pc->setfromoptionscalled = 0; 2507 } 2508 data->N += reused; 2509 if (data->share && swap) { 2510 /* swap back pointers so that variables follow the user-provided numbering */ 2511 std::swap(C, data->aux); 2512 std::swap(uis, data->is); 2513 PetscCall(MatDestroy(&C)); 2514 PetscCall(ISDestroy(&uis)); 2515 } 2516 if (algebraic) PetscCall(MatDestroy(&data->aux)); 2517 if (unsorted && unsorted != is[0]) { 2518 PetscCall(ISCopy(unsorted, data->is)); 2519 PetscCall(ISDestroy(&unsorted)); 2520 } 2521 #if PetscDefined(USE_DEBUG) 2522 PetscCheck((data->is && dis) || (!data->is && !dis), PETSC_COMM_SELF, PETSC_ERR_PLIB, "An IS pointer is NULL but not the other: input IS pointer (%p) v. output IS pointer (%p)", (void *)dis, (void *)data->is); 2523 if (data->is) { 2524 PetscCall(ISEqualUnsorted(data->is, dis, &flg)); 2525 PetscCall(ISDestroy(&dis)); 2526 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Input IS and output IS are not equal"); 2527 } 2528 PetscCheck((data->aux && daux) || (!data->aux && !daux), PETSC_COMM_SELF, PETSC_ERR_PLIB, "A Mat pointer is NULL but not the other: input Mat pointer (%p) v. output Mat pointer (%p)", (void *)daux, (void *)data->aux); 2529 if (data->aux) { 2530 PetscCall(MatMultEqual(data->aux, daux, 20, &flg)); 2531 PetscCall(MatDestroy(&daux)); 2532 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Input Mat and output Mat are not equal"); 2533 } 2534 #endif 2535 PetscFunctionReturn(PETSC_SUCCESS); 2536 } 2537 2538 /*@ 2539 PCHPDDMSetCoarseCorrectionType - Sets the coarse correction type. 2540 2541 Collective 2542 2543 Input Parameters: 2544 + pc - preconditioner context 2545 - type - `PC_HPDDM_COARSE_CORRECTION_DEFLATED`, `PC_HPDDM_COARSE_CORRECTION_ADDITIVE`, `PC_HPDDM_COARSE_CORRECTION_BALANCED`, or `PC_HPDDM_COARSE_CORRECTION_NONE` 2546 2547 Options Database Key: 2548 . -pc_hpddm_coarse_correction <deflated, additive, balanced, none> - type of coarse correction to apply 2549 2550 Level: intermediate 2551 2552 .seealso: [](ch_ksp), `PCHPDDMGetCoarseCorrectionType()`, `PCHPDDM`, `PCHPDDMCoarseCorrectionType` 2553 @*/ 2554 PetscErrorCode PCHPDDMSetCoarseCorrectionType(PC pc, PCHPDDMCoarseCorrectionType type) 2555 { 2556 PetscFunctionBegin; 2557 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2558 PetscValidLogicalCollectiveEnum(pc, type, 2); 2559 PetscTryMethod(pc, "PCHPDDMSetCoarseCorrectionType_C", (PC, PCHPDDMCoarseCorrectionType), (pc, type)); 2560 PetscFunctionReturn(PETSC_SUCCESS); 2561 } 2562 2563 /*@ 2564 PCHPDDMGetCoarseCorrectionType - Gets the coarse correction type. 2565 2566 Input Parameter: 2567 . pc - preconditioner context 2568 2569 Output Parameter: 2570 . type - `PC_HPDDM_COARSE_CORRECTION_DEFLATED`, `PC_HPDDM_COARSE_CORRECTION_ADDITIVE`, `PC_HPDDM_COARSE_CORRECTION_BALANCED`, or `PC_HPDDM_COARSE_CORRECTION_NONE` 2571 2572 Level: intermediate 2573 2574 .seealso: [](ch_ksp), `PCHPDDMSetCoarseCorrectionType()`, `PCHPDDM`, `PCHPDDMCoarseCorrectionType` 2575 @*/ 2576 PetscErrorCode PCHPDDMGetCoarseCorrectionType(PC pc, PCHPDDMCoarseCorrectionType *type) 2577 { 2578 PetscFunctionBegin; 2579 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2580 if (type) { 2581 PetscAssertPointer(type, 2); 2582 PetscUseMethod(pc, "PCHPDDMGetCoarseCorrectionType_C", (PC, PCHPDDMCoarseCorrectionType *), (pc, type)); 2583 } 2584 PetscFunctionReturn(PETSC_SUCCESS); 2585 } 2586 2587 static PetscErrorCode PCHPDDMSetCoarseCorrectionType_HPDDM(PC pc, PCHPDDMCoarseCorrectionType type) 2588 { 2589 PC_HPDDM *data = (PC_HPDDM *)pc->data; 2590 2591 PetscFunctionBegin; 2592 data->correction = type; 2593 PetscFunctionReturn(PETSC_SUCCESS); 2594 } 2595 2596 static PetscErrorCode PCHPDDMGetCoarseCorrectionType_HPDDM(PC pc, PCHPDDMCoarseCorrectionType *type) 2597 { 2598 PC_HPDDM *data = (PC_HPDDM *)pc->data; 2599 2600 PetscFunctionBegin; 2601 *type = data->correction; 2602 PetscFunctionReturn(PETSC_SUCCESS); 2603 } 2604 2605 /*@ 2606 PCHPDDMSetSTShareSubKSP - Sets whether the `KSP` in SLEPc `ST` and the fine-level subdomain solver should be shared. 2607 2608 Input Parameters: 2609 + pc - preconditioner context 2610 - share - whether the `KSP` should be shared or not 2611 2612 Note: 2613 This is not the same as `PCSetReusePreconditioner()`. Given certain conditions (visible using -info), a symbolic factorization can be skipped 2614 when using a subdomain `PCType` such as `PCLU` or `PCCHOLESKY`. 2615 2616 Level: advanced 2617 2618 .seealso: [](ch_ksp), `PCHPDDM`, `PCHPDDMGetSTShareSubKSP()` 2619 @*/ 2620 PetscErrorCode PCHPDDMSetSTShareSubKSP(PC pc, PetscBool share) 2621 { 2622 PetscFunctionBegin; 2623 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2624 PetscTryMethod(pc, "PCHPDDMSetSTShareSubKSP_C", (PC, PetscBool), (pc, share)); 2625 PetscFunctionReturn(PETSC_SUCCESS); 2626 } 2627 2628 /*@ 2629 PCHPDDMGetSTShareSubKSP - Gets whether the `KSP` in SLEPc `ST` and the fine-level subdomain solver is shared. 2630 2631 Input Parameter: 2632 . pc - preconditioner context 2633 2634 Output Parameter: 2635 . share - whether the `KSP` is shared or not 2636 2637 Note: 2638 This is not the same as `PCGetReusePreconditioner()`. The return value is unlikely to be true, but when it is, a symbolic factorization can be skipped 2639 when using a subdomain `PCType` such as `PCLU` or `PCCHOLESKY`. 2640 2641 Level: advanced 2642 2643 .seealso: [](ch_ksp), `PCHPDDM`, `PCHPDDMSetSTShareSubKSP()` 2644 @*/ 2645 PetscErrorCode PCHPDDMGetSTShareSubKSP(PC pc, PetscBool *share) 2646 { 2647 PetscFunctionBegin; 2648 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2649 if (share) { 2650 PetscAssertPointer(share, 2); 2651 PetscUseMethod(pc, "PCHPDDMGetSTShareSubKSP_C", (PC, PetscBool *), (pc, share)); 2652 } 2653 PetscFunctionReturn(PETSC_SUCCESS); 2654 } 2655 2656 static PetscErrorCode PCHPDDMSetSTShareSubKSP_HPDDM(PC pc, PetscBool share) 2657 { 2658 PC_HPDDM *data = (PC_HPDDM *)pc->data; 2659 2660 PetscFunctionBegin; 2661 data->share = share; 2662 PetscFunctionReturn(PETSC_SUCCESS); 2663 } 2664 2665 static PetscErrorCode PCHPDDMGetSTShareSubKSP_HPDDM(PC pc, PetscBool *share) 2666 { 2667 PC_HPDDM *data = (PC_HPDDM *)pc->data; 2668 2669 PetscFunctionBegin; 2670 *share = data->share; 2671 PetscFunctionReturn(PETSC_SUCCESS); 2672 } 2673 2674 /*@ 2675 PCHPDDMSetDeflationMat - Sets the deflation space used to assemble a coarser operator. 2676 2677 Input Parameters: 2678 + pc - preconditioner context 2679 . is - index set of the local deflation matrix 2680 - U - deflation sequential matrix stored as a `MATSEQDENSE` 2681 2682 Level: advanced 2683 2684 .seealso: [](ch_ksp), `PCHPDDM`, `PCDeflationSetSpace()`, `PCMGSetRestriction()` 2685 @*/ 2686 PetscErrorCode PCHPDDMSetDeflationMat(PC pc, IS is, Mat U) 2687 { 2688 PetscFunctionBegin; 2689 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2690 PetscValidHeaderSpecific(is, IS_CLASSID, 2); 2691 PetscValidHeaderSpecific(U, MAT_CLASSID, 3); 2692 PetscTryMethod(pc, "PCHPDDMSetDeflationMat_C", (PC, IS, Mat), (pc, is, U)); 2693 PetscFunctionReturn(PETSC_SUCCESS); 2694 } 2695 2696 static PetscErrorCode PCHPDDMSetDeflationMat_HPDDM(PC pc, IS is, Mat U) 2697 { 2698 PetscFunctionBegin; 2699 PetscCall(PCHPDDMSetAuxiliaryMat_Private(pc, is, U, PETSC_TRUE)); 2700 PetscFunctionReturn(PETSC_SUCCESS); 2701 } 2702 2703 PetscErrorCode HPDDMLoadDL_Private(PetscBool *found) 2704 { 2705 PetscBool flg; 2706 char lib[PETSC_MAX_PATH_LEN], dlib[PETSC_MAX_PATH_LEN], dir[PETSC_MAX_PATH_LEN]; 2707 2708 PetscFunctionBegin; 2709 PetscAssertPointer(found, 1); 2710 PetscCall(PetscStrncpy(dir, "${PETSC_LIB_DIR}", sizeof(dir))); 2711 PetscCall(PetscOptionsGetString(nullptr, nullptr, "-hpddm_dir", dir, sizeof(dir), nullptr)); 2712 PetscCall(PetscSNPrintf(lib, sizeof(lib), "%s/libhpddm_petsc", dir)); 2713 PetscCall(PetscDLLibraryRetrieve(PETSC_COMM_SELF, lib, dlib, 1024, found)); 2714 #if defined(SLEPC_LIB_DIR) /* this variable is passed during SLEPc ./configure when PETSc has not been configured */ 2715 if (!*found) { /* with --download-hpddm since slepcconf.h is not yet built (and thus can't be included) */ 2716 PetscCall(PetscStrncpy(dir, HPDDM_STR(SLEPC_LIB_DIR), sizeof(dir))); 2717 PetscCall(PetscSNPrintf(lib, sizeof(lib), "%s/libhpddm_petsc", dir)); 2718 PetscCall(PetscDLLibraryRetrieve(PETSC_COMM_SELF, lib, dlib, 1024, found)); 2719 } 2720 #endif 2721 if (!*found) { /* probable options for this to evaluate to PETSC_TRUE: system inconsistency (libhpddm_petsc moved by user?) or PETSc configured without --download-slepc */ 2722 PetscCall(PetscOptionsGetenv(PETSC_COMM_SELF, "SLEPC_DIR", dir, sizeof(dir), &flg)); 2723 if (flg) { /* if both PETSc and SLEPc are configured with --download-hpddm but PETSc has been configured without --download-slepc, one must ensure that libslepc is loaded before libhpddm_petsc */ 2724 PetscCall(PetscSNPrintf(lib, sizeof(lib), "%s/lib/libslepc", dir)); 2725 PetscCall(PetscDLLibraryRetrieve(PETSC_COMM_SELF, lib, dlib, 1024, found)); 2726 PetscCheck(*found, PETSC_COMM_SELF, PETSC_ERR_PLIB, "%s not found but SLEPC_DIR=%s", lib, dir); 2727 PetscCall(PetscDLLibraryAppend(PETSC_COMM_SELF, &PetscDLLibrariesLoaded, dlib)); 2728 PetscCall(PetscSNPrintf(lib, sizeof(lib), "%s/lib/libhpddm_petsc", dir)); /* libhpddm_petsc is always in the same directory as libslepc */ 2729 PetscCall(PetscDLLibraryRetrieve(PETSC_COMM_SELF, lib, dlib, 1024, found)); 2730 } 2731 } 2732 PetscCheck(*found, PETSC_COMM_SELF, PETSC_ERR_PLIB, "%s not found", lib); 2733 PetscCall(PetscDLLibraryAppend(PETSC_COMM_SELF, &PetscDLLibrariesLoaded, dlib)); 2734 PetscFunctionReturn(PETSC_SUCCESS); 2735 } 2736 2737 /*MC 2738 PCHPDDM - Interface with the HPDDM library. 2739 2740 This `PC` may be used to build multilevel spectral domain decomposition methods based on the GenEO framework {cite}`spillane2011robust` {cite}`al2021multilevel`. 2741 It may be viewed as an alternative to spectral 2742 AMGe or `PCBDDC` with adaptive selection of constraints. The interface is explained in details in {cite}`jolivetromanzampini2020` 2743 2744 The matrix to be preconditioned (Pmat) may be unassembled (`MATIS`), assembled (`MATAIJ`, `MATBAIJ`, or `MATSBAIJ`), hierarchical (`MATHTOOL`), or `MATNORMAL`. 2745 2746 For multilevel preconditioning, when using an assembled or hierarchical Pmat, one must provide an auxiliary local `Mat` (unassembled local operator for GenEO) using 2747 `PCHPDDMSetAuxiliaryMat()`. Calling this routine is not needed when using a `MATIS` Pmat, assembly is done internally using `MatConvert()`. 2748 2749 Options Database Keys: 2750 + -pc_hpddm_define_subdomains <true, default=false> - on the finest level, calls `PCASMSetLocalSubdomains()` with the `IS` supplied in `PCHPDDMSetAuxiliaryMat()` 2751 (not relevant with an unassembled Pmat) 2752 . -pc_hpddm_has_neumann <true, default=false> - on the finest level, informs the `PC` that the local Neumann matrix is supplied in `PCHPDDMSetAuxiliaryMat()` 2753 - -pc_hpddm_coarse_correction <type, default=deflated> - determines the `PCHPDDMCoarseCorrectionType` when calling `PCApply()` 2754 2755 Options for subdomain solvers, subdomain eigensolvers (for computing deflation vectors), and the coarse solver can be set using the following options database prefixes. 2756 .vb 2757 -pc_hpddm_levels_%d_pc_ 2758 -pc_hpddm_levels_%d_ksp_ 2759 -pc_hpddm_levels_%d_eps_ 2760 -pc_hpddm_levels_%d_p 2761 -pc_hpddm_levels_%d_mat_type 2762 -pc_hpddm_coarse_ 2763 -pc_hpddm_coarse_p 2764 -pc_hpddm_coarse_mat_type 2765 -pc_hpddm_coarse_mat_filter 2766 .ve 2767 2768 E.g., -pc_hpddm_levels_1_sub_pc_type lu -pc_hpddm_levels_1_eps_nev 10 -pc_hpddm_levels_2_p 4 -pc_hpddm_levels_2_sub_pc_type lu -pc_hpddm_levels_2_eps_nev 10 2769 -pc_hpddm_coarse_p 2 -pc_hpddm_coarse_mat_type baij will use 10 deflation vectors per subdomain on the fine "level 1", 2770 aggregate the fine subdomains into 4 "level 2" subdomains, then use 10 deflation vectors per subdomain on "level 2", 2771 and assemble the coarse matrix (of dimension 4 x 10 = 40) on two processes as a `MATBAIJ` (default is `MATSBAIJ`). 2772 2773 In order to activate a "level N+1" coarse correction, it is mandatory to call -pc_hpddm_levels_N_eps_nev <nu> or -pc_hpddm_levels_N_eps_threshold <val>. The default -pc_hpddm_coarse_p value is 1, meaning that the coarse operator is aggregated on a single process. 2774 2775 Level: intermediate 2776 2777 Notes: 2778 This preconditioner requires that PETSc is built with SLEPc (``--download-slepc``). 2779 2780 By default, the underlying concurrent eigenproblems 2781 are solved using SLEPc shift-and-invert spectral transformation. This is usually what gives the best performance for GenEO, cf. 2782 {cite}`spillane2011robust` {cite}`jolivet2013scalabledd`. As 2783 stated above, SLEPc options are available through -pc_hpddm_levels_%d_, e.g., -pc_hpddm_levels_1_eps_type arpack -pc_hpddm_levels_1_eps_nev 10 2784 -pc_hpddm_levels_1_st_type sinvert. There are furthermore three options related to the (subdomain-wise local) eigensolver that are not described in 2785 SLEPc documentation since they are specific to `PCHPDDM`. 2786 .vb 2787 -pc_hpddm_levels_1_st_share_sub_ksp 2788 -pc_hpddm_levels_%d_eps_threshold 2789 -pc_hpddm_levels_1_eps_use_inertia 2790 .ve 2791 2792 The first option from the list only applies to the fine-level eigensolver, see `PCHPDDMSetSTShareSubKSP()`. The second option from the list is 2793 used to filter eigenmodes retrieved after convergence of `EPSSolve()` at "level N" such that eigenvectors used to define a "level N+1" coarse 2794 correction are associated to eigenvalues whose magnitude are lower or equal than -pc_hpddm_levels_N_eps_threshold. When using an `EPS` which cannot 2795 determine a priori the proper -pc_hpddm_levels_N_eps_nev such that all wanted eigenmodes are retrieved, it is possible to get an estimation of the 2796 correct value using the third option from the list, -pc_hpddm_levels_1_eps_use_inertia, see `MatGetInertia()`. In that case, there is no need 2797 to supply -pc_hpddm_levels_1_eps_nev. This last option also only applies to the fine-level (N = 1) eigensolver. 2798 2799 See also {cite}`dolean2015introduction`, {cite}`al2022robust`, {cite}`al2022robustpd`, and {cite}`nataf2022recent` 2800 2801 .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCHPDDMSetAuxiliaryMat()`, `MATIS`, `PCBDDC`, `PCDEFLATION`, `PCTELESCOPE`, `PCASM`, 2802 `PCHPDDMSetCoarseCorrectionType()`, `PCHPDDMHasNeumannMat()`, `PCHPDDMSetRHSMat()`, `PCHPDDMSetDeflationMat()`, `PCHPDDMSetSTShareSubKSP()`, 2803 `PCHPDDMGetSTShareSubKSP()`, `PCHPDDMGetCoarseCorrectionType()`, `PCHPDDMGetComplexities()` 2804 M*/ 2805 PETSC_EXTERN PetscErrorCode PCCreate_HPDDM(PC pc) 2806 { 2807 PC_HPDDM *data; 2808 PetscBool found; 2809 2810 PetscFunctionBegin; 2811 if (!loadedSym) { 2812 PetscCall(HPDDMLoadDL_Private(&found)); 2813 if (found) PetscCall(PetscDLLibrarySym(PETSC_COMM_SELF, &PetscDLLibrariesLoaded, nullptr, "PCHPDDM_Internal", (void **)&loadedSym)); 2814 } 2815 PetscCheck(loadedSym, PETSC_COMM_SELF, PETSC_ERR_PLIB, "PCHPDDM_Internal symbol not found in loaded libhpddm_petsc"); 2816 PetscCall(PetscNew(&data)); 2817 pc->data = data; 2818 data->Neumann = PETSC_BOOL3_UNKNOWN; 2819 pc->ops->reset = PCReset_HPDDM; 2820 pc->ops->destroy = PCDestroy_HPDDM; 2821 pc->ops->setfromoptions = PCSetFromOptions_HPDDM; 2822 pc->ops->setup = PCSetUp_HPDDM; 2823 pc->ops->apply = PCApply_HPDDM; 2824 pc->ops->matapply = PCMatApply_HPDDM; 2825 pc->ops->view = PCView_HPDDM; 2826 pc->ops->presolve = PCPreSolve_HPDDM; 2827 2828 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetAuxiliaryMat_C", PCHPDDMSetAuxiliaryMat_HPDDM)); 2829 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMHasNeumannMat_C", PCHPDDMHasNeumannMat_HPDDM)); 2830 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetRHSMat_C", PCHPDDMSetRHSMat_HPDDM)); 2831 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetCoarseCorrectionType_C", PCHPDDMSetCoarseCorrectionType_HPDDM)); 2832 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMGetCoarseCorrectionType_C", PCHPDDMGetCoarseCorrectionType_HPDDM)); 2833 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetSTShareSubKSP_C", PCHPDDMSetSTShareSubKSP_HPDDM)); 2834 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMGetSTShareSubKSP_C", PCHPDDMGetSTShareSubKSP_HPDDM)); 2835 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetDeflationMat_C", PCHPDDMSetDeflationMat_HPDDM)); 2836 PetscFunctionReturn(PETSC_SUCCESS); 2837 } 2838 2839 /*@C 2840 PCHPDDMInitializePackage - This function initializes everything in the `PCHPDDM` package. It is called from `PCInitializePackage()`. 2841 2842 Level: developer 2843 2844 .seealso: [](ch_ksp), `PetscInitialize()` 2845 @*/ 2846 PetscErrorCode PCHPDDMInitializePackage(void) 2847 { 2848 char ename[32]; 2849 PetscInt i; 2850 2851 PetscFunctionBegin; 2852 if (PCHPDDMPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS); 2853 PCHPDDMPackageInitialized = PETSC_TRUE; 2854 PetscCall(PetscRegisterFinalize(PCHPDDMFinalizePackage)); 2855 /* general events registered once during package initialization */ 2856 /* some of these events are not triggered in libpetsc, */ 2857 /* but rather directly in libhpddm_petsc, */ 2858 /* which is in charge of performing the following operations */ 2859 2860 /* domain decomposition structure from Pmat sparsity pattern */ 2861 PetscCall(PetscLogEventRegister("PCHPDDMStrc", PC_CLASSID, &PC_HPDDM_Strc)); 2862 /* Galerkin product, redistribution, and setup (not triggered in libpetsc) */ 2863 PetscCall(PetscLogEventRegister("PCHPDDMPtAP", PC_CLASSID, &PC_HPDDM_PtAP)); 2864 /* Galerkin product with summation, redistribution, and setup (not triggered in libpetsc) */ 2865 PetscCall(PetscLogEventRegister("PCHPDDMPtBP", PC_CLASSID, &PC_HPDDM_PtBP)); 2866 /* next level construction using PtAP and PtBP (not triggered in libpetsc) */ 2867 PetscCall(PetscLogEventRegister("PCHPDDMNext", PC_CLASSID, &PC_HPDDM_Next)); 2868 static_assert(PETSC_PCHPDDM_MAXLEVELS <= 9, "PETSC_PCHPDDM_MAXLEVELS value is too high"); 2869 for (i = 1; i < PETSC_PCHPDDM_MAXLEVELS; ++i) { 2870 PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCHPDDMSetUp L%1" PetscInt_FMT, i)); 2871 /* events during a PCSetUp() at level #i _except_ the assembly */ 2872 /* of the Galerkin operator of the coarser level #(i + 1) */ 2873 PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_HPDDM_SetUp[i - 1])); 2874 PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCHPDDMSolve L%1" PetscInt_FMT, i)); 2875 /* events during a PCApply() at level #i _except_ */ 2876 /* the KSPSolve() of the coarser level #(i + 1) */ 2877 PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_HPDDM_Solve[i - 1])); 2878 } 2879 PetscFunctionReturn(PETSC_SUCCESS); 2880 } 2881 2882 /*@C 2883 PCHPDDMFinalizePackage - This function frees everything from the `PCHPDDM` package. It is called from `PetscFinalize()`. 2884 2885 Level: developer 2886 2887 .seealso: [](ch_ksp), `PetscFinalize()` 2888 @*/ 2889 PetscErrorCode PCHPDDMFinalizePackage(void) 2890 { 2891 PetscFunctionBegin; 2892 PCHPDDMPackageInitialized = PETSC_FALSE; 2893 PetscFunctionReturn(PETSC_SUCCESS); 2894 } 2895 2896 static PetscErrorCode MatMult_Harmonic(Mat A, Vec x, Vec y) 2897 { 2898 Harmonic h; /* [ A_00 A_01 ], furthermore, A_00 = [ A_loc,loc A_loc,ovl ], thus, A_01 = [ ] */ 2899 /* [ A_10 A_11 A_12 ] [ A_ovl,loc A_ovl,ovl ] [ A_ovl,1 ] */ 2900 Vec sub; /* y = A x = R_loc R_0 [ A_00 A_01 ]^-1 R_loc = [ I_loc ] */ 2901 /* [ A_10 A_11 ] R_1^T A_12 x [ ] */ 2902 PetscFunctionBegin; 2903 PetscCall(MatShellGetContext(A, &h)); 2904 PetscCall(VecSet(h->v, 0.0)); 2905 PetscCall(VecGetSubVector(h->v, h->is[0], &sub)); 2906 PetscCall(MatMult(h->A[0], x, sub)); 2907 PetscCall(VecRestoreSubVector(h->v, h->is[0], &sub)); 2908 PetscCall(KSPSolve(h->ksp, h->v, h->v)); 2909 PetscCall(VecISCopy(h->v, h->is[1], SCATTER_REVERSE, y)); 2910 PetscFunctionReturn(PETSC_SUCCESS); 2911 } 2912 2913 static PetscErrorCode MatMultTranspose_Harmonic(Mat A, Vec y, Vec x) 2914 { 2915 Harmonic h; /* x = A^T y = [ A_00 A_01 ]^-T R_0^T R_loc^T y */ 2916 Vec sub; /* A_12^T R_1 [ A_10 A_11 ] */ 2917 2918 PetscFunctionBegin; 2919 PetscCall(MatShellGetContext(A, &h)); 2920 PetscCall(VecSet(h->v, 0.0)); 2921 PetscCall(VecISCopy(h->v, h->is[1], SCATTER_FORWARD, y)); 2922 PetscCall(KSPSolveTranspose(h->ksp, h->v, h->v)); 2923 PetscCall(VecGetSubVector(h->v, h->is[0], &sub)); 2924 PetscCall(MatMultTranspose(h->A[0], sub, x)); 2925 PetscCall(VecRestoreSubVector(h->v, h->is[0], &sub)); 2926 PetscFunctionReturn(PETSC_SUCCESS); 2927 } 2928 2929 static PetscErrorCode MatProduct_AB_Harmonic(Mat S, Mat X, Mat Y, void *) 2930 { 2931 Harmonic h; 2932 Mat A, B; 2933 Vec a, b; 2934 2935 PetscFunctionBegin; 2936 PetscCall(MatShellGetContext(S, &h)); 2937 PetscCall(MatMatMult(h->A[0], X, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &A)); 2938 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, h->ksp->pc->mat->rmap->n, A->cmap->n, nullptr, &B)); 2939 for (PetscInt i = 0; i < A->cmap->n; ++i) { 2940 PetscCall(MatDenseGetColumnVecRead(A, i, &a)); 2941 PetscCall(MatDenseGetColumnVecWrite(B, i, &b)); 2942 PetscCall(VecISCopy(b, h->is[0], SCATTER_FORWARD, a)); 2943 PetscCall(MatDenseRestoreColumnVecWrite(B, i, &b)); 2944 PetscCall(MatDenseRestoreColumnVecRead(A, i, &a)); 2945 } 2946 PetscCall(MatDestroy(&A)); 2947 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, h->ksp->pc->mat->rmap->n, B->cmap->n, nullptr, &A)); 2948 PetscCall(KSPMatSolve(h->ksp, B, A)); 2949 PetscCall(MatDestroy(&B)); 2950 for (PetscInt i = 0; i < A->cmap->n; ++i) { 2951 PetscCall(MatDenseGetColumnVecRead(A, i, &a)); 2952 PetscCall(MatDenseGetColumnVecWrite(Y, i, &b)); 2953 PetscCall(VecISCopy(a, h->is[1], SCATTER_REVERSE, b)); 2954 PetscCall(MatDenseRestoreColumnVecWrite(Y, i, &b)); 2955 PetscCall(MatDenseRestoreColumnVecRead(A, i, &a)); 2956 } 2957 PetscCall(MatDestroy(&A)); 2958 PetscFunctionReturn(PETSC_SUCCESS); 2959 } 2960 2961 static PetscErrorCode MatDestroy_Harmonic(Mat A) 2962 { 2963 Harmonic h; 2964 2965 PetscFunctionBegin; 2966 PetscCall(MatShellGetContext(A, &h)); 2967 for (PetscInt i = 0; i < (h->A[1] ? 5 : 3); ++i) PetscCall(ISDestroy(h->is + i)); 2968 PetscCall(PetscFree(h->is)); 2969 PetscCall(VecDestroy(&h->v)); 2970 for (PetscInt i = 0; i < 2; ++i) PetscCall(MatDestroy(h->A + i)); 2971 PetscCall(PetscFree(h->A)); 2972 PetscCall(KSPDestroy(&h->ksp)); 2973 PetscCall(PetscFree(h)); 2974 PetscFunctionReturn(PETSC_SUCCESS); 2975 } 2976