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