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