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