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