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