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