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