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