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