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