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