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