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