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