1 #include <petscsf.h> 2 #include <petsc/private/vecimpl.h> 3 #include <petsc/private/matimpl.h> 4 #include <petsc/private/petschpddm.h> /*I "petscpc.h" I*/ 5 #include <petsc/private/pcimpl.h> 6 #include <petsc/private/dmimpl.h> /* this must be included after petschpddm.h so that DM_MAX_WORK_VECTORS is not defined */ 7 /* otherwise, it is assumed that one is compiling libhpddm_petsc => circular dependency */ 8 9 static PetscErrorCode (*loadedSym)(HPDDM::Schwarz<PetscScalar> *const, IS, Mat, Mat, Mat, std::vector<Vec>, PC_HPDDM_Level **const) = nullptr; 10 11 static PetscBool PCHPDDMPackageInitialized = PETSC_FALSE; 12 13 PetscLogEvent PC_HPDDM_Strc; 14 PetscLogEvent PC_HPDDM_PtAP; 15 PetscLogEvent PC_HPDDM_PtBP; 16 PetscLogEvent PC_HPDDM_Next; 17 PetscLogEvent PC_HPDDM_SetUp[PETSC_PCHPDDM_MAXLEVELS]; 18 PetscLogEvent PC_HPDDM_Solve[PETSC_PCHPDDM_MAXLEVELS]; 19 20 const char *const PCHPDDMCoarseCorrectionTypes[] = {"DEFLATED", "ADDITIVE", "BALANCED", "NONE", "PCHPDDMCoarseCorrectionType", "PC_HPDDM_COARSE_CORRECTION_", nullptr}; 21 const char *const PCHPDDMSchurPreTypes[] = {"LEAST_SQUARES", "GENEO", "PCHPDDMSchurPreType", "PC_HPDDM_SCHUR_PRE", nullptr}; 22 23 static PetscErrorCode PCReset_HPDDM(PC pc) 24 { 25 PC_HPDDM *data = (PC_HPDDM *)pc->data; 26 27 PetscFunctionBegin; 28 if (data->levels) { 29 for (PetscInt i = 0; i < PETSC_PCHPDDM_MAXLEVELS && data->levels[i]; ++i) { 30 PetscCall(KSPDestroy(&data->levels[i]->ksp)); 31 PetscCall(PCDestroy(&data->levels[i]->pc)); 32 PetscCall(PetscFree(data->levels[i])); 33 } 34 PetscCall(PetscFree(data->levels)); 35 } 36 PetscCall(ISDestroy(&data->is)); 37 PetscCall(MatDestroy(&data->aux)); 38 PetscCall(MatDestroy(&data->B)); 39 PetscCall(VecDestroy(&data->normal)); 40 data->correction = PC_HPDDM_COARSE_CORRECTION_DEFLATED; 41 data->Neumann = PETSC_BOOL3_UNKNOWN; 42 data->deflation = PETSC_FALSE; 43 data->setup = nullptr; 44 data->setup_ctx = nullptr; 45 PetscFunctionReturn(PETSC_SUCCESS); 46 } 47 48 static PetscErrorCode PCDestroy_HPDDM(PC pc) 49 { 50 PC_HPDDM *data = (PC_HPDDM *)pc->data; 51 52 PetscFunctionBegin; 53 PetscCall(PCReset_HPDDM(pc)); 54 PetscCall(PetscFree(data)); 55 PetscCall(PetscObjectChangeTypeName((PetscObject)pc, nullptr)); 56 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetAuxiliaryMat_C", nullptr)); 57 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMHasNeumannMat_C", nullptr)); 58 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetRHSMat_C", nullptr)); 59 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetCoarseCorrectionType_C", nullptr)); 60 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMGetCoarseCorrectionType_C", nullptr)); 61 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetSTShareSubKSP_C", nullptr)); 62 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMGetSTShareSubKSP_C", nullptr)); 63 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetDeflationMat_C", nullptr)); 64 PetscCall(PetscObjectCompose((PetscObject)pc, "_PCHPDDM_Schur", nullptr)); 65 PetscFunctionReturn(PETSC_SUCCESS); 66 } 67 68 static inline PetscErrorCode PCHPDDMSetAuxiliaryMat_Private(PC pc, IS is, Mat A, PetscBool deflation) 69 { 70 PC_HPDDM *data = (PC_HPDDM *)pc->data; 71 PCHPDDMCoarseCorrectionType type = data->correction; 72 73 PetscFunctionBegin; 74 PetscValidLogicalCollectiveBool(pc, deflation, 4); 75 if (is && A) { 76 PetscInt m[2]; 77 78 PetscCall(ISGetLocalSize(is, m)); 79 PetscCall(MatGetLocalSize(A, m + 1, nullptr)); 80 PetscCheck(m[0] == m[1], PETSC_COMM_SELF, PETSC_ERR_USER_INPUT, "Inconsistent IS and Mat sizes (%" PetscInt_FMT " v. %" PetscInt_FMT ")", m[0], m[1]); 81 } 82 if (is) { 83 PetscCall(PetscObjectReference((PetscObject)is)); 84 if (data->is) { /* new overlap definition resets the PC */ 85 PetscCall(PCReset_HPDDM(pc)); 86 pc->setfromoptionscalled = 0; 87 pc->setupcalled = 0; 88 data->correction = type; 89 } 90 PetscCall(ISDestroy(&data->is)); 91 data->is = is; 92 } 93 if (A) { 94 PetscCall(PetscObjectReference((PetscObject)A)); 95 PetscCall(MatDestroy(&data->aux)); 96 data->aux = A; 97 } 98 data->deflation = deflation; 99 PetscFunctionReturn(PETSC_SUCCESS); 100 } 101 102 static inline PetscErrorCode PCHPDDMSplittingMatNormal_Private(Mat A, IS *is, Mat *splitting[]) 103 { 104 Mat *sub; 105 IS zero; 106 107 PetscFunctionBegin; 108 PetscCall(MatSetOption(A, MAT_SUBMAT_SINGLEIS, PETSC_TRUE)); 109 PetscCall(MatCreateSubMatrices(A, 1, is + 2, is, MAT_INITIAL_MATRIX, splitting)); 110 PetscCall(MatCreateSubMatrices(**splitting, 1, is + 2, is + 1, MAT_INITIAL_MATRIX, &sub)); 111 PetscCall(MatFindZeroRows(*sub, &zero)); 112 PetscCall(MatDestroySubMatrices(1, &sub)); 113 PetscCall(MatSetOption(**splitting, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE)); 114 PetscCall(MatZeroRowsIS(**splitting, zero, 0.0, nullptr, nullptr)); 115 PetscCall(ISDestroy(&zero)); 116 PetscFunctionReturn(PETSC_SUCCESS); 117 } 118 119 static inline PetscErrorCode PCHPDDMSetAuxiliaryMatNormal_Private(PC pc, Mat A, Mat N, Mat *B, const char *pcpre, Vec *diagonal = nullptr, Mat B01 = nullptr) 120 { 121 PC_HPDDM *data = (PC_HPDDM *)pc->data; 122 Mat *splitting[2] = {}, aux; 123 Vec d; 124 IS is[3]; 125 PetscReal norm; 126 PetscBool flg; 127 char type[256] = {}; /* same size as in src/ksp/pc/interface/pcset.c */ 128 129 PetscFunctionBegin; 130 if (!B01) PetscCall(MatConvert(N, MATAIJ, MAT_INITIAL_MATRIX, B)); 131 else PetscCall(MatTransposeMatMult(B01, A, MAT_INITIAL_MATRIX, PETSC_DETERMINE, B)); 132 PetscCall(MatEliminateZeros(*B, PETSC_TRUE)); 133 PetscCall(ISCreateStride(PETSC_COMM_SELF, A->cmap->n, A->cmap->rstart, 1, is)); 134 PetscCall(MatIncreaseOverlap(*B, 1, is, 1)); 135 PetscCall(ISCreateStride(PETSC_COMM_SELF, A->cmap->n, A->cmap->rstart, 1, is + 2)); 136 PetscCall(ISEmbed(is[0], is[2], PETSC_TRUE, is + 1)); 137 PetscCall(ISDestroy(is + 2)); 138 PetscCall(ISCreateStride(PETSC_COMM_SELF, A->rmap->N, 0, 1, is + 2)); 139 PetscCall(PCHPDDMSplittingMatNormal_Private(A, is, &splitting[0])); 140 if (B01) { 141 PetscCall(PCHPDDMSplittingMatNormal_Private(B01, is, &splitting[1])); 142 PetscCall(MatDestroy(&B01)); 143 } 144 PetscCall(ISDestroy(is + 2)); 145 PetscCall(ISDestroy(is + 1)); 146 PetscCall(PetscOptionsGetString(nullptr, pcpre, "-pc_hpddm_levels_1_sub_pc_type", type, sizeof(type), nullptr)); 147 PetscCall(PetscStrcmp(type, PCQR, &flg)); 148 if (!flg) { 149 Mat conjugate = *splitting[splitting[1] ? 1 : 0]; 150 151 if (PetscDefined(USE_COMPLEX) && !splitting[1]) { 152 PetscCall(MatDuplicate(*splitting[0], MAT_COPY_VALUES, &conjugate)); 153 PetscCall(MatConjugate(conjugate)); 154 } 155 PetscCall(MatTransposeMatMult(conjugate, *splitting[0], MAT_INITIAL_MATRIX, PETSC_DETERMINE, &aux)); 156 if (PetscDefined(USE_COMPLEX) && !splitting[1]) PetscCall(MatDestroy(&conjugate)); 157 else if (splitting[1]) PetscCall(MatDestroySubMatrices(1, &splitting[1])); 158 PetscCall(MatNorm(aux, NORM_FROBENIUS, &norm)); 159 PetscCall(MatSetOption(aux, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE)); 160 if (diagonal) { 161 PetscReal norm; 162 163 PetscCall(VecScale(*diagonal, -1.0)); 164 PetscCall(VecNorm(*diagonal, NORM_INFINITY, &norm)); 165 if (norm > PETSC_SMALL) { 166 PetscSF scatter; 167 PetscInt n; 168 169 PetscCall(ISGetLocalSize(*is, &n)); 170 PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)pc), n, PETSC_DECIDE, &d)); 171 PetscCall(VecScatterCreate(*diagonal, *is, d, nullptr, &scatter)); 172 PetscCall(VecScatterBegin(scatter, *diagonal, d, INSERT_VALUES, SCATTER_FORWARD)); 173 PetscCall(VecScatterEnd(scatter, *diagonal, d, INSERT_VALUES, SCATTER_FORWARD)); 174 PetscCall(PetscSFDestroy(&scatter)); 175 PetscCall(MatDiagonalSet(aux, d, ADD_VALUES)); 176 PetscCall(VecDestroy(&d)); 177 } else PetscCall(VecDestroy(diagonal)); 178 } 179 if (!diagonal) PetscCall(MatShift(aux, PETSC_SMALL * norm)); 180 } else { 181 PetscBool flg; 182 183 PetscCheck(!splitting[1], PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Cannot use PCQR when A01 != A10^T"); 184 if (diagonal) { 185 PetscCall(VecNorm(*diagonal, NORM_INFINITY, &norm)); 186 PetscCheck(norm < PETSC_SMALL, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Nonzero diagonal A11 block"); 187 PetscCall(VecDestroy(diagonal)); 188 } 189 PetscCall(PetscObjectTypeCompare((PetscObject)N, MATNORMAL, &flg)); 190 if (flg) PetscCall(MatCreateNormal(*splitting[0], &aux)); 191 else PetscCall(MatCreateNormalHermitian(*splitting[0], &aux)); 192 } 193 PetscCall(MatDestroySubMatrices(1, &splitting[0])); 194 PetscCall(PCHPDDMSetAuxiliaryMat(pc, *is, aux, nullptr, nullptr)); 195 data->Neumann = PETSC_BOOL3_TRUE; 196 PetscCall(ISDestroy(is)); 197 PetscCall(MatDestroy(&aux)); 198 PetscFunctionReturn(PETSC_SUCCESS); 199 } 200 201 static PetscErrorCode PCHPDDMSetAuxiliaryMat_HPDDM(PC pc, IS is, Mat A, PetscErrorCode (*setup)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void *), void *setup_ctx) 202 { 203 PC_HPDDM *data = (PC_HPDDM *)pc->data; 204 205 PetscFunctionBegin; 206 PetscCall(PCHPDDMSetAuxiliaryMat_Private(pc, is, A, PETSC_FALSE)); 207 if (setup) { 208 data->setup = setup; 209 data->setup_ctx = setup_ctx; 210 } 211 PetscFunctionReturn(PETSC_SUCCESS); 212 } 213 214 /*@C 215 PCHPDDMSetAuxiliaryMat - Sets the auxiliary matrix used by `PCHPDDM` for the concurrent GenEO problems at the finest level. 216 217 Input Parameters: 218 + pc - preconditioner context 219 . is - index set of the local auxiliary, e.g., Neumann, matrix 220 . A - auxiliary sequential matrix 221 . setup - function for generating the auxiliary matrix entries, may be `NULL` 222 - ctx - context for `setup`, may be `NULL` 223 224 Calling sequence of `setup`: 225 + J - matrix whose values are to be set 226 . t - time 227 . X - linearization point 228 . X_t - time-derivative of the linearization point 229 . s - step 230 . ovl - index set of the local auxiliary, e.g., Neumann, matrix 231 - ctx - context for `setup`, may be `NULL` 232 233 Level: intermediate 234 235 Note: 236 As an example, in a finite element context with nonoverlapping subdomains plus (overlapping) ghost elements, this could be the unassembled (Neumann) 237 local overlapping operator. As opposed to the assembled (Dirichlet) local overlapping operator obtained by summing neighborhood contributions 238 at the interface of ghost elements. 239 240 Fortran Notes: 241 Only `PETSC_NULL_FUNCTION` is supported for `setup` and `ctx` is never accessed 242 243 .seealso: [](ch_ksp), `PCHPDDM`, `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCHPDDMSetRHSMat()`, `MATIS` 244 @*/ 245 PetscErrorCode PCHPDDMSetAuxiliaryMat(PC pc, IS is, Mat A, PetscErrorCode (*setup)(Mat J, PetscReal t, Vec X, Vec X_t, PetscReal s, IS ovl, void *ctx), void *ctx) 246 { 247 PetscFunctionBegin; 248 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 249 if (is) PetscValidHeaderSpecific(is, IS_CLASSID, 2); 250 if (A) PetscValidHeaderSpecific(A, MAT_CLASSID, 3); 251 PetscTryMethod(pc, "PCHPDDMSetAuxiliaryMat_C", (PC, IS, Mat, PetscErrorCode (*)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void *), void *), (pc, is, A, setup, ctx)); 252 PetscFunctionReturn(PETSC_SUCCESS); 253 } 254 255 static PetscErrorCode PCHPDDMHasNeumannMat_HPDDM(PC pc, PetscBool has) 256 { 257 PC_HPDDM *data = (PC_HPDDM *)pc->data; 258 259 PetscFunctionBegin; 260 data->Neumann = PetscBoolToBool3(has); 261 PetscFunctionReturn(PETSC_SUCCESS); 262 } 263 264 /*@ 265 PCHPDDMHasNeumannMat - Informs `PCHPDDM` that the `Mat` passed to `PCHPDDMSetAuxiliaryMat()` is the local Neumann matrix. 266 267 Input Parameters: 268 + pc - preconditioner context 269 - has - Boolean value 270 271 Level: intermediate 272 273 Notes: 274 This may be used to bypass a call to `MatCreateSubMatrices()` and to `MatConvert()` for `MATSBAIJ` matrices. 275 276 If a function is composed with DMCreateNeumannOverlap_C 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`. 277 278 .seealso: [](ch_ksp), `PCHPDDM`, `PCHPDDMSetAuxiliaryMat()` 279 @*/ 280 PetscErrorCode PCHPDDMHasNeumannMat(PC pc, PetscBool has) 281 { 282 PetscFunctionBegin; 283 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 284 PetscTryMethod(pc, "PCHPDDMHasNeumannMat_C", (PC, PetscBool), (pc, has)); 285 PetscFunctionReturn(PETSC_SUCCESS); 286 } 287 288 static PetscErrorCode PCHPDDMSetRHSMat_HPDDM(PC pc, Mat B) 289 { 290 PC_HPDDM *data = (PC_HPDDM *)pc->data; 291 292 PetscFunctionBegin; 293 PetscCall(PetscObjectReference((PetscObject)B)); 294 PetscCall(MatDestroy(&data->B)); 295 data->B = B; 296 PetscFunctionReturn(PETSC_SUCCESS); 297 } 298 299 /*@ 300 PCHPDDMSetRHSMat - Sets the right-hand side matrix used by `PCHPDDM` for the concurrent GenEO problems at the finest level. 301 302 Input Parameters: 303 + pc - preconditioner context 304 - B - right-hand side sequential matrix 305 306 Level: advanced 307 308 Note: 309 Must be used in conjunction with `PCHPDDMSetAuxiliaryMat`(N), so that Nv = lambda Bv is solved using `EPSSetOperators`(N, B). 310 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. 311 312 .seealso: [](ch_ksp), `PCHPDDMSetAuxiliaryMat()`, `PCHPDDM` 313 @*/ 314 PetscErrorCode PCHPDDMSetRHSMat(PC pc, Mat B) 315 { 316 PetscFunctionBegin; 317 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 318 if (B) { 319 PetscValidHeaderSpecific(B, MAT_CLASSID, 2); 320 PetscTryMethod(pc, "PCHPDDMSetRHSMat_C", (PC, Mat), (pc, B)); 321 } 322 PetscFunctionReturn(PETSC_SUCCESS); 323 } 324 325 static PetscErrorCode PCSetFromOptions_HPDDM(PC pc, PetscOptionItems PetscOptionsObject) 326 { 327 PC_HPDDM *data = (PC_HPDDM *)pc->data; 328 PC_HPDDM_Level **levels = data->levels; 329 char prefix[256]; 330 int i = 1; 331 PetscMPIInt size, previous; 332 PetscInt n, overlap = 1; 333 PCHPDDMCoarseCorrectionType type; 334 PetscBool flg = PETSC_TRUE, set; 335 336 PetscFunctionBegin; 337 if (!data->levels) { 338 PetscCall(PetscCalloc1(PETSC_PCHPDDM_MAXLEVELS, &levels)); 339 data->levels = levels; 340 } 341 PetscOptionsHeadBegin(PetscOptionsObject, "PCHPDDM options"); 342 PetscCall(PetscOptionsBoundedInt("-pc_hpddm_harmonic_overlap", "Overlap prior to computing local harmonic extensions", "PCHPDDM", overlap, &overlap, &set, 1)); 343 if (!set) overlap = -1; 344 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size)); 345 previous = size; 346 while (i < PETSC_PCHPDDM_MAXLEVELS) { 347 PetscInt p = 1; 348 349 if (!data->levels[i - 1]) PetscCall(PetscNew(data->levels + i - 1)); 350 data->levels[i - 1]->parent = data; 351 /* if the previous level has a single process, it is not possible to coarsen further */ 352 if (previous == 1 || !flg) break; 353 data->levels[i - 1]->nu = 0; 354 data->levels[i - 1]->threshold = -1.0; 355 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_%d_eps_nev", i)); 356 PetscCall(PetscOptionsBoundedInt(prefix, "Local number of deflation vectors computed by SLEPc", "EPSSetDimensions", data->levels[i - 1]->nu, &data->levels[i - 1]->nu, nullptr, 0)); 357 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_%d_eps_threshold", i)); 358 PetscCall(PetscOptionsReal(prefix, "Local threshold for selecting deflation vectors returned by SLEPc", "PCHPDDM", data->levels[i - 1]->threshold, &data->levels[i - 1]->threshold, nullptr)); 359 if (i == 1) { 360 PetscCheck(overlap == -1 || PetscAbsReal(data->levels[i - 1]->threshold + static_cast<PetscReal>(1.0)) < PETSC_MACHINE_EPSILON, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot supply both -pc_hpddm_levels_1_eps_threshold and -pc_hpddm_harmonic_overlap"); 361 if (overlap != -1) { 362 PetscInt nsv = 0; 363 364 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_%d_svd_nsv", i)); 365 PetscCall(PetscOptionsBoundedInt(prefix, "Local number of deflation vectors computed by SLEPc", "SVDSetDimensions", nsv, &nsv, nullptr, 0)); 366 PetscCheck(bool(data->levels[0]->nu) != bool(nsv), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot supply %s -pc_hpddm_levels_1_eps_nev %s -pc_hpddm_levels_1_svd_nsv", nsv ? "both" : "neither", nsv ? "and" : "nor"); 367 if (nsv) { 368 data->levels[0]->nu = nsv; 369 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_%d_svd_relative_threshold", i)); 370 } else PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_%d_eps_relative_threshold", i)); 371 PetscCall(PetscOptionsReal(prefix, "Local relative threshold for selecting deflation vectors returned by SLEPc", "PCHPDDM", data->levels[0]->threshold, &data->levels[0]->threshold, nullptr)); 372 } 373 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_1_st_share_sub_ksp")); 374 PetscCall(PetscOptionsBool(prefix, "Shared KSP between SLEPc ST and the fine-level subdomain solver", "PCHPDDMSetSTShareSubKSP", PETSC_FALSE, &data->share, nullptr)); 375 } 376 /* if there is no prescribed coarsening, just break out of the loop */ 377 if (data->levels[i - 1]->threshold <= PetscReal() && data->levels[i - 1]->nu <= 0 && !(data->deflation && i == 1)) break; 378 else { 379 ++i; 380 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_%d_eps_nev", i)); 381 PetscCall(PetscOptionsHasName(PetscOptionsObject->options, PetscOptionsObject->prefix, prefix, &flg)); 382 if (!flg) { 383 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_%d_eps_threshold", i)); 384 PetscCall(PetscOptionsHasName(PetscOptionsObject->options, PetscOptionsObject->prefix, prefix, &flg)); 385 } 386 if (flg) { 387 /* if there are coarsening options for the next level, then register it */ 388 /* otherwise, don't to avoid having both options levels_N_p and coarse_p */ 389 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_levels_%d_p", i)); 390 PetscCall(PetscOptionsRangeInt(prefix, "Number of processes used to assemble the coarse operator at this level", "PCHPDDM", p, &p, &flg, 1, PetscMax(1, previous / 2))); 391 previous = p; 392 } 393 } 394 } 395 data->N = i; 396 n = 1; 397 if (i > 1) { 398 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_coarse_p")); 399 PetscCall(PetscOptionsRangeInt(prefix, "Number of processes used to assemble the coarsest operator", "PCHPDDM", n, &n, nullptr, 1, PetscMax(1, previous / 2))); 400 #if PetscDefined(HAVE_MUMPS) 401 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "pc_hpddm_coarse_")); 402 PetscCall(PetscOptionsHasName(nullptr, prefix, "-mat_mumps_use_omp_threads", &flg)); 403 if (flg) { 404 char type[64]; /* same size as in src/ksp/pc/impls/factor/factimpl.c */ 405 406 PetscCall(PetscStrncpy(type, n > 1 && PetscDefined(HAVE_MUMPS) ? MATSOLVERMUMPS : MATSOLVERPETSC, sizeof(type))); /* default solver for a MatMPIAIJ or a MatSeqAIJ */ 407 PetscCall(PetscOptionsGetString(nullptr, prefix, "-pc_factor_mat_solver_type", type, sizeof(type), nullptr)); 408 PetscCall(PetscStrcmp(type, MATSOLVERMUMPS, &flg)); 409 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "-%smat_mumps_use_omp_threads and -%spc_factor_mat_solver_type != %s", prefix, prefix, MATSOLVERMUMPS); 410 size = n; 411 n = -1; 412 PetscCall(PetscOptionsGetInt(nullptr, prefix, "-mat_mumps_use_omp_threads", &n, nullptr)); 413 PetscCheck(n >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Need to specify a positive integer for -%smat_mumps_use_omp_threads", prefix); 414 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" : ""); 415 } 416 #endif 417 PetscCall(PetscOptionsEnum("-pc_hpddm_coarse_correction", "Type of coarse correction applied each iteration", "PCHPDDMSetCoarseCorrectionType", PCHPDDMCoarseCorrectionTypes, (PetscEnum)data->correction, (PetscEnum *)&type, &flg)); 418 if (flg) PetscCall(PCHPDDMSetCoarseCorrectionType(pc, type)); 419 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_has_neumann")); 420 PetscCall(PetscOptionsBool(prefix, "Is the auxiliary Mat the local Neumann matrix?", "PCHPDDMHasNeumannMat", PetscBool3ToBool(data->Neumann), &flg, &set)); 421 if (set) data->Neumann = PetscBoolToBool3(flg); 422 data->log_separate = PETSC_FALSE; 423 if (PetscDefined(USE_LOG)) { 424 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-pc_hpddm_log_separate")); 425 PetscCall(PetscOptionsBool(prefix, "Log events level by level instead of inside PCSetUp()/KSPSolve()", nullptr, data->log_separate, &data->log_separate, nullptr)); 426 } 427 } 428 PetscOptionsHeadEnd(); 429 while (i < PETSC_PCHPDDM_MAXLEVELS && data->levels[i]) PetscCall(PetscFree(data->levels[i++])); 430 PetscFunctionReturn(PETSC_SUCCESS); 431 } 432 433 static PetscErrorCode PCApply_HPDDM(PC pc, Vec x, Vec y) 434 { 435 PC_HPDDM *data = (PC_HPDDM *)pc->data; 436 437 PetscFunctionBegin; 438 PetscCall(PetscCitationsRegister(HPDDMCitation, &HPDDMCite)); 439 PetscCheck(data->levels[0]->ksp, PETSC_COMM_SELF, PETSC_ERR_PLIB, "No KSP attached to PCHPDDM"); 440 if (data->log_separate) PetscCall(PetscLogEventBegin(PC_HPDDM_Solve[0], data->levels[0]->ksp, nullptr, nullptr, nullptr)); /* coarser-level events are directly triggered in HPDDM */ 441 PetscCall(KSPSolve(data->levels[0]->ksp, x, y)); 442 if (data->log_separate) PetscCall(PetscLogEventEnd(PC_HPDDM_Solve[0], data->levels[0]->ksp, nullptr, nullptr, nullptr)); 443 PetscFunctionReturn(PETSC_SUCCESS); 444 } 445 446 static PetscErrorCode PCMatApply_HPDDM(PC pc, Mat X, Mat Y) 447 { 448 PC_HPDDM *data = (PC_HPDDM *)pc->data; 449 450 PetscFunctionBegin; 451 PetscCall(PetscCitationsRegister(HPDDMCitation, &HPDDMCite)); 452 PetscCheck(data->levels[0]->ksp, PETSC_COMM_SELF, PETSC_ERR_PLIB, "No KSP attached to PCHPDDM"); 453 PetscCall(KSPMatSolve(data->levels[0]->ksp, X, Y)); 454 PetscFunctionReturn(PETSC_SUCCESS); 455 } 456 457 /*@ 458 PCHPDDMGetComplexities - Computes the grid and operator complexities. 459 460 Collective 461 462 Input Parameter: 463 . pc - preconditioner context 464 465 Output Parameters: 466 + gc - grid complexity $ \sum_i m_i / m_1 $ 467 - oc - operator complexity $ \sum_i nnz_i / nnz_1 $ 468 469 Level: advanced 470 471 .seealso: [](ch_ksp), `PCMGGetGridComplexity()`, `PCHPDDM`, `PCHYPRE`, `PCGAMG` 472 @*/ 473 PetscErrorCode PCHPDDMGetComplexities(PC pc, PetscReal *gc, PetscReal *oc) 474 { 475 PC_HPDDM *data = (PC_HPDDM *)pc->data; 476 MatInfo info; 477 PetscLogDouble accumulate[2]{}, nnz1 = 1.0, m1 = 1.0; 478 479 PetscFunctionBegin; 480 if (gc) { 481 PetscAssertPointer(gc, 2); 482 *gc = 0; 483 } 484 if (oc) { 485 PetscAssertPointer(oc, 3); 486 *oc = 0; 487 } 488 for (PetscInt n = 0; n < data->N; ++n) { 489 if (data->levels[n]->ksp) { 490 Mat P, A = nullptr; 491 PetscInt m; 492 PetscBool flg = PETSC_FALSE; 493 494 PetscCall(KSPGetOperators(data->levels[n]->ksp, nullptr, &P)); 495 PetscCall(MatGetSize(P, &m, nullptr)); 496 accumulate[0] += m; 497 if (n == 0) { 498 PetscCall(PetscObjectTypeCompareAny((PetscObject)P, &flg, MATNORMAL, MATNORMALHERMITIAN, "")); 499 if (flg) { 500 PetscCall(MatConvert(P, MATAIJ, MAT_INITIAL_MATRIX, &A)); 501 P = A; 502 } else { 503 PetscCall(PetscObjectTypeCompare((PetscObject)P, MATSCHURCOMPLEMENT, &flg)); 504 PetscCall(PetscObjectReference((PetscObject)P)); 505 } 506 } 507 if (!A && flg) accumulate[1] += m * m; /* assumption that a MATSCHURCOMPLEMENT is dense if stored explicitly */ 508 else if (P->ops->getinfo) { 509 PetscCall(MatGetInfo(P, MAT_GLOBAL_SUM, &info)); 510 accumulate[1] += info.nz_used; 511 } 512 if (n == 0) { 513 m1 = m; 514 if (!A && flg) nnz1 = m * m; 515 else if (P->ops->getinfo) nnz1 = info.nz_used; 516 PetscCall(MatDestroy(&P)); 517 } 518 } 519 } 520 /* only process #0 has access to the full hierarchy by construction, so broadcast to ensure consistent outputs */ 521 PetscCallMPI(MPI_Bcast(accumulate, 2, MPIU_PETSCLOGDOUBLE, 0, PetscObjectComm((PetscObject)pc))); 522 if (gc) *gc = static_cast<PetscReal>(accumulate[0] / m1); 523 if (oc) *oc = static_cast<PetscReal>(accumulate[1] / nnz1); 524 PetscFunctionReturn(PETSC_SUCCESS); 525 } 526 527 static PetscErrorCode PCView_HPDDM(PC pc, PetscViewer viewer) 528 { 529 PC_HPDDM *data = (PC_HPDDM *)pc->data; 530 PetscViewer subviewer; 531 PetscViewerFormat format; 532 PetscSubcomm subcomm; 533 PetscReal oc, gc; 534 PetscInt tabs; 535 PetscMPIInt size, color, rank; 536 PetscBool flg; 537 const char *name; 538 539 PetscFunctionBegin; 540 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &flg)); 541 if (flg) { 542 PetscCall(PetscViewerASCIIPrintf(viewer, "level%s: %" PetscInt_FMT "\n", data->N > 1 ? "s" : "", data->N)); 543 PetscCall(PCHPDDMGetComplexities(pc, &gc, &oc)); 544 if (data->N > 1) { 545 if (!data->deflation) { 546 PetscCall(PetscViewerASCIIPrintf(viewer, "Neumann matrix attached? %s\n", PetscBools[PetscBool3ToBool(data->Neumann)])); 547 PetscCall(PetscViewerASCIIPrintf(viewer, "shared subdomain KSP between SLEPc and PETSc? %s\n", PetscBools[data->share])); 548 } else PetscCall(PetscViewerASCIIPrintf(viewer, "user-supplied deflation matrix\n")); 549 PetscCall(PetscViewerASCIIPrintf(viewer, "coarse correction: %s\n", PCHPDDMCoarseCorrectionTypes[data->correction])); 550 PetscCall(PetscViewerASCIIPrintf(viewer, "on process #0, value%s (+ threshold%s if available) for selecting deflation vectors:", data->N > 2 ? "s" : "", data->N > 2 ? "s" : "")); 551 PetscCall(PetscViewerASCIIGetTab(viewer, &tabs)); 552 PetscCall(PetscViewerASCIISetTab(viewer, 0)); 553 for (PetscInt i = 1; i < data->N; ++i) { 554 PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, data->levels[i - 1]->nu)); 555 if (data->levels[i - 1]->threshold > static_cast<PetscReal>(-0.1)) PetscCall(PetscViewerASCIIPrintf(viewer, " (%g)", (double)data->levels[i - 1]->threshold)); 556 } 557 PetscCall(PetscViewerASCIIPrintf(viewer, "\n")); 558 PetscCall(PetscViewerASCIISetTab(viewer, tabs)); 559 } 560 PetscCall(PetscViewerASCIIPrintf(viewer, "grid and operator complexities: %g %g\n", (double)gc, (double)oc)); 561 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size)); 562 if (data->levels[0]->ksp) { 563 PetscCall(KSPView(data->levels[0]->ksp, viewer)); 564 if (data->levels[0]->pc) PetscCall(PCView(data->levels[0]->pc, viewer)); 565 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank)); 566 for (PetscInt i = 1; i < data->N; ++i) { 567 if (data->levels[i]->ksp) color = 1; 568 else color = 0; 569 PetscCall(PetscSubcommCreate(PetscObjectComm((PetscObject)pc), &subcomm)); 570 PetscCall(PetscSubcommSetNumber(subcomm, PetscMin(size, 2))); 571 PetscCall(PetscSubcommSetTypeGeneral(subcomm, color, rank)); 572 PetscCall(PetscViewerASCIIPushTab(viewer)); 573 PetscCall(PetscViewerGetSubViewer(viewer, PetscSubcommChild(subcomm), &subviewer)); 574 if (color == 1) { 575 PetscCall(KSPView(data->levels[i]->ksp, subviewer)); 576 if (data->levels[i]->pc) PetscCall(PCView(data->levels[i]->pc, subviewer)); 577 PetscCall(PetscViewerFlush(subviewer)); 578 } 579 PetscCall(PetscViewerRestoreSubViewer(viewer, PetscSubcommChild(subcomm), &subviewer)); 580 PetscCall(PetscViewerASCIIPopTab(viewer)); 581 PetscCall(PetscSubcommDestroy(&subcomm)); 582 } 583 } 584 PetscCall(PetscViewerGetFormat(viewer, &format)); 585 if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 586 PetscCall(PetscViewerFileGetName(viewer, &name)); 587 if (name) { 588 Mat aux[2]; 589 IS is; 590 const PetscInt *indices; 591 PetscInt m, n, sizes[5] = {pc->mat->rmap->n, pc->mat->cmap->n, pc->mat->rmap->N, pc->mat->cmap->N, 0}; 592 char *tmp; 593 std::string prefix, suffix; 594 size_t pos; 595 596 PetscCall(PetscStrstr(name, ".", &tmp)); 597 if (tmp) { 598 pos = std::distance(const_cast<char *>(name), tmp); 599 prefix = std::string(name, pos); 600 suffix = std::string(name + pos + 1); 601 } else prefix = name; 602 if (data->aux) { 603 PetscCall(MatGetSize(data->aux, &m, &n)); 604 PetscCall(MatCreate(PetscObjectComm((PetscObject)pc), aux)); 605 PetscCall(MatSetSizes(aux[0], m, n, PETSC_DETERMINE, PETSC_DETERMINE)); 606 PetscCall(PetscObjectBaseTypeCompare((PetscObject)data->aux, MATSEQAIJ, &flg)); 607 if (flg) PetscCall(MatSetType(aux[0], MATMPIAIJ)); 608 else { 609 PetscCall(PetscObjectBaseTypeCompare((PetscObject)data->aux, MATSEQBAIJ, &flg)); 610 if (flg) PetscCall(MatSetType(aux[0], MATMPIBAIJ)); 611 else { 612 PetscCall(PetscObjectBaseTypeCompare((PetscObject)data->aux, MATSEQSBAIJ, &flg)); 613 PetscCheck(flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "MatType of auxiliary Mat (%s) is not any of the following: MATSEQAIJ, MATSEQBAIJ, or MATSEQSBAIJ", ((PetscObject)data->aux)->type_name); 614 PetscCall(MatSetType(aux[0], MATMPISBAIJ)); 615 } 616 } 617 PetscCall(MatSetBlockSizesFromMats(aux[0], data->aux, data->aux)); 618 PetscCall(MatAssemblyBegin(aux[0], MAT_FINAL_ASSEMBLY)); 619 PetscCall(MatAssemblyEnd(aux[0], MAT_FINAL_ASSEMBLY)); 620 PetscCall(MatGetDiagonalBlock(aux[0], aux + 1)); 621 PetscCall(MatCopy(data->aux, aux[1], DIFFERENT_NONZERO_PATTERN)); 622 PetscCall(PetscViewerBinaryOpen(PetscObjectComm((PetscObject)pc), std::string(prefix + "_aux_" + std::to_string(size) + (tmp ? ("." + suffix) : "")).c_str(), FILE_MODE_WRITE, &subviewer)); 623 PetscCall(MatView(aux[0], subviewer)); 624 PetscCall(PetscViewerDestroy(&subviewer)); 625 PetscCall(MatDestroy(aux)); 626 } 627 if (data->is) { 628 PetscCall(ISGetIndices(data->is, &indices)); 629 PetscCall(ISGetSize(data->is, sizes + 4)); 630 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pc), sizes[4], indices, PETSC_USE_POINTER, &is)); 631 PetscCall(PetscViewerBinaryOpen(PetscObjectComm((PetscObject)pc), std::string(prefix + "_is_" + std::to_string(size) + (tmp ? ("." + suffix) : "")).c_str(), FILE_MODE_WRITE, &subviewer)); 632 PetscCall(ISView(is, subviewer)); 633 PetscCall(PetscViewerDestroy(&subviewer)); 634 PetscCall(ISDestroy(&is)); 635 PetscCall(ISRestoreIndices(data->is, &indices)); 636 } 637 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pc), PETSC_STATIC_ARRAY_LENGTH(sizes), sizes, PETSC_USE_POINTER, &is)); 638 PetscCall(PetscViewerBinaryOpen(PetscObjectComm((PetscObject)pc), std::string(prefix + "_sizes_" + std::to_string(size) + (tmp ? ("." + suffix) : "")).c_str(), FILE_MODE_WRITE, &subviewer)); 639 PetscCall(ISView(is, subviewer)); 640 PetscCall(PetscViewerDestroy(&subviewer)); 641 PetscCall(ISDestroy(&is)); 642 } 643 } 644 } 645 PetscFunctionReturn(PETSC_SUCCESS); 646 } 647 648 static PetscErrorCode PCPreSolve_HPDDM(PC pc, KSP ksp, Vec, Vec) 649 { 650 PC_HPDDM *data = (PC_HPDDM *)pc->data; 651 Mat A; 652 PetscBool flg; 653 654 PetscFunctionBegin; 655 if (ksp) { 656 PetscCall(PetscObjectTypeCompare((PetscObject)ksp, KSPLSQR, &flg)); 657 if (flg && !data->normal) { 658 PetscCall(KSPGetOperators(ksp, &A, nullptr)); 659 PetscCall(MatCreateVecs(A, nullptr, &data->normal)); /* temporary Vec used in PCApply_HPDDMShell() for coarse grid corrections */ 660 } else if (!flg) { 661 PetscCall(PetscObjectTypeCompareAny((PetscObject)ksp, &flg, KSPCG, KSPGROPPCG, KSPPIPECG, KSPPIPECGRR, KSPPIPELCG, KSPPIPEPRCG, KSPPIPECG2, KSPSTCG, KSPFCG, KSPPIPEFCG, KSPMINRES, KSPNASH, KSPSYMMLQ, "")); 662 if (!flg) { 663 PetscCall(PetscObjectTypeCompare((PetscObject)ksp, KSPHPDDM, &flg)); 664 if (flg) { 665 KSPHPDDMType type; 666 667 PetscCall(KSPHPDDMGetType(ksp, &type)); 668 flg = (type == KSP_HPDDM_TYPE_CG || type == KSP_HPDDM_TYPE_BCG || type == KSP_HPDDM_TYPE_BFBCG ? PETSC_TRUE : PETSC_FALSE); 669 } 670 } 671 } 672 if (flg) { 673 if (data->correction == PC_HPDDM_COARSE_CORRECTION_DEFLATED) { 674 PetscCall(PetscOptionsHasName(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_hpddm_coarse_correction", &flg)); 675 PetscCheck(flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "PCHPDDMCoarseCorrectionType %s is known to be not symmetric, but KSPType %s requires a symmetric PC, if you insist on using this configuration, use the additional option -%spc_hpddm_coarse_correction %s, or alternatively, switch to a symmetric PCHPDDMCoarseCorrectionType such as %s", 676 PCHPDDMCoarseCorrectionTypes[data->correction], ((PetscObject)ksp)->type_name, ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", PCHPDDMCoarseCorrectionTypes[data->correction], PCHPDDMCoarseCorrectionTypes[PC_HPDDM_COARSE_CORRECTION_BALANCED]); 677 } 678 for (PetscInt n = 0; n < data->N; ++n) { 679 if (data->levels[n]->pc) { 680 PetscCall(PetscObjectTypeCompare((PetscObject)data->levels[n]->pc, PCASM, &flg)); 681 if (flg) { 682 PCASMType type; 683 684 PetscCall(PCASMGetType(data->levels[n]->pc, &type)); 685 if (type == PC_ASM_RESTRICT || type == PC_ASM_INTERPOLATE) { 686 PetscCall(PetscOptionsHasName(((PetscObject)data->levels[n]->pc)->options, ((PetscObject)data->levels[n]->pc)->prefix, "-pc_asm_type", &flg)); 687 PetscCheck(flg, PetscObjectComm((PetscObject)data->levels[n]->pc), PETSC_ERR_ARG_INCOMP, "PCASMType %s is known to be not symmetric, but KSPType %s requires a symmetric PC, if you insist on using this configuration, use the additional option -%spc_asm_type %s, or alternatively, switch to a symmetric PCASMType such as %s", PCASMTypes[type], 688 ((PetscObject)ksp)->type_name, ((PetscObject)data->levels[n]->pc)->prefix, PCASMTypes[type], PCASMTypes[PC_ASM_BASIC]); 689 } 690 } 691 } 692 } 693 } 694 } 695 PetscFunctionReturn(PETSC_SUCCESS); 696 } 697 698 static PetscErrorCode PCSetUp_HPDDMShell(PC pc) 699 { 700 PC_HPDDM_Level *ctx; 701 Mat A, P; 702 Vec x; 703 const char *pcpre; 704 705 PetscFunctionBegin; 706 PetscCall(PCShellGetContext(pc, &ctx)); 707 PetscCall(KSPGetOptionsPrefix(ctx->ksp, &pcpre)); 708 PetscCall(KSPGetOperators(ctx->ksp, &A, &P)); 709 /* smoother */ 710 PetscCall(PCSetOptionsPrefix(ctx->pc, pcpre)); 711 PetscCall(PCSetOperators(ctx->pc, A, P)); 712 if (!ctx->v[0]) { 713 PetscCall(VecDuplicateVecs(ctx->D, 1, &ctx->v[0])); 714 if (!std::is_same<PetscScalar, PetscReal>::value) PetscCall(VecDestroy(&ctx->D)); 715 PetscCall(MatCreateVecs(A, &x, nullptr)); 716 PetscCall(VecDuplicateVecs(x, 2, &ctx->v[1])); 717 PetscCall(VecDestroy(&x)); 718 } 719 std::fill_n(ctx->V, 3, nullptr); 720 PetscFunctionReturn(PETSC_SUCCESS); 721 } 722 723 template <class Type, typename std::enable_if<std::is_same<Type, Vec>::value>::type * = nullptr> 724 static inline PetscErrorCode PCHPDDMDeflate_Private(PC pc, Type x, Type y) 725 { 726 PC_HPDDM_Level *ctx; 727 PetscScalar *out; 728 729 PetscFunctionBegin; 730 PetscCall(PCShellGetContext(pc, &ctx)); 731 /* going from PETSc to HPDDM numbering */ 732 PetscCall(VecScatterBegin(ctx->scatter, x, ctx->v[0][0], INSERT_VALUES, SCATTER_FORWARD)); 733 PetscCall(VecScatterEnd(ctx->scatter, x, ctx->v[0][0], INSERT_VALUES, SCATTER_FORWARD)); 734 PetscCall(VecGetArrayWrite(ctx->v[0][0], &out)); 735 ctx->P->deflation<false>(nullptr, out, 1); /* y = Q x */ 736 PetscCall(VecRestoreArrayWrite(ctx->v[0][0], &out)); 737 /* going from HPDDM to PETSc numbering */ 738 PetscCall(VecScatterBegin(ctx->scatter, ctx->v[0][0], y, INSERT_VALUES, SCATTER_REVERSE)); 739 PetscCall(VecScatterEnd(ctx->scatter, ctx->v[0][0], y, INSERT_VALUES, SCATTER_REVERSE)); 740 PetscFunctionReturn(PETSC_SUCCESS); 741 } 742 743 template <class Type, typename std::enable_if<std::is_same<Type, Mat>::value>::type * = nullptr> 744 static inline PetscErrorCode PCHPDDMDeflate_Private(PC pc, Type X, Type Y) 745 { 746 PC_HPDDM_Level *ctx; 747 Vec vX, vY, vC; 748 PetscScalar *out; 749 PetscInt N; 750 751 PetscFunctionBegin; 752 PetscCall(PCShellGetContext(pc, &ctx)); 753 PetscCall(MatGetSize(X, nullptr, &N)); 754 /* going from PETSc to HPDDM numbering */ 755 for (PetscInt i = 0; i < N; ++i) { 756 PetscCall(MatDenseGetColumnVecRead(X, i, &vX)); 757 PetscCall(MatDenseGetColumnVecWrite(ctx->V[0], i, &vC)); 758 PetscCall(VecScatterBegin(ctx->scatter, vX, vC, INSERT_VALUES, SCATTER_FORWARD)); 759 PetscCall(VecScatterEnd(ctx->scatter, vX, vC, INSERT_VALUES, SCATTER_FORWARD)); 760 PetscCall(MatDenseRestoreColumnVecWrite(ctx->V[0], i, &vC)); 761 PetscCall(MatDenseRestoreColumnVecRead(X, i, &vX)); 762 } 763 PetscCall(MatDenseGetArrayWrite(ctx->V[0], &out)); 764 ctx->P->deflation<false>(nullptr, out, N); /* Y = Q X */ 765 PetscCall(MatDenseRestoreArrayWrite(ctx->V[0], &out)); 766 /* going from HPDDM to PETSc numbering */ 767 for (PetscInt i = 0; i < N; ++i) { 768 PetscCall(MatDenseGetColumnVecRead(ctx->V[0], i, &vC)); 769 PetscCall(MatDenseGetColumnVecWrite(Y, i, &vY)); 770 PetscCall(VecScatterBegin(ctx->scatter, vC, vY, INSERT_VALUES, SCATTER_REVERSE)); 771 PetscCall(VecScatterEnd(ctx->scatter, vC, vY, INSERT_VALUES, SCATTER_REVERSE)); 772 PetscCall(MatDenseRestoreColumnVecWrite(Y, i, &vY)); 773 PetscCall(MatDenseRestoreColumnVecRead(ctx->V[0], i, &vC)); 774 } 775 PetscFunctionReturn(PETSC_SUCCESS); 776 } 777 778 /* 779 PCApply_HPDDMShell - Applies a (2) deflated, (1) additive, (3) balanced, or (4) no coarse correction. In what follows, E = Z Pmat Z^T and Q = Z^T E^-1 Z. 780 781 .vb 782 (1) y = Pmat^-1 x + Q x, 783 (2) y = Pmat^-1 (I - Amat Q) x + Q x (default), 784 (3) y = (I - Q Amat^T) Pmat^-1 (I - Amat Q) x + Q x, 785 (4) y = Pmat^-1 x . 786 .ve 787 788 Input Parameters: 789 + pc - preconditioner context 790 - x - input vector 791 792 Output Parameter: 793 . y - output vector 794 795 Notes: 796 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. 797 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. 798 (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. 799 800 Level: advanced 801 802 Developer Note: 803 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 804 to trigger it. Likely the manual page is `PCHPDDM` 805 806 .seealso: [](ch_ksp), `PCHPDDM`, `PCHPDDMCoarseCorrectionType` 807 */ 808 static PetscErrorCode PCApply_HPDDMShell(PC pc, Vec x, Vec y) 809 { 810 PC_HPDDM_Level *ctx; 811 Mat A; 812 813 PetscFunctionBegin; 814 PetscCall(PCShellGetContext(pc, &ctx)); 815 PetscCheck(ctx->P, PETSC_COMM_SELF, PETSC_ERR_PLIB, "PCSHELL from PCHPDDM called with no HPDDM object"); 816 PetscCall(KSPGetOperators(ctx->ksp, &A, nullptr)); 817 if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_NONE) PetscCall(PCApply(ctx->pc, x, y)); /* y = M^-1 x */ 818 else { 819 PetscCall(PCHPDDMDeflate_Private(pc, x, y)); /* y = Q x */ 820 if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_DEFLATED || ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_BALANCED) { 821 if (!ctx->parent->normal || ctx != ctx->parent->levels[0]) PetscCall(MatMult(A, y, ctx->v[1][0])); /* y = A Q x */ 822 else { /* KSPLSQR and finest level */ PetscCall(MatMult(A, y, ctx->parent->normal)); /* y = A Q x */ 823 PetscCall(MatMultHermitianTranspose(A, ctx->parent->normal, ctx->v[1][0])); /* y = A^T A Q x */ 824 } 825 PetscCall(VecWAXPY(ctx->v[1][1], -1.0, ctx->v[1][0], x)); /* y = (I - A Q) x */ 826 PetscCall(PCApply(ctx->pc, ctx->v[1][1], ctx->v[1][0])); /* y = M^-1 (I - A Q) x */ 827 if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_BALANCED) { 828 if (!ctx->parent->normal || ctx != ctx->parent->levels[0]) PetscCall(MatMultHermitianTranspose(A, ctx->v[1][0], ctx->v[1][1])); /* z = A^T y */ 829 else { 830 PetscCall(MatMult(A, ctx->v[1][0], ctx->parent->normal)); 831 PetscCall(MatMultHermitianTranspose(A, ctx->parent->normal, ctx->v[1][1])); /* z = A^T A y */ 832 } 833 PetscCall(PCHPDDMDeflate_Private(pc, ctx->v[1][1], ctx->v[1][1])); 834 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 */ 835 } else PetscCall(VecAXPY(y, 1.0, ctx->v[1][0])); /* y = Q M^-1 (I - A Q) x + Q x */ 836 } else { 837 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); 838 PetscCall(PCApply(ctx->pc, x, ctx->v[1][0])); 839 PetscCall(VecAXPY(y, 1.0, ctx->v[1][0])); /* y = M^-1 x + Q x */ 840 } 841 } 842 PetscFunctionReturn(PETSC_SUCCESS); 843 } 844 845 /* 846 PCMatApply_HPDDMShell - Variant of PCApply_HPDDMShell() for blocks of vectors. 847 848 Input Parameters: 849 + pc - preconditioner context 850 - X - block of input vectors 851 852 Output Parameter: 853 . Y - block of output vectors 854 855 Level: advanced 856 857 .seealso: [](ch_ksp), `PCHPDDM`, `PCApply_HPDDMShell()`, `PCHPDDMCoarseCorrectionType` 858 */ 859 static PetscErrorCode PCMatApply_HPDDMShell(PC pc, Mat X, Mat Y) 860 { 861 PC_HPDDM_Level *ctx; 862 Mat A, *ptr; 863 PetscContainer container = nullptr; 864 PetscScalar *array; 865 PetscInt m, M, N, prev = 0; 866 PetscBool reset = PETSC_FALSE; 867 868 PetscFunctionBegin; 869 PetscCall(PCShellGetContext(pc, &ctx)); 870 PetscCheck(ctx->P, PETSC_COMM_SELF, PETSC_ERR_PLIB, "PCSHELL from PCHPDDM called with no HPDDM object"); 871 PetscCall(MatGetSize(X, nullptr, &N)); 872 PetscCall(KSPGetOperators(ctx->ksp, &A, nullptr)); 873 if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_NONE) PetscCall(PCMatApply(ctx->pc, X, Y)); 874 else { 875 PetscCall(PetscObjectQuery((PetscObject)A, "_HPDDM_MatProduct", (PetscObject *)&container)); 876 if (container) { /* MatProduct container already attached */ 877 PetscCall(PetscContainerGetPointer(container, (void **)&ptr)); 878 if (ptr[1] != ctx->V[2]) /* Mat has changed or may have been set first in KSPHPDDM */ 879 for (m = 0; m < 2; ++m) { 880 PetscCall(MatDestroy(ctx->V + m + 1)); 881 ctx->V[m + 1] = ptr[m]; 882 PetscCall(PetscObjectReference((PetscObject)ctx->V[m + 1])); 883 } 884 } 885 if (ctx->V[1]) PetscCall(MatGetSize(ctx->V[1], nullptr, &prev)); 886 if (N != prev || !ctx->V[0]) { 887 PetscCall(MatDestroy(ctx->V)); 888 PetscCall(VecGetLocalSize(ctx->v[0][0], &m)); 889 PetscCall(MatCreateDense(PetscObjectComm((PetscObject)pc), m, PETSC_DECIDE, PETSC_DECIDE, N, nullptr, ctx->V)); 890 if (N != prev) { 891 PetscCall(MatDestroy(ctx->V + 1)); 892 PetscCall(MatDestroy(ctx->V + 2)); 893 PetscCall(MatGetLocalSize(X, &m, nullptr)); 894 PetscCall(MatGetSize(X, &M, nullptr)); 895 if (ctx->parent->correction != PC_HPDDM_COARSE_CORRECTION_BALANCED) PetscCall(MatDenseGetArrayWrite(ctx->V[0], &array)); 896 else array = nullptr; 897 PetscCall(MatCreateDense(PetscObjectComm((PetscObject)pc), m, PETSC_DECIDE, M, N, array, ctx->V + 1)); 898 if (ctx->parent->correction != PC_HPDDM_COARSE_CORRECTION_BALANCED) PetscCall(MatDenseRestoreArrayWrite(ctx->V[0], &array)); 899 PetscCall(MatCreateDense(PetscObjectComm((PetscObject)pc), m, PETSC_DECIDE, M, N, nullptr, ctx->V + 2)); 900 PetscCall(MatProductCreateWithMat(A, Y, nullptr, ctx->V[1])); 901 PetscCall(MatProductSetType(ctx->V[1], MATPRODUCT_AB)); 902 PetscCall(MatProductSetFromOptions(ctx->V[1])); 903 PetscCall(MatProductSymbolic(ctx->V[1])); 904 if (!container) PetscCall(PetscObjectContainerCompose((PetscObject)A, "_HPDDM_MatProduct", ctx->V + 1, nullptr)); /* no MatProduct container attached, create one to be queried in KSPHPDDM or at the next call to PCMatApply() */ 905 else PetscCall(PetscContainerSetPointer(container, ctx->V + 1)); /* need to compose B and D from MatProductCreateWithMat(A, B, NULL, D), which are stored in the contiguous array ctx->V */ 906 } 907 if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_BALANCED) { 908 PetscCall(MatProductCreateWithMat(A, ctx->V[1], nullptr, ctx->V[2])); 909 PetscCall(MatProductSetType(ctx->V[2], MATPRODUCT_AtB)); 910 PetscCall(MatProductSetFromOptions(ctx->V[2])); 911 PetscCall(MatProductSymbolic(ctx->V[2])); 912 } 913 ctx->P->start(N); 914 } 915 if (N == prev || container) { /* when MatProduct container is attached, always need to MatProductReplaceMats() since KSPHPDDM may have replaced the Mat as well */ 916 PetscCall(MatProductReplaceMats(nullptr, Y, nullptr, ctx->V[1])); 917 if (container && ctx->parent->correction != PC_HPDDM_COARSE_CORRECTION_BALANCED) { 918 PetscCall(MatDenseGetArrayWrite(ctx->V[0], &array)); 919 PetscCall(MatDensePlaceArray(ctx->V[1], array)); 920 PetscCall(MatDenseRestoreArrayWrite(ctx->V[0], &array)); 921 reset = PETSC_TRUE; 922 } 923 } 924 PetscCall(PCHPDDMDeflate_Private(pc, X, Y)); 925 if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_DEFLATED || ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_BALANCED) { 926 PetscCall(MatProductNumeric(ctx->V[1])); 927 PetscCall(MatCopy(ctx->V[1], ctx->V[2], SAME_NONZERO_PATTERN)); 928 PetscCall(MatAXPY(ctx->V[2], -1.0, X, SAME_NONZERO_PATTERN)); 929 PetscCall(PCMatApply(ctx->pc, ctx->V[2], ctx->V[1])); 930 if (ctx->parent->correction == PC_HPDDM_COARSE_CORRECTION_BALANCED) { 931 PetscCall(MatProductNumeric(ctx->V[2])); 932 PetscCall(PCHPDDMDeflate_Private(pc, ctx->V[2], ctx->V[2])); 933 PetscCall(MatAXPY(ctx->V[1], -1.0, ctx->V[2], SAME_NONZERO_PATTERN)); 934 } 935 PetscCall(MatAXPY(Y, -1.0, ctx->V[1], SAME_NONZERO_PATTERN)); 936 } else { 937 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); 938 PetscCall(PCMatApply(ctx->pc, X, ctx->V[1])); 939 PetscCall(MatAXPY(Y, 1.0, ctx->V[1], SAME_NONZERO_PATTERN)); 940 } 941 if (reset) PetscCall(MatDenseResetArray(ctx->V[1])); 942 } 943 PetscFunctionReturn(PETSC_SUCCESS); 944 } 945 946 static PetscErrorCode PCDestroy_HPDDMShell(PC pc) 947 { 948 PC_HPDDM_Level *ctx; 949 950 PetscFunctionBegin; 951 PetscCall(PCShellGetContext(pc, &ctx)); 952 PetscCall(HPDDM::Schwarz<PetscScalar>::destroy(ctx, PETSC_TRUE)); 953 PetscCall(VecDestroyVecs(1, &ctx->v[0])); 954 PetscCall(VecDestroyVecs(2, &ctx->v[1])); 955 PetscCall(PetscObjectCompose((PetscObject)ctx->pc->mat, "_HPDDM_MatProduct", nullptr)); 956 PetscCall(MatDestroy(ctx->V)); 957 PetscCall(MatDestroy(ctx->V + 1)); 958 PetscCall(MatDestroy(ctx->V + 2)); 959 PetscCall(VecDestroy(&ctx->D)); 960 PetscCall(PetscSFDestroy(&ctx->scatter)); 961 PetscCall(PCDestroy(&ctx->pc)); 962 PetscFunctionReturn(PETSC_SUCCESS); 963 } 964 965 template <class Type, bool T = false, typename std::enable_if<std::is_same<Type, Vec>::value>::type * = nullptr> 966 static inline PetscErrorCode PCApply_Schur_Private(std::tuple<KSP, IS, Vec[2]> *p, PC factor, Type x, Type y) 967 { 968 PetscFunctionBegin; 969 PetscCall(VecISCopy(std::get<2>(*p)[0], std::get<1>(*p), SCATTER_FORWARD, x)); 970 if (!T) PetscCall(PCApply(factor, std::get<2>(*p)[0], std::get<2>(*p)[1])); 971 else PetscCall(PCApplyTranspose(factor, std::get<2>(*p)[0], std::get<2>(*p)[1])); 972 PetscCall(VecISCopy(std::get<2>(*p)[1], std::get<1>(*p), SCATTER_REVERSE, y)); 973 PetscFunctionReturn(PETSC_SUCCESS); 974 } 975 976 template <class Type, bool = false, typename std::enable_if<std::is_same<Type, Mat>::value>::type * = nullptr> 977 static inline PetscErrorCode PCApply_Schur_Private(std::tuple<KSP, IS, Vec[2]> *p, PC factor, Type X, Type Y) 978 { 979 Mat B[2]; 980 Vec x, y; 981 982 PetscFunctionBegin; 983 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, factor->mat->rmap->n, X->cmap->n, nullptr, B)); 984 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, factor->mat->rmap->n, X->cmap->n, nullptr, B + 1)); 985 for (PetscInt i = 0; i < X->cmap->n; ++i) { 986 PetscCall(MatDenseGetColumnVecRead(X, i, &x)); 987 PetscCall(MatDenseGetColumnVecWrite(B[0], i, &y)); 988 PetscCall(VecISCopy(y, std::get<1>(*p), SCATTER_FORWARD, x)); 989 PetscCall(MatDenseRestoreColumnVecWrite(B[0], i, &y)); 990 PetscCall(MatDenseRestoreColumnVecRead(X, i, &x)); 991 } 992 PetscCall(PCMatApply(factor, B[0], B[1])); 993 PetscCall(MatDestroy(B)); 994 for (PetscInt i = 0; i < X->cmap->n; ++i) { 995 PetscCall(MatDenseGetColumnVecRead(B[1], i, &x)); 996 PetscCall(MatDenseGetColumnVecWrite(Y, i, &y)); 997 PetscCall(VecISCopy(x, std::get<1>(*p), SCATTER_REVERSE, y)); 998 PetscCall(MatDenseRestoreColumnVecWrite(Y, i, &y)); 999 PetscCall(MatDenseRestoreColumnVecRead(B[1], i, &x)); 1000 } 1001 PetscCall(MatDestroy(B + 1)); 1002 PetscFunctionReturn(PETSC_SUCCESS); 1003 } 1004 1005 template <class Type = Vec, bool T = false> 1006 static PetscErrorCode PCApply_Schur(PC pc, Type x, Type y) 1007 { 1008 PC factor; 1009 Mat A; 1010 MatSolverType type; 1011 PetscBool flg; 1012 std::tuple<KSP, IS, Vec[2]> *p; 1013 1014 PetscFunctionBegin; 1015 PetscCall(PCShellGetContext(pc, &p)); 1016 PetscCall(KSPGetPC(std::get<0>(*p), &factor)); 1017 PetscCall(PCFactorGetMatSolverType(factor, &type)); 1018 PetscCall(PCFactorGetMatrix(factor, &A)); 1019 PetscCall(PetscStrcmp(type, MATSOLVERMUMPS, &flg)); 1020 if (flg) { 1021 PetscCheck(PetscDefined(HAVE_MUMPS), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent MatSolverType"); 1022 #if PetscDefined(HAVE_MUMPS) 1023 PetscCall(MatMumpsSetIcntl(A, 26, 0)); 1024 #endif 1025 } else { 1026 PetscCall(PetscStrcmp(type, MATSOLVERMKL_PARDISO, &flg)); 1027 PetscCheck(flg && PetscDefined(HAVE_MKL_PARDISO), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent MatSolverType"); 1028 flg = PETSC_FALSE; 1029 #if PetscDefined(HAVE_MKL_PARDISO) 1030 PetscCall(MatMkl_PardisoSetCntl(A, 70, 1)); 1031 #endif 1032 } 1033 PetscCall(PCApply_Schur_Private<Type, T>(p, factor, x, y)); 1034 if (flg) { 1035 #if PetscDefined(HAVE_MUMPS) 1036 PetscCall(MatMumpsSetIcntl(A, 26, -1)); 1037 #endif 1038 } else { 1039 #if PetscDefined(HAVE_MKL_PARDISO) 1040 PetscCall(MatMkl_PardisoSetCntl(A, 70, 0)); 1041 #endif 1042 } 1043 PetscFunctionReturn(PETSC_SUCCESS); 1044 } 1045 1046 static PetscErrorCode PCDestroy_Schur(PC pc) 1047 { 1048 std::tuple<KSP, IS, Vec[2]> *p; 1049 1050 PetscFunctionBegin; 1051 PetscCall(PCShellGetContext(pc, &p)); 1052 PetscCall(ISDestroy(&std::get<1>(*p))); 1053 PetscCall(VecDestroy(std::get<2>(*p))); 1054 PetscCall(VecDestroy(std::get<2>(*p) + 1)); 1055 PetscCall(PetscFree(p)); 1056 PetscFunctionReturn(PETSC_SUCCESS); 1057 } 1058 1059 static PetscErrorCode PCHPDDMSolve_Private(const PC_HPDDM_Level *ctx, PetscScalar *rhs, const unsigned short &mu) 1060 { 1061 Mat B, X; 1062 PetscInt n, N, j = 0; 1063 1064 PetscFunctionBegin; 1065 PetscCall(KSPGetOperators(ctx->ksp, &B, nullptr)); 1066 PetscCall(MatGetLocalSize(B, &n, nullptr)); 1067 PetscCall(MatGetSize(B, &N, nullptr)); 1068 if (ctx->parent->log_separate) { 1069 j = std::distance(ctx->parent->levels, std::find(ctx->parent->levels, ctx->parent->levels + ctx->parent->N, ctx)); 1070 PetscCall(PetscLogEventBegin(PC_HPDDM_Solve[j], ctx->ksp, nullptr, nullptr, nullptr)); 1071 } 1072 if (mu == 1) { 1073 if (!ctx->ksp->vec_rhs) { 1074 PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)ctx->ksp), 1, n, N, nullptr, &ctx->ksp->vec_rhs)); 1075 PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)ctx->ksp), n, N, &ctx->ksp->vec_sol)); 1076 } 1077 PetscCall(VecPlaceArray(ctx->ksp->vec_rhs, rhs)); 1078 PetscCall(KSPSolve(ctx->ksp, nullptr, nullptr)); 1079 PetscCall(VecCopy(ctx->ksp->vec_sol, ctx->ksp->vec_rhs)); 1080 PetscCall(VecResetArray(ctx->ksp->vec_rhs)); 1081 } else { 1082 PetscCall(MatCreateDense(PetscObjectComm((PetscObject)ctx->ksp), n, PETSC_DECIDE, N, mu, rhs, &B)); 1083 PetscCall(MatCreateDense(PetscObjectComm((PetscObject)ctx->ksp), n, PETSC_DECIDE, N, mu, nullptr, &X)); 1084 PetscCall(KSPMatSolve(ctx->ksp, B, X)); 1085 PetscCall(MatCopy(X, B, SAME_NONZERO_PATTERN)); 1086 PetscCall(MatDestroy(&X)); 1087 PetscCall(MatDestroy(&B)); 1088 } 1089 if (ctx->parent->log_separate) PetscCall(PetscLogEventEnd(PC_HPDDM_Solve[j], ctx->ksp, nullptr, nullptr, nullptr)); 1090 PetscFunctionReturn(PETSC_SUCCESS); 1091 } 1092 1093 static PetscErrorCode PCHPDDMSetUpNeumannOverlap_Private(PC pc) 1094 { 1095 PC_HPDDM *data = (PC_HPDDM *)pc->data; 1096 1097 PetscFunctionBegin; 1098 if (data->setup) { 1099 Mat P; 1100 Vec x, xt = nullptr; 1101 PetscReal t = 0.0, s = 0.0; 1102 1103 PetscCall(PCGetOperators(pc, nullptr, &P)); 1104 PetscCall(PetscObjectQuery((PetscObject)P, "__SNES_latest_X", (PetscObject *)&x)); 1105 PetscCallBack("PCHPDDM Neumann callback", (*data->setup)(data->aux, t, x, xt, s, data->is, data->setup_ctx)); 1106 } 1107 PetscFunctionReturn(PETSC_SUCCESS); 1108 } 1109 1110 static PetscErrorCode PCHPDDMCreateSubMatrices_Private(Mat mat, PetscInt n, const IS *, const IS *, MatReuse scall, Mat *submat[]) 1111 { 1112 Mat A; 1113 PetscBool flg; 1114 1115 PetscFunctionBegin; 1116 PetscCheck(n == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "MatCreateSubMatrices() called to extract %" PetscInt_FMT " submatrices, which is different than 1", n); 1117 /* previously composed Mat */ 1118 PetscCall(PetscObjectQuery((PetscObject)mat, "_PCHPDDM_SubMatrices", (PetscObject *)&A)); 1119 PetscCheck(A, PETSC_COMM_SELF, PETSC_ERR_PLIB, "SubMatrices not found in Mat"); 1120 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATSCHURCOMPLEMENT, &flg)); /* MATSCHURCOMPLEMENT has neither a MatDuplicate() nor a MatCopy() implementation */ 1121 if (scall == MAT_INITIAL_MATRIX) { 1122 PetscCall(PetscCalloc1(2, submat)); /* allocate an extra Mat to avoid errors in MatDestroySubMatrices_Dummy() */ 1123 if (!flg) PetscCall(MatDuplicate(A, MAT_COPY_VALUES, *submat)); 1124 } else if (!flg) PetscCall(MatCopy(A, (*submat)[0], SAME_NONZERO_PATTERN)); 1125 if (flg) { 1126 PetscCall(MatDestroy(*submat)); /* previously created Mat has to be destroyed */ 1127 (*submat)[0] = A; 1128 PetscCall(PetscObjectReference((PetscObject)A)); 1129 } 1130 PetscFunctionReturn(PETSC_SUCCESS); 1131 } 1132 1133 static PetscErrorCode PCHPDDMCommunicationAvoidingPCASM_Private(PC pc, Mat C, PetscBool sorted) 1134 { 1135 void (*op)(void); 1136 1137 PetscFunctionBegin; 1138 /* previously-composed Mat */ 1139 PetscCall(PetscObjectCompose((PetscObject)pc->pmat, "_PCHPDDM_SubMatrices", (PetscObject)C)); 1140 PetscCall(MatGetOperation(pc->pmat, MATOP_CREATE_SUBMATRICES, &op)); 1141 /* trick suggested by Barry https://lists.mcs.anl.gov/pipermail/petsc-dev/2020-January/025491.html */ 1142 PetscCall(MatSetOperation(pc->pmat, MATOP_CREATE_SUBMATRICES, (void (*)(void))PCHPDDMCreateSubMatrices_Private)); 1143 if (sorted) PetscCall(PCASMSetSortIndices(pc, PETSC_FALSE)); /* everything is already sorted */ 1144 PetscCall(PCSetFromOptions(pc)); /* otherwise -pc_hpddm_levels_1_pc_asm_sub_mat_type is not used */ 1145 PetscCall(PCSetUp(pc)); 1146 /* reset MatCreateSubMatrices() */ 1147 PetscCall(MatSetOperation(pc->pmat, MATOP_CREATE_SUBMATRICES, op)); 1148 PetscCall(PetscObjectCompose((PetscObject)pc->pmat, "_PCHPDDM_SubMatrices", nullptr)); 1149 PetscFunctionReturn(PETSC_SUCCESS); 1150 } 1151 1152 static PetscErrorCode PCHPDDMPermute_Private(IS is, IS in_is, IS *out_is, Mat in_C, Mat *out_C, IS *p) 1153 { 1154 IS perm; 1155 const PetscInt *ptr; 1156 PetscInt *concatenate, size, bs; 1157 std::map<PetscInt, PetscInt> order; 1158 PetscBool sorted; 1159 1160 PetscFunctionBegin; 1161 PetscValidHeaderSpecific(is, IS_CLASSID, 1); 1162 PetscValidHeaderSpecific(in_C, MAT_CLASSID, 4); 1163 PetscCall(ISSorted(is, &sorted)); 1164 if (!sorted) { 1165 PetscCall(ISGetLocalSize(is, &size)); 1166 PetscCall(ISGetIndices(is, &ptr)); 1167 PetscCall(ISGetBlockSize(is, &bs)); 1168 /* MatCreateSubMatrices(), called by PCASM, follows the global numbering of Pmat */ 1169 for (PetscInt n = 0; n < size; n += bs) order.insert(std::make_pair(ptr[n] / bs, n / bs)); 1170 PetscCall(ISRestoreIndices(is, &ptr)); 1171 size /= bs; 1172 if (out_C) { 1173 PetscCall(PetscMalloc1(size, &concatenate)); 1174 for (const std::pair<const PetscInt, PetscInt> &i : order) *concatenate++ = i.second; 1175 concatenate -= size; 1176 PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)in_C), bs, size, concatenate, PETSC_OWN_POINTER, &perm)); 1177 PetscCall(ISSetPermutation(perm)); 1178 /* permute user-provided Mat so that it matches with MatCreateSubMatrices() numbering */ 1179 PetscCall(MatPermute(in_C, perm, perm, out_C)); 1180 if (p) *p = perm; 1181 else PetscCall(ISDestroy(&perm)); /* no need to save the permutation */ 1182 } 1183 if (out_is) { 1184 PetscCall(PetscMalloc1(size, &concatenate)); 1185 for (const std::pair<const PetscInt, PetscInt> &i : order) *concatenate++ = i.first; 1186 concatenate -= size; 1187 /* permute user-provided IS so that it matches with MatCreateSubMatrices() numbering */ 1188 PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)in_is), bs, size, concatenate, PETSC_OWN_POINTER, out_is)); 1189 } 1190 } else { /* input IS is sorted, nothing to permute, simply duplicate inputs when needed */ 1191 if (out_C) PetscCall(MatDuplicate(in_C, MAT_COPY_VALUES, out_C)); 1192 if (out_is) PetscCall(ISDuplicate(in_is, out_is)); 1193 } 1194 PetscFunctionReturn(PETSC_SUCCESS); 1195 } 1196 1197 static PetscErrorCode PCHPDDMCheckSymmetry_Private(PC pc, Mat A01, Mat A10, Mat *B01 = nullptr) 1198 { 1199 Mat T, U = nullptr, B = nullptr; 1200 IS z; 1201 PetscBool flg, conjugate = PETSC_FALSE; 1202 1203 PetscFunctionBegin; 1204 PetscCall(PetscObjectTypeCompare((PetscObject)A10, MATTRANSPOSEVIRTUAL, &flg)); 1205 if (B01) *B01 = nullptr; 1206 if (flg) { 1207 PetscCall(MatShellGetScalingShifts(A10, (PetscScalar *)MAT_SHELL_NOT_ALLOWED, (PetscScalar *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Mat *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED)); 1208 PetscCall(MatTransposeGetMat(A10, &U)); 1209 } else { 1210 PetscCall(PetscObjectTypeCompare((PetscObject)A10, MATHERMITIANTRANSPOSEVIRTUAL, &flg)); 1211 if (flg) { 1212 PetscCall(MatShellGetScalingShifts(A10, (PetscScalar *)MAT_SHELL_NOT_ALLOWED, (PetscScalar *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Mat *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED)); 1213 PetscCall(MatHermitianTransposeGetMat(A10, &U)); 1214 conjugate = PETSC_TRUE; 1215 } 1216 } 1217 if (U) PetscCall(MatDuplicate(U, MAT_COPY_VALUES, &T)); 1218 else PetscCall(MatHermitianTranspose(A10, MAT_INITIAL_MATRIX, &T)); 1219 PetscCall(PetscObjectTypeCompare((PetscObject)A01, MATTRANSPOSEVIRTUAL, &flg)); 1220 if (flg) { 1221 PetscCall(MatShellGetScalingShifts(A01, (PetscScalar *)MAT_SHELL_NOT_ALLOWED, (PetscScalar *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Mat *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED)); 1222 PetscCall(MatTransposeGetMat(A01, &A01)); 1223 PetscCall(MatTranspose(A01, MAT_INITIAL_MATRIX, &B)); 1224 A01 = B; 1225 } else { 1226 PetscCall(PetscObjectTypeCompare((PetscObject)A01, MATHERMITIANTRANSPOSEVIRTUAL, &flg)); 1227 if (flg) { 1228 PetscCall(MatShellGetScalingShifts(A01, (PetscScalar *)MAT_SHELL_NOT_ALLOWED, (PetscScalar *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Mat *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED)); 1229 PetscCall(MatHermitianTransposeGetMat(A01, &A01)); 1230 PetscCall(MatHermitianTranspose(A01, MAT_INITIAL_MATRIX, &B)); 1231 A01 = B; 1232 } 1233 } 1234 PetscCall(PetscLayoutCompare(T->rmap, A01->rmap, &flg)); 1235 if (flg) { 1236 PetscCall(PetscLayoutCompare(T->cmap, A01->cmap, &flg)); 1237 if (flg) { 1238 PetscCall(MatFindZeroRows(A01, &z)); /* for essential boundary conditions, some implementations will */ 1239 if (z) { /* zero rows in [P00 A01] except for the diagonal of P00 */ 1240 if (B01) PetscCall(MatDuplicate(T, MAT_COPY_VALUES, B01)); 1241 PetscCall(MatSetOption(T, MAT_NO_OFF_PROC_ZERO_ROWS, PETSC_TRUE)); 1242 PetscCall(MatZeroRowsIS(T, z, 0.0, nullptr, nullptr)); /* corresponding zero rows from A01 */ 1243 } 1244 PetscCall(MatMultEqual(A01, T, 20, &flg)); 1245 if (!B01) PetscCheck(flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "A01 != A10^T"); 1246 else { 1247 PetscCall(PetscInfo(pc, "A01 and A10^T are equal? %s\n", PetscBools[flg])); 1248 if (!flg) { 1249 if (z) PetscCall(MatDestroy(&T)); 1250 else *B01 = T; 1251 flg = PETSC_TRUE; 1252 } else PetscCall(MatDestroy(B01)); 1253 } 1254 PetscCall(ISDestroy(&z)); 1255 } 1256 } 1257 if (!flg) PetscCall(PetscInfo(pc, "A01 and A10^T have non-congruent layouts, cannot test for equality\n")); 1258 if (!B01 || !*B01) PetscCall(MatDestroy(&T)); 1259 else if (conjugate) PetscCall(MatConjugate(T)); 1260 PetscCall(MatDestroy(&B)); 1261 PetscFunctionReturn(PETSC_SUCCESS); 1262 } 1263 1264 static PetscErrorCode PCHPDDMCheckInclusion_Private(PC pc, IS is, IS is_local, PetscBool check) 1265 { 1266 IS intersect; 1267 const char *str = "IS of the auxiliary Mat does not include all local rows of A"; 1268 PetscBool equal; 1269 1270 PetscFunctionBegin; 1271 PetscCall(ISIntersect(is, is_local, &intersect)); 1272 PetscCall(ISEqualUnsorted(is_local, intersect, &equal)); 1273 PetscCall(ISDestroy(&intersect)); 1274 if (check) PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "%s", str); 1275 else if (!equal) PetscCall(PetscInfo(pc, "%s\n", str)); 1276 PetscFunctionReturn(PETSC_SUCCESS); 1277 } 1278 1279 static PetscErrorCode PCHPDDMDestroySubMatrices_Private(PetscBool flg, PetscBool algebraic, Mat *sub) 1280 { 1281 IS is; 1282 1283 PetscFunctionBegin; 1284 if (!flg) { 1285 if (algebraic) { 1286 PetscCall(PetscObjectQuery((PetscObject)sub[0], "_PCHPDDM_Embed", (PetscObject *)&is)); 1287 PetscCall(ISDestroy(&is)); 1288 PetscCall(PetscObjectCompose((PetscObject)sub[0], "_PCHPDDM_Embed", nullptr)); 1289 PetscCall(PetscObjectCompose((PetscObject)sub[0], "_PCHPDDM_Compact", nullptr)); 1290 } 1291 PetscCall(MatDestroySubMatrices(algebraic ? 2 : 1, &sub)); 1292 } 1293 PetscFunctionReturn(PETSC_SUCCESS); 1294 } 1295 1296 static PetscErrorCode PCHPDDMAlgebraicAuxiliaryMat_Private(Mat P, IS *is, Mat *sub[], PetscBool block) 1297 { 1298 IS icol[3], irow[2]; 1299 Mat *M, Q; 1300 PetscReal *ptr; 1301 PetscInt *idx, p = 0, bs = P->cmap->bs; 1302 PetscBool flg; 1303 1304 PetscFunctionBegin; 1305 PetscCall(ISCreateStride(PETSC_COMM_SELF, P->cmap->N, 0, 1, icol + 2)); 1306 PetscCall(ISSetBlockSize(icol[2], bs)); 1307 PetscCall(ISSetIdentity(icol[2])); 1308 PetscCall(PetscObjectTypeCompare((PetscObject)P, MATMPISBAIJ, &flg)); 1309 if (flg) { 1310 /* MatCreateSubMatrices() does not handle MATMPISBAIJ properly when iscol != isrow, so convert first to MATMPIBAIJ */ 1311 PetscCall(MatConvert(P, MATMPIBAIJ, MAT_INITIAL_MATRIX, &Q)); 1312 std::swap(P, Q); 1313 } 1314 PetscCall(MatCreateSubMatrices(P, 1, is, icol + 2, MAT_INITIAL_MATRIX, &M)); 1315 if (flg) { 1316 std::swap(P, Q); 1317 PetscCall(MatDestroy(&Q)); 1318 } 1319 PetscCall(ISDestroy(icol + 2)); 1320 PetscCall(ISCreateStride(PETSC_COMM_SELF, M[0]->rmap->N, 0, 1, irow)); 1321 PetscCall(ISSetBlockSize(irow[0], bs)); 1322 PetscCall(ISSetIdentity(irow[0])); 1323 if (!block) { 1324 PetscCall(PetscMalloc2(P->cmap->N, &ptr, P->cmap->N / bs, &idx)); 1325 PetscCall(MatGetColumnNorms(M[0], NORM_INFINITY, ptr)); 1326 /* check for nonzero columns so that M[0] may be expressed in compact form */ 1327 for (PetscInt n = 0; n < P->cmap->N; n += bs) { 1328 if (std::find_if(ptr + n, ptr + n + bs, [](PetscReal v) { return v > PETSC_MACHINE_EPSILON; }) != ptr + n + bs) idx[p++] = n / bs; 1329 } 1330 PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, p, idx, PETSC_USE_POINTER, icol + 1)); 1331 PetscCall(ISSetInfo(icol[1], IS_SORTED, IS_GLOBAL, PETSC_TRUE, PETSC_TRUE)); 1332 PetscCall(ISEmbed(*is, icol[1], PETSC_FALSE, icol + 2)); 1333 irow[1] = irow[0]; 1334 /* 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 */ 1335 icol[0] = is[0]; 1336 PetscCall(MatCreateSubMatrices(M[0], 2, irow, icol, MAT_INITIAL_MATRIX, sub)); 1337 PetscCall(ISDestroy(icol + 1)); 1338 PetscCall(PetscFree2(ptr, idx)); 1339 /* IS used to go back and forth between the augmented and the original local linear system, see eq. (3.4) of [2022b] */ 1340 PetscCall(PetscObjectCompose((PetscObject)(*sub)[0], "_PCHPDDM_Embed", (PetscObject)icol[2])); 1341 /* Mat used in eq. (3.1) of [2022b] */ 1342 PetscCall(PetscObjectCompose((PetscObject)(*sub)[0], "_PCHPDDM_Compact", (PetscObject)(*sub)[1])); 1343 } else { 1344 Mat aux; 1345 1346 PetscCall(MatSetOption(M[0], MAT_SUBMAT_SINGLEIS, PETSC_TRUE)); 1347 /* diagonal block of the overlapping rows */ 1348 PetscCall(MatCreateSubMatrices(M[0], 1, irow, is, MAT_INITIAL_MATRIX, sub)); 1349 PetscCall(MatDuplicate((*sub)[0], MAT_COPY_VALUES, &aux)); 1350 PetscCall(MatSetOption(aux, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE)); 1351 if (bs == 1) { /* scalar case */ 1352 Vec sum[2]; 1353 1354 PetscCall(MatCreateVecs(aux, sum, sum + 1)); 1355 PetscCall(MatGetRowSum(M[0], sum[0])); 1356 PetscCall(MatGetRowSum(aux, sum[1])); 1357 /* off-diagonal block row sum (full rows - diagonal block rows) */ 1358 PetscCall(VecAXPY(sum[0], -1.0, sum[1])); 1359 /* subdomain matrix plus off-diagonal block row sum */ 1360 PetscCall(MatDiagonalSet(aux, sum[0], ADD_VALUES)); 1361 PetscCall(VecDestroy(sum)); 1362 PetscCall(VecDestroy(sum + 1)); 1363 } else { /* vectorial case */ 1364 /* TODO: missing MatGetValuesBlocked(), so the code below is */ 1365 /* an extension of the scalar case for when bs > 1, but it could */ 1366 /* be more efficient by avoiding all these MatMatMult() */ 1367 Mat sum[2], ones; 1368 PetscScalar *ptr; 1369 1370 PetscCall(PetscCalloc1(M[0]->cmap->n * bs, &ptr)); 1371 PetscCall(MatCreateDense(PETSC_COMM_SELF, M[0]->cmap->n, bs, M[0]->cmap->n, bs, ptr, &ones)); 1372 for (PetscInt n = 0; n < M[0]->cmap->n; n += bs) { 1373 for (p = 0; p < bs; ++p) ptr[n + p * (M[0]->cmap->n + 1)] = 1.0; 1374 } 1375 PetscCall(MatMatMult(M[0], ones, MAT_INITIAL_MATRIX, PETSC_CURRENT, sum)); 1376 PetscCall(MatDestroy(&ones)); 1377 PetscCall(MatCreateDense(PETSC_COMM_SELF, aux->cmap->n, bs, aux->cmap->n, bs, ptr, &ones)); 1378 PetscCall(MatDenseSetLDA(ones, M[0]->cmap->n)); 1379 PetscCall(MatMatMult(aux, ones, MAT_INITIAL_MATRIX, PETSC_CURRENT, sum + 1)); 1380 PetscCall(MatDestroy(&ones)); 1381 PetscCall(PetscFree(ptr)); 1382 /* off-diagonal block row sum (full rows - diagonal block rows) */ 1383 PetscCall(MatAXPY(sum[0], -1.0, sum[1], SAME_NONZERO_PATTERN)); 1384 PetscCall(MatDestroy(sum + 1)); 1385 /* re-order values to be consistent with MatSetValuesBlocked() */ 1386 /* equivalent to MatTranspose() which does not truly handle */ 1387 /* MAT_INPLACE_MATRIX in the rectangular case, as it calls PetscMalloc() */ 1388 PetscCall(MatDenseGetArrayWrite(sum[0], &ptr)); 1389 HPDDM::Wrapper<PetscScalar>::imatcopy<'T'>(bs, sum[0]->rmap->n, ptr, sum[0]->rmap->n, bs); 1390 /* subdomain matrix plus off-diagonal block row sum */ 1391 for (PetscInt n = 0; n < aux->cmap->n / bs; ++n) PetscCall(MatSetValuesBlocked(aux, 1, &n, 1, &n, ptr + n * bs * bs, ADD_VALUES)); 1392 PetscCall(MatAssemblyBegin(aux, MAT_FINAL_ASSEMBLY)); 1393 PetscCall(MatAssemblyEnd(aux, MAT_FINAL_ASSEMBLY)); 1394 PetscCall(MatDenseRestoreArrayWrite(sum[0], &ptr)); 1395 PetscCall(MatDestroy(sum)); 1396 } 1397 PetscCall(MatSetOption(aux, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE)); 1398 /* left-hand side of GenEO, with the same sparsity pattern as PCASM subdomain solvers */ 1399 PetscCall(PetscObjectCompose((PetscObject)(*sub)[0], "_PCHPDDM_Neumann_Mat", (PetscObject)aux)); 1400 } 1401 PetscCall(ISDestroy(irow)); 1402 PetscCall(MatDestroySubMatrices(1, &M)); 1403 PetscFunctionReturn(PETSC_SUCCESS); 1404 } 1405 1406 static PetscErrorCode PCApply_Nest(PC pc, Vec x, Vec y) 1407 { 1408 Mat A; 1409 MatSolverType type; 1410 IS is[2]; 1411 PetscBool flg; 1412 std::pair<PC, Vec[2]> *p; 1413 1414 PetscFunctionBegin; 1415 PetscCall(PCShellGetContext(pc, &p)); 1416 if (p->second[0]) { /* in case of a centralized Schur complement, some processes may have no local operator */ 1417 PetscCall(PCGetOperators(p->first, &A, nullptr)); 1418 PetscCall(MatNestGetISs(A, is, nullptr)); 1419 PetscCall(PetscObjectTypeCompareAny((PetscObject)p->first, &flg, PCLU, PCCHOLESKY, "")); 1420 if (flg) { /* partial solve currently only makes sense with exact factorizations */ 1421 PetscCall(PCFactorGetMatSolverType(p->first, &type)); 1422 PetscCall(PCFactorGetMatrix(p->first, &A)); 1423 if (A->schur) { 1424 PetscCall(PetscStrcmp(type, MATSOLVERMUMPS, &flg)); 1425 #if PetscDefined(HAVE_MUMPS) 1426 if (flg) PetscCall(MatMumpsSetIcntl(A, 26, 1)); /* reduction/condensation phase followed by Schur complement solve */ 1427 #endif 1428 } else flg = PETSC_FALSE; 1429 } 1430 PetscCall(VecISCopy(p->second[0], is[1], SCATTER_FORWARD, x)); /* assign the RHS associated to the Schur complement */ 1431 PetscCall(PCApply(p->first, p->second[0], p->second[1])); 1432 PetscCall(VecISCopy(p->second[1], is[1], SCATTER_REVERSE, y)); /* retrieve the partial solution associated to the Schur complement */ 1433 #if PetscDefined(HAVE_MUMPS) 1434 if (flg) PetscCall(MatMumpsSetIcntl(A, 26, -1)); /* default ICNTL(26) value in PETSc */ 1435 #endif 1436 } 1437 PetscFunctionReturn(PETSC_SUCCESS); 1438 } 1439 1440 static PetscErrorCode PCView_Nest(PC pc, PetscViewer viewer) 1441 { 1442 std::pair<PC, Vec[2]> *p; 1443 1444 PetscFunctionBegin; 1445 PetscCall(PCShellGetContext(pc, &p)); 1446 PetscCall(PCView(p->first, viewer)); 1447 PetscFunctionReturn(PETSC_SUCCESS); 1448 } 1449 1450 static PetscErrorCode PCDestroy_Nest(PC pc) 1451 { 1452 std::pair<PC, Vec[2]> *p; 1453 1454 PetscFunctionBegin; 1455 PetscCall(PCShellGetContext(pc, &p)); 1456 PetscCall(VecDestroy(p->second)); 1457 PetscCall(VecDestroy(p->second + 1)); 1458 PetscCall(PCDestroy(&p->first)); 1459 PetscCall(PetscFree(p)); 1460 PetscFunctionReturn(PETSC_SUCCESS); 1461 } 1462 1463 template <bool T = false> 1464 static PetscErrorCode MatMult_Schur(Mat A, Vec x, Vec y) 1465 { 1466 std::tuple<Mat, PetscSF, Vec[2]> *ctx; 1467 1468 PetscFunctionBegin; 1469 PetscCall(MatShellGetContext(A, &ctx)); 1470 PetscCall(VecScatterBegin(std::get<1>(*ctx), x, std::get<2>(*ctx)[0], INSERT_VALUES, SCATTER_FORWARD)); /* local Vec with overlap */ 1471 PetscCall(VecScatterEnd(std::get<1>(*ctx), x, std::get<2>(*ctx)[0], INSERT_VALUES, SCATTER_FORWARD)); 1472 if (!T) PetscCall(MatMult(std::get<0>(*ctx), std::get<2>(*ctx)[0], std::get<2>(*ctx)[1])); /* local Schur complement */ 1473 else PetscCall(MatMultTranspose(std::get<0>(*ctx), std::get<2>(*ctx)[0], std::get<2>(*ctx)[1])); 1474 PetscCall(VecSet(y, 0.0)); 1475 PetscCall(VecScatterBegin(std::get<1>(*ctx), std::get<2>(*ctx)[1], y, ADD_VALUES, SCATTER_REVERSE)); /* global Vec with summed up contributions on the overlap */ 1476 PetscCall(VecScatterEnd(std::get<1>(*ctx), std::get<2>(*ctx)[1], y, ADD_VALUES, SCATTER_REVERSE)); 1477 PetscFunctionReturn(PETSC_SUCCESS); 1478 } 1479 1480 static PetscErrorCode MatDestroy_Schur(Mat A) 1481 { 1482 std::tuple<Mat, PetscSF, Vec[2]> *ctx; 1483 1484 PetscFunctionBegin; 1485 PetscCall(MatShellGetContext(A, &ctx)); 1486 PetscCall(VecDestroy(std::get<2>(*ctx))); 1487 PetscCall(VecDestroy(std::get<2>(*ctx) + 1)); 1488 PetscCall(PetscFree(ctx)); 1489 PetscFunctionReturn(PETSC_SUCCESS); 1490 } 1491 1492 static PetscErrorCode MatMult_SchurCorrection(Mat A, Vec x, Vec y) 1493 { 1494 PC pc; 1495 std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *ctx; 1496 1497 PetscFunctionBegin; 1498 PetscCall(MatShellGetContext(A, &ctx)); 1499 pc = ((PC_HPDDM *)std::get<0>(*ctx)[0]->data)->levels[0]->ksp->pc; 1500 if (std::get<2>(*ctx) == PC_LEFT || std::get<2>(*ctx) == PC_SIDE_DEFAULT) { /* Q_0 is the coarse correction associated to the A00 block from PCFIELDSPLIT */ 1501 PetscCall(MatMult(std::get<1>(*ctx)[0], x, std::get<3>(*ctx)[1])); /* A_01 x */ 1502 PetscCall(PCHPDDMDeflate_Private(pc, std::get<3>(*ctx)[1], std::get<3>(*ctx)[1])); /* Q_0 A_01 x */ 1503 PetscCall(MatMult(std::get<1>(*ctx)[1], std::get<3>(*ctx)[1], std::get<3>(*ctx)[0])); /* A_10 Q_0 A_01 x */ 1504 PetscCall(PCApply(std::get<0>(*ctx)[1], std::get<3>(*ctx)[0], y)); /* y = M_S^-1 A_10 Q_0 A_01 x */ 1505 } else { 1506 PetscCall(PCApply(std::get<0>(*ctx)[1], x, std::get<3>(*ctx)[0])); /* M_S^-1 x */ 1507 PetscCall(MatMult(std::get<1>(*ctx)[0], std::get<3>(*ctx)[0], std::get<3>(*ctx)[1])); /* A_01 M_S^-1 x */ 1508 PetscCall(PCHPDDMDeflate_Private(pc, std::get<3>(*ctx)[1], std::get<3>(*ctx)[1])); /* Q_0 A_01 M_S^-1 x */ 1509 PetscCall(MatMult(std::get<1>(*ctx)[1], std::get<3>(*ctx)[1], y)); /* y = A_10 Q_0 A_01 M_S^-1 x */ 1510 } 1511 PetscCall(VecAXPY(y, -1.0, x)); /* y -= x, preconditioned eq. (24) of https://hal.science/hal-02343808v6/document (with a sign flip) */ 1512 PetscFunctionReturn(PETSC_SUCCESS); 1513 } 1514 1515 static PetscErrorCode MatView_SchurCorrection(Mat A, PetscViewer viewer) 1516 { 1517 PetscBool ascii; 1518 std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *ctx; 1519 1520 PetscFunctionBegin; 1521 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &ascii)); 1522 if (ascii) { 1523 PetscCall(MatShellGetContext(A, &ctx)); 1524 PetscCall(PetscViewerASCIIPrintf(viewer, "action of %s\n", std::get<2>(*ctx) == PC_LEFT || std::get<2>(*ctx) == PC_SIDE_DEFAULT ? "(I - M_S^-1 A_10 Q_0 A_01)" : "(I - A_10 Q_0 A_01 M_S^-1)")); 1525 PetscCall(PCView(std::get<0>(*ctx)[1], viewer)); /* no need to PCView(Q_0) since it will be done by PCFIELDSPLIT */ 1526 } 1527 PetscFunctionReturn(PETSC_SUCCESS); 1528 } 1529 1530 static PetscErrorCode MatDestroy_SchurCorrection(Mat A) 1531 { 1532 std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *ctx; 1533 1534 PetscFunctionBegin; 1535 PetscCall(MatShellGetContext(A, &ctx)); 1536 PetscCall(VecDestroy(std::get<3>(*ctx))); 1537 PetscCall(VecDestroy(std::get<3>(*ctx) + 1)); 1538 PetscCall(VecDestroy(std::get<3>(*ctx) + 2)); 1539 PetscCall(PCDestroy(std::get<0>(*ctx) + 1)); 1540 PetscCall(PetscFree(ctx)); 1541 PetscFunctionReturn(PETSC_SUCCESS); 1542 } 1543 1544 static PetscErrorCode PCPostSolve_SchurPreLeastSquares(PC, KSP, Vec, Vec x) 1545 { 1546 PetscFunctionBegin; 1547 PetscCall(VecScale(x, -1.0)); 1548 PetscFunctionReturn(PETSC_SUCCESS); 1549 } 1550 1551 static PetscErrorCode KSPPreSolve_SchurCorrection(KSP, Vec b, Vec, void *context) 1552 { 1553 std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *ctx = reinterpret_cast<std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *>(context); 1554 1555 PetscFunctionBegin; 1556 if (std::get<2>(*ctx) == PC_LEFT || std::get<2>(*ctx) == PC_SIDE_DEFAULT) { 1557 PetscCall(PCApply(std::get<0>(*ctx)[1], b, std::get<3>(*ctx)[2])); 1558 std::swap(*b, *std::get<3>(*ctx)[2]); /* replace b by M^-1 b, but need to keep a copy of the original RHS, so swap it with the work Vec */ 1559 } 1560 PetscFunctionReturn(PETSC_SUCCESS); 1561 } 1562 1563 static PetscErrorCode KSPPostSolve_SchurCorrection(KSP, Vec b, Vec x, void *context) 1564 { 1565 std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *ctx = reinterpret_cast<std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *>(context); 1566 1567 PetscFunctionBegin; 1568 if (std::get<2>(*ctx) == PC_LEFT || std::get<2>(*ctx) == PC_SIDE_DEFAULT) std::swap(*b, *std::get<3>(*ctx)[2]); /* put back the original RHS where it belongs */ 1569 else { 1570 PetscCall(PCApply(std::get<0>(*ctx)[1], x, std::get<3>(*ctx)[2])); 1571 PetscCall(VecCopy(std::get<3>(*ctx)[2], x)); /* replace x by M^-1 x */ 1572 } 1573 PetscFunctionReturn(PETSC_SUCCESS); 1574 } 1575 1576 static PetscErrorCode MatMult_Harmonic(Mat, Vec, Vec); 1577 static PetscErrorCode MatMultTranspose_Harmonic(Mat, Vec, Vec); 1578 static PetscErrorCode MatProduct_AB_Harmonic(Mat, Mat, Mat, void *); 1579 static PetscErrorCode MatDestroy_Harmonic(Mat); 1580 1581 static PetscErrorCode PCSetUp_HPDDM(PC pc) 1582 { 1583 PC_HPDDM *data = (PC_HPDDM *)pc->data; 1584 PC inner; 1585 KSP *ksp; 1586 Mat *sub, A, P, N, C = nullptr, uaux = nullptr, weighted, subA[2], S; 1587 Vec xin, v; 1588 std::vector<Vec> initial; 1589 IS is[1], loc, uis = data->is, unsorted = nullptr; 1590 ISLocalToGlobalMapping l2g; 1591 char prefix[256]; 1592 const char *pcpre; 1593 const PetscScalar *const *ev; 1594 PetscInt n, requested = data->N, reused = 0, overlap = -1; 1595 MatStructure structure = UNKNOWN_NONZERO_PATTERN; 1596 PetscBool subdomains = PETSC_FALSE, flg = PETSC_FALSE, ismatis, swap = PETSC_FALSE, algebraic = PETSC_FALSE, block = PETSC_FALSE; 1597 DM dm; 1598 std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *ctx = nullptr; 1599 #if PetscDefined(USE_DEBUG) 1600 IS dis = nullptr; 1601 Mat daux = nullptr; 1602 #endif 1603 1604 PetscFunctionBegin; 1605 PetscCheck(data->levels && data->levels[0], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Not a single level allocated"); 1606 PetscCall(PCGetOptionsPrefix(pc, &pcpre)); 1607 PetscCall(PCGetOperators(pc, &A, &P)); 1608 if (!data->levels[0]->ksp) { 1609 PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &data->levels[0]->ksp)); 1610 PetscCall(KSPSetNestLevel(data->levels[0]->ksp, pc->kspnestlevel)); 1611 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_%s_", pcpre ? pcpre : "", data->N > 1 ? "levels_1" : "coarse")); 1612 PetscCall(KSPSetOptionsPrefix(data->levels[0]->ksp, prefix)); 1613 PetscCall(KSPSetType(data->levels[0]->ksp, KSPPREONLY)); 1614 } else if (data->levels[0]->ksp->pc && data->levels[0]->ksp->pc->setupcalled == 1 && data->levels[0]->ksp->pc->reusepreconditioner) { 1615 /* if the fine-level PCSHELL exists, its setup has succeeded, and one wants to reuse it, */ 1616 /* then just propagate the appropriate flag to the coarser levels */ 1617 for (n = 0; n < PETSC_PCHPDDM_MAXLEVELS && data->levels[n]; ++n) { 1618 /* the following KSP and PC may be NULL for some processes, hence the check */ 1619 if (data->levels[n]->ksp) PetscCall(KSPSetReusePreconditioner(data->levels[n]->ksp, PETSC_TRUE)); 1620 if (data->levels[n]->pc) PetscCall(PCSetReusePreconditioner(data->levels[n]->pc, PETSC_TRUE)); 1621 } 1622 /* early bail out because there is nothing to do */ 1623 PetscFunctionReturn(PETSC_SUCCESS); 1624 } else { 1625 /* reset coarser levels */ 1626 for (n = 1; n < PETSC_PCHPDDM_MAXLEVELS && data->levels[n]; ++n) { 1627 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) { 1628 reused = data->N - n; 1629 break; 1630 } 1631 PetscCall(KSPDestroy(&data->levels[n]->ksp)); 1632 PetscCall(PCDestroy(&data->levels[n]->pc)); 1633 } 1634 /* check if some coarser levels are being reused */ 1635 PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &reused, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)pc))); 1636 const int *addr = data->levels[0]->P ? data->levels[0]->P->getAddrLocal() : &HPDDM::i__0; 1637 1638 if (addr != &HPDDM::i__0 && reused != data->N - 1) { 1639 /* reuse previously computed eigenvectors */ 1640 ev = data->levels[0]->P->getVectors(); 1641 if (ev) { 1642 initial.reserve(*addr); 1643 PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, data->levels[0]->P->getDof(), ev[0], &xin)); 1644 for (n = 0; n < *addr; ++n) { 1645 PetscCall(VecDuplicate(xin, &v)); 1646 PetscCall(VecPlaceArray(xin, ev[n])); 1647 PetscCall(VecCopy(xin, v)); 1648 initial.emplace_back(v); 1649 PetscCall(VecResetArray(xin)); 1650 } 1651 PetscCall(VecDestroy(&xin)); 1652 } 1653 } 1654 } 1655 data->N -= reused; 1656 PetscCall(KSPSetOperators(data->levels[0]->ksp, A, P)); 1657 1658 PetscCall(PetscObjectTypeCompare((PetscObject)P, MATIS, &ismatis)); 1659 if (!data->is && !ismatis) { 1660 PetscErrorCode (*create)(DM, IS *, Mat *, PetscErrorCode (**)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void *), void **) = nullptr; 1661 PetscErrorCode (*usetup)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void *) = nullptr; 1662 void *uctx = nullptr; 1663 1664 /* first see if we can get the data from the DM */ 1665 PetscCall(MatGetDM(P, &dm)); 1666 if (!dm) PetscCall(MatGetDM(A, &dm)); 1667 if (!dm) PetscCall(PCGetDM(pc, &dm)); 1668 if (dm) { /* this is the hook for DMPLEX for which the auxiliary Mat is the local Neumann matrix */ 1669 PetscCall(PetscObjectQueryFunction((PetscObject)dm, "DMCreateNeumannOverlap_C", &create)); 1670 if (create) { 1671 PetscCall((*create)(dm, &uis, &uaux, &usetup, &uctx)); 1672 if (data->Neumann == PETSC_BOOL3_UNKNOWN) data->Neumann = PETSC_BOOL3_TRUE; /* set the value only if it was not already provided by the user */ 1673 } 1674 } 1675 if (!create) { 1676 if (!uis) { 1677 PetscCall(PetscObjectQuery((PetscObject)pc, "_PCHPDDM_Neumann_IS", (PetscObject *)&uis)); 1678 PetscCall(PetscObjectReference((PetscObject)uis)); 1679 } 1680 if (!uaux) { 1681 PetscCall(PetscObjectQuery((PetscObject)pc, "_PCHPDDM_Neumann_Mat", (PetscObject *)&uaux)); 1682 PetscCall(PetscObjectReference((PetscObject)uaux)); 1683 } 1684 /* look inside the Pmat instead of the PC, needed for MatSchurComplementComputeExplicitOperator() */ 1685 if (!uis) { 1686 PetscCall(PetscObjectQuery((PetscObject)P, "_PCHPDDM_Neumann_IS", (PetscObject *)&uis)); 1687 PetscCall(PetscObjectReference((PetscObject)uis)); 1688 } 1689 if (!uaux) { 1690 PetscCall(PetscObjectQuery((PetscObject)P, "_PCHPDDM_Neumann_Mat", (PetscObject *)&uaux)); 1691 PetscCall(PetscObjectReference((PetscObject)uaux)); 1692 } 1693 } 1694 PetscCall(PCHPDDMSetAuxiliaryMat(pc, uis, uaux, usetup, uctx)); 1695 PetscCall(MatDestroy(&uaux)); 1696 PetscCall(ISDestroy(&uis)); 1697 } 1698 1699 if (!ismatis) { 1700 PetscCall(PCHPDDMSetUpNeumannOverlap_Private(pc)); 1701 PetscCall(PetscOptionsGetBool(nullptr, pcpre, "-pc_hpddm_block_splitting", &block, nullptr)); 1702 PetscCall(PetscOptionsGetInt(nullptr, pcpre, "-pc_hpddm_harmonic_overlap", &overlap, nullptr)); 1703 PetscCall(PetscObjectTypeCompare((PetscObject)P, MATSCHURCOMPLEMENT, &flg)); 1704 if (data->is || flg) { 1705 if (block || overlap != -1) { 1706 PetscCall(ISDestroy(&data->is)); 1707 PetscCall(MatDestroy(&data->aux)); 1708 } else if (flg) { 1709 PCHPDDMSchurPreType type = PC_HPDDM_SCHUR_PRE_GENEO; 1710 1711 PetscCall(PetscOptionsGetEnum(nullptr, pcpre, "-pc_hpddm_schur_precondition", PCHPDDMSchurPreTypes, (PetscEnum *)&type, &flg)); 1712 if (type == PC_HPDDM_SCHUR_PRE_LEAST_SQUARES) { 1713 PetscCall(ISDestroy(&data->is)); /* destroy any previously user-set objects since they will be set automatically */ 1714 PetscCall(MatDestroy(&data->aux)); 1715 } else if (type == PC_HPDDM_SCHUR_PRE_GENEO) { 1716 PetscContainer container = nullptr; 1717 1718 PetscCall(PetscObjectQuery((PetscObject)pc, "_PCHPDDM_Schur", (PetscObject *)&container)); 1719 if (!container) { /* first call to PCSetUp() on the PC associated to the Schur complement */ 1720 PC_HPDDM *data_00; 1721 KSP ksp, inner_ksp; 1722 PC pc_00; 1723 Mat A11 = nullptr; 1724 Vec d = nullptr; 1725 const PetscInt *ranges; 1726 PetscMPIInt size; 1727 char *prefix; 1728 1729 PetscCall(MatSchurComplementGetKSP(P, &ksp)); 1730 PetscCall(KSPGetPC(ksp, &pc_00)); 1731 PetscCall(PetscObjectTypeCompare((PetscObject)pc_00, PCHPDDM, &flg)); 1732 PetscCheck(flg, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-%spc_hpddm_schur_precondition %s and -%spc_type %s (!= %s)", pcpre ? pcpre : "", PCHPDDMSchurPreTypes[type], ((PetscObject)pc_00)->prefix ? ((PetscObject)pc_00)->prefix : "", 1733 ((PetscObject)pc_00)->type_name, PCHPDDM); 1734 data_00 = (PC_HPDDM *)pc_00->data; 1735 PetscCheck(data_00->N == 2, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-%spc_hpddm_schur_precondition %s and %" PetscInt_FMT " level%s instead of 2 for the A00 block -%s", pcpre ? pcpre : "", PCHPDDMSchurPreTypes[type], 1736 data_00->N, data_00->N > 1 ? "s" : "", ((PetscObject)pc_00)->prefix); 1737 PetscCheck(data_00->levels[0]->pc, PetscObjectComm((PetscObject)P), PETSC_ERR_ORDER, "PC of the first block%s not setup yet", ((PetscObject)pc_00)->prefix ? std::string(std::string(" (") + ((PetscObject)pc_00)->prefix + std::string(")")).c_str() : ""); 1738 PetscCall(PetscObjectTypeCompare((PetscObject)data_00->levels[0]->pc, PCASM, &flg)); 1739 PetscCheck(flg, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-%spc_hpddm_schur_precondition %s and -%spc_type %s (!= %s)", pcpre ? pcpre : "", PCHPDDMSchurPreTypes[type], ((PetscObject)data_00->levels[0]->pc)->prefix, 1740 ((PetscObject)data_00->levels[0]->pc)->type_name, PCASM); 1741 PetscCall(PetscNew(&ctx)); /* context to pass data around for the inner-most PC, which will be a proper PCHPDDM (or a dummy variable if the Schur complement is centralized on a single process) */ 1742 PetscCall(MatSchurComplementGetSubMatrices(P, nullptr, nullptr, nullptr, nullptr, &A11)); 1743 PetscCall(MatGetOwnershipRanges(A11, &ranges)); 1744 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A11), &size)); 1745 flg = PetscBool(std::find_if(ranges, ranges + size + 1, [&](PetscInt v) { return v != ranges[0] && v != ranges[size]; }) == ranges + size + 1); /* are all local matrices but one of dimension 0 (centralized Schur complement)? */ 1746 if (!flg) { 1747 if (PetscDefined(USE_DEBUG) || !data->is) { 1748 Mat A01, A10, B = nullptr, C = nullptr, *sub; 1749 1750 PetscCall(MatSchurComplementGetSubMatrices(P, &A, nullptr, &A01, &A10, nullptr)); 1751 PetscCall(PetscObjectTypeCompare((PetscObject)A10, MATTRANSPOSEVIRTUAL, &flg)); 1752 if (flg) { 1753 PetscCall(MatTransposeGetMat(A10, &C)); 1754 PetscCall(MatTranspose(C, MAT_INITIAL_MATRIX, &B)); 1755 } else { 1756 PetscCall(PetscObjectTypeCompare((PetscObject)A10, MATHERMITIANTRANSPOSEVIRTUAL, &flg)); 1757 if (flg) { 1758 PetscCall(MatHermitianTransposeGetMat(A10, &C)); 1759 PetscCall(MatHermitianTranspose(C, MAT_INITIAL_MATRIX, &B)); 1760 } 1761 } 1762 if (flg) 1763 PetscCall(MatShellGetScalingShifts(A10, (PetscScalar *)MAT_SHELL_NOT_ALLOWED, (PetscScalar *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Mat *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED)); 1764 if (!B) { 1765 B = A10; 1766 PetscCall(PetscObjectReference((PetscObject)B)); 1767 } else if (!data->is) { 1768 PetscCall(PetscObjectTypeCompareAny((PetscObject)A01, &flg, MATTRANSPOSEVIRTUAL, MATHERMITIANTRANSPOSEVIRTUAL, "")); 1769 if (!flg) C = A01; 1770 else 1771 PetscCall(MatShellGetScalingShifts(A01, (PetscScalar *)MAT_SHELL_NOT_ALLOWED, (PetscScalar *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Mat *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED)); 1772 } 1773 PetscCall(ISCreateStride(PETSC_COMM_SELF, B->rmap->N, 0, 1, &uis)); 1774 PetscCall(ISSetIdentity(uis)); 1775 if (!data->is) { 1776 if (C) PetscCall(PetscObjectReference((PetscObject)C)); 1777 else PetscCall(MatTranspose(B, MAT_INITIAL_MATRIX, &C)); 1778 PetscCall(ISDuplicate(data_00->is, is)); 1779 PetscCall(MatIncreaseOverlap(A, 1, is, 1)); 1780 PetscCall(MatSetOption(C, MAT_SUBMAT_SINGLEIS, PETSC_TRUE)); 1781 PetscCall(MatCreateSubMatrices(C, 1, is, &uis, MAT_INITIAL_MATRIX, &sub)); 1782 PetscCall(MatDestroy(&C)); 1783 PetscCall(MatTranspose(sub[0], MAT_INITIAL_MATRIX, &C)); 1784 PetscCall(MatDestroySubMatrices(1, &sub)); 1785 PetscCall(MatFindNonzeroRows(C, &data->is)); 1786 PetscCall(MatDestroy(&C)); 1787 PetscCall(ISDestroy(is)); 1788 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)data->is), A11->rmap->n, A11->rmap->rstart, 1, &loc)); 1789 if (PetscDefined(USE_DEBUG)) PetscCall(PCHPDDMCheckInclusion_Private(pc, data->is, loc, PETSC_FALSE)); 1790 PetscCall(ISExpand(data->is, loc, is)); 1791 PetscCall(ISDestroy(&loc)); 1792 PetscCall(ISDestroy(&data->is)); 1793 data->is = is[0]; 1794 is[0] = nullptr; 1795 } 1796 if (PetscDefined(USE_DEBUG)) { 1797 PetscCall(PCHPDDMCheckSymmetry_Private(pc, A01, A10)); 1798 PetscCall(MatCreateSubMatrices(B, 1, &uis, &data_00->is, MAT_INITIAL_MATRIX, &sub)); /* expensive check since all processes fetch all rows (but only some columns) of the constraint matrix */ 1799 PetscCall(ISDestroy(&uis)); 1800 PetscCall(ISDuplicate(data->is, &uis)); 1801 PetscCall(ISSort(uis)); 1802 PetscCall(ISComplement(uis, 0, B->rmap->N, is)); 1803 PetscCall(MatDuplicate(sub[0], MAT_COPY_VALUES, &C)); 1804 PetscCall(MatZeroRowsIS(C, is[0], 0.0, nullptr, nullptr)); 1805 PetscCall(ISDestroy(is)); 1806 PetscCall(MatMultEqual(sub[0], C, 20, &flg)); 1807 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "The image of A_10 (R_i^p)^T from the local primal (e.g., velocity) space to the full dual (e.g., pressure) space is not restricted to the local dual space: A_10 (R_i^p)^T != R_i^d (R_i^d)^T A_10 (R_i^p)^T"); /* cf. eq. (9) of https://hal.science/hal-02343808v6/document */ 1808 PetscCall(MatDestroy(&C)); 1809 PetscCall(MatDestroySubMatrices(1, &sub)); 1810 } 1811 PetscCall(ISDestroy(&uis)); 1812 PetscCall(MatDestroy(&B)); 1813 } 1814 flg = PETSC_FALSE; 1815 if (!data->aux) { 1816 Mat D; 1817 1818 PetscCall(MatCreateVecs(A11, &d, nullptr)); 1819 PetscCall(MatGetDiagonal(A11, d)); 1820 PetscCall(PetscObjectTypeCompareAny((PetscObject)A11, &flg, MATDIAGONAL, MATCONSTANTDIAGONAL, "")); 1821 if (!flg) { 1822 PetscCall(MatCreateDiagonal(d, &D)); 1823 PetscCall(MatMultEqual(A11, D, 20, &flg)); 1824 PetscCall(MatDestroy(&D)); 1825 } 1826 if (flg) PetscCall(PetscInfo(pc, "A11 block is likely diagonal so the PC will build an auxiliary Mat (which was not initially provided by the user)\n")); 1827 } 1828 if (data->Neumann != PETSC_BOOL3_TRUE && !flg && A11) { 1829 PetscReal norm; 1830 1831 PetscCall(MatNorm(A11, NORM_INFINITY, &norm)); 1832 PetscCheck(norm < PETSC_MACHINE_EPSILON * static_cast<PetscReal>(10.0), PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-%spc_hpddm_schur_precondition geneo and -%spc_hpddm_has_neumann != true with a nonzero or non-diagonal A11 block", pcpre ? pcpre : "", pcpre ? pcpre : ""); 1833 PetscCall(PetscInfo(pc, "A11 block is likely zero so the PC will build an auxiliary Mat (which was%s initially provided by the user)\n", data->aux ? "" : " not")); 1834 PetscCall(MatDestroy(&data->aux)); 1835 flg = PETSC_TRUE; 1836 } 1837 if (!data->aux) { /* if A11 is near zero, e.g., Stokes equation, or diagonal, build an auxiliary (Neumann) Mat which is a (possibly slightly shifted) diagonal weighted by the inverse of the multiplicity */ 1838 PetscSF scatter; 1839 const PetscScalar *read; 1840 PetscScalar *write, *diagonal = nullptr; 1841 1842 PetscCall(MatDestroy(&data->aux)); 1843 PetscCall(ISGetLocalSize(data->is, &n)); 1844 PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)P), n, PETSC_DECIDE, &xin)); 1845 PetscCall(VecDuplicate(xin, &v)); 1846 PetscCall(VecScatterCreate(xin, data->is, v, nullptr, &scatter)); 1847 PetscCall(VecSet(v, 1.0)); 1848 PetscCall(VecSet(xin, 1.0)); 1849 PetscCall(VecScatterBegin(scatter, v, xin, ADD_VALUES, SCATTER_REVERSE)); 1850 PetscCall(VecScatterEnd(scatter, v, xin, ADD_VALUES, SCATTER_REVERSE)); /* v has the multiplicity of all unknowns on the overlap */ 1851 PetscCall(PetscSFDestroy(&scatter)); 1852 if (d) { 1853 PetscCall(VecScatterCreate(d, data->is, v, nullptr, &scatter)); 1854 PetscCall(VecScatterBegin(scatter, d, v, INSERT_VALUES, SCATTER_FORWARD)); 1855 PetscCall(VecScatterEnd(scatter, d, v, INSERT_VALUES, SCATTER_FORWARD)); 1856 PetscCall(PetscSFDestroy(&scatter)); 1857 PetscCall(VecDestroy(&d)); 1858 PetscCall(PetscMalloc1(n, &diagonal)); 1859 PetscCall(VecGetArrayRead(v, &read)); 1860 PetscCallCXX(std::copy_n(read, n, diagonal)); 1861 PetscCall(VecRestoreArrayRead(v, &read)); 1862 } 1863 PetscCall(VecDestroy(&v)); 1864 PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &v)); 1865 PetscCall(VecGetArrayRead(xin, &read)); 1866 PetscCall(VecGetArrayWrite(v, &write)); 1867 for (PetscInt i = 0; i < n; ++i) write[i] = (!diagonal || std::abs(diagonal[i]) < PETSC_MACHINE_EPSILON) ? PETSC_SMALL / (static_cast<PetscReal>(1000.0) * read[i]) : diagonal[i] / read[i]; 1868 PetscCall(PetscFree(diagonal)); 1869 PetscCall(VecRestoreArrayRead(xin, &read)); 1870 PetscCall(VecRestoreArrayWrite(v, &write)); 1871 PetscCall(VecDestroy(&xin)); 1872 PetscCall(MatCreateDiagonal(v, &data->aux)); 1873 PetscCall(VecDestroy(&v)); 1874 } 1875 uis = data->is; 1876 uaux = data->aux; 1877 PetscCall(PetscObjectReference((PetscObject)uis)); 1878 PetscCall(PetscObjectReference((PetscObject)uaux)); 1879 PetscCall(PetscStrallocpy(pcpre, &prefix)); 1880 PetscCall(PCSetOptionsPrefix(pc, nullptr)); 1881 PetscCall(PCSetType(pc, PCKSP)); /* replace the PC associated to the Schur complement by PCKSP */ 1882 PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &inner_ksp)); /* new KSP that will be attached to the previously set PC */ 1883 PetscCall(PetscObjectGetTabLevel((PetscObject)pc, &n)); 1884 PetscCall(PetscObjectSetTabLevel((PetscObject)inner_ksp, n + 2)); 1885 PetscCall(KSPSetOperators(inner_ksp, pc->mat, pc->pmat)); 1886 PetscCall(KSPSetOptionsPrefix(inner_ksp, std::string(std::string(prefix) + "pc_hpddm_").c_str())); 1887 PetscCall(KSPSetSkipPCSetFromOptions(inner_ksp, PETSC_TRUE)); 1888 PetscCall(KSPSetFromOptions(inner_ksp)); 1889 PetscCall(KSPGetPC(inner_ksp, &inner)); 1890 PetscCall(PCSetOptionsPrefix(inner, nullptr)); 1891 PetscCall(PCSetType(inner, PCNONE)); /* no preconditioner since the action of M^-1 A or A M^-1 will be computed by the Amat */ 1892 PetscCall(PCKSPSetKSP(pc, inner_ksp)); 1893 std::get<0>(*ctx)[0] = pc_00; /* for coarse correction on the primal (e.g., velocity) space */ 1894 PetscCall(PCCreate(PetscObjectComm((PetscObject)pc), &std::get<0>(*ctx)[1])); 1895 PetscCall(PCSetOptionsPrefix(pc, prefix)); /* both PC share the same prefix so that the outer PC can be reset with PCSetFromOptions() */ 1896 PetscCall(PCSetOptionsPrefix(std::get<0>(*ctx)[1], prefix)); 1897 PetscCall(PetscFree(prefix)); 1898 PetscCall(PCSetOperators(std::get<0>(*ctx)[1], pc->mat, pc->pmat)); 1899 PetscCall(PCSetType(std::get<0>(*ctx)[1], PCHPDDM)); 1900 PetscCall(PCHPDDMSetAuxiliaryMat(std::get<0>(*ctx)[1], uis, uaux, nullptr, nullptr)); /* transfer ownership of the auxiliary inputs from the inner (PCKSP) to the inner-most (PCHPDDM) PC */ 1901 if (flg) static_cast<PC_HPDDM *>(std::get<0>(*ctx)[1]->data)->Neumann = PETSC_BOOL3_TRUE; 1902 PetscCall(PCSetFromOptions(std::get<0>(*ctx)[1])); 1903 PetscCall(PetscObjectDereference((PetscObject)uis)); 1904 PetscCall(PetscObjectDereference((PetscObject)uaux)); 1905 PetscCall(MatCreateShell(PetscObjectComm((PetscObject)pc), inner->mat->rmap->n, inner->mat->cmap->n, inner->mat->rmap->N, inner->mat->cmap->N, ctx, &S)); /* MatShell computing the action of M^-1 A or A M^-1 */ 1906 PetscCall(MatShellSetOperation(S, MATOP_MULT, (void (*)(void))MatMult_SchurCorrection)); 1907 PetscCall(MatShellSetOperation(S, MATOP_VIEW, (void (*)(void))MatView_SchurCorrection)); 1908 PetscCall(MatShellSetOperation(S, MATOP_DESTROY, (void (*)(void))MatDestroy_SchurCorrection)); 1909 PetscCall(KSPGetPCSide(inner_ksp, &(std::get<2>(*ctx)))); 1910 if (std::get<2>(*ctx) == PC_LEFT || std::get<2>(*ctx) == PC_SIDE_DEFAULT) { 1911 PetscCall(KSPSetPreSolve(inner_ksp, KSPPreSolve_SchurCorrection, ctx)); 1912 } else { /* no support for PC_SYMMETRIC */ 1913 PetscCheck(std::get<2>(*ctx) == PC_RIGHT, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "PCSide %s (!= %s or %s or %s)", PCSides[std::get<2>(*ctx)], PCSides[PC_SIDE_DEFAULT], PCSides[PC_LEFT], PCSides[PC_RIGHT]); 1914 } 1915 PetscCall(KSPSetPostSolve(inner_ksp, KSPPostSolve_SchurCorrection, ctx)); 1916 PetscCall(PetscObjectContainerCompose((PetscObject)std::get<0>(*ctx)[1], "_PCHPDDM_Schur", ctx, nullptr)); 1917 PetscCall(PCSetUp(std::get<0>(*ctx)[1])); 1918 PetscCall(KSPSetOperators(inner_ksp, S, S)); 1919 PetscCall(MatCreateVecs(std::get<1>(*ctx)[0], std::get<3>(*ctx), std::get<3>(*ctx) + 1)); 1920 PetscCall(VecDuplicate(std::get<3>(*ctx)[0], std::get<3>(*ctx) + 2)); 1921 PetscCall(PetscObjectDereference((PetscObject)inner_ksp)); 1922 PetscCall(PetscObjectDereference((PetscObject)S)); 1923 } else { 1924 std::get<0>(*ctx)[0] = pc_00; 1925 PetscCall(PetscObjectContainerCompose((PetscObject)pc, "_PCHPDDM_Schur", ctx, nullptr)); 1926 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)data_00->is), A11->rmap->n, A11->rmap->rstart, 1, &data->is)); /* dummy variables in the case of a centralized Schur complement */ 1927 PetscCall(MatGetDiagonalBlock(A11, &data->aux)); 1928 PetscCall(PetscObjectReference((PetscObject)data->aux)); 1929 PetscCall(PCSetUp(pc)); 1930 } 1931 for (std::vector<Vec>::iterator it = initial.begin(); it != initial.end(); ++it) PetscCall(VecDestroy(&*it)); 1932 PetscFunctionReturn(PETSC_SUCCESS); 1933 } else { /* second call to PCSetUp() on the PC associated to the Schur complement, retrieve previously set context */ 1934 PetscCall(PetscContainerGetPointer(container, (void **)&ctx)); 1935 } 1936 } 1937 } 1938 } 1939 if (!data->is && data->N > 1) { 1940 char type[256] = {}; /* same size as in src/ksp/pc/interface/pcset.c */ 1941 1942 PetscCall(PetscObjectTypeCompareAny((PetscObject)P, &flg, MATNORMAL, MATNORMALHERMITIAN, "")); 1943 if (flg || (A->rmap->N != A->cmap->N && P->rmap->N == P->cmap->N && P->rmap->N == A->cmap->N)) { 1944 Mat B; 1945 1946 PetscCall(PCHPDDMSetAuxiliaryMatNormal_Private(pc, A, P, &B, pcpre)); 1947 if (data->correction == PC_HPDDM_COARSE_CORRECTION_DEFLATED) data->correction = PC_HPDDM_COARSE_CORRECTION_BALANCED; 1948 PetscCall(MatDestroy(&B)); 1949 } else { 1950 PetscCall(PetscObjectTypeCompare((PetscObject)P, MATSCHURCOMPLEMENT, &flg)); 1951 if (flg) { 1952 Mat A00, P00, A01, A10, A11, B, N; 1953 PCHPDDMSchurPreType type = PC_HPDDM_SCHUR_PRE_LEAST_SQUARES; 1954 1955 PetscCall(MatSchurComplementGetSubMatrices(P, &A00, &P00, &A01, &A10, &A11)); 1956 PetscCall(PetscOptionsGetEnum(nullptr, pcpre, "-pc_hpddm_schur_precondition", PCHPDDMSchurPreTypes, (PetscEnum *)&type, &flg)); 1957 if (type == PC_HPDDM_SCHUR_PRE_LEAST_SQUARES) { 1958 Mat B01; 1959 Vec diagonal = nullptr; 1960 const PetscScalar *array; 1961 MatSchurComplementAinvType type; 1962 1963 PetscCall(PCHPDDMCheckSymmetry_Private(pc, A01, A10, &B01)); 1964 if (A11) { 1965 PetscCall(MatCreateVecs(A11, &diagonal, nullptr)); 1966 PetscCall(MatGetDiagonal(A11, diagonal)); 1967 } 1968 PetscCall(MatCreateVecs(P00, &v, nullptr)); 1969 PetscCall(MatSchurComplementGetAinvType(P, &type)); 1970 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]); 1971 if (type == MAT_SCHUR_COMPLEMENT_AINV_LUMP) { 1972 PetscCall(MatGetRowSum(P00, v)); 1973 if (A00 == P00) PetscCall(PetscObjectReference((PetscObject)A00)); 1974 PetscCall(MatDestroy(&P00)); 1975 PetscCall(VecGetArrayRead(v, &array)); 1976 PetscCall(MatCreateAIJ(PetscObjectComm((PetscObject)A00), A00->rmap->n, A00->cmap->n, A00->rmap->N, A00->cmap->N, 1, nullptr, 0, nullptr, &P00)); 1977 PetscCall(MatSetOption(P00, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE)); 1978 for (n = A00->rmap->rstart; n < A00->rmap->rend; ++n) PetscCall(MatSetValue(P00, n, n, array[n - A00->rmap->rstart], INSERT_VALUES)); 1979 PetscCall(MatAssemblyBegin(P00, MAT_FINAL_ASSEMBLY)); 1980 PetscCall(MatAssemblyEnd(P00, MAT_FINAL_ASSEMBLY)); 1981 PetscCall(VecRestoreArrayRead(v, &array)); 1982 PetscCall(MatSchurComplementUpdateSubMatrices(P, A00, P00, A01, A10, A11)); /* replace P00 by diag(sum of each row of P00) */ 1983 PetscCall(MatDestroy(&P00)); 1984 } else PetscCall(MatGetDiagonal(P00, v)); 1985 PetscCall(VecReciprocal(v)); /* inv(diag(P00)) */ 1986 PetscCall(VecSqrtAbs(v)); /* sqrt(inv(diag(P00))) */ 1987 PetscCall(MatDuplicate(A01, MAT_COPY_VALUES, &B)); 1988 PetscCall(MatDiagonalScale(B, v, nullptr)); 1989 if (B01) PetscCall(MatDiagonalScale(B01, v, nullptr)); 1990 PetscCall(VecDestroy(&v)); 1991 PetscCall(MatCreateNormalHermitian(B, &N)); 1992 PetscCall(PCHPDDMSetAuxiliaryMatNormal_Private(pc, B, N, &P, pcpre, &diagonal, B01)); 1993 PetscCall(PetscObjectTypeCompare((PetscObject)data->aux, MATSEQAIJ, &flg)); 1994 if (!flg) { 1995 PetscCall(MatDestroy(&P)); 1996 P = N; 1997 PetscCall(PetscObjectReference((PetscObject)P)); 1998 } 1999 if (diagonal) { 2000 PetscCall(MatDiagonalSet(P, diagonal, ADD_VALUES)); 2001 PetscCall(PCSetOperators(pc, P, P)); /* replace P by A01^T inv(diag(P00)) A01 - diag(P11) */ 2002 PetscCall(VecDestroy(&diagonal)); 2003 } else PetscCall(PCSetOperators(pc, B01 ? P : N, P)); /* replace P by A01^T inv(diag(P00)) A01 */ 2004 pc->ops->postsolve = PCPostSolve_SchurPreLeastSquares; /* PCFIELDSPLIT expect a KSP for (P11 - A10 inv(diag(P00)) A01) */ 2005 PetscCall(MatDestroy(&N)); /* but a PC for (A10 inv(diag(P00)) A10 - P11) is setup instead */ 2006 PetscCall(MatDestroy(&P)); /* so the sign of the solution must be flipped */ 2007 PetscCall(MatDestroy(&B)); 2008 } else 2009 PetscCheck(type != PC_HPDDM_SCHUR_PRE_GENEO, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-%spc_hpddm_schur_precondition %s without a prior call to PCHPDDMSetAuxiliaryMat() on the A11 block%s%s", pcpre ? pcpre : "", PCHPDDMSchurPreTypes[type], pcpre ? " -" : "", pcpre ? pcpre : ""); 2010 for (std::vector<Vec>::iterator it = initial.begin(); it != initial.end(); ++it) PetscCall(VecDestroy(&*it)); 2011 PetscFunctionReturn(PETSC_SUCCESS); 2012 } else { 2013 PetscCall(PetscOptionsGetString(nullptr, pcpre, "-pc_hpddm_levels_1_st_pc_type", type, sizeof(type), nullptr)); 2014 PetscCall(PetscStrcmp(type, PCMAT, &algebraic)); 2015 PetscCheck(!algebraic || !block, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-pc_hpddm_levels_1_st_pc_type mat and -pc_hpddm_block_splitting"); 2016 if (overlap != -1) { 2017 PetscCheck(!block && !algebraic, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-pc_hpddm_%s and -pc_hpddm_harmonic_overlap", block ? "block_splitting" : "levels_1_st_pc_type mat"); 2018 PetscCheck(overlap >= 1, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_WRONG, "-pc_hpddm_harmonic_overlap %" PetscInt_FMT " < 1", overlap); 2019 } 2020 if (block || overlap != -1) algebraic = PETSC_TRUE; 2021 if (algebraic) { 2022 PetscCall(ISCreateStride(PETSC_COMM_SELF, P->rmap->n, P->rmap->rstart, 1, &data->is)); 2023 PetscCall(MatIncreaseOverlap(P, 1, &data->is, 1)); 2024 PetscCall(ISSort(data->is)); 2025 } else 2026 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 and -%spc_hpddm_harmonic_overlap < 1\n", pcpre ? pcpre : "", pcpre ? pcpre : "", pcpre ? pcpre : "")); 2027 } 2028 } 2029 } 2030 } 2031 #if PetscDefined(USE_DEBUG) 2032 if (data->is) PetscCall(ISDuplicate(data->is, &dis)); 2033 if (data->aux) PetscCall(MatDuplicate(data->aux, MAT_COPY_VALUES, &daux)); 2034 #endif 2035 if (data->is || (ismatis && data->N > 1)) { 2036 if (ismatis) { 2037 std::initializer_list<std::string> list = {MATSEQBAIJ, MATSEQSBAIJ}; 2038 PetscCall(MatISGetLocalMat(P, &N)); 2039 std::initializer_list<std::string>::const_iterator it = std::find(list.begin(), list.end(), ((PetscObject)N)->type_name); 2040 PetscCall(MatISRestoreLocalMat(P, &N)); 2041 switch (std::distance(list.begin(), it)) { 2042 case 0: 2043 PetscCall(MatConvert(P, MATMPIBAIJ, MAT_INITIAL_MATRIX, &C)); 2044 break; 2045 case 1: 2046 /* MatCreateSubMatrices() does not work with MATSBAIJ and unsorted ISes, so convert to MPIBAIJ */ 2047 PetscCall(MatConvert(P, MATMPIBAIJ, MAT_INITIAL_MATRIX, &C)); 2048 PetscCall(MatSetOption(C, MAT_SYMMETRIC, PETSC_TRUE)); 2049 break; 2050 default: 2051 PetscCall(MatConvert(P, MATMPIAIJ, MAT_INITIAL_MATRIX, &C)); 2052 } 2053 PetscCall(MatISGetLocalToGlobalMapping(P, &l2g, nullptr)); 2054 PetscCall(PetscObjectReference((PetscObject)P)); 2055 PetscCall(KSPSetOperators(data->levels[0]->ksp, A, C)); 2056 std::swap(C, P); 2057 PetscCall(ISLocalToGlobalMappingGetSize(l2g, &n)); 2058 PetscCall(ISCreateStride(PETSC_COMM_SELF, n, 0, 1, &loc)); 2059 PetscCall(ISLocalToGlobalMappingApplyIS(l2g, loc, &is[0])); 2060 PetscCall(ISDestroy(&loc)); 2061 /* the auxiliary Mat is _not_ the local Neumann matrix */ 2062 /* it is the local Neumann matrix augmented (with zeros) through MatIncreaseOverlap() */ 2063 data->Neumann = PETSC_BOOL3_FALSE; 2064 structure = SAME_NONZERO_PATTERN; 2065 } else { 2066 is[0] = data->is; 2067 if (algebraic || ctx) subdomains = PETSC_TRUE; 2068 PetscCall(PetscOptionsGetBool(nullptr, pcpre, "-pc_hpddm_define_subdomains", &subdomains, nullptr)); 2069 if (ctx) PetscCheck(subdomains, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-%spc_hpddm_schur_precondition geneo and -%spc_hpddm_define_subdomains false", pcpre, pcpre); 2070 if (PetscBool3ToBool(data->Neumann)) { 2071 PetscCheck(!block, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-pc_hpddm_block_splitting and -pc_hpddm_has_neumann"); 2072 PetscCheck(overlap == -1, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-pc_hpddm_harmonic_overlap %" PetscInt_FMT " and -pc_hpddm_has_neumann", overlap); 2073 PetscCheck(!algebraic, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-pc_hpddm_levels_1_st_pc_type mat and -pc_hpddm_has_neumann"); 2074 } 2075 if (PetscBool3ToBool(data->Neumann) || block) structure = SAME_NONZERO_PATTERN; 2076 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)data->is), P->rmap->n, P->rmap->rstart, 1, &loc)); 2077 } 2078 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_levels_1_", pcpre ? pcpre : "")); 2079 PetscCall(PetscOptionsGetEnum(nullptr, prefix, "-st_matstructure", MatStructures, (PetscEnum *)&structure, &flg)); /* if not user-provided, force its value when possible */ 2080 if (!flg && structure == SAME_NONZERO_PATTERN) { /* cannot call STSetMatStructure() yet, insert the appropriate option in the database, parsed by STSetFromOptions() */ 2081 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-%spc_hpddm_levels_1_st_matstructure", pcpre ? pcpre : "")); 2082 PetscCall(PetscOptionsSetValue(nullptr, prefix, MatStructures[structure])); 2083 } 2084 flg = PETSC_FALSE; 2085 if (data->share) { 2086 data->share = PETSC_FALSE; /* will be reset to PETSC_TRUE if none of the conditions below are true */ 2087 if (!subdomains) PetscCall(PetscInfo(pc, "Cannot share subdomain KSP between SLEPc and PETSc since -%spc_hpddm_define_subdomains is not true\n", pcpre ? pcpre : "")); 2088 else if (data->deflation) PetscCall(PetscInfo(pc, "Nothing to share since PCHPDDMSetDeflationMat() has been called\n")); 2089 else if (ismatis) PetscCall(PetscInfo(pc, "Cannot share subdomain KSP between SLEPc and PETSc with a Pmat of type MATIS\n")); 2090 else if (!algebraic && structure != SAME_NONZERO_PATTERN) 2091 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])); 2092 else data->share = PETSC_TRUE; 2093 } 2094 if (!ismatis) { 2095 if (data->share || (!PetscBool3ToBool(data->Neumann) && subdomains)) PetscCall(ISDuplicate(is[0], &unsorted)); 2096 else unsorted = is[0]; 2097 } 2098 if ((ctx || data->N > 1) && (data->aux || ismatis || algebraic)) { 2099 PetscCheck(loadedSym, PETSC_COMM_SELF, PETSC_ERR_PLIB, "HPDDM library not loaded, cannot use more than one level"); 2100 PetscCall(MatSetOption(P, MAT_SUBMAT_SINGLEIS, PETSC_TRUE)); 2101 if (ismatis) { 2102 /* needed by HPDDM (currently) so that the partition of unity is 0 on subdomain interfaces */ 2103 PetscCall(MatIncreaseOverlap(P, 1, is, 1)); 2104 PetscCall(ISDestroy(&data->is)); 2105 data->is = is[0]; 2106 } else { 2107 if (PetscDefined(USE_DEBUG)) PetscCall(PCHPDDMCheckInclusion_Private(pc, data->is, loc, PETSC_TRUE)); 2108 if (!ctx && overlap == -1) PetscCall(PetscObjectComposeFunction((PetscObject)pc->pmat, "PCHPDDMAlgebraicAuxiliaryMat_Private_C", PCHPDDMAlgebraicAuxiliaryMat_Private)); 2109 if (!PetscBool3ToBool(data->Neumann) && (!algebraic || overlap != -1)) { 2110 PetscCall(PetscObjectTypeCompare((PetscObject)P, MATMPISBAIJ, &flg)); 2111 if (flg) { 2112 /* maybe better to ISSort(is[0]), MatCreateSubMatrices(), and then MatPermute() */ 2113 /* but there is no MatPermute_SeqSBAIJ(), so as before, just use MATMPIBAIJ */ 2114 PetscCall(MatConvert(P, MATMPIBAIJ, MAT_INITIAL_MATRIX, &uaux)); 2115 flg = PETSC_FALSE; 2116 } 2117 } 2118 } 2119 if (algebraic && overlap == -1) { 2120 PetscUseMethod(pc->pmat, "PCHPDDMAlgebraicAuxiliaryMat_Private_C", (Mat, IS *, Mat *[], PetscBool), (P, is, &sub, block)); 2121 if (block) { 2122 PetscCall(PetscObjectQuery((PetscObject)sub[0], "_PCHPDDM_Neumann_Mat", (PetscObject *)&data->aux)); 2123 PetscCall(PetscObjectCompose((PetscObject)sub[0], "_PCHPDDM_Neumann_Mat", nullptr)); 2124 } 2125 } else if (!uaux || overlap != -1) { 2126 if (!ctx) { 2127 if (PetscBool3ToBool(data->Neumann)) sub = &data->aux; 2128 else { 2129 PetscBool flg; 2130 if (overlap != -1) { 2131 Harmonic h; 2132 Mat A0, *a; /* with an SVD: [ A_00 A_01 ] */ 2133 IS ov[2], rows, cols, stride; /* [ A_10 A_11 A_12 ] */ 2134 const PetscInt *i[2], bs = P->cmap->bs; /* with a GEVP: [ A_00 A_01 ] */ 2135 PetscInt n[2]; /* [ A_10 A_11 A_12 ] */ 2136 std::vector<PetscInt> v[2]; /* [ A_21 A_22 ] */ 2137 2138 do { 2139 PetscCall(ISDuplicate(data->is, ov)); 2140 if (overlap > 1) PetscCall(MatIncreaseOverlap(P, 1, ov, overlap - 1)); 2141 PetscCall(ISDuplicate(ov[0], ov + 1)); 2142 PetscCall(MatIncreaseOverlap(P, 1, ov + 1, 1)); 2143 PetscCall(ISGetLocalSize(ov[0], n)); 2144 PetscCall(ISGetLocalSize(ov[1], n + 1)); 2145 flg = PetscBool(n[0] == n[1]); 2146 PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &flg, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)pc))); 2147 if (flg) { 2148 PetscCall(ISDestroy(ov)); 2149 PetscCall(ISDestroy(ov + 1)); 2150 PetscCheck(--overlap, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "No oversampling possible"); 2151 PetscCall(PetscInfo(pc, "Supplied -%spc_hpddm_harmonic_overlap parameter is too large, it has been decreased to %" PetscInt_FMT "\n", pcpre ? pcpre : "", overlap)); 2152 } else break; 2153 } while (1); 2154 PetscCall(PetscNew(&h)); 2155 h->ksp = nullptr; 2156 PetscCall(PetscCalloc1(2, &h->A)); 2157 PetscCall(PetscOptionsHasName(nullptr, prefix, "-svd_nsv", &flg)); 2158 if (!flg) PetscCall(PetscOptionsHasName(nullptr, prefix, "-svd_relative_threshold", &flg)); 2159 PetscCall(ISSort(ov[0])); 2160 if (!flg) PetscCall(ISSort(ov[1])); 2161 PetscCall(PetscCalloc1(5, &h->is)); 2162 PetscCall(MatCreateSubMatrices(uaux ? uaux : P, 1, ov + !flg, ov + 1, MAT_INITIAL_MATRIX, &a)); /* submatrix from above, either square (!flg) or rectangular (flg) */ 2163 for (PetscInt j = 0; j < 2; ++j) PetscCall(ISGetIndices(ov[j], i + j)); 2164 v[1].reserve((n[1] - n[0]) / bs); 2165 for (PetscInt j = 0; j < n[1]; j += bs) { /* indices of the (2,2) block */ 2166 PetscInt location; 2167 PetscCall(ISLocate(ov[0], i[1][j], &location)); 2168 if (location < 0) v[1].emplace_back(j / bs); 2169 } 2170 if (!flg) { 2171 h->A[1] = a[0]; 2172 PetscCall(PetscObjectReference((PetscObject)h->A[1])); 2173 v[0].reserve((n[0] - P->rmap->n) / bs); 2174 for (PetscInt j = 0; j < n[1]; j += bs) { /* row indices of the (1,2) block */ 2175 PetscInt location; 2176 PetscCall(ISLocate(loc, i[1][j], &location)); 2177 if (location < 0) { 2178 PetscCall(ISLocate(ov[0], i[1][j], &location)); 2179 if (location >= 0) v[0].emplace_back(j / bs); 2180 } 2181 } 2182 PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, v[0].size(), v[0].data(), PETSC_USE_POINTER, &rows)); 2183 PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, v[1].size(), v[1].data(), PETSC_COPY_VALUES, h->is + 4)); 2184 PetscCall(MatCreateSubMatrix(a[0], rows, h->is[4], MAT_INITIAL_MATRIX, h->A)); /* A_12 submatrix from above */ 2185 PetscCall(ISDestroy(&rows)); 2186 if (uaux) PetscCall(MatConvert(a[0], MATSEQSBAIJ, MAT_INPLACE_MATRIX, a)); /* initial Pmat was MATSBAIJ, convert back to the same format since the rectangular A_12 submatrix has been created */ 2187 PetscCall(ISEmbed(ov[0], ov[1], PETSC_TRUE, &rows)); 2188 PetscCall(MatCreateSubMatrix(a[0], rows, cols = rows, MAT_INITIAL_MATRIX, &A0)); /* [ A_00 A_01 ; A_10 A_11 ] submatrix from above */ 2189 PetscCall(ISDestroy(&rows)); 2190 v[0].clear(); 2191 PetscCall(ISEmbed(loc, ov[1], PETSC_TRUE, h->is + 3)); 2192 PetscCall(ISEmbed(data->is, ov[1], PETSC_TRUE, h->is + 2)); 2193 } 2194 v[0].reserve((n[0] - P->rmap->n) / bs); 2195 for (PetscInt j = 0; j < n[0]; j += bs) { 2196 PetscInt location; 2197 PetscCall(ISLocate(loc, i[0][j], &location)); 2198 if (location < 0) v[0].emplace_back(j / bs); 2199 } 2200 PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, v[0].size(), v[0].data(), PETSC_USE_POINTER, &rows)); 2201 for (PetscInt j = 0; j < 2; ++j) PetscCall(ISRestoreIndices(ov[j], i + j)); 2202 if (flg) { 2203 IS is; 2204 PetscCall(ISCreateStride(PETSC_COMM_SELF, a[0]->rmap->n, 0, 1, &is)); 2205 PetscCall(ISEmbed(ov[0], ov[1], PETSC_TRUE, &cols)); 2206 PetscCall(MatCreateSubMatrix(a[0], is, cols, MAT_INITIAL_MATRIX, &A0)); /* [ A_00 A_01 ; A_10 A_11 ] submatrix from above */ 2207 PetscCall(ISDestroy(&cols)); 2208 PetscCall(ISDestroy(&is)); 2209 if (uaux) PetscCall(MatConvert(A0, MATSEQSBAIJ, MAT_INPLACE_MATRIX, &A0)); /* initial Pmat was MATSBAIJ, convert back to the same format since this submatrix is square */ 2210 PetscCall(ISEmbed(loc, data->is, PETSC_TRUE, h->is + 2)); 2211 PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, v[1].size(), v[1].data(), PETSC_USE_POINTER, &cols)); 2212 PetscCall(MatCreateSubMatrix(a[0], rows, cols, MAT_INITIAL_MATRIX, h->A)); /* A_12 submatrix from above */ 2213 PetscCall(ISDestroy(&cols)); 2214 } 2215 PetscCall(ISCreateStride(PETSC_COMM_SELF, A0->rmap->n, 0, 1, &stride)); 2216 PetscCall(ISEmbed(rows, stride, PETSC_TRUE, h->is)); 2217 PetscCall(ISDestroy(&stride)); 2218 PetscCall(ISDestroy(&rows)); 2219 PetscCall(ISEmbed(loc, ov[0], PETSC_TRUE, h->is + 1)); 2220 if (subdomains) { 2221 if (!data->levels[0]->pc) { 2222 PetscCall(PCCreate(PetscObjectComm((PetscObject)pc), &data->levels[0]->pc)); 2223 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_levels_1_", pcpre ? pcpre : "")); 2224 PetscCall(PCSetOptionsPrefix(data->levels[0]->pc, prefix)); 2225 PetscCall(PCSetOperators(data->levels[0]->pc, A, P)); 2226 } 2227 PetscCall(PCSetType(data->levels[0]->pc, PCASM)); 2228 if (!data->levels[0]->pc->setupcalled) PetscCall(PCASMSetLocalSubdomains(data->levels[0]->pc, 1, ov + !flg, &loc)); 2229 PetscCall(PCHPDDMCommunicationAvoidingPCASM_Private(data->levels[0]->pc, flg ? A0 : a[0], PETSC_TRUE)); 2230 if (!flg) ++overlap; 2231 if (data->share) { 2232 PetscInt n = -1; 2233 PetscTryMethod(data->levels[0]->pc, "PCASMGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (data->levels[0]->pc, &n, nullptr, &ksp)); 2234 PetscCheck(n == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of subdomain solver %" PetscInt_FMT " != 1", n); 2235 if (flg) { 2236 h->ksp = ksp[0]; 2237 PetscCall(PetscObjectReference((PetscObject)h->ksp)); 2238 } 2239 } 2240 } 2241 if (!h->ksp) { 2242 PetscBool share = data->share; 2243 PetscCall(KSPCreate(PETSC_COMM_SELF, &h->ksp)); 2244 PetscCall(KSPSetType(h->ksp, KSPPREONLY)); 2245 PetscCall(KSPSetOperators(h->ksp, A0, A0)); 2246 do { 2247 if (!data->share) { 2248 share = PETSC_FALSE; 2249 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_levels_1_%s", pcpre ? pcpre : "", flg ? "svd_" : "eps_")); 2250 PetscCall(KSPSetOptionsPrefix(h->ksp, prefix)); 2251 PetscCall(KSPSetFromOptions(h->ksp)); 2252 } else { 2253 MatSolverType type; 2254 PetscCall(KSPGetPC(ksp[0], &pc)); 2255 PetscCall(PetscObjectTypeCompareAny((PetscObject)pc, &data->share, PCLU, PCCHOLESKY, "")); 2256 if (data->share) { 2257 PetscCall(PCFactorGetMatSolverType(pc, &type)); 2258 if (!type) { 2259 if (PetscDefined(HAVE_MUMPS)) PetscCall(PCFactorSetMatSolverType(pc, MATSOLVERMUMPS)); 2260 else if (PetscDefined(HAVE_MKL_PARDISO)) PetscCall(PCFactorSetMatSolverType(pc, MATSOLVERMKL_PARDISO)); 2261 else data->share = PETSC_FALSE; 2262 if (data->share) PetscCall(PCSetFromOptions(pc)); 2263 } else { 2264 PetscCall(PetscStrcmp(type, MATSOLVERMUMPS, &data->share)); 2265 if (!data->share) PetscCall(PetscStrcmp(type, MATSOLVERMKL_PARDISO, &data->share)); 2266 } 2267 if (data->share) { 2268 std::tuple<KSP, IS, Vec[2]> *p; 2269 PetscCall(PCFactorGetMatrix(pc, &A)); 2270 PetscCall(MatFactorSetSchurIS(A, h->is[4])); 2271 PetscCall(KSPSetUp(ksp[0])); 2272 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_levels_1_eps_shell_", pcpre ? pcpre : "")); 2273 PetscCall(KSPSetOptionsPrefix(h->ksp, prefix)); 2274 PetscCall(KSPSetFromOptions(h->ksp)); 2275 PetscCall(KSPGetPC(h->ksp, &pc)); 2276 PetscCall(PCSetType(pc, PCSHELL)); 2277 PetscCall(PetscNew(&p)); 2278 std::get<0>(*p) = ksp[0]; 2279 PetscCall(ISEmbed(ov[0], ov[1], PETSC_TRUE, &std::get<1>(*p))); 2280 PetscCall(MatCreateVecs(A, std::get<2>(*p), std::get<2>(*p) + 1)); 2281 PetscCall(PCShellSetContext(pc, p)); 2282 PetscCall(PCShellSetApply(pc, PCApply_Schur)); 2283 PetscCall(PCShellSetApplyTranspose(pc, PCApply_Schur<Vec, true>)); 2284 PetscCall(PCShellSetMatApply(pc, PCApply_Schur<Mat>)); 2285 PetscCall(PCShellSetDestroy(pc, PCDestroy_Schur)); 2286 } 2287 } 2288 if (!data->share) PetscCall(PetscInfo(pc, "Cannot share subdomain KSP between SLEPc and PETSc since neither MUMPS nor MKL PARDISO is used\n")); 2289 } 2290 } while (!share != !data->share); /* if data->share is initially PETSC_TRUE, but then reset to PETSC_FALSE, then go back to the beginning of the do loop */ 2291 } 2292 PetscCall(ISDestroy(ov)); 2293 PetscCall(ISDestroy(ov + 1)); 2294 if (overlap == 1 && subdomains && flg) { 2295 *subA = A0; 2296 sub = subA; 2297 if (uaux) PetscCall(MatDestroy(&uaux)); 2298 } else PetscCall(MatDestroy(&A0)); 2299 PetscCall(MatCreateShell(PETSC_COMM_SELF, P->rmap->n, n[1] - n[0], P->rmap->n, n[1] - n[0], h, &data->aux)); 2300 PetscCall(KSPSetErrorIfNotConverged(h->ksp, PETSC_TRUE)); /* bail out as early as possible to avoid (apparently) unrelated error messages */ 2301 PetscCall(MatCreateVecs(h->ksp->pc->pmat, &h->v, nullptr)); 2302 PetscCall(MatShellSetOperation(data->aux, MATOP_MULT, (void (*)(void))MatMult_Harmonic)); 2303 PetscCall(MatShellSetOperation(data->aux, MATOP_MULT_TRANSPOSE, (void (*)(void))MatMultTranspose_Harmonic)); 2304 PetscCall(MatShellSetMatProductOperation(data->aux, MATPRODUCT_AB, nullptr, MatProduct_AB_Harmonic, nullptr, MATDENSE, MATDENSE)); 2305 PetscCall(MatShellSetOperation(data->aux, MATOP_DESTROY, (void (*)(void))MatDestroy_Harmonic)); 2306 PetscCall(MatDestroySubMatrices(1, &a)); 2307 } 2308 if (overlap != 1 || !subdomains) { 2309 PetscCall(MatCreateSubMatrices(uaux ? uaux : P, 1, is, is, MAT_INITIAL_MATRIX, &sub)); 2310 if (ismatis) { 2311 PetscCall(MatISGetLocalMat(C, &N)); 2312 PetscCall(PetscObjectTypeCompare((PetscObject)N, MATSEQSBAIJ, &flg)); 2313 if (flg) PetscCall(MatConvert(sub[0], MATSEQSBAIJ, MAT_INPLACE_MATRIX, sub)); 2314 PetscCall(MatISRestoreLocalMat(C, &N)); 2315 } 2316 } 2317 if (uaux) { 2318 PetscCall(MatDestroy(&uaux)); 2319 PetscCall(MatConvert(sub[0], MATSEQSBAIJ, MAT_INPLACE_MATRIX, sub)); 2320 } 2321 } 2322 } 2323 } else if (!ctx) { 2324 PetscCall(MatCreateSubMatrices(uaux, 1, is, is, MAT_INITIAL_MATRIX, &sub)); 2325 PetscCall(MatDestroy(&uaux)); 2326 PetscCall(MatConvert(sub[0], MATSEQSBAIJ, MAT_INPLACE_MATRIX, sub)); 2327 } 2328 if (data->N > 1) { 2329 /* Vec holding the partition of unity */ 2330 if (!data->levels[0]->D) { 2331 PetscCall(ISGetLocalSize(data->is, &n)); 2332 PetscCall(VecCreateMPI(PETSC_COMM_SELF, n, PETSC_DETERMINE, &data->levels[0]->D)); 2333 } 2334 if (data->share && overlap == -1) { 2335 Mat D; 2336 IS perm = nullptr; 2337 PetscInt size = -1; 2338 2339 if (!data->levels[0]->pc) { 2340 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_levels_1_", pcpre ? pcpre : "")); 2341 PetscCall(PCCreate(PetscObjectComm((PetscObject)pc), &data->levels[0]->pc)); 2342 PetscCall(PCSetOptionsPrefix(data->levels[0]->pc, prefix)); 2343 PetscCall(PCSetOperators(data->levels[0]->pc, A, P)); 2344 } 2345 PetscCall(PCSetType(data->levels[0]->pc, PCASM)); 2346 if (!ctx) { 2347 if (!data->levels[0]->pc->setupcalled) { 2348 IS sorted; /* PCASM will sort the input IS, duplicate it to return an unmodified (PCHPDDM) input IS */ 2349 PetscCall(ISDuplicate(is[0], &sorted)); 2350 PetscCall(PCASMSetLocalSubdomains(data->levels[0]->pc, 1, &sorted, &loc)); 2351 PetscCall(PetscObjectDereference((PetscObject)sorted)); 2352 } 2353 PetscCall(PCSetFromOptions(data->levels[0]->pc)); 2354 if (block) { 2355 PetscCall(PCHPDDMPermute_Private(unsorted, data->is, &uis, sub[0], &C, &perm)); 2356 PetscCall(PCHPDDMCommunicationAvoidingPCASM_Private(data->levels[0]->pc, C, algebraic)); 2357 } else PetscCall(PCSetUp(data->levels[0]->pc)); 2358 PetscTryMethod(data->levels[0]->pc, "PCASMGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (data->levels[0]->pc, &size, nullptr, &ksp)); 2359 if (size != 1) { 2360 data->share = PETSC_FALSE; 2361 PetscCheck(size == -1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of subdomain solver %" PetscInt_FMT " != 1", size); 2362 PetscCall(PetscInfo(pc, "Cannot share subdomain KSP between SLEPc and PETSc since PCASMGetSubKSP() not found in fine-level PC\n")); 2363 PetscCall(ISDestroy(&unsorted)); 2364 unsorted = is[0]; 2365 } else { 2366 const char *matpre; 2367 PetscBool cmp[4]; 2368 2369 if (!block && !ctx) PetscCall(PCHPDDMPermute_Private(unsorted, data->is, &uis, PetscBool3ToBool(data->Neumann) ? sub[0] : data->aux, &C, &perm)); 2370 if (perm) { /* unsorted input IS */ 2371 if (!PetscBool3ToBool(data->Neumann) && !block) { 2372 PetscCall(MatPermute(sub[0], perm, perm, &D)); /* permute since PCASM will call ISSort() */ 2373 PetscCall(MatHeaderReplace(sub[0], &D)); 2374 } 2375 if (data->B) { /* see PCHPDDMSetRHSMat() */ 2376 PetscCall(MatPermute(data->B, perm, perm, &D)); 2377 PetscCall(MatHeaderReplace(data->B, &D)); 2378 } 2379 PetscCall(ISDestroy(&perm)); 2380 } 2381 PetscCall(KSPGetOperators(ksp[0], subA, subA + 1)); 2382 PetscCall(PetscObjectReference((PetscObject)subA[0])); 2383 PetscCall(MatDuplicate(subA[1], MAT_SHARE_NONZERO_PATTERN, &D)); 2384 PetscCall(MatGetOptionsPrefix(subA[1], &matpre)); 2385 PetscCall(MatSetOptionsPrefix(D, matpre)); 2386 PetscCall(PetscObjectTypeCompare((PetscObject)D, MATNORMAL, cmp)); 2387 PetscCall(PetscObjectTypeCompare((PetscObject)C, MATNORMAL, cmp + 1)); 2388 if (!cmp[0]) PetscCall(PetscObjectTypeCompare((PetscObject)D, MATNORMALHERMITIAN, cmp + 2)); 2389 else cmp[2] = PETSC_FALSE; 2390 if (!cmp[1]) PetscCall(PetscObjectTypeCompare((PetscObject)C, MATNORMALHERMITIAN, cmp + 3)); 2391 else cmp[3] = PETSC_FALSE; 2392 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); 2393 if (!cmp[0] && !cmp[2]) { 2394 if (!block) PetscCall(MatAXPY(D, 1.0, C, SUBSET_NONZERO_PATTERN)); 2395 else { 2396 PetscCall(MatMissingDiagonal(D, cmp, nullptr)); 2397 if (cmp[0]) structure = DIFFERENT_NONZERO_PATTERN; /* data->aux has no missing diagonal entry */ 2398 PetscCall(MatAXPY(D, 1.0, data->aux, structure)); 2399 } 2400 } else { 2401 Mat mat[2]; 2402 2403 if (cmp[0]) { 2404 PetscCall(MatNormalGetMat(D, mat)); 2405 PetscCall(MatNormalGetMat(C, mat + 1)); 2406 } else { 2407 PetscCall(MatNormalHermitianGetMat(D, mat)); 2408 PetscCall(MatNormalHermitianGetMat(C, mat + 1)); 2409 } 2410 PetscCall(MatAXPY(mat[0], 1.0, mat[1], SUBSET_NONZERO_PATTERN)); 2411 } 2412 PetscCall(MatPropagateSymmetryOptions(C, D)); 2413 PetscCall(MatDestroy(&C)); 2414 C = D; 2415 /* swap pointers so that variables stay consistent throughout PCSetUp() */ 2416 std::swap(C, data->aux); 2417 std::swap(uis, data->is); 2418 swap = PETSC_TRUE; 2419 } 2420 } 2421 } 2422 } 2423 if (ctx) { 2424 PC_HPDDM *data_00 = (PC_HPDDM *)std::get<0>(*ctx)[0]->data; 2425 PC s; 2426 Mat A00, P00, A01 = nullptr, A10, A11, N, b[4]; 2427 IS sorted, is[2], *is_00; 2428 MatSolverType type; 2429 std::pair<PC, Vec[2]> *p; 2430 2431 n = -1; 2432 PetscTryMethod(data_00->levels[0]->pc, "PCASMGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (data_00->levels[0]->pc, &n, nullptr, &ksp)); 2433 PetscCheck(n == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of subdomain solver %" PetscInt_FMT " != 1", n); 2434 PetscCall(KSPGetOperators(ksp[0], subA, subA + 1)); 2435 PetscCall(ISGetLocalSize(data_00->is, &n)); 2436 if (n != subA[0]->rmap->n || n != subA[0]->cmap->n) { 2437 PetscCall(PCASMGetLocalSubdomains(data_00->levels[0]->pc, &n, &is_00, nullptr)); 2438 PetscCall(ISGetLocalSize(*is_00, &n)); 2439 PetscCheck(n == subA[0]->rmap->n && n == subA[0]->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "-%spc_hpddm_schur_precondition geneo and -%spc_hpddm_define_subdomains false", pcpre ? pcpre : "", ((PetscObject)pc)->prefix); 2440 } else is_00 = &data_00->is; 2441 PetscCall(PCHPDDMPermute_Private(unsorted, data->is, &uis, data->aux, &C, nullptr)); /* permute since PCASM works with a sorted IS */ 2442 std::swap(C, data->aux); 2443 std::swap(uis, data->is); 2444 swap = PETSC_TRUE; 2445 PetscCall(MatSchurComplementGetSubMatrices(P, &A00, &P00, std::get<1>(*ctx), &A10, &A11)); 2446 std::get<1>(*ctx)[1] = A10; 2447 PetscCall(PetscObjectTypeCompare((PetscObject)A10, MATTRANSPOSEVIRTUAL, &flg)); 2448 if (flg) PetscCall(MatTransposeGetMat(A10, &A01)); 2449 else { 2450 PetscBool flg; 2451 2452 PetscCall(PetscObjectTypeCompare((PetscObject)A10, MATHERMITIANTRANSPOSEVIRTUAL, &flg)); 2453 if (flg) PetscCall(MatHermitianTransposeGetMat(A10, &A01)); 2454 } 2455 PetscCall(ISDuplicate(*is_00, &sorted)); /* during setup of the PC associated to the A00 block, this IS has already been sorted, but it's put back to its original state at the end of PCSetUp_HPDDM(), which may be unsorted */ 2456 PetscCall(ISSort(sorted)); /* this is to avoid changing users inputs, but it requires a new call to ISSort() here */ 2457 if (!A01) { 2458 PetscCall(MatSetOption(A10, MAT_SUBMAT_SINGLEIS, PETSC_TRUE)); 2459 PetscCall(MatCreateSubMatrices(A10, 1, &data->is, &sorted, MAT_INITIAL_MATRIX, &sub)); 2460 b[2] = sub[0]; 2461 PetscCall(PetscObjectReference((PetscObject)sub[0])); 2462 PetscCall(MatDestroySubMatrices(1, &sub)); 2463 PetscCall(PetscObjectTypeCompare((PetscObject)std::get<1>(*ctx)[0], MATTRANSPOSEVIRTUAL, &flg)); 2464 A10 = nullptr; 2465 if (flg) PetscCall(MatTransposeGetMat(std::get<1>(*ctx)[0], &A10)); 2466 else { 2467 PetscBool flg; 2468 2469 PetscCall(PetscObjectTypeCompare((PetscObject)std::get<1>(*ctx)[0], MATHERMITIANTRANSPOSEVIRTUAL, &flg)); 2470 if (flg) PetscCall(MatHermitianTransposeGetMat(std::get<1>(*ctx)[0], &A10)); 2471 } 2472 if (!A10) PetscCall(MatCreateSubMatrices(std::get<1>(*ctx)[0], 1, &sorted, &data->is, MAT_INITIAL_MATRIX, &sub)); 2473 else { 2474 if (flg) PetscCall(MatCreateTranspose(b[2], b + 1)); 2475 else PetscCall(MatCreateHermitianTranspose(b[2], b + 1)); 2476 } 2477 } else { 2478 PetscCall(MatSetOption(A01, MAT_SUBMAT_SINGLEIS, PETSC_TRUE)); 2479 PetscCall(MatCreateSubMatrices(A01, 1, &sorted, &data->is, MAT_INITIAL_MATRIX, &sub)); 2480 if (flg) PetscCall(MatCreateTranspose(*sub, b + 2)); 2481 else PetscCall(MatCreateHermitianTranspose(*sub, b + 2)); 2482 } 2483 if (A01 || !A10) { 2484 b[1] = sub[0]; 2485 PetscCall(PetscObjectReference((PetscObject)sub[0])); 2486 } 2487 PetscCall(MatDestroySubMatrices(1, &sub)); 2488 PetscCall(ISDestroy(&sorted)); 2489 b[3] = data->aux; 2490 PetscCall(MatCreateSchurComplement(subA[0], subA[1], b[1], b[2], b[3], &S)); 2491 PetscCall(MatSchurComplementSetKSP(S, ksp[0])); 2492 if (data->N != 1) { 2493 PetscCall(PCASMSetType(data->levels[0]->pc, PC_ASM_NONE)); /* "Neumann--Neumann" preconditioning with overlap and a Boolean partition of unity */ 2494 PetscCall(PCASMSetLocalSubdomains(data->levels[0]->pc, 1, &data->is, &loc)); 2495 PetscCall(PCSetFromOptions(data->levels[0]->pc)); /* action of eq. (15) of https://hal.science/hal-02343808v6/document (with a sign flip) */ 2496 s = data->levels[0]->pc; 2497 } else { 2498 is[0] = data->is; 2499 PetscCall(PetscObjectReference((PetscObject)is[0])); 2500 PetscCall(PetscObjectReference((PetscObject)b[3])); 2501 PetscCall(PCSetType(pc, PCASM)); /* change the type of the current PC */ 2502 data = nullptr; /* destroyed in the previous PCSetType(), so reset to NULL to avoid any faulty use */ 2503 PetscCall(PCAppendOptionsPrefix(pc, "pc_hpddm_coarse_")); /* same prefix as when using PCHPDDM with a single level */ 2504 PetscCall(PCASMSetLocalSubdomains(pc, 1, is, &loc)); 2505 PetscCall(ISDestroy(is)); 2506 PetscCall(ISDestroy(&loc)); 2507 s = pc; 2508 } 2509 PetscCall(PCHPDDMCommunicationAvoidingPCASM_Private(s, S, PETSC_TRUE)); /* the subdomain Mat is already known and the input IS of PCASMSetLocalSubdomains() is already sorted */ 2510 PetscTryMethod(s, "PCASMGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (s, &n, nullptr, &ksp)); 2511 PetscCheck(n == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of subdomain solver %" PetscInt_FMT " != 1", n); 2512 PetscCall(KSPGetPC(ksp[0], &inner)); 2513 PetscCall(PCSetType(inner, PCSHELL)); /* compute the action of the inverse of the local Schur complement with a PCSHELL */ 2514 b[0] = subA[0]; 2515 PetscCall(MatCreateNest(PETSC_COMM_SELF, 2, nullptr, 2, nullptr, b, &N)); /* instead of computing inv(A11 - A10 inv(A00) A01), compute inv([A00, A01; A10, A11]) followed by a partial solution associated to the A11 block */ 2516 if (!data) PetscCall(PetscObjectDereference((PetscObject)b[3])); 2517 PetscCall(PetscObjectDereference((PetscObject)b[1])); 2518 PetscCall(PetscObjectDereference((PetscObject)b[2])); 2519 PetscCall(PCCreate(PETSC_COMM_SELF, &s)); 2520 PetscCall(PCSetOptionsPrefix(s, ((PetscObject)inner)->prefix)); 2521 PetscCall(PCSetOptionsPrefix(inner, nullptr)); 2522 PetscCall(KSPSetSkipPCSetFromOptions(ksp[0], PETSC_TRUE)); 2523 PetscCall(PCSetType(s, PCLU)); 2524 if (PetscDefined(HAVE_MUMPS)) PetscCall(PCFactorSetMatSolverType(s, MATSOLVERMUMPS)); /* only MATSOLVERMUMPS handles MATNEST, so for the others, e.g., MATSOLVERPETSC or MATSOLVERMKL_PARDISO, convert to plain MATAIJ */ 2525 PetscCall(PCSetFromOptions(s)); 2526 PetscCall(PCFactorGetMatSolverType(s, &type)); 2527 PetscCall(PetscStrcmp(type, MATSOLVERMUMPS, &flg)); 2528 PetscCall(MatGetLocalSize(A11, &n, nullptr)); 2529 if (flg || n == 0) { 2530 PetscCall(PCSetOperators(s, N, N)); 2531 if (n) { 2532 PetscCall(PCFactorGetMatrix(s, b)); 2533 PetscCall(MatSetOptionsPrefix(*b, ((PetscObject)s)->prefix)); 2534 n = -1; 2535 PetscCall(PetscOptionsGetInt(nullptr, ((PetscObject)s)->prefix, "-mat_mumps_icntl_26", &n, nullptr)); 2536 if (n == 1) { /* allocates a square MatDense of size is[1]->map->n, so one */ 2537 PetscCall(MatNestGetISs(N, is, nullptr)); /* needs to be able to deactivate this path when dealing */ 2538 PetscCall(MatFactorSetSchurIS(*b, is[1])); /* with a large constraint space in order to avoid OOM */ 2539 } 2540 } else PetscCall(PCSetType(s, PCNONE)); /* empty local Schur complement (e.g., centralized on another process) */ 2541 } else { 2542 PetscCall(MatConvert(N, MATAIJ, MAT_INITIAL_MATRIX, b)); 2543 PetscCall(PCSetOperators(s, N, *b)); 2544 PetscCall(PetscObjectDereference((PetscObject)*b)); 2545 PetscCall(PetscObjectTypeCompareAny((PetscObject)s, &flg, PCLU, PCCHOLESKY, PCILU, PCICC, PCQR, "")); 2546 if (flg) PetscCall(PCFactorGetMatrix(s, b)); /* MATSOLVERMKL_PARDISO cannot compute in PETSc (yet) a partial solution associated to the A11 block, only partial solution associated to the A00 block or full solution */ 2547 } 2548 PetscCall(PetscNew(&p)); 2549 p->first = s; 2550 if (n != 0) PetscCall(MatCreateVecs(*b, p->second, p->second + 1)); 2551 else p->second[0] = p->second[1] = nullptr; 2552 PetscCall(PCShellSetContext(inner, p)); 2553 PetscCall(PCShellSetApply(inner, PCApply_Nest)); 2554 PetscCall(PCShellSetView(inner, PCView_Nest)); 2555 PetscCall(PCShellSetDestroy(inner, PCDestroy_Nest)); 2556 PetscCall(PetscObjectDereference((PetscObject)N)); 2557 if (!data) { 2558 PetscCall(MatDestroy(&S)); 2559 PetscCall(ISDestroy(&unsorted)); 2560 PetscCall(MatDestroy(&C)); 2561 PetscCall(ISDestroy(&uis)); 2562 PetscCall(PetscFree(ctx)); 2563 #if PetscDefined(USE_DEBUG) 2564 PetscCall(ISDestroy(&dis)); 2565 PetscCall(MatDestroy(&daux)); 2566 #endif 2567 PetscFunctionReturn(PETSC_SUCCESS); 2568 } 2569 } 2570 if (!data->levels[0]->scatter) { 2571 PetscCall(MatCreateVecs(P, &xin, nullptr)); 2572 if (ismatis) PetscCall(MatDestroy(&P)); 2573 PetscCall(VecScatterCreate(xin, data->is, data->levels[0]->D, nullptr, &data->levels[0]->scatter)); 2574 PetscCall(VecDestroy(&xin)); 2575 } 2576 if (data->levels[0]->P) { 2577 /* if the pattern is the same and PCSetUp() has previously succeeded, reuse HPDDM buffers and connectivity */ 2578 PetscCall(HPDDM::Schwarz<PetscScalar>::destroy(data->levels[0], pc->setupcalled < 1 || pc->flag == DIFFERENT_NONZERO_PATTERN ? PETSC_TRUE : PETSC_FALSE)); 2579 } 2580 if (!data->levels[0]->P) data->levels[0]->P = new HPDDM::Schwarz<PetscScalar>(); 2581 if (data->log_separate) PetscCall(PetscLogEventBegin(PC_HPDDM_SetUp[0], data->levels[0]->ksp, nullptr, nullptr, nullptr)); 2582 else PetscCall(PetscLogEventBegin(PC_HPDDM_Strc, data->levels[0]->ksp, nullptr, nullptr, nullptr)); 2583 /* HPDDM internal data structure */ 2584 PetscCall(data->levels[0]->P->structure(loc, data->is, !ctx ? sub[0] : nullptr, ismatis ? C : data->aux, data->levels)); 2585 if (!data->log_separate) PetscCall(PetscLogEventEnd(PC_HPDDM_Strc, data->levels[0]->ksp, nullptr, nullptr, nullptr)); 2586 /* matrix pencil of the generalized eigenvalue problem on the overlap (GenEO) */ 2587 if (!ctx) { 2588 if (data->deflation || overlap != -1) weighted = data->aux; 2589 else if (!data->B) { 2590 PetscBool cmp; 2591 2592 PetscCall(MatDuplicate(sub[0], MAT_COPY_VALUES, &weighted)); 2593 PetscCall(PetscObjectTypeCompareAny((PetscObject)weighted, &cmp, MATNORMAL, MATNORMALHERMITIAN, "")); 2594 if (cmp) flg = PETSC_FALSE; 2595 PetscCall(MatDiagonalScale(weighted, data->levels[0]->D, data->levels[0]->D)); 2596 /* neither MatDuplicate() nor MatDiagonalScale() handles the symmetry options, so propagate the options explicitly */ 2597 /* only useful for -mat_type baij -pc_hpddm_levels_1_st_pc_type cholesky (no problem with MATAIJ or MATSBAIJ) */ 2598 PetscCall(MatPropagateSymmetryOptions(sub[0], weighted)); 2599 } else weighted = data->B; 2600 } else weighted = nullptr; 2601 /* SLEPc is used inside the loaded symbol */ 2602 PetscCall((*loadedSym)(data->levels[0]->P, data->is, ismatis ? C : (algebraic && !block && overlap == -1 ? sub[0] : (!ctx ? data->aux : S)), weighted, data->B, initial, data->levels)); 2603 if (!ctx && data->share && overlap == -1) { 2604 Mat st[2]; 2605 PetscCall(KSPGetOperators(ksp[0], st, st + 1)); 2606 PetscCall(MatCopy(subA[0], st[0], structure)); 2607 if (subA[1] != subA[0] || st[1] != st[0]) PetscCall(MatCopy(subA[1], st[1], SAME_NONZERO_PATTERN)); 2608 PetscCall(PetscObjectDereference((PetscObject)subA[0])); 2609 } 2610 if (data->log_separate) PetscCall(PetscLogEventEnd(PC_HPDDM_SetUp[0], data->levels[0]->ksp, nullptr, nullptr, nullptr)); 2611 if (ismatis) PetscCall(MatISGetLocalMat(C, &N)); 2612 else N = data->aux; 2613 if (!ctx) P = sub[0]; 2614 else P = S; 2615 /* going through the grid hierarchy */ 2616 for (n = 1; n < data->N; ++n) { 2617 if (data->log_separate) PetscCall(PetscLogEventBegin(PC_HPDDM_SetUp[n], data->levels[n]->ksp, nullptr, nullptr, nullptr)); 2618 /* method composed in the loaded symbol since there, SLEPc is used as well */ 2619 PetscTryMethod(data->levels[0]->ksp, "PCHPDDMSetUp_Private_C", (Mat *, Mat *, PetscInt, PetscInt *const, PC_HPDDM_Level **const), (&P, &N, n, &data->N, data->levels)); 2620 if (data->log_separate) PetscCall(PetscLogEventEnd(PC_HPDDM_SetUp[n], data->levels[n]->ksp, nullptr, nullptr, nullptr)); 2621 } 2622 /* reset to NULL to avoid any faulty use */ 2623 PetscCall(PetscObjectComposeFunction((PetscObject)data->levels[0]->ksp, "PCHPDDMSetUp_Private_C", nullptr)); 2624 if (!ismatis) PetscCall(PetscObjectComposeFunction((PetscObject)pc->pmat, "PCHPDDMAlgebraicAuxiliaryMat_C", nullptr)); 2625 else PetscCall(PetscObjectDereference((PetscObject)C)); /* matching PetscObjectReference() above */ 2626 for (n = 0; n < data->N - 1; ++n) 2627 if (data->levels[n]->P) { 2628 /* HPDDM internal work buffers */ 2629 data->levels[n]->P->setBuffer(); 2630 data->levels[n]->P->super::start(); 2631 } 2632 if (ismatis || !subdomains) PetscCall(PCHPDDMDestroySubMatrices_Private(PetscBool3ToBool(data->Neumann), PetscBool(algebraic && !block && overlap == -1), sub)); 2633 if (ismatis) data->is = nullptr; 2634 for (n = 0; n < data->N - 1 + (reused > 0); ++n) { 2635 if (data->levels[n]->P) { 2636 PC spc; 2637 2638 /* force the PC to be PCSHELL to do the coarse grid corrections */ 2639 PetscCall(KSPSetSkipPCSetFromOptions(data->levels[n]->ksp, PETSC_TRUE)); 2640 PetscCall(KSPGetPC(data->levels[n]->ksp, &spc)); 2641 PetscCall(PCSetType(spc, PCSHELL)); 2642 PetscCall(PCShellSetContext(spc, data->levels[n])); 2643 PetscCall(PCShellSetSetUp(spc, PCSetUp_HPDDMShell)); 2644 PetscCall(PCShellSetApply(spc, PCApply_HPDDMShell)); 2645 PetscCall(PCShellSetMatApply(spc, PCMatApply_HPDDMShell)); 2646 if (ctx && n == 0) { 2647 Mat Amat, Pmat; 2648 PetscInt m, M; 2649 std::tuple<Mat, PetscSF, Vec[2]> *ctx; 2650 2651 PetscCall(KSPGetOperators(data->levels[n]->ksp, nullptr, &Pmat)); 2652 PetscCall(MatGetLocalSize(Pmat, &m, nullptr)); 2653 PetscCall(MatGetSize(Pmat, &M, nullptr)); 2654 PetscCall(PetscNew(&ctx)); 2655 std::get<0>(*ctx) = S; 2656 std::get<1>(*ctx) = data->levels[n]->scatter; 2657 PetscCall(MatCreateShell(PetscObjectComm((PetscObject)data->levels[n]->ksp), m, m, M, M, ctx, &Amat)); 2658 PetscCall(MatShellSetOperation(Amat, MATOP_MULT, (void (*)(void))MatMult_Schur<false>)); 2659 PetscCall(MatShellSetOperation(Amat, MATOP_MULT_TRANSPOSE, (void (*)(void))MatMult_Schur<true>)); 2660 PetscCall(MatShellSetOperation(Amat, MATOP_DESTROY, (void (*)(void))MatDestroy_Schur)); 2661 PetscCall(MatCreateVecs(S, std::get<2>(*ctx), std::get<2>(*ctx) + 1)); 2662 PetscCall(KSPSetOperators(data->levels[n]->ksp, Amat, Pmat)); 2663 PetscCall(PetscObjectDereference((PetscObject)Amat)); 2664 } 2665 PetscCall(PCShellSetDestroy(spc, PCDestroy_HPDDMShell)); 2666 if (!data->levels[n]->pc) PetscCall(PCCreate(PetscObjectComm((PetscObject)data->levels[n]->ksp), &data->levels[n]->pc)); 2667 if (n < reused) { 2668 PetscCall(PCSetReusePreconditioner(spc, PETSC_TRUE)); 2669 PetscCall(PCSetReusePreconditioner(data->levels[n]->pc, PETSC_TRUE)); 2670 } 2671 PetscCall(PCSetUp(spc)); 2672 } 2673 } 2674 if (ctx) PetscCall(MatDestroy(&S)); 2675 if (overlap == -1) PetscCall(PetscObjectComposeFunction((PetscObject)pc->pmat, "PCHPDDMAlgebraicAuxiliaryMat_Private_C", nullptr)); 2676 } else flg = reused ? PETSC_FALSE : PETSC_TRUE; 2677 if (!ismatis && subdomains) { 2678 if (flg) PetscCall(KSPGetPC(data->levels[0]->ksp, &inner)); 2679 else inner = data->levels[0]->pc; 2680 if (inner) { 2681 if (!inner->setupcalled) PetscCall(PCSetType(inner, PCASM)); 2682 PetscCall(PCSetFromOptions(inner)); 2683 PetscCall(PetscStrcmp(((PetscObject)inner)->type_name, PCASM, &flg)); 2684 if (flg) { 2685 if (!inner->setupcalled) { /* evaluates to PETSC_FALSE when -pc_hpddm_block_splitting */ 2686 IS sorted; /* PCASM will sort the input IS, duplicate it to return an unmodified (PCHPDDM) input IS */ 2687 2688 PetscCall(ISDuplicate(is[0], &sorted)); 2689 PetscCall(PCASMSetLocalSubdomains(inner, 1, &sorted, &loc)); 2690 PetscCall(PetscObjectDereference((PetscObject)sorted)); 2691 } 2692 if (!PetscBool3ToBool(data->Neumann) && data->N > 1) { /* subdomain matrices are already created for the eigenproblem, reuse them for the fine-level PC */ 2693 PetscCall(PCHPDDMPermute_Private(*is, nullptr, nullptr, sub[0], &P, nullptr)); 2694 PetscCall(PCHPDDMCommunicationAvoidingPCASM_Private(inner, P, algebraic)); 2695 PetscCall(PetscObjectDereference((PetscObject)P)); 2696 } 2697 } 2698 } 2699 if (data->N > 1) { 2700 if (overlap != 1) PetscCall(PCHPDDMDestroySubMatrices_Private(PetscBool3ToBool(data->Neumann), PetscBool(algebraic && !block && overlap == -1), sub)); 2701 if (overlap == 1) PetscCall(MatDestroy(subA)); 2702 } 2703 } 2704 PetscCall(ISDestroy(&loc)); 2705 } else data->N = 1 + reused; /* enforce this value to 1 + reused if there is no way to build another level */ 2706 if (requested != data->N + reused) { 2707 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, 2708 data->N, pcpre ? pcpre : "")); 2709 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)); 2710 /* cannot use PCDestroy_HPDDMShell() because PCSHELL not set for unassembled levels */ 2711 for (n = data->N - 1; n < requested - 1; ++n) { 2712 if (data->levels[n]->P) { 2713 PetscCall(HPDDM::Schwarz<PetscScalar>::destroy(data->levels[n], PETSC_TRUE)); 2714 PetscCall(VecDestroyVecs(1, &data->levels[n]->v[0])); 2715 PetscCall(VecDestroyVecs(2, &data->levels[n]->v[1])); 2716 PetscCall(MatDestroy(data->levels[n]->V)); 2717 PetscCall(MatDestroy(data->levels[n]->V + 1)); 2718 PetscCall(MatDestroy(data->levels[n]->V + 2)); 2719 PetscCall(VecDestroy(&data->levels[n]->D)); 2720 PetscCall(PetscSFDestroy(&data->levels[n]->scatter)); 2721 } 2722 } 2723 if (reused) { 2724 for (n = reused; n < PETSC_PCHPDDM_MAXLEVELS && data->levels[n]; ++n) { 2725 PetscCall(KSPDestroy(&data->levels[n]->ksp)); 2726 PetscCall(PCDestroy(&data->levels[n]->pc)); 2727 } 2728 } 2729 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, 2730 data->N, reused, data->N, pcpre ? pcpre : "", pcpre ? pcpre : "", data->N); 2731 } 2732 /* these solvers are created after PCSetFromOptions() is called */ 2733 if (pc->setfromoptionscalled) { 2734 for (n = 0; n < data->N; ++n) { 2735 if (data->levels[n]->ksp) PetscCall(KSPSetFromOptions(data->levels[n]->ksp)); 2736 if (data->levels[n]->pc) PetscCall(PCSetFromOptions(data->levels[n]->pc)); 2737 } 2738 pc->setfromoptionscalled = 0; 2739 } 2740 data->N += reused; 2741 if (data->share && swap) { 2742 /* swap back pointers so that variables follow the user-provided numbering */ 2743 std::swap(C, data->aux); 2744 std::swap(uis, data->is); 2745 PetscCall(MatDestroy(&C)); 2746 PetscCall(ISDestroy(&uis)); 2747 } 2748 if (algebraic) PetscCall(MatDestroy(&data->aux)); 2749 if (unsorted && unsorted != is[0]) { 2750 PetscCall(ISCopy(unsorted, data->is)); 2751 PetscCall(ISDestroy(&unsorted)); 2752 } 2753 #if PetscDefined(USE_DEBUG) 2754 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); 2755 if (data->is) { 2756 PetscCall(ISEqualUnsorted(data->is, dis, &flg)); 2757 PetscCall(ISDestroy(&dis)); 2758 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Input IS and output IS are not equal"); 2759 } 2760 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); 2761 if (data->aux) { 2762 PetscCall(MatMultEqual(data->aux, daux, 20, &flg)); 2763 PetscCall(MatDestroy(&daux)); 2764 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Input Mat and output Mat are not equal"); 2765 } 2766 #endif 2767 PetscFunctionReturn(PETSC_SUCCESS); 2768 } 2769 2770 /*@ 2771 PCHPDDMSetCoarseCorrectionType - Sets the coarse correction type. 2772 2773 Collective 2774 2775 Input Parameters: 2776 + pc - preconditioner context 2777 - type - `PC_HPDDM_COARSE_CORRECTION_DEFLATED`, `PC_HPDDM_COARSE_CORRECTION_ADDITIVE`, `PC_HPDDM_COARSE_CORRECTION_BALANCED`, or `PC_HPDDM_COARSE_CORRECTION_NONE` 2778 2779 Options Database Key: 2780 . -pc_hpddm_coarse_correction <deflated, additive, balanced, none> - type of coarse correction to apply 2781 2782 Level: intermediate 2783 2784 .seealso: [](ch_ksp), `PCHPDDMGetCoarseCorrectionType()`, `PCHPDDM`, `PCHPDDMCoarseCorrectionType` 2785 @*/ 2786 PetscErrorCode PCHPDDMSetCoarseCorrectionType(PC pc, PCHPDDMCoarseCorrectionType type) 2787 { 2788 PetscFunctionBegin; 2789 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2790 PetscValidLogicalCollectiveEnum(pc, type, 2); 2791 PetscTryMethod(pc, "PCHPDDMSetCoarseCorrectionType_C", (PC, PCHPDDMCoarseCorrectionType), (pc, type)); 2792 PetscFunctionReturn(PETSC_SUCCESS); 2793 } 2794 2795 /*@ 2796 PCHPDDMGetCoarseCorrectionType - Gets the coarse correction type. 2797 2798 Input Parameter: 2799 . pc - preconditioner context 2800 2801 Output Parameter: 2802 . type - `PC_HPDDM_COARSE_CORRECTION_DEFLATED`, `PC_HPDDM_COARSE_CORRECTION_ADDITIVE`, `PC_HPDDM_COARSE_CORRECTION_BALANCED`, or `PC_HPDDM_COARSE_CORRECTION_NONE` 2803 2804 Level: intermediate 2805 2806 .seealso: [](ch_ksp), `PCHPDDMSetCoarseCorrectionType()`, `PCHPDDM`, `PCHPDDMCoarseCorrectionType` 2807 @*/ 2808 PetscErrorCode PCHPDDMGetCoarseCorrectionType(PC pc, PCHPDDMCoarseCorrectionType *type) 2809 { 2810 PetscFunctionBegin; 2811 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2812 if (type) { 2813 PetscAssertPointer(type, 2); 2814 PetscUseMethod(pc, "PCHPDDMGetCoarseCorrectionType_C", (PC, PCHPDDMCoarseCorrectionType *), (pc, type)); 2815 } 2816 PetscFunctionReturn(PETSC_SUCCESS); 2817 } 2818 2819 static PetscErrorCode PCHPDDMSetCoarseCorrectionType_HPDDM(PC pc, PCHPDDMCoarseCorrectionType type) 2820 { 2821 PC_HPDDM *data = (PC_HPDDM *)pc->data; 2822 2823 PetscFunctionBegin; 2824 data->correction = type; 2825 PetscFunctionReturn(PETSC_SUCCESS); 2826 } 2827 2828 static PetscErrorCode PCHPDDMGetCoarseCorrectionType_HPDDM(PC pc, PCHPDDMCoarseCorrectionType *type) 2829 { 2830 PC_HPDDM *data = (PC_HPDDM *)pc->data; 2831 2832 PetscFunctionBegin; 2833 *type = data->correction; 2834 PetscFunctionReturn(PETSC_SUCCESS); 2835 } 2836 2837 /*@ 2838 PCHPDDMSetSTShareSubKSP - Sets whether the `KSP` in SLEPc `ST` and the fine-level subdomain solver should be shared. 2839 2840 Input Parameters: 2841 + pc - preconditioner context 2842 - share - whether the `KSP` should be shared or not 2843 2844 Note: 2845 This is not the same as `PCSetReusePreconditioner()`. Given certain conditions (visible using -info), a symbolic factorization can be skipped 2846 when using a subdomain `PCType` such as `PCLU` or `PCCHOLESKY`. 2847 2848 Level: advanced 2849 2850 .seealso: [](ch_ksp), `PCHPDDM`, `PCHPDDMGetSTShareSubKSP()` 2851 @*/ 2852 PetscErrorCode PCHPDDMSetSTShareSubKSP(PC pc, PetscBool share) 2853 { 2854 PetscFunctionBegin; 2855 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2856 PetscTryMethod(pc, "PCHPDDMSetSTShareSubKSP_C", (PC, PetscBool), (pc, share)); 2857 PetscFunctionReturn(PETSC_SUCCESS); 2858 } 2859 2860 /*@ 2861 PCHPDDMGetSTShareSubKSP - Gets whether the `KSP` in SLEPc `ST` and the fine-level subdomain solver is shared. 2862 2863 Input Parameter: 2864 . pc - preconditioner context 2865 2866 Output Parameter: 2867 . share - whether the `KSP` is shared or not 2868 2869 Note: 2870 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 2871 when using a subdomain `PCType` such as `PCLU` or `PCCHOLESKY`. 2872 2873 Level: advanced 2874 2875 .seealso: [](ch_ksp), `PCHPDDM`, `PCHPDDMSetSTShareSubKSP()` 2876 @*/ 2877 PetscErrorCode PCHPDDMGetSTShareSubKSP(PC pc, PetscBool *share) 2878 { 2879 PetscFunctionBegin; 2880 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2881 if (share) { 2882 PetscAssertPointer(share, 2); 2883 PetscUseMethod(pc, "PCHPDDMGetSTShareSubKSP_C", (PC, PetscBool *), (pc, share)); 2884 } 2885 PetscFunctionReturn(PETSC_SUCCESS); 2886 } 2887 2888 static PetscErrorCode PCHPDDMSetSTShareSubKSP_HPDDM(PC pc, PetscBool share) 2889 { 2890 PC_HPDDM *data = (PC_HPDDM *)pc->data; 2891 2892 PetscFunctionBegin; 2893 data->share = share; 2894 PetscFunctionReturn(PETSC_SUCCESS); 2895 } 2896 2897 static PetscErrorCode PCHPDDMGetSTShareSubKSP_HPDDM(PC pc, PetscBool *share) 2898 { 2899 PC_HPDDM *data = (PC_HPDDM *)pc->data; 2900 2901 PetscFunctionBegin; 2902 *share = data->share; 2903 PetscFunctionReturn(PETSC_SUCCESS); 2904 } 2905 2906 /*@ 2907 PCHPDDMSetDeflationMat - Sets the deflation space used to assemble a coarser operator. 2908 2909 Input Parameters: 2910 + pc - preconditioner context 2911 . is - index set of the local deflation matrix 2912 - U - deflation sequential matrix stored as a `MATSEQDENSE` 2913 2914 Level: advanced 2915 2916 .seealso: [](ch_ksp), `PCHPDDM`, `PCDeflationSetSpace()`, `PCMGSetRestriction()` 2917 @*/ 2918 PetscErrorCode PCHPDDMSetDeflationMat(PC pc, IS is, Mat U) 2919 { 2920 PetscFunctionBegin; 2921 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2922 PetscValidHeaderSpecific(is, IS_CLASSID, 2); 2923 PetscValidHeaderSpecific(U, MAT_CLASSID, 3); 2924 PetscTryMethod(pc, "PCHPDDMSetDeflationMat_C", (PC, IS, Mat), (pc, is, U)); 2925 PetscFunctionReturn(PETSC_SUCCESS); 2926 } 2927 2928 static PetscErrorCode PCHPDDMSetDeflationMat_HPDDM(PC pc, IS is, Mat U) 2929 { 2930 PetscFunctionBegin; 2931 PetscCall(PCHPDDMSetAuxiliaryMat_Private(pc, is, U, PETSC_TRUE)); 2932 PetscFunctionReturn(PETSC_SUCCESS); 2933 } 2934 2935 PetscErrorCode HPDDMLoadDL_Private(PetscBool *found) 2936 { 2937 PetscBool flg; 2938 char lib[PETSC_MAX_PATH_LEN], dlib[PETSC_MAX_PATH_LEN], dir[PETSC_MAX_PATH_LEN]; 2939 2940 PetscFunctionBegin; 2941 PetscAssertPointer(found, 1); 2942 PetscCall(PetscStrncpy(dir, "${PETSC_LIB_DIR}", sizeof(dir))); 2943 PetscCall(PetscOptionsGetString(nullptr, nullptr, "-hpddm_dir", dir, sizeof(dir), nullptr)); 2944 PetscCall(PetscSNPrintf(lib, sizeof(lib), "%s/libhpddm_petsc", dir)); 2945 PetscCall(PetscDLLibraryRetrieve(PETSC_COMM_SELF, lib, dlib, 1024, found)); 2946 #if defined(SLEPC_LIB_DIR) /* this variable is passed during SLEPc ./configure when PETSc has not been configured */ 2947 if (!*found) { /* with --download-hpddm since slepcconf.h is not yet built (and thus can't be included) */ 2948 PetscCall(PetscStrncpy(dir, HPDDM_STR(SLEPC_LIB_DIR), sizeof(dir))); 2949 PetscCall(PetscSNPrintf(lib, sizeof(lib), "%s/libhpddm_petsc", dir)); 2950 PetscCall(PetscDLLibraryRetrieve(PETSC_COMM_SELF, lib, dlib, 1024, found)); 2951 } 2952 #endif 2953 if (!*found) { /* probable options for this to evaluate to PETSC_TRUE: system inconsistency (libhpddm_petsc moved by user?) or PETSc configured without --download-slepc */ 2954 PetscCall(PetscOptionsGetenv(PETSC_COMM_SELF, "SLEPC_DIR", dir, sizeof(dir), &flg)); 2955 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 */ 2956 PetscCall(PetscSNPrintf(lib, sizeof(lib), "%s/lib/libslepc", dir)); 2957 PetscCall(PetscDLLibraryRetrieve(PETSC_COMM_SELF, lib, dlib, 1024, found)); 2958 PetscCheck(*found, PETSC_COMM_SELF, PETSC_ERR_PLIB, "%s not found but SLEPC_DIR=%s", lib, dir); 2959 PetscCall(PetscDLLibraryAppend(PETSC_COMM_SELF, &PetscDLLibrariesLoaded, dlib)); 2960 PetscCall(PetscSNPrintf(lib, sizeof(lib), "%s/lib/libhpddm_petsc", dir)); /* libhpddm_petsc is always in the same directory as libslepc */ 2961 PetscCall(PetscDLLibraryRetrieve(PETSC_COMM_SELF, lib, dlib, 1024, found)); 2962 } 2963 } 2964 PetscCheck(*found, PETSC_COMM_SELF, PETSC_ERR_PLIB, "%s not found", lib); 2965 PetscCall(PetscDLLibraryAppend(PETSC_COMM_SELF, &PetscDLLibrariesLoaded, dlib)); 2966 PetscFunctionReturn(PETSC_SUCCESS); 2967 } 2968 2969 /*MC 2970 PCHPDDM - Interface with the HPDDM library. 2971 2972 This `PC` may be used to build multilevel spectral domain decomposition methods based on the GenEO framework {cite}`spillane2011robust` {cite}`al2021multilevel`. 2973 It may be viewed as an alternative to spectral 2974 AMGe or `PCBDDC` with adaptive selection of constraints. The interface is explained in details in {cite}`jolivetromanzampini2020` 2975 2976 The matrix used for building the preconditioner (Pmat) may be unassembled (`MATIS`), assembled (`MATAIJ`, `MATBAIJ`, or `MATSBAIJ`), hierarchical (`MATHTOOL`), `MATNORMAL`, `MATNORMALHERMITIAN`, or `MATSCHURCOMPLEMENT` (when `PCHPDDM` is used as part of an outer `PCFIELDSPLIT`). 2977 2978 For multilevel preconditioning, when using an assembled or hierarchical Pmat, one must provide an auxiliary local `Mat` (unassembled local operator for GenEO) using 2979 `PCHPDDMSetAuxiliaryMat()`. Calling this routine is not needed when using a `MATIS` Pmat, assembly is done internally using `MatConvert()`. 2980 2981 Options Database Keys: 2982 + -pc_hpddm_define_subdomains <true, default=false> - on the finest level, calls `PCASMSetLocalSubdomains()` with the `IS` supplied in `PCHPDDMSetAuxiliaryMat()` 2983 (not relevant with an unassembled Pmat) 2984 . -pc_hpddm_has_neumann <true, default=false> - on the finest level, informs the `PC` that the local Neumann matrix is supplied in `PCHPDDMSetAuxiliaryMat()` 2985 - -pc_hpddm_coarse_correction <type, default=deflated> - determines the `PCHPDDMCoarseCorrectionType` when calling `PCApply()` 2986 2987 Options for subdomain solvers, subdomain eigensolvers (for computing deflation vectors), and the coarse solver can be set using the following options database prefixes. 2988 .vb 2989 -pc_hpddm_levels_%d_pc_ 2990 -pc_hpddm_levels_%d_ksp_ 2991 -pc_hpddm_levels_%d_eps_ 2992 -pc_hpddm_levels_%d_p 2993 -pc_hpddm_levels_%d_mat_type 2994 -pc_hpddm_coarse_ 2995 -pc_hpddm_coarse_p 2996 -pc_hpddm_coarse_mat_type 2997 -pc_hpddm_coarse_mat_filter 2998 .ve 2999 3000 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 3001 -pc_hpddm_coarse_p 2 -pc_hpddm_coarse_mat_type baij will use 10 deflation vectors per subdomain on the fine "level 1", 3002 aggregate the fine subdomains into 4 "level 2" subdomains, then use 10 deflation vectors per subdomain on "level 2", 3003 and assemble the coarse matrix (of dimension 4 x 10 = 40) on two processes as a `MATBAIJ` (default is `MATSBAIJ`). 3004 3005 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. 3006 3007 Level: intermediate 3008 3009 Notes: 3010 This preconditioner requires that PETSc is built with SLEPc (``--download-slepc``). 3011 3012 By default, the underlying concurrent eigenproblems 3013 are solved using SLEPc shift-and-invert spectral transformation. This is usually what gives the best performance for GenEO, cf. 3014 {cite}`spillane2011robust` {cite}`jolivet2013scalabledd`. As 3015 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 3016 -pc_hpddm_levels_1_st_type sinvert. There are furthermore three options related to the (subdomain-wise local) eigensolver that are not described in 3017 SLEPc documentation since they are specific to `PCHPDDM`. 3018 .vb 3019 -pc_hpddm_levels_1_st_share_sub_ksp 3020 -pc_hpddm_levels_%d_eps_threshold 3021 -pc_hpddm_levels_1_eps_use_inertia 3022 .ve 3023 3024 The first option from the list only applies to the fine-level eigensolver, see `PCHPDDMSetSTShareSubKSP()`. The second option from the list is 3025 used to filter eigenmodes retrieved after convergence of `EPSSolve()` at "level N" such that eigenvectors used to define a "level N+1" coarse 3026 correction are associated to eigenvalues whose magnitude are lower or equal than -pc_hpddm_levels_N_eps_threshold. When using an `EPS` which cannot 3027 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 3028 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 3029 to supply -pc_hpddm_levels_1_eps_nev. This last option also only applies to the fine-level (N = 1) eigensolver. 3030 3031 See also {cite}`dolean2015introduction`, {cite}`al2022robust`, {cite}`al2022robustpd`, and {cite}`nataf2022recent` 3032 3033 .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCHPDDMSetAuxiliaryMat()`, `MATIS`, `PCBDDC`, `PCDEFLATION`, `PCTELESCOPE`, `PCASM`, 3034 `PCHPDDMSetCoarseCorrectionType()`, `PCHPDDMHasNeumannMat()`, `PCHPDDMSetRHSMat()`, `PCHPDDMSetDeflationMat()`, `PCHPDDMSetSTShareSubKSP()`, 3035 `PCHPDDMGetSTShareSubKSP()`, `PCHPDDMGetCoarseCorrectionType()`, `PCHPDDMGetComplexities()` 3036 M*/ 3037 PETSC_EXTERN PetscErrorCode PCCreate_HPDDM(PC pc) 3038 { 3039 PC_HPDDM *data; 3040 PetscBool found; 3041 3042 PetscFunctionBegin; 3043 if (!loadedSym) { 3044 PetscCall(HPDDMLoadDL_Private(&found)); 3045 if (found) PetscCall(PetscDLLibrarySym(PETSC_COMM_SELF, &PetscDLLibrariesLoaded, nullptr, "PCHPDDM_Internal", (void **)&loadedSym)); 3046 } 3047 PetscCheck(loadedSym, PETSC_COMM_SELF, PETSC_ERR_PLIB, "PCHPDDM_Internal symbol not found in loaded libhpddm_petsc"); 3048 PetscCall(PetscNew(&data)); 3049 pc->data = data; 3050 data->Neumann = PETSC_BOOL3_UNKNOWN; 3051 pc->ops->reset = PCReset_HPDDM; 3052 pc->ops->destroy = PCDestroy_HPDDM; 3053 pc->ops->setfromoptions = PCSetFromOptions_HPDDM; 3054 pc->ops->setup = PCSetUp_HPDDM; 3055 pc->ops->apply = PCApply_HPDDM; 3056 pc->ops->matapply = PCMatApply_HPDDM; 3057 pc->ops->view = PCView_HPDDM; 3058 pc->ops->presolve = PCPreSolve_HPDDM; 3059 3060 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetAuxiliaryMat_C", PCHPDDMSetAuxiliaryMat_HPDDM)); 3061 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMHasNeumannMat_C", PCHPDDMHasNeumannMat_HPDDM)); 3062 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetRHSMat_C", PCHPDDMSetRHSMat_HPDDM)); 3063 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetCoarseCorrectionType_C", PCHPDDMSetCoarseCorrectionType_HPDDM)); 3064 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMGetCoarseCorrectionType_C", PCHPDDMGetCoarseCorrectionType_HPDDM)); 3065 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetSTShareSubKSP_C", PCHPDDMSetSTShareSubKSP_HPDDM)); 3066 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMGetSTShareSubKSP_C", PCHPDDMGetSTShareSubKSP_HPDDM)); 3067 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetDeflationMat_C", PCHPDDMSetDeflationMat_HPDDM)); 3068 PetscFunctionReturn(PETSC_SUCCESS); 3069 } 3070 3071 /*@C 3072 PCHPDDMInitializePackage - This function initializes everything in the `PCHPDDM` package. It is called from `PCInitializePackage()`. 3073 3074 Level: developer 3075 3076 .seealso: [](ch_ksp), `PetscInitialize()` 3077 @*/ 3078 PetscErrorCode PCHPDDMInitializePackage(void) 3079 { 3080 char ename[32]; 3081 3082 PetscFunctionBegin; 3083 if (PCHPDDMPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS); 3084 PCHPDDMPackageInitialized = PETSC_TRUE; 3085 PetscCall(PetscRegisterFinalize(PCHPDDMFinalizePackage)); 3086 /* general events registered once during package initialization */ 3087 /* some of these events are not triggered in libpetsc, */ 3088 /* but rather directly in libhpddm_petsc, */ 3089 /* which is in charge of performing the following operations */ 3090 3091 /* domain decomposition structure from Pmat sparsity pattern */ 3092 PetscCall(PetscLogEventRegister("PCHPDDMStrc", PC_CLASSID, &PC_HPDDM_Strc)); 3093 /* Galerkin product, redistribution, and setup (not triggered in libpetsc) */ 3094 PetscCall(PetscLogEventRegister("PCHPDDMPtAP", PC_CLASSID, &PC_HPDDM_PtAP)); 3095 /* Galerkin product with summation, redistribution, and setup (not triggered in libpetsc) */ 3096 PetscCall(PetscLogEventRegister("PCHPDDMPtBP", PC_CLASSID, &PC_HPDDM_PtBP)); 3097 /* next level construction using PtAP and PtBP (not triggered in libpetsc) */ 3098 PetscCall(PetscLogEventRegister("PCHPDDMNext", PC_CLASSID, &PC_HPDDM_Next)); 3099 static_assert(PETSC_PCHPDDM_MAXLEVELS <= 9, "PETSC_PCHPDDM_MAXLEVELS value is too high"); 3100 for (PetscInt i = 1; i < PETSC_PCHPDDM_MAXLEVELS; ++i) { 3101 PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCHPDDMSetUp L%1" PetscInt_FMT, i)); 3102 /* events during a PCSetUp() at level #i _except_ the assembly */ 3103 /* of the Galerkin operator of the coarser level #(i + 1) */ 3104 PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_HPDDM_SetUp[i - 1])); 3105 PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCHPDDMSolve L%1" PetscInt_FMT, i)); 3106 /* events during a PCApply() at level #i _except_ */ 3107 /* the KSPSolve() of the coarser level #(i + 1) */ 3108 PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_HPDDM_Solve[i - 1])); 3109 } 3110 PetscFunctionReturn(PETSC_SUCCESS); 3111 } 3112 3113 /*@C 3114 PCHPDDMFinalizePackage - This function frees everything from the `PCHPDDM` package. It is called from `PetscFinalize()`. 3115 3116 Level: developer 3117 3118 .seealso: [](ch_ksp), `PetscFinalize()` 3119 @*/ 3120 PetscErrorCode PCHPDDMFinalizePackage(void) 3121 { 3122 PetscFunctionBegin; 3123 PCHPDDMPackageInitialized = PETSC_FALSE; 3124 PetscFunctionReturn(PETSC_SUCCESS); 3125 } 3126 3127 static PetscErrorCode MatMult_Harmonic(Mat A, Vec x, Vec y) 3128 { 3129 Harmonic h; /* [ A_00 A_01 ], furthermore, A_00 = [ A_loc,loc A_loc,ovl ], thus, A_01 = [ ] */ 3130 /* [ A_10 A_11 A_12 ] [ A_ovl,loc A_ovl,ovl ] [ A_ovl,1 ] */ 3131 Vec sub; /* y = A x = R_loc R_0 [ A_00 A_01 ]^-1 R_loc = [ I_loc ] */ 3132 /* [ A_10 A_11 ] R_1^T A_12 x [ ] */ 3133 PetscFunctionBegin; 3134 PetscCall(MatShellGetContext(A, &h)); 3135 PetscCall(VecSet(h->v, 0.0)); 3136 PetscCall(VecGetSubVector(h->v, h->is[0], &sub)); 3137 PetscCall(MatMult(h->A[0], x, sub)); 3138 PetscCall(VecRestoreSubVector(h->v, h->is[0], &sub)); 3139 PetscCall(KSPSolve(h->ksp, h->v, h->v)); 3140 PetscCall(VecISCopy(h->v, h->is[1], SCATTER_REVERSE, y)); 3141 PetscFunctionReturn(PETSC_SUCCESS); 3142 } 3143 3144 static PetscErrorCode MatMultTranspose_Harmonic(Mat A, Vec y, Vec x) 3145 { 3146 Harmonic h; /* x = A^T y = [ A_00 A_01 ]^-T R_0^T R_loc^T y */ 3147 Vec sub; /* A_12^T R_1 [ A_10 A_11 ] */ 3148 3149 PetscFunctionBegin; 3150 PetscCall(MatShellGetContext(A, &h)); 3151 PetscCall(VecSet(h->v, 0.0)); 3152 PetscCall(VecISCopy(h->v, h->is[1], SCATTER_FORWARD, y)); 3153 PetscCall(KSPSolveTranspose(h->ksp, h->v, h->v)); 3154 PetscCall(VecGetSubVector(h->v, h->is[0], &sub)); 3155 PetscCall(MatMultTranspose(h->A[0], sub, x)); 3156 PetscCall(VecRestoreSubVector(h->v, h->is[0], &sub)); 3157 PetscFunctionReturn(PETSC_SUCCESS); 3158 } 3159 3160 static PetscErrorCode MatProduct_AB_Harmonic(Mat S, Mat X, Mat Y, void *) 3161 { 3162 Harmonic h; 3163 Mat A, B; 3164 Vec a, b; 3165 3166 PetscFunctionBegin; 3167 PetscCall(MatShellGetContext(S, &h)); 3168 PetscCall(MatMatMult(h->A[0], X, MAT_INITIAL_MATRIX, PETSC_CURRENT, &A)); 3169 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, h->ksp->pc->mat->rmap->n, A->cmap->n, nullptr, &B)); 3170 for (PetscInt i = 0; i < A->cmap->n; ++i) { 3171 PetscCall(MatDenseGetColumnVecRead(A, i, &a)); 3172 PetscCall(MatDenseGetColumnVecWrite(B, i, &b)); 3173 PetscCall(VecISCopy(b, h->is[0], SCATTER_FORWARD, a)); 3174 PetscCall(MatDenseRestoreColumnVecWrite(B, i, &b)); 3175 PetscCall(MatDenseRestoreColumnVecRead(A, i, &a)); 3176 } 3177 PetscCall(MatDestroy(&A)); 3178 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, h->ksp->pc->mat->rmap->n, B->cmap->n, nullptr, &A)); 3179 PetscCall(KSPMatSolve(h->ksp, B, A)); 3180 PetscCall(MatDestroy(&B)); 3181 for (PetscInt i = 0; i < A->cmap->n; ++i) { 3182 PetscCall(MatDenseGetColumnVecRead(A, i, &a)); 3183 PetscCall(MatDenseGetColumnVecWrite(Y, i, &b)); 3184 PetscCall(VecISCopy(a, h->is[1], SCATTER_REVERSE, b)); 3185 PetscCall(MatDenseRestoreColumnVecWrite(Y, i, &b)); 3186 PetscCall(MatDenseRestoreColumnVecRead(A, i, &a)); 3187 } 3188 PetscCall(MatDestroy(&A)); 3189 PetscFunctionReturn(PETSC_SUCCESS); 3190 } 3191 3192 static PetscErrorCode MatDestroy_Harmonic(Mat A) 3193 { 3194 Harmonic h; 3195 3196 PetscFunctionBegin; 3197 PetscCall(MatShellGetContext(A, &h)); 3198 for (PetscInt i = 0; i < 5; ++i) PetscCall(ISDestroy(h->is + i)); 3199 PetscCall(PetscFree(h->is)); 3200 PetscCall(VecDestroy(&h->v)); 3201 for (PetscInt i = 0; i < 2; ++i) PetscCall(MatDestroy(h->A + i)); 3202 PetscCall(PetscFree(h->A)); 3203 PetscCall(KSPDestroy(&h->ksp)); 3204 PetscCall(PetscFree(h)); 3205 PetscFunctionReturn(PETSC_SUCCESS); 3206 } 3207