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