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