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