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