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