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