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