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