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 PetscCall(MatSetOption(C, MAT_SYMMETRIC, PETSC_TRUE)); 2222 break; 2223 default: 2224 PetscCall(MatConvert(P, MATMPIAIJ, MAT_INITIAL_MATRIX, &C)); 2225 } 2226 PetscCall(MatISGetLocalToGlobalMapping(P, &l2g, nullptr)); 2227 PetscCall(PetscObjectReference((PetscObject)P)); 2228 PetscCall(KSPSetOperators(data->levels[0]->ksp, A, C)); 2229 std::swap(C, P); 2230 PetscCall(ISLocalToGlobalMappingGetSize(l2g, &n)); 2231 PetscCall(ISCreateStride(PETSC_COMM_SELF, n, 0, 1, &loc)); 2232 PetscCall(ISLocalToGlobalMappingApplyIS(l2g, loc, &is[0])); 2233 PetscCall(ISDestroy(&loc)); 2234 /* the auxiliary Mat is _not_ the local Neumann matrix */ 2235 /* it is the local Neumann matrix augmented (with zeros) through MatIncreaseOverlap() */ 2236 data->Neumann = PETSC_BOOL3_FALSE; 2237 structure = SAME_NONZERO_PATTERN; 2238 } else { 2239 is[0] = data->is; 2240 if (algebraic || ctx) subdomains = PETSC_TRUE; 2241 PetscCall(PetscOptionsGetBool(nullptr, pcpre, "-pc_hpddm_define_subdomains", &subdomains, nullptr)); 2242 if (ctx) PetscCheck(subdomains, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-%spc_hpddm_schur_precondition geneo and -%spc_hpddm_define_subdomains false", pcpre, pcpre); 2243 if (PetscBool3ToBool(data->Neumann)) { 2244 PetscCheck(!block, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-pc_hpddm_block_splitting and -pc_hpddm_has_neumann"); 2245 PetscCheck(overlap == -1, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-pc_hpddm_harmonic_overlap %" PetscInt_FMT " and -pc_hpddm_has_neumann", overlap); 2246 PetscCheck(!algebraic, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-pc_hpddm_levels_1_st_pc_type mat and -pc_hpddm_has_neumann"); 2247 } 2248 if (PetscBool3ToBool(data->Neumann) || block) structure = SAME_NONZERO_PATTERN; 2249 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)data->is), P->rmap->n, P->rmap->rstart, 1, &loc)); 2250 } 2251 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_levels_1_", pcpre ? pcpre : "")); 2252 PetscCall(PetscOptionsGetEnum(nullptr, prefix, "-st_matstructure", MatStructures, (PetscEnum *)&structure, &flg)); /* if not user-provided, force its value when possible */ 2253 if (!flg && structure == SAME_NONZERO_PATTERN) { /* cannot call STSetMatStructure() yet, insert the appropriate option in the database, parsed by STSetFromOptions() */ 2254 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-%spc_hpddm_levels_1_st_matstructure", pcpre ? pcpre : "")); 2255 PetscCall(PetscOptionsSetValue(nullptr, prefix, MatStructures[structure])); 2256 } 2257 flg = PETSC_FALSE; 2258 if (data->share) { 2259 data->share = PETSC_FALSE; /* will be reset to PETSC_TRUE if none of the conditions below are true */ 2260 if (!subdomains) PetscCall(PetscInfo(pc, "Cannot share subdomain KSP between SLEPc and PETSc since -%spc_hpddm_define_subdomains is not true\n", pcpre ? pcpre : "")); 2261 else if (data->deflation) PetscCall(PetscInfo(pc, "Nothing to share since PCHPDDMSetDeflationMat() has been called\n")); 2262 else if (ismatis) PetscCall(PetscInfo(pc, "Cannot share subdomain KSP between SLEPc and PETSc with a Pmat of type MATIS\n")); 2263 else if (!algebraic && structure != SAME_NONZERO_PATTERN) 2264 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])); 2265 else data->share = PETSC_TRUE; 2266 } 2267 if (!ismatis) { 2268 if (data->share || (!PetscBool3ToBool(data->Neumann) && subdomains)) PetscCall(ISDuplicate(is[0], &unsorted)); 2269 else unsorted = is[0]; 2270 } 2271 if ((ctx || data->N > 1) && (data->aux || ismatis || algebraic)) { 2272 PetscCheck(loadedSym, PETSC_COMM_SELF, PETSC_ERR_PLIB, "HPDDM library not loaded, cannot use more than one level"); 2273 PetscCall(MatSetOption(P, MAT_SUBMAT_SINGLEIS, PETSC_TRUE)); 2274 if (ismatis) { 2275 /* needed by HPDDM (currently) so that the partition of unity is 0 on subdomain interfaces */ 2276 PetscCall(MatIncreaseOverlap(P, 1, is, 1)); 2277 PetscCall(ISDestroy(&data->is)); 2278 data->is = is[0]; 2279 } else { 2280 if (PetscDefined(USE_DEBUG)) PetscCall(PCHPDDMCheckInclusion_Private(pc, data->is, loc, PETSC_TRUE)); 2281 if (!ctx && overlap == -1) PetscCall(PetscObjectComposeFunction((PetscObject)pc->pmat, "PCHPDDMAlgebraicAuxiliaryMat_Private_C", PCHPDDMAlgebraicAuxiliaryMat_Private)); 2282 if (!PetscBool3ToBool(data->Neumann) && (!algebraic || overlap != -1)) { 2283 PetscCall(PetscObjectTypeCompare((PetscObject)P, MATMPISBAIJ, &flg)); 2284 if (flg) { 2285 /* maybe better to ISSort(is[0]), MatCreateSubMatrices(), and then MatPermute() */ 2286 /* but there is no MatPermute_SeqSBAIJ(), so as before, just use MATMPIBAIJ */ 2287 PetscCall(MatConvert(P, MATMPIBAIJ, MAT_INITIAL_MATRIX, &uaux)); 2288 flg = PETSC_FALSE; 2289 } 2290 } 2291 } 2292 if (algebraic && overlap == -1) { 2293 PetscUseMethod(pc->pmat, "PCHPDDMAlgebraicAuxiliaryMat_Private_C", (Mat, IS *, Mat *[], PetscBool), (P, is, &sub, block)); 2294 if (block) { 2295 PetscCall(PetscObjectQuery((PetscObject)sub[0], "_PCHPDDM_Neumann_Mat", (PetscObject *)&data->aux)); 2296 PetscCall(PetscObjectCompose((PetscObject)sub[0], "_PCHPDDM_Neumann_Mat", nullptr)); 2297 } 2298 } else if (!uaux || overlap != -1) { 2299 if (!ctx) { 2300 if (PetscBool3ToBool(data->Neumann)) sub = &data->aux; 2301 else { 2302 PetscBool flg; 2303 if (overlap != -1) { 2304 Harmonic h; 2305 Mat A0, *a; /* with an SVD: [ A_00 A_01 ] */ 2306 IS ov[2], rows, cols, stride; /* [ A_10 A_11 A_12 ] */ 2307 const PetscInt *i[2], bs = P->cmap->bs; /* with a GEVP: [ A_00 A_01 ] */ 2308 PetscInt n[2]; /* [ A_10 A_11 A_12 ] */ 2309 std::vector<PetscInt> v[2]; /* [ A_21 A_22 ] */ 2310 2311 do { 2312 PetscCall(ISDuplicate(data->is, ov)); 2313 if (overlap > 1) PetscCall(MatIncreaseOverlap(P, 1, ov, overlap - 1)); 2314 PetscCall(ISDuplicate(ov[0], ov + 1)); 2315 PetscCall(MatIncreaseOverlap(P, 1, ov + 1, 1)); 2316 PetscCall(ISGetLocalSize(ov[0], n)); 2317 PetscCall(ISGetLocalSize(ov[1], n + 1)); 2318 flg = PetscBool(n[0] == n[1] && n[0] != P->rmap->n); 2319 PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &flg, 1, MPI_C_BOOL, MPI_LOR, PetscObjectComm((PetscObject)pc))); 2320 if (flg) { 2321 PetscCall(ISDestroy(ov)); 2322 PetscCall(ISDestroy(ov + 1)); 2323 PetscCheck(--overlap, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "No oversampling possible"); 2324 PetscCall(PetscInfo(pc, "Supplied -%spc_hpddm_harmonic_overlap parameter is too large, it has been decreased to %" PetscInt_FMT "\n", pcpre ? pcpre : "", overlap)); 2325 } else break; 2326 } while (1); 2327 PetscCall(PetscNew(&h)); 2328 h->ksp = nullptr; 2329 PetscCall(PetscCalloc1(2, &h->A)); 2330 PetscCall(PetscOptionsHasName(nullptr, prefix, "-eps_nev", &flg)); 2331 if (!flg) { 2332 PetscCall(PetscOptionsHasName(nullptr, prefix, "-svd_nsv", &flg)); 2333 if (!flg) PetscCall(PetscOptionsHasName(nullptr, prefix, "-svd_threshold_relative", &flg)); 2334 } else flg = PETSC_FALSE; 2335 PetscCall(ISSort(ov[0])); 2336 if (!flg) PetscCall(ISSort(ov[1])); 2337 PetscCall(PetscCalloc1(5, &h->is)); 2338 PetscCall(MatCreateSubMatrices(uaux ? uaux : P, 1, ov + !flg, ov + 1, MAT_INITIAL_MATRIX, &a)); /* submatrix from above, either square (!flg) or rectangular (flg) */ 2339 for (PetscInt j = 0; j < 2; ++j) PetscCall(ISGetIndices(ov[j], i + j)); 2340 v[1].reserve((n[1] - n[0]) / bs); 2341 for (PetscInt j = 0; j < n[1]; j += bs) { /* indices of the (2,2) block */ 2342 PetscInt location; 2343 PetscCall(ISLocate(ov[0], i[1][j], &location)); 2344 if (location < 0) v[1].emplace_back(j / bs); 2345 } 2346 if (!flg) { 2347 h->A[1] = a[0]; 2348 PetscCall(PetscObjectReference((PetscObject)h->A[1])); 2349 v[0].reserve((n[0] - P->rmap->n) / bs); 2350 for (PetscInt j = 0; j < n[1]; j += bs) { /* row indices of the (1,2) block */ 2351 PetscInt location; 2352 PetscCall(ISLocate(loc, i[1][j], &location)); 2353 if (location < 0) { 2354 PetscCall(ISLocate(ov[0], i[1][j], &location)); 2355 if (location >= 0) v[0].emplace_back(j / bs); 2356 } 2357 } 2358 PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, v[0].size(), v[0].data(), PETSC_USE_POINTER, &rows)); 2359 PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, v[1].size(), v[1].data(), PETSC_COPY_VALUES, h->is + 4)); 2360 PetscCall(MatCreateSubMatrix(a[0], rows, h->is[4], MAT_INITIAL_MATRIX, h->A)); /* A_12 submatrix from above */ 2361 PetscCall(ISDestroy(&rows)); 2362 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 */ 2363 PetscCall(ISEmbed(ov[0], ov[1], PETSC_TRUE, &rows)); 2364 PetscCall(MatCreateSubMatrix(a[0], rows, cols = rows, MAT_INITIAL_MATRIX, &A0)); /* [ A_00 A_01 ; A_10 A_11 ] submatrix from above */ 2365 PetscCall(ISDestroy(&rows)); 2366 v[0].clear(); 2367 PetscCall(ISEmbed(loc, ov[1], PETSC_TRUE, h->is + 3)); 2368 PetscCall(ISEmbed(data->is, ov[1], PETSC_TRUE, h->is + 2)); 2369 } 2370 v[0].reserve((n[0] - P->rmap->n) / bs); 2371 for (PetscInt j = 0; j < n[0]; j += bs) { 2372 PetscInt location; 2373 PetscCall(ISLocate(loc, i[0][j], &location)); 2374 if (location < 0) v[0].emplace_back(j / bs); 2375 } 2376 PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, v[0].size(), v[0].data(), PETSC_USE_POINTER, &rows)); 2377 for (PetscInt j = 0; j < 2; ++j) PetscCall(ISRestoreIndices(ov[j], i + j)); 2378 if (flg) { 2379 IS is; 2380 PetscCall(ISCreateStride(PETSC_COMM_SELF, a[0]->rmap->n, 0, 1, &is)); 2381 PetscCall(ISEmbed(ov[0], ov[1], PETSC_TRUE, &cols)); 2382 PetscCall(MatCreateSubMatrix(a[0], is, cols, MAT_INITIAL_MATRIX, &A0)); /* [ A_00 A_01 ; A_10 A_11 ] submatrix from above */ 2383 PetscCall(ISDestroy(&cols)); 2384 PetscCall(ISDestroy(&is)); 2385 if (uaux) PetscCall(MatConvert(A0, MATSEQSBAIJ, MAT_INPLACE_MATRIX, &A0)); /* initial Pmat was MATSBAIJ, convert back to the same format since this submatrix is square */ 2386 PetscCall(ISEmbed(loc, data->is, PETSC_TRUE, h->is + 2)); 2387 PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, v[1].size(), v[1].data(), PETSC_USE_POINTER, &cols)); 2388 PetscCall(MatCreateSubMatrix(a[0], rows, cols, MAT_INITIAL_MATRIX, h->A)); /* A_12 submatrix from above */ 2389 PetscCall(ISDestroy(&cols)); 2390 } 2391 PetscCall(ISCreateStride(PETSC_COMM_SELF, A0->rmap->n, 0, 1, &stride)); 2392 PetscCall(ISEmbed(rows, stride, PETSC_TRUE, h->is)); 2393 PetscCall(ISDestroy(&stride)); 2394 PetscCall(ISDestroy(&rows)); 2395 PetscCall(ISEmbed(loc, ov[0], PETSC_TRUE, h->is + 1)); 2396 if (subdomains) { 2397 if (!data->levels[0]->pc) { 2398 PetscCall(PCCreate(PetscObjectComm((PetscObject)pc), &data->levels[0]->pc)); 2399 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_levels_1_", pcpre ? pcpre : "")); 2400 PetscCall(PCSetOptionsPrefix(data->levels[0]->pc, prefix)); 2401 PetscCall(PCSetOperators(data->levels[0]->pc, A, P)); 2402 } 2403 PetscCall(PCSetType(data->levels[0]->pc, PCASM)); 2404 if (!data->levels[0]->pc->setupcalled) PetscCall(PCASMSetLocalSubdomains(data->levels[0]->pc, 1, ov + !flg, &loc)); 2405 PetscCall(PCHPDDMCommunicationAvoidingPCASM_Private(data->levels[0]->pc, flg ? A0 : a[0], PETSC_TRUE)); 2406 if (!flg) ++overlap; 2407 if (data->share) { 2408 PetscInt n = -1; 2409 PetscTryMethod(data->levels[0]->pc, "PCASMGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (data->levels[0]->pc, &n, nullptr, &ksp)); 2410 PetscCheck(n == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of subdomain solver %" PetscInt_FMT " != 1", n); 2411 if (flg) { 2412 h->ksp = ksp[0]; 2413 PetscCall(PetscObjectReference((PetscObject)h->ksp)); 2414 } 2415 } 2416 } 2417 if (!h->ksp) { 2418 PetscBool share = data->share; 2419 PetscCall(KSPCreate(PETSC_COMM_SELF, &h->ksp)); 2420 PetscCall(KSPSetType(h->ksp, KSPPREONLY)); 2421 PetscCall(KSPSetOperators(h->ksp, A0, A0)); 2422 do { 2423 if (!data->share) { 2424 share = PETSC_FALSE; 2425 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_levels_1_%s", pcpre ? pcpre : "", flg ? "svd_" : "eps_")); 2426 PetscCall(KSPSetOptionsPrefix(h->ksp, prefix)); 2427 PetscCall(KSPSetFromOptions(h->ksp)); 2428 } else { 2429 MatSolverType type; 2430 PetscCall(KSPGetPC(ksp[0], &pc)); 2431 PetscCall(PetscObjectTypeCompareAny((PetscObject)pc, &data->share, PCLU, PCCHOLESKY, "")); 2432 if (data->share) { 2433 PetscCall(PCFactorGetMatSolverType(pc, &type)); 2434 if (!type) { 2435 if (PetscDefined(HAVE_MUMPS)) PetscCall(PCFactorSetMatSolverType(pc, MATSOLVERMUMPS)); 2436 else if (PetscDefined(HAVE_MKL_PARDISO)) PetscCall(PCFactorSetMatSolverType(pc, MATSOLVERMKL_PARDISO)); 2437 else data->share = PETSC_FALSE; 2438 if (data->share) PetscCall(PCSetFromOptions(pc)); 2439 } else { 2440 PetscCall(PetscStrcmp(type, MATSOLVERMUMPS, &data->share)); 2441 if (!data->share) PetscCall(PetscStrcmp(type, MATSOLVERMKL_PARDISO, &data->share)); 2442 } 2443 if (data->share) { 2444 std::tuple<KSP, IS, Vec[2]> *p; 2445 PetscCall(PCFactorGetMatrix(pc, &A)); 2446 PetscCall(MatFactorSetSchurIS(A, h->is[4])); 2447 PetscCall(KSPSetUp(ksp[0])); 2448 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_levels_1_eps_shell_", pcpre ? pcpre : "")); 2449 PetscCall(KSPSetOptionsPrefix(h->ksp, prefix)); 2450 PetscCall(KSPSetFromOptions(h->ksp)); 2451 PetscCall(KSPGetPC(h->ksp, &pc)); 2452 PetscCall(PCSetType(pc, PCSHELL)); 2453 PetscCall(PetscNew(&p)); 2454 std::get<0>(*p) = ksp[0]; 2455 PetscCall(ISEmbed(ov[0], ov[1], PETSC_TRUE, &std::get<1>(*p))); 2456 PetscCall(MatCreateVecs(A, std::get<2>(*p), std::get<2>(*p) + 1)); 2457 PetscCall(PCShellSetContext(pc, p)); 2458 PetscCall(PCShellSetApply(pc, PCApply_Schur)); 2459 PetscCall(PCShellSetApplyTranspose(pc, PCApply_Schur<Vec, true>)); 2460 PetscCall(PCShellSetMatApply(pc, PCApply_Schur<Mat>)); 2461 PetscCall(PCShellSetDestroy(pc, PCDestroy_Schur)); 2462 } 2463 } 2464 if (!data->share) PetscCall(PetscInfo(pc, "Cannot share subdomain KSP between SLEPc and PETSc since neither MUMPS nor MKL PARDISO is used\n")); 2465 } 2466 } 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 */ 2467 } 2468 PetscCall(ISDestroy(ov)); 2469 PetscCall(ISDestroy(ov + 1)); 2470 if (overlap == 1 && subdomains && flg) { 2471 *subA = A0; 2472 sub = subA; 2473 if (uaux) PetscCall(MatDestroy(&uaux)); 2474 } else PetscCall(MatDestroy(&A0)); 2475 PetscCall(MatCreateShell(PETSC_COMM_SELF, P->rmap->n, n[1] - n[0], P->rmap->n, n[1] - n[0], h, &data->aux)); 2476 PetscCall(KSPSetErrorIfNotConverged(h->ksp, PETSC_TRUE)); /* bail out as early as possible to avoid (apparently) unrelated error messages */ 2477 PetscCall(MatCreateVecs(h->ksp->pc->pmat, &h->v, nullptr)); 2478 PetscCall(MatShellSetOperation(data->aux, MATOP_MULT, (PetscErrorCodeFn *)MatMult_Harmonic)); 2479 PetscCall(MatShellSetOperation(data->aux, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)MatMultTranspose_Harmonic)); 2480 PetscCall(MatShellSetMatProductOperation(data->aux, MATPRODUCT_AB, nullptr, MatProduct_AB_Harmonic, nullptr, MATDENSE, MATDENSE)); 2481 PetscCall(MatShellSetMatProductOperation(data->aux, MATPRODUCT_AtB, nullptr, MatProduct_AtB_Harmonic, nullptr, MATDENSE, MATDENSE)); 2482 PetscCall(MatShellSetOperation(data->aux, MATOP_DESTROY, (PetscErrorCodeFn *)MatDestroy_Harmonic)); 2483 PetscCall(MatDestroySubMatrices(1, &a)); 2484 } 2485 if (overlap != 1 || !subdomains) { 2486 PetscCall(MatCreateSubMatrices(uaux ? uaux : P, 1, is, is, MAT_INITIAL_MATRIX, &sub)); 2487 if (ismatis) { 2488 PetscCall(MatISGetLocalMat(C, &N)); 2489 PetscCall(PetscObjectTypeCompare((PetscObject)N, MATSEQSBAIJ, &flg)); 2490 if (flg) PetscCall(MatConvert(sub[0], MATSEQSBAIJ, MAT_INPLACE_MATRIX, sub)); 2491 PetscCall(MatISRestoreLocalMat(C, &N)); 2492 } 2493 } 2494 if (uaux) { 2495 PetscCall(MatDestroy(&uaux)); 2496 PetscCall(MatConvert(sub[0], MATSEQSBAIJ, MAT_INPLACE_MATRIX, sub)); 2497 } 2498 } 2499 } 2500 } else if (!ctx) { 2501 PetscCall(MatCreateSubMatrices(uaux, 1, is, is, MAT_INITIAL_MATRIX, &sub)); 2502 PetscCall(MatDestroy(&uaux)); 2503 PetscCall(MatConvert(sub[0], MATSEQSBAIJ, MAT_INPLACE_MATRIX, sub)); 2504 } 2505 if (data->N > 1) { 2506 /* Vec holding the partition of unity */ 2507 if (!data->levels[0]->D) { 2508 PetscCall(ISGetLocalSize(data->is, &n)); 2509 PetscCall(VecCreateMPI(PETSC_COMM_SELF, n, PETSC_DETERMINE, &data->levels[0]->D)); 2510 } 2511 if (data->share && overlap == -1) { 2512 Mat D; 2513 IS perm = nullptr; 2514 PetscInt size = -1; 2515 2516 if (!data->levels[0]->pc) { 2517 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_levels_1_", pcpre ? pcpre : "")); 2518 PetscCall(PCCreate(PetscObjectComm((PetscObject)pc), &data->levels[0]->pc)); 2519 PetscCall(PCSetOptionsPrefix(data->levels[0]->pc, prefix)); 2520 PetscCall(PCSetOperators(data->levels[0]->pc, A, P)); 2521 } 2522 PetscCall(PCSetType(data->levels[0]->pc, PCASM)); 2523 if (!ctx) { 2524 if (!data->levels[0]->pc->setupcalled) { 2525 IS sorted; /* PCASM will sort the input IS, duplicate it to return an unmodified (PCHPDDM) input IS */ 2526 PetscCall(ISDuplicate(is[0], &sorted)); 2527 PetscCall(PCASMSetLocalSubdomains(data->levels[0]->pc, 1, &sorted, &loc)); 2528 PetscCall(PetscObjectDereference((PetscObject)sorted)); 2529 } 2530 PetscCall(PCSetFromOptions(data->levels[0]->pc)); 2531 if (block) { 2532 PetscCall(PCHPDDMPermute_Private(unsorted, data->is, &uis, sub[0], &C, &perm)); 2533 PetscCall(PCHPDDMCommunicationAvoidingPCASM_Private(data->levels[0]->pc, C, algebraic)); 2534 } else PetscCall(PCSetUp(data->levels[0]->pc)); 2535 PetscTryMethod(data->levels[0]->pc, "PCASMGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (data->levels[0]->pc, &size, nullptr, &ksp)); 2536 if (size != 1) { 2537 data->share = PETSC_FALSE; 2538 PetscCheck(size == -1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of subdomain solver %" PetscInt_FMT " != 1", size); 2539 PetscCall(PetscInfo(pc, "Cannot share subdomain KSP between SLEPc and PETSc since PCASMGetSubKSP() not found in fine-level PC\n")); 2540 PetscCall(ISDestroy(&unsorted)); 2541 unsorted = is[0]; 2542 } else { 2543 const char *matpre; 2544 PetscBool cmp[4]; 2545 2546 if (!block && !ctx) PetscCall(PCHPDDMPermute_Private(unsorted, data->is, &uis, PetscBool3ToBool(data->Neumann) ? sub[0] : data->aux, &C, &perm)); 2547 if (perm) { /* unsorted input IS */ 2548 if (!PetscBool3ToBool(data->Neumann) && !block) { 2549 PetscCall(MatPermute(sub[0], perm, perm, &D)); /* permute since PCASM will call ISSort() */ 2550 PetscCall(MatHeaderReplace(sub[0], &D)); 2551 } 2552 if (data->B) { /* see PCHPDDMSetRHSMat() */ 2553 PetscCall(MatPermute(data->B, perm, perm, &D)); 2554 PetscCall(MatHeaderReplace(data->B, &D)); 2555 } 2556 PetscCall(ISDestroy(&perm)); 2557 } 2558 PetscCall(KSPGetOperators(ksp[0], subA, subA + 1)); 2559 PetscCall(PetscObjectReference((PetscObject)subA[0])); 2560 PetscCall(MatDuplicate(subA[1], MAT_SHARE_NONZERO_PATTERN, &D)); 2561 PetscCall(MatGetOptionsPrefix(subA[1], &matpre)); 2562 PetscCall(MatSetOptionsPrefix(D, matpre)); 2563 PetscCall(PetscObjectTypeCompare((PetscObject)D, MATNORMAL, cmp)); 2564 PetscCall(PetscObjectTypeCompare((PetscObject)C, MATNORMAL, cmp + 1)); 2565 if (!cmp[0]) PetscCall(PetscObjectTypeCompare((PetscObject)D, MATNORMALHERMITIAN, cmp + 2)); 2566 else cmp[2] = PETSC_FALSE; 2567 if (!cmp[1]) PetscCall(PetscObjectTypeCompare((PetscObject)C, MATNORMALHERMITIAN, cmp + 3)); 2568 else cmp[3] = PETSC_FALSE; 2569 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); 2570 if (!cmp[0] && !cmp[2]) { 2571 if (!block) { 2572 if (PetscDefined(USE_DEBUG)) PetscCall(PCHPDDMCheckMatStructure_Private(pc, D, C)); 2573 PetscCall(MatAXPY(D, 1.0, C, SUBSET_NONZERO_PATTERN)); 2574 } else { 2575 structure = DIFFERENT_NONZERO_PATTERN; 2576 PetscCall(MatAXPY(D, 1.0, data->aux, structure)); 2577 } 2578 } else { 2579 Mat mat[2]; 2580 2581 if (cmp[0]) { 2582 PetscCall(MatNormalGetMat(D, mat)); 2583 PetscCall(MatNormalGetMat(C, mat + 1)); 2584 } else { 2585 PetscCall(MatNormalHermitianGetMat(D, mat)); 2586 PetscCall(MatNormalHermitianGetMat(C, mat + 1)); 2587 } 2588 PetscCall(MatAXPY(mat[0], 1.0, mat[1], SUBSET_NONZERO_PATTERN)); 2589 } 2590 PetscCall(MatPropagateSymmetryOptions(C, D)); 2591 PetscCall(MatDestroy(&C)); 2592 C = D; 2593 /* swap pointers so that variables stay consistent throughout PCSetUp() */ 2594 std::swap(C, data->aux); 2595 std::swap(uis, data->is); 2596 swap = PETSC_TRUE; 2597 } 2598 } 2599 } 2600 } 2601 if (ctx) { 2602 PC_HPDDM *data_00 = (PC_HPDDM *)std::get<0>(*ctx)[0]->data; 2603 PC s; 2604 Mat A00, P00, A01 = nullptr, A10, A11, N, b[4]; 2605 IS sorted, is[2], *is_00; 2606 MatSolverType type; 2607 std::pair<PC, Vec[2]> *p; 2608 2609 n = -1; 2610 PetscTryMethod(data_00->levels[0]->pc, "PCASMGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (data_00->levels[0]->pc, &n, nullptr, &ksp)); 2611 PetscCheck(n == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of subdomain solver %" PetscInt_FMT " != 1", n); 2612 PetscCall(KSPGetOperators(ksp[0], subA, subA + 1)); 2613 PetscCall(ISGetLocalSize(data_00->is, &n)); 2614 if (n != subA[0]->rmap->n || n != subA[0]->cmap->n) { 2615 PetscCall(PCASMGetLocalSubdomains(data_00->levels[0]->pc, &n, &is_00, nullptr)); 2616 PetscCall(ISGetLocalSize(*is_00, &n)); 2617 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); 2618 } else is_00 = &data_00->is; 2619 PetscCall(PCHPDDMPermute_Private(unsorted, data->is, &uis, data->aux, &C, nullptr)); /* permute since PCASM works with a sorted IS */ 2620 std::swap(C, data->aux); 2621 std::swap(uis, data->is); 2622 swap = PETSC_TRUE; 2623 PetscCall(MatSchurComplementGetSubMatrices(P, &A00, &P00, std::get<1>(*ctx), &A10, &A11)); 2624 std::get<1>(*ctx)[1] = A10; 2625 PetscCall(PetscObjectTypeCompare((PetscObject)A10, MATTRANSPOSEVIRTUAL, &flg)); 2626 if (flg) PetscCall(MatTransposeGetMat(A10, &A01)); 2627 else { 2628 PetscBool flg; 2629 2630 PetscCall(PetscObjectTypeCompare((PetscObject)A10, MATHERMITIANTRANSPOSEVIRTUAL, &flg)); 2631 if (flg) PetscCall(MatHermitianTransposeGetMat(A10, &A01)); 2632 } 2633 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 */ 2634 PetscCall(ISSort(sorted)); /* this is to avoid changing users inputs, but it requires a new call to ISSort() here */ 2635 if (!A01) { 2636 PetscCall(MatSetOption(A10, MAT_SUBMAT_SINGLEIS, PETSC_TRUE)); 2637 PetscCall(MatCreateSubMatrices(A10, 1, &data->is, &sorted, MAT_INITIAL_MATRIX, &sub)); 2638 b[2] = sub[0]; 2639 PetscCall(PetscObjectReference((PetscObject)sub[0])); 2640 PetscCall(MatDestroySubMatrices(1, &sub)); 2641 PetscCall(PetscObjectTypeCompare((PetscObject)std::get<1>(*ctx)[0], MATTRANSPOSEVIRTUAL, &flg)); 2642 A10 = nullptr; 2643 if (flg) PetscCall(MatTransposeGetMat(std::get<1>(*ctx)[0], &A10)); 2644 else { 2645 PetscBool flg; 2646 2647 PetscCall(PetscObjectTypeCompare((PetscObject)std::get<1>(*ctx)[0], MATHERMITIANTRANSPOSEVIRTUAL, &flg)); 2648 if (flg) PetscCall(MatHermitianTransposeGetMat(std::get<1>(*ctx)[0], &A10)); 2649 } 2650 if (!A10) PetscCall(MatCreateSubMatrices(std::get<1>(*ctx)[0], 1, &sorted, &data->is, MAT_INITIAL_MATRIX, &sub)); 2651 else { 2652 if (flg) PetscCall(MatCreateTranspose(b[2], b + 1)); 2653 else PetscCall(MatCreateHermitianTranspose(b[2], b + 1)); 2654 } 2655 } else { 2656 PetscCall(MatSetOption(A01, MAT_SUBMAT_SINGLEIS, PETSC_TRUE)); 2657 PetscCall(MatCreateSubMatrices(A01, 1, &sorted, &data->is, MAT_INITIAL_MATRIX, &sub)); 2658 if (flg) PetscCall(MatCreateTranspose(*sub, b + 2)); 2659 else PetscCall(MatCreateHermitianTranspose(*sub, b + 2)); 2660 } 2661 if (A01 || !A10) { 2662 b[1] = sub[0]; 2663 PetscCall(PetscObjectReference((PetscObject)sub[0])); 2664 } 2665 PetscCall(MatDestroySubMatrices(1, &sub)); 2666 PetscCall(ISDestroy(&sorted)); 2667 b[3] = data->aux; 2668 PetscCall(MatCreateSchurComplement(subA[0], subA[1], b[1], b[2], b[3], &S)); 2669 PetscCall(MatSchurComplementSetKSP(S, ksp[0])); 2670 if (data->N != 1) { 2671 PetscCall(PCASMSetType(data->levels[0]->pc, PC_ASM_NONE)); /* "Neumann--Neumann" preconditioning with overlap and a Boolean partition of unity */ 2672 PetscCall(PCASMSetLocalSubdomains(data->levels[0]->pc, 1, &data->is, &loc)); 2673 PetscCall(PCSetFromOptions(data->levels[0]->pc)); /* action of eq. (15) of https://hal.science/hal-02343808v6/document (with a sign flip) */ 2674 s = data->levels[0]->pc; 2675 } else { 2676 is[0] = data->is; 2677 PetscCall(PetscObjectReference((PetscObject)is[0])); 2678 PetscCall(PetscObjectReference((PetscObject)b[3])); 2679 PetscCall(PCSetType(pc, PCASM)); /* change the type of the current PC */ 2680 data = nullptr; /* destroyed in the previous PCSetType(), so reset to NULL to avoid any faulty use */ 2681 PetscCall(PCAppendOptionsPrefix(pc, "pc_hpddm_coarse_")); /* same prefix as when using PCHPDDM with a single level */ 2682 PetscCall(PCASMSetLocalSubdomains(pc, 1, is, &loc)); 2683 PetscCall(ISDestroy(is)); 2684 PetscCall(ISDestroy(&loc)); 2685 s = pc; 2686 } 2687 PetscCall(PCHPDDMCommunicationAvoidingPCASM_Private(s, S, PETSC_TRUE)); /* the subdomain Mat is already known and the input IS of PCASMSetLocalSubdomains() is already sorted */ 2688 PetscTryMethod(s, "PCASMGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (s, &n, nullptr, &ksp)); 2689 PetscCheck(n == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of subdomain solver %" PetscInt_FMT " != 1", n); 2690 PetscCall(KSPGetPC(ksp[0], &inner)); 2691 PetscCall(PCSetType(inner, PCSHELL)); /* compute the action of the inverse of the local Schur complement with a PCSHELL */ 2692 b[0] = subA[0]; 2693 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 */ 2694 if (!data) PetscCall(PetscObjectDereference((PetscObject)b[3])); 2695 PetscCall(PetscObjectDereference((PetscObject)b[1])); 2696 PetscCall(PetscObjectDereference((PetscObject)b[2])); 2697 PetscCall(PCCreate(PETSC_COMM_SELF, &s)); 2698 PetscCall(PCSetOptionsPrefix(s, ((PetscObject)inner)->prefix)); 2699 PetscCall(PCSetOptionsPrefix(inner, nullptr)); 2700 PetscCall(KSPSetSkipPCSetFromOptions(ksp[0], PETSC_TRUE)); 2701 PetscCall(PCSetType(s, PCLU)); 2702 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 */ 2703 PetscCall(PCSetFromOptions(s)); 2704 PetscCall(PCFactorGetMatSolverType(s, &type)); 2705 PetscCall(PetscStrcmp(type, MATSOLVERMUMPS, &flg)); 2706 PetscCall(MatGetLocalSize(A11, &n, nullptr)); 2707 if (flg || n == 0) { 2708 PetscCall(PCSetOperators(s, N, N)); 2709 if (n) { 2710 PetscCall(PCFactorGetMatrix(s, b)); 2711 PetscCall(MatSetOptionsPrefix(*b, ((PetscObject)s)->prefix)); 2712 n = -1; 2713 PetscCall(PetscOptionsGetInt(nullptr, ((PetscObject)s)->prefix, "-mat_mumps_icntl_26", &n, nullptr)); 2714 if (n == 1) { /* allocates a square MatDense of size is[1]->map->n, so one */ 2715 PetscCall(MatNestGetISs(N, is, nullptr)); /* needs to be able to deactivate this path when dealing */ 2716 PetscCall(MatFactorSetSchurIS(*b, is[1])); /* with a large constraint space in order to avoid OOM */ 2717 } 2718 } else PetscCall(PCSetType(s, PCNONE)); /* empty local Schur complement (e.g., centralized on another process) */ 2719 } else { 2720 PetscCall(MatConvert(N, MATAIJ, MAT_INITIAL_MATRIX, b)); 2721 PetscCall(PCSetOperators(s, N, *b)); 2722 PetscCall(PetscObjectDereference((PetscObject)*b)); 2723 PetscCall(PetscObjectTypeCompareAny((PetscObject)s, &flg, PCLU, PCCHOLESKY, PCILU, PCICC, PCQR, "")); 2724 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 */ 2725 } 2726 PetscCall(PetscNew(&p)); 2727 p->first = s; 2728 if (n != 0) PetscCall(MatCreateVecs(*b, p->second, p->second + 1)); 2729 else p->second[0] = p->second[1] = nullptr; 2730 PetscCall(PCShellSetContext(inner, p)); 2731 PetscCall(PCShellSetApply(inner, PCApply_Nest)); 2732 PetscCall(PCShellSetView(inner, PCView_Nest)); 2733 PetscCall(PCShellSetDestroy(inner, PCDestroy_Nest)); 2734 PetscCall(PetscObjectDereference((PetscObject)N)); 2735 if (!data) { 2736 PetscCall(MatDestroy(&S)); 2737 PetscCall(ISDestroy(&unsorted)); 2738 PetscCall(MatDestroy(&C)); 2739 PetscCall(ISDestroy(&uis)); 2740 PetscCall(PetscFree(ctx)); 2741 #if PetscDefined(USE_DEBUG) 2742 PetscCall(ISDestroy(&dis)); 2743 PetscCall(MatDestroy(&daux)); 2744 #endif 2745 PetscFunctionReturn(PETSC_SUCCESS); 2746 } 2747 } 2748 if (!data->levels[0]->scatter) { 2749 PetscCall(MatCreateVecs(P, &xin, nullptr)); 2750 if (ismatis) PetscCall(MatDestroy(&P)); 2751 PetscCall(VecScatterCreate(xin, data->is, data->levels[0]->D, nullptr, &data->levels[0]->scatter)); 2752 PetscCall(VecDestroy(&xin)); 2753 } 2754 if (data->levels[0]->P) { 2755 /* if the pattern is the same and PCSetUp() has previously succeeded, reuse HPDDM buffers and connectivity */ 2756 PetscCall(HPDDM::Schwarz<PetscScalar>::destroy(data->levels[0], !pc->setupcalled || pc->flag == DIFFERENT_NONZERO_PATTERN ? PETSC_TRUE : PETSC_FALSE)); 2757 } 2758 if (!data->levels[0]->P) data->levels[0]->P = new HPDDM::Schwarz<PetscScalar>(); 2759 if (data->log_separate) PetscCall(PetscLogEventBegin(PC_HPDDM_SetUp[0], data->levels[0]->ksp, nullptr, nullptr, nullptr)); 2760 else PetscCall(PetscLogEventBegin(PC_HPDDM_Strc, data->levels[0]->ksp, nullptr, nullptr, nullptr)); 2761 /* HPDDM internal data structure */ 2762 PetscCall(data->levels[0]->P->structure(loc, data->is, !ctx ? sub[0] : nullptr, ismatis ? C : data->aux, data->levels)); 2763 if (!data->log_separate) PetscCall(PetscLogEventEnd(PC_HPDDM_Strc, data->levels[0]->ksp, nullptr, nullptr, nullptr)); 2764 /* matrix pencil of the generalized eigenvalue problem on the overlap (GenEO) */ 2765 if (!ctx) { 2766 if (data->deflation || overlap != -1) weighted = data->aux; 2767 else if (!data->B) { 2768 PetscBool cmp; 2769 2770 PetscCall(MatDuplicate(sub[0], MAT_COPY_VALUES, &weighted)); 2771 PetscCall(PetscObjectTypeCompareAny((PetscObject)weighted, &cmp, MATNORMAL, MATNORMALHERMITIAN, "")); 2772 if (cmp) flg = PETSC_FALSE; 2773 PetscCall(MatDiagonalScale(weighted, data->levels[0]->D, data->levels[0]->D)); 2774 /* neither MatDuplicate() nor MatDiagonalScale() handles the symmetry options, so propagate the options explicitly */ 2775 /* only useful for -mat_type baij -pc_hpddm_levels_1_st_pc_type cholesky (no problem with MATAIJ or MATSBAIJ) */ 2776 PetscCall(MatPropagateSymmetryOptions(sub[0], weighted)); 2777 if (PetscDefined(USE_DEBUG) && PetscBool3ToBool(data->Neumann)) { 2778 Mat *sub, A[3]; 2779 PetscReal norm[2]; 2780 PetscBool flg; 2781 2782 PetscCall(PetscObjectTypeCompare((PetscObject)P, MATMPISBAIJ, &flg)); /* MatCreateSubMatrices() does not work with MATSBAIJ and unsorted ISes, so convert to MPIAIJ */ 2783 if (flg) PetscCall(MatConvert(P, MATMPIAIJ, MAT_INITIAL_MATRIX, A)); 2784 else { 2785 A[0] = P; 2786 PetscCall(PetscObjectReference((PetscObject)P)); 2787 } 2788 PetscCall(MatCreateSubMatrices(A[0], 1, &data->is, &data->is, MAT_INITIAL_MATRIX, &sub)); 2789 PetscCall(MatDiagonalScale(sub[0], data->levels[0]->D, data->levels[0]->D)); 2790 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 */ 2791 PetscCall(MatConvert(weighted, MATSEQAIJ, MAT_INITIAL_MATRIX, A + 2)); 2792 PetscCall(MatAXPY(A[1], -1.0, A[2], UNKNOWN_NONZERO_PATTERN)); 2793 PetscCall(MatNorm(A[1], NORM_FROBENIUS, norm)); 2794 if (norm[0]) { 2795 PetscCall(MatNorm(A[2], NORM_FROBENIUS, norm + 1)); 2796 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 : ""); 2797 } 2798 PetscCall(MatDestroySubMatrices(1, &sub)); 2799 for (PetscInt i = 0; i < 3; ++i) PetscCall(MatDestroy(A + i)); 2800 } 2801 } else weighted = data->B; 2802 } else weighted = nullptr; 2803 /* SLEPc is used inside the loaded symbol */ 2804 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)); 2805 if (!ctx && data->share && overlap == -1) { 2806 Mat st[2]; 2807 2808 PetscCall(KSPGetOperators(ksp[0], st, st + 1)); 2809 PetscCall(MatCopy(subA[0], st[0], structure)); 2810 if (subA[1] != subA[0] || st[1] != st[0]) PetscCall(MatCopy(subA[1], st[1], SAME_NONZERO_PATTERN)); 2811 PetscCall(PetscObjectDereference((PetscObject)subA[0])); 2812 } 2813 if (data->log_separate) PetscCall(PetscLogEventEnd(PC_HPDDM_SetUp[0], data->levels[0]->ksp, nullptr, nullptr, nullptr)); 2814 if (ismatis) PetscCall(MatISGetLocalMat(C, &N)); 2815 else N = data->aux; 2816 if (!ctx) P = sub[0]; 2817 else P = S; 2818 /* going through the grid hierarchy */ 2819 for (n = 1; n < data->N; ++n) { 2820 if (data->log_separate) PetscCall(PetscLogEventBegin(PC_HPDDM_SetUp[n], data->levels[n]->ksp, nullptr, nullptr, nullptr)); 2821 /* method composed in the loaded symbol since there, SLEPc is used as well */ 2822 PetscTryMethod(data->levels[0]->ksp, "PCHPDDMSetUp_Private_C", (Mat *, Mat *, PetscInt, PetscInt *const, PC_HPDDM_Level **const), (&P, &N, n, &data->N, data->levels)); 2823 if (data->log_separate) PetscCall(PetscLogEventEnd(PC_HPDDM_SetUp[n], data->levels[n]->ksp, nullptr, nullptr, nullptr)); 2824 } 2825 /* reset to NULL to avoid any faulty use */ 2826 PetscCall(PetscObjectComposeFunction((PetscObject)data->levels[0]->ksp, "PCHPDDMSetUp_Private_C", nullptr)); 2827 if (!ismatis) PetscCall(PetscObjectComposeFunction((PetscObject)pc->pmat, "PCHPDDMAlgebraicAuxiliaryMat_C", nullptr)); 2828 else PetscCall(PetscObjectDereference((PetscObject)C)); /* matching PetscObjectReference() above */ 2829 for (n = 0; n < data->N - 1; ++n) 2830 if (data->levels[n]->P) { 2831 /* HPDDM internal work buffers */ 2832 data->levels[n]->P->setBuffer(); 2833 data->levels[n]->P->super::start(); 2834 } 2835 if (ismatis || !subdomains) PetscCall(PCHPDDMDestroySubMatrices_Private(PetscBool3ToBool(data->Neumann), PetscBool(algebraic && !block && overlap == -1), sub)); 2836 if (ismatis) data->is = nullptr; 2837 for (n = 0; n < data->N - 1 + (reused > 0); ++n) { 2838 if (data->levels[n]->P) { 2839 PC spc; 2840 2841 /* force the PC to be PCSHELL to do the coarse grid corrections */ 2842 PetscCall(KSPSetSkipPCSetFromOptions(data->levels[n]->ksp, PETSC_TRUE)); 2843 PetscCall(KSPGetPC(data->levels[n]->ksp, &spc)); 2844 PetscCall(PCSetType(spc, PCSHELL)); 2845 PetscCall(PCShellSetContext(spc, data->levels[n])); 2846 PetscCall(PCShellSetSetUp(spc, PCSetUp_HPDDMShell)); 2847 PetscCall(PCShellSetApply(spc, PCApply_HPDDMShell)); 2848 PetscCall(PCShellSetMatApply(spc, PCMatApply_HPDDMShell)); 2849 PetscCall(PCShellSetApplyTranspose(spc, PCApplyTranspose_HPDDMShell)); 2850 PetscCall(PCShellSetMatApplyTranspose(spc, PCMatApplyTranspose_HPDDMShell)); 2851 if (ctx && n == 0) { 2852 Mat Amat, Pmat; 2853 PetscInt m, M; 2854 std::tuple<Mat, PetscSF, Vec[2]> *ctx; 2855 2856 PetscCall(KSPGetOperators(data->levels[n]->ksp, nullptr, &Pmat)); 2857 PetscCall(MatGetLocalSize(Pmat, &m, nullptr)); 2858 PetscCall(MatGetSize(Pmat, &M, nullptr)); 2859 PetscCall(PetscNew(&ctx)); 2860 std::get<0>(*ctx) = S; 2861 std::get<1>(*ctx) = data->levels[n]->scatter; 2862 PetscCall(MatCreateShell(PetscObjectComm((PetscObject)data->levels[n]->ksp), m, m, M, M, ctx, &Amat)); 2863 PetscCall(MatShellSetOperation(Amat, MATOP_MULT, (PetscErrorCodeFn *)MatMult_Schur<false>)); 2864 PetscCall(MatShellSetOperation(Amat, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)MatMult_Schur<true>)); 2865 PetscCall(MatShellSetOperation(Amat, MATOP_DESTROY, (PetscErrorCodeFn *)MatDestroy_Schur)); 2866 PetscCall(MatCreateVecs(S, std::get<2>(*ctx), std::get<2>(*ctx) + 1)); 2867 PetscCall(KSPSetOperators(data->levels[n]->ksp, Amat, Pmat)); 2868 PetscCall(PetscObjectDereference((PetscObject)Amat)); 2869 } 2870 PetscCall(PCShellSetDestroy(spc, PCDestroy_HPDDMShell)); 2871 if (!data->levels[n]->pc) PetscCall(PCCreate(PetscObjectComm((PetscObject)data->levels[n]->ksp), &data->levels[n]->pc)); 2872 if (n < reused) { 2873 PetscCall(PCSetReusePreconditioner(spc, PETSC_TRUE)); 2874 PetscCall(PCSetReusePreconditioner(data->levels[n]->pc, PETSC_TRUE)); 2875 } 2876 PetscCall(PCSetUp(spc)); 2877 } 2878 } 2879 if (ctx) PetscCall(MatDestroy(&S)); 2880 if (overlap == -1) PetscCall(PetscObjectComposeFunction((PetscObject)pc->pmat, "PCHPDDMAlgebraicAuxiliaryMat_Private_C", nullptr)); 2881 } else flg = reused ? PETSC_FALSE : PETSC_TRUE; 2882 if (!ismatis && subdomains) { 2883 if (flg) PetscCall(KSPGetPC(data->levels[0]->ksp, &inner)); 2884 else inner = data->levels[0]->pc; 2885 if (inner) { 2886 if (!inner->setupcalled) PetscCall(PCSetType(inner, PCASM)); 2887 PetscCall(PCSetFromOptions(inner)); 2888 PetscCall(PetscStrcmp(((PetscObject)inner)->type_name, PCASM, &flg)); 2889 if (flg) { 2890 if (!inner->setupcalled) { /* evaluates to PETSC_FALSE when -pc_hpddm_block_splitting */ 2891 IS sorted; /* PCASM will sort the input IS, duplicate it to return an unmodified (PCHPDDM) input IS */ 2892 2893 PetscCall(ISDuplicate(is[0], &sorted)); 2894 PetscCall(PCASMSetLocalSubdomains(inner, 1, &sorted, &loc)); 2895 PetscCall(PetscObjectDereference((PetscObject)sorted)); 2896 } 2897 if (!PetscBool3ToBool(data->Neumann) && data->N > 1) { /* subdomain matrices are already created for the eigenproblem, reuse them for the fine-level PC */ 2898 PetscCall(PCHPDDMPermute_Private(*is, nullptr, nullptr, sub[0], &P, nullptr)); 2899 PetscCall(PCHPDDMCommunicationAvoidingPCASM_Private(inner, P, algebraic)); 2900 PetscCall(PetscObjectDereference((PetscObject)P)); 2901 } 2902 } 2903 } 2904 if (data->N > 1) { 2905 if (overlap != 1) PetscCall(PCHPDDMDestroySubMatrices_Private(PetscBool3ToBool(data->Neumann), PetscBool(algebraic && !block && overlap == -1), sub)); 2906 if (overlap == 1) PetscCall(MatDestroy(subA)); 2907 } 2908 } 2909 PetscCall(ISDestroy(&loc)); 2910 } else data->N = 1 + reused; /* enforce this value to 1 + reused if there is no way to build another level */ 2911 if (requested != data->N + reused) { 2912 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, 2913 data->N, pcpre ? pcpre : "")); 2914 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 : "", 2915 data->N, pcpre ? pcpre : "", data->N)); 2916 /* cannot use PCDestroy_HPDDMShell() because PCSHELL not set for unassembled levels */ 2917 for (n = data->N - 1; n < requested - 1; ++n) { 2918 if (data->levels[n]->P) { 2919 PetscCall(HPDDM::Schwarz<PetscScalar>::destroy(data->levels[n], PETSC_TRUE)); 2920 PetscCall(VecDestroyVecs(1, &data->levels[n]->v[0])); 2921 PetscCall(VecDestroyVecs(2, &data->levels[n]->v[1])); 2922 PetscCall(MatDestroy(data->levels[n]->V)); 2923 PetscCall(MatDestroy(data->levels[n]->V + 1)); 2924 PetscCall(MatDestroy(data->levels[n]->V + 2)); 2925 PetscCall(VecDestroy(&data->levels[n]->D)); 2926 PetscCall(PetscSFDestroy(&data->levels[n]->scatter)); 2927 } 2928 } 2929 if (reused) { 2930 for (n = reused; n < PETSC_PCHPDDM_MAXLEVELS && data->levels[n]; ++n) { 2931 PetscCall(KSPDestroy(&data->levels[n]->ksp)); 2932 PetscCall(PCDestroy(&data->levels[n]->pc)); 2933 } 2934 } 2935 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, 2936 data->N, reused, data->N, pcpre ? pcpre : "", pcpre ? pcpre : "", data->N, pcpre ? pcpre : "", data->N); 2937 } 2938 /* these solvers are created after PCSetFromOptions() is called */ 2939 if (pc->setfromoptionscalled) { 2940 for (n = 0; n < data->N; ++n) { 2941 if (data->levels[n]->ksp) PetscCall(KSPSetFromOptions(data->levels[n]->ksp)); 2942 if (data->levels[n]->pc) PetscCall(PCSetFromOptions(data->levels[n]->pc)); 2943 } 2944 pc->setfromoptionscalled = 0; 2945 } 2946 data->N += reused; 2947 if (data->share && swap) { 2948 /* swap back pointers so that variables follow the user-provided numbering */ 2949 std::swap(C, data->aux); 2950 std::swap(uis, data->is); 2951 PetscCall(MatDestroy(&C)); 2952 PetscCall(ISDestroy(&uis)); 2953 } 2954 if (algebraic) PetscCall(MatDestroy(&data->aux)); 2955 if (unsorted && unsorted != is[0]) { 2956 PetscCall(ISCopy(unsorted, data->is)); 2957 PetscCall(ISDestroy(&unsorted)); 2958 } 2959 #if PetscDefined(USE_DEBUG) 2960 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); 2961 if (data->is) { 2962 PetscCall(ISEqualUnsorted(data->is, dis, &flg)); 2963 PetscCall(ISDestroy(&dis)); 2964 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Input IS and output IS are not equal"); 2965 } 2966 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); 2967 if (data->aux) { 2968 PetscCall(MatMultEqual(data->aux, daux, 20, &flg)); 2969 PetscCall(MatDestroy(&daux)); 2970 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Input Mat and output Mat are not equal"); 2971 } 2972 #endif 2973 PetscFunctionReturn(PETSC_SUCCESS); 2974 } 2975 2976 /*@ 2977 PCHPDDMSetCoarseCorrectionType - Sets the coarse correction type. 2978 2979 Collective 2980 2981 Input Parameters: 2982 + pc - preconditioner context 2983 - type - `PC_HPDDM_COARSE_CORRECTION_DEFLATED`, `PC_HPDDM_COARSE_CORRECTION_ADDITIVE`, `PC_HPDDM_COARSE_CORRECTION_BALANCED`, or `PC_HPDDM_COARSE_CORRECTION_NONE` 2984 2985 Options Database Key: 2986 . -pc_hpddm_coarse_correction <deflated, additive, balanced, none> - type of coarse correction to apply 2987 2988 Level: intermediate 2989 2990 .seealso: [](ch_ksp), `PCHPDDMGetCoarseCorrectionType()`, `PCHPDDM`, `PCHPDDMCoarseCorrectionType` 2991 @*/ 2992 PetscErrorCode PCHPDDMSetCoarseCorrectionType(PC pc, PCHPDDMCoarseCorrectionType type) 2993 { 2994 PetscFunctionBegin; 2995 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2996 PetscValidLogicalCollectiveEnum(pc, type, 2); 2997 PetscTryMethod(pc, "PCHPDDMSetCoarseCorrectionType_C", (PC, PCHPDDMCoarseCorrectionType), (pc, type)); 2998 PetscFunctionReturn(PETSC_SUCCESS); 2999 } 3000 3001 /*@ 3002 PCHPDDMGetCoarseCorrectionType - Gets the coarse correction type. 3003 3004 Input Parameter: 3005 . pc - preconditioner context 3006 3007 Output Parameter: 3008 . type - `PC_HPDDM_COARSE_CORRECTION_DEFLATED`, `PC_HPDDM_COARSE_CORRECTION_ADDITIVE`, `PC_HPDDM_COARSE_CORRECTION_BALANCED`, or `PC_HPDDM_COARSE_CORRECTION_NONE` 3009 3010 Level: intermediate 3011 3012 .seealso: [](ch_ksp), `PCHPDDMSetCoarseCorrectionType()`, `PCHPDDM`, `PCHPDDMCoarseCorrectionType` 3013 @*/ 3014 PetscErrorCode PCHPDDMGetCoarseCorrectionType(PC pc, PCHPDDMCoarseCorrectionType *type) 3015 { 3016 PetscFunctionBegin; 3017 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 3018 if (type) { 3019 PetscAssertPointer(type, 2); 3020 PetscUseMethod(pc, "PCHPDDMGetCoarseCorrectionType_C", (PC, PCHPDDMCoarseCorrectionType *), (pc, type)); 3021 } 3022 PetscFunctionReturn(PETSC_SUCCESS); 3023 } 3024 3025 static PetscErrorCode PCHPDDMSetCoarseCorrectionType_HPDDM(PC pc, PCHPDDMCoarseCorrectionType type) 3026 { 3027 PC_HPDDM *data = (PC_HPDDM *)pc->data; 3028 3029 PetscFunctionBegin; 3030 data->correction = type; 3031 PetscFunctionReturn(PETSC_SUCCESS); 3032 } 3033 3034 static PetscErrorCode PCHPDDMGetCoarseCorrectionType_HPDDM(PC pc, PCHPDDMCoarseCorrectionType *type) 3035 { 3036 PC_HPDDM *data = (PC_HPDDM *)pc->data; 3037 3038 PetscFunctionBegin; 3039 *type = data->correction; 3040 PetscFunctionReturn(PETSC_SUCCESS); 3041 } 3042 3043 /*@ 3044 PCHPDDMSetSTShareSubKSP - Sets whether the `KSP` in SLEPc `ST` and the fine-level subdomain solver should be shared. 3045 3046 Input Parameters: 3047 + pc - preconditioner context 3048 - share - whether the `KSP` should be shared or not 3049 3050 Note: 3051 This is not the same as `PCSetReusePreconditioner()`. Given certain conditions (visible using -info), a symbolic factorization can be skipped 3052 when using a subdomain `PCType` such as `PCLU` or `PCCHOLESKY`. 3053 3054 Level: advanced 3055 3056 .seealso: [](ch_ksp), `PCHPDDM`, `PCHPDDMGetSTShareSubKSP()` 3057 @*/ 3058 PetscErrorCode PCHPDDMSetSTShareSubKSP(PC pc, PetscBool share) 3059 { 3060 PetscFunctionBegin; 3061 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 3062 PetscTryMethod(pc, "PCHPDDMSetSTShareSubKSP_C", (PC, PetscBool), (pc, share)); 3063 PetscFunctionReturn(PETSC_SUCCESS); 3064 } 3065 3066 /*@ 3067 PCHPDDMGetSTShareSubKSP - Gets whether the `KSP` in SLEPc `ST` and the fine-level subdomain solver is shared. 3068 3069 Input Parameter: 3070 . pc - preconditioner context 3071 3072 Output Parameter: 3073 . share - whether the `KSP` is shared or not 3074 3075 Note: 3076 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 3077 when using a subdomain `PCType` such as `PCLU` or `PCCHOLESKY`. 3078 3079 Level: advanced 3080 3081 .seealso: [](ch_ksp), `PCHPDDM`, `PCHPDDMSetSTShareSubKSP()` 3082 @*/ 3083 PetscErrorCode PCHPDDMGetSTShareSubKSP(PC pc, PetscBool *share) 3084 { 3085 PetscFunctionBegin; 3086 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 3087 if (share) { 3088 PetscAssertPointer(share, 2); 3089 PetscUseMethod(pc, "PCHPDDMGetSTShareSubKSP_C", (PC, PetscBool *), (pc, share)); 3090 } 3091 PetscFunctionReturn(PETSC_SUCCESS); 3092 } 3093 3094 static PetscErrorCode PCHPDDMSetSTShareSubKSP_HPDDM(PC pc, PetscBool share) 3095 { 3096 PC_HPDDM *data = (PC_HPDDM *)pc->data; 3097 3098 PetscFunctionBegin; 3099 data->share = share; 3100 PetscFunctionReturn(PETSC_SUCCESS); 3101 } 3102 3103 static PetscErrorCode PCHPDDMGetSTShareSubKSP_HPDDM(PC pc, PetscBool *share) 3104 { 3105 PC_HPDDM *data = (PC_HPDDM *)pc->data; 3106 3107 PetscFunctionBegin; 3108 *share = data->share; 3109 PetscFunctionReturn(PETSC_SUCCESS); 3110 } 3111 3112 /*@ 3113 PCHPDDMSetDeflationMat - Sets the deflation space used to assemble a coarser operator. 3114 3115 Input Parameters: 3116 + pc - preconditioner context 3117 . is - index set of the local deflation matrix 3118 - U - deflation sequential matrix stored as a `MATSEQDENSE` 3119 3120 Level: advanced 3121 3122 .seealso: [](ch_ksp), `PCHPDDM`, `PCDeflationSetSpace()`, `PCMGSetRestriction()` 3123 @*/ 3124 PetscErrorCode PCHPDDMSetDeflationMat(PC pc, IS is, Mat U) 3125 { 3126 PetscFunctionBegin; 3127 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 3128 PetscValidHeaderSpecific(is, IS_CLASSID, 2); 3129 PetscValidHeaderSpecific(U, MAT_CLASSID, 3); 3130 PetscTryMethod(pc, "PCHPDDMSetDeflationMat_C", (PC, IS, Mat), (pc, is, U)); 3131 PetscFunctionReturn(PETSC_SUCCESS); 3132 } 3133 3134 static PetscErrorCode PCHPDDMSetDeflationMat_HPDDM(PC pc, IS is, Mat U) 3135 { 3136 PetscFunctionBegin; 3137 PetscCall(PCHPDDMSetAuxiliaryMat_Private(pc, is, U, PETSC_TRUE)); 3138 PetscFunctionReturn(PETSC_SUCCESS); 3139 } 3140 3141 PetscErrorCode HPDDMLoadDL_Private(PetscBool *found) 3142 { 3143 PetscBool flg; 3144 char lib[PETSC_MAX_PATH_LEN], dlib[PETSC_MAX_PATH_LEN], dir[PETSC_MAX_PATH_LEN]; 3145 3146 PetscFunctionBegin; 3147 PetscAssertPointer(found, 1); 3148 PetscCall(PetscStrncpy(dir, "${PETSC_LIB_DIR}", sizeof(dir))); 3149 PetscCall(PetscOptionsGetString(nullptr, nullptr, "-hpddm_dir", dir, sizeof(dir), nullptr)); 3150 PetscCall(PetscSNPrintf(lib, sizeof(lib), "%s/libhpddm_petsc", dir)); 3151 PetscCall(PetscDLLibraryRetrieve(PETSC_COMM_SELF, lib, dlib, 1024, found)); 3152 #if defined(SLEPC_LIB_DIR) /* this variable is passed during SLEPc ./configure when PETSc has not been configured */ 3153 if (!*found) { /* with --download-hpddm since slepcconf.h is not yet built (and thus can't be included) */ 3154 PetscCall(PetscStrncpy(dir, HPDDM_STR(SLEPC_LIB_DIR), sizeof(dir))); 3155 PetscCall(PetscSNPrintf(lib, sizeof(lib), "%s/libhpddm_petsc", dir)); 3156 PetscCall(PetscDLLibraryRetrieve(PETSC_COMM_SELF, lib, dlib, 1024, found)); 3157 } 3158 #endif 3159 if (!*found) { /* probable options for this to evaluate to PETSC_TRUE: system inconsistency (libhpddm_petsc moved by user?) or PETSc configured without --download-slepc */ 3160 PetscCall(PetscOptionsGetenv(PETSC_COMM_SELF, "SLEPC_DIR", dir, sizeof(dir), &flg)); 3161 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 */ 3162 PetscCall(PetscSNPrintf(lib, sizeof(lib), "%s/lib/libslepc", dir)); 3163 PetscCall(PetscDLLibraryRetrieve(PETSC_COMM_SELF, lib, dlib, 1024, found)); 3164 PetscCheck(*found, PETSC_COMM_SELF, PETSC_ERR_PLIB, "%s not found but SLEPC_DIR=%s", lib, dir); 3165 PetscCall(PetscDLLibraryAppend(PETSC_COMM_SELF, &PetscDLLibrariesLoaded, dlib)); 3166 PetscCall(PetscSNPrintf(lib, sizeof(lib), "%s/lib/libhpddm_petsc", dir)); /* libhpddm_petsc is always in the same directory as libslepc */ 3167 PetscCall(PetscDLLibraryRetrieve(PETSC_COMM_SELF, lib, dlib, 1024, found)); 3168 } 3169 } 3170 PetscCheck(*found, PETSC_COMM_SELF, PETSC_ERR_PLIB, "%s not found", lib); 3171 PetscCall(PetscDLLibraryAppend(PETSC_COMM_SELF, &PetscDLLibrariesLoaded, dlib)); 3172 PetscFunctionReturn(PETSC_SUCCESS); 3173 } 3174 3175 /*MC 3176 PCHPDDM - Interface with the HPDDM library. 3177 3178 This `PC` may be used to build multilevel spectral domain decomposition methods based on the GenEO framework {cite}`spillane2011robust` {cite}`al2021multilevel`. 3179 It may be viewed as an alternative to spectral 3180 AMGe or `PCBDDC` with adaptive selection of constraints. The interface is explained in details in {cite}`jolivetromanzampini2020` 3181 3182 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`). 3183 3184 For multilevel preconditioning, when using an assembled or hierarchical Pmat, one must provide an auxiliary local `Mat` (unassembled local operator for GenEO) using 3185 `PCHPDDMSetAuxiliaryMat()`. Calling this routine is not needed when using a `MATIS` Pmat, assembly is done internally using `MatConvert()`. 3186 3187 Options Database Keys: 3188 + -pc_hpddm_define_subdomains <true, default=false> - on the finest level, calls `PCASMSetLocalSubdomains()` with the `IS` supplied in `PCHPDDMSetAuxiliaryMat()` 3189 (not relevant with an unassembled Pmat) 3190 . -pc_hpddm_has_neumann <true, default=false> - on the finest level, informs the `PC` that the local Neumann matrix is supplied in `PCHPDDMSetAuxiliaryMat()` 3191 - -pc_hpddm_coarse_correction <type, default=deflated> - determines the `PCHPDDMCoarseCorrectionType` when calling `PCApply()` 3192 3193 Options for subdomain solvers, subdomain eigensolvers (for computing deflation vectors), and the coarse solver can be set using the following options database prefixes. 3194 .vb 3195 -pc_hpddm_levels_%d_pc_ 3196 -pc_hpddm_levels_%d_ksp_ 3197 -pc_hpddm_levels_%d_eps_ 3198 -pc_hpddm_levels_%d_p 3199 -pc_hpddm_levels_%d_mat_type 3200 -pc_hpddm_coarse_ 3201 -pc_hpddm_coarse_p 3202 -pc_hpddm_coarse_mat_type 3203 -pc_hpddm_coarse_mat_filter 3204 .ve 3205 3206 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 3207 -pc_hpddm_coarse_p 2 -pc_hpddm_coarse_mat_type baij will use 10 deflation vectors per subdomain on the fine "level 1", 3208 aggregate the fine subdomains into 4 "level 2" subdomains, then use 10 deflation vectors per subdomain on "level 2", 3209 and assemble the coarse matrix (of dimension 4 x 10 = 40) on two processes as a `MATBAIJ` (default is `MATSBAIJ`). 3210 3211 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. 3212 3213 Level: intermediate 3214 3215 Notes: 3216 This preconditioner requires that PETSc is built with SLEPc (``--download-slepc``). 3217 3218 By default, the underlying concurrent eigenproblems 3219 are solved using SLEPc shift-and-invert spectral transformation. This is usually what gives the best performance for GenEO, cf. 3220 {cite}`spillane2011robust` {cite}`jolivet2013scalabledd`. As 3221 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 3222 -pc_hpddm_levels_1_st_type sinvert. There are furthermore three options related to the (subdomain-wise local) eigensolver that are not described in 3223 SLEPc documentation since they are specific to `PCHPDDM`. 3224 .vb 3225 -pc_hpddm_levels_1_st_share_sub_ksp 3226 -pc_hpddm_levels_%d_eps_threshold_absolute 3227 -pc_hpddm_levels_1_eps_use_inertia 3228 .ve 3229 3230 The first option from the list only applies to the fine-level eigensolver, see `PCHPDDMSetSTShareSubKSP()`. The second option from the list is 3231 used to filter eigenmodes retrieved after convergence of `EPSSolve()` at "level N" such that eigenvectors used to define a "level N+1" coarse 3232 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 3233 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 3234 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 3235 to supply -pc_hpddm_levels_1_eps_nev. This last option also only applies to the fine-level (N = 1) eigensolver. 3236 3237 See also {cite}`dolean2015introduction`, {cite}`al2022robust`, {cite}`al2022robustpd`, and {cite}`nataf2022recent` 3238 3239 .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCHPDDMSetAuxiliaryMat()`, `MATIS`, `PCBDDC`, `PCDEFLATION`, `PCTELESCOPE`, `PCASM`, 3240 `PCHPDDMSetCoarseCorrectionType()`, `PCHPDDMHasNeumannMat()`, `PCHPDDMSetRHSMat()`, `PCHPDDMSetDeflationMat()`, `PCHPDDMSetSTShareSubKSP()`, 3241 `PCHPDDMGetSTShareSubKSP()`, `PCHPDDMGetCoarseCorrectionType()`, `PCHPDDMGetComplexities()` 3242 M*/ 3243 PETSC_EXTERN PetscErrorCode PCCreate_HPDDM(PC pc) 3244 { 3245 PC_HPDDM *data; 3246 PetscBool found; 3247 3248 PetscFunctionBegin; 3249 if (!loadedSym) { 3250 PetscCall(HPDDMLoadDL_Private(&found)); 3251 if (found) PetscCall(PetscDLLibrarySym(PETSC_COMM_SELF, &PetscDLLibrariesLoaded, nullptr, "PCHPDDM_Internal", (void **)&loadedSym)); 3252 } 3253 PetscCheck(loadedSym, PETSC_COMM_SELF, PETSC_ERR_PLIB, "PCHPDDM_Internal symbol not found in loaded libhpddm_petsc"); 3254 PetscCall(PetscNew(&data)); 3255 pc->data = data; 3256 data->Neumann = PETSC_BOOL3_UNKNOWN; 3257 pc->ops->reset = PCReset_HPDDM; 3258 pc->ops->destroy = PCDestroy_HPDDM; 3259 pc->ops->setfromoptions = PCSetFromOptions_HPDDM; 3260 pc->ops->setup = PCSetUp_HPDDM; 3261 pc->ops->apply = PCApply_HPDDM<false>; 3262 pc->ops->matapply = PCMatApply_HPDDM<false>; 3263 pc->ops->applytranspose = PCApply_HPDDM<true>; 3264 pc->ops->matapplytranspose = PCMatApply_HPDDM<true>; 3265 pc->ops->view = PCView_HPDDM; 3266 pc->ops->presolve = PCPreSolve_HPDDM; 3267 3268 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetAuxiliaryMat_C", PCHPDDMSetAuxiliaryMat_HPDDM)); 3269 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMHasNeumannMat_C", PCHPDDMHasNeumannMat_HPDDM)); 3270 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetRHSMat_C", PCHPDDMSetRHSMat_HPDDM)); 3271 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetCoarseCorrectionType_C", PCHPDDMSetCoarseCorrectionType_HPDDM)); 3272 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMGetCoarseCorrectionType_C", PCHPDDMGetCoarseCorrectionType_HPDDM)); 3273 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetSTShareSubKSP_C", PCHPDDMSetSTShareSubKSP_HPDDM)); 3274 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMGetSTShareSubKSP_C", PCHPDDMGetSTShareSubKSP_HPDDM)); 3275 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetDeflationMat_C", PCHPDDMSetDeflationMat_HPDDM)); 3276 PetscFunctionReturn(PETSC_SUCCESS); 3277 } 3278 3279 /*@C 3280 PCHPDDMInitializePackage - This function initializes everything in the `PCHPDDM` package. It is called from `PCInitializePackage()`. 3281 3282 Level: developer 3283 3284 .seealso: [](ch_ksp), `PetscInitialize()` 3285 @*/ 3286 PetscErrorCode PCHPDDMInitializePackage(void) 3287 { 3288 char ename[32]; 3289 3290 PetscFunctionBegin; 3291 if (PCHPDDMPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS); 3292 PCHPDDMPackageInitialized = PETSC_TRUE; 3293 PetscCall(PetscRegisterFinalize(PCHPDDMFinalizePackage)); 3294 /* general events registered once during package initialization */ 3295 /* some of these events are not triggered in libpetsc, */ 3296 /* but rather directly in libhpddm_petsc, */ 3297 /* which is in charge of performing the following operations */ 3298 3299 /* domain decomposition structure from Pmat sparsity pattern */ 3300 PetscCall(PetscLogEventRegister("PCHPDDMStrc", PC_CLASSID, &PC_HPDDM_Strc)); 3301 /* Galerkin product, redistribution, and setup (not triggered in libpetsc) */ 3302 PetscCall(PetscLogEventRegister("PCHPDDMPtAP", PC_CLASSID, &PC_HPDDM_PtAP)); 3303 /* Galerkin product with summation, redistribution, and setup (not triggered in libpetsc) */ 3304 PetscCall(PetscLogEventRegister("PCHPDDMPtBP", PC_CLASSID, &PC_HPDDM_PtBP)); 3305 /* next level construction using PtAP and PtBP (not triggered in libpetsc) */ 3306 PetscCall(PetscLogEventRegister("PCHPDDMNext", PC_CLASSID, &PC_HPDDM_Next)); 3307 static_assert(PETSC_PCHPDDM_MAXLEVELS <= 9, "PETSC_PCHPDDM_MAXLEVELS value is too high"); 3308 for (PetscInt i = 1; i < PETSC_PCHPDDM_MAXLEVELS; ++i) { 3309 PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCHPDDMSetUp L%1" PetscInt_FMT, i)); 3310 /* events during a PCSetUp() at level #i _except_ the assembly */ 3311 /* of the Galerkin operator of the coarser level #(i + 1) */ 3312 PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_HPDDM_SetUp[i - 1])); 3313 PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCHPDDMSolve L%1" PetscInt_FMT, i)); 3314 /* events during a PCApply() at level #i _except_ */ 3315 /* the KSPSolve() of the coarser level #(i + 1) */ 3316 PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_HPDDM_Solve[i - 1])); 3317 } 3318 PetscFunctionReturn(PETSC_SUCCESS); 3319 } 3320 3321 /*@C 3322 PCHPDDMFinalizePackage - This function frees everything from the `PCHPDDM` package. It is called from `PetscFinalize()`. 3323 3324 Level: developer 3325 3326 .seealso: [](ch_ksp), `PetscFinalize()` 3327 @*/ 3328 PetscErrorCode PCHPDDMFinalizePackage(void) 3329 { 3330 PetscFunctionBegin; 3331 PCHPDDMPackageInitialized = PETSC_FALSE; 3332 PetscFunctionReturn(PETSC_SUCCESS); 3333 } 3334 3335 static PetscErrorCode MatMult_Harmonic(Mat A, Vec x, Vec y) 3336 { 3337 Harmonic h; /* [ A_00 A_01 ], furthermore, A_00 = [ A_loc,loc A_loc,ovl ], thus, A_01 = [ ] */ 3338 /* [ A_10 A_11 A_12 ] [ A_ovl,loc A_ovl,ovl ] [ A_ovl,1 ] */ 3339 Vec sub; /* y = A x = R_loc R_0 [ A_00 A_01 ]^-1 R_loc = [ I_loc ] */ 3340 /* [ A_10 A_11 ] R_1^T A_12 x [ ] */ 3341 PetscFunctionBegin; 3342 PetscCall(MatShellGetContext(A, &h)); 3343 PetscCall(VecSet(h->v, 0.0)); 3344 PetscCall(VecGetSubVector(h->v, h->is[0], &sub)); 3345 PetscCall(MatMult(h->A[0], x, sub)); 3346 PetscCall(VecRestoreSubVector(h->v, h->is[0], &sub)); 3347 PetscCall(KSPSolve(h->ksp, h->v, h->v)); 3348 PetscCall(VecISCopy(h->v, h->is[1], SCATTER_REVERSE, y)); 3349 PetscFunctionReturn(PETSC_SUCCESS); 3350 } 3351 3352 static PetscErrorCode MatMultTranspose_Harmonic(Mat A, Vec y, Vec x) 3353 { 3354 Harmonic h; /* x = A^T y = [ A_00 A_01 ]^-T R_0^T R_loc^T y */ 3355 Vec sub; /* A_12^T R_1 [ A_10 A_11 ] */ 3356 3357 PetscFunctionBegin; 3358 PetscCall(MatShellGetContext(A, &h)); 3359 PetscCall(VecSet(h->v, 0.0)); 3360 PetscCall(VecISCopy(h->v, h->is[1], SCATTER_FORWARD, y)); 3361 PetscCall(KSPSolveTranspose(h->ksp, h->v, h->v)); 3362 PetscCall(VecGetSubVector(h->v, h->is[0], &sub)); 3363 PetscCall(MatMultTranspose(h->A[0], sub, x)); 3364 PetscCall(VecRestoreSubVector(h->v, h->is[0], &sub)); 3365 PetscFunctionReturn(PETSC_SUCCESS); 3366 } 3367 3368 static PetscErrorCode MatProduct_AB_Harmonic(Mat S, Mat X, Mat Y, void *) 3369 { 3370 Harmonic h; 3371 Mat A, B; 3372 Vec a, b; 3373 3374 PetscFunctionBegin; 3375 PetscCall(MatShellGetContext(S, &h)); 3376 PetscCall(MatMatMult(h->A[0], X, MAT_INITIAL_MATRIX, PETSC_CURRENT, &A)); 3377 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, h->ksp->pc->mat->rmap->n, A->cmap->n, nullptr, &B)); 3378 for (PetscInt i = 0; i < A->cmap->n; ++i) { 3379 PetscCall(MatDenseGetColumnVecRead(A, i, &a)); 3380 PetscCall(MatDenseGetColumnVecWrite(B, i, &b)); 3381 PetscCall(VecISCopy(b, h->is[0], SCATTER_FORWARD, a)); 3382 PetscCall(MatDenseRestoreColumnVecWrite(B, i, &b)); 3383 PetscCall(MatDenseRestoreColumnVecRead(A, i, &a)); 3384 } 3385 PetscCall(MatDestroy(&A)); 3386 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, h->ksp->pc->mat->rmap->n, B->cmap->n, nullptr, &A)); 3387 PetscCall(KSPMatSolve(h->ksp, B, A)); 3388 PetscCall(MatDestroy(&B)); 3389 for (PetscInt i = 0; i < A->cmap->n; ++i) { 3390 PetscCall(MatDenseGetColumnVecRead(A, i, &a)); 3391 PetscCall(MatDenseGetColumnVecWrite(Y, i, &b)); 3392 PetscCall(VecISCopy(a, h->is[1], SCATTER_REVERSE, b)); 3393 PetscCall(MatDenseRestoreColumnVecWrite(Y, i, &b)); 3394 PetscCall(MatDenseRestoreColumnVecRead(A, i, &a)); 3395 } 3396 PetscCall(MatDestroy(&A)); 3397 PetscFunctionReturn(PETSC_SUCCESS); 3398 } 3399 3400 static PetscErrorCode MatProduct_AtB_Harmonic(Mat S, Mat Y, Mat X, void *) 3401 { 3402 Harmonic h; 3403 Mat A, B; 3404 Vec a, b; 3405 3406 PetscFunctionBegin; 3407 PetscCall(MatShellGetContext(S, &h)); 3408 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, h->ksp->pc->mat->rmap->n, Y->cmap->n, nullptr, &A)); 3409 for (PetscInt i = 0; i < A->cmap->n; ++i) { 3410 PetscCall(MatDenseGetColumnVecRead(Y, i, &b)); 3411 PetscCall(MatDenseGetColumnVecWrite(A, i, &a)); 3412 PetscCall(VecISCopy(a, h->is[1], SCATTER_FORWARD, b)); 3413 PetscCall(MatDenseRestoreColumnVecWrite(A, i, &a)); 3414 PetscCall(MatDenseRestoreColumnVecRead(Y, i, &b)); 3415 } 3416 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, h->ksp->pc->mat->rmap->n, A->cmap->n, nullptr, &B)); 3417 PetscCall(KSPMatSolveTranspose(h->ksp, A, B)); 3418 PetscCall(MatDestroy(&A)); 3419 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, h->A[0]->rmap->n, B->cmap->n, nullptr, &A)); 3420 for (PetscInt i = 0; i < A->cmap->n; ++i) { 3421 PetscCall(MatDenseGetColumnVecRead(B, i, &b)); 3422 PetscCall(MatDenseGetColumnVecWrite(A, i, &a)); 3423 PetscCall(VecISCopy(b, h->is[0], SCATTER_REVERSE, a)); 3424 PetscCall(MatDenseRestoreColumnVecWrite(A, i, &a)); 3425 PetscCall(MatDenseRestoreColumnVecRead(B, i, &b)); 3426 } 3427 PetscCall(MatDestroy(&B)); 3428 PetscCall(MatTransposeMatMult(h->A[0], A, MAT_REUSE_MATRIX, PETSC_CURRENT, &X)); 3429 PetscCall(MatDestroy(&A)); 3430 PetscFunctionReturn(PETSC_SUCCESS); 3431 } 3432 3433 static PetscErrorCode MatDestroy_Harmonic(Mat A) 3434 { 3435 Harmonic h; 3436 3437 PetscFunctionBegin; 3438 PetscCall(MatShellGetContext(A, &h)); 3439 for (PetscInt i = 0; i < 5; ++i) PetscCall(ISDestroy(h->is + i)); 3440 PetscCall(PetscFree(h->is)); 3441 PetscCall(VecDestroy(&h->v)); 3442 for (PetscInt i = 0; i < 2; ++i) PetscCall(MatDestroy(h->A + i)); 3443 PetscCall(PetscFree(h->A)); 3444 PetscCall(KSPDestroy(&h->ksp)); 3445 PetscCall(PetscFree(h)); 3446 PetscFunctionReturn(PETSC_SUCCESS); 3447 } 3448