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