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 PetscCall(ISDuplicate(data->is, ov)); 2139 if (overlap > 1) PetscCall(MatIncreaseOverlap(P, 1, ov, overlap - 1)); 2140 PetscCall(ISDuplicate(ov[0], ov + 1)); 2141 PetscCall(MatIncreaseOverlap(P, 1, ov + 1, 1)); 2142 PetscCall(PetscNew(&h)); 2143 h->ksp = nullptr; 2144 PetscCall(PetscCalloc1(2, &h->A)); 2145 PetscCall(PetscOptionsHasName(nullptr, prefix, "-svd_nsv", &flg)); 2146 if (!flg) PetscCall(PetscOptionsHasName(nullptr, prefix, "-svd_relative_threshold", &flg)); 2147 PetscCall(ISSort(ov[0])); 2148 if (!flg) PetscCall(ISSort(ov[1])); 2149 PetscCall(PetscCalloc1(5, &h->is)); 2150 PetscCall(MatCreateSubMatrices(uaux ? uaux : P, 1, ov + !flg, ov + 1, MAT_INITIAL_MATRIX, &a)); /* submatrix from above, either square (!flg) or rectangular (flg) */ 2151 for (PetscInt j = 0; j < 2; ++j) { 2152 PetscCall(ISGetIndices(ov[j], i + j)); 2153 PetscCall(ISGetLocalSize(ov[j], n + j)); 2154 } 2155 v[1].reserve((n[1] - n[0]) / bs); 2156 for (PetscInt j = 0; j < n[1]; j += bs) { /* indices of the (2,2) block */ 2157 PetscInt location; 2158 PetscCall(ISLocate(ov[0], i[1][j], &location)); 2159 if (location < 0) v[1].emplace_back(j / bs); 2160 } 2161 if (!flg) { 2162 h->A[1] = a[0]; 2163 PetscCall(PetscObjectReference((PetscObject)h->A[1])); 2164 v[0].reserve((n[0] - P->rmap->n) / bs); 2165 for (PetscInt j = 0; j < n[1]; j += bs) { /* row indices of the (1,2) block */ 2166 PetscInt location; 2167 PetscCall(ISLocate(loc, i[1][j], &location)); 2168 if (location < 0) { 2169 PetscCall(ISLocate(ov[0], i[1][j], &location)); 2170 if (location >= 0) v[0].emplace_back(j / bs); 2171 } 2172 } 2173 PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, v[0].size(), v[0].data(), PETSC_USE_POINTER, &rows)); 2174 PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, v[1].size(), v[1].data(), PETSC_COPY_VALUES, h->is + 4)); 2175 PetscCall(MatCreateSubMatrix(a[0], rows, h->is[4], MAT_INITIAL_MATRIX, h->A)); /* A_12 submatrix from above */ 2176 PetscCall(ISDestroy(&rows)); 2177 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 */ 2178 PetscCall(ISEmbed(ov[0], ov[1], PETSC_TRUE, &rows)); 2179 PetscCall(MatCreateSubMatrix(a[0], rows, cols = rows, MAT_INITIAL_MATRIX, &A0)); /* [ A_00 A_01 ; A_10 A_11 ] submatrix from above */ 2180 PetscCall(ISDestroy(&rows)); 2181 v[0].clear(); 2182 PetscCall(ISEmbed(loc, ov[1], PETSC_TRUE, h->is + 3)); 2183 PetscCall(ISEmbed(data->is, ov[1], PETSC_TRUE, h->is + 2)); 2184 } 2185 v[0].reserve((n[0] - P->rmap->n) / bs); 2186 for (PetscInt j = 0; j < n[0]; j += bs) { 2187 PetscInt location; 2188 PetscCall(ISLocate(loc, i[0][j], &location)); 2189 if (location < 0) v[0].emplace_back(j / bs); 2190 } 2191 PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, v[0].size(), v[0].data(), PETSC_USE_POINTER, &rows)); 2192 for (PetscInt j = 0; j < 2; ++j) PetscCall(ISRestoreIndices(ov[j], i + j)); 2193 if (flg) { 2194 IS is; 2195 PetscCall(ISCreateStride(PETSC_COMM_SELF, a[0]->rmap->n, 0, 1, &is)); 2196 PetscCall(ISEmbed(ov[0], ov[1], PETSC_TRUE, &cols)); 2197 PetscCall(MatCreateSubMatrix(a[0], is, cols, MAT_INITIAL_MATRIX, &A0)); /* [ A_00 A_01 ; A_10 A_11 ] submatrix from above */ 2198 PetscCall(ISDestroy(&cols)); 2199 PetscCall(ISDestroy(&is)); 2200 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 */ 2201 PetscCall(ISEmbed(loc, data->is, PETSC_TRUE, h->is + 2)); 2202 PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, v[1].size(), v[1].data(), PETSC_USE_POINTER, &cols)); 2203 PetscCall(MatCreateSubMatrix(a[0], rows, cols, MAT_INITIAL_MATRIX, h->A)); /* A_12 submatrix from above */ 2204 PetscCall(ISDestroy(&cols)); 2205 } 2206 PetscCall(ISCreateStride(PETSC_COMM_SELF, A0->rmap->n, 0, 1, &stride)); 2207 PetscCall(ISEmbed(rows, stride, PETSC_TRUE, h->is)); 2208 PetscCall(ISDestroy(&stride)); 2209 PetscCall(ISDestroy(&rows)); 2210 PetscCall(ISEmbed(loc, ov[0], PETSC_TRUE, h->is + 1)); 2211 if (subdomains) { 2212 if (!data->levels[0]->pc) { 2213 PetscCall(PCCreate(PetscObjectComm((PetscObject)pc), &data->levels[0]->pc)); 2214 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_levels_1_", pcpre ? pcpre : "")); 2215 PetscCall(PCSetOptionsPrefix(data->levels[0]->pc, prefix)); 2216 PetscCall(PCSetOperators(data->levels[0]->pc, A, P)); 2217 } 2218 PetscCall(PCSetType(data->levels[0]->pc, PCASM)); 2219 if (!data->levels[0]->pc->setupcalled) PetscCall(PCASMSetLocalSubdomains(data->levels[0]->pc, 1, ov + !flg, &loc)); 2220 PetscCall(PCHPDDMCommunicationAvoidingPCASM_Private(data->levels[0]->pc, flg ? A0 : a[0], PETSC_TRUE)); 2221 if (!flg) ++overlap; 2222 if (data->share) { 2223 PetscInt n = -1; 2224 PetscTryMethod(data->levels[0]->pc, "PCASMGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (data->levels[0]->pc, &n, nullptr, &ksp)); 2225 PetscCheck(n == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of subdomain solver %" PetscInt_FMT " != 1", n); 2226 if (flg) { 2227 h->ksp = ksp[0]; 2228 PetscCall(PetscObjectReference((PetscObject)h->ksp)); 2229 } 2230 } 2231 } 2232 if (!h->ksp) { 2233 PetscBool share = data->share; 2234 PetscCall(KSPCreate(PETSC_COMM_SELF, &h->ksp)); 2235 PetscCall(KSPSetType(h->ksp, KSPPREONLY)); 2236 PetscCall(KSPSetOperators(h->ksp, A0, A0)); 2237 do { 2238 if (!data->share) { 2239 share = PETSC_FALSE; 2240 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_levels_1_%s", pcpre ? pcpre : "", flg ? "svd_" : "eps_")); 2241 PetscCall(KSPSetOptionsPrefix(h->ksp, prefix)); 2242 PetscCall(KSPSetFromOptions(h->ksp)); 2243 } else { 2244 MatSolverType type; 2245 PetscCall(KSPGetPC(ksp[0], &pc)); 2246 PetscCall(PetscObjectTypeCompareAny((PetscObject)pc, &data->share, PCLU, PCCHOLESKY, "")); 2247 if (data->share) { 2248 PetscCall(PCFactorGetMatSolverType(pc, &type)); 2249 if (!type) { 2250 if (PetscDefined(HAVE_MUMPS)) PetscCall(PCFactorSetMatSolverType(pc, MATSOLVERMUMPS)); 2251 else if (PetscDefined(HAVE_MKL_PARDISO)) PetscCall(PCFactorSetMatSolverType(pc, MATSOLVERMKL_PARDISO)); 2252 else data->share = PETSC_FALSE; 2253 if (data->share) PetscCall(PCSetFromOptions(pc)); 2254 } else { 2255 PetscCall(PetscStrcmp(type, MATSOLVERMUMPS, &data->share)); 2256 if (!data->share) PetscCall(PetscStrcmp(type, MATSOLVERMKL_PARDISO, &data->share)); 2257 } 2258 if (data->share) { 2259 std::tuple<KSP, IS, Vec[2]> *p; 2260 PetscCall(PCFactorGetMatrix(pc, &A)); 2261 PetscCall(MatFactorSetSchurIS(A, h->is[4])); 2262 PetscCall(KSPSetUp(ksp[0])); 2263 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_levels_1_eps_shell_", pcpre ? pcpre : "")); 2264 PetscCall(KSPSetOptionsPrefix(h->ksp, prefix)); 2265 PetscCall(KSPSetFromOptions(h->ksp)); 2266 PetscCall(KSPGetPC(h->ksp, &pc)); 2267 PetscCall(PCSetType(pc, PCSHELL)); 2268 PetscCall(PetscNew(&p)); 2269 std::get<0>(*p) = ksp[0]; 2270 PetscCall(ISEmbed(ov[0], ov[1], PETSC_TRUE, &std::get<1>(*p))); 2271 PetscCall(MatCreateVecs(A, std::get<2>(*p), std::get<2>(*p) + 1)); 2272 PetscCall(PCShellSetContext(pc, p)); 2273 PetscCall(PCShellSetApply(pc, PCApply_Schur)); 2274 PetscCall(PCShellSetApplyTranspose(pc, PCApply_Schur<Vec, true>)); 2275 PetscCall(PCShellSetMatApply(pc, PCApply_Schur<Mat>)); 2276 PetscCall(PCShellSetDestroy(pc, PCDestroy_Schur)); 2277 } 2278 } 2279 if (!data->share) PetscCall(PetscInfo(pc, "Cannot share subdomain KSP between SLEPc and PETSc since neither MUMPS nor MKL PARDISO is used\n")); 2280 } 2281 } 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 */ 2282 } 2283 PetscCall(ISDestroy(ov)); 2284 PetscCall(ISDestroy(ov + 1)); 2285 if (overlap == 1 && subdomains && flg) { 2286 *subA = A0; 2287 sub = subA; 2288 if (uaux) PetscCall(MatDestroy(&uaux)); 2289 } else PetscCall(MatDestroy(&A0)); 2290 PetscCall(MatCreateShell(PETSC_COMM_SELF, P->rmap->n, n[1] - n[0], P->rmap->n, n[1] - n[0], h, &data->aux)); 2291 PetscCall(KSPSetErrorIfNotConverged(h->ksp, PETSC_TRUE)); /* bail out as early as possible to avoid (apparently) unrelated error messages */ 2292 PetscCall(MatCreateVecs(h->ksp->pc->pmat, &h->v, nullptr)); 2293 PetscCall(MatShellSetOperation(data->aux, MATOP_MULT, (void (*)(void))MatMult_Harmonic)); 2294 PetscCall(MatShellSetOperation(data->aux, MATOP_MULT_TRANSPOSE, (void (*)(void))MatMultTranspose_Harmonic)); 2295 PetscCall(MatShellSetMatProductOperation(data->aux, MATPRODUCT_AB, nullptr, MatProduct_AB_Harmonic, nullptr, MATDENSE, MATDENSE)); 2296 PetscCall(MatShellSetOperation(data->aux, MATOP_DESTROY, (void (*)(void))MatDestroy_Harmonic)); 2297 PetscCall(MatDestroySubMatrices(1, &a)); 2298 } 2299 if (overlap != 1 || !subdomains) { 2300 PetscCall(MatCreateSubMatrices(uaux ? uaux : P, 1, is, is, MAT_INITIAL_MATRIX, &sub)); 2301 if (ismatis) { 2302 PetscCall(MatISGetLocalMat(C, &N)); 2303 PetscCall(PetscObjectTypeCompare((PetscObject)N, MATSEQSBAIJ, &flg)); 2304 if (flg) PetscCall(MatConvert(sub[0], MATSEQSBAIJ, MAT_INPLACE_MATRIX, sub)); 2305 PetscCall(MatISRestoreLocalMat(C, &N)); 2306 } 2307 } 2308 if (uaux) { 2309 PetscCall(MatDestroy(&uaux)); 2310 PetscCall(MatConvert(sub[0], MATSEQSBAIJ, MAT_INPLACE_MATRIX, sub)); 2311 } 2312 } 2313 } 2314 } else if (!ctx) { 2315 PetscCall(MatCreateSubMatrices(uaux, 1, is, is, MAT_INITIAL_MATRIX, &sub)); 2316 PetscCall(MatDestroy(&uaux)); 2317 PetscCall(MatConvert(sub[0], MATSEQSBAIJ, MAT_INPLACE_MATRIX, sub)); 2318 } 2319 if (data->N > 1) { 2320 /* Vec holding the partition of unity */ 2321 if (!data->levels[0]->D) { 2322 PetscCall(ISGetLocalSize(data->is, &n)); 2323 PetscCall(VecCreateMPI(PETSC_COMM_SELF, n, PETSC_DETERMINE, &data->levels[0]->D)); 2324 } 2325 if (data->share && overlap == -1) { 2326 Mat D; 2327 IS perm = nullptr; 2328 PetscInt size = -1; 2329 2330 if (!data->levels[0]->pc) { 2331 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_levels_1_", pcpre ? pcpre : "")); 2332 PetscCall(PCCreate(PetscObjectComm((PetscObject)pc), &data->levels[0]->pc)); 2333 PetscCall(PCSetOptionsPrefix(data->levels[0]->pc, prefix)); 2334 PetscCall(PCSetOperators(data->levels[0]->pc, A, P)); 2335 } 2336 PetscCall(PCSetType(data->levels[0]->pc, PCASM)); 2337 if (!ctx) { 2338 if (!data->levels[0]->pc->setupcalled) { 2339 IS sorted; /* PCASM will sort the input IS, duplicate it to return an unmodified (PCHPDDM) input IS */ 2340 PetscCall(ISDuplicate(is[0], &sorted)); 2341 PetscCall(PCASMSetLocalSubdomains(data->levels[0]->pc, 1, &sorted, &loc)); 2342 PetscCall(PetscObjectDereference((PetscObject)sorted)); 2343 } 2344 PetscCall(PCSetFromOptions(data->levels[0]->pc)); 2345 if (block) { 2346 PetscCall(PCHPDDMPermute_Private(unsorted, data->is, &uis, sub[0], &C, &perm)); 2347 PetscCall(PCHPDDMCommunicationAvoidingPCASM_Private(data->levels[0]->pc, C, algebraic)); 2348 } else PetscCall(PCSetUp(data->levels[0]->pc)); 2349 PetscTryMethod(data->levels[0]->pc, "PCASMGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (data->levels[0]->pc, &size, nullptr, &ksp)); 2350 if (size != 1) { 2351 data->share = PETSC_FALSE; 2352 PetscCheck(size == -1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of subdomain solver %" PetscInt_FMT " != 1", size); 2353 PetscCall(PetscInfo(pc, "Cannot share subdomain KSP between SLEPc and PETSc since PCASMGetSubKSP() not found in fine-level PC\n")); 2354 PetscCall(ISDestroy(&unsorted)); 2355 unsorted = is[0]; 2356 } else { 2357 const char *matpre; 2358 PetscBool cmp[4]; 2359 2360 if (!block && !ctx) PetscCall(PCHPDDMPermute_Private(unsorted, data->is, &uis, PetscBool3ToBool(data->Neumann) ? sub[0] : data->aux, &C, &perm)); 2361 if (perm) { /* unsorted input IS */ 2362 if (!PetscBool3ToBool(data->Neumann) && !block) { 2363 PetscCall(MatPermute(sub[0], perm, perm, &D)); /* permute since PCASM will call ISSort() */ 2364 PetscCall(MatHeaderReplace(sub[0], &D)); 2365 } 2366 if (data->B) { /* see PCHPDDMSetRHSMat() */ 2367 PetscCall(MatPermute(data->B, perm, perm, &D)); 2368 PetscCall(MatHeaderReplace(data->B, &D)); 2369 } 2370 PetscCall(ISDestroy(&perm)); 2371 } 2372 PetscCall(KSPGetOperators(ksp[0], subA, subA + 1)); 2373 PetscCall(PetscObjectReference((PetscObject)subA[0])); 2374 PetscCall(MatDuplicate(subA[1], MAT_SHARE_NONZERO_PATTERN, &D)); 2375 PetscCall(MatGetOptionsPrefix(subA[1], &matpre)); 2376 PetscCall(MatSetOptionsPrefix(D, matpre)); 2377 PetscCall(PetscObjectTypeCompare((PetscObject)D, MATNORMAL, cmp)); 2378 PetscCall(PetscObjectTypeCompare((PetscObject)C, MATNORMAL, cmp + 1)); 2379 if (!cmp[0]) PetscCall(PetscObjectTypeCompare((PetscObject)D, MATNORMALHERMITIAN, cmp + 2)); 2380 else cmp[2] = PETSC_FALSE; 2381 if (!cmp[1]) PetscCall(PetscObjectTypeCompare((PetscObject)C, MATNORMALHERMITIAN, cmp + 3)); 2382 else cmp[3] = PETSC_FALSE; 2383 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); 2384 if (!cmp[0] && !cmp[2]) { 2385 if (!block) PetscCall(MatAXPY(D, 1.0, C, SUBSET_NONZERO_PATTERN)); 2386 else { 2387 PetscCall(MatMissingDiagonal(D, cmp, nullptr)); 2388 if (cmp[0]) structure = DIFFERENT_NONZERO_PATTERN; /* data->aux has no missing diagonal entry */ 2389 PetscCall(MatAXPY(D, 1.0, data->aux, structure)); 2390 } 2391 } else { 2392 Mat mat[2]; 2393 2394 if (cmp[0]) { 2395 PetscCall(MatNormalGetMat(D, mat)); 2396 PetscCall(MatNormalGetMat(C, mat + 1)); 2397 } else { 2398 PetscCall(MatNormalHermitianGetMat(D, mat)); 2399 PetscCall(MatNormalHermitianGetMat(C, mat + 1)); 2400 } 2401 PetscCall(MatAXPY(mat[0], 1.0, mat[1], SUBSET_NONZERO_PATTERN)); 2402 } 2403 PetscCall(MatPropagateSymmetryOptions(C, D)); 2404 PetscCall(MatDestroy(&C)); 2405 C = D; 2406 /* swap pointers so that variables stay consistent throughout PCSetUp() */ 2407 std::swap(C, data->aux); 2408 std::swap(uis, data->is); 2409 swap = PETSC_TRUE; 2410 } 2411 } 2412 } 2413 } 2414 if (ctx) { 2415 PC_HPDDM *data_00 = (PC_HPDDM *)std::get<0>(*ctx)[0]->data; 2416 PC s; 2417 Mat A00, P00, A01 = nullptr, A10, A11, N, b[4]; 2418 IS sorted, is[2], *is_00; 2419 MatSolverType type; 2420 std::pair<PC, Vec[2]> *p; 2421 2422 n = -1; 2423 PetscTryMethod(data_00->levels[0]->pc, "PCASMGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (data_00->levels[0]->pc, &n, nullptr, &ksp)); 2424 PetscCheck(n == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of subdomain solver %" PetscInt_FMT " != 1", n); 2425 PetscCall(KSPGetOperators(ksp[0], subA, subA + 1)); 2426 PetscCall(ISGetLocalSize(data_00->is, &n)); 2427 if (n != subA[0]->rmap->n || n != subA[0]->cmap->n) { 2428 PetscCall(PCASMGetLocalSubdomains(data_00->levels[0]->pc, &n, &is_00, nullptr)); 2429 PetscCall(ISGetLocalSize(*is_00, &n)); 2430 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); 2431 } else is_00 = &data_00->is; 2432 PetscCall(PCHPDDMPermute_Private(unsorted, data->is, &uis, data->aux, &C, nullptr)); /* permute since PCASM works with a sorted IS */ 2433 std::swap(C, data->aux); 2434 std::swap(uis, data->is); 2435 swap = PETSC_TRUE; 2436 PetscCall(MatSchurComplementGetSubMatrices(P, &A00, &P00, std::get<1>(*ctx), &A10, &A11)); 2437 std::get<1>(*ctx)[1] = A10; 2438 PetscCall(PetscObjectTypeCompare((PetscObject)A10, MATTRANSPOSEVIRTUAL, &flg)); 2439 if (flg) PetscCall(MatTransposeGetMat(A10, &A01)); 2440 else { 2441 PetscBool flg; 2442 2443 PetscCall(PetscObjectTypeCompare((PetscObject)A10, MATHERMITIANTRANSPOSEVIRTUAL, &flg)); 2444 if (flg) PetscCall(MatHermitianTransposeGetMat(A10, &A01)); 2445 } 2446 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 */ 2447 PetscCall(ISSort(sorted)); /* this is to avoid changing users inputs, but it requires a new call to ISSort() here */ 2448 if (!A01) { 2449 PetscCall(MatSetOption(A10, MAT_SUBMAT_SINGLEIS, PETSC_TRUE)); 2450 PetscCall(MatCreateSubMatrices(A10, 1, &data->is, &sorted, MAT_INITIAL_MATRIX, &sub)); 2451 b[2] = sub[0]; 2452 PetscCall(PetscObjectReference((PetscObject)sub[0])); 2453 PetscCall(MatDestroySubMatrices(1, &sub)); 2454 PetscCall(PetscObjectTypeCompare((PetscObject)std::get<1>(*ctx)[0], MATTRANSPOSEVIRTUAL, &flg)); 2455 A10 = nullptr; 2456 if (flg) PetscCall(MatTransposeGetMat(std::get<1>(*ctx)[0], &A10)); 2457 else { 2458 PetscBool flg; 2459 2460 PetscCall(PetscObjectTypeCompare((PetscObject)std::get<1>(*ctx)[0], MATHERMITIANTRANSPOSEVIRTUAL, &flg)); 2461 if (flg) PetscCall(MatHermitianTransposeGetMat(std::get<1>(*ctx)[0], &A10)); 2462 } 2463 if (!A10) PetscCall(MatCreateSubMatrices(std::get<1>(*ctx)[0], 1, &sorted, &data->is, MAT_INITIAL_MATRIX, &sub)); 2464 else { 2465 if (flg) PetscCall(MatCreateTranspose(b[2], b + 1)); 2466 else PetscCall(MatCreateHermitianTranspose(b[2], b + 1)); 2467 } 2468 } else { 2469 PetscCall(MatSetOption(A01, MAT_SUBMAT_SINGLEIS, PETSC_TRUE)); 2470 PetscCall(MatCreateSubMatrices(A01, 1, &sorted, &data->is, MAT_INITIAL_MATRIX, &sub)); 2471 if (flg) PetscCall(MatCreateTranspose(*sub, b + 2)); 2472 else PetscCall(MatCreateHermitianTranspose(*sub, b + 2)); 2473 } 2474 if (A01 || !A10) { 2475 b[1] = sub[0]; 2476 PetscCall(PetscObjectReference((PetscObject)sub[0])); 2477 } 2478 PetscCall(MatDestroySubMatrices(1, &sub)); 2479 PetscCall(ISDestroy(&sorted)); 2480 b[3] = data->aux; 2481 PetscCall(MatCreateSchurComplement(subA[0], subA[1], b[1], b[2], b[3], &S)); 2482 PetscCall(MatSchurComplementSetKSP(S, ksp[0])); 2483 if (data->N != 1) { 2484 PetscCall(PCASMSetType(data->levels[0]->pc, PC_ASM_NONE)); /* "Neumann--Neumann" preconditioning with overlap and a Boolean partition of unity */ 2485 PetscCall(PCASMSetLocalSubdomains(data->levels[0]->pc, 1, &data->is, &loc)); 2486 PetscCall(PCSetFromOptions(data->levels[0]->pc)); /* action of eq. (15) of https://hal.science/hal-02343808v6/document (with a sign flip) */ 2487 s = data->levels[0]->pc; 2488 } else { 2489 is[0] = data->is; 2490 PetscCall(PetscObjectReference((PetscObject)is[0])); 2491 PetscCall(PetscObjectReference((PetscObject)b[3])); 2492 PetscCall(PCSetType(pc, PCASM)); /* change the type of the current PC */ 2493 data = nullptr; /* destroyed in the previous PCSetType(), so reset to NULL to avoid any faulty use */ 2494 PetscCall(PCAppendOptionsPrefix(pc, "pc_hpddm_coarse_")); /* same prefix as when using PCHPDDM with a single level */ 2495 PetscCall(PCASMSetLocalSubdomains(pc, 1, is, &loc)); 2496 PetscCall(ISDestroy(is)); 2497 PetscCall(ISDestroy(&loc)); 2498 s = pc; 2499 } 2500 PetscCall(PCHPDDMCommunicationAvoidingPCASM_Private(s, S, PETSC_TRUE)); /* the subdomain Mat is already known and the input IS of PCASMSetLocalSubdomains() is already sorted */ 2501 PetscTryMethod(s, "PCASMGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (s, &n, nullptr, &ksp)); 2502 PetscCheck(n == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of subdomain solver %" PetscInt_FMT " != 1", n); 2503 PetscCall(KSPGetPC(ksp[0], &inner)); 2504 PetscCall(PCSetType(inner, PCSHELL)); /* compute the action of the inverse of the local Schur complement with a PCSHELL */ 2505 b[0] = subA[0]; 2506 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 */ 2507 if (!data) PetscCall(PetscObjectDereference((PetscObject)b[3])); 2508 PetscCall(PetscObjectDereference((PetscObject)b[1])); 2509 PetscCall(PetscObjectDereference((PetscObject)b[2])); 2510 PetscCall(PCCreate(PETSC_COMM_SELF, &s)); 2511 PetscCall(PCSetOptionsPrefix(s, ((PetscObject)inner)->prefix)); 2512 PetscCall(PCSetOptionsPrefix(inner, nullptr)); 2513 PetscCall(KSPSetSkipPCSetFromOptions(ksp[0], PETSC_TRUE)); 2514 PetscCall(PCSetType(s, PCLU)); 2515 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 */ 2516 PetscCall(PCSetFromOptions(s)); 2517 PetscCall(PCFactorGetMatSolverType(s, &type)); 2518 PetscCall(PetscStrcmp(type, MATSOLVERMUMPS, &flg)); 2519 PetscCall(MatGetLocalSize(A11, &n, nullptr)); 2520 if (flg || n == 0) { 2521 PetscCall(PCSetOperators(s, N, N)); 2522 if (n) { 2523 PetscCall(PCFactorGetMatrix(s, b)); 2524 PetscCall(MatSetOptionsPrefix(*b, ((PetscObject)s)->prefix)); 2525 n = -1; 2526 PetscCall(PetscOptionsGetInt(nullptr, ((PetscObject)s)->prefix, "-mat_mumps_icntl_26", &n, nullptr)); 2527 if (n == 1) { /* allocates a square MatDense of size is[1]->map->n, so one */ 2528 PetscCall(MatNestGetISs(N, is, nullptr)); /* needs to be able to deactivate this path when dealing */ 2529 PetscCall(MatFactorSetSchurIS(*b, is[1])); /* with a large constraint space in order to avoid OOM */ 2530 } 2531 } else PetscCall(PCSetType(s, PCNONE)); /* empty local Schur complement (e.g., centralized on another process) */ 2532 } else { 2533 PetscCall(MatConvert(N, MATAIJ, MAT_INITIAL_MATRIX, b)); 2534 PetscCall(PCSetOperators(s, N, *b)); 2535 PetscCall(PetscObjectDereference((PetscObject)*b)); 2536 PetscCall(PetscObjectTypeCompareAny((PetscObject)s, &flg, PCLU, PCCHOLESKY, PCILU, PCICC, PCQR, "")); 2537 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 */ 2538 } 2539 PetscCall(PetscNew(&p)); 2540 p->first = s; 2541 if (n != 0) PetscCall(MatCreateVecs(*b, p->second, p->second + 1)); 2542 else p->second[0] = p->second[1] = nullptr; 2543 PetscCall(PCShellSetContext(inner, p)); 2544 PetscCall(PCShellSetApply(inner, PCApply_Nest)); 2545 PetscCall(PCShellSetView(inner, PCView_Nest)); 2546 PetscCall(PCShellSetDestroy(inner, PCDestroy_Nest)); 2547 PetscCall(PetscObjectDereference((PetscObject)N)); 2548 if (!data) { 2549 PetscCall(MatDestroy(&S)); 2550 PetscCall(ISDestroy(&unsorted)); 2551 PetscCall(MatDestroy(&C)); 2552 PetscCall(ISDestroy(&uis)); 2553 PetscCall(PetscFree(ctx)); 2554 #if PetscDefined(USE_DEBUG) 2555 PetscCall(ISDestroy(&dis)); 2556 PetscCall(MatDestroy(&daux)); 2557 #endif 2558 PetscFunctionReturn(PETSC_SUCCESS); 2559 } 2560 } 2561 if (!data->levels[0]->scatter) { 2562 PetscCall(MatCreateVecs(P, &xin, nullptr)); 2563 if (ismatis) PetscCall(MatDestroy(&P)); 2564 PetscCall(VecScatterCreate(xin, data->is, data->levels[0]->D, nullptr, &data->levels[0]->scatter)); 2565 PetscCall(VecDestroy(&xin)); 2566 } 2567 if (data->levels[0]->P) { 2568 /* if the pattern is the same and PCSetUp() has previously succeeded, reuse HPDDM buffers and connectivity */ 2569 PetscCall(HPDDM::Schwarz<PetscScalar>::destroy(data->levels[0], pc->setupcalled < 1 || pc->flag == DIFFERENT_NONZERO_PATTERN ? PETSC_TRUE : PETSC_FALSE)); 2570 } 2571 if (!data->levels[0]->P) data->levels[0]->P = new HPDDM::Schwarz<PetscScalar>(); 2572 if (data->log_separate) PetscCall(PetscLogEventBegin(PC_HPDDM_SetUp[0], data->levels[0]->ksp, nullptr, nullptr, nullptr)); 2573 else PetscCall(PetscLogEventBegin(PC_HPDDM_Strc, data->levels[0]->ksp, nullptr, nullptr, nullptr)); 2574 /* HPDDM internal data structure */ 2575 PetscCall(data->levels[0]->P->structure(loc, data->is, !ctx ? sub[0] : nullptr, ismatis ? C : data->aux, data->levels)); 2576 if (!data->log_separate) PetscCall(PetscLogEventEnd(PC_HPDDM_Strc, data->levels[0]->ksp, nullptr, nullptr, nullptr)); 2577 /* matrix pencil of the generalized eigenvalue problem on the overlap (GenEO) */ 2578 if (!ctx) { 2579 if (data->deflation || overlap != -1) weighted = data->aux; 2580 else if (!data->B) { 2581 PetscBool cmp; 2582 2583 PetscCall(MatDuplicate(sub[0], MAT_COPY_VALUES, &weighted)); 2584 PetscCall(PetscObjectTypeCompareAny((PetscObject)weighted, &cmp, MATNORMAL, MATNORMALHERMITIAN, "")); 2585 if (cmp) flg = PETSC_FALSE; 2586 PetscCall(MatDiagonalScale(weighted, data->levels[0]->D, data->levels[0]->D)); 2587 /* neither MatDuplicate() nor MatDiagonalScale() handles the symmetry options, so propagate the options explicitly */ 2588 /* only useful for -mat_type baij -pc_hpddm_levels_1_st_pc_type cholesky (no problem with MATAIJ or MATSBAIJ) */ 2589 PetscCall(MatPropagateSymmetryOptions(sub[0], weighted)); 2590 } else weighted = data->B; 2591 } else weighted = nullptr; 2592 /* SLEPc is used inside the loaded symbol */ 2593 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)); 2594 if (!ctx && data->share && overlap == -1) { 2595 Mat st[2]; 2596 PetscCall(KSPGetOperators(ksp[0], st, st + 1)); 2597 PetscCall(MatCopy(subA[0], st[0], structure)); 2598 if (subA[1] != subA[0] || st[1] != st[0]) PetscCall(MatCopy(subA[1], st[1], SAME_NONZERO_PATTERN)); 2599 PetscCall(PetscObjectDereference((PetscObject)subA[0])); 2600 } 2601 if (data->log_separate) PetscCall(PetscLogEventEnd(PC_HPDDM_SetUp[0], data->levels[0]->ksp, nullptr, nullptr, nullptr)); 2602 if (ismatis) PetscCall(MatISGetLocalMat(C, &N)); 2603 else N = data->aux; 2604 if (!ctx) P = sub[0]; 2605 else P = S; 2606 /* going through the grid hierarchy */ 2607 for (n = 1; n < data->N; ++n) { 2608 if (data->log_separate) PetscCall(PetscLogEventBegin(PC_HPDDM_SetUp[n], data->levels[n]->ksp, nullptr, nullptr, nullptr)); 2609 /* method composed in the loaded symbol since there, SLEPc is used as well */ 2610 PetscTryMethod(data->levels[0]->ksp, "PCHPDDMSetUp_Private_C", (Mat *, Mat *, PetscInt, PetscInt *const, PC_HPDDM_Level **const), (&P, &N, n, &data->N, data->levels)); 2611 if (data->log_separate) PetscCall(PetscLogEventEnd(PC_HPDDM_SetUp[n], data->levels[n]->ksp, nullptr, nullptr, nullptr)); 2612 } 2613 /* reset to NULL to avoid any faulty use */ 2614 PetscCall(PetscObjectComposeFunction((PetscObject)data->levels[0]->ksp, "PCHPDDMSetUp_Private_C", nullptr)); 2615 if (!ismatis) PetscCall(PetscObjectComposeFunction((PetscObject)pc->pmat, "PCHPDDMAlgebraicAuxiliaryMat_C", nullptr)); 2616 else PetscCall(PetscObjectDereference((PetscObject)C)); /* matching PetscObjectReference() above */ 2617 for (n = 0; n < data->N - 1; ++n) 2618 if (data->levels[n]->P) { 2619 /* HPDDM internal work buffers */ 2620 data->levels[n]->P->setBuffer(); 2621 data->levels[n]->P->super::start(); 2622 } 2623 if (ismatis || !subdomains) PetscCall(PCHPDDMDestroySubMatrices_Private(PetscBool3ToBool(data->Neumann), PetscBool(algebraic && !block && overlap == -1), sub)); 2624 if (ismatis) data->is = nullptr; 2625 for (n = 0; n < data->N - 1 + (reused > 0); ++n) { 2626 if (data->levels[n]->P) { 2627 PC spc; 2628 2629 /* force the PC to be PCSHELL to do the coarse grid corrections */ 2630 PetscCall(KSPSetSkipPCSetFromOptions(data->levels[n]->ksp, PETSC_TRUE)); 2631 PetscCall(KSPGetPC(data->levels[n]->ksp, &spc)); 2632 PetscCall(PCSetType(spc, PCSHELL)); 2633 PetscCall(PCShellSetContext(spc, data->levels[n])); 2634 PetscCall(PCShellSetSetUp(spc, PCSetUp_HPDDMShell)); 2635 PetscCall(PCShellSetApply(spc, PCApply_HPDDMShell)); 2636 PetscCall(PCShellSetMatApply(spc, PCMatApply_HPDDMShell)); 2637 if (ctx && n == 0) { 2638 Mat Amat, Pmat; 2639 PetscInt m, M; 2640 std::tuple<Mat, PetscSF, Vec[2]> *ctx; 2641 2642 PetscCall(KSPGetOperators(data->levels[n]->ksp, nullptr, &Pmat)); 2643 PetscCall(MatGetLocalSize(Pmat, &m, nullptr)); 2644 PetscCall(MatGetSize(Pmat, &M, nullptr)); 2645 PetscCall(PetscNew(&ctx)); 2646 std::get<0>(*ctx) = S; 2647 std::get<1>(*ctx) = data->levels[n]->scatter; 2648 PetscCall(MatCreateShell(PetscObjectComm((PetscObject)data->levels[n]->ksp), m, m, M, M, ctx, &Amat)); 2649 PetscCall(MatShellSetOperation(Amat, MATOP_MULT, (void (*)(void))MatMult_Schur<false>)); 2650 PetscCall(MatShellSetOperation(Amat, MATOP_MULT_TRANSPOSE, (void (*)(void))MatMult_Schur<true>)); 2651 PetscCall(MatShellSetOperation(Amat, MATOP_DESTROY, (void (*)(void))MatDestroy_Schur)); 2652 PetscCall(MatCreateVecs(S, std::get<2>(*ctx), std::get<2>(*ctx) + 1)); 2653 PetscCall(KSPSetOperators(data->levels[n]->ksp, Amat, Pmat)); 2654 PetscCall(PetscObjectDereference((PetscObject)Amat)); 2655 } 2656 PetscCall(PCShellSetDestroy(spc, PCDestroy_HPDDMShell)); 2657 if (!data->levels[n]->pc) PetscCall(PCCreate(PetscObjectComm((PetscObject)data->levels[n]->ksp), &data->levels[n]->pc)); 2658 if (n < reused) { 2659 PetscCall(PCSetReusePreconditioner(spc, PETSC_TRUE)); 2660 PetscCall(PCSetReusePreconditioner(data->levels[n]->pc, PETSC_TRUE)); 2661 } 2662 PetscCall(PCSetUp(spc)); 2663 } 2664 } 2665 if (ctx) PetscCall(MatDestroy(&S)); 2666 if (overlap == -1) PetscCall(PetscObjectComposeFunction((PetscObject)pc->pmat, "PCHPDDMAlgebraicAuxiliaryMat_Private_C", nullptr)); 2667 } else flg = reused ? PETSC_FALSE : PETSC_TRUE; 2668 if (!ismatis && subdomains) { 2669 if (flg) PetscCall(KSPGetPC(data->levels[0]->ksp, &inner)); 2670 else inner = data->levels[0]->pc; 2671 if (inner) { 2672 if (!inner->setupcalled) PetscCall(PCSetType(inner, PCASM)); 2673 PetscCall(PCSetFromOptions(inner)); 2674 PetscCall(PetscStrcmp(((PetscObject)inner)->type_name, PCASM, &flg)); 2675 if (flg) { 2676 if (!inner->setupcalled) { /* evaluates to PETSC_FALSE when -pc_hpddm_block_splitting */ 2677 IS sorted; /* PCASM will sort the input IS, duplicate it to return an unmodified (PCHPDDM) input IS */ 2678 2679 PetscCall(ISDuplicate(is[0], &sorted)); 2680 PetscCall(PCASMSetLocalSubdomains(inner, 1, &sorted, &loc)); 2681 PetscCall(PetscObjectDereference((PetscObject)sorted)); 2682 } 2683 if (!PetscBool3ToBool(data->Neumann) && data->N > 1) { /* subdomain matrices are already created for the eigenproblem, reuse them for the fine-level PC */ 2684 PetscCall(PCHPDDMPermute_Private(*is, nullptr, nullptr, sub[0], &P, nullptr)); 2685 PetscCall(PCHPDDMCommunicationAvoidingPCASM_Private(inner, P, algebraic)); 2686 PetscCall(PetscObjectDereference((PetscObject)P)); 2687 } 2688 } 2689 } 2690 if (data->N > 1) { 2691 if (overlap != 1) PetscCall(PCHPDDMDestroySubMatrices_Private(PetscBool3ToBool(data->Neumann), PetscBool(algebraic && !block && overlap == -1), sub)); 2692 if (overlap == 1) PetscCall(MatDestroy(subA)); 2693 } 2694 } 2695 PetscCall(ISDestroy(&loc)); 2696 } else data->N = 1 + reused; /* enforce this value to 1 + reused if there is no way to build another level */ 2697 if (requested != data->N + reused) { 2698 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, 2699 data->N, pcpre ? pcpre : "")); 2700 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)); 2701 /* cannot use PCDestroy_HPDDMShell() because PCSHELL not set for unassembled levels */ 2702 for (n = data->N - 1; n < requested - 1; ++n) { 2703 if (data->levels[n]->P) { 2704 PetscCall(HPDDM::Schwarz<PetscScalar>::destroy(data->levels[n], PETSC_TRUE)); 2705 PetscCall(VecDestroyVecs(1, &data->levels[n]->v[0])); 2706 PetscCall(VecDestroyVecs(2, &data->levels[n]->v[1])); 2707 PetscCall(MatDestroy(data->levels[n]->V)); 2708 PetscCall(MatDestroy(data->levels[n]->V + 1)); 2709 PetscCall(MatDestroy(data->levels[n]->V + 2)); 2710 PetscCall(VecDestroy(&data->levels[n]->D)); 2711 PetscCall(PetscSFDestroy(&data->levels[n]->scatter)); 2712 } 2713 } 2714 if (reused) { 2715 for (n = reused; n < PETSC_PCHPDDM_MAXLEVELS && data->levels[n]; ++n) { 2716 PetscCall(KSPDestroy(&data->levels[n]->ksp)); 2717 PetscCall(PCDestroy(&data->levels[n]->pc)); 2718 } 2719 } 2720 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, 2721 data->N, reused, data->N, pcpre ? pcpre : "", pcpre ? pcpre : "", data->N); 2722 } 2723 /* these solvers are created after PCSetFromOptions() is called */ 2724 if (pc->setfromoptionscalled) { 2725 for (n = 0; n < data->N; ++n) { 2726 if (data->levels[n]->ksp) PetscCall(KSPSetFromOptions(data->levels[n]->ksp)); 2727 if (data->levels[n]->pc) PetscCall(PCSetFromOptions(data->levels[n]->pc)); 2728 } 2729 pc->setfromoptionscalled = 0; 2730 } 2731 data->N += reused; 2732 if (data->share && swap) { 2733 /* swap back pointers so that variables follow the user-provided numbering */ 2734 std::swap(C, data->aux); 2735 std::swap(uis, data->is); 2736 PetscCall(MatDestroy(&C)); 2737 PetscCall(ISDestroy(&uis)); 2738 } 2739 if (algebraic) PetscCall(MatDestroy(&data->aux)); 2740 if (unsorted && unsorted != is[0]) { 2741 PetscCall(ISCopy(unsorted, data->is)); 2742 PetscCall(ISDestroy(&unsorted)); 2743 } 2744 #if PetscDefined(USE_DEBUG) 2745 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); 2746 if (data->is) { 2747 PetscCall(ISEqualUnsorted(data->is, dis, &flg)); 2748 PetscCall(ISDestroy(&dis)); 2749 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Input IS and output IS are not equal"); 2750 } 2751 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); 2752 if (data->aux) { 2753 PetscCall(MatMultEqual(data->aux, daux, 20, &flg)); 2754 PetscCall(MatDestroy(&daux)); 2755 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Input Mat and output Mat are not equal"); 2756 } 2757 #endif 2758 PetscFunctionReturn(PETSC_SUCCESS); 2759 } 2760 2761 /*@ 2762 PCHPDDMSetCoarseCorrectionType - Sets the coarse correction type. 2763 2764 Collective 2765 2766 Input Parameters: 2767 + pc - preconditioner context 2768 - type - `PC_HPDDM_COARSE_CORRECTION_DEFLATED`, `PC_HPDDM_COARSE_CORRECTION_ADDITIVE`, `PC_HPDDM_COARSE_CORRECTION_BALANCED`, or `PC_HPDDM_COARSE_CORRECTION_NONE` 2769 2770 Options Database Key: 2771 . -pc_hpddm_coarse_correction <deflated, additive, balanced, none> - type of coarse correction to apply 2772 2773 Level: intermediate 2774 2775 .seealso: [](ch_ksp), `PCHPDDMGetCoarseCorrectionType()`, `PCHPDDM`, `PCHPDDMCoarseCorrectionType` 2776 @*/ 2777 PetscErrorCode PCHPDDMSetCoarseCorrectionType(PC pc, PCHPDDMCoarseCorrectionType type) 2778 { 2779 PetscFunctionBegin; 2780 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2781 PetscValidLogicalCollectiveEnum(pc, type, 2); 2782 PetscTryMethod(pc, "PCHPDDMSetCoarseCorrectionType_C", (PC, PCHPDDMCoarseCorrectionType), (pc, type)); 2783 PetscFunctionReturn(PETSC_SUCCESS); 2784 } 2785 2786 /*@ 2787 PCHPDDMGetCoarseCorrectionType - Gets the coarse correction type. 2788 2789 Input Parameter: 2790 . pc - preconditioner context 2791 2792 Output Parameter: 2793 . type - `PC_HPDDM_COARSE_CORRECTION_DEFLATED`, `PC_HPDDM_COARSE_CORRECTION_ADDITIVE`, `PC_HPDDM_COARSE_CORRECTION_BALANCED`, or `PC_HPDDM_COARSE_CORRECTION_NONE` 2794 2795 Level: intermediate 2796 2797 .seealso: [](ch_ksp), `PCHPDDMSetCoarseCorrectionType()`, `PCHPDDM`, `PCHPDDMCoarseCorrectionType` 2798 @*/ 2799 PetscErrorCode PCHPDDMGetCoarseCorrectionType(PC pc, PCHPDDMCoarseCorrectionType *type) 2800 { 2801 PetscFunctionBegin; 2802 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2803 if (type) { 2804 PetscAssertPointer(type, 2); 2805 PetscUseMethod(pc, "PCHPDDMGetCoarseCorrectionType_C", (PC, PCHPDDMCoarseCorrectionType *), (pc, type)); 2806 } 2807 PetscFunctionReturn(PETSC_SUCCESS); 2808 } 2809 2810 static PetscErrorCode PCHPDDMSetCoarseCorrectionType_HPDDM(PC pc, PCHPDDMCoarseCorrectionType type) 2811 { 2812 PC_HPDDM *data = (PC_HPDDM *)pc->data; 2813 2814 PetscFunctionBegin; 2815 data->correction = type; 2816 PetscFunctionReturn(PETSC_SUCCESS); 2817 } 2818 2819 static PetscErrorCode PCHPDDMGetCoarseCorrectionType_HPDDM(PC pc, PCHPDDMCoarseCorrectionType *type) 2820 { 2821 PC_HPDDM *data = (PC_HPDDM *)pc->data; 2822 2823 PetscFunctionBegin; 2824 *type = data->correction; 2825 PetscFunctionReturn(PETSC_SUCCESS); 2826 } 2827 2828 /*@ 2829 PCHPDDMSetSTShareSubKSP - Sets whether the `KSP` in SLEPc `ST` and the fine-level subdomain solver should be shared. 2830 2831 Input Parameters: 2832 + pc - preconditioner context 2833 - share - whether the `KSP` should be shared or not 2834 2835 Note: 2836 This is not the same as `PCSetReusePreconditioner()`. Given certain conditions (visible using -info), a symbolic factorization can be skipped 2837 when using a subdomain `PCType` such as `PCLU` or `PCCHOLESKY`. 2838 2839 Level: advanced 2840 2841 .seealso: [](ch_ksp), `PCHPDDM`, `PCHPDDMGetSTShareSubKSP()` 2842 @*/ 2843 PetscErrorCode PCHPDDMSetSTShareSubKSP(PC pc, PetscBool share) 2844 { 2845 PetscFunctionBegin; 2846 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2847 PetscTryMethod(pc, "PCHPDDMSetSTShareSubKSP_C", (PC, PetscBool), (pc, share)); 2848 PetscFunctionReturn(PETSC_SUCCESS); 2849 } 2850 2851 /*@ 2852 PCHPDDMGetSTShareSubKSP - Gets whether the `KSP` in SLEPc `ST` and the fine-level subdomain solver is shared. 2853 2854 Input Parameter: 2855 . pc - preconditioner context 2856 2857 Output Parameter: 2858 . share - whether the `KSP` is shared or not 2859 2860 Note: 2861 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 2862 when using a subdomain `PCType` such as `PCLU` or `PCCHOLESKY`. 2863 2864 Level: advanced 2865 2866 .seealso: [](ch_ksp), `PCHPDDM`, `PCHPDDMSetSTShareSubKSP()` 2867 @*/ 2868 PetscErrorCode PCHPDDMGetSTShareSubKSP(PC pc, PetscBool *share) 2869 { 2870 PetscFunctionBegin; 2871 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2872 if (share) { 2873 PetscAssertPointer(share, 2); 2874 PetscUseMethod(pc, "PCHPDDMGetSTShareSubKSP_C", (PC, PetscBool *), (pc, share)); 2875 } 2876 PetscFunctionReturn(PETSC_SUCCESS); 2877 } 2878 2879 static PetscErrorCode PCHPDDMSetSTShareSubKSP_HPDDM(PC pc, PetscBool share) 2880 { 2881 PC_HPDDM *data = (PC_HPDDM *)pc->data; 2882 2883 PetscFunctionBegin; 2884 data->share = share; 2885 PetscFunctionReturn(PETSC_SUCCESS); 2886 } 2887 2888 static PetscErrorCode PCHPDDMGetSTShareSubKSP_HPDDM(PC pc, PetscBool *share) 2889 { 2890 PC_HPDDM *data = (PC_HPDDM *)pc->data; 2891 2892 PetscFunctionBegin; 2893 *share = data->share; 2894 PetscFunctionReturn(PETSC_SUCCESS); 2895 } 2896 2897 /*@ 2898 PCHPDDMSetDeflationMat - Sets the deflation space used to assemble a coarser operator. 2899 2900 Input Parameters: 2901 + pc - preconditioner context 2902 . is - index set of the local deflation matrix 2903 - U - deflation sequential matrix stored as a `MATSEQDENSE` 2904 2905 Level: advanced 2906 2907 .seealso: [](ch_ksp), `PCHPDDM`, `PCDeflationSetSpace()`, `PCMGSetRestriction()` 2908 @*/ 2909 PetscErrorCode PCHPDDMSetDeflationMat(PC pc, IS is, Mat U) 2910 { 2911 PetscFunctionBegin; 2912 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2913 PetscValidHeaderSpecific(is, IS_CLASSID, 2); 2914 PetscValidHeaderSpecific(U, MAT_CLASSID, 3); 2915 PetscTryMethod(pc, "PCHPDDMSetDeflationMat_C", (PC, IS, Mat), (pc, is, U)); 2916 PetscFunctionReturn(PETSC_SUCCESS); 2917 } 2918 2919 static PetscErrorCode PCHPDDMSetDeflationMat_HPDDM(PC pc, IS is, Mat U) 2920 { 2921 PetscFunctionBegin; 2922 PetscCall(PCHPDDMSetAuxiliaryMat_Private(pc, is, U, PETSC_TRUE)); 2923 PetscFunctionReturn(PETSC_SUCCESS); 2924 } 2925 2926 PetscErrorCode HPDDMLoadDL_Private(PetscBool *found) 2927 { 2928 PetscBool flg; 2929 char lib[PETSC_MAX_PATH_LEN], dlib[PETSC_MAX_PATH_LEN], dir[PETSC_MAX_PATH_LEN]; 2930 2931 PetscFunctionBegin; 2932 PetscAssertPointer(found, 1); 2933 PetscCall(PetscStrncpy(dir, "${PETSC_LIB_DIR}", sizeof(dir))); 2934 PetscCall(PetscOptionsGetString(nullptr, nullptr, "-hpddm_dir", dir, sizeof(dir), nullptr)); 2935 PetscCall(PetscSNPrintf(lib, sizeof(lib), "%s/libhpddm_petsc", dir)); 2936 PetscCall(PetscDLLibraryRetrieve(PETSC_COMM_SELF, lib, dlib, 1024, found)); 2937 #if defined(SLEPC_LIB_DIR) /* this variable is passed during SLEPc ./configure when PETSc has not been configured */ 2938 if (!*found) { /* with --download-hpddm since slepcconf.h is not yet built (and thus can't be included) */ 2939 PetscCall(PetscStrncpy(dir, HPDDM_STR(SLEPC_LIB_DIR), sizeof(dir))); 2940 PetscCall(PetscSNPrintf(lib, sizeof(lib), "%s/libhpddm_petsc", dir)); 2941 PetscCall(PetscDLLibraryRetrieve(PETSC_COMM_SELF, lib, dlib, 1024, found)); 2942 } 2943 #endif 2944 if (!*found) { /* probable options for this to evaluate to PETSC_TRUE: system inconsistency (libhpddm_petsc moved by user?) or PETSc configured without --download-slepc */ 2945 PetscCall(PetscOptionsGetenv(PETSC_COMM_SELF, "SLEPC_DIR", dir, sizeof(dir), &flg)); 2946 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 */ 2947 PetscCall(PetscSNPrintf(lib, sizeof(lib), "%s/lib/libslepc", dir)); 2948 PetscCall(PetscDLLibraryRetrieve(PETSC_COMM_SELF, lib, dlib, 1024, found)); 2949 PetscCheck(*found, PETSC_COMM_SELF, PETSC_ERR_PLIB, "%s not found but SLEPC_DIR=%s", lib, dir); 2950 PetscCall(PetscDLLibraryAppend(PETSC_COMM_SELF, &PetscDLLibrariesLoaded, dlib)); 2951 PetscCall(PetscSNPrintf(lib, sizeof(lib), "%s/lib/libhpddm_petsc", dir)); /* libhpddm_petsc is always in the same directory as libslepc */ 2952 PetscCall(PetscDLLibraryRetrieve(PETSC_COMM_SELF, lib, dlib, 1024, found)); 2953 } 2954 } 2955 PetscCheck(*found, PETSC_COMM_SELF, PETSC_ERR_PLIB, "%s not found", lib); 2956 PetscCall(PetscDLLibraryAppend(PETSC_COMM_SELF, &PetscDLLibrariesLoaded, dlib)); 2957 PetscFunctionReturn(PETSC_SUCCESS); 2958 } 2959 2960 /*MC 2961 PCHPDDM - Interface with the HPDDM library. 2962 2963 This `PC` may be used to build multilevel spectral domain decomposition methods based on the GenEO framework {cite}`spillane2011robust` {cite}`al2021multilevel`. 2964 It may be viewed as an alternative to spectral 2965 AMGe or `PCBDDC` with adaptive selection of constraints. The interface is explained in details in {cite}`jolivetromanzampini2020` 2966 2967 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`). 2968 2969 For multilevel preconditioning, when using an assembled or hierarchical Pmat, one must provide an auxiliary local `Mat` (unassembled local operator for GenEO) using 2970 `PCHPDDMSetAuxiliaryMat()`. Calling this routine is not needed when using a `MATIS` Pmat, assembly is done internally using `MatConvert()`. 2971 2972 Options Database Keys: 2973 + -pc_hpddm_define_subdomains <true, default=false> - on the finest level, calls `PCASMSetLocalSubdomains()` with the `IS` supplied in `PCHPDDMSetAuxiliaryMat()` 2974 (not relevant with an unassembled Pmat) 2975 . -pc_hpddm_has_neumann <true, default=false> - on the finest level, informs the `PC` that the local Neumann matrix is supplied in `PCHPDDMSetAuxiliaryMat()` 2976 - -pc_hpddm_coarse_correction <type, default=deflated> - determines the `PCHPDDMCoarseCorrectionType` when calling `PCApply()` 2977 2978 Options for subdomain solvers, subdomain eigensolvers (for computing deflation vectors), and the coarse solver can be set using the following options database prefixes. 2979 .vb 2980 -pc_hpddm_levels_%d_pc_ 2981 -pc_hpddm_levels_%d_ksp_ 2982 -pc_hpddm_levels_%d_eps_ 2983 -pc_hpddm_levels_%d_p 2984 -pc_hpddm_levels_%d_mat_type 2985 -pc_hpddm_coarse_ 2986 -pc_hpddm_coarse_p 2987 -pc_hpddm_coarse_mat_type 2988 -pc_hpddm_coarse_mat_filter 2989 .ve 2990 2991 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 2992 -pc_hpddm_coarse_p 2 -pc_hpddm_coarse_mat_type baij will use 10 deflation vectors per subdomain on the fine "level 1", 2993 aggregate the fine subdomains into 4 "level 2" subdomains, then use 10 deflation vectors per subdomain on "level 2", 2994 and assemble the coarse matrix (of dimension 4 x 10 = 40) on two processes as a `MATBAIJ` (default is `MATSBAIJ`). 2995 2996 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. 2997 2998 Level: intermediate 2999 3000 Notes: 3001 This preconditioner requires that PETSc is built with SLEPc (``--download-slepc``). 3002 3003 By default, the underlying concurrent eigenproblems 3004 are solved using SLEPc shift-and-invert spectral transformation. This is usually what gives the best performance for GenEO, cf. 3005 {cite}`spillane2011robust` {cite}`jolivet2013scalabledd`. As 3006 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 3007 -pc_hpddm_levels_1_st_type sinvert. There are furthermore three options related to the (subdomain-wise local) eigensolver that are not described in 3008 SLEPc documentation since they are specific to `PCHPDDM`. 3009 .vb 3010 -pc_hpddm_levels_1_st_share_sub_ksp 3011 -pc_hpddm_levels_%d_eps_threshold 3012 -pc_hpddm_levels_1_eps_use_inertia 3013 .ve 3014 3015 The first option from the list only applies to the fine-level eigensolver, see `PCHPDDMSetSTShareSubKSP()`. The second option from the list is 3016 used to filter eigenmodes retrieved after convergence of `EPSSolve()` at "level N" such that eigenvectors used to define a "level N+1" coarse 3017 correction are associated to eigenvalues whose magnitude are lower or equal than -pc_hpddm_levels_N_eps_threshold. When using an `EPS` which cannot 3018 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 3019 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 3020 to supply -pc_hpddm_levels_1_eps_nev. This last option also only applies to the fine-level (N = 1) eigensolver. 3021 3022 See also {cite}`dolean2015introduction`, {cite}`al2022robust`, {cite}`al2022robustpd`, and {cite}`nataf2022recent` 3023 3024 .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCHPDDMSetAuxiliaryMat()`, `MATIS`, `PCBDDC`, `PCDEFLATION`, `PCTELESCOPE`, `PCASM`, 3025 `PCHPDDMSetCoarseCorrectionType()`, `PCHPDDMHasNeumannMat()`, `PCHPDDMSetRHSMat()`, `PCHPDDMSetDeflationMat()`, `PCHPDDMSetSTShareSubKSP()`, 3026 `PCHPDDMGetSTShareSubKSP()`, `PCHPDDMGetCoarseCorrectionType()`, `PCHPDDMGetComplexities()` 3027 M*/ 3028 PETSC_EXTERN PetscErrorCode PCCreate_HPDDM(PC pc) 3029 { 3030 PC_HPDDM *data; 3031 PetscBool found; 3032 3033 PetscFunctionBegin; 3034 if (!loadedSym) { 3035 PetscCall(HPDDMLoadDL_Private(&found)); 3036 if (found) PetscCall(PetscDLLibrarySym(PETSC_COMM_SELF, &PetscDLLibrariesLoaded, nullptr, "PCHPDDM_Internal", (void **)&loadedSym)); 3037 } 3038 PetscCheck(loadedSym, PETSC_COMM_SELF, PETSC_ERR_PLIB, "PCHPDDM_Internal symbol not found in loaded libhpddm_petsc"); 3039 PetscCall(PetscNew(&data)); 3040 pc->data = data; 3041 data->Neumann = PETSC_BOOL3_UNKNOWN; 3042 pc->ops->reset = PCReset_HPDDM; 3043 pc->ops->destroy = PCDestroy_HPDDM; 3044 pc->ops->setfromoptions = PCSetFromOptions_HPDDM; 3045 pc->ops->setup = PCSetUp_HPDDM; 3046 pc->ops->apply = PCApply_HPDDM; 3047 pc->ops->matapply = PCMatApply_HPDDM; 3048 pc->ops->view = PCView_HPDDM; 3049 pc->ops->presolve = PCPreSolve_HPDDM; 3050 3051 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetAuxiliaryMat_C", PCHPDDMSetAuxiliaryMat_HPDDM)); 3052 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMHasNeumannMat_C", PCHPDDMHasNeumannMat_HPDDM)); 3053 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetRHSMat_C", PCHPDDMSetRHSMat_HPDDM)); 3054 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetCoarseCorrectionType_C", PCHPDDMSetCoarseCorrectionType_HPDDM)); 3055 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMGetCoarseCorrectionType_C", PCHPDDMGetCoarseCorrectionType_HPDDM)); 3056 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetSTShareSubKSP_C", PCHPDDMSetSTShareSubKSP_HPDDM)); 3057 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMGetSTShareSubKSP_C", PCHPDDMGetSTShareSubKSP_HPDDM)); 3058 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetDeflationMat_C", PCHPDDMSetDeflationMat_HPDDM)); 3059 PetscFunctionReturn(PETSC_SUCCESS); 3060 } 3061 3062 /*@C 3063 PCHPDDMInitializePackage - This function initializes everything in the `PCHPDDM` package. It is called from `PCInitializePackage()`. 3064 3065 Level: developer 3066 3067 .seealso: [](ch_ksp), `PetscInitialize()` 3068 @*/ 3069 PetscErrorCode PCHPDDMInitializePackage(void) 3070 { 3071 char ename[32]; 3072 3073 PetscFunctionBegin; 3074 if (PCHPDDMPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS); 3075 PCHPDDMPackageInitialized = PETSC_TRUE; 3076 PetscCall(PetscRegisterFinalize(PCHPDDMFinalizePackage)); 3077 /* general events registered once during package initialization */ 3078 /* some of these events are not triggered in libpetsc, */ 3079 /* but rather directly in libhpddm_petsc, */ 3080 /* which is in charge of performing the following operations */ 3081 3082 /* domain decomposition structure from Pmat sparsity pattern */ 3083 PetscCall(PetscLogEventRegister("PCHPDDMStrc", PC_CLASSID, &PC_HPDDM_Strc)); 3084 /* Galerkin product, redistribution, and setup (not triggered in libpetsc) */ 3085 PetscCall(PetscLogEventRegister("PCHPDDMPtAP", PC_CLASSID, &PC_HPDDM_PtAP)); 3086 /* Galerkin product with summation, redistribution, and setup (not triggered in libpetsc) */ 3087 PetscCall(PetscLogEventRegister("PCHPDDMPtBP", PC_CLASSID, &PC_HPDDM_PtBP)); 3088 /* next level construction using PtAP and PtBP (not triggered in libpetsc) */ 3089 PetscCall(PetscLogEventRegister("PCHPDDMNext", PC_CLASSID, &PC_HPDDM_Next)); 3090 static_assert(PETSC_PCHPDDM_MAXLEVELS <= 9, "PETSC_PCHPDDM_MAXLEVELS value is too high"); 3091 for (PetscInt i = 1; i < PETSC_PCHPDDM_MAXLEVELS; ++i) { 3092 PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCHPDDMSetUp L%1" PetscInt_FMT, i)); 3093 /* events during a PCSetUp() at level #i _except_ the assembly */ 3094 /* of the Galerkin operator of the coarser level #(i + 1) */ 3095 PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_HPDDM_SetUp[i - 1])); 3096 PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCHPDDMSolve L%1" PetscInt_FMT, i)); 3097 /* events during a PCApply() at level #i _except_ */ 3098 /* the KSPSolve() of the coarser level #(i + 1) */ 3099 PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_HPDDM_Solve[i - 1])); 3100 } 3101 PetscFunctionReturn(PETSC_SUCCESS); 3102 } 3103 3104 /*@C 3105 PCHPDDMFinalizePackage - This function frees everything from the `PCHPDDM` package. It is called from `PetscFinalize()`. 3106 3107 Level: developer 3108 3109 .seealso: [](ch_ksp), `PetscFinalize()` 3110 @*/ 3111 PetscErrorCode PCHPDDMFinalizePackage(void) 3112 { 3113 PetscFunctionBegin; 3114 PCHPDDMPackageInitialized = PETSC_FALSE; 3115 PetscFunctionReturn(PETSC_SUCCESS); 3116 } 3117 3118 static PetscErrorCode MatMult_Harmonic(Mat A, Vec x, Vec y) 3119 { 3120 Harmonic h; /* [ A_00 A_01 ], furthermore, A_00 = [ A_loc,loc A_loc,ovl ], thus, A_01 = [ ] */ 3121 /* [ A_10 A_11 A_12 ] [ A_ovl,loc A_ovl,ovl ] [ A_ovl,1 ] */ 3122 Vec sub; /* y = A x = R_loc R_0 [ A_00 A_01 ]^-1 R_loc = [ I_loc ] */ 3123 /* [ A_10 A_11 ] R_1^T A_12 x [ ] */ 3124 PetscFunctionBegin; 3125 PetscCall(MatShellGetContext(A, &h)); 3126 PetscCall(VecSet(h->v, 0.0)); 3127 PetscCall(VecGetSubVector(h->v, h->is[0], &sub)); 3128 PetscCall(MatMult(h->A[0], x, sub)); 3129 PetscCall(VecRestoreSubVector(h->v, h->is[0], &sub)); 3130 PetscCall(KSPSolve(h->ksp, h->v, h->v)); 3131 PetscCall(VecISCopy(h->v, h->is[1], SCATTER_REVERSE, y)); 3132 PetscFunctionReturn(PETSC_SUCCESS); 3133 } 3134 3135 static PetscErrorCode MatMultTranspose_Harmonic(Mat A, Vec y, Vec x) 3136 { 3137 Harmonic h; /* x = A^T y = [ A_00 A_01 ]^-T R_0^T R_loc^T y */ 3138 Vec sub; /* A_12^T R_1 [ A_10 A_11 ] */ 3139 3140 PetscFunctionBegin; 3141 PetscCall(MatShellGetContext(A, &h)); 3142 PetscCall(VecSet(h->v, 0.0)); 3143 PetscCall(VecISCopy(h->v, h->is[1], SCATTER_FORWARD, y)); 3144 PetscCall(KSPSolveTranspose(h->ksp, h->v, h->v)); 3145 PetscCall(VecGetSubVector(h->v, h->is[0], &sub)); 3146 PetscCall(MatMultTranspose(h->A[0], sub, x)); 3147 PetscCall(VecRestoreSubVector(h->v, h->is[0], &sub)); 3148 PetscFunctionReturn(PETSC_SUCCESS); 3149 } 3150 3151 static PetscErrorCode MatProduct_AB_Harmonic(Mat S, Mat X, Mat Y, void *) 3152 { 3153 Harmonic h; 3154 Mat A, B; 3155 Vec a, b; 3156 3157 PetscFunctionBegin; 3158 PetscCall(MatShellGetContext(S, &h)); 3159 PetscCall(MatMatMult(h->A[0], X, MAT_INITIAL_MATRIX, PETSC_CURRENT, &A)); 3160 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, h->ksp->pc->mat->rmap->n, A->cmap->n, nullptr, &B)); 3161 for (PetscInt i = 0; i < A->cmap->n; ++i) { 3162 PetscCall(MatDenseGetColumnVecRead(A, i, &a)); 3163 PetscCall(MatDenseGetColumnVecWrite(B, i, &b)); 3164 PetscCall(VecISCopy(b, h->is[0], SCATTER_FORWARD, a)); 3165 PetscCall(MatDenseRestoreColumnVecWrite(B, i, &b)); 3166 PetscCall(MatDenseRestoreColumnVecRead(A, i, &a)); 3167 } 3168 PetscCall(MatDestroy(&A)); 3169 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, h->ksp->pc->mat->rmap->n, B->cmap->n, nullptr, &A)); 3170 PetscCall(KSPMatSolve(h->ksp, B, A)); 3171 PetscCall(MatDestroy(&B)); 3172 for (PetscInt i = 0; i < A->cmap->n; ++i) { 3173 PetscCall(MatDenseGetColumnVecRead(A, i, &a)); 3174 PetscCall(MatDenseGetColumnVecWrite(Y, i, &b)); 3175 PetscCall(VecISCopy(a, h->is[1], SCATTER_REVERSE, b)); 3176 PetscCall(MatDenseRestoreColumnVecWrite(Y, i, &b)); 3177 PetscCall(MatDenseRestoreColumnVecRead(A, i, &a)); 3178 } 3179 PetscCall(MatDestroy(&A)); 3180 PetscFunctionReturn(PETSC_SUCCESS); 3181 } 3182 3183 static PetscErrorCode MatDestroy_Harmonic(Mat A) 3184 { 3185 Harmonic h; 3186 3187 PetscFunctionBegin; 3188 PetscCall(MatShellGetContext(A, &h)); 3189 for (PetscInt i = 0; i < 5; ++i) PetscCall(ISDestroy(h->is + i)); 3190 PetscCall(PetscFree(h->is)); 3191 PetscCall(VecDestroy(&h->v)); 3192 for (PetscInt i = 0; i < 2; ++i) PetscCall(MatDestroy(h->A + i)); 3193 PetscCall(PetscFree(h->A)); 3194 PetscCall(KSPDestroy(&h->ksp)); 3195 PetscCall(PetscFree(h)); 3196 PetscFunctionReturn(PETSC_SUCCESS); 3197 } 3198