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