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