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