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