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