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