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 PetscValidHeaderSpecific(is, IS_CLASSID, 1); 1125 PetscValidHeaderSpecific(in_C, MAT_CLASSID, 4); 1126 PetscCall(ISSorted(is, &sorted)); 1127 if (!sorted) { 1128 PetscCall(ISGetLocalSize(is, &size)); 1129 PetscCall(ISGetIndices(is, &ptr)); 1130 PetscCall(ISGetBlockSize(is, &bs)); 1131 /* MatCreateSubMatrices(), called by PCASM, follows the global numbering of Pmat */ 1132 for (n = 0; n < size; n += bs) order.insert(std::make_pair(ptr[n] / bs, n / bs)); 1133 PetscCall(ISRestoreIndices(is, &ptr)); 1134 size /= bs; 1135 if (out_C) { 1136 PetscCall(PetscMalloc1(size, &concatenate)); 1137 for (const std::pair<const PetscInt, PetscInt> &i : order) *concatenate++ = i.second; 1138 concatenate -= size; 1139 PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)in_C), bs, size, concatenate, PETSC_OWN_POINTER, &perm)); 1140 PetscCall(ISSetPermutation(perm)); 1141 /* permute user-provided Mat so that it matches with MatCreateSubMatrices() numbering */ 1142 PetscCall(MatPermute(in_C, perm, perm, out_C)); 1143 if (p) *p = perm; 1144 else PetscCall(ISDestroy(&perm)); /* no need to save the permutation */ 1145 } 1146 if (out_is) { 1147 PetscCall(PetscMalloc1(size, &concatenate)); 1148 for (const std::pair<const PetscInt, PetscInt> &i : order) *concatenate++ = i.first; 1149 concatenate -= size; 1150 /* permute user-provided IS so that it matches with MatCreateSubMatrices() numbering */ 1151 PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)in_is), bs, size, concatenate, PETSC_OWN_POINTER, out_is)); 1152 } 1153 } else { /* input IS is sorted, nothing to permute, simply duplicate inputs when needed */ 1154 if (out_C) PetscCall(MatDuplicate(in_C, MAT_COPY_VALUES, out_C)); 1155 if (out_is) PetscCall(ISDuplicate(in_is, out_is)); 1156 } 1157 PetscFunctionReturn(PETSC_SUCCESS); 1158 } 1159 1160 static PetscErrorCode PCHPDDMCheckSymmetry_Private(PC pc, Mat A01, Mat A10) 1161 { 1162 Mat T, U = nullptr, B = nullptr; 1163 IS z; 1164 PetscBool flg[2]; 1165 1166 PetscFunctionBegin; 1167 PetscCall(PetscObjectTypeCompare((PetscObject)A10, MATTRANSPOSEVIRTUAL, flg)); 1168 if (flg[0]) PetscCall(MatTransposeGetMat(A10, &U)); 1169 else { 1170 PetscCall(PetscObjectTypeCompare((PetscObject)A10, MATHERMITIANTRANSPOSEVIRTUAL, flg + 1)); 1171 if (flg[1]) PetscCall(MatHermitianTransposeGetMat(A10, &U)); 1172 } 1173 if (U) PetscCall(MatDuplicate(U, MAT_COPY_VALUES, &T)); 1174 else PetscCall(MatHermitianTranspose(A10, MAT_INITIAL_MATRIX, &T)); 1175 PetscCall(PetscObjectTypeCompare((PetscObject)A01, MATTRANSPOSEVIRTUAL, flg)); 1176 if (flg[0]) { 1177 PetscCall(MatTransposeGetMat(A01, &A01)); 1178 PetscCall(MatTranspose(A01, MAT_INITIAL_MATRIX, &B)); 1179 A01 = B; 1180 } else { 1181 PetscCall(PetscObjectTypeCompare((PetscObject)A01, MATHERMITIANTRANSPOSEVIRTUAL, flg)); 1182 if (flg[0]) { 1183 PetscCall(MatHermitianTransposeGetMat(A01, &A01)); 1184 PetscCall(MatHermitianTranspose(A01, MAT_INITIAL_MATRIX, &B)); 1185 A01 = B; 1186 } 1187 } 1188 PetscCall(PetscLayoutCompare(T->rmap, A01->rmap, flg)); 1189 if (flg[0]) { 1190 PetscCall(PetscLayoutCompare(T->cmap, A01->cmap, flg)); 1191 if (flg[0]) { 1192 PetscCall(MatFindZeroRows(A01, &z)); /* for essential boundary conditions, some implementations will */ 1193 if (z) { /* zero rows in [P00 A01] except for the diagonal of P00 */ 1194 PetscCall(MatSetOption(T, MAT_NO_OFF_PROC_ZERO_ROWS, PETSC_TRUE)); 1195 PetscCall(MatZeroRowsIS(T, z, 0.0, nullptr, nullptr)); /* corresponding zero rows from A01 */ 1196 PetscCall(ISDestroy(&z)); 1197 } 1198 PetscCall(MatMultEqual(A01, T, 20, flg)); 1199 PetscCheck(flg[0], PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "A01 != A10^T"); 1200 } else PetscCall(PetscInfo(pc, "A01 and A10^T have non-congruent column layouts, cannot test for equality\n")); 1201 } 1202 PetscCall(MatDestroy(&B)); 1203 PetscCall(MatDestroy(&T)); 1204 PetscFunctionReturn(PETSC_SUCCESS); 1205 } 1206 1207 static PetscErrorCode PCHPDDMDestroySubMatrices_Private(PetscBool flg, PetscBool algebraic, Mat *sub) 1208 { 1209 IS is; 1210 1211 PetscFunctionBegin; 1212 if (!flg) { 1213 if (algebraic) { 1214 PetscCall(PetscObjectQuery((PetscObject)sub[0], "_PCHPDDM_Embed", (PetscObject *)&is)); 1215 PetscCall(ISDestroy(&is)); 1216 PetscCall(PetscObjectCompose((PetscObject)sub[0], "_PCHPDDM_Embed", nullptr)); 1217 PetscCall(PetscObjectCompose((PetscObject)sub[0], "_PCHPDDM_Compact", nullptr)); 1218 } 1219 PetscCall(MatDestroySubMatrices(algebraic ? 2 : 1, &sub)); 1220 } 1221 PetscFunctionReturn(PETSC_SUCCESS); 1222 } 1223 1224 static PetscErrorCode PCHPDDMAlgebraicAuxiliaryMat_Private(Mat P, IS *is, Mat *sub[], PetscBool block) 1225 { 1226 IS icol[3], irow[2]; 1227 Mat *M, Q; 1228 PetscReal *ptr; 1229 PetscInt *idx, p = 0, n, bs = PetscAbs(P->cmap->bs); 1230 PetscBool flg; 1231 1232 PetscFunctionBegin; 1233 PetscCall(ISCreateStride(PETSC_COMM_SELF, P->cmap->N, 0, 1, icol + 2)); 1234 PetscCall(ISSetBlockSize(icol[2], bs)); 1235 PetscCall(ISSetIdentity(icol[2])); 1236 PetscCall(PetscObjectTypeCompare((PetscObject)P, MATMPISBAIJ, &flg)); 1237 if (flg) { 1238 /* MatCreateSubMatrices() does not handle MATMPISBAIJ properly when iscol != isrow, so convert first to MATMPIBAIJ */ 1239 PetscCall(MatConvert(P, MATMPIBAIJ, MAT_INITIAL_MATRIX, &Q)); 1240 std::swap(P, Q); 1241 } 1242 PetscCall(MatCreateSubMatrices(P, 1, is, icol + 2, MAT_INITIAL_MATRIX, &M)); 1243 if (flg) { 1244 std::swap(P, Q); 1245 PetscCall(MatDestroy(&Q)); 1246 } 1247 PetscCall(ISDestroy(icol + 2)); 1248 PetscCall(ISCreateStride(PETSC_COMM_SELF, M[0]->rmap->N, 0, 1, irow)); 1249 PetscCall(ISSetBlockSize(irow[0], bs)); 1250 PetscCall(ISSetIdentity(irow[0])); 1251 if (!block) { 1252 PetscCall(PetscMalloc2(P->cmap->N, &ptr, P->cmap->N / bs, &idx)); 1253 PetscCall(MatGetColumnNorms(M[0], NORM_INFINITY, ptr)); 1254 /* check for nonzero columns so that M[0] may be expressed in compact form */ 1255 for (n = 0; n < P->cmap->N; n += bs) { 1256 if (std::find_if(ptr + n, ptr + n + bs, [](PetscReal v) { return v > PETSC_MACHINE_EPSILON; }) != ptr + n + bs) idx[p++] = n / bs; 1257 } 1258 PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, p, idx, PETSC_USE_POINTER, icol + 1)); 1259 PetscCall(ISSetInfo(icol[1], IS_SORTED, IS_GLOBAL, PETSC_TRUE, PETSC_TRUE)); 1260 PetscCall(ISEmbed(*is, icol[1], PETSC_FALSE, icol + 2)); 1261 irow[1] = irow[0]; 1262 /* first Mat will be used in PCASM (if it is used as a PC on this level) and as the left-hand side of GenEO */ 1263 icol[0] = is[0]; 1264 PetscCall(MatCreateSubMatrices(M[0], 2, irow, icol, MAT_INITIAL_MATRIX, sub)); 1265 PetscCall(ISDestroy(icol + 1)); 1266 PetscCall(PetscFree2(ptr, idx)); 1267 /* IS used to go back and forth between the augmented and the original local linear system, see eq. (3.4) of [2022b] */ 1268 PetscCall(PetscObjectCompose((PetscObject)(*sub)[0], "_PCHPDDM_Embed", (PetscObject)icol[2])); 1269 /* Mat used in eq. (3.1) of [2022b] */ 1270 PetscCall(PetscObjectCompose((PetscObject)(*sub)[0], "_PCHPDDM_Compact", (PetscObject)(*sub)[1])); 1271 } else { 1272 Mat aux; 1273 PetscCall(MatSetOption(M[0], MAT_SUBMAT_SINGLEIS, PETSC_TRUE)); 1274 /* diagonal block of the overlapping rows */ 1275 PetscCall(MatCreateSubMatrices(M[0], 1, irow, is, MAT_INITIAL_MATRIX, sub)); 1276 PetscCall(MatDuplicate((*sub)[0], MAT_COPY_VALUES, &aux)); 1277 PetscCall(MatSetOption(aux, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE)); 1278 if (bs == 1) { /* scalar case */ 1279 Vec sum[2]; 1280 PetscCall(MatCreateVecs(aux, sum, sum + 1)); 1281 PetscCall(MatGetRowSum(M[0], sum[0])); 1282 PetscCall(MatGetRowSum(aux, sum[1])); 1283 /* off-diagonal block row sum (full rows - diagonal block rows) */ 1284 PetscCall(VecAXPY(sum[0], -1.0, sum[1])); 1285 /* subdomain matrix plus off-diagonal block row sum */ 1286 PetscCall(MatDiagonalSet(aux, sum[0], ADD_VALUES)); 1287 PetscCall(VecDestroy(sum)); 1288 PetscCall(VecDestroy(sum + 1)); 1289 } else { /* vectorial case */ 1290 /* TODO: missing MatGetValuesBlocked(), so the code below is */ 1291 /* an extension of the scalar case for when bs > 1, but it could */ 1292 /* be more efficient by avoiding all these MatMatMult() */ 1293 Mat sum[2], ones; 1294 PetscScalar *ptr; 1295 PetscCall(PetscCalloc1(M[0]->cmap->n * bs, &ptr)); 1296 PetscCall(MatCreateDense(PETSC_COMM_SELF, M[0]->cmap->n, bs, M[0]->cmap->n, bs, ptr, &ones)); 1297 for (n = 0; n < M[0]->cmap->n; n += bs) { 1298 for (p = 0; p < bs; ++p) ptr[n + p * (M[0]->cmap->n + 1)] = 1.0; 1299 } 1300 PetscCall(MatMatMult(M[0], ones, MAT_INITIAL_MATRIX, PETSC_DEFAULT, sum)); 1301 PetscCall(MatDestroy(&ones)); 1302 PetscCall(MatCreateDense(PETSC_COMM_SELF, aux->cmap->n, bs, aux->cmap->n, bs, ptr, &ones)); 1303 PetscCall(MatDenseSetLDA(ones, M[0]->cmap->n)); 1304 PetscCall(MatMatMult(aux, ones, MAT_INITIAL_MATRIX, PETSC_DEFAULT, sum + 1)); 1305 PetscCall(MatDestroy(&ones)); 1306 PetscCall(PetscFree(ptr)); 1307 /* off-diagonal block row sum (full rows - diagonal block rows) */ 1308 PetscCall(MatAXPY(sum[0], -1.0, sum[1], SAME_NONZERO_PATTERN)); 1309 PetscCall(MatDestroy(sum + 1)); 1310 /* re-order values to be consistent with MatSetValuesBlocked() */ 1311 /* equivalent to MatTranspose() which does not truly handle */ 1312 /* MAT_INPLACE_MATRIX in the rectangular case, as it calls PetscMalloc() */ 1313 PetscCall(MatDenseGetArrayWrite(sum[0], &ptr)); 1314 HPDDM::Wrapper<PetscScalar>::imatcopy<'T'>(bs, sum[0]->rmap->n, ptr, sum[0]->rmap->n, bs); 1315 /* subdomain matrix plus off-diagonal block row sum */ 1316 for (n = 0; n < aux->cmap->n / bs; ++n) PetscCall(MatSetValuesBlocked(aux, 1, &n, 1, &n, ptr + n * bs * bs, ADD_VALUES)); 1317 PetscCall(MatAssemblyBegin(aux, MAT_FINAL_ASSEMBLY)); 1318 PetscCall(MatAssemblyEnd(aux, MAT_FINAL_ASSEMBLY)); 1319 PetscCall(MatDenseRestoreArrayWrite(sum[0], &ptr)); 1320 PetscCall(MatDestroy(sum)); 1321 } 1322 PetscCall(MatSetOption(aux, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE)); 1323 /* left-hand side of GenEO, with the same sparsity pattern as PCASM subdomain solvers */ 1324 PetscCall(PetscObjectCompose((PetscObject)(*sub)[0], "_PCHPDDM_Neumann_Mat", (PetscObject)aux)); 1325 } 1326 PetscCall(ISDestroy(irow)); 1327 PetscCall(MatDestroySubMatrices(1, &M)); 1328 PetscFunctionReturn(PETSC_SUCCESS); 1329 } 1330 1331 static PetscErrorCode PCApply_Nest(PC pc, Vec x, Vec y) 1332 { 1333 Mat A; 1334 MatSolverType type; 1335 IS is[2]; 1336 PetscBool flg; 1337 std::pair<PC, Vec[2]> *p; 1338 1339 PetscFunctionBegin; 1340 PetscCall(PCShellGetContext(pc, &p)); 1341 PetscCall(PCGetOperators(p->first, &A, nullptr)); 1342 PetscCall(MatNestGetISs(A, is, nullptr)); 1343 PetscCall(PCFactorGetMatSolverType(p->first, &type)); 1344 PetscCall(PCFactorGetMatrix(p->first, &A)); 1345 PetscCall(PetscStrcmp(type, MATSOLVERMUMPS, &flg)); 1346 if (flg && A->schur) { 1347 #if PetscDefined(HAVE_MUMPS) 1348 PetscCall(MatMumpsSetIcntl(A, 26, 1)); /* reduction/condensation phase followed by Schur complement solve */ 1349 #endif 1350 } 1351 PetscCall(VecISCopy(p->second[0], is[1], SCATTER_FORWARD, x)); /* assign the RHS associated to the Schur complement */ 1352 PetscCall(PCApply(p->first, p->second[0], p->second[1])); 1353 PetscCall(VecISCopy(p->second[1], is[1], SCATTER_REVERSE, y)); /* retrieve the partial solution associated to the Schur complement */ 1354 if (flg) { 1355 #if PetscDefined(HAVE_MUMPS) 1356 PetscCall(MatMumpsSetIcntl(A, 26, -1)); /* default ICNTL(26) value in PETSc */ 1357 #endif 1358 } 1359 PetscFunctionReturn(PETSC_SUCCESS); 1360 } 1361 1362 static PetscErrorCode PCView_Nest(PC pc, PetscViewer viewer) 1363 { 1364 std::pair<PC, Vec[2]> *p; 1365 1366 PetscFunctionBegin; 1367 PetscCall(PCShellGetContext(pc, &p)); 1368 PetscCall(PCView(p->first, viewer)); 1369 PetscFunctionReturn(PETSC_SUCCESS); 1370 } 1371 1372 static PetscErrorCode PCDestroy_Nest(PC pc) 1373 { 1374 std::pair<PC, Vec[2]> *p; 1375 1376 PetscFunctionBegin; 1377 PetscCall(PCShellGetContext(pc, &p)); 1378 PetscCall(VecDestroy(p->second)); 1379 PetscCall(VecDestroy(p->second + 1)); 1380 PetscCall(PCDestroy(&p->first)); 1381 PetscCall(PetscFree(p)); 1382 PetscFunctionReturn(PETSC_SUCCESS); 1383 } 1384 1385 template <bool T = false> 1386 static PetscErrorCode MatMult_Schur(Mat A, Vec x, Vec y) 1387 { 1388 std::tuple<Mat, PetscSF, Vec[2]> *ctx; 1389 1390 PetscFunctionBegin; 1391 PetscCall(MatShellGetContext(A, &ctx)); 1392 PetscCall(VecScatterBegin(std::get<1>(*ctx), x, std::get<2>(*ctx)[0], INSERT_VALUES, SCATTER_FORWARD)); /* local Vec with overlap */ 1393 PetscCall(VecScatterEnd(std::get<1>(*ctx), x, std::get<2>(*ctx)[0], INSERT_VALUES, SCATTER_FORWARD)); 1394 if (!T) PetscCall(MatMult(std::get<0>(*ctx), std::get<2>(*ctx)[0], std::get<2>(*ctx)[1])); /* local Schur complement */ 1395 else PetscCall(MatMultTranspose(std::get<0>(*ctx), std::get<2>(*ctx)[0], std::get<2>(*ctx)[1])); 1396 PetscCall(VecSet(y, 0.0)); 1397 PetscCall(VecScatterBegin(std::get<1>(*ctx), std::get<2>(*ctx)[1], y, ADD_VALUES, SCATTER_REVERSE)); /* global Vec with summed up contributions on the overlap */ 1398 PetscCall(VecScatterEnd(std::get<1>(*ctx), std::get<2>(*ctx)[1], y, ADD_VALUES, SCATTER_REVERSE)); 1399 PetscFunctionReturn(PETSC_SUCCESS); 1400 } 1401 1402 static PetscErrorCode MatDestroy_Schur(Mat A) 1403 { 1404 std::tuple<Mat, PetscSF, Vec[2]> *ctx; 1405 1406 PetscFunctionBegin; 1407 PetscCall(MatShellGetContext(A, &ctx)); 1408 PetscCall(VecDestroy(std::get<2>(*ctx))); 1409 PetscCall(VecDestroy(std::get<2>(*ctx) + 1)); 1410 PetscCall(PetscFree(ctx)); 1411 PetscFunctionReturn(PETSC_SUCCESS); 1412 } 1413 1414 static PetscErrorCode MatMult_SchurCorrection(Mat A, Vec x, Vec y) 1415 { 1416 PC pc; 1417 std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *ctx; 1418 1419 PetscFunctionBegin; 1420 PetscCall(MatShellGetContext(A, &ctx)); 1421 pc = ((PC_HPDDM *)std::get<0>(*ctx)[0]->data)->levels[0]->ksp->pc; 1422 if (std::get<2>(*ctx) == PC_LEFT || std::get<2>(*ctx) == PC_SIDE_DEFAULT) { /* Q_0 is the coarse correction associated to the A00 block from PCFIELDSPLIT */ 1423 PetscCall(MatMult(std::get<1>(*ctx)[0], x, std::get<3>(*ctx)[1])); /* A_01 x */ 1424 PetscCall(PCHPDDMDeflate_Private(pc, std::get<3>(*ctx)[1], std::get<3>(*ctx)[1])); /* Q_0 A_01 x */ 1425 PetscCall(MatMult(std::get<1>(*ctx)[1], std::get<3>(*ctx)[1], std::get<3>(*ctx)[0])); /* A_10 Q_0 A_01 x */ 1426 PetscCall(PCApply(std::get<0>(*ctx)[1], std::get<3>(*ctx)[0], y)); /* y = M_S^-1 A_10 Q_0 A_01 x */ 1427 } else { 1428 PetscCall(PCApply(std::get<0>(*ctx)[1], x, std::get<3>(*ctx)[0])); /* M_S^-1 x */ 1429 PetscCall(MatMult(std::get<1>(*ctx)[0], std::get<3>(*ctx)[0], std::get<3>(*ctx)[1])); /* A_01 M_S^-1 x */ 1430 PetscCall(PCHPDDMDeflate_Private(pc, std::get<3>(*ctx)[1], std::get<3>(*ctx)[1])); /* Q_0 A_01 M_S^-1 x */ 1431 PetscCall(MatMult(std::get<1>(*ctx)[1], std::get<3>(*ctx)[1], y)); /* y = A_10 Q_0 A_01 M_S^-1 x */ 1432 } 1433 PetscCall(VecAXPY(y, -1.0, x)); /* y -= x, preconditioned eq. (24) of https://hal.science/hal-02343808v6/document (with a sign flip) */ 1434 PetscFunctionReturn(PETSC_SUCCESS); 1435 } 1436 1437 static PetscErrorCode MatView_SchurCorrection(Mat A, PetscViewer viewer) 1438 { 1439 PetscBool ascii; 1440 std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *ctx; 1441 1442 PetscFunctionBegin; 1443 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &ascii)); 1444 if (ascii) { 1445 PetscCall(MatShellGetContext(A, &ctx)); 1446 PetscCall(PetscViewerASCIIPrintf(viewer, "action of %s\n", std::get<2>(*ctx) == PC_LEFT || std::get<2>(*ctx) == PC_SIDE_DEFAULT ? "(I - M_S^-1 A_10 Q_0 A_01)" : "(I - A_10 Q_0 A_01 M_S^-1)")); 1447 PetscCall(PCView(std::get<0>(*ctx)[1], viewer)); /* no need to PCView(Q_0) since it will be done by PCFIELDSPLIT */ 1448 } 1449 PetscFunctionReturn(PETSC_SUCCESS); 1450 } 1451 1452 static PetscErrorCode MatDestroy_SchurCorrection(Mat A) 1453 { 1454 std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *ctx; 1455 1456 PetscFunctionBegin; 1457 PetscCall(MatShellGetContext(A, &ctx)); 1458 PetscCall(VecDestroy(std::get<3>(*ctx))); 1459 PetscCall(VecDestroy(std::get<3>(*ctx) + 1)); 1460 PetscCall(VecDestroy(std::get<3>(*ctx) + 2)); 1461 PetscCall(PCDestroy(std::get<0>(*ctx) + 1)); 1462 PetscCall(PetscFree(ctx)); 1463 PetscFunctionReturn(PETSC_SUCCESS); 1464 } 1465 1466 static PetscErrorCode KSPPreSolve_SchurCorrection(KSP, Vec b, Vec, void *context) 1467 { 1468 std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *ctx = reinterpret_cast<std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *>(context); 1469 1470 PetscFunctionBegin; 1471 if (std::get<2>(*ctx) == PC_LEFT || std::get<2>(*ctx) == PC_SIDE_DEFAULT) { 1472 PetscCall(PCApply(std::get<0>(*ctx)[1], b, std::get<3>(*ctx)[2])); 1473 std::swap(*b, *std::get<3>(*ctx)[2]); /* replace b by M^-1 b, but need to keep a copy of the original RHS, so swap it with the work Vec */ 1474 } 1475 PetscFunctionReturn(PETSC_SUCCESS); 1476 } 1477 1478 static PetscErrorCode KSPPostSolve_SchurCorrection(KSP, Vec b, Vec x, void *context) 1479 { 1480 std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *ctx = reinterpret_cast<std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *>(context); 1481 1482 PetscFunctionBegin; 1483 if (std::get<2>(*ctx) == PC_LEFT || std::get<2>(*ctx) == PC_SIDE_DEFAULT) std::swap(*b, *std::get<3>(*ctx)[2]); /* put back the original RHS where it belongs */ 1484 else { 1485 PetscCall(PCApply(std::get<0>(*ctx)[1], x, std::get<3>(*ctx)[2])); 1486 PetscCall(VecCopy(std::get<3>(*ctx)[2], x)); /* replace x by M^-1 x */ 1487 } 1488 PetscFunctionReturn(PETSC_SUCCESS); 1489 } 1490 1491 static PetscErrorCode MatMult_Harmonic(Mat, Vec, Vec); 1492 static PetscErrorCode MatMultTranspose_Harmonic(Mat, Vec, Vec); 1493 static PetscErrorCode MatProduct_AB_Harmonic(Mat, Mat, Mat, void *); 1494 static PetscErrorCode MatDestroy_Harmonic(Mat); 1495 1496 static PetscErrorCode PCSetUp_HPDDM(PC pc) 1497 { 1498 PC_HPDDM *data = (PC_HPDDM *)pc->data; 1499 PC inner; 1500 KSP *ksp; 1501 Mat *sub, A, P, N, C = nullptr, uaux = nullptr, weighted, subA[2], S; 1502 Vec xin, v; 1503 std::vector<Vec> initial; 1504 IS is[1], loc, uis = data->is, unsorted = nullptr; 1505 ISLocalToGlobalMapping l2g; 1506 char prefix[256]; 1507 const char *pcpre; 1508 const PetscScalar *const *ev; 1509 PetscInt n, requested = data->N, reused = 0, overlap = -1; 1510 MatStructure structure = UNKNOWN_NONZERO_PATTERN; 1511 PetscBool subdomains = PETSC_FALSE, flg = PETSC_FALSE, ismatis, swap = PETSC_FALSE, algebraic = PETSC_FALSE, block = PETSC_FALSE; 1512 DM dm; 1513 std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *ctx = nullptr; 1514 #if PetscDefined(USE_DEBUG) 1515 IS dis = nullptr; 1516 Mat daux = nullptr; 1517 #endif 1518 1519 PetscFunctionBegin; 1520 PetscCheck(data->levels && data->levels[0], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Not a single level allocated"); 1521 PetscCall(PCGetOptionsPrefix(pc, &pcpre)); 1522 PetscCall(PCGetOperators(pc, &A, &P)); 1523 if (!data->levels[0]->ksp) { 1524 PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &data->levels[0]->ksp)); 1525 PetscCall(KSPSetNestLevel(data->levels[0]->ksp, pc->kspnestlevel)); 1526 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_%s_", pcpre ? pcpre : "", data->N > 1 ? "levels_1" : "coarse")); 1527 PetscCall(KSPSetOptionsPrefix(data->levels[0]->ksp, prefix)); 1528 PetscCall(KSPSetType(data->levels[0]->ksp, KSPPREONLY)); 1529 } else if (data->levels[0]->ksp->pc && data->levels[0]->ksp->pc->setupcalled == 1 && data->levels[0]->ksp->pc->reusepreconditioner) { 1530 /* if the fine-level PCSHELL exists, its setup has succeeded, and one wants to reuse it, */ 1531 /* then just propagate the appropriate flag to the coarser levels */ 1532 for (n = 0; n < PETSC_PCHPDDM_MAXLEVELS && data->levels[n]; ++n) { 1533 /* the following KSP and PC may be NULL for some processes, hence the check */ 1534 if (data->levels[n]->ksp) PetscCall(KSPSetReusePreconditioner(data->levels[n]->ksp, PETSC_TRUE)); 1535 if (data->levels[n]->pc) PetscCall(PCSetReusePreconditioner(data->levels[n]->pc, PETSC_TRUE)); 1536 } 1537 /* early bail out because there is nothing to do */ 1538 PetscFunctionReturn(PETSC_SUCCESS); 1539 } else { 1540 /* reset coarser levels */ 1541 for (n = 1; n < PETSC_PCHPDDM_MAXLEVELS && data->levels[n]; ++n) { 1542 if (data->levels[n]->ksp && data->levels[n]->ksp->pc && data->levels[n]->ksp->pc->setupcalled == 1 && data->levels[n]->ksp->pc->reusepreconditioner && n < data->N) { 1543 reused = data->N - n; 1544 break; 1545 } 1546 PetscCall(KSPDestroy(&data->levels[n]->ksp)); 1547 PetscCall(PCDestroy(&data->levels[n]->pc)); 1548 } 1549 /* check if some coarser levels are being reused */ 1550 PetscCall(MPIU_Allreduce(MPI_IN_PLACE, &reused, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)pc))); 1551 const int *addr = data->levels[0]->P ? data->levels[0]->P->getAddrLocal() : &HPDDM::i__0; 1552 1553 if (addr != &HPDDM::i__0 && reused != data->N - 1) { 1554 /* reuse previously computed eigenvectors */ 1555 ev = data->levels[0]->P->getVectors(); 1556 if (ev) { 1557 initial.reserve(*addr); 1558 PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, data->levels[0]->P->getDof(), ev[0], &xin)); 1559 for (n = 0; n < *addr; ++n) { 1560 PetscCall(VecDuplicate(xin, &v)); 1561 PetscCall(VecPlaceArray(xin, ev[n])); 1562 PetscCall(VecCopy(xin, v)); 1563 initial.emplace_back(v); 1564 PetscCall(VecResetArray(xin)); 1565 } 1566 PetscCall(VecDestroy(&xin)); 1567 } 1568 } 1569 } 1570 data->N -= reused; 1571 PetscCall(KSPSetOperators(data->levels[0]->ksp, A, P)); 1572 1573 PetscCall(PetscObjectTypeCompare((PetscObject)P, MATIS, &ismatis)); 1574 if (!data->is && !ismatis) { 1575 PetscErrorCode (*create)(DM, IS *, Mat *, PetscErrorCode (**)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void *), void **) = nullptr; 1576 PetscErrorCode (*usetup)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void *) = nullptr; 1577 void *uctx = nullptr; 1578 1579 /* first see if we can get the data from the DM */ 1580 PetscCall(MatGetDM(P, &dm)); 1581 if (!dm) PetscCall(MatGetDM(A, &dm)); 1582 if (!dm) PetscCall(PCGetDM(pc, &dm)); 1583 if (dm) { /* this is the hook for DMPLEX for which the auxiliary Mat is the local Neumann matrix */ 1584 PetscCall(PetscObjectQueryFunction((PetscObject)dm, "DMCreateNeumannOverlap_C", &create)); 1585 if (create) { 1586 PetscCall((*create)(dm, &uis, &uaux, &usetup, &uctx)); 1587 if (data->Neumann == PETSC_BOOL3_UNKNOWN) data->Neumann = PETSC_BOOL3_TRUE; /* set the value only if it was not already provided by the user */ 1588 } 1589 } 1590 if (!create) { 1591 if (!uis) { 1592 PetscCall(PetscObjectQuery((PetscObject)pc, "_PCHPDDM_Neumann_IS", (PetscObject *)&uis)); 1593 PetscCall(PetscObjectReference((PetscObject)uis)); 1594 } 1595 if (!uaux) { 1596 PetscCall(PetscObjectQuery((PetscObject)pc, "_PCHPDDM_Neumann_Mat", (PetscObject *)&uaux)); 1597 PetscCall(PetscObjectReference((PetscObject)uaux)); 1598 } 1599 /* look inside the Pmat instead of the PC, needed for MatSchurComplementComputeExplicitOperator() */ 1600 if (!uis) { 1601 PetscCall(PetscObjectQuery((PetscObject)P, "_PCHPDDM_Neumann_IS", (PetscObject *)&uis)); 1602 PetscCall(PetscObjectReference((PetscObject)uis)); 1603 } 1604 if (!uaux) { 1605 PetscCall(PetscObjectQuery((PetscObject)P, "_PCHPDDM_Neumann_Mat", (PetscObject *)&uaux)); 1606 PetscCall(PetscObjectReference((PetscObject)uaux)); 1607 } 1608 } 1609 PetscCall(PCHPDDMSetAuxiliaryMat(pc, uis, uaux, usetup, uctx)); 1610 PetscCall(MatDestroy(&uaux)); 1611 PetscCall(ISDestroy(&uis)); 1612 } 1613 1614 if (!ismatis) { 1615 PetscCall(PCHPDDMSetUpNeumannOverlap_Private(pc)); 1616 PetscCall(PetscOptionsGetBool(nullptr, pcpre, "-pc_hpddm_block_splitting", &block, nullptr)); 1617 PetscCall(PetscOptionsGetInt(nullptr, pcpre, "-pc_hpddm_harmonic_overlap", &overlap, nullptr)); 1618 PetscCall(PetscObjectTypeCompare((PetscObject)P, MATSCHURCOMPLEMENT, &flg)); 1619 if (data->is || (data->N > 1 && flg)) { 1620 if (block || overlap != -1) { 1621 PetscCall(ISDestroy(&data->is)); 1622 PetscCall(MatDestroy(&data->aux)); 1623 } else if (data->N > 1 && flg) { 1624 PCHPDDMSchurPreType type = PC_HPDDM_SCHUR_PRE_GENEO; 1625 1626 PetscCall(PetscOptionsGetEnum(nullptr, pcpre, "-pc_hpddm_schur_precondition", PCHPDDMSchurPreTypes, (PetscEnum *)&type, &flg)); 1627 if (type == PC_HPDDM_SCHUR_PRE_LEAST_SQUARES) { 1628 PetscCall(ISDestroy(&data->is)); /* destroy any previously user-set objects since they will be set automatically */ 1629 PetscCall(MatDestroy(&data->aux)); 1630 } else if (type == PC_HPDDM_SCHUR_PRE_GENEO) { 1631 PetscContainer container = nullptr; 1632 1633 PetscCall(PetscObjectQuery((PetscObject)pc, "_PCHPDDM_Schur", (PetscObject *)&container)); 1634 if (!container) { /* first call to PCSetUp() on the PC associated to the Schur complement */ 1635 PC_HPDDM *data_00; 1636 KSP ksp, inner_ksp; 1637 PC pc_00; 1638 Mat A11 = nullptr; 1639 Vec d = nullptr; 1640 char *prefix; 1641 1642 PetscCall(MatSchurComplementGetKSP(P, &ksp)); 1643 PetscCall(KSPGetPC(ksp, &pc_00)); 1644 PetscCall(PetscObjectTypeCompare((PetscObject)pc_00, PCHPDDM, &flg)); 1645 PetscCheck(flg, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-%spc_hpddm_schur_precondition %s and -%spc_type %s (!= %s)", pcpre ? pcpre : "", PCHPDDMSchurPreTypes[type], ((PetscObject)pc_00)->prefix ? ((PetscObject)pc_00)->prefix : "", 1646 ((PetscObject)pc_00)->type_name, PCHPDDM); 1647 data_00 = (PC_HPDDM *)pc_00->data; 1648 PetscCheck(data_00->N == 2, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-%spc_hpddm_schur_precondition %s and %" PetscInt_FMT " level%s instead of 2 for the A00 block -%s", pcpre ? pcpre : "", PCHPDDMSchurPreTypes[type], 1649 data_00->N, data_00->N > 1 ? "s" : "", ((PetscObject)pc_00)->prefix); 1650 PetscCall(PetscObjectTypeCompare((PetscObject)data_00->levels[0]->pc, PCASM, &flg)); 1651 PetscCheck(flg, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-%spc_hpddm_schur_precondition %s and -%spc_type %s (!= %s)", pcpre ? pcpre : "", PCHPDDMSchurPreTypes[type], ((PetscObject)data_00->levels[0]->pc)->prefix, 1652 ((PetscObject)data_00->levels[0]->pc)->type_name, PCASM); 1653 if (PetscDefined(USE_DEBUG) || !data->is) { 1654 Mat A01, A10, B = nullptr, C = nullptr, *sub; 1655 1656 PetscCall(MatSchurComplementGetSubMatrices(P, &A, nullptr, &A01, &A10, nullptr)); 1657 PetscCall(PetscObjectTypeCompare((PetscObject)A10, MATTRANSPOSEVIRTUAL, &flg)); 1658 if (flg) { 1659 PetscCall(MatTransposeGetMat(A10, &C)); 1660 PetscCall(MatTranspose(C, MAT_INITIAL_MATRIX, &B)); 1661 } else { 1662 PetscCall(PetscObjectTypeCompare((PetscObject)A10, MATHERMITIANTRANSPOSEVIRTUAL, &flg)); 1663 if (flg) { 1664 PetscCall(MatHermitianTransposeGetMat(A10, &C)); 1665 PetscCall(MatHermitianTranspose(C, MAT_INITIAL_MATRIX, &B)); 1666 } 1667 } 1668 if (!B) { 1669 B = A10; 1670 PetscCall(PetscObjectReference((PetscObject)B)); 1671 } else if (!data->is) { 1672 PetscCall(PetscObjectTypeCompareAny((PetscObject)A01, &flg, MATTRANSPOSEVIRTUAL, MATHERMITIANTRANSPOSEVIRTUAL, "")); 1673 if (!flg) C = A01; 1674 } 1675 PetscCall(ISCreateStride(PETSC_COMM_SELF, B->rmap->N, 0, 1, &uis)); 1676 PetscCall(ISSetIdentity(uis)); 1677 if (!data->is) { 1678 if (C) PetscCall(PetscObjectReference((PetscObject)C)); 1679 else PetscCall(MatTranspose(B, MAT_INITIAL_MATRIX, &C)); 1680 PetscCall(ISDuplicate(data_00->is, is)); 1681 PetscCall(MatIncreaseOverlap(A, 1, is, 1)); 1682 PetscCall(MatSetOption(C, MAT_SUBMAT_SINGLEIS, PETSC_TRUE)); 1683 PetscCall(MatCreateSubMatrices(C, 1, is, &uis, MAT_INITIAL_MATRIX, &sub)); 1684 PetscCall(MatDestroy(&C)); 1685 PetscCall(MatTranspose(sub[0], MAT_INITIAL_MATRIX, &C)); 1686 PetscCall(MatDestroySubMatrices(1, &sub)); 1687 PetscCall(MatFindNonzeroRows(C, &data->is)); 1688 PetscCall(MatDestroy(&C)); 1689 PetscCall(ISDestroy(is)); 1690 } 1691 if (PetscDefined(USE_DEBUG)) { 1692 PetscCall(PCHPDDMCheckSymmetry_Private(pc, A01, A10)); 1693 PetscCall(MatCreateSubMatrices(B, 1, &uis, &data_00->is, MAT_INITIAL_MATRIX, &sub)); /* expensive check since all processes fetch all rows (but only some columns) of the constraint matrix */ 1694 PetscCall(ISDestroy(&uis)); 1695 PetscCall(ISDuplicate(data->is, &uis)); 1696 PetscCall(ISSort(uis)); 1697 PetscCall(ISComplement(uis, 0, B->rmap->N, is)); 1698 PetscCall(MatDuplicate(sub[0], MAT_COPY_VALUES, &C)); 1699 PetscCall(MatZeroRowsIS(C, is[0], 0.0, nullptr, nullptr)); 1700 PetscCall(ISDestroy(is)); 1701 PetscCall(MatMultEqual(sub[0], C, 20, &flg)); 1702 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "The image of A_10 (R_i^p)^T from the local primal (e.g., velocity) space to the full dual (e.g., pressure) space is not restricted to the local dual space: A_10 (R_i^p)^T != R_i^d (R_i^d)^T A_10 (R_i^p)^T"); /* cf. eq. (9) of https://hal.science/hal-02343808v6/document */ 1703 PetscCall(MatDestroy(&C)); 1704 PetscCall(MatDestroySubMatrices(1, &sub)); 1705 } 1706 PetscCall(ISDestroy(&uis)); 1707 PetscCall(MatDestroy(&B)); 1708 } 1709 PetscCall(MatSchurComplementGetSubMatrices(P, nullptr, nullptr, nullptr, nullptr, &A11)); 1710 flg = PETSC_FALSE; 1711 if (!data->aux) { 1712 Mat D; 1713 1714 PetscCall(MatCreateVecs(A11, &d, nullptr)); 1715 PetscCall(MatGetDiagonal(A11, d)); 1716 PetscCall(PetscObjectTypeCompareAny((PetscObject)A11, &flg, MATDIAGONAL, MATCONSTANTDIAGONAL, "")); 1717 if (!flg) { 1718 PetscCall(MatCreateDiagonal(d, &D)); 1719 PetscCall(MatMultEqual(A11, D, 20, &flg)); 1720 PetscCall(MatDestroy(&D)); 1721 } 1722 if (flg) PetscCall(PetscInfo(pc, "A11 block is likely diagonal so the PC will build an auxiliary Mat (which was not initially provided by the user)\n")); 1723 } 1724 if (data->Neumann != PETSC_BOOL3_TRUE && !flg && A11) { 1725 PetscReal norm; 1726 1727 PetscCall(MatNorm(A11, NORM_INFINITY, &norm)); 1728 PetscCheck(norm < PETSC_MACHINE_EPSILON * static_cast<PetscReal>(10.0), PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-%spc_hpddm_schur_precondition geneo and -%spc_hpddm_has_neumann != true with a nonzero or non-diagonal A11 block", pcpre ? pcpre : "", pcpre ? pcpre : ""); 1729 PetscCall(PetscInfo(pc, "A11 block is likely zero so the PC will build an auxiliary Mat (which was%s initially provided by the user)\n", data->aux ? "" : " not")); 1730 PetscCall(MatDestroy(&data->aux)); 1731 flg = PETSC_TRUE; 1732 } 1733 if (!data->aux) { /* if A11 is near zero, e.g., Stokes equation, or diagonal, build an auxiliary (Neumann) Mat which is a (possibly slightly shifted) diagonal weighted by the inverse of the multiplicity */ 1734 PetscSF scatter; 1735 const PetscScalar *read; 1736 PetscScalar *write, *diagonal = nullptr; 1737 1738 PetscCall(MatDestroy(&data->aux)); 1739 PetscCall(ISGetLocalSize(data->is, &n)); 1740 PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)P), n, PETSC_DECIDE, &xin)); 1741 PetscCall(VecDuplicate(xin, &v)); 1742 PetscCall(VecScatterCreate(xin, data->is, v, nullptr, &scatter)); 1743 PetscCall(VecSet(v, 1.0)); 1744 PetscCall(VecSet(xin, 1.0)); 1745 PetscCall(VecScatterBegin(scatter, v, xin, ADD_VALUES, SCATTER_REVERSE)); 1746 PetscCall(VecScatterEnd(scatter, v, xin, ADD_VALUES, SCATTER_REVERSE)); /* v has the multiplicity of all unknowns on the overlap */ 1747 PetscCall(PetscSFDestroy(&scatter)); 1748 if (d) { 1749 PetscCall(VecScatterCreate(d, data->is, v, nullptr, &scatter)); 1750 PetscCall(VecScatterBegin(scatter, d, v, INSERT_VALUES, SCATTER_FORWARD)); 1751 PetscCall(VecScatterEnd(scatter, d, v, INSERT_VALUES, SCATTER_FORWARD)); 1752 PetscCall(PetscSFDestroy(&scatter)); 1753 PetscCall(VecDestroy(&d)); 1754 PetscCall(PetscMalloc1(n, &diagonal)); 1755 PetscCall(VecGetArrayRead(v, &read)); 1756 PetscCallCXX(std::copy_n(read, n, diagonal)); 1757 PetscCall(VecRestoreArrayRead(v, &read)); 1758 } 1759 PetscCall(VecDestroy(&v)); 1760 PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &v)); 1761 PetscCall(VecGetArrayRead(xin, &read)); 1762 PetscCall(VecGetArrayWrite(v, &write)); 1763 for (PetscInt i = 0; i < n; ++i) write[i] = (!diagonal || std::abs(diagonal[i]) < PETSC_MACHINE_EPSILON) ? PETSC_SMALL / (static_cast<PetscReal>(1000.0) * read[i]) : diagonal[i] / read[i]; 1764 PetscCall(PetscFree(diagonal)); 1765 PetscCall(VecRestoreArrayRead(xin, &read)); 1766 PetscCall(VecRestoreArrayWrite(v, &write)); 1767 PetscCall(VecDestroy(&xin)); 1768 PetscCall(MatCreateDiagonal(v, &data->aux)); 1769 PetscCall(VecDestroy(&v)); 1770 } 1771 uis = data->is; 1772 uaux = data->aux; 1773 PetscCall(PetscObjectReference((PetscObject)uis)); 1774 PetscCall(PetscObjectReference((PetscObject)uaux)); 1775 PetscCall(PetscStrallocpy(pcpre, &prefix)); 1776 PetscCall(PCSetOptionsPrefix(pc, nullptr)); 1777 PetscCall(PCSetType(pc, PCKSP)); /* replace the PC associated to the Schur complement by PCKSP */ 1778 PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &inner_ksp)); /* new KSP that will be attached to the previously set PC */ 1779 pc->ops->setup = PCSetUp_KSP; 1780 PetscCall(PetscObjectGetTabLevel((PetscObject)pc, &n)); 1781 PetscCall(PetscObjectSetTabLevel((PetscObject)inner_ksp, n + 2)); 1782 PetscCall(KSPSetOperators(inner_ksp, pc->mat, pc->pmat)); 1783 PetscCall(KSPSetOptionsPrefix(inner_ksp, std::string(std::string(prefix) + "pc_hpddm_").c_str())); 1784 PetscCall(KSPSetSkipPCSetFromOptions(inner_ksp, PETSC_TRUE)); 1785 PetscCall(KSPSetFromOptions(inner_ksp)); 1786 PetscCall(KSPGetPC(inner_ksp, &inner)); 1787 PetscCall(PCSetOptionsPrefix(inner, nullptr)); 1788 PetscCall(PCSetType(inner, PCNONE)); /* no preconditioner since the action of M^-1 A or A M^-1 will be computed by the Amat */ 1789 PetscCall(PCKSPSetKSP(pc, inner_ksp)); 1790 PetscCall(PetscContainerCreate(PetscObjectComm((PetscObject)pc), &container)); 1791 PetscCall(PetscNew(&ctx)); /* context to pass data around for the inner-most PC, which will be a proper PCHPDDM */ 1792 PetscCall(PetscContainerSetPointer(container, ctx)); 1793 std::get<0>(*ctx)[0] = pc_00; /* for coarse correction on the primal (e.g., velocity) space */ 1794 PetscCall(PCCreate(PetscObjectComm((PetscObject)pc), &std::get<0>(*ctx)[1])); 1795 PetscCall(PCSetOptionsPrefix(std::get<0>(*ctx)[1], prefix)); 1796 PetscCall(PetscFree(prefix)); 1797 PetscCall(PCSetOperators(std::get<0>(*ctx)[1], pc->mat, pc->pmat)); 1798 PetscCall(PCSetType(std::get<0>(*ctx)[1], PCHPDDM)); 1799 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 */ 1800 if (flg) static_cast<PC_HPDDM *>(std::get<0>(*ctx)[1]->data)->Neumann = PETSC_BOOL3_TRUE; 1801 PetscCall(PCSetFromOptions(std::get<0>(*ctx)[1])); 1802 PetscCall(PetscObjectDereference((PetscObject)uis)); 1803 PetscCall(PetscObjectDereference((PetscObject)uaux)); 1804 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 */ 1805 PetscCall(MatShellSetOperation(S, MATOP_MULT, (void (*)(void))MatMult_SchurCorrection)); 1806 PetscCall(MatShellSetOperation(S, MATOP_VIEW, (void (*)(void))MatView_SchurCorrection)); 1807 PetscCall(MatShellSetOperation(S, MATOP_DESTROY, (void (*)(void))MatDestroy_SchurCorrection)); 1808 PetscCall(KSPGetPCSide(inner_ksp, &(std::get<2>(*ctx)))); 1809 if (std::get<2>(*ctx) == PC_LEFT || std::get<2>(*ctx) == PC_SIDE_DEFAULT) { 1810 PetscCall(KSPSetPreSolve(inner_ksp, KSPPreSolve_SchurCorrection, ctx)); 1811 } else { /* no support for PC_SYMMETRIC */ 1812 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]); 1813 } 1814 PetscCall(KSPSetPostSolve(inner_ksp, KSPPostSolve_SchurCorrection, ctx)); 1815 PetscCall(PetscObjectCompose((PetscObject)(std::get<0>(*ctx)[1]), "_PCHPDDM_Schur", (PetscObject)container)); 1816 PetscCall(PetscObjectDereference((PetscObject)container)); 1817 PetscCall(PCSetUp(std::get<0>(*ctx)[1])); 1818 PetscCall(KSPSetOperators(inner_ksp, S, S)); 1819 PetscCall(MatCreateVecs(std::get<1>(*ctx)[0], std::get<3>(*ctx), std::get<3>(*ctx) + 1)); 1820 PetscCall(VecDuplicate(std::get<3>(*ctx)[0], std::get<3>(*ctx) + 2)); 1821 PetscCall(PetscObjectDereference((PetscObject)inner_ksp)); 1822 PetscCall(PetscObjectDereference((PetscObject)S)); 1823 for (std::vector<Vec>::iterator it = initial.begin(); it != initial.end(); ++it) PetscCall(VecDestroy(&*it)); 1824 PetscFunctionReturn(PETSC_SUCCESS); 1825 } else { /* second call to PCSetUp() on the PC associated to the Schur complement, retrieve previously set context */ 1826 PetscCall(PetscContainerGetPointer(container, (void **)&ctx)); 1827 } 1828 } 1829 } 1830 } 1831 if (!data->is && data->N > 1) { 1832 char type[256] = {}; /* same size as in src/ksp/pc/interface/pcset.c */ 1833 PetscCall(PetscObjectTypeCompareAny((PetscObject)P, &flg, MATNORMAL, MATNORMALHERMITIAN, "")); 1834 if (flg || (A->rmap->N != A->cmap->N && P->rmap->N == P->cmap->N && P->rmap->N == A->cmap->N)) { 1835 Mat B; 1836 PetscCall(PCHPDDMSetAuxiliaryMatNormal_Private(pc, A, P, &B, pcpre)); 1837 if (data->correction == PC_HPDDM_COARSE_CORRECTION_DEFLATED) data->correction = PC_HPDDM_COARSE_CORRECTION_BALANCED; 1838 PetscCall(MatDestroy(&B)); 1839 } else { 1840 PetscCall(PetscObjectTypeCompare((PetscObject)P, MATSCHURCOMPLEMENT, &flg)); 1841 if (flg) { 1842 Mat A00, P00, A01, A10, A11, B, N; 1843 PCHPDDMSchurPreType type = PC_HPDDM_SCHUR_PRE_LEAST_SQUARES; 1844 1845 PetscCall(MatSchurComplementGetSubMatrices(P, &A00, &P00, &A01, &A10, &A11)); 1846 if (PetscDefined(USE_DEBUG)) PetscCall(PCHPDDMCheckSymmetry_Private(pc, A01, A10)); 1847 PetscCall(PetscOptionsGetEnum(nullptr, pcpre, "-pc_hpddm_schur_precondition", PCHPDDMSchurPreTypes, (PetscEnum *)&type, &flg)); 1848 if (type == PC_HPDDM_SCHUR_PRE_LEAST_SQUARES) { 1849 Vec diagonal = nullptr; 1850 const PetscScalar *array; 1851 MatSchurComplementAinvType type; 1852 1853 if (A11) { 1854 PetscCall(MatCreateVecs(A11, &diagonal, nullptr)); 1855 PetscCall(MatGetDiagonal(A11, diagonal)); 1856 } 1857 PetscCall(MatCreateVecs(P00, &v, nullptr)); 1858 PetscCall(MatSchurComplementGetAinvType(P, &type)); 1859 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]); 1860 if (type == MAT_SCHUR_COMPLEMENT_AINV_LUMP) { 1861 PetscCall(MatGetRowSum(P00, v)); 1862 if (A00 == P00) PetscCall(PetscObjectReference((PetscObject)A00)); 1863 PetscCall(MatDestroy(&P00)); 1864 PetscCall(VecGetArrayRead(v, &array)); 1865 PetscCall(MatCreateAIJ(PetscObjectComm((PetscObject)A00), A00->rmap->n, A00->cmap->n, A00->rmap->N, A00->cmap->N, 1, nullptr, 0, nullptr, &P00)); 1866 PetscCall(MatSetOption(P00, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE)); 1867 for (n = A00->rmap->rstart; n < A00->rmap->rend; ++n) PetscCall(MatSetValue(P00, n, n, array[n - A00->rmap->rstart], INSERT_VALUES)); 1868 PetscCall(MatAssemblyBegin(P00, MAT_FINAL_ASSEMBLY)); 1869 PetscCall(MatAssemblyEnd(P00, MAT_FINAL_ASSEMBLY)); 1870 PetscCall(VecRestoreArrayRead(v, &array)); 1871 PetscCall(MatSchurComplementUpdateSubMatrices(P, A00, P00, A01, A10, A11)); /* replace P00 by diag(sum of each row of P00) */ 1872 PetscCall(MatDestroy(&P00)); 1873 } else PetscCall(MatGetDiagonal(P00, v)); 1874 PetscCall(VecReciprocal(v)); /* inv(diag(P00)) */ 1875 PetscCall(VecSqrtAbs(v)); /* sqrt(inv(diag(P00))) */ 1876 PetscCall(MatDuplicate(A01, MAT_COPY_VALUES, &B)); 1877 PetscCall(MatDiagonalScale(B, v, nullptr)); 1878 PetscCall(VecDestroy(&v)); 1879 PetscCall(MatCreateNormalHermitian(B, &N)); 1880 PetscCall(PCHPDDMSetAuxiliaryMatNormal_Private(pc, B, N, &P, pcpre, &diagonal)); 1881 PetscCall(PetscObjectTypeCompare((PetscObject)data->aux, MATSEQAIJ, &flg)); 1882 if (!flg) { 1883 PetscCall(MatDestroy(&P)); 1884 P = N; 1885 PetscCall(PetscObjectReference((PetscObject)P)); 1886 } else PetscCall(MatScale(P, -1.0)); 1887 if (diagonal) { 1888 PetscCall(MatDiagonalSet(P, diagonal, ADD_VALUES)); 1889 PetscCall(PCSetOperators(pc, P, P)); /* replace P by diag(P11) - A01^T inv(diag(P00)) A01 */ 1890 PetscCall(VecDestroy(&diagonal)); 1891 } else { 1892 PetscCall(MatScale(N, -1.0)); 1893 PetscCall(PCSetOperators(pc, N, P)); /* replace P by - A01^T inv(diag(P00)) A01 */ 1894 } 1895 PetscCall(MatDestroy(&N)); 1896 PetscCall(MatDestroy(&P)); 1897 PetscCall(MatDestroy(&B)); 1898 } else 1899 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 : ""); 1900 for (std::vector<Vec>::iterator it = initial.begin(); it != initial.end(); ++it) PetscCall(VecDestroy(&*it)); 1901 PetscFunctionReturn(PETSC_SUCCESS); 1902 } else { 1903 PetscCall(PetscOptionsGetString(nullptr, pcpre, "-pc_hpddm_levels_1_st_pc_type", type, sizeof(type), nullptr)); 1904 PetscCall(PetscStrcmp(type, PCMAT, &algebraic)); 1905 PetscCheck(!algebraic || !block, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-pc_hpddm_levels_1_st_pc_type mat and -pc_hpddm_block_splitting"); 1906 if (overlap != -1) { 1907 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"); 1908 PetscCheck(overlap >= 1, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_WRONG, "-pc_hpddm_harmonic_overlap %" PetscInt_FMT " < 1", overlap); 1909 } 1910 if (block || overlap != -1) algebraic = PETSC_TRUE; 1911 if (algebraic) { 1912 PetscCall(ISCreateStride(PETSC_COMM_SELF, P->rmap->n, P->rmap->rstart, 1, &data->is)); 1913 PetscCall(MatIncreaseOverlap(P, 1, &data->is, 1)); 1914 PetscCall(ISSort(data->is)); 1915 } else 1916 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 : "")); 1917 } 1918 } 1919 } 1920 } 1921 #if PetscDefined(USE_DEBUG) 1922 if (data->is) PetscCall(ISDuplicate(data->is, &dis)); 1923 if (data->aux) PetscCall(MatDuplicate(data->aux, MAT_COPY_VALUES, &daux)); 1924 #endif 1925 if (data->is || (ismatis && data->N > 1)) { 1926 if (ismatis) { 1927 std::initializer_list<std::string> list = {MATSEQBAIJ, MATSEQSBAIJ}; 1928 PetscCall(MatISGetLocalMat(P, &N)); 1929 std::initializer_list<std::string>::const_iterator it = std::find(list.begin(), list.end(), ((PetscObject)N)->type_name); 1930 PetscCall(MatISRestoreLocalMat(P, &N)); 1931 switch (std::distance(list.begin(), it)) { 1932 case 0: 1933 PetscCall(MatConvert(P, MATMPIBAIJ, MAT_INITIAL_MATRIX, &C)); 1934 break; 1935 case 1: 1936 /* MatCreateSubMatrices() does not work with MATSBAIJ and unsorted ISes, so convert to MPIBAIJ */ 1937 PetscCall(MatConvert(P, MATMPIBAIJ, MAT_INITIAL_MATRIX, &C)); 1938 PetscCall(MatSetOption(C, MAT_SYMMETRIC, PETSC_TRUE)); 1939 break; 1940 default: 1941 PetscCall(MatConvert(P, MATMPIAIJ, MAT_INITIAL_MATRIX, &C)); 1942 } 1943 PetscCall(MatISGetLocalToGlobalMapping(P, &l2g, nullptr)); 1944 PetscCall(PetscObjectReference((PetscObject)P)); 1945 PetscCall(KSPSetOperators(data->levels[0]->ksp, A, C)); 1946 std::swap(C, P); 1947 PetscCall(ISLocalToGlobalMappingGetSize(l2g, &n)); 1948 PetscCall(ISCreateStride(PETSC_COMM_SELF, n, 0, 1, &loc)); 1949 PetscCall(ISLocalToGlobalMappingApplyIS(l2g, loc, &is[0])); 1950 PetscCall(ISDestroy(&loc)); 1951 /* the auxiliary Mat is _not_ the local Neumann matrix */ 1952 /* it is the local Neumann matrix augmented (with zeros) through MatIncreaseOverlap() */ 1953 data->Neumann = PETSC_BOOL3_FALSE; 1954 structure = SAME_NONZERO_PATTERN; 1955 } else { 1956 is[0] = data->is; 1957 if (algebraic || ctx) subdomains = PETSC_TRUE; 1958 PetscCall(PetscOptionsGetBool(nullptr, pcpre, "-pc_hpddm_define_subdomains", &subdomains, nullptr)); 1959 if (ctx) PetscCheck(subdomains, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-%spc_hpddm_schur_precondition geneo and -%spc_hpddm_define_subdomains false", pcpre, pcpre); 1960 if (PetscBool3ToBool(data->Neumann)) { 1961 PetscCheck(!block, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-pc_hpddm_block_splitting and -pc_hpddm_has_neumann"); 1962 PetscCheck(overlap == -1, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-pc_hpddm_harmonic_overlap %" PetscInt_FMT " and -pc_hpddm_has_neumann", overlap); 1963 PetscCheck(!algebraic, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-pc_hpddm_levels_1_st_pc_type mat and -pc_hpddm_has_neumann"); 1964 } 1965 if (PetscBool3ToBool(data->Neumann) || block) structure = SAME_NONZERO_PATTERN; 1966 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)data->is), P->rmap->n, P->rmap->rstart, 1, &loc)); 1967 } 1968 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_levels_1_", pcpre ? pcpre : "")); 1969 PetscCall(PetscOptionsGetEnum(nullptr, prefix, "-st_matstructure", MatStructures, (PetscEnum *)&structure, &flg)); /* if not user-provided, force its value when possible */ 1970 if (!flg && structure == SAME_NONZERO_PATTERN) { /* cannot call STSetMatStructure() yet, insert the appropriate option in the database, parsed by STSetFromOptions() */ 1971 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-%spc_hpddm_levels_1_st_matstructure", pcpre ? pcpre : "")); 1972 PetscCall(PetscOptionsSetValue(nullptr, prefix, MatStructures[structure])); 1973 } 1974 flg = PETSC_FALSE; 1975 if (data->share) { 1976 data->share = PETSC_FALSE; /* will be reset to PETSC_TRUE if none of the conditions below are true */ 1977 if (!subdomains) PetscCall(PetscInfo(pc, "Cannot share subdomain KSP between SLEPc and PETSc since -%spc_hpddm_define_subdomains is not true\n", pcpre ? pcpre : "")); 1978 else if (data->deflation) PetscCall(PetscInfo(pc, "Nothing to share since PCHPDDMSetDeflationMat() has been called\n")); 1979 else if (ismatis) PetscCall(PetscInfo(pc, "Cannot share subdomain KSP between SLEPc and PETSc with a Pmat of type MATIS\n")); 1980 else if (!algebraic && structure != SAME_NONZERO_PATTERN) 1981 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])); 1982 else data->share = PETSC_TRUE; 1983 } 1984 if (!ismatis) { 1985 if (data->share || (!PetscBool3ToBool(data->Neumann) && subdomains)) PetscCall(ISDuplicate(is[0], &unsorted)); 1986 else unsorted = is[0]; 1987 } 1988 if (data->N > 1 && (data->aux || ismatis || algebraic)) { 1989 PetscCheck(loadedSym, PETSC_COMM_SELF, PETSC_ERR_PLIB, "HPDDM library not loaded, cannot use more than one level"); 1990 PetscCall(MatSetOption(P, MAT_SUBMAT_SINGLEIS, PETSC_TRUE)); 1991 if (ismatis) { 1992 /* needed by HPDDM (currently) so that the partition of unity is 0 on subdomain interfaces */ 1993 PetscCall(MatIncreaseOverlap(P, 1, is, 1)); 1994 PetscCall(ISDestroy(&data->is)); 1995 data->is = is[0]; 1996 } else { 1997 if (PetscDefined(USE_DEBUG)) { 1998 PetscBool equal; 1999 IS intersect; 2000 2001 PetscCall(ISIntersect(data->is, loc, &intersect)); 2002 PetscCall(ISEqualUnsorted(loc, intersect, &equal)); 2003 PetscCall(ISDestroy(&intersect)); 2004 PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "IS of the auxiliary Mat does not include all local rows of A"); 2005 } 2006 if (overlap == -1) PetscCall(PetscObjectComposeFunction((PetscObject)pc->pmat, "PCHPDDMAlgebraicAuxiliaryMat_Private_C", PCHPDDMAlgebraicAuxiliaryMat_Private)); 2007 if (!PetscBool3ToBool(data->Neumann) && (!algebraic || overlap != -1)) { 2008 PetscCall(PetscObjectTypeCompare((PetscObject)P, MATMPISBAIJ, &flg)); 2009 if (flg) { 2010 /* maybe better to ISSort(is[0]), MatCreateSubMatrices(), and then MatPermute() */ 2011 /* but there is no MatPermute_SeqSBAIJ(), so as before, just use MATMPIBAIJ */ 2012 PetscCall(MatConvert(P, MATMPIBAIJ, MAT_INITIAL_MATRIX, &uaux)); 2013 flg = PETSC_FALSE; 2014 } 2015 } 2016 } 2017 if (algebraic && overlap == -1) { 2018 PetscUseMethod(pc->pmat, "PCHPDDMAlgebraicAuxiliaryMat_Private_C", (Mat, IS *, Mat *[], PetscBool), (P, is, &sub, block)); 2019 if (block) { 2020 PetscCall(PetscObjectQuery((PetscObject)sub[0], "_PCHPDDM_Neumann_Mat", (PetscObject *)&data->aux)); 2021 PetscCall(PetscObjectCompose((PetscObject)sub[0], "_PCHPDDM_Neumann_Mat", nullptr)); 2022 } 2023 } else if (!uaux || overlap != -1) { 2024 if (!ctx) { 2025 if (PetscBool3ToBool(data->Neumann)) sub = &data->aux; 2026 else { 2027 PetscBool flg; 2028 if (overlap != -1) { 2029 Harmonic h; 2030 Mat A0, *a; /* with an SVD: [ A_00 A_01 ] */ 2031 IS ov[2], rows, cols, stride; /* [ A_10 A_11 A_12 ] */ 2032 const PetscInt *i[2], bs = PetscAbs(P->cmap->bs); /* with a GEVP: [ A_00 A_01 ] */ 2033 PetscInt n[2]; /* [ A_10 A_11 A_12 ] */ 2034 std::vector<PetscInt> v[2]; /* [ A_21 A_22 ] */ 2035 2036 PetscCall(ISDuplicate(data->is, ov)); 2037 if (overlap > 1) PetscCall(MatIncreaseOverlap(P, 1, ov, overlap - 1)); 2038 PetscCall(ISDuplicate(ov[0], ov + 1)); 2039 PetscCall(MatIncreaseOverlap(P, 1, ov + 1, 1)); 2040 PetscCall(PetscNew(&h)); 2041 h->ksp = nullptr; 2042 PetscCall(PetscCalloc1(2, &h->A)); 2043 PetscCall(PetscOptionsHasName(nullptr, prefix, "-svd_nsv", &flg)); 2044 if (!flg) PetscCall(PetscOptionsHasName(nullptr, prefix, "-svd_relative_threshold", &flg)); 2045 PetscCall(ISSort(ov[0])); 2046 if (!flg) PetscCall(ISSort(ov[1])); 2047 PetscCall(PetscMalloc1(!flg ? 5 : 3, &h->is)); 2048 PetscCall(MatCreateSubMatrices(uaux ? uaux : P, 1, ov + !flg, ov + 1, MAT_INITIAL_MATRIX, &a)); /* submatrix from above, either square (!flg) or rectangular (flg) */ 2049 for (PetscInt j = 0; j < 2; ++j) { 2050 PetscCall(ISGetIndices(ov[j], i + j)); 2051 PetscCall(ISGetLocalSize(ov[j], n + j)); 2052 } 2053 v[1].reserve((n[1] - n[0]) / bs); 2054 for (PetscInt j = 0; j < n[1]; j += bs) { /* indices of the (2,2) block */ 2055 PetscInt location; 2056 PetscCall(ISLocate(ov[0], i[1][j], &location)); 2057 if (location < 0) v[1].emplace_back(j / bs); 2058 } 2059 if (!flg) { 2060 h->A[1] = a[0]; 2061 PetscCall(PetscObjectReference((PetscObject)h->A[1])); 2062 v[0].reserve((n[0] - P->rmap->n) / bs); 2063 for (PetscInt j = 0; j < n[1]; j += bs) { /* row indices of the (1,2) block */ 2064 PetscInt location; 2065 PetscCall(ISLocate(loc, i[1][j], &location)); 2066 if (location < 0) { 2067 PetscCall(ISLocate(ov[0], i[1][j], &location)); 2068 if (location >= 0) v[0].emplace_back(j / bs); 2069 } 2070 } 2071 PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, v[0].size(), v[0].data(), PETSC_USE_POINTER, &rows)); 2072 PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, v[1].size(), v[1].data(), PETSC_COPY_VALUES, h->is + 4)); 2073 PetscCall(MatCreateSubMatrix(a[0], rows, h->is[4], MAT_INITIAL_MATRIX, h->A)); /* A_12 submatrix from above */ 2074 PetscCall(ISDestroy(&rows)); 2075 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 */ 2076 PetscCall(ISEmbed(ov[0], ov[1], PETSC_TRUE, &rows)); 2077 PetscCall(MatCreateSubMatrix(a[0], rows, cols = rows, MAT_INITIAL_MATRIX, &A0)); /* [ A_00 A_01 ; A_10 A_11 ] submatrix from above */ 2078 PetscCall(ISDestroy(&rows)); 2079 v[0].clear(); 2080 PetscCall(ISEmbed(loc, ov[1], PETSC_TRUE, h->is + 3)); 2081 PetscCall(ISEmbed(data->is, ov[1], PETSC_TRUE, h->is + 2)); 2082 } 2083 v[0].reserve((n[0] - P->rmap->n) / bs); 2084 for (PetscInt j = 0; j < n[0]; j += bs) { 2085 PetscInt location; 2086 PetscCall(ISLocate(loc, i[0][j], &location)); 2087 if (location < 0) v[0].emplace_back(j / bs); 2088 } 2089 PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, v[0].size(), v[0].data(), PETSC_USE_POINTER, &rows)); 2090 for (PetscInt j = 0; j < 2; ++j) PetscCall(ISRestoreIndices(ov[j], i + j)); 2091 if (flg) { 2092 IS is; 2093 PetscCall(ISCreateStride(PETSC_COMM_SELF, a[0]->rmap->n, 0, 1, &is)); 2094 PetscCall(ISEmbed(ov[0], ov[1], PETSC_TRUE, &cols)); 2095 PetscCall(MatCreateSubMatrix(a[0], is, cols, MAT_INITIAL_MATRIX, &A0)); /* [ A_00 A_01 ; A_10 A_11 ] submatrix from above */ 2096 PetscCall(ISDestroy(&cols)); 2097 PetscCall(ISDestroy(&is)); 2098 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 */ 2099 PetscCall(ISEmbed(loc, data->is, PETSC_TRUE, h->is + 2)); 2100 PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, v[1].size(), v[1].data(), PETSC_USE_POINTER, &cols)); 2101 PetscCall(MatCreateSubMatrix(a[0], rows, cols, MAT_INITIAL_MATRIX, h->A)); /* A_12 submatrix from above */ 2102 PetscCall(ISDestroy(&cols)); 2103 } 2104 PetscCall(ISCreateStride(PETSC_COMM_SELF, A0->rmap->n, 0, 1, &stride)); 2105 PetscCall(ISEmbed(rows, stride, PETSC_TRUE, h->is)); 2106 PetscCall(ISDestroy(&stride)); 2107 PetscCall(ISDestroy(&rows)); 2108 PetscCall(ISEmbed(loc, ov[0], PETSC_TRUE, h->is + 1)); 2109 if (subdomains) { 2110 if (!data->levels[0]->pc) { 2111 PetscCall(PCCreate(PetscObjectComm((PetscObject)pc), &data->levels[0]->pc)); 2112 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_levels_1_", pcpre ? pcpre : "")); 2113 PetscCall(PCSetOptionsPrefix(data->levels[0]->pc, prefix)); 2114 PetscCall(PCSetOperators(data->levels[0]->pc, A, P)); 2115 } 2116 PetscCall(PCSetType(data->levels[0]->pc, PCASM)); 2117 if (!data->levels[0]->pc->setupcalled) PetscCall(PCASMSetLocalSubdomains(data->levels[0]->pc, 1, ov + !flg, &loc)); 2118 PetscCall(PCHPDDMCommunicationAvoidingPCASM_Private(data->levels[0]->pc, flg ? A0 : a[0], PETSC_TRUE)); 2119 if (!flg) ++overlap; 2120 if (data->share) { 2121 PetscInt n = -1; 2122 PetscTryMethod(data->levels[0]->pc, "PCASMGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (data->levels[0]->pc, &n, nullptr, &ksp)); 2123 PetscCheck(n == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of subdomain solver %" PetscInt_FMT " != 1", n); 2124 if (flg) { 2125 h->ksp = ksp[0]; 2126 PetscCall(PetscObjectReference((PetscObject)h->ksp)); 2127 } 2128 } 2129 } 2130 if (!h->ksp) { 2131 PetscBool share = data->share; 2132 PetscCall(KSPCreate(PETSC_COMM_SELF, &h->ksp)); 2133 PetscCall(KSPSetType(h->ksp, KSPPREONLY)); 2134 PetscCall(KSPSetOperators(h->ksp, A0, A0)); 2135 do { 2136 if (!data->share) { 2137 share = PETSC_FALSE; 2138 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_levels_1_%s", pcpre ? pcpre : "", flg ? "svd_" : "eps_")); 2139 PetscCall(KSPSetOptionsPrefix(h->ksp, prefix)); 2140 PetscCall(KSPSetFromOptions(h->ksp)); 2141 } else { 2142 MatSolverType type; 2143 PetscCall(KSPGetPC(ksp[0], &pc)); 2144 PetscCall(PetscObjectTypeCompareAny((PetscObject)pc, &data->share, PCLU, PCCHOLESKY, "")); 2145 if (data->share) { 2146 PetscCall(PCFactorGetMatSolverType(pc, &type)); 2147 if (!type) { 2148 if (PetscDefined(HAVE_MUMPS)) PetscCall(PCFactorSetMatSolverType(pc, MATSOLVERMUMPS)); 2149 else if (PetscDefined(HAVE_MKL_PARDISO)) PetscCall(PCFactorSetMatSolverType(pc, MATSOLVERMKL_PARDISO)); 2150 else data->share = PETSC_FALSE; 2151 if (data->share) PetscCall(PCSetFromOptions(pc)); 2152 } else { 2153 PetscCall(PetscStrcmp(type, MATSOLVERMUMPS, &data->share)); 2154 if (!data->share) PetscCall(PetscStrcmp(type, MATSOLVERMKL_PARDISO, &data->share)); 2155 } 2156 if (data->share) { 2157 std::tuple<KSP, IS, Vec[2]> *p; 2158 PetscCall(PCFactorGetMatrix(pc, &A)); 2159 PetscCall(MatFactorSetSchurIS(A, h->is[4])); 2160 PetscCall(KSPSetUp(ksp[0])); 2161 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_levels_1_eps_shell_", pcpre ? pcpre : "")); 2162 PetscCall(KSPSetOptionsPrefix(h->ksp, prefix)); 2163 PetscCall(KSPSetFromOptions(h->ksp)); 2164 PetscCall(KSPGetPC(h->ksp, &pc)); 2165 PetscCall(PCSetType(pc, PCSHELL)); 2166 PetscCall(PetscNew(&p)); 2167 std::get<0>(*p) = ksp[0]; 2168 PetscCall(ISEmbed(ov[0], ov[1], PETSC_TRUE, &std::get<1>(*p))); 2169 PetscCall(MatCreateVecs(A, std::get<2>(*p), std::get<2>(*p) + 1)); 2170 PetscCall(PCShellSetContext(pc, p)); 2171 PetscCall(PCShellSetApply(pc, PCApply_Schur)); 2172 PetscCall(PCShellSetApplyTranspose(pc, PCApply_Schur<Vec, true>)); 2173 PetscCall(PCShellSetMatApply(pc, PCApply_Schur<Mat>)); 2174 PetscCall(PCShellSetDestroy(pc, PCDestroy_Schur)); 2175 } 2176 } 2177 if (!data->share) PetscCall(PetscInfo(pc, "Cannot share subdomain KSP between SLEPc and PETSc since neither MUMPS nor MKL PARDISO is used\n")); 2178 } 2179 } 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 */ 2180 } 2181 PetscCall(ISDestroy(ov)); 2182 PetscCall(ISDestroy(ov + 1)); 2183 if (overlap == 1 && subdomains && flg) { 2184 *subA = A0; 2185 sub = subA; 2186 if (uaux) PetscCall(MatDestroy(&uaux)); 2187 } else PetscCall(MatDestroy(&A0)); 2188 PetscCall(MatCreateShell(PETSC_COMM_SELF, P->rmap->n, n[1] - n[0], P->rmap->n, n[1] - n[0], h, &data->aux)); 2189 PetscCall(MatCreateVecs(h->ksp->pc->pmat, &h->v, nullptr)); 2190 PetscCall(MatShellSetOperation(data->aux, MATOP_MULT, (void (*)(void))MatMult_Harmonic)); 2191 PetscCall(MatShellSetOperation(data->aux, MATOP_MULT_TRANSPOSE, (void (*)(void))MatMultTranspose_Harmonic)); 2192 PetscCall(MatShellSetMatProductOperation(data->aux, MATPRODUCT_AB, nullptr, MatProduct_AB_Harmonic, nullptr, MATDENSE, MATDENSE)); 2193 PetscCall(MatShellSetOperation(data->aux, MATOP_DESTROY, (void (*)(void))MatDestroy_Harmonic)); 2194 PetscCall(MatDestroySubMatrices(1, &a)); 2195 } 2196 if (overlap != 1 || !subdomains) { 2197 PetscCall(MatCreateSubMatrices(uaux ? uaux : P, 1, is, is, MAT_INITIAL_MATRIX, &sub)); 2198 if (ismatis) { 2199 PetscCall(MatISGetLocalMat(C, &N)); 2200 PetscCall(PetscObjectTypeCompare((PetscObject)N, MATSEQSBAIJ, &flg)); 2201 if (flg) PetscCall(MatConvert(sub[0], MATSEQSBAIJ, MAT_INPLACE_MATRIX, sub)); 2202 PetscCall(MatISRestoreLocalMat(C, &N)); 2203 } 2204 } 2205 if (uaux) { 2206 PetscCall(MatDestroy(&uaux)); 2207 PetscCall(MatConvert(sub[0], MATSEQSBAIJ, MAT_INPLACE_MATRIX, sub)); 2208 } 2209 } 2210 } 2211 } else { 2212 PetscCall(MatCreateSubMatrices(uaux, 1, is, is, MAT_INITIAL_MATRIX, &sub)); 2213 PetscCall(MatDestroy(&uaux)); 2214 PetscCall(MatConvert(sub[0], MATSEQSBAIJ, MAT_INPLACE_MATRIX, sub)); 2215 } 2216 /* Vec holding the partition of unity */ 2217 if (!data->levels[0]->D) { 2218 PetscCall(ISGetLocalSize(data->is, &n)); 2219 PetscCall(VecCreateMPI(PETSC_COMM_SELF, n, PETSC_DETERMINE, &data->levels[0]->D)); 2220 } 2221 if (data->share && overlap == -1) { 2222 Mat D; 2223 IS perm = nullptr; 2224 PetscInt size = -1; 2225 if (!data->levels[0]->pc) { 2226 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_levels_1_", pcpre ? pcpre : "")); 2227 PetscCall(PCCreate(PetscObjectComm((PetscObject)pc), &data->levels[0]->pc)); 2228 PetscCall(PCSetOptionsPrefix(data->levels[0]->pc, prefix)); 2229 PetscCall(PCSetOperators(data->levels[0]->pc, A, P)); 2230 } 2231 PetscCall(PCSetType(data->levels[0]->pc, PCASM)); 2232 if (!ctx) { 2233 if (!data->levels[0]->pc->setupcalled) { 2234 IS sorted; /* PCASM will sort the input IS, duplicate it to return an unmodified (PCHPDDM) input IS */ 2235 PetscCall(ISDuplicate(is[0], &sorted)); 2236 PetscCall(PCASMSetLocalSubdomains(data->levels[0]->pc, 1, &sorted, &loc)); 2237 PetscCall(PetscObjectDereference((PetscObject)sorted)); 2238 } 2239 PetscCall(PCSetFromOptions(data->levels[0]->pc)); 2240 if (block) { 2241 PetscCall(PCHPDDMPermute_Private(unsorted, data->is, &uis, sub[0], &C, &perm)); 2242 PetscCall(PCHPDDMCommunicationAvoidingPCASM_Private(data->levels[0]->pc, C, algebraic)); 2243 } else PetscCall(PCSetUp(data->levels[0]->pc)); 2244 PetscTryMethod(data->levels[0]->pc, "PCASMGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (data->levels[0]->pc, &size, nullptr, &ksp)); 2245 if (size != 1) { 2246 data->share = PETSC_FALSE; 2247 PetscCheck(size == -1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of subdomain solver %" PetscInt_FMT " != 1", size); 2248 PetscCall(PetscInfo(pc, "Cannot share subdomain KSP between SLEPc and PETSc since PCASMGetSubKSP() not found in fine-level PC\n")); 2249 PetscCall(ISDestroy(&unsorted)); 2250 unsorted = is[0]; 2251 } else { 2252 if (!block && !ctx) PetscCall(PCHPDDMPermute_Private(unsorted, data->is, &uis, PetscBool3ToBool(data->Neumann) ? sub[0] : data->aux, &C, &perm)); 2253 if (!PetscBool3ToBool(data->Neumann) && !block) { 2254 PetscCall(MatPermute(sub[0], perm, perm, &D)); /* permute since PCASM will call ISSort() */ 2255 PetscCall(MatHeaderReplace(sub[0], &D)); 2256 } 2257 if (data->B) { /* see PCHPDDMSetRHSMat() */ 2258 PetscCall(MatPermute(data->B, perm, perm, &D)); 2259 PetscCall(MatHeaderReplace(data->B, &D)); 2260 } 2261 PetscCall(ISDestroy(&perm)); 2262 const char *matpre; 2263 PetscBool cmp[4]; 2264 PetscCall(KSPGetOperators(ksp[0], subA, subA + 1)); 2265 PetscCall(PetscObjectReference((PetscObject)subA[0])); 2266 PetscCall(MatDuplicate(subA[1], MAT_SHARE_NONZERO_PATTERN, &D)); 2267 PetscCall(MatGetOptionsPrefix(subA[1], &matpre)); 2268 PetscCall(MatSetOptionsPrefix(D, matpre)); 2269 PetscCall(PetscObjectTypeCompare((PetscObject)D, MATNORMAL, cmp)); 2270 PetscCall(PetscObjectTypeCompare((PetscObject)C, MATNORMAL, cmp + 1)); 2271 if (!cmp[0]) PetscCall(PetscObjectTypeCompare((PetscObject)D, MATNORMALHERMITIAN, cmp + 2)); 2272 else cmp[2] = PETSC_FALSE; 2273 if (!cmp[1]) PetscCall(PetscObjectTypeCompare((PetscObject)C, MATNORMALHERMITIAN, cmp + 3)); 2274 else cmp[3] = PETSC_FALSE; 2275 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); 2276 if (!cmp[0] && !cmp[2]) { 2277 if (!block) PetscCall(MatAXPY(D, 1.0, C, SUBSET_NONZERO_PATTERN)); 2278 else { 2279 PetscCall(MatMissingDiagonal(D, cmp, nullptr)); 2280 if (cmp[0]) structure = DIFFERENT_NONZERO_PATTERN; /* data->aux has no missing diagonal entry */ 2281 PetscCall(MatAXPY(D, 1.0, data->aux, structure)); 2282 } 2283 } else { 2284 Mat mat[2]; 2285 if (cmp[0]) { 2286 PetscCall(MatNormalGetMat(D, mat)); 2287 PetscCall(MatNormalGetMat(C, mat + 1)); 2288 } else { 2289 PetscCall(MatNormalHermitianGetMat(D, mat)); 2290 PetscCall(MatNormalHermitianGetMat(C, mat + 1)); 2291 } 2292 PetscCall(MatAXPY(mat[0], 1.0, mat[1], SUBSET_NONZERO_PATTERN)); 2293 } 2294 PetscCall(MatPropagateSymmetryOptions(C, D)); 2295 PetscCall(MatDestroy(&C)); 2296 C = D; 2297 /* swap pointers so that variables stay consistent throughout PCSetUp() */ 2298 std::swap(C, data->aux); 2299 std::swap(uis, data->is); 2300 swap = PETSC_TRUE; 2301 } 2302 } 2303 } 2304 if (ctx) { 2305 PC_HPDDM *data_00 = (PC_HPDDM *)std::get<0>(*ctx)[0]->data; 2306 PC s; 2307 Mat A00, P00, A01 = nullptr, A10, A11, N, b[4]; 2308 IS sorted, is[2], *is_00; 2309 MatSolverType type; 2310 std::pair<PC, Vec[2]> *p; 2311 2312 n = -1; 2313 PetscTryMethod(data_00->levels[0]->pc, "PCASMGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (data_00->levels[0]->pc, &n, nullptr, &ksp)); 2314 PetscCheck(n == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of subdomain solver %" PetscInt_FMT " != 1", n); 2315 PetscCall(KSPGetOperators(ksp[0], subA, subA + 1)); 2316 PetscCall(ISGetLocalSize(data_00->is, &n)); 2317 if (n != subA[0]->rmap->n || n != subA[0]->cmap->n) { 2318 PetscCall(PCASMGetLocalSubdomains(data_00->levels[0]->pc, &n, &is_00, nullptr)); 2319 PetscCall(ISGetLocalSize(*is_00, &n)); 2320 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); 2321 } else is_00 = &data_00->is; 2322 PetscCall(PCHPDDMPermute_Private(unsorted, data->is, &uis, data->aux, &C, nullptr)); /* permute since PCASM works with a sorted IS */ 2323 std::swap(C, data->aux); 2324 std::swap(uis, data->is); 2325 swap = PETSC_TRUE; 2326 PetscCall(PCASMSetType(data->levels[0]->pc, PC_ASM_NONE)); /* "Neumann--Neumann" preconditioning with overlap and a Boolean partition of unity */ 2327 PetscCall(PCASMSetLocalSubdomains(data->levels[0]->pc, 1, &data->is, &loc)); 2328 PetscCall(PCSetFromOptions(data->levels[0]->pc)); /* action of eq. (15) of https://hal.science/hal-02343808v6/document (with a sign flip) */ 2329 PetscCall(MatSchurComplementGetSubMatrices(P, &A00, &P00, std::get<1>(*ctx), &A10, &A11)); 2330 std::get<1>(*ctx)[1] = A10; 2331 PetscCall(PetscObjectTypeCompare((PetscObject)A10, MATTRANSPOSEVIRTUAL, &flg)); 2332 if (flg) PetscCall(MatTransposeGetMat(A10, &A01)); 2333 else { 2334 PetscBool flg; 2335 2336 PetscCall(PetscObjectTypeCompare((PetscObject)A10, MATHERMITIANTRANSPOSEVIRTUAL, &flg)); 2337 if (flg) PetscCall(MatHermitianTransposeGetMat(A10, &A01)); 2338 } 2339 PetscCall(ISDuplicate(*is_00, &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 */ 2340 PetscCall(ISSort(sorted)); /* this is to avoid changing users inputs, but it requires a new call to ISSort() here */ 2341 if (!A01) { 2342 PetscCall(MatSetOption(A10, MAT_SUBMAT_SINGLEIS, PETSC_TRUE)); 2343 PetscCall(MatCreateSubMatrices(A10, 1, &data->is, &sorted, MAT_INITIAL_MATRIX, &sub)); 2344 b[2] = sub[0]; 2345 PetscCall(PetscObjectReference((PetscObject)sub[0])); 2346 PetscCall(MatDestroySubMatrices(1, &sub)); 2347 PetscCall(PetscObjectTypeCompare((PetscObject)std::get<1>(*ctx)[0], MATTRANSPOSEVIRTUAL, &flg)); 2348 A10 = nullptr; 2349 if (flg) PetscCall(MatTransposeGetMat(std::get<1>(*ctx)[0], &A10)); 2350 else { 2351 PetscBool flg; 2352 2353 PetscCall(PetscObjectTypeCompare((PetscObject)(PetscObject)std::get<1>(*ctx)[0], MATHERMITIANTRANSPOSEVIRTUAL, &flg)); 2354 if (flg) PetscCall(MatHermitianTransposeGetMat(std::get<1>(*ctx)[0], &A10)); 2355 } 2356 if (!A10) { 2357 PetscCall(MatCreateSubMatrices(std::get<1>(*ctx)[0], 1, &sorted, &data->is, MAT_INITIAL_MATRIX, &sub)); 2358 b[1] = sub[0]; 2359 PetscCall(PetscObjectReference((PetscObject)sub[0])); 2360 } 2361 } else { 2362 PetscCall(MatSetOption(A01, MAT_SUBMAT_SINGLEIS, PETSC_TRUE)); 2363 PetscCall(MatCreateSubMatrices(A01, 1, &sorted, &data->is, MAT_INITIAL_MATRIX, &sub)); 2364 if (flg) PetscCall(MatTranspose(*sub, MAT_INITIAL_MATRIX, b + 2)); 2365 else PetscCall(MatHermitianTranspose(*sub, MAT_INITIAL_MATRIX, b + 2)); 2366 } 2367 PetscCall(MatDestroySubMatrices(1, &sub)); 2368 PetscCall(ISDestroy(&sorted)); 2369 if (A01 || A10) { 2370 if (flg) PetscCall(MatTranspose(b[2], MAT_INITIAL_MATRIX, b + 1)); 2371 else PetscCall(MatHermitianTranspose(b[2], MAT_INITIAL_MATRIX, b + 1)); 2372 } 2373 PetscCall(MatCreateSchurComplement(subA[0], subA[1], b[1], b[2], data->aux, &S)); 2374 PetscCall(MatSchurComplementSetKSP(S, ksp[0])); 2375 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 */ 2376 PetscTryMethod(data->levels[0]->pc, "PCASMGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (data->levels[0]->pc, &n, nullptr, &ksp)); 2377 PetscCheck(n == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of subdomain solver %" PetscInt_FMT " != 1", n); 2378 PetscCall(KSPGetPC(ksp[0], &inner)); 2379 PetscCall(PCSetType(inner, PCSHELL)); /* compute the action of the inverse of the local Schur complement with a PCSHELL */ 2380 b[0] = subA[0]; 2381 b[3] = data->aux; 2382 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 */ 2383 PetscCall(PetscObjectDereference((PetscObject)b[1])); 2384 PetscCall(PetscObjectDereference((PetscObject)b[2])); 2385 PetscCall(PCCreate(PETSC_COMM_SELF, &s)); 2386 PetscCall(PCSetOptionsPrefix(s, ((PetscObject)inner)->prefix)); 2387 PetscCall(PCSetOptionsPrefix(inner, nullptr)); 2388 PetscCall(KSPSetSkipPCSetFromOptions(ksp[0], PETSC_TRUE)); 2389 PetscCall(PCSetType(s, PCLU)); 2390 if (PetscDefined(HAVE_MUMPS)) { /* only MATSOLVERMUMPS handles MATNEST, so for the others, e.g., MATSOLVERPETSC or MATSOLVERMKL_PARDISO, convert to plain MATAIJ */ 2391 PetscCall(PCFactorSetMatSolverType(s, MATSOLVERMUMPS)); 2392 } 2393 PetscCall(PCSetFromOptions(s)); 2394 PetscCall(PCFactorGetMatSolverType(s, &type)); 2395 PetscCall(PetscStrcmp(type, MATSOLVERMUMPS, &flg)); 2396 if (flg) { 2397 PetscCall(PCSetOperators(s, N, N)); 2398 PetscCall(PCFactorGetMatrix(s, b)); 2399 PetscCall(MatSetOptionsPrefix(*b, ((PetscObject)s)->prefix)); 2400 n = -1; 2401 PetscCall(PetscOptionsGetInt(nullptr, ((PetscObject)s)->prefix, "-mat_mumps_icntl_26", &n, nullptr)); 2402 if (n == 1) { /* allocates a square MatDense of size is[1]->map->n, so one */ 2403 PetscCall(MatNestGetISs(N, is, nullptr)); /* needs to be able to deactivate this path when dealing */ 2404 PetscCall(MatFactorSetSchurIS(*b, is[1])); /* with a large constraint space in order to avoid OOM */ 2405 } 2406 } else { 2407 PetscCall(MatConvert(N, MATAIJ, MAT_INITIAL_MATRIX, b)); 2408 PetscCall(PCSetOperators(s, N, *b)); 2409 PetscCall(PetscObjectDereference((PetscObject)*b)); 2410 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 */ 2411 } 2412 PetscCall(PetscNew(&p)); 2413 p->first = s; 2414 PetscCall(MatCreateVecs(*b, p->second, p->second + 1)); 2415 PetscCall(PCShellSetContext(inner, p)); 2416 PetscCall(PCShellSetApply(inner, PCApply_Nest)); 2417 PetscCall(PCShellSetView(inner, PCView_Nest)); 2418 PetscCall(PCShellSetDestroy(inner, PCDestroy_Nest)); 2419 PetscCall(PetscObjectDereference((PetscObject)N)); 2420 } 2421 if (!data->levels[0]->scatter) { 2422 PetscCall(MatCreateVecs(P, &xin, nullptr)); 2423 if (ismatis) PetscCall(MatDestroy(&P)); 2424 PetscCall(VecScatterCreate(xin, data->is, data->levels[0]->D, nullptr, &data->levels[0]->scatter)); 2425 PetscCall(VecDestroy(&xin)); 2426 } 2427 if (data->levels[0]->P) { 2428 /* if the pattern is the same and PCSetUp() has previously succeeded, reuse HPDDM buffers and connectivity */ 2429 PetscCall(HPDDM::Schwarz<PetscScalar>::destroy(data->levels[0], pc->setupcalled < 1 || pc->flag == DIFFERENT_NONZERO_PATTERN ? PETSC_TRUE : PETSC_FALSE)); 2430 } 2431 if (!data->levels[0]->P) data->levels[0]->P = new HPDDM::Schwarz<PetscScalar>(); 2432 if (data->log_separate) PetscCall(PetscLogEventBegin(PC_HPDDM_SetUp[0], data->levels[0]->ksp, nullptr, nullptr, nullptr)); 2433 else PetscCall(PetscLogEventBegin(PC_HPDDM_Strc, data->levels[0]->ksp, nullptr, nullptr, nullptr)); 2434 /* HPDDM internal data structure */ 2435 PetscCall(data->levels[0]->P->structure(loc, data->is, !ctx ? sub[0] : nullptr, ismatis ? C : data->aux, data->levels)); 2436 if (!data->log_separate) PetscCall(PetscLogEventEnd(PC_HPDDM_Strc, data->levels[0]->ksp, nullptr, nullptr, nullptr)); 2437 /* matrix pencil of the generalized eigenvalue problem on the overlap (GenEO) */ 2438 if (!ctx) { 2439 if (data->deflation || overlap != -1) weighted = data->aux; 2440 else if (!data->B) { 2441 PetscBool cmp; 2442 PetscCall(MatDuplicate(sub[0], MAT_COPY_VALUES, &weighted)); 2443 PetscCall(PetscObjectTypeCompareAny((PetscObject)weighted, &cmp, MATNORMAL, MATNORMALHERMITIAN, "")); 2444 if (cmp) flg = PETSC_FALSE; 2445 PetscCall(MatDiagonalScale(weighted, data->levels[0]->D, data->levels[0]->D)); 2446 /* neither MatDuplicate() nor MatDiagonalScale() handles the symmetry options, so propagate the options explicitly */ 2447 /* only useful for -mat_type baij -pc_hpddm_levels_1_st_pc_type cholesky (no problem with MATAIJ or MATSBAIJ) */ 2448 PetscCall(MatPropagateSymmetryOptions(sub[0], weighted)); 2449 } else weighted = data->B; 2450 } else weighted = nullptr; 2451 /* SLEPc is used inside the loaded symbol */ 2452 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)); 2453 if (!ctx && data->share && overlap == -1) { 2454 Mat st[2]; 2455 PetscCall(KSPGetOperators(ksp[0], st, st + 1)); 2456 PetscCall(MatCopy(subA[0], st[0], structure)); 2457 if (subA[1] != subA[0] || st[1] != st[0]) PetscCall(MatCopy(subA[1], st[1], SAME_NONZERO_PATTERN)); 2458 PetscCall(PetscObjectDereference((PetscObject)subA[0])); 2459 } 2460 if (data->log_separate) PetscCall(PetscLogEventEnd(PC_HPDDM_SetUp[0], data->levels[0]->ksp, nullptr, nullptr, nullptr)); 2461 if (ismatis) PetscCall(MatISGetLocalMat(C, &N)); 2462 else N = data->aux; 2463 if (!ctx) P = sub[0]; 2464 else P = S; 2465 /* going through the grid hierarchy */ 2466 for (n = 1; n < data->N; ++n) { 2467 if (data->log_separate) PetscCall(PetscLogEventBegin(PC_HPDDM_SetUp[n], data->levels[n]->ksp, nullptr, nullptr, nullptr)); 2468 /* method composed in the loaded symbol since there, SLEPc is used as well */ 2469 PetscTryMethod(data->levels[0]->ksp, "PCHPDDMSetUp_Private_C", (Mat *, Mat *, PetscInt, PetscInt *const, PC_HPDDM_Level **const), (&P, &N, n, &data->N, data->levels)); 2470 if (data->log_separate) PetscCall(PetscLogEventEnd(PC_HPDDM_SetUp[n], data->levels[n]->ksp, nullptr, nullptr, nullptr)); 2471 } 2472 /* reset to NULL to avoid any faulty use */ 2473 PetscCall(PetscObjectComposeFunction((PetscObject)data->levels[0]->ksp, "PCHPDDMSetUp_Private_C", nullptr)); 2474 if (!ismatis) PetscCall(PetscObjectComposeFunction((PetscObject)pc->pmat, "PCHPDDMAlgebraicAuxiliaryMat_C", nullptr)); 2475 else PetscCall(PetscObjectDereference((PetscObject)C)); /* matching PetscObjectReference() above */ 2476 for (n = 0; n < data->N - 1; ++n) 2477 if (data->levels[n]->P) { 2478 /* HPDDM internal work buffers */ 2479 data->levels[n]->P->setBuffer(); 2480 data->levels[n]->P->super::start(); 2481 } 2482 if (ismatis || !subdomains) PetscCall(PCHPDDMDestroySubMatrices_Private(PetscBool3ToBool(data->Neumann), PetscBool(algebraic && !block && overlap == -1), sub)); 2483 if (ismatis) data->is = nullptr; 2484 for (n = 0; n < data->N - 1 + (reused > 0); ++n) { 2485 if (data->levels[n]->P) { 2486 PC spc; 2487 2488 /* force the PC to be PCSHELL to do the coarse grid corrections */ 2489 PetscCall(KSPSetSkipPCSetFromOptions(data->levels[n]->ksp, PETSC_TRUE)); 2490 PetscCall(KSPGetPC(data->levels[n]->ksp, &spc)); 2491 PetscCall(PCSetType(spc, PCSHELL)); 2492 PetscCall(PCShellSetContext(spc, data->levels[n])); 2493 PetscCall(PCShellSetSetUp(spc, PCSetUp_HPDDMShell)); 2494 PetscCall(PCShellSetApply(spc, PCApply_HPDDMShell)); 2495 PetscCall(PCShellSetMatApply(spc, PCMatApply_HPDDMShell)); 2496 if (ctx && n == 0) { 2497 Mat Amat, Pmat; 2498 PetscInt m, M; 2499 std::tuple<Mat, PetscSF, Vec[2]> *ctx; 2500 2501 PetscCall(KSPGetOperators(data->levels[n]->ksp, nullptr, &Pmat)); 2502 PetscCall(MatGetLocalSize(Pmat, &m, nullptr)); 2503 PetscCall(MatGetSize(Pmat, &M, nullptr)); 2504 PetscCall(PetscNew(&ctx)); 2505 std::get<0>(*ctx) = S; 2506 std::get<1>(*ctx) = data->levels[n]->scatter; 2507 PetscCall(MatCreateShell(PetscObjectComm((PetscObject)data->levels[n]->ksp), m, m, M, M, ctx, &Amat)); 2508 PetscCall(MatShellSetOperation(Amat, MATOP_MULT, (void (*)(void))MatMult_Schur<false>)); 2509 PetscCall(MatShellSetOperation(Amat, MATOP_MULT_TRANSPOSE, (void (*)(void))MatMult_Schur<true>)); 2510 PetscCall(MatShellSetOperation(Amat, MATOP_DESTROY, (void (*)(void))MatDestroy_Schur)); 2511 PetscCall(MatCreateVecs(S, std::get<2>(*ctx), std::get<2>(*ctx) + 1)); 2512 PetscCall(KSPSetOperators(data->levels[n]->ksp, Amat, Pmat)); 2513 PetscCall(PetscObjectDereference((PetscObject)Amat)); 2514 } 2515 PetscCall(PCShellSetDestroy(spc, PCDestroy_HPDDMShell)); 2516 if (!data->levels[n]->pc) PetscCall(PCCreate(PetscObjectComm((PetscObject)data->levels[n]->ksp), &data->levels[n]->pc)); 2517 if (n < reused) { 2518 PetscCall(PCSetReusePreconditioner(spc, PETSC_TRUE)); 2519 PetscCall(PCSetReusePreconditioner(data->levels[n]->pc, PETSC_TRUE)); 2520 } 2521 PetscCall(PCSetUp(spc)); 2522 } 2523 } 2524 if (ctx) PetscCall(MatDestroy(&S)); 2525 if (overlap == -1) PetscCall(PetscObjectComposeFunction((PetscObject)pc->pmat, "PCHPDDMAlgebraicAuxiliaryMat_Private_C", nullptr)); 2526 } else flg = reused ? PETSC_FALSE : PETSC_TRUE; 2527 if (!ismatis && subdomains) { 2528 if (flg) PetscCall(KSPGetPC(data->levels[0]->ksp, &inner)); 2529 else inner = data->levels[0]->pc; 2530 if (inner) { 2531 if (!inner->setupcalled) PetscCall(PCSetType(inner, PCASM)); 2532 PetscCall(PCSetFromOptions(inner)); 2533 PetscCall(PetscStrcmp(((PetscObject)inner)->type_name, PCASM, &flg)); 2534 if (flg) { 2535 if (!inner->setupcalled) { /* evaluates to PETSC_FALSE when -pc_hpddm_block_splitting */ 2536 IS sorted; /* PCASM will sort the input IS, duplicate it to return an unmodified (PCHPDDM) input IS */ 2537 PetscCall(ISDuplicate(is[0], &sorted)); 2538 PetscCall(PCASMSetLocalSubdomains(inner, 1, &sorted, &loc)); 2539 PetscCall(PetscObjectDereference((PetscObject)sorted)); 2540 } 2541 if (!PetscBool3ToBool(data->Neumann) && data->N > 1) { /* subdomain matrices are already created for the eigenproblem, reuse them for the fine-level PC */ 2542 PetscCall(PCHPDDMPermute_Private(*is, nullptr, nullptr, sub[0], &P, nullptr)); 2543 PetscCall(PCHPDDMCommunicationAvoidingPCASM_Private(inner, P, algebraic)); 2544 PetscCall(PetscObjectDereference((PetscObject)P)); 2545 } 2546 } 2547 } 2548 if (data->N > 1) { 2549 if (overlap != 1) PetscCall(PCHPDDMDestroySubMatrices_Private(PetscBool3ToBool(data->Neumann), PetscBool(algebraic && !block && overlap == -1), sub)); 2550 if (overlap == 1) PetscCall(MatDestroy(subA)); 2551 } 2552 } 2553 PetscCall(ISDestroy(&loc)); 2554 } else data->N = 1 + reused; /* enforce this value to 1 + reused if there is no way to build another level */ 2555 if (requested != data->N + reused) { 2556 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, 2557 data->N, pcpre ? pcpre : "")); 2558 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)); 2559 /* cannot use PCDestroy_HPDDMShell() because PCSHELL not set for unassembled levels */ 2560 for (n = data->N - 1; n < requested - 1; ++n) { 2561 if (data->levels[n]->P) { 2562 PetscCall(HPDDM::Schwarz<PetscScalar>::destroy(data->levels[n], PETSC_TRUE)); 2563 PetscCall(VecDestroyVecs(1, &data->levels[n]->v[0])); 2564 PetscCall(VecDestroyVecs(2, &data->levels[n]->v[1])); 2565 PetscCall(MatDestroy(data->levels[n]->V)); 2566 PetscCall(MatDestroy(data->levels[n]->V + 1)); 2567 PetscCall(MatDestroy(data->levels[n]->V + 2)); 2568 PetscCall(VecDestroy(&data->levels[n]->D)); 2569 PetscCall(PetscSFDestroy(&data->levels[n]->scatter)); 2570 } 2571 } 2572 if (reused) { 2573 for (n = reused; n < PETSC_PCHPDDM_MAXLEVELS && data->levels[n]; ++n) { 2574 PetscCall(KSPDestroy(&data->levels[n]->ksp)); 2575 PetscCall(PCDestroy(&data->levels[n]->pc)); 2576 } 2577 } 2578 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, 2579 data->N, reused, data->N, pcpre ? pcpre : "", pcpre ? pcpre : "", data->N); 2580 } 2581 /* these solvers are created after PCSetFromOptions() is called */ 2582 if (pc->setfromoptionscalled) { 2583 for (n = 0; n < data->N; ++n) { 2584 if (data->levels[n]->ksp) PetscCall(KSPSetFromOptions(data->levels[n]->ksp)); 2585 if (data->levels[n]->pc) PetscCall(PCSetFromOptions(data->levels[n]->pc)); 2586 } 2587 pc->setfromoptionscalled = 0; 2588 } 2589 data->N += reused; 2590 if (data->share && swap) { 2591 /* swap back pointers so that variables follow the user-provided numbering */ 2592 std::swap(C, data->aux); 2593 std::swap(uis, data->is); 2594 PetscCall(MatDestroy(&C)); 2595 PetscCall(ISDestroy(&uis)); 2596 } 2597 if (algebraic) PetscCall(MatDestroy(&data->aux)); 2598 if (unsorted && unsorted != is[0]) { 2599 PetscCall(ISCopy(unsorted, data->is)); 2600 PetscCall(ISDestroy(&unsorted)); 2601 } 2602 #if PetscDefined(USE_DEBUG) 2603 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); 2604 if (data->is) { 2605 PetscCall(ISEqualUnsorted(data->is, dis, &flg)); 2606 PetscCall(ISDestroy(&dis)); 2607 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Input IS and output IS are not equal"); 2608 } 2609 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); 2610 if (data->aux) { 2611 PetscCall(MatMultEqual(data->aux, daux, 20, &flg)); 2612 PetscCall(MatDestroy(&daux)); 2613 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Input Mat and output Mat are not equal"); 2614 } 2615 #endif 2616 PetscFunctionReturn(PETSC_SUCCESS); 2617 } 2618 2619 /*@ 2620 PCHPDDMSetCoarseCorrectionType - Sets the coarse correction type. 2621 2622 Collective 2623 2624 Input Parameters: 2625 + pc - preconditioner context 2626 - type - `PC_HPDDM_COARSE_CORRECTION_DEFLATED`, `PC_HPDDM_COARSE_CORRECTION_ADDITIVE`, `PC_HPDDM_COARSE_CORRECTION_BALANCED`, or `PC_HPDDM_COARSE_CORRECTION_NONE` 2627 2628 Options Database Key: 2629 . -pc_hpddm_coarse_correction <deflated, additive, balanced, none> - type of coarse correction to apply 2630 2631 Level: intermediate 2632 2633 .seealso: [](ch_ksp), `PCHPDDMGetCoarseCorrectionType()`, `PCHPDDM`, `PCHPDDMCoarseCorrectionType` 2634 @*/ 2635 PetscErrorCode PCHPDDMSetCoarseCorrectionType(PC pc, PCHPDDMCoarseCorrectionType type) 2636 { 2637 PetscFunctionBegin; 2638 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2639 PetscValidLogicalCollectiveEnum(pc, type, 2); 2640 PetscTryMethod(pc, "PCHPDDMSetCoarseCorrectionType_C", (PC, PCHPDDMCoarseCorrectionType), (pc, type)); 2641 PetscFunctionReturn(PETSC_SUCCESS); 2642 } 2643 2644 /*@ 2645 PCHPDDMGetCoarseCorrectionType - Gets the coarse correction type. 2646 2647 Input Parameter: 2648 . pc - preconditioner context 2649 2650 Output Parameter: 2651 . type - `PC_HPDDM_COARSE_CORRECTION_DEFLATED`, `PC_HPDDM_COARSE_CORRECTION_ADDITIVE`, `PC_HPDDM_COARSE_CORRECTION_BALANCED`, or `PC_HPDDM_COARSE_CORRECTION_NONE` 2652 2653 Level: intermediate 2654 2655 .seealso: [](ch_ksp), `PCHPDDMSetCoarseCorrectionType()`, `PCHPDDM`, `PCHPDDMCoarseCorrectionType` 2656 @*/ 2657 PetscErrorCode PCHPDDMGetCoarseCorrectionType(PC pc, PCHPDDMCoarseCorrectionType *type) 2658 { 2659 PetscFunctionBegin; 2660 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2661 if (type) { 2662 PetscAssertPointer(type, 2); 2663 PetscUseMethod(pc, "PCHPDDMGetCoarseCorrectionType_C", (PC, PCHPDDMCoarseCorrectionType *), (pc, type)); 2664 } 2665 PetscFunctionReturn(PETSC_SUCCESS); 2666 } 2667 2668 static PetscErrorCode PCHPDDMSetCoarseCorrectionType_HPDDM(PC pc, PCHPDDMCoarseCorrectionType type) 2669 { 2670 PC_HPDDM *data = (PC_HPDDM *)pc->data; 2671 2672 PetscFunctionBegin; 2673 data->correction = type; 2674 PetscFunctionReturn(PETSC_SUCCESS); 2675 } 2676 2677 static PetscErrorCode PCHPDDMGetCoarseCorrectionType_HPDDM(PC pc, PCHPDDMCoarseCorrectionType *type) 2678 { 2679 PC_HPDDM *data = (PC_HPDDM *)pc->data; 2680 2681 PetscFunctionBegin; 2682 *type = data->correction; 2683 PetscFunctionReturn(PETSC_SUCCESS); 2684 } 2685 2686 /*@ 2687 PCHPDDMSetSTShareSubKSP - Sets whether the `KSP` in SLEPc `ST` and the fine-level subdomain solver should be shared. 2688 2689 Input Parameters: 2690 + pc - preconditioner context 2691 - share - whether the `KSP` should be shared or not 2692 2693 Note: 2694 This is not the same as `PCSetReusePreconditioner()`. Given certain conditions (visible using -info), a symbolic factorization can be skipped 2695 when using a subdomain `PCType` such as `PCLU` or `PCCHOLESKY`. 2696 2697 Level: advanced 2698 2699 .seealso: [](ch_ksp), `PCHPDDM`, `PCHPDDMGetSTShareSubKSP()` 2700 @*/ 2701 PetscErrorCode PCHPDDMSetSTShareSubKSP(PC pc, PetscBool share) 2702 { 2703 PetscFunctionBegin; 2704 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2705 PetscTryMethod(pc, "PCHPDDMSetSTShareSubKSP_C", (PC, PetscBool), (pc, share)); 2706 PetscFunctionReturn(PETSC_SUCCESS); 2707 } 2708 2709 /*@ 2710 PCHPDDMGetSTShareSubKSP - Gets whether the `KSP` in SLEPc `ST` and the fine-level subdomain solver is shared. 2711 2712 Input Parameter: 2713 . pc - preconditioner context 2714 2715 Output Parameter: 2716 . share - whether the `KSP` is shared or not 2717 2718 Note: 2719 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 2720 when using a subdomain `PCType` such as `PCLU` or `PCCHOLESKY`. 2721 2722 Level: advanced 2723 2724 .seealso: [](ch_ksp), `PCHPDDM`, `PCHPDDMSetSTShareSubKSP()` 2725 @*/ 2726 PetscErrorCode PCHPDDMGetSTShareSubKSP(PC pc, PetscBool *share) 2727 { 2728 PetscFunctionBegin; 2729 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2730 if (share) { 2731 PetscAssertPointer(share, 2); 2732 PetscUseMethod(pc, "PCHPDDMGetSTShareSubKSP_C", (PC, PetscBool *), (pc, share)); 2733 } 2734 PetscFunctionReturn(PETSC_SUCCESS); 2735 } 2736 2737 static PetscErrorCode PCHPDDMSetSTShareSubKSP_HPDDM(PC pc, PetscBool share) 2738 { 2739 PC_HPDDM *data = (PC_HPDDM *)pc->data; 2740 2741 PetscFunctionBegin; 2742 data->share = share; 2743 PetscFunctionReturn(PETSC_SUCCESS); 2744 } 2745 2746 static PetscErrorCode PCHPDDMGetSTShareSubKSP_HPDDM(PC pc, PetscBool *share) 2747 { 2748 PC_HPDDM *data = (PC_HPDDM *)pc->data; 2749 2750 PetscFunctionBegin; 2751 *share = data->share; 2752 PetscFunctionReturn(PETSC_SUCCESS); 2753 } 2754 2755 /*@ 2756 PCHPDDMSetDeflationMat - Sets the deflation space used to assemble a coarser operator. 2757 2758 Input Parameters: 2759 + pc - preconditioner context 2760 . is - index set of the local deflation matrix 2761 - U - deflation sequential matrix stored as a `MATSEQDENSE` 2762 2763 Level: advanced 2764 2765 .seealso: [](ch_ksp), `PCHPDDM`, `PCDeflationSetSpace()`, `PCMGSetRestriction()` 2766 @*/ 2767 PetscErrorCode PCHPDDMSetDeflationMat(PC pc, IS is, Mat U) 2768 { 2769 PetscFunctionBegin; 2770 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2771 PetscValidHeaderSpecific(is, IS_CLASSID, 2); 2772 PetscValidHeaderSpecific(U, MAT_CLASSID, 3); 2773 PetscTryMethod(pc, "PCHPDDMSetDeflationMat_C", (PC, IS, Mat), (pc, is, U)); 2774 PetscFunctionReturn(PETSC_SUCCESS); 2775 } 2776 2777 static PetscErrorCode PCHPDDMSetDeflationMat_HPDDM(PC pc, IS is, Mat U) 2778 { 2779 PetscFunctionBegin; 2780 PetscCall(PCHPDDMSetAuxiliaryMat_Private(pc, is, U, PETSC_TRUE)); 2781 PetscFunctionReturn(PETSC_SUCCESS); 2782 } 2783 2784 PetscErrorCode HPDDMLoadDL_Private(PetscBool *found) 2785 { 2786 PetscBool flg; 2787 char lib[PETSC_MAX_PATH_LEN], dlib[PETSC_MAX_PATH_LEN], dir[PETSC_MAX_PATH_LEN]; 2788 2789 PetscFunctionBegin; 2790 PetscAssertPointer(found, 1); 2791 PetscCall(PetscStrncpy(dir, "${PETSC_LIB_DIR}", sizeof(dir))); 2792 PetscCall(PetscOptionsGetString(nullptr, nullptr, "-hpddm_dir", dir, sizeof(dir), nullptr)); 2793 PetscCall(PetscSNPrintf(lib, sizeof(lib), "%s/libhpddm_petsc", dir)); 2794 PetscCall(PetscDLLibraryRetrieve(PETSC_COMM_SELF, lib, dlib, 1024, found)); 2795 #if defined(SLEPC_LIB_DIR) /* this variable is passed during SLEPc ./configure when PETSc has not been configured */ 2796 if (!*found) { /* with --download-hpddm since slepcconf.h is not yet built (and thus can't be included) */ 2797 PetscCall(PetscStrncpy(dir, HPDDM_STR(SLEPC_LIB_DIR), sizeof(dir))); 2798 PetscCall(PetscSNPrintf(lib, sizeof(lib), "%s/libhpddm_petsc", dir)); 2799 PetscCall(PetscDLLibraryRetrieve(PETSC_COMM_SELF, lib, dlib, 1024, found)); 2800 } 2801 #endif 2802 if (!*found) { /* probable options for this to evaluate to PETSC_TRUE: system inconsistency (libhpddm_petsc moved by user?) or PETSc configured without --download-slepc */ 2803 PetscCall(PetscOptionsGetenv(PETSC_COMM_SELF, "SLEPC_DIR", dir, sizeof(dir), &flg)); 2804 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 */ 2805 PetscCall(PetscSNPrintf(lib, sizeof(lib), "%s/lib/libslepc", dir)); 2806 PetscCall(PetscDLLibraryRetrieve(PETSC_COMM_SELF, lib, dlib, 1024, found)); 2807 PetscCheck(*found, PETSC_COMM_SELF, PETSC_ERR_PLIB, "%s not found but SLEPC_DIR=%s", lib, dir); 2808 PetscCall(PetscDLLibraryAppend(PETSC_COMM_SELF, &PetscDLLibrariesLoaded, dlib)); 2809 PetscCall(PetscSNPrintf(lib, sizeof(lib), "%s/lib/libhpddm_petsc", dir)); /* libhpddm_petsc is always in the same directory as libslepc */ 2810 PetscCall(PetscDLLibraryRetrieve(PETSC_COMM_SELF, lib, dlib, 1024, found)); 2811 } 2812 } 2813 PetscCheck(*found, PETSC_COMM_SELF, PETSC_ERR_PLIB, "%s not found", lib); 2814 PetscCall(PetscDLLibraryAppend(PETSC_COMM_SELF, &PetscDLLibrariesLoaded, dlib)); 2815 PetscFunctionReturn(PETSC_SUCCESS); 2816 } 2817 2818 /*MC 2819 PCHPDDM - Interface with the HPDDM library. 2820 2821 This `PC` may be used to build multilevel spectral domain decomposition methods based on the GenEO framework {cite}`spillane2011robust` {cite}`al2021multilevel`. 2822 It may be viewed as an alternative to spectral 2823 AMGe or `PCBDDC` with adaptive selection of constraints. The interface is explained in details in {cite}`jolivetromanzampini2020` 2824 2825 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`). 2826 2827 For multilevel preconditioning, when using an assembled or hierarchical Pmat, one must provide an auxiliary local `Mat` (unassembled local operator for GenEO) using 2828 `PCHPDDMSetAuxiliaryMat()`. Calling this routine is not needed when using a `MATIS` Pmat, assembly is done internally using `MatConvert()`. 2829 2830 Options Database Keys: 2831 + -pc_hpddm_define_subdomains <true, default=false> - on the finest level, calls `PCASMSetLocalSubdomains()` with the `IS` supplied in `PCHPDDMSetAuxiliaryMat()` 2832 (not relevant with an unassembled Pmat) 2833 . -pc_hpddm_has_neumann <true, default=false> - on the finest level, informs the `PC` that the local Neumann matrix is supplied in `PCHPDDMSetAuxiliaryMat()` 2834 - -pc_hpddm_coarse_correction <type, default=deflated> - determines the `PCHPDDMCoarseCorrectionType` when calling `PCApply()` 2835 2836 Options for subdomain solvers, subdomain eigensolvers (for computing deflation vectors), and the coarse solver can be set using the following options database prefixes. 2837 .vb 2838 -pc_hpddm_levels_%d_pc_ 2839 -pc_hpddm_levels_%d_ksp_ 2840 -pc_hpddm_levels_%d_eps_ 2841 -pc_hpddm_levels_%d_p 2842 -pc_hpddm_levels_%d_mat_type 2843 -pc_hpddm_coarse_ 2844 -pc_hpddm_coarse_p 2845 -pc_hpddm_coarse_mat_type 2846 -pc_hpddm_coarse_mat_filter 2847 .ve 2848 2849 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 2850 -pc_hpddm_coarse_p 2 -pc_hpddm_coarse_mat_type baij will use 10 deflation vectors per subdomain on the fine "level 1", 2851 aggregate the fine subdomains into 4 "level 2" subdomains, then use 10 deflation vectors per subdomain on "level 2", 2852 and assemble the coarse matrix (of dimension 4 x 10 = 40) on two processes as a `MATBAIJ` (default is `MATSBAIJ`). 2853 2854 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. 2855 2856 Level: intermediate 2857 2858 Notes: 2859 This preconditioner requires that PETSc is built with SLEPc (``--download-slepc``). 2860 2861 By default, the underlying concurrent eigenproblems 2862 are solved using SLEPc shift-and-invert spectral transformation. This is usually what gives the best performance for GenEO, cf. 2863 {cite}`spillane2011robust` {cite}`jolivet2013scalabledd`. As 2864 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 2865 -pc_hpddm_levels_1_st_type sinvert. There are furthermore three options related to the (subdomain-wise local) eigensolver that are not described in 2866 SLEPc documentation since they are specific to `PCHPDDM`. 2867 .vb 2868 -pc_hpddm_levels_1_st_share_sub_ksp 2869 -pc_hpddm_levels_%d_eps_threshold 2870 -pc_hpddm_levels_1_eps_use_inertia 2871 .ve 2872 2873 The first option from the list only applies to the fine-level eigensolver, see `PCHPDDMSetSTShareSubKSP()`. The second option from the list is 2874 used to filter eigenmodes retrieved after convergence of `EPSSolve()` at "level N" such that eigenvectors used to define a "level N+1" coarse 2875 correction are associated to eigenvalues whose magnitude are lower or equal than -pc_hpddm_levels_N_eps_threshold. When using an `EPS` which cannot 2876 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 2877 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 2878 to supply -pc_hpddm_levels_1_eps_nev. This last option also only applies to the fine-level (N = 1) eigensolver. 2879 2880 See also {cite}`dolean2015introduction`, {cite}`al2022robust`, {cite}`al2022robustpd`, and {cite}`nataf2022recent` 2881 2882 .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCHPDDMSetAuxiliaryMat()`, `MATIS`, `PCBDDC`, `PCDEFLATION`, `PCTELESCOPE`, `PCASM`, 2883 `PCHPDDMSetCoarseCorrectionType()`, `PCHPDDMHasNeumannMat()`, `PCHPDDMSetRHSMat()`, `PCHPDDMSetDeflationMat()`, `PCHPDDMSetSTShareSubKSP()`, 2884 `PCHPDDMGetSTShareSubKSP()`, `PCHPDDMGetCoarseCorrectionType()`, `PCHPDDMGetComplexities()` 2885 M*/ 2886 PETSC_EXTERN PetscErrorCode PCCreate_HPDDM(PC pc) 2887 { 2888 PC_HPDDM *data; 2889 PetscBool found; 2890 2891 PetscFunctionBegin; 2892 if (!loadedSym) { 2893 PetscCall(HPDDMLoadDL_Private(&found)); 2894 if (found) PetscCall(PetscDLLibrarySym(PETSC_COMM_SELF, &PetscDLLibrariesLoaded, nullptr, "PCHPDDM_Internal", (void **)&loadedSym)); 2895 } 2896 PetscCheck(loadedSym, PETSC_COMM_SELF, PETSC_ERR_PLIB, "PCHPDDM_Internal symbol not found in loaded libhpddm_petsc"); 2897 PetscCall(PetscNew(&data)); 2898 pc->data = data; 2899 data->Neumann = PETSC_BOOL3_UNKNOWN; 2900 pc->ops->reset = PCReset_HPDDM; 2901 pc->ops->destroy = PCDestroy_HPDDM; 2902 pc->ops->setfromoptions = PCSetFromOptions_HPDDM; 2903 pc->ops->setup = PCSetUp_HPDDM; 2904 pc->ops->apply = PCApply_HPDDM; 2905 pc->ops->matapply = PCMatApply_HPDDM; 2906 pc->ops->view = PCView_HPDDM; 2907 pc->ops->presolve = PCPreSolve_HPDDM; 2908 2909 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetAuxiliaryMat_C", PCHPDDMSetAuxiliaryMat_HPDDM)); 2910 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMHasNeumannMat_C", PCHPDDMHasNeumannMat_HPDDM)); 2911 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetRHSMat_C", PCHPDDMSetRHSMat_HPDDM)); 2912 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetCoarseCorrectionType_C", PCHPDDMSetCoarseCorrectionType_HPDDM)); 2913 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMGetCoarseCorrectionType_C", PCHPDDMGetCoarseCorrectionType_HPDDM)); 2914 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetSTShareSubKSP_C", PCHPDDMSetSTShareSubKSP_HPDDM)); 2915 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMGetSTShareSubKSP_C", PCHPDDMGetSTShareSubKSP_HPDDM)); 2916 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetDeflationMat_C", PCHPDDMSetDeflationMat_HPDDM)); 2917 PetscFunctionReturn(PETSC_SUCCESS); 2918 } 2919 2920 /*@C 2921 PCHPDDMInitializePackage - This function initializes everything in the `PCHPDDM` package. It is called from `PCInitializePackage()`. 2922 2923 Level: developer 2924 2925 .seealso: [](ch_ksp), `PetscInitialize()` 2926 @*/ 2927 PetscErrorCode PCHPDDMInitializePackage(void) 2928 { 2929 char ename[32]; 2930 PetscInt i; 2931 2932 PetscFunctionBegin; 2933 if (PCHPDDMPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS); 2934 PCHPDDMPackageInitialized = PETSC_TRUE; 2935 PetscCall(PetscRegisterFinalize(PCHPDDMFinalizePackage)); 2936 /* general events registered once during package initialization */ 2937 /* some of these events are not triggered in libpetsc, */ 2938 /* but rather directly in libhpddm_petsc, */ 2939 /* which is in charge of performing the following operations */ 2940 2941 /* domain decomposition structure from Pmat sparsity pattern */ 2942 PetscCall(PetscLogEventRegister("PCHPDDMStrc", PC_CLASSID, &PC_HPDDM_Strc)); 2943 /* Galerkin product, redistribution, and setup (not triggered in libpetsc) */ 2944 PetscCall(PetscLogEventRegister("PCHPDDMPtAP", PC_CLASSID, &PC_HPDDM_PtAP)); 2945 /* Galerkin product with summation, redistribution, and setup (not triggered in libpetsc) */ 2946 PetscCall(PetscLogEventRegister("PCHPDDMPtBP", PC_CLASSID, &PC_HPDDM_PtBP)); 2947 /* next level construction using PtAP and PtBP (not triggered in libpetsc) */ 2948 PetscCall(PetscLogEventRegister("PCHPDDMNext", PC_CLASSID, &PC_HPDDM_Next)); 2949 static_assert(PETSC_PCHPDDM_MAXLEVELS <= 9, "PETSC_PCHPDDM_MAXLEVELS value is too high"); 2950 for (i = 1; i < PETSC_PCHPDDM_MAXLEVELS; ++i) { 2951 PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCHPDDMSetUp L%1" PetscInt_FMT, i)); 2952 /* events during a PCSetUp() at level #i _except_ the assembly */ 2953 /* of the Galerkin operator of the coarser level #(i + 1) */ 2954 PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_HPDDM_SetUp[i - 1])); 2955 PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCHPDDMSolve L%1" PetscInt_FMT, i)); 2956 /* events during a PCApply() at level #i _except_ */ 2957 /* the KSPSolve() of the coarser level #(i + 1) */ 2958 PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_HPDDM_Solve[i - 1])); 2959 } 2960 PetscFunctionReturn(PETSC_SUCCESS); 2961 } 2962 2963 /*@C 2964 PCHPDDMFinalizePackage - This function frees everything from the `PCHPDDM` package. It is called from `PetscFinalize()`. 2965 2966 Level: developer 2967 2968 .seealso: [](ch_ksp), `PetscFinalize()` 2969 @*/ 2970 PetscErrorCode PCHPDDMFinalizePackage(void) 2971 { 2972 PetscFunctionBegin; 2973 PCHPDDMPackageInitialized = PETSC_FALSE; 2974 PetscFunctionReturn(PETSC_SUCCESS); 2975 } 2976 2977 static PetscErrorCode MatMult_Harmonic(Mat A, Vec x, Vec y) 2978 { 2979 Harmonic h; /* [ A_00 A_01 ], furthermore, A_00 = [ A_loc,loc A_loc,ovl ], thus, A_01 = [ ] */ 2980 /* [ A_10 A_11 A_12 ] [ A_ovl,loc A_ovl,ovl ] [ A_ovl,1 ] */ 2981 Vec sub; /* y = A x = R_loc R_0 [ A_00 A_01 ]^-1 R_loc = [ I_loc ] */ 2982 /* [ A_10 A_11 ] R_1^T A_12 x [ ] */ 2983 PetscFunctionBegin; 2984 PetscCall(MatShellGetContext(A, &h)); 2985 PetscCall(VecSet(h->v, 0.0)); 2986 PetscCall(VecGetSubVector(h->v, h->is[0], &sub)); 2987 PetscCall(MatMult(h->A[0], x, sub)); 2988 PetscCall(VecRestoreSubVector(h->v, h->is[0], &sub)); 2989 PetscCall(KSPSolve(h->ksp, h->v, h->v)); 2990 PetscCall(VecISCopy(h->v, h->is[1], SCATTER_REVERSE, y)); 2991 PetscFunctionReturn(PETSC_SUCCESS); 2992 } 2993 2994 static PetscErrorCode MatMultTranspose_Harmonic(Mat A, Vec y, Vec x) 2995 { 2996 Harmonic h; /* x = A^T y = [ A_00 A_01 ]^-T R_0^T R_loc^T y */ 2997 Vec sub; /* A_12^T R_1 [ A_10 A_11 ] */ 2998 2999 PetscFunctionBegin; 3000 PetscCall(MatShellGetContext(A, &h)); 3001 PetscCall(VecSet(h->v, 0.0)); 3002 PetscCall(VecISCopy(h->v, h->is[1], SCATTER_FORWARD, y)); 3003 PetscCall(KSPSolveTranspose(h->ksp, h->v, h->v)); 3004 PetscCall(VecGetSubVector(h->v, h->is[0], &sub)); 3005 PetscCall(MatMultTranspose(h->A[0], sub, x)); 3006 PetscCall(VecRestoreSubVector(h->v, h->is[0], &sub)); 3007 PetscFunctionReturn(PETSC_SUCCESS); 3008 } 3009 3010 static PetscErrorCode MatProduct_AB_Harmonic(Mat S, Mat X, Mat Y, void *) 3011 { 3012 Harmonic h; 3013 Mat A, B; 3014 Vec a, b; 3015 3016 PetscFunctionBegin; 3017 PetscCall(MatShellGetContext(S, &h)); 3018 PetscCall(MatMatMult(h->A[0], X, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &A)); 3019 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, h->ksp->pc->mat->rmap->n, A->cmap->n, nullptr, &B)); 3020 for (PetscInt i = 0; i < A->cmap->n; ++i) { 3021 PetscCall(MatDenseGetColumnVecRead(A, i, &a)); 3022 PetscCall(MatDenseGetColumnVecWrite(B, i, &b)); 3023 PetscCall(VecISCopy(b, h->is[0], SCATTER_FORWARD, a)); 3024 PetscCall(MatDenseRestoreColumnVecWrite(B, i, &b)); 3025 PetscCall(MatDenseRestoreColumnVecRead(A, i, &a)); 3026 } 3027 PetscCall(MatDestroy(&A)); 3028 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, h->ksp->pc->mat->rmap->n, B->cmap->n, nullptr, &A)); 3029 PetscCall(KSPMatSolve(h->ksp, B, A)); 3030 PetscCall(MatDestroy(&B)); 3031 for (PetscInt i = 0; i < A->cmap->n; ++i) { 3032 PetscCall(MatDenseGetColumnVecRead(A, i, &a)); 3033 PetscCall(MatDenseGetColumnVecWrite(Y, i, &b)); 3034 PetscCall(VecISCopy(a, h->is[1], SCATTER_REVERSE, b)); 3035 PetscCall(MatDenseRestoreColumnVecWrite(Y, i, &b)); 3036 PetscCall(MatDenseRestoreColumnVecRead(A, i, &a)); 3037 } 3038 PetscCall(MatDestroy(&A)); 3039 PetscFunctionReturn(PETSC_SUCCESS); 3040 } 3041 3042 static PetscErrorCode MatDestroy_Harmonic(Mat A) 3043 { 3044 Harmonic h; 3045 3046 PetscFunctionBegin; 3047 PetscCall(MatShellGetContext(A, &h)); 3048 for (PetscInt i = 0; i < (h->A[1] ? 5 : 3); ++i) PetscCall(ISDestroy(h->is + i)); 3049 PetscCall(PetscFree(h->is)); 3050 PetscCall(VecDestroy(&h->v)); 3051 for (PetscInt i = 0; i < 2; ++i) PetscCall(MatDestroy(h->A + i)); 3052 PetscCall(PetscFree(h->A)); 3053 PetscCall(KSPDestroy(&h->ksp)); 3054 PetscCall(PetscFree(h)); 3055 PetscFunctionReturn(PETSC_SUCCESS); 3056 } 3057