1 #include <petsc/private/dmimpl.h> 2 #include <petsc/private/matimpl.h> 3 #include <petsc/private/petschpddm.h> /*I "petscpc.h" I*/ 4 #include <petsc/private/pcimpl.h> /* this must be included after petschpddm.h so that _PCIMPL_H is not defined */ 5 /* otherwise, it is assumed that one is compiling libhpddm_petsc => circular dependency */ 6 #if PetscDefined(USE_FORTRAN_BINDINGS) 7 #include <petsc/private/fortranimpl.h> 8 #endif 9 10 static PetscErrorCode (*loadedSym)(HPDDM::Schwarz<PetscScalar> *const, IS, Mat, Mat, Mat, std::vector<Vec>, PC_HPDDM_Level **const) = nullptr; 11 12 static PetscBool PCHPDDMPackageInitialized = PETSC_FALSE; 13 14 PetscLogEvent PC_HPDDM_Strc; 15 PetscLogEvent PC_HPDDM_PtAP; 16 PetscLogEvent PC_HPDDM_PtBP; 17 PetscLogEvent PC_HPDDM_Next; 18 PetscLogEvent PC_HPDDM_SetUp[PETSC_PCHPDDM_MAXLEVELS]; 19 PetscLogEvent PC_HPDDM_Solve[PETSC_PCHPDDM_MAXLEVELS]; 20 21 const char *const PCHPDDMCoarseCorrectionTypes[] = {"DEFLATED", "ADDITIVE", "BALANCED", "PCHPDDMCoarseCorrectionType", "PC_HPDDM_COARSE_CORRECTION_", nullptr}; 22 23 static PetscErrorCode PCReset_HPDDM(PC pc) 24 { 25 PC_HPDDM *data = (PC_HPDDM *)pc->data; 26 PetscInt i; 27 28 PetscFunctionBegin; 29 if (data->levels) { 30 for (i = 0; i < PETSC_PCHPDDM_MAXLEVELS && data->levels[i]; ++i) { 31 PetscCall(KSPDestroy(&data->levels[i]->ksp)); 32 PetscCall(PCDestroy(&data->levels[i]->pc)); 33 PetscCall(PetscFree(data->levels[i])); 34 } 35 PetscCall(PetscFree(data->levels)); 36 } 37 38 PetscCall(ISDestroy(&data->is)); 39 PetscCall(MatDestroy(&data->aux)); 40 PetscCall(MatDestroy(&data->B)); 41 PetscCall(VecDestroy(&data->normal)); 42 data->correction = PC_HPDDM_COARSE_CORRECTION_DEFLATED; 43 data->Neumann = PETSC_BOOL3_UNKNOWN; 44 data->deflation = PETSC_FALSE; 45 data->setup = nullptr; 46 data->setup_ctx = nullptr; 47 PetscFunctionReturn(PETSC_SUCCESS); 48 } 49 50 static PetscErrorCode PCDestroy_HPDDM(PC pc) 51 { 52 PC_HPDDM *data = (PC_HPDDM *)pc->data; 53 54 PetscFunctionBegin; 55 PetscCall(PCReset_HPDDM(pc)); 56 PetscCall(PetscFree(data)); 57 PetscCall(PetscObjectChangeTypeName((PetscObject)pc, nullptr)); 58 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetAuxiliaryMat_C", nullptr)); 59 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMHasNeumannMat_C", nullptr)); 60 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetRHSMat_C", nullptr)); 61 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetCoarseCorrectionType_C", nullptr)); 62 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMGetCoarseCorrectionType_C", nullptr)); 63 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetSTShareSubKSP_C", nullptr)); 64 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMGetSTShareSubKSP_C", nullptr)); 65 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetDeflationMat_C", nullptr)); 66 PetscFunctionReturn(PETSC_SUCCESS); 67 } 68 69 static inline PetscErrorCode PCHPDDMSetAuxiliaryMat_Private(PC pc, IS is, Mat A, PetscBool deflation) 70 { 71 PC_HPDDM *data = (PC_HPDDM *)pc->data; 72 73 PetscFunctionBegin; 74 if (is) { 75 PetscCall(PetscObjectReference((PetscObject)is)); 76 if (data->is) { /* new overlap definition resets the PC */ 77 PetscCall(PCReset_HPDDM(pc)); 78 pc->setfromoptionscalled = 0; 79 pc->setupcalled = 0; 80 } 81 PetscCall(ISDestroy(&data->is)); 82 data->is = is; 83 } 84 if (A) { 85 PetscCall(PetscObjectReference((PetscObject)A)); 86 PetscCall(MatDestroy(&data->aux)); 87 data->aux = A; 88 } 89 data->deflation = deflation; 90 PetscFunctionReturn(PETSC_SUCCESS); 91 } 92 93 static inline PetscErrorCode PCHPDDMSetAuxiliaryMatNormal_Private(PC pc, Mat A, Mat N, Mat *B, const char *pcpre, Vec *diagonal = nullptr) 94 { 95 PC_HPDDM *data = (PC_HPDDM *)pc->data; 96 Mat *splitting, *sub, aux; 97 Vec d; 98 IS is, cols[2], rows; 99 PetscReal norm; 100 PetscBool flg; 101 char type[256] = {}; /* same size as in src/ksp/pc/interface/pcset.c */ 102 103 PetscFunctionBegin; 104 PetscCall(MatConvert(N, MATAIJ, MAT_INITIAL_MATRIX, B)); 105 PetscCall(ISCreateStride(PETSC_COMM_SELF, A->cmap->n, A->cmap->rstart, 1, cols)); 106 PetscCall(ISCreateStride(PETSC_COMM_SELF, A->rmap->N, 0, 1, &rows)); 107 PetscCall(MatSetOption(A, MAT_SUBMAT_SINGLEIS, PETSC_TRUE)); 108 PetscCall(MatIncreaseOverlap(*B, 1, cols, 1)); 109 PetscCall(MatCreateSubMatrices(A, 1, &rows, cols, MAT_INITIAL_MATRIX, &splitting)); 110 PetscCall(ISCreateStride(PETSC_COMM_SELF, A->cmap->n, A->cmap->rstart, 1, &is)); 111 PetscCall(ISEmbed(*cols, is, PETSC_TRUE, cols + 1)); 112 PetscCall(ISDestroy(&is)); 113 PetscCall(MatCreateSubMatrices(*splitting, 1, &rows, cols + 1, MAT_INITIAL_MATRIX, &sub)); 114 PetscCall(ISDestroy(cols + 1)); 115 PetscCall(MatFindZeroRows(*sub, &is)); 116 PetscCall(MatDestroySubMatrices(1, &sub)); 117 PetscCall(ISDestroy(&rows)); 118 PetscCall(MatSetOption(*splitting, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE)); 119 PetscCall(MatZeroRowsIS(*splitting, is, 0.0, nullptr, nullptr)); 120 PetscCall(ISDestroy(&is)); 121 PetscCall(PetscOptionsGetString(nullptr, pcpre, "-pc_hpddm_levels_1_sub_pc_type", type, sizeof(type), nullptr)); 122 PetscCall(PetscStrcmp(type, PCQR, &flg)); 123 if (!flg) { 124 Mat conjugate = *splitting; 125 if (PetscDefined(USE_COMPLEX)) { 126 PetscCall(MatDuplicate(*splitting, MAT_COPY_VALUES, &conjugate)); 127 PetscCall(MatConjugate(conjugate)); 128 } 129 PetscCall(MatTransposeMatMult(conjugate, *splitting, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &aux)); 130 if (PetscDefined(USE_COMPLEX)) PetscCall(MatDestroy(&conjugate)); 131 PetscCall(MatNorm(aux, NORM_FROBENIUS, &norm)); 132 PetscCall(MatSetOption(aux, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE)); 133 if (diagonal) { 134 PetscCall(VecNorm(*diagonal, NORM_INFINITY, &norm)); 135 if (norm > PETSC_SMALL) { 136 VecScatter scatter; 137 PetscInt n; 138 PetscCall(ISGetLocalSize(*cols, &n)); 139 PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)pc), n, PETSC_DECIDE, &d)); 140 PetscCall(VecScatterCreate(*diagonal, *cols, d, nullptr, &scatter)); 141 PetscCall(VecScatterBegin(scatter, *diagonal, d, INSERT_VALUES, SCATTER_FORWARD)); 142 PetscCall(VecScatterEnd(scatter, *diagonal, d, INSERT_VALUES, SCATTER_FORWARD)); 143 PetscCall(VecScatterDestroy(&scatter)); 144 PetscCall(MatScale(aux, -1.0)); 145 PetscCall(MatDiagonalSet(aux, d, ADD_VALUES)); 146 PetscCall(VecDestroy(&d)); 147 } else PetscCall(VecDestroy(diagonal)); 148 } 149 if (!diagonal) PetscCall(MatShift(aux, PETSC_SMALL * norm)); 150 } else { 151 PetscBool flg; 152 if (diagonal) { 153 PetscCall(VecNorm(*diagonal, NORM_INFINITY, &norm)); 154 PetscCheck(norm < PETSC_SMALL, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Nonzero diagonal A11 block"); 155 PetscCall(VecDestroy(diagonal)); 156 } 157 PetscCall(PetscObjectTypeCompare((PetscObject)N, MATNORMAL, &flg)); 158 if (flg) PetscCall(MatCreateNormal(*splitting, &aux)); 159 else PetscCall(MatCreateNormalHermitian(*splitting, &aux)); 160 } 161 PetscCall(MatDestroySubMatrices(1, &splitting)); 162 PetscCall(PCHPDDMSetAuxiliaryMat(pc, *cols, aux, nullptr, nullptr)); 163 data->Neumann = PETSC_BOOL3_TRUE; 164 PetscCall(ISDestroy(cols)); 165 PetscCall(MatDestroy(&aux)); 166 PetscFunctionReturn(PETSC_SUCCESS); 167 } 168 169 static PetscErrorCode PCHPDDMSetAuxiliaryMat_HPDDM(PC pc, IS is, Mat A, PetscErrorCode (*setup)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void *), void *setup_ctx) 170 { 171 PC_HPDDM *data = (PC_HPDDM *)pc->data; 172 173 PetscFunctionBegin; 174 PetscCall(PCHPDDMSetAuxiliaryMat_Private(pc, is, A, PETSC_FALSE)); 175 if (setup) { 176 data->setup = setup; 177 data->setup_ctx = setup_ctx; 178 } 179 PetscFunctionReturn(PETSC_SUCCESS); 180 } 181 182 /*@C 183 PCHPDDMSetAuxiliaryMat - Sets the auxiliary matrix used by `PCHPDDM` for the concurrent GenEO problems at the finest level. 184 185 Input Parameters: 186 + pc - preconditioner context 187 . is - index set of the local auxiliary, e.g., Neumann, matrix 188 . A - auxiliary sequential matrix 189 . setup - function for generating the auxiliary matrix entries, may be `NULL` 190 - ctx - context for `setup`, may be `NULL` 191 192 Calling sequence of `setup`: 193 + J - matrix whose values are to be set 194 . t - time 195 . X - linearization point 196 . X_t - time-derivative of the linearization point 197 . s - step 198 . ovl - index set of the local auxiliary, e.g., Neumann, matrix 199 - ctx - context for `setup`, may be `NULL` 200 201 Level: intermediate 202 203 Note: 204 As an example, in a finite element context with nonoverlapping subdomains plus (overlapping) ghost elements, this could be the unassembled (Neumann) 205 local overlapping operator. As opposed to the assembled (Dirichlet) local overlapping operator obtained by summing neighborhood contributions 206 at the interface of ghost elements. 207 208 Fortran Notes: 209 Only `PETSC_NULL_FUNCTION` is supported for `setup` and `ctx` is never accessed 210 211 .seealso: `PCHPDDM`, `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCHPDDMSetRHSMat()`, `MATIS` 212 @*/ 213 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) 214 { 215 PetscFunctionBegin; 216 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 217 if (is) PetscValidHeaderSpecific(is, IS_CLASSID, 2); 218 if (A) PetscValidHeaderSpecific(A, MAT_CLASSID, 3); 219 PetscTryMethod(pc, "PCHPDDMSetAuxiliaryMat_C", (PC, IS, Mat, PetscErrorCode(*)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void *), void *), (pc, is, A, setup, ctx)); 220 PetscFunctionReturn(PETSC_SUCCESS); 221 } 222 223 static PetscErrorCode PCHPDDMHasNeumannMat_HPDDM(PC pc, PetscBool has) 224 { 225 PC_HPDDM *data = (PC_HPDDM *)pc->data; 226 227 PetscFunctionBegin; 228 data->Neumann = PetscBoolToBool3(has); 229 PetscFunctionReturn(PETSC_SUCCESS); 230 } 231 232 /*@ 233 PCHPDDMHasNeumannMat - Informs `PCHPDDM` that the `Mat` passed to `PCHPDDMSetAuxiliaryMat()` is the local Neumann matrix. 234 235 Input Parameters: 236 + pc - preconditioner context 237 - has - Boolean value 238 239 Level: intermediate 240 241 Notes: 242 This may be used to bypass a call to `MatCreateSubMatrices()` and to `MatConvert()` for `MATSBAIJ` matrices. 243 244 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`. 245 246 .seealso: `PCHPDDM`, `PCHPDDMSetAuxiliaryMat()` 247 @*/ 248 PetscErrorCode PCHPDDMHasNeumannMat(PC pc, PetscBool has) 249 { 250 PetscFunctionBegin; 251 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 252 PetscTryMethod(pc, "PCHPDDMHasNeumannMat_C", (PC, PetscBool), (pc, has)); 253 PetscFunctionReturn(PETSC_SUCCESS); 254 } 255 256 static PetscErrorCode PCHPDDMSetRHSMat_HPDDM(PC pc, Mat B) 257 { 258 PC_HPDDM *data = (PC_HPDDM *)pc->data; 259 260 PetscFunctionBegin; 261 PetscCall(PetscObjectReference((PetscObject)B)); 262 PetscCall(MatDestroy(&data->B)); 263 data->B = B; 264 PetscFunctionReturn(PETSC_SUCCESS); 265 } 266 267 /*@ 268 PCHPDDMSetRHSMat - Sets the right-hand side matrix used by `PCHPDDM` for the concurrent GenEO problems at the finest level. 269 270 Input Parameters: 271 + pc - preconditioner context 272 - B - right-hand side sequential matrix 273 274 Level: advanced 275 276 Note: 277 Must be used in conjunction with `PCHPDDMSetAuxiliaryMat`(N), so that Nv = lambda Bv is solved using `EPSSetOperators`(N, B). 278 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. 279 280 .seealso: `PCHPDDMSetAuxiliaryMat()`, `PCHPDDM` 281 @*/ 282 PetscErrorCode PCHPDDMSetRHSMat(PC pc, Mat B) 283 { 284 PetscFunctionBegin; 285 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 286 if (B) { 287 PetscValidHeaderSpecific(B, MAT_CLASSID, 2); 288 PetscTryMethod(pc, "PCHPDDMSetRHSMat_C", (PC, Mat), (pc, B)); 289 } 290 PetscFunctionReturn(PETSC_SUCCESS); 291 } 292 293 static PetscErrorCode PCSetFromOptions_HPDDM(PC pc, PetscOptionItems *PetscOptionsObject) 294 { 295 PC_HPDDM *data = (PC_HPDDM *)pc->data; 296 PC_HPDDM_Level **levels = data->levels; 297 char prefix[256]; 298 int i = 1; 299 PetscMPIInt size, previous; 300 PetscInt n; 301 PCHPDDMCoarseCorrectionType type; 302 PetscBool flg = PETSC_TRUE, set; 303 304 PetscFunctionBegin; 305 if (!data->levels) { 306 PetscCall(PetscCalloc1(PETSC_PCHPDDM_MAXLEVELS, &levels)); 307 data->levels = levels; 308 } 309 PetscOptionsHeadBegin(PetscOptionsObject, "PCHPDDM options"); 310 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size)); 311 previous = size; 312 while (i < PETSC_PCHPDDM_MAXLEVELS) { 313 PetscInt p = 1; 314 315 if (!data->levels[i - 1]) PetscCall(PetscNew(data->levels + i - 1)); 316 data->levels[i - 1]->parent = data; 317 /* if the previous level has a single process, it is not possible to coarsen further */ 318 if (previous == 1 || !flg) break; 319 data->levels[i - 1]->nu = 0; 320 data->levels[i - 1]->threshold = -1.0; 321 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_%d_eps_nev", i)); 322 PetscCall(PetscOptionsBoundedInt(prefix, "Local number of deflation vectors computed by SLEPc", "EPSSetDimensions", data->levels[i - 1]->nu, &data->levels[i - 1]->nu, nullptr, 0)); 323 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_%d_eps_threshold", i)); 324 PetscCall(PetscOptionsReal(prefix, "Local threshold for selecting deflation vectors returned by SLEPc", "PCHPDDM", data->levels[i - 1]->threshold, &data->levels[i - 1]->threshold, nullptr)); 325 if (i == 1) { 326 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_1_st_share_sub_ksp")); 327 PetscCall(PetscOptionsBool(prefix, "Shared KSP between SLEPc ST and the fine-level subdomain solver", "PCHPDDMSetSTShareSubKSP", PETSC_FALSE, &data->share, nullptr)); 328 } 329 /* if there is no prescribed coarsening, just break out of the loop */ 330 if (data->levels[i - 1]->threshold <= 0.0 && data->levels[i - 1]->nu <= 0 && !(data->deflation && i == 1)) break; 331 else { 332 ++i; 333 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_%d_eps_nev", i)); 334 PetscCall(PetscOptionsHasName(PetscOptionsObject->options, PetscOptionsObject->prefix, prefix, &flg)); 335 if (!flg) { 336 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_%d_eps_threshold", i)); 337 PetscCall(PetscOptionsHasName(PetscOptionsObject->options, PetscOptionsObject->prefix, prefix, &flg)); 338 } 339 if (flg) { 340 /* if there are coarsening options for the next level, then register it */ 341 /* otherwise, don't to avoid having both options levels_N_p and coarse_p */ 342 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_%d_p", i)); 343 PetscCall(PetscOptionsRangeInt(prefix, "Number of processes used to assemble the coarse operator at this level", "PCHPDDM", p, &p, &flg, 1, PetscMax(1, previous / 2))); 344 previous = p; 345 } 346 } 347 } 348 data->N = i; 349 n = 1; 350 if (i > 1) { 351 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_coarse_p")); 352 PetscCall(PetscOptionsRangeInt(prefix, "Number of processes used to assemble the coarsest operator", "PCHPDDM", n, &n, nullptr, 1, PetscMax(1, previous / 2))); 353 #if PetscDefined(HAVE_MUMPS) 354 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "pc_hpddm_coarse_")); 355 PetscCall(PetscOptionsHasName(nullptr, prefix, "-mat_mumps_use_omp_threads", &flg)); 356 if (flg) { 357 char type[64]; /* same size as in src/ksp/pc/impls/factor/factimpl.c */ 358 PetscCall(PetscStrncpy(type, n > 1 && PetscDefined(HAVE_MUMPS) ? MATSOLVERMUMPS : MATSOLVERPETSC, sizeof(type))); /* default solver for a MatMPIAIJ or a MatSeqAIJ */ 359 PetscCall(PetscOptionsGetString(nullptr, prefix, "-pc_factor_mat_solver_type", type, sizeof(type), nullptr)); 360 PetscCall(PetscStrcmp(type, MATSOLVERMUMPS, &flg)); 361 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "-%smat_mumps_use_omp_threads and -%spc_factor_mat_solver_type != %s", prefix, prefix, MATSOLVERMUMPS); 362 size = n; 363 n = -1; 364 PetscCall(PetscOptionsGetInt(nullptr, prefix, "-mat_mumps_use_omp_threads", &n, nullptr)); 365 PetscCheck(n >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Need to specify a positive integer for -%smat_mumps_use_omp_threads", prefix); 366 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" : ""); 367 } 368 #endif 369 PetscCall(PetscOptionsEnum("-pc_hpddm_coarse_correction", "Type of coarse correction applied each iteration", "PCHPDDMSetCoarseCorrectionType", PCHPDDMCoarseCorrectionTypes, (PetscEnum)data->correction, (PetscEnum *)&type, &flg)); 370 if (flg) PetscCall(PCHPDDMSetCoarseCorrectionType(pc, type)); 371 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_has_neumann")); 372 PetscCall(PetscOptionsBool(prefix, "Is the auxiliary Mat the local Neumann matrix?", "PCHPDDMHasNeumannMat", PetscBool3ToBool(data->Neumann), &flg, &set)); 373 if (set) data->Neumann = PetscBoolToBool3(flg); 374 data->log_separate = PETSC_FALSE; 375 if (PetscDefined(USE_LOG)) { 376 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_log_separate")); 377 PetscCall(PetscOptionsBool(prefix, "Log events level by level instead of inside PCSetUp()/KSPSolve()", nullptr, data->log_separate, &data->log_separate, nullptr)); 378 } 379 } 380 PetscOptionsHeadEnd(); 381 while (i < PETSC_PCHPDDM_MAXLEVELS && data->levels[i]) PetscCall(PetscFree(data->levels[i++])); 382 PetscFunctionReturn(PETSC_SUCCESS); 383 } 384 385 static PetscErrorCode PCApply_HPDDM(PC pc, Vec x, Vec y) 386 { 387 PC_HPDDM *data = (PC_HPDDM *)pc->data; 388 389 PetscFunctionBegin; 390 PetscCall(PetscCitationsRegister(HPDDMCitation, &HPDDMCite)); 391 PetscCheck(data->levels[0]->ksp, PETSC_COMM_SELF, PETSC_ERR_PLIB, "No KSP attached to PCHPDDM"); 392 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 */ 393 PetscCall(KSPSolve(data->levels[0]->ksp, x, y)); 394 if (data->log_separate) PetscCall(PetscLogEventEnd(PC_HPDDM_Solve[0], data->levels[0]->ksp, nullptr, nullptr, nullptr)); 395 PetscFunctionReturn(PETSC_SUCCESS); 396 } 397 398 static PetscErrorCode PCMatApply_HPDDM(PC pc, Mat X, Mat Y) 399 { 400 PC_HPDDM *data = (PC_HPDDM *)pc->data; 401 402 PetscFunctionBegin; 403 PetscCall(PetscCitationsRegister(HPDDMCitation, &HPDDMCite)); 404 PetscCheck(data->levels[0]->ksp, PETSC_COMM_SELF, PETSC_ERR_PLIB, "No KSP attached to PCHPDDM"); 405 PetscCall(KSPMatSolve(data->levels[0]->ksp, X, Y)); 406 PetscFunctionReturn(PETSC_SUCCESS); 407 } 408 409 // PetscClangLinter pragma disable: -fdoc-internal-linkage 410 /*@C 411 PCHPDDMGetComplexities - Computes the grid and operator complexities. 412 413 Input Parameter: 414 . pc - preconditioner context 415 416 Output Parameters: 417 + gc - grid complexity = sum_i(m_i) / m_1 418 - oc - operator complexity = sum_i(nnz_i) / nnz_1 419 420 Level: advanced 421 422 .seealso: `PCMGGetGridComplexity()`, `PCHPDDM`, `PCHYPRE`, `PCGAMG` 423 @*/ 424 static PetscErrorCode PCHPDDMGetComplexities(PC pc, PetscReal *gc, PetscReal *oc) 425 { 426 PC_HPDDM *data = (PC_HPDDM *)pc->data; 427 MatInfo info; 428 PetscInt n, m; 429 PetscLogDouble accumulate[2]{}, nnz1 = 1.0, m1 = 1.0; 430 431 PetscFunctionBegin; 432 for (n = 0, *gc = 0, *oc = 0; n < data->N; ++n) { 433 if (data->levels[n]->ksp) { 434 Mat P, A; 435 PetscCall(KSPGetOperators(data->levels[n]->ksp, nullptr, &P)); 436 PetscCall(MatGetSize(P, &m, nullptr)); 437 accumulate[0] += m; 438 if (n == 0) { 439 PetscBool flg; 440 PetscCall(PetscObjectTypeCompareAny((PetscObject)P, &flg, MATNORMAL, MATNORMALHERMITIAN, "")); 441 if (flg) { 442 PetscCall(MatConvert(P, MATAIJ, MAT_INITIAL_MATRIX, &A)); 443 P = A; 444 } else PetscCall(PetscObjectReference((PetscObject)P)); 445 } 446 if (P->ops->getinfo) { 447 PetscCall(MatGetInfo(P, MAT_GLOBAL_SUM, &info)); 448 accumulate[1] += info.nz_used; 449 } 450 if (n == 0) { 451 m1 = m; 452 if (P->ops->getinfo) nnz1 = info.nz_used; 453 PetscCall(MatDestroy(&P)); 454 } 455 } 456 } 457 *gc = static_cast<PetscReal>(accumulate[0] / m1); 458 *oc = static_cast<PetscReal>(accumulate[1] / nnz1); 459 PetscFunctionReturn(PETSC_SUCCESS); 460 } 461 462 static PetscErrorCode PCView_HPDDM(PC pc, PetscViewer viewer) 463 { 464 PC_HPDDM *data = (PC_HPDDM *)pc->data; 465 PetscViewer subviewer; 466 PetscSubcomm subcomm; 467 PetscReal oc, gc; 468 PetscInt i, tabs; 469 PetscMPIInt size, color, rank; 470 PetscBool ascii; 471 472 PetscFunctionBegin; 473 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &ascii)); 474 if (ascii) { 475 PetscCall(PetscViewerASCIIPrintf(viewer, "level%s: %" PetscInt_FMT "\n", data->N > 1 ? "s" : "", data->N)); 476 PetscCall(PCHPDDMGetComplexities(pc, &gc, &oc)); 477 if (data->N > 1) { 478 if (!data->deflation) { 479 PetscCall(PetscViewerASCIIPrintf(viewer, "Neumann matrix attached? %s\n", PetscBools[PetscBool3ToBool(data->Neumann)])); 480 PetscCall(PetscViewerASCIIPrintf(viewer, "shared subdomain KSP between SLEPc and PETSc? %s\n", PetscBools[data->share])); 481 } else PetscCall(PetscViewerASCIIPrintf(viewer, "user-supplied deflation matrix\n")); 482 PetscCall(PetscViewerASCIIPrintf(viewer, "coarse correction: %s\n", PCHPDDMCoarseCorrectionTypes[data->correction])); 483 PetscCall(PetscViewerASCIIPrintf(viewer, "on process #0, value%s (+ threshold%s if available) for selecting deflation vectors:", data->N > 2 ? "s" : "", data->N > 2 ? "s" : "")); 484 PetscCall(PetscViewerASCIIGetTab(viewer, &tabs)); 485 PetscCall(PetscViewerASCIISetTab(viewer, 0)); 486 for (i = 1; i < data->N; ++i) { 487 PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, data->levels[i - 1]->nu)); 488 if (data->levels[i - 1]->threshold > -0.1) PetscCall(PetscViewerASCIIPrintf(viewer, " (%g)", (double)data->levels[i - 1]->threshold)); 489 } 490 PetscCall(PetscViewerASCIIPrintf(viewer, "\n")); 491 PetscCall(PetscViewerASCIISetTab(viewer, tabs)); 492 } 493 PetscCall(PetscViewerASCIIPrintf(viewer, "grid and operator complexities: %g %g\n", (double)gc, (double)oc)); 494 if (data->levels[0]->ksp) { 495 PetscCall(KSPView(data->levels[0]->ksp, viewer)); 496 if (data->levels[0]->pc) PetscCall(PCView(data->levels[0]->pc, viewer)); 497 for (i = 1; i < data->N; ++i) { 498 if (data->levels[i]->ksp) color = 1; 499 else color = 0; 500 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size)); 501 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank)); 502 PetscCall(PetscSubcommCreate(PetscObjectComm((PetscObject)pc), &subcomm)); 503 PetscCall(PetscSubcommSetNumber(subcomm, PetscMin(size, 2))); 504 PetscCall(PetscSubcommSetTypeGeneral(subcomm, color, rank)); 505 PetscCall(PetscViewerASCIIPushTab(viewer)); 506 PetscCall(PetscViewerGetSubViewer(viewer, PetscSubcommChild(subcomm), &subviewer)); 507 if (color == 1) { 508 PetscCall(KSPView(data->levels[i]->ksp, subviewer)); 509 if (data->levels[i]->pc) PetscCall(PCView(data->levels[i]->pc, subviewer)); 510 PetscCall(PetscViewerFlush(subviewer)); 511 } 512 PetscCall(PetscViewerRestoreSubViewer(viewer, PetscSubcommChild(subcomm), &subviewer)); 513 PetscCall(PetscViewerASCIIPopTab(viewer)); 514 PetscCall(PetscSubcommDestroy(&subcomm)); 515 PetscCall(PetscViewerFlush(viewer)); 516 } 517 } 518 } 519 PetscFunctionReturn(PETSC_SUCCESS); 520 } 521 522 static PetscErrorCode PCPreSolve_HPDDM(PC pc, KSP ksp, Vec, Vec) 523 { 524 PC_HPDDM *data = (PC_HPDDM *)pc->data; 525 PetscBool flg; 526 Mat A; 527 528 PetscFunctionBegin; 529 if (ksp) { 530 PetscCall(PetscObjectTypeCompare((PetscObject)ksp, KSPLSQR, &flg)); 531 if (flg && !data->normal) { 532 PetscCall(KSPGetOperators(ksp, &A, nullptr)); 533 PetscCall(MatCreateVecs(A, nullptr, &data->normal)); /* temporary Vec used in PCApply_HPDDMShell() for coarse grid corrections */ 534 } 535 } 536 PetscFunctionReturn(PETSC_SUCCESS); 537 } 538 539 static PetscErrorCode PCSetUp_HPDDMShell(PC pc) 540 { 541 PC_HPDDM_Level *ctx; 542 Mat A, P; 543 Vec x; 544 const char *pcpre; 545 546 PetscFunctionBegin; 547 PetscCall(PCShellGetContext(pc, &ctx)); 548 PetscCall(KSPGetOptionsPrefix(ctx->ksp, &pcpre)); 549 PetscCall(KSPGetOperators(ctx->ksp, &A, &P)); 550 /* smoother */ 551 PetscCall(PCSetOptionsPrefix(ctx->pc, pcpre)); 552 PetscCall(PCSetOperators(ctx->pc, A, P)); 553 if (!ctx->v[0]) { 554 PetscCall(VecDuplicateVecs(ctx->D, 1, &ctx->v[0])); 555 if (!std::is_same<PetscScalar, PetscReal>::value) PetscCall(VecDestroy(&ctx->D)); 556 PetscCall(MatCreateVecs(A, &x, nullptr)); 557 PetscCall(VecDuplicateVecs(x, 2, &ctx->v[1])); 558 PetscCall(VecDestroy(&x)); 559 } 560 std::fill_n(ctx->V, 3, nullptr); 561 PetscFunctionReturn(PETSC_SUCCESS); 562 } 563 564 template <class Type, typename std::enable_if<std::is_same<Type, Vec>::value>::type * = nullptr> 565 static inline PetscErrorCode PCHPDDMDeflate_Private(PC pc, Type x, Type y) 566 { 567 PC_HPDDM_Level *ctx; 568 PetscScalar *out; 569 570 PetscFunctionBegin; 571 PetscCall(PCShellGetContext(pc, &ctx)); 572 /* going from PETSc to HPDDM numbering */ 573 PetscCall(VecScatterBegin(ctx->scatter, x, ctx->v[0][0], INSERT_VALUES, SCATTER_FORWARD)); 574 PetscCall(VecScatterEnd(ctx->scatter, x, ctx->v[0][0], INSERT_VALUES, SCATTER_FORWARD)); 575 PetscCall(VecGetArrayWrite(ctx->v[0][0], &out)); 576 ctx->P->deflation<false>(nullptr, out, 1); /* y = Q x */ 577 PetscCall(VecRestoreArrayWrite(ctx->v[0][0], &out)); 578 /* going from HPDDM to PETSc numbering */ 579 PetscCall(VecScatterBegin(ctx->scatter, ctx->v[0][0], y, INSERT_VALUES, SCATTER_REVERSE)); 580 PetscCall(VecScatterEnd(ctx->scatter, ctx->v[0][0], y, INSERT_VALUES, SCATTER_REVERSE)); 581 PetscFunctionReturn(PETSC_SUCCESS); 582 } 583 584 template <class Type, typename std::enable_if<std::is_same<Type, Mat>::value>::type * = nullptr> 585 static inline PetscErrorCode PCHPDDMDeflate_Private(PC pc, Type X, Type Y) 586 { 587 PC_HPDDM_Level *ctx; 588 Vec vX, vY, vC; 589 PetscScalar *out; 590 PetscInt i, N; 591 592 PetscFunctionBegin; 593 PetscCall(PCShellGetContext(pc, &ctx)); 594 PetscCall(MatGetSize(X, nullptr, &N)); 595 /* going from PETSc to HPDDM numbering */ 596 for (i = 0; i < N; ++i) { 597 PetscCall(MatDenseGetColumnVecRead(X, i, &vX)); 598 PetscCall(MatDenseGetColumnVecWrite(ctx->V[0], i, &vC)); 599 PetscCall(VecScatterBegin(ctx->scatter, vX, vC, INSERT_VALUES, SCATTER_FORWARD)); 600 PetscCall(VecScatterEnd(ctx->scatter, vX, vC, INSERT_VALUES, SCATTER_FORWARD)); 601 PetscCall(MatDenseRestoreColumnVecWrite(ctx->V[0], i, &vC)); 602 PetscCall(MatDenseRestoreColumnVecRead(X, i, &vX)); 603 } 604 PetscCall(MatDenseGetArrayWrite(ctx->V[0], &out)); 605 ctx->P->deflation<false>(nullptr, out, N); /* Y = Q X */ 606 PetscCall(MatDenseRestoreArrayWrite(ctx->V[0], &out)); 607 /* going from HPDDM to PETSc numbering */ 608 for (i = 0; i < N; ++i) { 609 PetscCall(MatDenseGetColumnVecRead(ctx->V[0], i, &vC)); 610 PetscCall(MatDenseGetColumnVecWrite(Y, i, &vY)); 611 PetscCall(VecScatterBegin(ctx->scatter, vC, vY, INSERT_VALUES, SCATTER_REVERSE)); 612 PetscCall(VecScatterEnd(ctx->scatter, vC, vY, INSERT_VALUES, SCATTER_REVERSE)); 613 PetscCall(MatDenseRestoreColumnVecWrite(Y, i, &vY)); 614 PetscCall(MatDenseRestoreColumnVecRead(ctx->V[0], i, &vC)); 615 } 616 PetscFunctionReturn(PETSC_SUCCESS); 617 } 618 619 /* 620 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. 621 622 .vb 623 (1) y = Pmat^-1 x + Q x, 624 (2) y = Pmat^-1 (I - Amat Q) x + Q x (default), 625 (3) y = (I - Q Amat^T) Pmat^-1 (I - Amat Q) x + Q x. 626 .ve 627 628 Input Parameters: 629 + pc - preconditioner context 630 - x - input vector 631 632 Output Parameter: 633 . y - output vector 634 635 Notes: 636 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. 637 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. 638 (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. 639 640 Level: advanced 641 642 Developer Note: 643 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 644 to trigger it. Likely the manual page is `PCHPDDM` 645 646 .seealso: `PCHPDDM`, `PCHPDDMCoarseCorrectionType` 647 */ 648 static PetscErrorCode PCApply_HPDDMShell(PC pc, Vec x, Vec y) 649 { 650 PC_HPDDM_Level *ctx; 651 Mat A; 652 653 PetscFunctionBegin; 654 PetscCall(PCShellGetContext(pc, &ctx)); 655 PetscCheck(ctx->P, PETSC_COMM_SELF, PETSC_ERR_PLIB, "PCSHELL from PCHPDDM called with no HPDDM object"); 656 PetscCall(KSPGetOperators(ctx->ksp, &A, nullptr)); 657 PetscCall(PCHPDDMDeflate_Private(pc, x, y)); /* y = Q x */ 658 if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_DEFLATED || ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_BALANCED) { 659 if (!ctx->parent->normal || ctx != ctx->parent->levels[0]) PetscCall(MatMult(A, y, ctx->v[1][0])); /* y = A Q x */ 660 else { /* KSPLSQR and finest level */ PetscCall(MatMult(A, y, ctx->parent->normal)); /* y = A Q x */ 661 PetscCall(MatMultHermitianTranspose(A, ctx->parent->normal, ctx->v[1][0])); /* y = A^T A Q x */ 662 } 663 PetscCall(VecWAXPY(ctx->v[1][1], -1.0, ctx->v[1][0], x)); /* y = (I - A Q) x */ 664 PetscCall(PCApply(ctx->pc, ctx->v[1][1], ctx->v[1][0])); /* y = M^-1 (I - A Q) x */ 665 if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_BALANCED) { 666 if (!ctx->parent->normal || ctx != ctx->parent->levels[0]) PetscCall(MatMultHermitianTranspose(A, ctx->v[1][0], ctx->v[1][1])); /* z = A^T y */ 667 else { 668 PetscCall(MatMult(A, ctx->v[1][0], ctx->parent->normal)); 669 PetscCall(MatMultHermitianTranspose(A, ctx->parent->normal, ctx->v[1][1])); /* z = A^T A y */ 670 } 671 PetscCall(PCHPDDMDeflate_Private(pc, ctx->v[1][1], ctx->v[1][1])); 672 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 */ 673 } else PetscCall(VecAXPY(y, 1.0, ctx->v[1][0])); /* y = Q M^-1 (I - A Q) x + Q x */ 674 } else { 675 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); 676 PetscCall(PCApply(ctx->pc, x, ctx->v[1][0])); 677 PetscCall(VecAXPY(y, 1.0, ctx->v[1][0])); /* y = M^-1 x + Q x */ 678 } 679 PetscFunctionReturn(PETSC_SUCCESS); 680 } 681 682 /* 683 PCMatApply_HPDDMShell - Variant of PCApply_HPDDMShell() for blocks of vectors. 684 685 Input Parameters: 686 + pc - preconditioner context 687 - X - block of input vectors 688 689 Output Parameter: 690 . Y - block of output vectors 691 692 Level: advanced 693 694 .seealso: `PCHPDDM`, `PCApply_HPDDMShell()`, `PCHPDDMCoarseCorrectionType` 695 */ 696 static PetscErrorCode PCMatApply_HPDDMShell(PC pc, Mat X, Mat Y) 697 { 698 PC_HPDDM_Level *ctx; 699 Mat A, *ptr; 700 PetscContainer container = nullptr; 701 PetscScalar *array; 702 PetscInt m, M, N, prev = 0; 703 PetscBool reset = PETSC_FALSE; 704 705 PetscFunctionBegin; 706 PetscCall(PCShellGetContext(pc, &ctx)); 707 PetscCheck(ctx->P, PETSC_COMM_SELF, PETSC_ERR_PLIB, "PCSHELL from PCHPDDM called with no HPDDM object"); 708 PetscCall(MatGetSize(X, nullptr, &N)); 709 PetscCall(KSPGetOperators(ctx->ksp, &A, nullptr)); 710 PetscCall(PetscObjectQuery((PetscObject)A, "_HPDDM_MatProduct", (PetscObject *)&container)); 711 if (container) { /* MatProduct container already attached */ 712 PetscCall(PetscContainerGetPointer(container, (void **)&ptr)); 713 if (ptr[1] != ctx->V[2]) /* Mat has changed or may have been set first in KSPHPDDM */ 714 for (m = 0; m < 2; ++m) { 715 PetscCall(MatDestroy(ctx->V + m + 1)); 716 ctx->V[m + 1] = ptr[m]; 717 PetscCall(PetscObjectReference((PetscObject)ctx->V[m + 1])); 718 } 719 } 720 if (ctx->V[1]) PetscCall(MatGetSize(ctx->V[1], nullptr, &prev)); 721 if (N != prev || !ctx->V[0]) { 722 PetscCall(MatDestroy(ctx->V)); 723 PetscCall(VecGetLocalSize(ctx->v[0][0], &m)); 724 PetscCall(MatCreateDense(PetscObjectComm((PetscObject)pc), m, PETSC_DECIDE, PETSC_DECIDE, N, nullptr, ctx->V)); 725 if (N != prev) { 726 PetscCall(MatDestroy(ctx->V + 1)); 727 PetscCall(MatDestroy(ctx->V + 2)); 728 PetscCall(MatGetLocalSize(X, &m, nullptr)); 729 PetscCall(MatGetSize(X, &M, nullptr)); 730 if (ctx->parent->correction != PC_HPDDM_COARSE_CORRECTION_BALANCED) PetscCall(MatDenseGetArrayWrite(ctx->V[0], &array)); 731 else array = nullptr; 732 PetscCall(MatCreateDense(PetscObjectComm((PetscObject)pc), m, PETSC_DECIDE, M, N, array, ctx->V + 1)); 733 if (ctx->parent->correction != PC_HPDDM_COARSE_CORRECTION_BALANCED) PetscCall(MatDenseRestoreArrayWrite(ctx->V[0], &array)); 734 PetscCall(MatCreateDense(PetscObjectComm((PetscObject)pc), m, PETSC_DECIDE, M, N, nullptr, ctx->V + 2)); 735 PetscCall(MatProductCreateWithMat(A, Y, nullptr, ctx->V[1])); 736 PetscCall(MatProductSetType(ctx->V[1], MATPRODUCT_AB)); 737 PetscCall(MatProductSetFromOptions(ctx->V[1])); 738 PetscCall(MatProductSymbolic(ctx->V[1])); 739 if (!container) { /* no MatProduct container attached, create one to be queried in KSPHPDDM or at the next call to PCMatApply() */ 740 PetscCall(PetscContainerCreate(PetscObjectComm((PetscObject)A), &container)); 741 PetscCall(PetscObjectCompose((PetscObject)A, "_HPDDM_MatProduct", (PetscObject)container)); 742 } 743 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 */ 744 } 745 if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_BALANCED) { 746 PetscCall(MatProductCreateWithMat(A, ctx->V[1], nullptr, ctx->V[2])); 747 PetscCall(MatProductSetType(ctx->V[2], MATPRODUCT_AtB)); 748 PetscCall(MatProductSetFromOptions(ctx->V[2])); 749 PetscCall(MatProductSymbolic(ctx->V[2])); 750 } 751 ctx->P->start(N); 752 } 753 if (N == prev || container) { /* when MatProduct container is attached, always need to MatProductReplaceMats() since KSPHPDDM may have replaced the Mat as well */ 754 PetscCall(MatProductReplaceMats(nullptr, Y, nullptr, ctx->V[1])); 755 if (container && ctx->parent->correction != PC_HPDDM_COARSE_CORRECTION_BALANCED) { 756 PetscCall(MatDenseGetArrayWrite(ctx->V[0], &array)); 757 PetscCall(MatDensePlaceArray(ctx->V[1], array)); 758 PetscCall(MatDenseRestoreArrayWrite(ctx->V[0], &array)); 759 reset = PETSC_TRUE; 760 } 761 } 762 PetscCall(PCHPDDMDeflate_Private(pc, X, Y)); 763 if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_DEFLATED || ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_BALANCED) { 764 PetscCall(MatProductNumeric(ctx->V[1])); 765 PetscCall(MatCopy(ctx->V[1], ctx->V[2], SAME_NONZERO_PATTERN)); 766 PetscCall(MatAXPY(ctx->V[2], -1.0, X, SAME_NONZERO_PATTERN)); 767 PetscCall(PCMatApply(ctx->pc, ctx->V[2], ctx->V[1])); 768 if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_BALANCED) { 769 PetscCall(MatProductNumeric(ctx->V[2])); 770 PetscCall(PCHPDDMDeflate_Private(pc, ctx->V[2], ctx->V[2])); 771 PetscCall(MatAXPY(ctx->V[1], -1.0, ctx->V[2], SAME_NONZERO_PATTERN)); 772 } 773 PetscCall(MatAXPY(Y, -1.0, ctx->V[1], SAME_NONZERO_PATTERN)); 774 } else { 775 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); 776 PetscCall(PCMatApply(ctx->pc, X, ctx->V[1])); 777 PetscCall(MatAXPY(Y, 1.0, ctx->V[1], SAME_NONZERO_PATTERN)); 778 } 779 if (reset) PetscCall(MatDenseResetArray(ctx->V[1])); 780 PetscFunctionReturn(PETSC_SUCCESS); 781 } 782 783 static PetscErrorCode PCDestroy_HPDDMShell(PC pc) 784 { 785 PC_HPDDM_Level *ctx; 786 PetscContainer container; 787 788 PetscFunctionBegin; 789 PetscCall(PCShellGetContext(pc, &ctx)); 790 PetscCall(HPDDM::Schwarz<PetscScalar>::destroy(ctx, PETSC_TRUE)); 791 PetscCall(VecDestroyVecs(1, &ctx->v[0])); 792 PetscCall(VecDestroyVecs(2, &ctx->v[1])); 793 PetscCall(PetscObjectQuery((PetscObject)(ctx->pc)->mat, "_HPDDM_MatProduct", (PetscObject *)&container)); 794 PetscCall(PetscContainerDestroy(&container)); 795 PetscCall(PetscObjectCompose((PetscObject)(ctx->pc)->mat, "_HPDDM_MatProduct", nullptr)); 796 PetscCall(MatDestroy(ctx->V)); 797 PetscCall(MatDestroy(ctx->V + 1)); 798 PetscCall(MatDestroy(ctx->V + 2)); 799 PetscCall(VecDestroy(&ctx->D)); 800 PetscCall(VecScatterDestroy(&ctx->scatter)); 801 PetscCall(PCDestroy(&ctx->pc)); 802 PetscFunctionReturn(PETSC_SUCCESS); 803 } 804 805 static PetscErrorCode PCHPDDMSolve_Private(const PC_HPDDM_Level *ctx, PetscScalar *rhs, const unsigned short &mu) 806 { 807 Mat B, X; 808 PetscInt n, N, j = 0; 809 810 PetscFunctionBegin; 811 PetscCall(KSPGetOperators(ctx->ksp, &B, nullptr)); 812 PetscCall(MatGetLocalSize(B, &n, nullptr)); 813 PetscCall(MatGetSize(B, &N, nullptr)); 814 if (ctx->parent->log_separate) { 815 j = std::distance(ctx->parent->levels, std::find(ctx->parent->levels, ctx->parent->levels + ctx->parent->N, ctx)); 816 PetscCall(PetscLogEventBegin(PC_HPDDM_Solve[j], ctx->ksp, nullptr, nullptr, nullptr)); 817 } 818 if (mu == 1) { 819 if (!ctx->ksp->vec_rhs) { 820 PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)ctx->ksp), 1, n, N, nullptr, &ctx->ksp->vec_rhs)); 821 PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)ctx->ksp), n, N, &ctx->ksp->vec_sol)); 822 } 823 PetscCall(VecPlaceArray(ctx->ksp->vec_rhs, rhs)); 824 PetscCall(KSPSolve(ctx->ksp, nullptr, nullptr)); 825 PetscCall(VecCopy(ctx->ksp->vec_sol, ctx->ksp->vec_rhs)); 826 PetscCall(VecResetArray(ctx->ksp->vec_rhs)); 827 } else { 828 PetscCall(MatCreateDense(PetscObjectComm((PetscObject)ctx->ksp), n, PETSC_DECIDE, N, mu, rhs, &B)); 829 PetscCall(MatCreateDense(PetscObjectComm((PetscObject)ctx->ksp), n, PETSC_DECIDE, N, mu, nullptr, &X)); 830 PetscCall(KSPMatSolve(ctx->ksp, B, X)); 831 PetscCall(MatCopy(X, B, SAME_NONZERO_PATTERN)); 832 PetscCall(MatDestroy(&X)); 833 PetscCall(MatDestroy(&B)); 834 } 835 if (ctx->parent->log_separate) PetscCall(PetscLogEventEnd(PC_HPDDM_Solve[j], ctx->ksp, nullptr, nullptr, nullptr)); 836 PetscFunctionReturn(PETSC_SUCCESS); 837 } 838 839 static PetscErrorCode PCHPDDMSetUpNeumannOverlap_Private(PC pc) 840 { 841 PC_HPDDM *data = (PC_HPDDM *)pc->data; 842 843 PetscFunctionBegin; 844 if (data->setup) { 845 Mat P; 846 Vec x, xt = nullptr; 847 PetscReal t = 0.0, s = 0.0; 848 849 PetscCall(PCGetOperators(pc, nullptr, &P)); 850 PetscCall(PetscObjectQuery((PetscObject)P, "__SNES_latest_X", (PetscObject *)&x)); 851 PetscCallBack("PCHPDDM Neumann callback", (*data->setup)(data->aux, t, x, xt, s, data->is, data->setup_ctx)); 852 } 853 PetscFunctionReturn(PETSC_SUCCESS); 854 } 855 856 static PetscErrorCode PCHPDDMCreateSubMatrices_Private(Mat mat, PetscInt n, const IS *, const IS *, MatReuse scall, Mat *submat[]) 857 { 858 Mat A; 859 860 PetscFunctionBegin; 861 PetscCheck(n == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "MatCreateSubMatrices() called to extract %" PetscInt_FMT " submatrices, which is different than 1", n); 862 /* previously composed Mat */ 863 PetscCall(PetscObjectQuery((PetscObject)mat, "_PCHPDDM_SubMatrices", (PetscObject *)&A)); 864 PetscCheck(A, PETSC_COMM_SELF, PETSC_ERR_PLIB, "SubMatrices not found in Mat"); 865 if (scall == MAT_INITIAL_MATRIX) { 866 PetscCall(PetscCalloc1(2, submat)); /* allocate an extra Mat to avoid errors in MatDestroySubMatrices_Dummy() */ 867 PetscCall(MatDuplicate(A, MAT_COPY_VALUES, *submat)); 868 } else PetscCall(MatCopy(A, (*submat)[0], SAME_NONZERO_PATTERN)); 869 PetscFunctionReturn(PETSC_SUCCESS); 870 } 871 872 static PetscErrorCode PCHPDDMCommunicationAvoidingPCASM_Private(PC pc, Mat C, PetscBool sorted) 873 { 874 void (*op)(void); 875 876 PetscFunctionBegin; 877 /* previously-composed Mat */ 878 PetscCall(PetscObjectCompose((PetscObject)pc->pmat, "_PCHPDDM_SubMatrices", (PetscObject)C)); 879 PetscCall(MatGetOperation(pc->pmat, MATOP_CREATE_SUBMATRICES, &op)); 880 /* trick suggested by Barry https://lists.mcs.anl.gov/pipermail/petsc-dev/2020-January/025491.html */ 881 PetscCall(MatSetOperation(pc->pmat, MATOP_CREATE_SUBMATRICES, (void (*)(void))PCHPDDMCreateSubMatrices_Private)); 882 if (sorted) PetscCall(PCASMSetSortIndices(pc, PETSC_FALSE)); /* everything is already sorted */ 883 PetscCall(PCSetFromOptions(pc)); /* otherwise -pc_hpddm_levels_1_pc_asm_sub_mat_type is not used */ 884 PetscCall(PCSetUp(pc)); 885 /* reset MatCreateSubMatrices() */ 886 PetscCall(MatSetOperation(pc->pmat, MATOP_CREATE_SUBMATRICES, op)); 887 PetscCall(PetscObjectCompose((PetscObject)pc->pmat, "_PCHPDDM_SubMatrices", nullptr)); 888 PetscFunctionReturn(PETSC_SUCCESS); 889 } 890 891 static PetscErrorCode PCHPDDMPermute_Private(IS is, IS in_is, IS *out_is, Mat in_C, Mat *out_C, IS *p) 892 { 893 IS perm; 894 const PetscInt *ptr; 895 PetscInt *concatenate, size, n, bs; 896 std::map<PetscInt, PetscInt> order; 897 PetscBool sorted; 898 899 PetscFunctionBegin; 900 PetscCall(ISSorted(is, &sorted)); 901 if (!sorted) { 902 PetscCall(ISGetLocalSize(is, &size)); 903 PetscCall(ISGetIndices(is, &ptr)); 904 PetscCall(ISGetBlockSize(is, &bs)); 905 /* MatCreateSubMatrices(), called by PCASM, follows the global numbering of Pmat */ 906 for (n = 0; n < size; n += bs) order.insert(std::make_pair(ptr[n] / bs, n / bs)); 907 PetscCall(ISRestoreIndices(is, &ptr)); 908 size /= bs; 909 if (out_C) { 910 PetscCall(PetscMalloc1(size, &concatenate)); 911 for (const std::pair<const PetscInt, PetscInt> &i : order) *concatenate++ = i.second; 912 concatenate -= size; 913 PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)in_C), bs, size, concatenate, PETSC_OWN_POINTER, &perm)); 914 PetscCall(ISSetPermutation(perm)); 915 /* permute user-provided Mat so that it matches with MatCreateSubMatrices() numbering */ 916 PetscCall(MatPermute(in_C, perm, perm, out_C)); 917 if (p) *p = perm; 918 else PetscCall(ISDestroy(&perm)); /* no need to save the permutation */ 919 } 920 if (out_is) { 921 PetscCall(PetscMalloc1(size, &concatenate)); 922 for (const std::pair<const PetscInt, PetscInt> &i : order) *concatenate++ = i.first; 923 concatenate -= size; 924 /* permute user-provided IS so that it matches with MatCreateSubMatrices() numbering */ 925 PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)in_is), bs, size, concatenate, PETSC_OWN_POINTER, out_is)); 926 } 927 } else { /* input IS is sorted, nothing to permute, simply duplicate inputs when needed */ 928 if (out_C) PetscCall(MatDuplicate(in_C, MAT_COPY_VALUES, out_C)); 929 if (out_is) PetscCall(ISDuplicate(in_is, out_is)); 930 } 931 PetscFunctionReturn(PETSC_SUCCESS); 932 } 933 934 static PetscErrorCode PCHPDDMDestroySubMatrices_Private(PetscBool flg, PetscBool algebraic, Mat *sub) 935 { 936 IS is; 937 938 PetscFunctionBegin; 939 if (!flg) { 940 if (algebraic) { 941 PetscCall(PetscObjectQuery((PetscObject)sub[0], "_PCHPDDM_Embed", (PetscObject *)&is)); 942 PetscCall(ISDestroy(&is)); 943 PetscCall(PetscObjectCompose((PetscObject)sub[0], "_PCHPDDM_Embed", nullptr)); 944 PetscCall(PetscObjectCompose((PetscObject)sub[0], "_PCHPDDM_Compact", nullptr)); 945 } 946 PetscCall(MatDestroySubMatrices(algebraic ? 2 : 1, &sub)); 947 } 948 PetscFunctionReturn(PETSC_SUCCESS); 949 } 950 951 static PetscErrorCode PCHPDDMAlgebraicAuxiliaryMat_Private(Mat P, IS *is, Mat *sub[], PetscBool block) 952 { 953 IS icol[3], irow[2]; 954 Mat *M, Q; 955 PetscReal *ptr; 956 PetscInt *idx, p = 0, n, bs = PetscAbs(P->cmap->bs); 957 PetscBool flg; 958 959 PetscFunctionBegin; 960 PetscCall(ISCreateStride(PETSC_COMM_SELF, P->cmap->N, 0, 1, icol + 2)); 961 PetscCall(ISSetBlockSize(icol[2], bs)); 962 PetscCall(ISSetIdentity(icol[2])); 963 PetscCall(PetscObjectTypeCompare((PetscObject)P, MATMPISBAIJ, &flg)); 964 if (flg) { 965 /* MatCreateSubMatrices() does not handle MATMPISBAIJ properly when iscol != isrow, so convert first to MATMPIBAIJ */ 966 PetscCall(MatConvert(P, MATMPIBAIJ, MAT_INITIAL_MATRIX, &Q)); 967 std::swap(P, Q); 968 } 969 PetscCall(MatCreateSubMatrices(P, 1, is, icol + 2, MAT_INITIAL_MATRIX, &M)); 970 if (flg) { 971 std::swap(P, Q); 972 PetscCall(MatDestroy(&Q)); 973 } 974 PetscCall(ISDestroy(icol + 2)); 975 PetscCall(ISCreateStride(PETSC_COMM_SELF, M[0]->rmap->N, 0, 1, irow)); 976 PetscCall(ISSetBlockSize(irow[0], bs)); 977 PetscCall(ISSetIdentity(irow[0])); 978 if (!block) { 979 PetscCall(PetscMalloc2(P->cmap->N, &ptr, P->cmap->N / bs, &idx)); 980 PetscCall(MatGetColumnNorms(M[0], NORM_INFINITY, ptr)); 981 /* check for nonzero columns so that M[0] may be expressed in compact form */ 982 for (n = 0; n < P->cmap->N; n += bs) { 983 if (std::find_if(ptr + n, ptr + n + bs, [](PetscReal v) { return v > PETSC_MACHINE_EPSILON; }) != ptr + n + bs) idx[p++] = n / bs; 984 } 985 PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, p, idx, PETSC_USE_POINTER, icol + 1)); 986 PetscCall(ISSetInfo(icol[1], IS_SORTED, IS_GLOBAL, PETSC_TRUE, PETSC_TRUE)); 987 PetscCall(ISEmbed(*is, icol[1], PETSC_FALSE, icol + 2)); 988 irow[1] = irow[0]; 989 /* 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 */ 990 icol[0] = is[0]; 991 PetscCall(MatCreateSubMatrices(M[0], 2, irow, icol, MAT_INITIAL_MATRIX, sub)); 992 PetscCall(ISDestroy(icol + 1)); 993 PetscCall(PetscFree2(ptr, idx)); 994 /* IS used to go back and forth between the augmented and the original local linear system, see eq. (3.4) of [2022b] */ 995 PetscCall(PetscObjectCompose((PetscObject)(*sub)[0], "_PCHPDDM_Embed", (PetscObject)icol[2])); 996 /* Mat used in eq. (3.1) of [2022b] */ 997 PetscCall(PetscObjectCompose((PetscObject)(*sub)[0], "_PCHPDDM_Compact", (PetscObject)(*sub)[1])); 998 } else { 999 Mat aux; 1000 PetscCall(MatSetOption(M[0], MAT_SUBMAT_SINGLEIS, PETSC_TRUE)); 1001 /* diagonal block of the overlapping rows */ 1002 PetscCall(MatCreateSubMatrices(M[0], 1, irow, is, MAT_INITIAL_MATRIX, sub)); 1003 PetscCall(MatDuplicate((*sub)[0], MAT_COPY_VALUES, &aux)); 1004 PetscCall(MatSetOption(aux, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE)); 1005 if (bs == 1) { /* scalar case */ 1006 Vec sum[2]; 1007 PetscCall(MatCreateVecs(aux, sum, sum + 1)); 1008 PetscCall(MatGetRowSum(M[0], sum[0])); 1009 PetscCall(MatGetRowSum(aux, sum[1])); 1010 /* off-diagonal block row sum (full rows - diagonal block rows) */ 1011 PetscCall(VecAXPY(sum[0], -1.0, sum[1])); 1012 /* subdomain matrix plus off-diagonal block row sum */ 1013 PetscCall(MatDiagonalSet(aux, sum[0], ADD_VALUES)); 1014 PetscCall(VecDestroy(sum)); 1015 PetscCall(VecDestroy(sum + 1)); 1016 } else { /* vectorial case */ 1017 /* TODO: missing MatGetValuesBlocked(), so the code below is */ 1018 /* an extension of the scalar case for when bs > 1, but it could */ 1019 /* be more efficient by avoiding all these MatMatMult() */ 1020 Mat sum[2], ones; 1021 PetscScalar *ptr; 1022 PetscCall(PetscCalloc1(M[0]->cmap->n * bs, &ptr)); 1023 PetscCall(MatCreateDense(PETSC_COMM_SELF, M[0]->cmap->n, bs, M[0]->cmap->n, bs, ptr, &ones)); 1024 for (n = 0; n < M[0]->cmap->n; n += bs) { 1025 for (p = 0; p < bs; ++p) ptr[n + p * (M[0]->cmap->n + 1)] = 1.0; 1026 } 1027 PetscCall(MatMatMult(M[0], ones, MAT_INITIAL_MATRIX, PETSC_DEFAULT, sum)); 1028 PetscCall(MatDestroy(&ones)); 1029 PetscCall(MatCreateDense(PETSC_COMM_SELF, aux->cmap->n, bs, aux->cmap->n, bs, ptr, &ones)); 1030 PetscCall(MatDenseSetLDA(ones, M[0]->cmap->n)); 1031 PetscCall(MatMatMult(aux, ones, MAT_INITIAL_MATRIX, PETSC_DEFAULT, sum + 1)); 1032 PetscCall(MatDestroy(&ones)); 1033 PetscCall(PetscFree(ptr)); 1034 /* off-diagonal block row sum (full rows - diagonal block rows) */ 1035 PetscCall(MatAXPY(sum[0], -1.0, sum[1], SAME_NONZERO_PATTERN)); 1036 PetscCall(MatDestroy(sum + 1)); 1037 /* re-order values to be consistent with MatSetValuesBlocked() */ 1038 /* equivalent to MatTranspose() which does not truly handle */ 1039 /* MAT_INPLACE_MATRIX in the rectangular case, as it calls PetscMalloc() */ 1040 PetscCall(MatDenseGetArrayWrite(sum[0], &ptr)); 1041 HPDDM::Wrapper<PetscScalar>::imatcopy<'T'>(bs, sum[0]->rmap->n, ptr, sum[0]->rmap->n, bs); 1042 /* subdomain matrix plus off-diagonal block row sum */ 1043 for (n = 0; n < aux->cmap->n / bs; ++n) PetscCall(MatSetValuesBlocked(aux, 1, &n, 1, &n, ptr + n * bs * bs, ADD_VALUES)); 1044 PetscCall(MatAssemblyBegin(aux, MAT_FINAL_ASSEMBLY)); 1045 PetscCall(MatAssemblyEnd(aux, MAT_FINAL_ASSEMBLY)); 1046 PetscCall(MatDenseRestoreArrayWrite(sum[0], &ptr)); 1047 PetscCall(MatDestroy(sum)); 1048 } 1049 PetscCall(MatSetOption(aux, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE)); 1050 /* left-hand side of GenEO, with the same sparsity pattern as PCASM subdomain solvers */ 1051 PetscCall(PetscObjectCompose((PetscObject)(*sub)[0], "_PCHPDDM_Neumann_Mat", (PetscObject)aux)); 1052 } 1053 PetscCall(ISDestroy(irow)); 1054 PetscCall(MatDestroySubMatrices(1, &M)); 1055 PetscFunctionReturn(PETSC_SUCCESS); 1056 } 1057 1058 static PetscErrorCode PCSetUp_HPDDM(PC pc) 1059 { 1060 PC_HPDDM *data = (PC_HPDDM *)pc->data; 1061 PC inner; 1062 KSP *ksp; 1063 Mat *sub, A, P, N, C = nullptr, uaux = nullptr, weighted, subA[2]; 1064 Vec xin, v; 1065 std::vector<Vec> initial; 1066 IS is[1], loc, uis = data->is, unsorted = nullptr; 1067 ISLocalToGlobalMapping l2g; 1068 char prefix[256]; 1069 const char *pcpre; 1070 const PetscScalar *const *ev; 1071 PetscInt n, requested = data->N, reused = 0; 1072 MatStructure structure = UNKNOWN_NONZERO_PATTERN; 1073 PetscBool subdomains = PETSC_FALSE, flg = PETSC_FALSE, ismatis, swap = PETSC_FALSE, algebraic = PETSC_FALSE, block = PETSC_FALSE; 1074 DM dm; 1075 #if PetscDefined(USE_DEBUG) 1076 IS dis = nullptr; 1077 Mat daux = nullptr; 1078 #endif 1079 1080 PetscFunctionBegin; 1081 PetscCheck(data->levels && data->levels[0], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Not a single level allocated"); 1082 PetscCall(PCGetOptionsPrefix(pc, &pcpre)); 1083 PetscCall(PCGetOperators(pc, &A, &P)); 1084 if (!data->levels[0]->ksp) { 1085 PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &data->levels[0]->ksp)); 1086 PetscCall(KSPSetNestLevel(data->levels[0]->ksp, pc->kspnestlevel)); 1087 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_%s_", pcpre ? pcpre : "", data->N > 1 ? "levels_1" : "coarse")); 1088 PetscCall(KSPSetOptionsPrefix(data->levels[0]->ksp, prefix)); 1089 PetscCall(KSPSetType(data->levels[0]->ksp, KSPPREONLY)); 1090 } else if (data->levels[0]->ksp->pc && data->levels[0]->ksp->pc->setupcalled == 1 && data->levels[0]->ksp->pc->reusepreconditioner) { 1091 /* if the fine-level PCSHELL exists, its setup has succeeded, and one wants to reuse it, */ 1092 /* then just propagate the appropriate flag to the coarser levels */ 1093 for (n = 0; n < PETSC_PCHPDDM_MAXLEVELS && data->levels[n]; ++n) { 1094 /* the following KSP and PC may be NULL for some processes, hence the check */ 1095 if (data->levels[n]->ksp) PetscCall(KSPSetReusePreconditioner(data->levels[n]->ksp, PETSC_TRUE)); 1096 if (data->levels[n]->pc) PetscCall(PCSetReusePreconditioner(data->levels[n]->pc, PETSC_TRUE)); 1097 } 1098 /* early bail out because there is nothing to do */ 1099 PetscFunctionReturn(PETSC_SUCCESS); 1100 } else { 1101 /* reset coarser levels */ 1102 for (n = 1; n < PETSC_PCHPDDM_MAXLEVELS && data->levels[n]; ++n) { 1103 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) { 1104 reused = data->N - n; 1105 break; 1106 } 1107 PetscCall(KSPDestroy(&data->levels[n]->ksp)); 1108 PetscCall(PCDestroy(&data->levels[n]->pc)); 1109 } 1110 /* check if some coarser levels are being reused */ 1111 PetscCall(MPIU_Allreduce(MPI_IN_PLACE, &reused, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)pc))); 1112 const int *addr = data->levels[0]->P ? data->levels[0]->P->getAddrLocal() : &HPDDM::i__0; 1113 1114 if (addr != &HPDDM::i__0 && reused != data->N - 1) { 1115 /* reuse previously computed eigenvectors */ 1116 ev = data->levels[0]->P->getVectors(); 1117 if (ev) { 1118 initial.reserve(*addr); 1119 PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, data->levels[0]->P->getDof(), ev[0], &xin)); 1120 for (n = 0; n < *addr; ++n) { 1121 PetscCall(VecDuplicate(xin, &v)); 1122 PetscCall(VecPlaceArray(xin, ev[n])); 1123 PetscCall(VecCopy(xin, v)); 1124 initial.emplace_back(v); 1125 PetscCall(VecResetArray(xin)); 1126 } 1127 PetscCall(VecDestroy(&xin)); 1128 } 1129 } 1130 } 1131 data->N -= reused; 1132 PetscCall(KSPSetOperators(data->levels[0]->ksp, A, P)); 1133 1134 PetscCall(PetscObjectTypeCompare((PetscObject)P, MATIS, &ismatis)); 1135 if (!data->is && !ismatis) { 1136 PetscErrorCode (*create)(DM, IS *, Mat *, PetscErrorCode (**)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void *), void **) = nullptr; 1137 PetscErrorCode (*usetup)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void *) = nullptr; 1138 void *uctx = nullptr; 1139 1140 /* first see if we can get the data from the DM */ 1141 PetscCall(MatGetDM(P, &dm)); 1142 if (!dm) PetscCall(MatGetDM(A, &dm)); 1143 if (!dm) PetscCall(PCGetDM(pc, &dm)); 1144 if (dm) { /* this is the hook for DMPLEX and DMDA for which the auxiliary Mat is the local Neumann matrix */ 1145 PetscCall(PetscObjectQueryFunction((PetscObject)dm, "DMCreateNeumannOverlap_C", &create)); 1146 if (create) { 1147 PetscCall((*create)(dm, &uis, &uaux, &usetup, &uctx)); 1148 if (data->Neumann == PETSC_BOOL3_UNKNOWN) data->Neumann = PETSC_BOOL3_TRUE; /* set the value only if it was not already provided by the user */ 1149 } 1150 } 1151 if (!create) { 1152 if (!uis) { 1153 PetscCall(PetscObjectQuery((PetscObject)pc, "_PCHPDDM_Neumann_IS", (PetscObject *)&uis)); 1154 PetscCall(PetscObjectReference((PetscObject)uis)); 1155 } 1156 if (!uaux) { 1157 PetscCall(PetscObjectQuery((PetscObject)pc, "_PCHPDDM_Neumann_Mat", (PetscObject *)&uaux)); 1158 PetscCall(PetscObjectReference((PetscObject)uaux)); 1159 } 1160 /* look inside the Pmat instead of the PC, needed for MatSchurComplementComputeExplicitOperator() */ 1161 if (!uis) { 1162 PetscCall(PetscObjectQuery((PetscObject)P, "_PCHPDDM_Neumann_IS", (PetscObject *)&uis)); 1163 PetscCall(PetscObjectReference((PetscObject)uis)); 1164 } 1165 if (!uaux) { 1166 PetscCall(PetscObjectQuery((PetscObject)P, "_PCHPDDM_Neumann_Mat", (PetscObject *)&uaux)); 1167 PetscCall(PetscObjectReference((PetscObject)uaux)); 1168 } 1169 } 1170 PetscCall(PCHPDDMSetAuxiliaryMat(pc, uis, uaux, usetup, uctx)); 1171 PetscCall(MatDestroy(&uaux)); 1172 PetscCall(ISDestroy(&uis)); 1173 } 1174 1175 if (!ismatis) { 1176 PetscCall(PCHPDDMSetUpNeumannOverlap_Private(pc)); 1177 PetscCall(PetscOptionsGetBool(nullptr, pcpre, "-pc_hpddm_block_splitting", &block, nullptr)); 1178 if (data->is && block) { 1179 PetscCall(ISDestroy(&data->is)); 1180 PetscCall(MatDestroy(&data->aux)); 1181 } 1182 if (!data->is && data->N > 1) { 1183 char type[256] = {}; /* same size as in src/ksp/pc/interface/pcset.c */ 1184 PetscCall(PetscObjectTypeCompareAny((PetscObject)P, &flg, MATNORMAL, MATNORMALHERMITIAN, "")); 1185 if (flg || (A->rmap->N != A->cmap->N && P->rmap->N == P->cmap->N && P->rmap->N == A->cmap->N)) { 1186 Mat B; 1187 PetscCall(PCHPDDMSetAuxiliaryMatNormal_Private(pc, A, P, &B, pcpre)); 1188 if (data->correction == PC_HPDDM_COARSE_CORRECTION_DEFLATED) data->correction = PC_HPDDM_COARSE_CORRECTION_BALANCED; 1189 PetscCall(MatDestroy(&B)); 1190 } else { 1191 PetscCall(PetscObjectTypeCompare((PetscObject)P, MATSCHURCOMPLEMENT, &flg)); 1192 if (flg) { 1193 Mat A00, P00, A01, A10, A11, B, N; 1194 Vec diagonal = nullptr; 1195 const PetscScalar *array; 1196 MatSchurComplementAinvType type; 1197 1198 PetscCall(MatSchurComplementGetSubMatrices(P, &A00, &P00, &A01, &A10, &A11)); 1199 if (A11) { 1200 PetscCall(MatCreateVecs(A11, &diagonal, nullptr)); 1201 PetscCall(MatGetDiagonal(A11, diagonal)); 1202 } 1203 if (PetscDefined(USE_DEBUG)) { 1204 Mat T, U = nullptr; 1205 IS z; 1206 PetscCall(PetscObjectTypeCompare((PetscObject)A10, MATTRANSPOSEVIRTUAL, &flg)); 1207 if (flg) PetscCall(MatTransposeGetMat(A10, &U)); 1208 else { 1209 PetscCall(PetscObjectTypeCompare((PetscObject)A10, MATHERMITIANTRANSPOSEVIRTUAL, &flg)); 1210 if (flg) PetscCall(MatHermitianTransposeGetMat(A10, &U)); 1211 } 1212 if (U) PetscCall(MatDuplicate(U, MAT_COPY_VALUES, &T)); 1213 else PetscCall(MatHermitianTranspose(A10, MAT_INITIAL_MATRIX, &T)); 1214 PetscCall(PetscLayoutCompare(T->rmap, A01->rmap, &flg)); 1215 if (flg) { 1216 PetscCall(PetscLayoutCompare(T->cmap, A01->cmap, &flg)); 1217 if (flg) { 1218 PetscCall(MatFindZeroRows(A01, &z)); /* for essential boundary conditions, some implementations will */ 1219 if (z) { /* zero rows in [P00 A01] except for the diagonal of P00 */ 1220 PetscCall(MatSetOption(T, MAT_NO_OFF_PROC_ZERO_ROWS, PETSC_TRUE)); 1221 PetscCall(MatZeroRowsIS(T, z, 0.0, nullptr, nullptr)); /* corresponding zero rows from A01 */ 1222 PetscCall(ISDestroy(&z)); 1223 } 1224 PetscCall(MatMultEqual(A01, T, 10, &flg)); 1225 PetscCheck(flg, PetscObjectComm((PetscObject)P), PETSC_ERR_SUP, "A01 != A10^T"); 1226 } else PetscCall(PetscInfo(pc, "A01 and A10^T have non-congruent column layouts, cannot test for equality\n")); 1227 } 1228 PetscCall(MatDestroy(&T)); 1229 } 1230 PetscCall(MatCreateVecs(P00, &v, nullptr)); 1231 PetscCall(MatSchurComplementGetAinvType(P, &type)); 1232 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]); 1233 if (type == MAT_SCHUR_COMPLEMENT_AINV_LUMP) { 1234 PetscCall(MatGetRowSum(P00, v)); 1235 if (A00 == P00) PetscCall(PetscObjectReference((PetscObject)A00)); 1236 PetscCall(MatDestroy(&P00)); 1237 PetscCall(VecGetArrayRead(v, &array)); 1238 PetscCall(MatCreateAIJ(PetscObjectComm((PetscObject)A00), A00->rmap->n, A00->cmap->n, A00->rmap->N, A00->cmap->N, 1, nullptr, 0, nullptr, &P00)); 1239 PetscCall(MatSetOption(P00, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE)); 1240 for (n = A00->rmap->rstart; n < A00->rmap->rend; ++n) PetscCall(MatSetValue(P00, n, n, array[n - A00->rmap->rstart], INSERT_VALUES)); 1241 PetscCall(MatAssemblyBegin(P00, MAT_FINAL_ASSEMBLY)); 1242 PetscCall(MatAssemblyEnd(P00, MAT_FINAL_ASSEMBLY)); 1243 PetscCall(VecRestoreArrayRead(v, &array)); 1244 PetscCall(MatSchurComplementUpdateSubMatrices(P, A00, P00, A01, A10, A11)); /* replace P00 by diag(sum of each row of P00) */ 1245 PetscCall(MatDestroy(&P00)); 1246 } else PetscCall(MatGetDiagonal(P00, v)); 1247 PetscCall(VecReciprocal(v)); /* inv(diag(P00)) */ 1248 PetscCall(VecSqrtAbs(v)); /* sqrt(inv(diag(P00))) */ 1249 PetscCall(MatDuplicate(A01, MAT_COPY_VALUES, &B)); 1250 PetscCall(MatDiagonalScale(B, v, nullptr)); 1251 PetscCall(VecDestroy(&v)); 1252 PetscCall(MatCreateNormalHermitian(B, &N)); 1253 PetscCall(PCHPDDMSetAuxiliaryMatNormal_Private(pc, B, N, &P, pcpre, &diagonal)); 1254 PetscCall(PetscObjectTypeCompare((PetscObject)data->aux, MATSEQAIJ, &flg)); 1255 if (!flg) { 1256 PetscCall(MatDestroy(&P)); 1257 P = N; 1258 PetscCall(PetscObjectReference((PetscObject)P)); 1259 } else PetscCall(MatScale(P, -1.0)); 1260 if (diagonal) { 1261 PetscCall(MatDiagonalSet(P, diagonal, ADD_VALUES)); 1262 PetscCall(PCSetOperators(pc, P, P)); /* replace P by diag(P11) - A01' inv(diag(P00)) A01 */ 1263 PetscCall(VecDestroy(&diagonal)); 1264 } else { 1265 PetscCall(MatScale(N, -1.0)); 1266 PetscCall(PCSetOperators(pc, N, P)); /* replace P by - A01' inv(diag(P00)) A01 */ 1267 } 1268 PetscCall(MatDestroy(&N)); 1269 PetscCall(MatDestroy(&P)); 1270 PetscCall(MatDestroy(&B)); 1271 PetscFunctionReturn(PETSC_SUCCESS); 1272 } else { 1273 PetscCall(PetscOptionsGetString(nullptr, pcpre, "-pc_hpddm_levels_1_st_pc_type", type, sizeof(type), nullptr)); 1274 PetscCall(PetscStrcmp(type, PCMAT, &algebraic)); 1275 PetscCall(PetscOptionsGetBool(nullptr, pcpre, "-pc_hpddm_block_splitting", &block, nullptr)); 1276 PetscCheck(!algebraic || !block, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-pc_hpddm_levels_1_st_pc_type mat and -pc_hpddm_block_splitting"); 1277 if (block) algebraic = PETSC_TRUE; 1278 if (algebraic) { 1279 PetscCall(ISCreateStride(PETSC_COMM_SELF, P->rmap->n, P->rmap->rstart, 1, &data->is)); 1280 PetscCall(MatIncreaseOverlap(P, 1, &data->is, 1)); 1281 PetscCall(ISSort(data->is)); 1282 } 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 : "")); 1283 } 1284 } 1285 } 1286 } 1287 #if PetscDefined(USE_DEBUG) 1288 if (data->is) PetscCall(ISDuplicate(data->is, &dis)); 1289 if (data->aux) PetscCall(MatDuplicate(data->aux, MAT_COPY_VALUES, &daux)); 1290 #endif 1291 if (data->is || (ismatis && data->N > 1)) { 1292 if (ismatis) { 1293 std::initializer_list<std::string> list = {MATSEQBAIJ, MATSEQSBAIJ}; 1294 PetscCall(MatISGetLocalMat(P, &N)); 1295 std::initializer_list<std::string>::const_iterator it = std::find(list.begin(), list.end(), ((PetscObject)N)->type_name); 1296 PetscCall(MatISRestoreLocalMat(P, &N)); 1297 switch (std::distance(list.begin(), it)) { 1298 case 0: 1299 PetscCall(MatConvert(P, MATMPIBAIJ, MAT_INITIAL_MATRIX, &C)); 1300 break; 1301 case 1: 1302 /* MatCreateSubMatrices() does not work with MATSBAIJ and unsorted ISes, so convert to MPIBAIJ */ 1303 PetscCall(MatConvert(P, MATMPIBAIJ, MAT_INITIAL_MATRIX, &C)); 1304 PetscCall(MatSetOption(C, MAT_SYMMETRIC, PETSC_TRUE)); 1305 break; 1306 default: 1307 PetscCall(MatConvert(P, MATMPIAIJ, MAT_INITIAL_MATRIX, &C)); 1308 } 1309 PetscCall(MatISGetLocalToGlobalMapping(P, &l2g, nullptr)); 1310 PetscCall(PetscObjectReference((PetscObject)P)); 1311 PetscCall(KSPSetOperators(data->levels[0]->ksp, A, C)); 1312 std::swap(C, P); 1313 PetscCall(ISLocalToGlobalMappingGetSize(l2g, &n)); 1314 PetscCall(ISCreateStride(PETSC_COMM_SELF, n, 0, 1, &loc)); 1315 PetscCall(ISLocalToGlobalMappingApplyIS(l2g, loc, &is[0])); 1316 PetscCall(ISDestroy(&loc)); 1317 /* the auxiliary Mat is _not_ the local Neumann matrix */ 1318 /* it is the local Neumann matrix augmented (with zeros) through MatIncreaseOverlap() */ 1319 data->Neumann = PETSC_BOOL3_FALSE; 1320 structure = SAME_NONZERO_PATTERN; 1321 } else { 1322 is[0] = data->is; 1323 if (algebraic) subdomains = PETSC_TRUE; 1324 PetscCall(PetscOptionsGetBool(nullptr, pcpre, "-pc_hpddm_define_subdomains", &subdomains, nullptr)); 1325 if (PetscBool3ToBool(data->Neumann)) { 1326 PetscCheck(!block, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-pc_hpddm_block_splitting and -pc_hpddm_has_neumann"); 1327 PetscCheck(!algebraic, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-pc_hpddm_levels_1_st_pc_type mat and -pc_hpddm_has_neumann"); 1328 } 1329 if (PetscBool3ToBool(data->Neumann) || block) structure = SAME_NONZERO_PATTERN; 1330 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)data->is), P->rmap->n, P->rmap->rstart, 1, &loc)); 1331 } 1332 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_levels_1_", pcpre ? pcpre : "")); 1333 PetscCall(PetscOptionsGetEnum(nullptr, prefix, "-st_matstructure", MatStructures, (PetscEnum *)&structure, &flg)); /* if not user-provided, force its value when possible */ 1334 if (!flg && structure == SAME_NONZERO_PATTERN) { /* cannot call STSetMatStructure() yet, insert the appropriate option in the database, parsed by STSetFromOptions() */ 1335 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-%spc_hpddm_levels_1_st_matstructure", pcpre ? pcpre : "")); 1336 PetscCall(PetscOptionsSetValue(nullptr, prefix, MatStructures[structure])); 1337 } 1338 flg = PETSC_FALSE; 1339 if (data->share) { 1340 data->share = PETSC_FALSE; /* will be reset to PETSC_TRUE if none of the conditions below are true */ 1341 if (!subdomains) PetscCall(PetscInfo(pc, "Cannot share subdomain KSP between SLEPc and PETSc since -%spc_hpddm_define_subdomains is not true\n", pcpre ? pcpre : "")); 1342 else if (data->deflation) PetscCall(PetscInfo(pc, "Nothing to share since PCHPDDMSetDeflationMat() has been called\n")); 1343 else if (ismatis) PetscCall(PetscInfo(pc, "Cannot share subdomain KSP between SLEPc and PETSc with a Pmat of type MATIS\n")); 1344 else if (!algebraic && structure != SAME_NONZERO_PATTERN) 1345 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])); 1346 else data->share = PETSC_TRUE; 1347 } 1348 if (!ismatis) { 1349 if (data->share || (!PetscBool3ToBool(data->Neumann) && subdomains)) PetscCall(ISDuplicate(is[0], &unsorted)); 1350 else unsorted = is[0]; 1351 } 1352 if (data->N > 1 && (data->aux || ismatis || algebraic)) { 1353 PetscCheck(loadedSym, PETSC_COMM_SELF, PETSC_ERR_PLIB, "HPDDM library not loaded, cannot use more than one level"); 1354 PetscCall(MatSetOption(P, MAT_SUBMAT_SINGLEIS, PETSC_TRUE)); 1355 if (ismatis) { 1356 /* needed by HPDDM (currently) so that the partition of unity is 0 on subdomain interfaces */ 1357 PetscCall(MatIncreaseOverlap(P, 1, is, 1)); 1358 PetscCall(ISDestroy(&data->is)); 1359 data->is = is[0]; 1360 } else { 1361 if (PetscDefined(USE_DEBUG)) { 1362 PetscBool equal; 1363 IS intersect; 1364 1365 PetscCall(ISIntersect(data->is, loc, &intersect)); 1366 PetscCall(ISEqualUnsorted(loc, intersect, &equal)); 1367 PetscCall(ISDestroy(&intersect)); 1368 PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "IS of the auxiliary Mat does not include all local rows of A"); 1369 } 1370 PetscCall(PetscObjectComposeFunction((PetscObject)pc->pmat, "PCHPDDMAlgebraicAuxiliaryMat_Private_C", PCHPDDMAlgebraicAuxiliaryMat_Private)); 1371 if (!PetscBool3ToBool(data->Neumann) && !algebraic) { 1372 PetscCall(PetscObjectTypeCompare((PetscObject)P, MATMPISBAIJ, &flg)); 1373 if (flg) { 1374 /* maybe better to ISSort(is[0]), MatCreateSubMatrices(), and then MatPermute() */ 1375 /* but there is no MatPermute_SeqSBAIJ(), so as before, just use MATMPIBAIJ */ 1376 PetscCall(MatConvert(P, MATMPIBAIJ, MAT_INITIAL_MATRIX, &uaux)); 1377 flg = PETSC_FALSE; 1378 } 1379 } 1380 } 1381 if (algebraic) { 1382 PetscUseMethod(pc->pmat, "PCHPDDMAlgebraicAuxiliaryMat_Private_C", (Mat, IS *, Mat *[], PetscBool), (P, is, &sub, block)); 1383 if (block) { 1384 PetscCall(PetscObjectQuery((PetscObject)sub[0], "_PCHPDDM_Neumann_Mat", (PetscObject *)&data->aux)); 1385 PetscCall(PetscObjectCompose((PetscObject)sub[0], "_PCHPDDM_Neumann_Mat", nullptr)); 1386 } 1387 } else if (!uaux) { 1388 if (PetscBool3ToBool(data->Neumann)) sub = &data->aux; 1389 else PetscCall(MatCreateSubMatrices(P, 1, is, is, MAT_INITIAL_MATRIX, &sub)); 1390 } else { 1391 PetscCall(MatCreateSubMatrices(uaux, 1, is, is, MAT_INITIAL_MATRIX, &sub)); 1392 PetscCall(MatDestroy(&uaux)); 1393 PetscCall(MatConvert(sub[0], MATSEQSBAIJ, MAT_INPLACE_MATRIX, sub)); 1394 } 1395 /* Vec holding the partition of unity */ 1396 if (!data->levels[0]->D) { 1397 PetscCall(ISGetLocalSize(data->is, &n)); 1398 PetscCall(VecCreateMPI(PETSC_COMM_SELF, n, PETSC_DETERMINE, &data->levels[0]->D)); 1399 } 1400 if (data->share) { 1401 Mat D; 1402 IS perm = nullptr; 1403 PetscInt size = -1; 1404 if (!data->levels[0]->pc) { 1405 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_levels_1_", pcpre ? pcpre : "")); 1406 PetscCall(PCCreate(PetscObjectComm((PetscObject)pc), &data->levels[0]->pc)); 1407 PetscCall(PCSetOptionsPrefix(data->levels[0]->pc, prefix)); 1408 PetscCall(PCSetOperators(data->levels[0]->pc, A, P)); 1409 } 1410 PetscCall(PCSetType(data->levels[0]->pc, PCASM)); 1411 if (!data->levels[0]->pc->setupcalled) { 1412 IS sorted; /* PCASM will sort the input IS, duplicate it to return an unmodified (PCHPDDM) input IS */ 1413 PetscCall(ISDuplicate(is[0], &sorted)); 1414 PetscCall(PCASMSetLocalSubdomains(data->levels[0]->pc, 1, &sorted, &loc)); 1415 PetscCall(PetscObjectDereference((PetscObject)sorted)); 1416 } 1417 PetscCall(PCSetFromOptions(data->levels[0]->pc)); 1418 if (block) { 1419 PetscCall(PCHPDDMPermute_Private(unsorted, data->is, &uis, sub[0], &C, &perm)); 1420 PetscCall(PCHPDDMCommunicationAvoidingPCASM_Private(data->levels[0]->pc, C, algebraic)); 1421 } else PetscCall(PCSetUp(data->levels[0]->pc)); 1422 PetscTryMethod(data->levels[0]->pc, "PCASMGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (data->levels[0]->pc, &size, nullptr, &ksp)); 1423 if (size != 1) { 1424 data->share = PETSC_FALSE; 1425 PetscCheck(size == -1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of subdomain solver %" PetscInt_FMT " != 1", size); 1426 PetscCall(PetscInfo(pc, "Cannot share subdomain KSP between SLEPc and PETSc since PCASMGetSubKSP() not found in fine-level PC\n")); 1427 PetscCall(ISDestroy(&unsorted)); 1428 unsorted = is[0]; 1429 } else { 1430 if (!block) PetscCall(PCHPDDMPermute_Private(unsorted, data->is, &uis, PetscBool3ToBool(data->Neumann) ? sub[0] : data->aux, &C, &perm)); 1431 if (!PetscBool3ToBool(data->Neumann) && !block) { 1432 PetscCall(MatPermute(sub[0], perm, perm, &D)); /* permute since PCASM will call ISSort() */ 1433 PetscCall(MatHeaderReplace(sub[0], &D)); 1434 } 1435 if (data->B) { /* see PCHPDDMSetRHSMat() */ 1436 PetscCall(MatPermute(data->B, perm, perm, &D)); 1437 PetscCall(MatHeaderReplace(data->B, &D)); 1438 } 1439 PetscCall(ISDestroy(&perm)); 1440 const char *matpre; 1441 PetscBool cmp[4]; 1442 PetscCall(KSPGetOperators(ksp[0], subA, subA + 1)); 1443 PetscCall(MatDuplicate(subA[1], MAT_SHARE_NONZERO_PATTERN, &D)); 1444 PetscCall(MatGetOptionsPrefix(subA[1], &matpre)); 1445 PetscCall(MatSetOptionsPrefix(D, matpre)); 1446 PetscCall(PetscObjectTypeCompare((PetscObject)D, MATNORMAL, cmp)); 1447 PetscCall(PetscObjectTypeCompare((PetscObject)C, MATNORMAL, cmp + 1)); 1448 if (!cmp[0]) PetscCall(PetscObjectTypeCompare((PetscObject)D, MATNORMALHERMITIAN, cmp + 2)); 1449 else cmp[2] = PETSC_FALSE; 1450 if (!cmp[1]) PetscCall(PetscObjectTypeCompare((PetscObject)C, MATNORMALHERMITIAN, cmp + 3)); 1451 else cmp[3] = PETSC_FALSE; 1452 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); 1453 if (!cmp[0] && !cmp[2]) { 1454 if (!block) PetscCall(MatAXPY(D, 1.0, C, SUBSET_NONZERO_PATTERN)); 1455 else { 1456 PetscCall(MatMissingDiagonal(D, cmp, nullptr)); 1457 if (cmp[0]) structure = DIFFERENT_NONZERO_PATTERN; /* data->aux has no missing diagonal entry */ 1458 PetscCall(MatAXPY(D, 1.0, data->aux, structure)); 1459 } 1460 } else { 1461 Mat mat[2]; 1462 if (cmp[0]) { 1463 PetscCall(MatNormalGetMat(D, mat)); 1464 PetscCall(MatNormalGetMat(C, mat + 1)); 1465 } else { 1466 PetscCall(MatNormalHermitianGetMat(D, mat)); 1467 PetscCall(MatNormalHermitianGetMat(C, mat + 1)); 1468 } 1469 PetscCall(MatAXPY(mat[0], 1.0, mat[1], SUBSET_NONZERO_PATTERN)); 1470 } 1471 PetscCall(MatPropagateSymmetryOptions(C, D)); 1472 PetscCall(MatDestroy(&C)); 1473 C = D; 1474 /* swap pointers so that variables stay consistent throughout PCSetUp() */ 1475 std::swap(C, data->aux); 1476 std::swap(uis, data->is); 1477 swap = PETSC_TRUE; 1478 } 1479 } 1480 if (!data->levels[0]->scatter) { 1481 PetscCall(MatCreateVecs(P, &xin, nullptr)); 1482 if (ismatis) PetscCall(MatDestroy(&P)); 1483 PetscCall(VecScatterCreate(xin, data->is, data->levels[0]->D, nullptr, &data->levels[0]->scatter)); 1484 PetscCall(VecDestroy(&xin)); 1485 } 1486 if (data->levels[0]->P) { 1487 /* if the pattern is the same and PCSetUp() has previously succeeded, reuse HPDDM buffers and connectivity */ 1488 PetscCall(HPDDM::Schwarz<PetscScalar>::destroy(data->levels[0], pc->setupcalled < 1 || pc->flag == DIFFERENT_NONZERO_PATTERN ? PETSC_TRUE : PETSC_FALSE)); 1489 } 1490 if (!data->levels[0]->P) data->levels[0]->P = new HPDDM::Schwarz<PetscScalar>(); 1491 if (data->log_separate) PetscCall(PetscLogEventBegin(PC_HPDDM_SetUp[0], data->levels[0]->ksp, nullptr, nullptr, nullptr)); 1492 else PetscCall(PetscLogEventBegin(PC_HPDDM_Strc, data->levels[0]->ksp, nullptr, nullptr, nullptr)); 1493 /* HPDDM internal data structure */ 1494 PetscCall(data->levels[0]->P->structure(loc, data->is, sub[0], ismatis ? C : data->aux, data->levels)); 1495 if (!data->log_separate) PetscCall(PetscLogEventEnd(PC_HPDDM_Strc, data->levels[0]->ksp, nullptr, nullptr, nullptr)); 1496 /* matrix pencil of the generalized eigenvalue problem on the overlap (GenEO) */ 1497 if (data->deflation) weighted = data->aux; 1498 else if (!data->B) { 1499 PetscBool cmp[2]; 1500 PetscCall(MatDuplicate(sub[0], MAT_COPY_VALUES, &weighted)); 1501 PetscCall(PetscObjectTypeCompare((PetscObject)weighted, MATNORMAL, cmp)); 1502 if (!cmp[0]) PetscCall(PetscObjectTypeCompare((PetscObject)weighted, MATNORMALHERMITIAN, cmp + 1)); 1503 else cmp[1] = PETSC_FALSE; 1504 if (!cmp[0] && !cmp[1]) PetscCall(MatDiagonalScale(weighted, data->levels[0]->D, data->levels[0]->D)); 1505 else { /* MATNORMAL applies MatDiagonalScale() in a matrix-free fashion, not what is needed since this won't be passed to SLEPc during the eigensolve */ 1506 if (cmp[0]) PetscCall(MatNormalGetMat(weighted, &data->B)); 1507 else PetscCall(MatNormalHermitianGetMat(weighted, &data->B)); 1508 PetscCall(MatDiagonalScale(data->B, nullptr, data->levels[0]->D)); 1509 data->B = nullptr; 1510 flg = PETSC_FALSE; 1511 } 1512 /* neither MatDuplicate() nor MatDiagonaleScale() handles the symmetry options, so propagate the options explicitly */ 1513 /* only useful for -mat_type baij -pc_hpddm_levels_1_st_pc_type cholesky (no problem with MATAIJ or MATSBAIJ) */ 1514 PetscCall(MatPropagateSymmetryOptions(sub[0], weighted)); 1515 } else weighted = data->B; 1516 /* SLEPc is used inside the loaded symbol */ 1517 PetscCall((*loadedSym)(data->levels[0]->P, data->is, ismatis ? C : (algebraic && !block ? sub[0] : data->aux), weighted, data->B, initial, data->levels)); 1518 if (data->share) { 1519 Mat st[2]; 1520 PetscCall(KSPGetOperators(ksp[0], st, st + 1)); 1521 PetscCall(MatCopy(subA[0], st[0], structure)); 1522 if (subA[1] != subA[0] || st[1] != st[0]) PetscCall(MatCopy(subA[1], st[1], SAME_NONZERO_PATTERN)); 1523 } 1524 if (data->log_separate) PetscCall(PetscLogEventEnd(PC_HPDDM_SetUp[0], data->levels[0]->ksp, nullptr, nullptr, nullptr)); 1525 if (ismatis) PetscCall(MatISGetLocalMat(C, &N)); 1526 else N = data->aux; 1527 P = sub[0]; 1528 /* going through the grid hierarchy */ 1529 for (n = 1; n < data->N; ++n) { 1530 if (data->log_separate) PetscCall(PetscLogEventBegin(PC_HPDDM_SetUp[n], data->levels[n]->ksp, nullptr, nullptr, nullptr)); 1531 /* method composed in the loaded symbol since there, SLEPc is used as well */ 1532 PetscTryMethod(data->levels[0]->ksp, "PCHPDDMSetUp_Private_C", (Mat *, Mat *, PetscInt, PetscInt *const, PC_HPDDM_Level **const), (&P, &N, n, &data->N, data->levels)); 1533 if (data->log_separate) PetscCall(PetscLogEventEnd(PC_HPDDM_SetUp[n], data->levels[n]->ksp, nullptr, nullptr, nullptr)); 1534 } 1535 /* reset to NULL to avoid any faulty use */ 1536 PetscCall(PetscObjectComposeFunction((PetscObject)data->levels[0]->ksp, "PCHPDDMSetUp_Private_C", nullptr)); 1537 if (!ismatis) PetscCall(PetscObjectComposeFunction((PetscObject)pc->pmat, "PCHPDDMAlgebraicAuxiliaryMat_C", nullptr)); 1538 else PetscCall(PetscObjectDereference((PetscObject)C)); /* matching PetscObjectReference() above */ 1539 for (n = 0; n < data->N - 1; ++n) 1540 if (data->levels[n]->P) { 1541 /* HPDDM internal work buffers */ 1542 data->levels[n]->P->setBuffer(); 1543 data->levels[n]->P->super::start(); 1544 } 1545 if (ismatis || !subdomains) PetscCall(PCHPDDMDestroySubMatrices_Private(PetscBool3ToBool(data->Neumann), PetscBool(algebraic && !block), sub)); 1546 if (ismatis) data->is = nullptr; 1547 for (n = 0; n < data->N - 1 + (reused > 0); ++n) { 1548 if (data->levels[n]->P) { 1549 PC spc; 1550 1551 /* force the PC to be PCSHELL to do the coarse grid corrections */ 1552 PetscCall(KSPSetSkipPCSetFromOptions(data->levels[n]->ksp, PETSC_TRUE)); 1553 PetscCall(KSPGetPC(data->levels[n]->ksp, &spc)); 1554 PetscCall(PCSetType(spc, PCSHELL)); 1555 PetscCall(PCShellSetContext(spc, data->levels[n])); 1556 PetscCall(PCShellSetSetUp(spc, PCSetUp_HPDDMShell)); 1557 PetscCall(PCShellSetApply(spc, PCApply_HPDDMShell)); 1558 PetscCall(PCShellSetMatApply(spc, PCMatApply_HPDDMShell)); 1559 PetscCall(PCShellSetDestroy(spc, PCDestroy_HPDDMShell)); 1560 if (!data->levels[n]->pc) PetscCall(PCCreate(PetscObjectComm((PetscObject)data->levels[n]->ksp), &data->levels[n]->pc)); 1561 if (n < reused) { 1562 PetscCall(PCSetReusePreconditioner(spc, PETSC_TRUE)); 1563 PetscCall(PCSetReusePreconditioner(data->levels[n]->pc, PETSC_TRUE)); 1564 } 1565 PetscCall(PCSetUp(spc)); 1566 } 1567 } 1568 PetscCall(PetscObjectComposeFunction((PetscObject)pc->pmat, "PCHPDDMAlgebraicAuxiliaryMat_Private_C", nullptr)); 1569 } else flg = reused ? PETSC_FALSE : PETSC_TRUE; 1570 if (!ismatis && subdomains) { 1571 if (flg) PetscCall(KSPGetPC(data->levels[0]->ksp, &inner)); 1572 else inner = data->levels[0]->pc; 1573 if (inner) { 1574 if (!inner->setupcalled) PetscCall(PCSetType(inner, PCASM)); 1575 PetscCall(PCSetFromOptions(inner)); 1576 PetscCall(PetscStrcmp(((PetscObject)inner)->type_name, PCASM, &flg)); 1577 if (flg) { 1578 if (!inner->setupcalled) { /* evaluates to PETSC_FALSE when -pc_hpddm_block_splitting */ 1579 IS sorted; /* PCASM will sort the input IS, duplicate it to return an unmodified (PCHPDDM) input IS */ 1580 PetscCall(ISDuplicate(is[0], &sorted)); 1581 PetscCall(PCASMSetLocalSubdomains(inner, 1, &sorted, &loc)); 1582 PetscCall(PetscObjectDereference((PetscObject)sorted)); 1583 } 1584 if (!PetscBool3ToBool(data->Neumann) && data->N > 1) { /* subdomain matrices are already created for the eigenproblem, reuse them for the fine-level PC */ 1585 PetscCall(PCHPDDMPermute_Private(*is, nullptr, nullptr, sub[0], &P, nullptr)); 1586 PetscCall(PCHPDDMCommunicationAvoidingPCASM_Private(inner, P, algebraic)); 1587 PetscCall(PetscObjectDereference((PetscObject)P)); 1588 } 1589 } 1590 } 1591 if (data->N > 1) PetscCall(PCHPDDMDestroySubMatrices_Private(PetscBool3ToBool(data->Neumann), PetscBool(algebraic && !block), sub)); 1592 } 1593 PetscCall(ISDestroy(&loc)); 1594 } else data->N = 1 + reused; /* enforce this value to 1 + reused if there is no way to build another level */ 1595 if (requested != data->N + reused) { 1596 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, 1597 data->N, pcpre ? pcpre : "")); 1598 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)); 1599 /* cannot use PCDestroy_HPDDMShell() because PCSHELL not set for unassembled levels */ 1600 for (n = data->N - 1; n < requested - 1; ++n) { 1601 if (data->levels[n]->P) { 1602 PetscCall(HPDDM::Schwarz<PetscScalar>::destroy(data->levels[n], PETSC_TRUE)); 1603 PetscCall(VecDestroyVecs(1, &data->levels[n]->v[0])); 1604 PetscCall(VecDestroyVecs(2, &data->levels[n]->v[1])); 1605 PetscCall(MatDestroy(data->levels[n]->V)); 1606 PetscCall(MatDestroy(data->levels[n]->V + 1)); 1607 PetscCall(MatDestroy(data->levels[n]->V + 2)); 1608 PetscCall(VecDestroy(&data->levels[n]->D)); 1609 PetscCall(VecScatterDestroy(&data->levels[n]->scatter)); 1610 } 1611 } 1612 if (reused) { 1613 for (n = reused; n < PETSC_PCHPDDM_MAXLEVELS && data->levels[n]; ++n) { 1614 PetscCall(KSPDestroy(&data->levels[n]->ksp)); 1615 PetscCall(PCDestroy(&data->levels[n]->pc)); 1616 } 1617 } 1618 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, 1619 data->N, reused, data->N, pcpre ? pcpre : "", pcpre ? pcpre : "", data->N); 1620 } 1621 /* these solvers are created after PCSetFromOptions() is called */ 1622 if (pc->setfromoptionscalled) { 1623 for (n = 0; n < data->N; ++n) { 1624 if (data->levels[n]->ksp) PetscCall(KSPSetFromOptions(data->levels[n]->ksp)); 1625 if (data->levels[n]->pc) PetscCall(PCSetFromOptions(data->levels[n]->pc)); 1626 } 1627 pc->setfromoptionscalled = 0; 1628 } 1629 data->N += reused; 1630 if (data->share && swap) { 1631 /* swap back pointers so that variables follow the user-provided numbering */ 1632 std::swap(C, data->aux); 1633 std::swap(uis, data->is); 1634 PetscCall(MatDestroy(&C)); 1635 PetscCall(ISDestroy(&uis)); 1636 } 1637 if (algebraic) PetscCall(MatDestroy(&data->aux)); 1638 if (unsorted && unsorted != is[0]) { 1639 PetscCall(ISCopy(unsorted, data->is)); 1640 PetscCall(ISDestroy(&unsorted)); 1641 } 1642 #if PetscDefined(USE_DEBUG) 1643 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); 1644 if (data->is) { 1645 PetscCall(ISEqualUnsorted(data->is, dis, &flg)); 1646 PetscCall(ISDestroy(&dis)); 1647 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Input IS and output IS are not equal"); 1648 } 1649 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); 1650 if (data->aux) { 1651 PetscCall(MatMultEqual(data->aux, daux, 20, &flg)); 1652 PetscCall(MatDestroy(&daux)); 1653 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Input Mat and output Mat are not equal"); 1654 } 1655 #endif 1656 PetscFunctionReturn(PETSC_SUCCESS); 1657 } 1658 1659 /*@ 1660 PCHPDDMSetCoarseCorrectionType - Sets the coarse correction type. 1661 1662 Collective 1663 1664 Input Parameters: 1665 + pc - preconditioner context 1666 - type - `PC_HPDDM_COARSE_CORRECTION_DEFLATED`, `PC_HPDDM_COARSE_CORRECTION_ADDITIVE`, or `PC_HPDDM_COARSE_CORRECTION_BALANCED` 1667 1668 Options Database Key: 1669 . -pc_hpddm_coarse_correction <deflated, additive, balanced> - type of coarse correction to apply 1670 1671 Level: intermediate 1672 1673 .seealso: `PCHPDDMGetCoarseCorrectionType()`, `PCHPDDM`, `PCHPDDMCoarseCorrectionType` 1674 @*/ 1675 PetscErrorCode PCHPDDMSetCoarseCorrectionType(PC pc, PCHPDDMCoarseCorrectionType type) 1676 { 1677 PetscFunctionBegin; 1678 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1679 PetscValidLogicalCollectiveEnum(pc, type, 2); 1680 PetscTryMethod(pc, "PCHPDDMSetCoarseCorrectionType_C", (PC, PCHPDDMCoarseCorrectionType), (pc, type)); 1681 PetscFunctionReturn(PETSC_SUCCESS); 1682 } 1683 1684 /*@ 1685 PCHPDDMGetCoarseCorrectionType - Gets the coarse correction type. 1686 1687 Input Parameter: 1688 . pc - preconditioner context 1689 1690 Output Parameter: 1691 . type - `PC_HPDDM_COARSE_CORRECTION_DEFLATED`, `PC_HPDDM_COARSE_CORRECTION_ADDITIVE`, or `PC_HPDDM_COARSE_CORRECTION_BALANCED` 1692 1693 Level: intermediate 1694 1695 .seealso: `PCHPDDMSetCoarseCorrectionType()`, `PCHPDDM`, `PCHPDDMCoarseCorrectionType` 1696 @*/ 1697 PetscErrorCode PCHPDDMGetCoarseCorrectionType(PC pc, PCHPDDMCoarseCorrectionType *type) 1698 { 1699 PetscFunctionBegin; 1700 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1701 if (type) { 1702 PetscAssertPointer(type, 2); 1703 PetscUseMethod(pc, "PCHPDDMGetCoarseCorrectionType_C", (PC, PCHPDDMCoarseCorrectionType *), (pc, type)); 1704 } 1705 PetscFunctionReturn(PETSC_SUCCESS); 1706 } 1707 1708 static PetscErrorCode PCHPDDMSetCoarseCorrectionType_HPDDM(PC pc, PCHPDDMCoarseCorrectionType type) 1709 { 1710 PC_HPDDM *data = (PC_HPDDM *)pc->data; 1711 1712 PetscFunctionBegin; 1713 PetscCheck(type >= 0 && type <= 2, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PCHPDDMCoarseCorrectionType %d", type); 1714 data->correction = type; 1715 PetscFunctionReturn(PETSC_SUCCESS); 1716 } 1717 1718 static PetscErrorCode PCHPDDMGetCoarseCorrectionType_HPDDM(PC pc, PCHPDDMCoarseCorrectionType *type) 1719 { 1720 PC_HPDDM *data = (PC_HPDDM *)pc->data; 1721 1722 PetscFunctionBegin; 1723 *type = data->correction; 1724 PetscFunctionReturn(PETSC_SUCCESS); 1725 } 1726 1727 /*@ 1728 PCHPDDMSetSTShareSubKSP - Sets whether the `KSP` in SLEPc `ST` and the fine-level subdomain solver should be shared. 1729 1730 Input Parameters: 1731 + pc - preconditioner context 1732 - share - whether the `KSP` should be shared or not 1733 1734 Note: 1735 This is not the same as `PCSetReusePreconditioner()`. Given certain conditions (visible using -info), a symbolic factorization can be skipped 1736 when using a subdomain `PCType` such as `PCLU` or `PCCHOLESKY`. 1737 1738 Level: advanced 1739 1740 .seealso: `PCHPDDM`, `PCHPDDMGetSTShareSubKSP()` 1741 @*/ 1742 PetscErrorCode PCHPDDMSetSTShareSubKSP(PC pc, PetscBool share) 1743 { 1744 PetscFunctionBegin; 1745 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1746 PetscTryMethod(pc, "PCHPDDMSetSTShareSubKSP_C", (PC, PetscBool), (pc, share)); 1747 PetscFunctionReturn(PETSC_SUCCESS); 1748 } 1749 1750 /*@ 1751 PCHPDDMGetSTShareSubKSP - Gets whether the `KSP` in SLEPc `ST` and the fine-level subdomain solver is shared. 1752 1753 Input Parameter: 1754 . pc - preconditioner context 1755 1756 Output Parameter: 1757 . share - whether the `KSP` is shared or not 1758 1759 Note: 1760 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 1761 when using a subdomain `PCType` such as `PCLU` or `PCCHOLESKY`. 1762 1763 Level: advanced 1764 1765 .seealso: `PCHPDDM`, `PCHPDDMSetSTShareSubKSP()` 1766 @*/ 1767 PetscErrorCode PCHPDDMGetSTShareSubKSP(PC pc, PetscBool *share) 1768 { 1769 PetscFunctionBegin; 1770 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1771 if (share) { 1772 PetscAssertPointer(share, 2); 1773 PetscUseMethod(pc, "PCHPDDMGetSTShareSubKSP_C", (PC, PetscBool *), (pc, share)); 1774 } 1775 PetscFunctionReturn(PETSC_SUCCESS); 1776 } 1777 1778 static PetscErrorCode PCHPDDMSetSTShareSubKSP_HPDDM(PC pc, PetscBool share) 1779 { 1780 PC_HPDDM *data = (PC_HPDDM *)pc->data; 1781 1782 PetscFunctionBegin; 1783 data->share = share; 1784 PetscFunctionReturn(PETSC_SUCCESS); 1785 } 1786 1787 static PetscErrorCode PCHPDDMGetSTShareSubKSP_HPDDM(PC pc, PetscBool *share) 1788 { 1789 PC_HPDDM *data = (PC_HPDDM *)pc->data; 1790 1791 PetscFunctionBegin; 1792 *share = data->share; 1793 PetscFunctionReturn(PETSC_SUCCESS); 1794 } 1795 1796 /*@ 1797 PCHPDDMSetDeflationMat - Sets the deflation space used to assemble a coarser operator. 1798 1799 Input Parameters: 1800 + pc - preconditioner context 1801 . is - index set of the local deflation matrix 1802 - U - deflation sequential matrix stored as a `MATSEQDENSE` 1803 1804 Level: advanced 1805 1806 .seealso: `PCHPDDM`, `PCDeflationSetSpace()`, `PCMGSetRestriction()` 1807 @*/ 1808 PetscErrorCode PCHPDDMSetDeflationMat(PC pc, IS is, Mat U) 1809 { 1810 PetscFunctionBegin; 1811 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1812 PetscValidHeaderSpecific(is, IS_CLASSID, 2); 1813 PetscValidHeaderSpecific(U, MAT_CLASSID, 3); 1814 PetscTryMethod(pc, "PCHPDDMSetDeflationMat_C", (PC, IS, Mat), (pc, is, U)); 1815 PetscFunctionReturn(PETSC_SUCCESS); 1816 } 1817 1818 static PetscErrorCode PCHPDDMSetDeflationMat_HPDDM(PC pc, IS is, Mat U) 1819 { 1820 PetscFunctionBegin; 1821 PetscCall(PCHPDDMSetAuxiliaryMat_Private(pc, is, U, PETSC_TRUE)); 1822 PetscFunctionReturn(PETSC_SUCCESS); 1823 } 1824 1825 PetscErrorCode HPDDMLoadDL_Private(PetscBool *found) 1826 { 1827 PetscBool flg; 1828 char lib[PETSC_MAX_PATH_LEN], dlib[PETSC_MAX_PATH_LEN], dir[PETSC_MAX_PATH_LEN]; 1829 1830 PetscFunctionBegin; 1831 PetscAssertPointer(found, 1); 1832 PetscCall(PetscStrncpy(dir, "${PETSC_LIB_DIR}", sizeof(dir))); 1833 PetscCall(PetscOptionsGetString(nullptr, nullptr, "-hpddm_dir", dir, sizeof(dir), nullptr)); 1834 PetscCall(PetscSNPrintf(lib, sizeof(lib), "%s/libhpddm_petsc", dir)); 1835 PetscCall(PetscDLLibraryRetrieve(PETSC_COMM_SELF, lib, dlib, 1024, found)); 1836 #if defined(SLEPC_LIB_DIR) /* this variable is passed during SLEPc ./configure when PETSc has not been configured */ 1837 if (!*found) { /* with --download-hpddm since slepcconf.h is not yet built (and thus can't be included) */ 1838 PetscCall(PetscStrncpy(dir, HPDDM_STR(SLEPC_LIB_DIR), sizeof(dir))); 1839 PetscCall(PetscSNPrintf(lib, sizeof(lib), "%s/libhpddm_petsc", dir)); 1840 PetscCall(PetscDLLibraryRetrieve(PETSC_COMM_SELF, lib, dlib, 1024, found)); 1841 } 1842 #endif 1843 if (!*found) { /* probable options for this to evaluate to PETSC_TRUE: system inconsistency (libhpddm_petsc moved by user?) or PETSc configured without --download-slepc */ 1844 PetscCall(PetscOptionsGetenv(PETSC_COMM_SELF, "SLEPC_DIR", dir, sizeof(dir), &flg)); 1845 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 */ 1846 PetscCall(PetscSNPrintf(lib, sizeof(lib), "%s/lib/libslepc", dir)); 1847 PetscCall(PetscDLLibraryRetrieve(PETSC_COMM_SELF, lib, dlib, 1024, found)); 1848 PetscCheck(*found, PETSC_COMM_SELF, PETSC_ERR_PLIB, "%s not found but SLEPC_DIR=%s", lib, dir); 1849 PetscCall(PetscDLLibraryAppend(PETSC_COMM_SELF, &PetscDLLibrariesLoaded, dlib)); 1850 PetscCall(PetscSNPrintf(lib, sizeof(lib), "%s/lib/libhpddm_petsc", dir)); /* libhpddm_petsc is always in the same directory as libslepc */ 1851 PetscCall(PetscDLLibraryRetrieve(PETSC_COMM_SELF, lib, dlib, 1024, found)); 1852 } 1853 } 1854 PetscCheck(*found, PETSC_COMM_SELF, PETSC_ERR_PLIB, "%s not found", lib); 1855 PetscCall(PetscDLLibraryAppend(PETSC_COMM_SELF, &PetscDLLibrariesLoaded, dlib)); 1856 PetscFunctionReturn(PETSC_SUCCESS); 1857 } 1858 1859 /*MC 1860 PCHPDDM - Interface with the HPDDM library. 1861 1862 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 1863 AMGe or `PCBDDC` with adaptive selection of constraints. The interface is explained in details in [2021] (see references below) 1864 1865 The matrix to be preconditioned (Pmat) may be unassembled (`MATIS`), assembled (`MATAIJ`, `MATBAIJ`, or `MATSBAIJ`), hierarchical (`MATHTOOL`), or `MATNORMAL`. 1866 1867 For multilevel preconditioning, when using an assembled or hierarchical Pmat, one must provide an auxiliary local `Mat` (unassembled local operator for GenEO) using 1868 `PCHPDDMSetAuxiliaryMat()`. Calling this routine is not needed when using a `MATIS` Pmat, assembly is done internally using `MatConvert()`. 1869 1870 Options Database Keys: 1871 + -pc_hpddm_define_subdomains <true, default=false> - on the finest level, calls `PCASMSetLocalSubdomains()` with the `IS` supplied in `PCHPDDMSetAuxiliaryMat()` 1872 (not relevant with an unassembled Pmat) 1873 . -pc_hpddm_has_neumann <true, default=false> - on the finest level, informs the `PC` that the local Neumann matrix is supplied in `PCHPDDMSetAuxiliaryMat()` 1874 - -pc_hpddm_coarse_correction <type, default=deflated> - determines the `PCHPDDMCoarseCorrectionType` when calling `PCApply()` 1875 1876 Options for subdomain solvers, subdomain eigensolvers (for computing deflation vectors), and the coarse solver can be set using the following options database prefixes. 1877 .vb 1878 -pc_hpddm_levels_%d_pc_ 1879 -pc_hpddm_levels_%d_ksp_ 1880 -pc_hpddm_levels_%d_eps_ 1881 -pc_hpddm_levels_%d_p 1882 -pc_hpddm_levels_%d_mat_type 1883 -pc_hpddm_coarse_ 1884 -pc_hpddm_coarse_p 1885 -pc_hpddm_coarse_mat_type 1886 -pc_hpddm_coarse_mat_filter 1887 .ve 1888 1889 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 1890 -pc_hpddm_coarse_p 2 -pc_hpddm_coarse_mat_type baij will use 10 deflation vectors per subdomain on the fine "level 1", 1891 aggregate the fine subdomains into 4 "level 2" subdomains, then use 10 deflation vectors per subdomain on "level 2", 1892 and assemble the coarse matrix (of dimension 4 x 10 = 40) on two processes as a `MATBAIJ` (default is `MATSBAIJ`). 1893 1894 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. 1895 1896 This preconditioner requires that you build PETSc with SLEPc (``--download-slepc``). By default, the underlying concurrent eigenproblems 1897 are solved using SLEPc shift-and-invert spectral transformation. This is usually what gives the best performance for GenEO, cf. [2011, 2013]. As 1898 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 1899 -pc_hpddm_levels_1_st_type sinvert. There are furthermore three options related to the (subdomain-wise local) eigensolver that are not described in 1900 SLEPc documentation since they are specific to `PCHPDDM`. 1901 .vb 1902 -pc_hpddm_levels_1_st_share_sub_ksp 1903 -pc_hpddm_levels_%d_eps_threshold 1904 -pc_hpddm_levels_1_eps_use_inertia 1905 .ve 1906 1907 The first option from the list only applies to the fine-level eigensolver, see `PCHPDDMSetSTShareSubKSP()`. The second option from the list is 1908 used to filter eigenmodes retrieved after convergence of `EPSSolve()` at "level N" such that eigenvectors used to define a "level N+1" coarse 1909 correction are associated to eigenvalues whose magnitude are lower or equal than -pc_hpddm_levels_N_eps_threshold. When using an `EPS` which cannot 1910 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 1911 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 1912 to supply -pc_hpddm_levels_1_eps_nev. This last option also only applies to the fine-level (N = 1) eigensolver. 1913 1914 References: 1915 + 2011 - A robust two-level domain decomposition preconditioner for systems of PDEs. Spillane, Dolean, Hauret, Nataf, Pechstein, and Scheichl. Comptes Rendus Mathematique. 1916 . 2013 - Scalable domain decomposition preconditioners for heterogeneous elliptic problems. Jolivet, Hecht, Nataf, and Prud'homme. SC13. 1917 . 2015 - An introduction to domain decomposition methods: algorithms, theory, and parallel implementation. Dolean, Jolivet, and Nataf. SIAM. 1918 . 2019 - A multilevel Schwarz preconditioner based on a hierarchy of robust coarse spaces. Al Daas, Grigori, Jolivet, and Tournier. SIAM Journal on Scientific Computing. 1919 . 2021 - KSPHPDDM and PCHPDDM: extending PETSc with advanced Krylov methods and robust multilevel overlapping Schwarz preconditioners. Jolivet, Roman, and Zampini. Computer & Mathematics with Applications. 1920 . 2022a - A robust algebraic domain decomposition preconditioner for sparse normal equations. Al Daas, Jolivet, and Scott. SIAM Journal on Scientific Computing. 1921 - 2022b - A robust algebraic multilevel domain decomposition preconditioner for sparse symmetric positive definite matrices. Al Daas and Jolivet. 1922 1923 Level: intermediate 1924 1925 .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCHPDDMSetAuxiliaryMat()`, `MATIS`, `PCBDDC`, `PCDEFLATION`, `PCTELESCOPE`, `PCASM`, 1926 `PCHPDDMSetCoarseCorrectionType()`, `PCHPDDMHasNeumannMat()`, `PCHPDDMSetRHSMat()`, `PCHPDDMSetDeflationMat()`, `PCHPDDMSetSTShareSubKSP()`, 1927 `PCHPDDMGetSTShareSubKSP()`, `PCHPDDMGetCoarseCorrectionType()`, `PCHPDDMGetComplexities()` 1928 M*/ 1929 PETSC_EXTERN PetscErrorCode PCCreate_HPDDM(PC pc) 1930 { 1931 PC_HPDDM *data; 1932 PetscBool found; 1933 1934 PetscFunctionBegin; 1935 if (!loadedSym) { 1936 PetscCall(HPDDMLoadDL_Private(&found)); 1937 if (found) PetscCall(PetscDLLibrarySym(PETSC_COMM_SELF, &PetscDLLibrariesLoaded, nullptr, "PCHPDDM_Internal", (void **)&loadedSym)); 1938 } 1939 PetscCheck(loadedSym, PETSC_COMM_SELF, PETSC_ERR_PLIB, "PCHPDDM_Internal symbol not found in loaded libhpddm_petsc"); 1940 PetscCall(PetscNew(&data)); 1941 pc->data = data; 1942 data->Neumann = PETSC_BOOL3_UNKNOWN; 1943 pc->ops->reset = PCReset_HPDDM; 1944 pc->ops->destroy = PCDestroy_HPDDM; 1945 pc->ops->setfromoptions = PCSetFromOptions_HPDDM; 1946 pc->ops->setup = PCSetUp_HPDDM; 1947 pc->ops->apply = PCApply_HPDDM; 1948 pc->ops->matapply = PCMatApply_HPDDM; 1949 pc->ops->view = PCView_HPDDM; 1950 pc->ops->presolve = PCPreSolve_HPDDM; 1951 1952 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetAuxiliaryMat_C", PCHPDDMSetAuxiliaryMat_HPDDM)); 1953 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMHasNeumannMat_C", PCHPDDMHasNeumannMat_HPDDM)); 1954 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetRHSMat_C", PCHPDDMSetRHSMat_HPDDM)); 1955 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetCoarseCorrectionType_C", PCHPDDMSetCoarseCorrectionType_HPDDM)); 1956 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMGetCoarseCorrectionType_C", PCHPDDMGetCoarseCorrectionType_HPDDM)); 1957 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetSTShareSubKSP_C", PCHPDDMSetSTShareSubKSP_HPDDM)); 1958 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMGetSTShareSubKSP_C", PCHPDDMGetSTShareSubKSP_HPDDM)); 1959 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetDeflationMat_C", PCHPDDMSetDeflationMat_HPDDM)); 1960 PetscFunctionReturn(PETSC_SUCCESS); 1961 } 1962 1963 /*@C 1964 PCHPDDMInitializePackage - This function initializes everything in the `PCHPDDM` package. It is called from `PCInitializePackage()`. 1965 1966 Level: developer 1967 1968 .seealso: `PetscInitialize()` 1969 @*/ 1970 PetscErrorCode PCHPDDMInitializePackage(void) 1971 { 1972 char ename[32]; 1973 PetscInt i; 1974 1975 PetscFunctionBegin; 1976 if (PCHPDDMPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS); 1977 PCHPDDMPackageInitialized = PETSC_TRUE; 1978 PetscCall(PetscRegisterFinalize(PCHPDDMFinalizePackage)); 1979 /* general events registered once during package initialization */ 1980 /* some of these events are not triggered in libpetsc, */ 1981 /* but rather directly in libhpddm_petsc, */ 1982 /* which is in charge of performing the following operations */ 1983 1984 /* domain decomposition structure from Pmat sparsity pattern */ 1985 PetscCall(PetscLogEventRegister("PCHPDDMStrc", PC_CLASSID, &PC_HPDDM_Strc)); 1986 /* Galerkin product, redistribution, and setup (not triggered in libpetsc) */ 1987 PetscCall(PetscLogEventRegister("PCHPDDMPtAP", PC_CLASSID, &PC_HPDDM_PtAP)); 1988 /* Galerkin product with summation, redistribution, and setup (not triggered in libpetsc) */ 1989 PetscCall(PetscLogEventRegister("PCHPDDMPtBP", PC_CLASSID, &PC_HPDDM_PtBP)); 1990 /* next level construction using PtAP and PtBP (not triggered in libpetsc) */ 1991 PetscCall(PetscLogEventRegister("PCHPDDMNext", PC_CLASSID, &PC_HPDDM_Next)); 1992 static_assert(PETSC_PCHPDDM_MAXLEVELS <= 9, "PETSC_PCHPDDM_MAXLEVELS value is too high"); 1993 for (i = 1; i < PETSC_PCHPDDM_MAXLEVELS; ++i) { 1994 PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCHPDDMSetUp L%1" PetscInt_FMT, i)); 1995 /* events during a PCSetUp() at level #i _except_ the assembly */ 1996 /* of the Galerkin operator of the coarser level #(i + 1) */ 1997 PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_HPDDM_SetUp[i - 1])); 1998 PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCHPDDMSolve L%1" PetscInt_FMT, i)); 1999 /* events during a PCApply() at level #i _except_ */ 2000 /* the KSPSolve() of the coarser level #(i + 1) */ 2001 PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_HPDDM_Solve[i - 1])); 2002 } 2003 PetscFunctionReturn(PETSC_SUCCESS); 2004 } 2005 2006 /*@C 2007 PCHPDDMFinalizePackage - This function frees everything from the `PCHPDDM` package. It is called from `PetscFinalize()`. 2008 2009 Level: developer 2010 2011 .seealso: `PetscFinalize()` 2012 @*/ 2013 PetscErrorCode PCHPDDMFinalizePackage(void) 2014 { 2015 PetscFunctionBegin; 2016 PCHPDDMPackageInitialized = PETSC_FALSE; 2017 PetscFunctionReturn(PETSC_SUCCESS); 2018 } 2019