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