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