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