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