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