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