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 = PETSC_FALSE; 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 PetscCall(MatMumpsSetIcntl(A, 26, 0)); 1023 } else { 1024 PetscCall(PetscStrcmp(type, MATSOLVERMKL_PARDISO, &flg)); 1025 PetscCheck(flg && PetscDefined(HAVE_MKL_PARDISO), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent MatSolverType"); 1026 flg = PETSC_FALSE; 1027 #if PetscDefined(HAVE_MKL_PARDISO) 1028 PetscCall(MatMkl_PardisoSetCntl(A, 70, 1)); 1029 #endif 1030 } 1031 PetscCall(PCApply_Schur_Private<Type, T>(p, factor, x, y)); 1032 if (flg) { 1033 PetscCall(MatMumpsSetIcntl(A, 26, -1)); 1034 } else { 1035 #if PetscDefined(HAVE_MKL_PARDISO) 1036 PetscCall(MatMkl_PardisoSetCntl(A, 70, 0)); 1037 #endif 1038 } 1039 PetscFunctionReturn(PETSC_SUCCESS); 1040 } 1041 1042 static PetscErrorCode PCDestroy_Schur(PC pc) 1043 { 1044 std::tuple<KSP, IS, Vec[2]> *p; 1045 1046 PetscFunctionBegin; 1047 PetscCall(PCShellGetContext(pc, &p)); 1048 PetscCall(ISDestroy(&std::get<1>(*p))); 1049 PetscCall(VecDestroy(std::get<2>(*p))); 1050 PetscCall(VecDestroy(std::get<2>(*p) + 1)); 1051 PetscCall(PetscFree(p)); 1052 PetscFunctionReturn(PETSC_SUCCESS); 1053 } 1054 1055 static PetscErrorCode PCHPDDMSolve_Private(const PC_HPDDM_Level *ctx, PetscScalar *rhs, const unsigned short &mu) 1056 { 1057 Mat B, X; 1058 PetscInt n, N, j = 0; 1059 1060 PetscFunctionBegin; 1061 PetscCall(KSPGetOperators(ctx->ksp, &B, nullptr)); 1062 PetscCall(MatGetLocalSize(B, &n, nullptr)); 1063 PetscCall(MatGetSize(B, &N, nullptr)); 1064 if (ctx->parent->log_separate) { 1065 j = std::distance(ctx->parent->levels, std::find(ctx->parent->levels, ctx->parent->levels + ctx->parent->N, ctx)); 1066 PetscCall(PetscLogEventBegin(PC_HPDDM_Solve[j], ctx->ksp, nullptr, nullptr, nullptr)); 1067 } 1068 if (mu == 1) { 1069 if (!ctx->ksp->vec_rhs) { 1070 PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)ctx->ksp), 1, n, N, nullptr, &ctx->ksp->vec_rhs)); 1071 PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)ctx->ksp), n, N, &ctx->ksp->vec_sol)); 1072 } 1073 PetscCall(VecPlaceArray(ctx->ksp->vec_rhs, rhs)); 1074 PetscCall(KSPSolve(ctx->ksp, nullptr, nullptr)); 1075 PetscCall(VecCopy(ctx->ksp->vec_sol, ctx->ksp->vec_rhs)); 1076 PetscCall(VecResetArray(ctx->ksp->vec_rhs)); 1077 } else { 1078 PetscCall(MatCreateDense(PetscObjectComm((PetscObject)ctx->ksp), n, PETSC_DECIDE, N, mu, rhs, &B)); 1079 PetscCall(MatCreateDense(PetscObjectComm((PetscObject)ctx->ksp), n, PETSC_DECIDE, N, mu, nullptr, &X)); 1080 PetscCall(KSPMatSolve(ctx->ksp, B, X)); 1081 PetscCall(MatCopy(X, B, SAME_NONZERO_PATTERN)); 1082 PetscCall(MatDestroy(&X)); 1083 PetscCall(MatDestroy(&B)); 1084 } 1085 if (ctx->parent->log_separate) PetscCall(PetscLogEventEnd(PC_HPDDM_Solve[j], ctx->ksp, nullptr, nullptr, nullptr)); 1086 PetscFunctionReturn(PETSC_SUCCESS); 1087 } 1088 1089 static PetscErrorCode PCHPDDMSetUpNeumannOverlap_Private(PC pc) 1090 { 1091 PC_HPDDM *data = (PC_HPDDM *)pc->data; 1092 1093 PetscFunctionBegin; 1094 if (data->setup) { 1095 Mat P; 1096 Vec x, xt = nullptr; 1097 PetscReal t = 0.0, s = 0.0; 1098 1099 PetscCall(PCGetOperators(pc, nullptr, &P)); 1100 PetscCall(PetscObjectQuery((PetscObject)P, "__SNES_latest_X", (PetscObject *)&x)); 1101 PetscCallBack("PCHPDDM Neumann callback", (*data->setup)(data->aux, t, x, xt, s, data->is, data->setup_ctx)); 1102 } 1103 PetscFunctionReturn(PETSC_SUCCESS); 1104 } 1105 1106 static PetscErrorCode PCHPDDMCreateSubMatrices_Private(Mat mat, PetscInt n, const IS *, const IS *, MatReuse scall, Mat *submat[]) 1107 { 1108 Mat A; 1109 PetscBool flg; 1110 1111 PetscFunctionBegin; 1112 PetscCheck(n == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "MatCreateSubMatrices() called to extract %" PetscInt_FMT " submatrices, which is different than 1", n); 1113 /* previously composed Mat */ 1114 PetscCall(PetscObjectQuery((PetscObject)mat, "_PCHPDDM_SubMatrices", (PetscObject *)&A)); 1115 PetscCheck(A, PETSC_COMM_SELF, PETSC_ERR_PLIB, "SubMatrices not found in Mat"); 1116 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATSCHURCOMPLEMENT, &flg)); /* MATSCHURCOMPLEMENT has neither a MatDuplicate() nor a MatCopy() implementation */ 1117 if (scall == MAT_INITIAL_MATRIX) { 1118 PetscCall(PetscCalloc1(2, submat)); /* allocate an extra Mat to avoid errors in MatDestroySubMatrices_Dummy() */ 1119 if (!flg) PetscCall(MatDuplicate(A, MAT_COPY_VALUES, *submat)); 1120 } else if (!flg) PetscCall(MatCopy(A, (*submat)[0], SAME_NONZERO_PATTERN)); 1121 if (flg) { 1122 PetscCall(MatDestroy(*submat)); /* previously created Mat has to be destroyed */ 1123 (*submat)[0] = A; 1124 PetscCall(PetscObjectReference((PetscObject)A)); 1125 } 1126 PetscFunctionReturn(PETSC_SUCCESS); 1127 } 1128 1129 static PetscErrorCode PCHPDDMCommunicationAvoidingPCASM_Private(PC pc, Mat C, PetscBool sorted) 1130 { 1131 void (*op)(void); 1132 1133 PetscFunctionBegin; 1134 /* previously-composed Mat */ 1135 PetscCall(PetscObjectCompose((PetscObject)pc->pmat, "_PCHPDDM_SubMatrices", (PetscObject)C)); 1136 PetscCall(MatGetOperation(pc->pmat, MATOP_CREATE_SUBMATRICES, &op)); 1137 /* trick suggested by Barry https://lists.mcs.anl.gov/pipermail/petsc-dev/2020-January/025491.html */ 1138 PetscCall(MatSetOperation(pc->pmat, MATOP_CREATE_SUBMATRICES, (void (*)(void))PCHPDDMCreateSubMatrices_Private)); 1139 if (sorted) PetscCall(PCASMSetSortIndices(pc, PETSC_FALSE)); /* everything is already sorted */ 1140 PetscCall(PCSetFromOptions(pc)); /* otherwise -pc_hpddm_levels_1_pc_asm_sub_mat_type is not used */ 1141 PetscCall(PCSetUp(pc)); 1142 /* reset MatCreateSubMatrices() */ 1143 PetscCall(MatSetOperation(pc->pmat, MATOP_CREATE_SUBMATRICES, op)); 1144 PetscCall(PetscObjectCompose((PetscObject)pc->pmat, "_PCHPDDM_SubMatrices", nullptr)); 1145 PetscFunctionReturn(PETSC_SUCCESS); 1146 } 1147 1148 static PetscErrorCode PCHPDDMPermute_Private(IS is, IS in_is, IS *out_is, Mat in_C, Mat *out_C, IS *p) 1149 { 1150 IS perm; 1151 const PetscInt *ptr; 1152 PetscInt *concatenate, size, bs; 1153 std::map<PetscInt, PetscInt> order; 1154 PetscBool sorted; 1155 1156 PetscFunctionBegin; 1157 PetscValidHeaderSpecific(is, IS_CLASSID, 1); 1158 PetscValidHeaderSpecific(in_C, MAT_CLASSID, 4); 1159 PetscCall(ISSorted(is, &sorted)); 1160 if (!sorted) { 1161 PetscCall(ISGetLocalSize(is, &size)); 1162 PetscCall(ISGetIndices(is, &ptr)); 1163 PetscCall(ISGetBlockSize(is, &bs)); 1164 /* MatCreateSubMatrices(), called by PCASM, follows the global numbering of Pmat */ 1165 for (PetscInt n = 0; n < size; n += bs) order.insert(std::make_pair(ptr[n] / bs, n / bs)); 1166 PetscCall(ISRestoreIndices(is, &ptr)); 1167 size /= bs; 1168 if (out_C) { 1169 PetscCall(PetscMalloc1(size, &concatenate)); 1170 for (const std::pair<const PetscInt, PetscInt> &i : order) *concatenate++ = i.second; 1171 concatenate -= size; 1172 PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)in_C), bs, size, concatenate, PETSC_OWN_POINTER, &perm)); 1173 PetscCall(ISSetPermutation(perm)); 1174 /* permute user-provided Mat so that it matches with MatCreateSubMatrices() numbering */ 1175 PetscCall(MatPermute(in_C, perm, perm, out_C)); 1176 if (p) *p = perm; 1177 else PetscCall(ISDestroy(&perm)); /* no need to save the permutation */ 1178 } 1179 if (out_is) { 1180 PetscCall(PetscMalloc1(size, &concatenate)); 1181 for (const std::pair<const PetscInt, PetscInt> &i : order) *concatenate++ = i.first; 1182 concatenate -= size; 1183 /* permute user-provided IS so that it matches with MatCreateSubMatrices() numbering */ 1184 PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)in_is), bs, size, concatenate, PETSC_OWN_POINTER, out_is)); 1185 } 1186 } else { /* input IS is sorted, nothing to permute, simply duplicate inputs when needed */ 1187 if (out_C) PetscCall(MatDuplicate(in_C, MAT_COPY_VALUES, out_C)); 1188 if (out_is) PetscCall(ISDuplicate(in_is, out_is)); 1189 } 1190 PetscFunctionReturn(PETSC_SUCCESS); 1191 } 1192 1193 static PetscErrorCode PCHPDDMCheckSymmetry_Private(PC pc, Mat A01, Mat A10, Mat *B01 = nullptr) 1194 { 1195 Mat T, U = nullptr, B = nullptr; 1196 IS z; 1197 PetscBool flg, conjugate = PETSC_FALSE; 1198 1199 PetscFunctionBegin; 1200 PetscCall(PetscObjectTypeCompare((PetscObject)A10, MATTRANSPOSEVIRTUAL, &flg)); 1201 if (B01) *B01 = nullptr; 1202 if (flg) { 1203 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)); 1204 PetscCall(MatTransposeGetMat(A10, &U)); 1205 } else { 1206 PetscCall(PetscObjectTypeCompare((PetscObject)A10, MATHERMITIANTRANSPOSEVIRTUAL, &flg)); 1207 if (flg) { 1208 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)); 1209 PetscCall(MatHermitianTransposeGetMat(A10, &U)); 1210 conjugate = PETSC_TRUE; 1211 } 1212 } 1213 if (U) PetscCall(MatDuplicate(U, MAT_COPY_VALUES, &T)); 1214 else PetscCall(MatHermitianTranspose(A10, MAT_INITIAL_MATRIX, &T)); 1215 PetscCall(PetscObjectTypeCompare((PetscObject)A01, MATTRANSPOSEVIRTUAL, &flg)); 1216 if (flg) { 1217 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)); 1218 PetscCall(MatTransposeGetMat(A01, &A01)); 1219 PetscCall(MatTranspose(A01, MAT_INITIAL_MATRIX, &B)); 1220 A01 = B; 1221 } else { 1222 PetscCall(PetscObjectTypeCompare((PetscObject)A01, MATHERMITIANTRANSPOSEVIRTUAL, &flg)); 1223 if (flg) { 1224 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)); 1225 PetscCall(MatHermitianTransposeGetMat(A01, &A01)); 1226 PetscCall(MatHermitianTranspose(A01, MAT_INITIAL_MATRIX, &B)); 1227 A01 = B; 1228 } 1229 } 1230 PetscCall(PetscLayoutCompare(T->rmap, A01->rmap, &flg)); 1231 if (flg) { 1232 PetscCall(PetscLayoutCompare(T->cmap, A01->cmap, &flg)); 1233 if (flg) { 1234 PetscCall(MatFindZeroRows(A01, &z)); /* for essential boundary conditions, some implementations will */ 1235 if (z) { /* zero rows in [P00 A01] except for the diagonal of P00 */ 1236 if (B01) PetscCall(MatDuplicate(T, MAT_COPY_VALUES, B01)); 1237 PetscCall(MatSetOption(T, MAT_NO_OFF_PROC_ZERO_ROWS, PETSC_TRUE)); 1238 PetscCall(MatZeroRowsIS(T, z, 0.0, nullptr, nullptr)); /* corresponding zero rows from A01 */ 1239 } 1240 PetscCall(MatMultEqual(A01, T, 20, &flg)); 1241 if (!B01) PetscCheck(flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "A01 != A10^T"); 1242 else { 1243 PetscCall(PetscInfo(pc, "A01 and A10^T are equal? %s\n", PetscBools[flg])); 1244 if (!flg) { 1245 if (z) PetscCall(MatDestroy(&T)); 1246 else *B01 = T; 1247 flg = PETSC_TRUE; 1248 } else PetscCall(MatDestroy(B01)); 1249 } 1250 PetscCall(ISDestroy(&z)); 1251 } 1252 } 1253 if (!flg) PetscCall(PetscInfo(pc, "A01 and A10^T have non-congruent layouts, cannot test for equality\n")); 1254 if (!B01 || !*B01) PetscCall(MatDestroy(&T)); 1255 else if (conjugate) PetscCall(MatConjugate(T)); 1256 PetscCall(MatDestroy(&B)); 1257 PetscFunctionReturn(PETSC_SUCCESS); 1258 } 1259 1260 static PetscErrorCode PCHPDDMCheckInclusion_Private(PC pc, IS is, IS is_local, PetscBool check) 1261 { 1262 IS intersect; 1263 const char *str = "IS of the auxiliary Mat does not include all local rows of A"; 1264 PetscBool equal; 1265 1266 PetscFunctionBegin; 1267 PetscCall(ISIntersect(is, is_local, &intersect)); 1268 PetscCall(ISEqualUnsorted(is_local, intersect, &equal)); 1269 PetscCall(ISDestroy(&intersect)); 1270 if (check) PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "%s", str); 1271 else if (!equal) PetscCall(PetscInfo(pc, "%s\n", str)); 1272 PetscFunctionReturn(PETSC_SUCCESS); 1273 } 1274 1275 static PetscErrorCode PCHPDDMDestroySubMatrices_Private(PetscBool flg, PetscBool algebraic, Mat *sub) 1276 { 1277 IS is; 1278 1279 PetscFunctionBegin; 1280 if (!flg) { 1281 if (algebraic) { 1282 PetscCall(PetscObjectQuery((PetscObject)sub[0], "_PCHPDDM_Embed", (PetscObject *)&is)); 1283 PetscCall(ISDestroy(&is)); 1284 PetscCall(PetscObjectCompose((PetscObject)sub[0], "_PCHPDDM_Embed", nullptr)); 1285 PetscCall(PetscObjectCompose((PetscObject)sub[0], "_PCHPDDM_Compact", nullptr)); 1286 } 1287 PetscCall(MatDestroySubMatrices(algebraic ? 2 : 1, &sub)); 1288 } 1289 PetscFunctionReturn(PETSC_SUCCESS); 1290 } 1291 1292 static PetscErrorCode PCHPDDMAlgebraicAuxiliaryMat_Private(Mat P, IS *is, Mat *sub[], PetscBool block) 1293 { 1294 IS icol[3], irow[2]; 1295 Mat *M, Q; 1296 PetscReal *ptr; 1297 PetscInt *idx, p = 0, bs = P->cmap->bs; 1298 PetscBool flg; 1299 1300 PetscFunctionBegin; 1301 PetscCall(ISCreateStride(PETSC_COMM_SELF, P->cmap->N, 0, 1, icol + 2)); 1302 PetscCall(ISSetBlockSize(icol[2], bs)); 1303 PetscCall(ISSetIdentity(icol[2])); 1304 PetscCall(PetscObjectTypeCompare((PetscObject)P, MATMPISBAIJ, &flg)); 1305 if (flg) { 1306 /* MatCreateSubMatrices() does not handle MATMPISBAIJ properly when iscol != isrow, so convert first to MATMPIBAIJ */ 1307 PetscCall(MatConvert(P, MATMPIBAIJ, MAT_INITIAL_MATRIX, &Q)); 1308 std::swap(P, Q); 1309 } 1310 PetscCall(MatCreateSubMatrices(P, 1, is, icol + 2, MAT_INITIAL_MATRIX, &M)); 1311 if (flg) { 1312 std::swap(P, Q); 1313 PetscCall(MatDestroy(&Q)); 1314 } 1315 PetscCall(ISDestroy(icol + 2)); 1316 PetscCall(ISCreateStride(PETSC_COMM_SELF, M[0]->rmap->N, 0, 1, irow)); 1317 PetscCall(ISSetBlockSize(irow[0], bs)); 1318 PetscCall(ISSetIdentity(irow[0])); 1319 if (!block) { 1320 PetscCall(PetscMalloc2(P->cmap->N, &ptr, P->cmap->N / bs, &idx)); 1321 PetscCall(MatGetColumnNorms(M[0], NORM_INFINITY, ptr)); 1322 /* check for nonzero columns so that M[0] may be expressed in compact form */ 1323 for (PetscInt n = 0; n < P->cmap->N; n += bs) { 1324 if (std::find_if(ptr + n, ptr + n + bs, [](PetscReal v) { return v > PETSC_MACHINE_EPSILON; }) != ptr + n + bs) idx[p++] = n / bs; 1325 } 1326 PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, p, idx, PETSC_USE_POINTER, icol + 1)); 1327 PetscCall(ISSetInfo(icol[1], IS_SORTED, IS_GLOBAL, PETSC_TRUE, PETSC_TRUE)); 1328 PetscCall(ISEmbed(*is, icol[1], PETSC_FALSE, icol + 2)); 1329 irow[1] = irow[0]; 1330 /* 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 */ 1331 icol[0] = is[0]; 1332 PetscCall(MatCreateSubMatrices(M[0], 2, irow, icol, MAT_INITIAL_MATRIX, sub)); 1333 PetscCall(ISDestroy(icol + 1)); 1334 PetscCall(PetscFree2(ptr, idx)); 1335 /* IS used to go back and forth between the augmented and the original local linear system, see eq. (3.4) of [2022b] */ 1336 PetscCall(PetscObjectCompose((PetscObject)(*sub)[0], "_PCHPDDM_Embed", (PetscObject)icol[2])); 1337 /* Mat used in eq. (3.1) of [2022b] */ 1338 PetscCall(PetscObjectCompose((PetscObject)(*sub)[0], "_PCHPDDM_Compact", (PetscObject)(*sub)[1])); 1339 } else { 1340 Mat aux; 1341 1342 PetscCall(MatSetOption(M[0], MAT_SUBMAT_SINGLEIS, PETSC_TRUE)); 1343 /* diagonal block of the overlapping rows */ 1344 PetscCall(MatCreateSubMatrices(M[0], 1, irow, is, MAT_INITIAL_MATRIX, sub)); 1345 PetscCall(MatDuplicate((*sub)[0], MAT_COPY_VALUES, &aux)); 1346 PetscCall(MatSetOption(aux, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE)); 1347 if (bs == 1) { /* scalar case */ 1348 Vec sum[2]; 1349 1350 PetscCall(MatCreateVecs(aux, sum, sum + 1)); 1351 PetscCall(MatGetRowSum(M[0], sum[0])); 1352 PetscCall(MatGetRowSum(aux, sum[1])); 1353 /* off-diagonal block row sum (full rows - diagonal block rows) */ 1354 PetscCall(VecAXPY(sum[0], -1.0, sum[1])); 1355 /* subdomain matrix plus off-diagonal block row sum */ 1356 PetscCall(MatDiagonalSet(aux, sum[0], ADD_VALUES)); 1357 PetscCall(VecDestroy(sum)); 1358 PetscCall(VecDestroy(sum + 1)); 1359 } else { /* vectorial case */ 1360 /* TODO: missing MatGetValuesBlocked(), so the code below is */ 1361 /* an extension of the scalar case for when bs > 1, but it could */ 1362 /* be more efficient by avoiding all these MatMatMult() */ 1363 Mat sum[2], ones; 1364 PetscScalar *ptr; 1365 1366 PetscCall(PetscCalloc1(M[0]->cmap->n * bs, &ptr)); 1367 PetscCall(MatCreateDense(PETSC_COMM_SELF, M[0]->cmap->n, bs, M[0]->cmap->n, bs, ptr, &ones)); 1368 for (PetscInt n = 0; n < M[0]->cmap->n; n += bs) { 1369 for (p = 0; p < bs; ++p) ptr[n + p * (M[0]->cmap->n + 1)] = 1.0; 1370 } 1371 PetscCall(MatMatMult(M[0], ones, MAT_INITIAL_MATRIX, PETSC_CURRENT, sum)); 1372 PetscCall(MatDestroy(&ones)); 1373 PetscCall(MatCreateDense(PETSC_COMM_SELF, aux->cmap->n, bs, aux->cmap->n, bs, ptr, &ones)); 1374 PetscCall(MatDenseSetLDA(ones, M[0]->cmap->n)); 1375 PetscCall(MatMatMult(aux, ones, MAT_INITIAL_MATRIX, PETSC_CURRENT, sum + 1)); 1376 PetscCall(MatDestroy(&ones)); 1377 PetscCall(PetscFree(ptr)); 1378 /* off-diagonal block row sum (full rows - diagonal block rows) */ 1379 PetscCall(MatAXPY(sum[0], -1.0, sum[1], SAME_NONZERO_PATTERN)); 1380 PetscCall(MatDestroy(sum + 1)); 1381 /* re-order values to be consistent with MatSetValuesBlocked() */ 1382 /* equivalent to MatTranspose() which does not truly handle */ 1383 /* MAT_INPLACE_MATRIX in the rectangular case, as it calls PetscMalloc() */ 1384 PetscCall(MatDenseGetArrayWrite(sum[0], &ptr)); 1385 HPDDM::Wrapper<PetscScalar>::imatcopy<'T'>(bs, sum[0]->rmap->n, ptr, sum[0]->rmap->n, bs); 1386 /* subdomain matrix plus off-diagonal block row sum */ 1387 for (PetscInt n = 0; n < aux->cmap->n / bs; ++n) PetscCall(MatSetValuesBlocked(aux, 1, &n, 1, &n, ptr + n * bs * bs, ADD_VALUES)); 1388 PetscCall(MatAssemblyBegin(aux, MAT_FINAL_ASSEMBLY)); 1389 PetscCall(MatAssemblyEnd(aux, MAT_FINAL_ASSEMBLY)); 1390 PetscCall(MatDenseRestoreArrayWrite(sum[0], &ptr)); 1391 PetscCall(MatDestroy(sum)); 1392 } 1393 PetscCall(MatSetOption(aux, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE)); 1394 /* left-hand side of GenEO, with the same sparsity pattern as PCASM subdomain solvers */ 1395 PetscCall(PetscObjectCompose((PetscObject)(*sub)[0], "_PCHPDDM_Neumann_Mat", (PetscObject)aux)); 1396 } 1397 PetscCall(ISDestroy(irow)); 1398 PetscCall(MatDestroySubMatrices(1, &M)); 1399 PetscFunctionReturn(PETSC_SUCCESS); 1400 } 1401 1402 static PetscErrorCode PCApply_Nest(PC pc, Vec x, Vec y) 1403 { 1404 Mat A; 1405 MatSolverType type; 1406 IS is[2]; 1407 PetscBool flg; 1408 std::pair<PC, Vec[2]> *p; 1409 1410 PetscFunctionBegin; 1411 PetscCall(PCShellGetContext(pc, &p)); 1412 if (p->second[0]) { /* in case of a centralized Schur complement, some processes may have no local operator */ 1413 PetscCall(PCGetOperators(p->first, &A, nullptr)); 1414 PetscCall(MatNestGetISs(A, is, nullptr)); 1415 PetscCall(PetscObjectTypeCompareAny((PetscObject)p->first, &flg, PCLU, PCCHOLESKY, "")); 1416 if (flg) { /* partial solve currently only makes sense with exact factorizations */ 1417 PetscCall(PCFactorGetMatSolverType(p->first, &type)); 1418 PetscCall(PCFactorGetMatrix(p->first, &A)); 1419 if (A->schur) { 1420 PetscCall(PetscStrcmp(type, MATSOLVERMUMPS, &flg)); 1421 if (flg) PetscCall(MatMumpsSetIcntl(A, 26, 1)); /* reduction/condensation phase followed by Schur complement solve */ 1422 } else flg = PETSC_FALSE; 1423 } 1424 PetscCall(VecISCopy(p->second[0], is[1], SCATTER_FORWARD, x)); /* assign the RHS associated to the Schur complement */ 1425 PetscCall(PCApply(p->first, p->second[0], p->second[1])); 1426 PetscCall(VecISCopy(p->second[1], is[1], SCATTER_REVERSE, y)); /* retrieve the partial solution associated to the Schur complement */ 1427 if (flg) PetscCall(MatMumpsSetIcntl(A, 26, -1)); /* default ICNTL(26) value in PETSc */ 1428 } 1429 PetscFunctionReturn(PETSC_SUCCESS); 1430 } 1431 1432 static PetscErrorCode PCView_Nest(PC pc, PetscViewer viewer) 1433 { 1434 std::pair<PC, Vec[2]> *p; 1435 1436 PetscFunctionBegin; 1437 PetscCall(PCShellGetContext(pc, &p)); 1438 PetscCall(PCView(p->first, viewer)); 1439 PetscFunctionReturn(PETSC_SUCCESS); 1440 } 1441 1442 static PetscErrorCode PCDestroy_Nest(PC pc) 1443 { 1444 std::pair<PC, Vec[2]> *p; 1445 1446 PetscFunctionBegin; 1447 PetscCall(PCShellGetContext(pc, &p)); 1448 PetscCall(VecDestroy(p->second)); 1449 PetscCall(VecDestroy(p->second + 1)); 1450 PetscCall(PCDestroy(&p->first)); 1451 PetscCall(PetscFree(p)); 1452 PetscFunctionReturn(PETSC_SUCCESS); 1453 } 1454 1455 template <bool T = false> 1456 static PetscErrorCode MatMult_Schur(Mat A, Vec x, Vec y) 1457 { 1458 std::tuple<Mat, PetscSF, Vec[2]> *ctx; 1459 1460 PetscFunctionBegin; 1461 PetscCall(MatShellGetContext(A, &ctx)); 1462 PetscCall(VecScatterBegin(std::get<1>(*ctx), x, std::get<2>(*ctx)[0], INSERT_VALUES, SCATTER_FORWARD)); /* local Vec with overlap */ 1463 PetscCall(VecScatterEnd(std::get<1>(*ctx), x, std::get<2>(*ctx)[0], INSERT_VALUES, SCATTER_FORWARD)); 1464 if (!T) PetscCall(MatMult(std::get<0>(*ctx), std::get<2>(*ctx)[0], std::get<2>(*ctx)[1])); /* local Schur complement */ 1465 else PetscCall(MatMultTranspose(std::get<0>(*ctx), std::get<2>(*ctx)[0], std::get<2>(*ctx)[1])); 1466 PetscCall(VecSet(y, 0.0)); 1467 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 */ 1468 PetscCall(VecScatterEnd(std::get<1>(*ctx), std::get<2>(*ctx)[1], y, ADD_VALUES, SCATTER_REVERSE)); 1469 PetscFunctionReturn(PETSC_SUCCESS); 1470 } 1471 1472 static PetscErrorCode MatDestroy_Schur(Mat A) 1473 { 1474 std::tuple<Mat, PetscSF, Vec[2]> *ctx; 1475 1476 PetscFunctionBegin; 1477 PetscCall(MatShellGetContext(A, &ctx)); 1478 PetscCall(VecDestroy(std::get<2>(*ctx))); 1479 PetscCall(VecDestroy(std::get<2>(*ctx) + 1)); 1480 PetscCall(PetscFree(ctx)); 1481 PetscFunctionReturn(PETSC_SUCCESS); 1482 } 1483 1484 static PetscErrorCode MatMult_SchurCorrection(Mat A, Vec x, Vec y) 1485 { 1486 PC pc; 1487 std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *ctx; 1488 1489 PetscFunctionBegin; 1490 PetscCall(MatShellGetContext(A, &ctx)); 1491 pc = ((PC_HPDDM *)std::get<0>(*ctx)[0]->data)->levels[0]->ksp->pc; 1492 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 */ 1493 PetscCall(MatMult(std::get<1>(*ctx)[0], x, std::get<3>(*ctx)[1])); /* A_01 x */ 1494 PetscCall(PCHPDDMDeflate_Private(pc, std::get<3>(*ctx)[1], std::get<3>(*ctx)[1])); /* Q_0 A_01 x */ 1495 PetscCall(MatMult(std::get<1>(*ctx)[1], std::get<3>(*ctx)[1], std::get<3>(*ctx)[0])); /* A_10 Q_0 A_01 x */ 1496 PetscCall(PCApply(std::get<0>(*ctx)[1], std::get<3>(*ctx)[0], y)); /* y = M_S^-1 A_10 Q_0 A_01 x */ 1497 } else { 1498 PetscCall(PCApply(std::get<0>(*ctx)[1], x, std::get<3>(*ctx)[0])); /* M_S^-1 x */ 1499 PetscCall(MatMult(std::get<1>(*ctx)[0], std::get<3>(*ctx)[0], std::get<3>(*ctx)[1])); /* A_01 M_S^-1 x */ 1500 PetscCall(PCHPDDMDeflate_Private(pc, std::get<3>(*ctx)[1], std::get<3>(*ctx)[1])); /* Q_0 A_01 M_S^-1 x */ 1501 PetscCall(MatMult(std::get<1>(*ctx)[1], std::get<3>(*ctx)[1], y)); /* y = A_10 Q_0 A_01 M_S^-1 x */ 1502 } 1503 PetscCall(VecAXPY(y, -1.0, x)); /* y -= x, preconditioned eq. (24) of https://hal.science/hal-02343808v6/document (with a sign flip) */ 1504 PetscFunctionReturn(PETSC_SUCCESS); 1505 } 1506 1507 static PetscErrorCode MatView_SchurCorrection(Mat A, PetscViewer viewer) 1508 { 1509 PetscBool ascii; 1510 std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *ctx; 1511 1512 PetscFunctionBegin; 1513 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &ascii)); 1514 if (ascii) { 1515 PetscCall(MatShellGetContext(A, &ctx)); 1516 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)")); 1517 PetscCall(PCView(std::get<0>(*ctx)[1], viewer)); /* no need to PCView(Q_0) since it will be done by PCFIELDSPLIT */ 1518 } 1519 PetscFunctionReturn(PETSC_SUCCESS); 1520 } 1521 1522 static PetscErrorCode MatDestroy_SchurCorrection(Mat A) 1523 { 1524 std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *ctx; 1525 1526 PetscFunctionBegin; 1527 PetscCall(MatShellGetContext(A, &ctx)); 1528 PetscCall(VecDestroy(std::get<3>(*ctx))); 1529 PetscCall(VecDestroy(std::get<3>(*ctx) + 1)); 1530 PetscCall(VecDestroy(std::get<3>(*ctx) + 2)); 1531 PetscCall(PCDestroy(std::get<0>(*ctx) + 1)); 1532 PetscCall(PetscFree(ctx)); 1533 PetscFunctionReturn(PETSC_SUCCESS); 1534 } 1535 1536 static PetscErrorCode PCPostSolve_SchurPreLeastSquares(PC, KSP, Vec, Vec x) 1537 { 1538 PetscFunctionBegin; 1539 PetscCall(VecScale(x, -1.0)); 1540 PetscFunctionReturn(PETSC_SUCCESS); 1541 } 1542 1543 static PetscErrorCode KSPPreSolve_SchurCorrection(KSP, Vec b, Vec, void *context) 1544 { 1545 std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *ctx = reinterpret_cast<std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *>(context); 1546 1547 PetscFunctionBegin; 1548 if (std::get<2>(*ctx) == PC_LEFT || std::get<2>(*ctx) == PC_SIDE_DEFAULT) { 1549 PetscCall(PCApply(std::get<0>(*ctx)[1], b, std::get<3>(*ctx)[2])); 1550 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 */ 1551 } 1552 PetscFunctionReturn(PETSC_SUCCESS); 1553 } 1554 1555 static PetscErrorCode KSPPostSolve_SchurCorrection(KSP, Vec b, Vec x, void *context) 1556 { 1557 std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *ctx = reinterpret_cast<std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *>(context); 1558 1559 PetscFunctionBegin; 1560 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 */ 1561 else { 1562 PetscCall(PCApply(std::get<0>(*ctx)[1], x, std::get<3>(*ctx)[2])); 1563 PetscCall(VecCopy(std::get<3>(*ctx)[2], x)); /* replace x by M^-1 x */ 1564 } 1565 PetscFunctionReturn(PETSC_SUCCESS); 1566 } 1567 1568 static PetscErrorCode MatMult_Harmonic(Mat, Vec, Vec); 1569 static PetscErrorCode MatMultTranspose_Harmonic(Mat, Vec, Vec); 1570 static PetscErrorCode MatProduct_AB_Harmonic(Mat, Mat, Mat, void *); 1571 static PetscErrorCode MatDestroy_Harmonic(Mat); 1572 1573 static PetscErrorCode PCSetUp_HPDDM(PC pc) 1574 { 1575 PC_HPDDM *data = (PC_HPDDM *)pc->data; 1576 PC inner; 1577 KSP *ksp; 1578 Mat *sub, A, P, N, C = nullptr, uaux = nullptr, weighted, subA[2], S; 1579 Vec xin, v; 1580 std::vector<Vec> initial; 1581 IS is[1], loc, uis = data->is, unsorted = nullptr; 1582 ISLocalToGlobalMapping l2g; 1583 char prefix[256]; 1584 const char *pcpre; 1585 const PetscScalar *const *ev; 1586 PetscInt n, requested = data->N, reused = 0, overlap = -1; 1587 MatStructure structure = UNKNOWN_NONZERO_PATTERN; 1588 PetscBool subdomains = PETSC_FALSE, flg = PETSC_FALSE, ismatis, swap = PETSC_FALSE, algebraic = PETSC_FALSE, block = PETSC_FALSE; 1589 DM dm; 1590 std::tuple<PC[2], Mat[2], PCSide, Vec[3]> *ctx = nullptr; 1591 #if PetscDefined(USE_DEBUG) 1592 IS dis = nullptr; 1593 Mat daux = nullptr; 1594 #endif 1595 1596 PetscFunctionBegin; 1597 PetscCheck(data->levels && data->levels[0], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Not a single level allocated"); 1598 PetscCall(PCGetOptionsPrefix(pc, &pcpre)); 1599 PetscCall(PCGetOperators(pc, &A, &P)); 1600 if (!data->levels[0]->ksp) { 1601 PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &data->levels[0]->ksp)); 1602 PetscCall(KSPSetNestLevel(data->levels[0]->ksp, pc->kspnestlevel)); 1603 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_%s_", pcpre ? pcpre : "", data->N > 1 ? "levels_1" : "coarse")); 1604 PetscCall(KSPSetOptionsPrefix(data->levels[0]->ksp, prefix)); 1605 PetscCall(KSPSetType(data->levels[0]->ksp, KSPPREONLY)); 1606 } else if (data->levels[0]->ksp->pc && data->levels[0]->ksp->pc->setupcalled && data->levels[0]->ksp->pc->reusepreconditioner) { 1607 /* if the fine-level PCSHELL exists, its setup has succeeded, and one wants to reuse it, */ 1608 /* then just propagate the appropriate flag to the coarser levels */ 1609 for (n = 0; n < PETSC_PCHPDDM_MAXLEVELS && data->levels[n]; ++n) { 1610 /* the following KSP and PC may be NULL for some processes, hence the check */ 1611 if (data->levels[n]->ksp) PetscCall(KSPSetReusePreconditioner(data->levels[n]->ksp, PETSC_TRUE)); 1612 if (data->levels[n]->pc) PetscCall(PCSetReusePreconditioner(data->levels[n]->pc, PETSC_TRUE)); 1613 } 1614 /* early bail out because there is nothing to do */ 1615 PetscFunctionReturn(PETSC_SUCCESS); 1616 } else { 1617 /* reset coarser levels */ 1618 for (n = 1; n < PETSC_PCHPDDM_MAXLEVELS && data->levels[n]; ++n) { 1619 if (data->levels[n]->ksp && data->levels[n]->ksp->pc && data->levels[n]->ksp->pc->setupcalled && data->levels[n]->ksp->pc->reusepreconditioner && n < data->N) { 1620 reused = data->N - n; 1621 break; 1622 } 1623 PetscCall(KSPDestroy(&data->levels[n]->ksp)); 1624 PetscCall(PCDestroy(&data->levels[n]->pc)); 1625 } 1626 /* check if some coarser levels are being reused */ 1627 PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &reused, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)pc))); 1628 const int *addr = data->levels[0]->P ? data->levels[0]->P->getAddrLocal() : &HPDDM::i__0; 1629 1630 if (addr != &HPDDM::i__0 && reused != data->N - 1) { 1631 /* reuse previously computed eigenvectors */ 1632 ev = data->levels[0]->P->getVectors(); 1633 if (ev) { 1634 initial.reserve(*addr); 1635 PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, data->levels[0]->P->getDof(), ev[0], &xin)); 1636 for (n = 0; n < *addr; ++n) { 1637 PetscCall(VecDuplicate(xin, &v)); 1638 PetscCall(VecPlaceArray(xin, ev[n])); 1639 PetscCall(VecCopy(xin, v)); 1640 initial.emplace_back(v); 1641 PetscCall(VecResetArray(xin)); 1642 } 1643 PetscCall(VecDestroy(&xin)); 1644 } 1645 } 1646 } 1647 data->N -= reused; 1648 PetscCall(KSPSetOperators(data->levels[0]->ksp, A, P)); 1649 1650 PetscCall(PetscObjectTypeCompare((PetscObject)P, MATIS, &ismatis)); 1651 if (!data->is && !ismatis) { 1652 PetscErrorCode (*create)(DM, IS *, Mat *, PetscErrorCode (**)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void *), void **) = nullptr; 1653 PetscErrorCode (*usetup)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void *) = nullptr; 1654 void *uctx = nullptr; 1655 1656 /* first see if we can get the data from the DM */ 1657 PetscCall(MatGetDM(P, &dm)); 1658 if (!dm) PetscCall(MatGetDM(A, &dm)); 1659 if (!dm) PetscCall(PCGetDM(pc, &dm)); 1660 if (dm) { /* this is the hook for DMPLEX for which the auxiliary Mat is the local Neumann matrix */ 1661 PetscCall(PetscObjectQueryFunction((PetscObject)dm, "DMCreateNeumannOverlap_C", &create)); 1662 if (create) { 1663 PetscCall((*create)(dm, &uis, &uaux, &usetup, &uctx)); 1664 if (data->Neumann == PETSC_BOOL3_UNKNOWN) data->Neumann = PETSC_BOOL3_TRUE; /* set the value only if it was not already provided by the user */ 1665 } 1666 } 1667 if (!create) { 1668 if (!uis) { 1669 PetscCall(PetscObjectQuery((PetscObject)pc, "_PCHPDDM_Neumann_IS", (PetscObject *)&uis)); 1670 PetscCall(PetscObjectReference((PetscObject)uis)); 1671 } 1672 if (!uaux) { 1673 PetscCall(PetscObjectQuery((PetscObject)pc, "_PCHPDDM_Neumann_Mat", (PetscObject *)&uaux)); 1674 PetscCall(PetscObjectReference((PetscObject)uaux)); 1675 } 1676 /* look inside the Pmat instead of the PC, needed for MatSchurComplementComputeExplicitOperator() */ 1677 if (!uis) { 1678 PetscCall(PetscObjectQuery((PetscObject)P, "_PCHPDDM_Neumann_IS", (PetscObject *)&uis)); 1679 PetscCall(PetscObjectReference((PetscObject)uis)); 1680 } 1681 if (!uaux) { 1682 PetscCall(PetscObjectQuery((PetscObject)P, "_PCHPDDM_Neumann_Mat", (PetscObject *)&uaux)); 1683 PetscCall(PetscObjectReference((PetscObject)uaux)); 1684 } 1685 } 1686 PetscCall(PCHPDDMSetAuxiliaryMat(pc, uis, uaux, usetup, uctx)); 1687 PetscCall(MatDestroy(&uaux)); 1688 PetscCall(ISDestroy(&uis)); 1689 } 1690 1691 if (!ismatis) { 1692 PetscCall(PCHPDDMSetUpNeumannOverlap_Private(pc)); 1693 PetscCall(PetscOptionsGetBool(nullptr, pcpre, "-pc_hpddm_block_splitting", &block, nullptr)); 1694 PetscCall(PetscOptionsGetInt(nullptr, pcpre, "-pc_hpddm_harmonic_overlap", &overlap, nullptr)); 1695 PetscCall(PetscObjectTypeCompare((PetscObject)P, MATSCHURCOMPLEMENT, &flg)); 1696 if (data->is || flg) { 1697 if (block || overlap != -1) { 1698 PetscCall(ISDestroy(&data->is)); 1699 PetscCall(MatDestroy(&data->aux)); 1700 } else if (flg) { 1701 PCHPDDMSchurPreType type = PC_HPDDM_SCHUR_PRE_GENEO; 1702 1703 PetscCall(PetscOptionsGetEnum(nullptr, pcpre, "-pc_hpddm_schur_precondition", PCHPDDMSchurPreTypes, (PetscEnum *)&type, &flg)); 1704 if (type == PC_HPDDM_SCHUR_PRE_LEAST_SQUARES) { 1705 PetscCall(ISDestroy(&data->is)); /* destroy any previously user-set objects since they will be set automatically */ 1706 PetscCall(MatDestroy(&data->aux)); 1707 } else if (type == PC_HPDDM_SCHUR_PRE_GENEO) { 1708 PetscContainer container = nullptr; 1709 1710 PetscCall(PetscObjectQuery((PetscObject)pc, "_PCHPDDM_Schur", (PetscObject *)&container)); 1711 if (!container) { /* first call to PCSetUp() on the PC associated to the Schur complement */ 1712 PC_HPDDM *data_00; 1713 KSP ksp, inner_ksp; 1714 PC pc_00; 1715 Mat A11 = nullptr; 1716 Vec d = nullptr; 1717 const PetscInt *ranges; 1718 PetscMPIInt size; 1719 char *prefix; 1720 1721 PetscCall(MatSchurComplementGetKSP(P, &ksp)); 1722 PetscCall(KSPGetPC(ksp, &pc_00)); 1723 PetscCall(PetscObjectTypeCompare((PetscObject)pc_00, PCHPDDM, &flg)); 1724 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 : "", 1725 ((PetscObject)pc_00)->type_name, PCHPDDM); 1726 data_00 = (PC_HPDDM *)pc_00->data; 1727 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], 1728 data_00->N, data_00->N > 1 ? "s" : "", ((PetscObject)pc_00)->prefix); 1729 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() : ""); 1730 PetscCall(PetscObjectTypeCompare((PetscObject)data_00->levels[0]->pc, PCASM, &flg)); 1731 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, 1732 ((PetscObject)data_00->levels[0]->pc)->type_name, PCASM); 1733 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) */ 1734 PetscCall(MatSchurComplementGetSubMatrices(P, nullptr, nullptr, nullptr, nullptr, &A11)); 1735 PetscCall(MatGetOwnershipRanges(A11, &ranges)); 1736 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A11), &size)); 1737 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)? */ 1738 if (!flg) { 1739 if (PetscDefined(USE_DEBUG) || !data->is) { 1740 Mat A01, A10, B = nullptr, C = nullptr, *sub; 1741 1742 PetscCall(MatSchurComplementGetSubMatrices(P, &A, nullptr, &A01, &A10, nullptr)); 1743 PetscCall(PetscObjectTypeCompare((PetscObject)A10, MATTRANSPOSEVIRTUAL, &flg)); 1744 if (flg) { 1745 PetscCall(MatTransposeGetMat(A10, &C)); 1746 PetscCall(MatTranspose(C, MAT_INITIAL_MATRIX, &B)); 1747 } else { 1748 PetscCall(PetscObjectTypeCompare((PetscObject)A10, MATHERMITIANTRANSPOSEVIRTUAL, &flg)); 1749 if (flg) { 1750 PetscCall(MatHermitianTransposeGetMat(A10, &C)); 1751 PetscCall(MatHermitianTranspose(C, MAT_INITIAL_MATRIX, &B)); 1752 } 1753 } 1754 if (flg) 1755 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)); 1756 if (!B) { 1757 B = A10; 1758 PetscCall(PetscObjectReference((PetscObject)B)); 1759 } else if (!data->is) { 1760 PetscCall(PetscObjectTypeCompareAny((PetscObject)A01, &flg, MATTRANSPOSEVIRTUAL, MATHERMITIANTRANSPOSEVIRTUAL, "")); 1761 if (!flg) C = A01; 1762 else 1763 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)); 1764 } 1765 PetscCall(ISCreateStride(PETSC_COMM_SELF, B->rmap->N, 0, 1, &uis)); 1766 PetscCall(ISSetIdentity(uis)); 1767 if (!data->is) { 1768 if (C) PetscCall(PetscObjectReference((PetscObject)C)); 1769 else PetscCall(MatTranspose(B, MAT_INITIAL_MATRIX, &C)); 1770 PetscCall(ISDuplicate(data_00->is, is)); 1771 PetscCall(MatIncreaseOverlap(A, 1, is, 1)); 1772 PetscCall(MatSetOption(C, MAT_SUBMAT_SINGLEIS, PETSC_TRUE)); 1773 PetscCall(MatCreateSubMatrices(C, 1, is, &uis, MAT_INITIAL_MATRIX, &sub)); 1774 PetscCall(MatDestroy(&C)); 1775 PetscCall(MatTranspose(sub[0], MAT_INITIAL_MATRIX, &C)); 1776 PetscCall(MatDestroySubMatrices(1, &sub)); 1777 PetscCall(MatFindNonzeroRows(C, &data->is)); 1778 PetscCall(MatDestroy(&C)); 1779 PetscCall(ISDestroy(is)); 1780 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)data->is), A11->rmap->n, A11->rmap->rstart, 1, &loc)); 1781 if (PetscDefined(USE_DEBUG)) PetscCall(PCHPDDMCheckInclusion_Private(pc, data->is, loc, PETSC_FALSE)); 1782 PetscCall(ISExpand(data->is, loc, is)); 1783 PetscCall(ISDestroy(&loc)); 1784 PetscCall(ISDestroy(&data->is)); 1785 data->is = is[0]; 1786 is[0] = nullptr; 1787 } 1788 if (PetscDefined(USE_DEBUG)) { 1789 PetscCall(PCHPDDMCheckSymmetry_Private(pc, A01, A10)); 1790 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 */ 1791 PetscCall(ISDestroy(&uis)); 1792 PetscCall(ISDuplicate(data->is, &uis)); 1793 PetscCall(ISSort(uis)); 1794 PetscCall(ISComplement(uis, 0, B->rmap->N, is)); 1795 PetscCall(MatDuplicate(sub[0], MAT_COPY_VALUES, &C)); 1796 PetscCall(MatZeroRowsIS(C, is[0], 0.0, nullptr, nullptr)); 1797 PetscCall(ISDestroy(is)); 1798 PetscCall(MatMultEqual(sub[0], C, 20, &flg)); 1799 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 */ 1800 PetscCall(MatDestroy(&C)); 1801 PetscCall(MatDestroySubMatrices(1, &sub)); 1802 } 1803 PetscCall(ISDestroy(&uis)); 1804 PetscCall(MatDestroy(&B)); 1805 } 1806 flg = PETSC_FALSE; 1807 if (!data->aux) { 1808 Mat D; 1809 1810 PetscCall(MatCreateVecs(A11, &d, nullptr)); 1811 PetscCall(MatGetDiagonal(A11, d)); 1812 PetscCall(PetscObjectTypeCompareAny((PetscObject)A11, &flg, MATDIAGONAL, MATCONSTANTDIAGONAL, "")); 1813 if (!flg) { 1814 PetscCall(MatCreateDiagonal(d, &D)); 1815 PetscCall(MatMultEqual(A11, D, 20, &flg)); 1816 PetscCall(MatDestroy(&D)); 1817 } 1818 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")); 1819 } 1820 if (data->Neumann != PETSC_BOOL3_TRUE && !flg && A11) { 1821 PetscReal norm; 1822 1823 PetscCall(MatNorm(A11, NORM_INFINITY, &norm)); 1824 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 : ""); 1825 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")); 1826 PetscCall(MatDestroy(&data->aux)); 1827 flg = PETSC_TRUE; 1828 } 1829 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 */ 1830 PetscSF scatter; 1831 const PetscScalar *read; 1832 PetscScalar *write, *diagonal = nullptr; 1833 1834 PetscCall(MatDestroy(&data->aux)); 1835 PetscCall(ISGetLocalSize(data->is, &n)); 1836 PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)P), n, PETSC_DECIDE, &xin)); 1837 PetscCall(VecDuplicate(xin, &v)); 1838 PetscCall(VecScatterCreate(xin, data->is, v, nullptr, &scatter)); 1839 PetscCall(VecSet(v, 1.0)); 1840 PetscCall(VecSet(xin, 1.0)); 1841 PetscCall(VecScatterBegin(scatter, v, xin, ADD_VALUES, SCATTER_REVERSE)); 1842 PetscCall(VecScatterEnd(scatter, v, xin, ADD_VALUES, SCATTER_REVERSE)); /* v has the multiplicity of all unknowns on the overlap */ 1843 PetscCall(PetscSFDestroy(&scatter)); 1844 if (d) { 1845 PetscCall(VecScatterCreate(d, data->is, v, nullptr, &scatter)); 1846 PetscCall(VecScatterBegin(scatter, d, v, INSERT_VALUES, SCATTER_FORWARD)); 1847 PetscCall(VecScatterEnd(scatter, d, v, INSERT_VALUES, SCATTER_FORWARD)); 1848 PetscCall(PetscSFDestroy(&scatter)); 1849 PetscCall(VecDestroy(&d)); 1850 PetscCall(PetscMalloc1(n, &diagonal)); 1851 PetscCall(VecGetArrayRead(v, &read)); 1852 PetscCallCXX(std::copy_n(read, n, diagonal)); 1853 PetscCall(VecRestoreArrayRead(v, &read)); 1854 } 1855 PetscCall(VecDestroy(&v)); 1856 PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &v)); 1857 PetscCall(VecGetArrayRead(xin, &read)); 1858 PetscCall(VecGetArrayWrite(v, &write)); 1859 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]; 1860 PetscCall(PetscFree(diagonal)); 1861 PetscCall(VecRestoreArrayRead(xin, &read)); 1862 PetscCall(VecRestoreArrayWrite(v, &write)); 1863 PetscCall(VecDestroy(&xin)); 1864 PetscCall(MatCreateDiagonal(v, &data->aux)); 1865 PetscCall(VecDestroy(&v)); 1866 } 1867 uis = data->is; 1868 uaux = data->aux; 1869 PetscCall(PetscObjectReference((PetscObject)uis)); 1870 PetscCall(PetscObjectReference((PetscObject)uaux)); 1871 PetscCall(PetscStrallocpy(pcpre, &prefix)); 1872 PetscCall(PCSetOptionsPrefix(pc, nullptr)); 1873 PetscCall(PCSetType(pc, PCKSP)); /* replace the PC associated to the Schur complement by PCKSP */ 1874 PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &inner_ksp)); /* new KSP that will be attached to the previously set PC */ 1875 PetscCall(PetscObjectGetTabLevel((PetscObject)pc, &n)); 1876 PetscCall(PetscObjectSetTabLevel((PetscObject)inner_ksp, n + 2)); 1877 PetscCall(KSPSetOperators(inner_ksp, pc->mat, pc->pmat)); 1878 PetscCall(KSPSetOptionsPrefix(inner_ksp, std::string(std::string(prefix) + "pc_hpddm_").c_str())); 1879 PetscCall(KSPSetSkipPCSetFromOptions(inner_ksp, PETSC_TRUE)); 1880 PetscCall(KSPSetFromOptions(inner_ksp)); 1881 PetscCall(KSPGetPC(inner_ksp, &inner)); 1882 PetscCall(PCSetOptionsPrefix(inner, nullptr)); 1883 PetscCall(PCSetType(inner, PCNONE)); /* no preconditioner since the action of M^-1 A or A M^-1 will be computed by the Amat */ 1884 PetscCall(PCKSPSetKSP(pc, inner_ksp)); 1885 std::get<0>(*ctx)[0] = pc_00; /* for coarse correction on the primal (e.g., velocity) space */ 1886 PetscCall(PCCreate(PetscObjectComm((PetscObject)pc), &std::get<0>(*ctx)[1])); 1887 PetscCall(PCSetOptionsPrefix(pc, prefix)); /* both PC share the same prefix so that the outer PC can be reset with PCSetFromOptions() */ 1888 PetscCall(PCSetOptionsPrefix(std::get<0>(*ctx)[1], prefix)); 1889 PetscCall(PetscFree(prefix)); 1890 PetscCall(PCSetOperators(std::get<0>(*ctx)[1], pc->mat, pc->pmat)); 1891 PetscCall(PCSetType(std::get<0>(*ctx)[1], PCHPDDM)); 1892 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 */ 1893 if (flg) static_cast<PC_HPDDM *>(std::get<0>(*ctx)[1]->data)->Neumann = PETSC_BOOL3_TRUE; 1894 PetscCall(PCSetFromOptions(std::get<0>(*ctx)[1])); 1895 PetscCall(PetscObjectDereference((PetscObject)uis)); 1896 PetscCall(PetscObjectDereference((PetscObject)uaux)); 1897 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 */ 1898 PetscCall(MatShellSetOperation(S, MATOP_MULT, (void (*)(void))MatMult_SchurCorrection)); 1899 PetscCall(MatShellSetOperation(S, MATOP_VIEW, (void (*)(void))MatView_SchurCorrection)); 1900 PetscCall(MatShellSetOperation(S, MATOP_DESTROY, (void (*)(void))MatDestroy_SchurCorrection)); 1901 PetscCall(KSPGetPCSide(inner_ksp, &(std::get<2>(*ctx)))); 1902 if (std::get<2>(*ctx) == PC_LEFT || std::get<2>(*ctx) == PC_SIDE_DEFAULT) { 1903 PetscCall(KSPSetPreSolve(inner_ksp, KSPPreSolve_SchurCorrection, ctx)); 1904 } else { /* no support for PC_SYMMETRIC */ 1905 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]); 1906 } 1907 PetscCall(KSPSetPostSolve(inner_ksp, KSPPostSolve_SchurCorrection, ctx)); 1908 PetscCall(PetscObjectContainerCompose((PetscObject)std::get<0>(*ctx)[1], "_PCHPDDM_Schur", ctx, nullptr)); 1909 PetscCall(PCSetUp(std::get<0>(*ctx)[1])); 1910 PetscCall(KSPSetOperators(inner_ksp, S, S)); 1911 PetscCall(MatCreateVecs(std::get<1>(*ctx)[0], std::get<3>(*ctx), std::get<3>(*ctx) + 1)); 1912 PetscCall(VecDuplicate(std::get<3>(*ctx)[0], std::get<3>(*ctx) + 2)); 1913 PetscCall(PetscObjectDereference((PetscObject)inner_ksp)); 1914 PetscCall(PetscObjectDereference((PetscObject)S)); 1915 } else { 1916 std::get<0>(*ctx)[0] = pc_00; 1917 PetscCall(PetscObjectContainerCompose((PetscObject)pc, "_PCHPDDM_Schur", ctx, nullptr)); 1918 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 */ 1919 PetscCall(MatGetDiagonalBlock(A11, &data->aux)); 1920 PetscCall(PetscObjectReference((PetscObject)data->aux)); 1921 PetscCall(PCSetUp(pc)); 1922 } 1923 for (std::vector<Vec>::iterator it = initial.begin(); it != initial.end(); ++it) PetscCall(VecDestroy(&*it)); 1924 PetscFunctionReturn(PETSC_SUCCESS); 1925 } else { /* second call to PCSetUp() on the PC associated to the Schur complement, retrieve previously set context */ 1926 PetscCall(PetscContainerGetPointer(container, (void **)&ctx)); 1927 } 1928 } 1929 } 1930 } 1931 if (!data->is && data->N > 1) { 1932 char type[256] = {}; /* same size as in src/ksp/pc/interface/pcset.c */ 1933 1934 PetscCall(PetscObjectTypeCompareAny((PetscObject)P, &flg, MATNORMAL, MATNORMALHERMITIAN, "")); 1935 if (flg || (A->rmap->N != A->cmap->N && P->rmap->N == P->cmap->N && P->rmap->N == A->cmap->N)) { 1936 Mat B; 1937 1938 PetscCall(PCHPDDMSetAuxiliaryMatNormal_Private(pc, A, P, &B, pcpre)); 1939 if (data->correction == PC_HPDDM_COARSE_CORRECTION_DEFLATED) data->correction = PC_HPDDM_COARSE_CORRECTION_BALANCED; 1940 PetscCall(MatDestroy(&B)); 1941 } else { 1942 PetscCall(PetscObjectTypeCompare((PetscObject)P, MATSCHURCOMPLEMENT, &flg)); 1943 if (flg) { 1944 Mat A00, P00, A01, A10, A11, B, N; 1945 PCHPDDMSchurPreType type = PC_HPDDM_SCHUR_PRE_LEAST_SQUARES; 1946 1947 PetscCall(MatSchurComplementGetSubMatrices(P, &A00, &P00, &A01, &A10, &A11)); 1948 PetscCall(PetscOptionsGetEnum(nullptr, pcpre, "-pc_hpddm_schur_precondition", PCHPDDMSchurPreTypes, (PetscEnum *)&type, &flg)); 1949 if (type == PC_HPDDM_SCHUR_PRE_LEAST_SQUARES) { 1950 Mat B01; 1951 Vec diagonal = nullptr; 1952 const PetscScalar *array; 1953 MatSchurComplementAinvType type; 1954 1955 PetscCall(PCHPDDMCheckSymmetry_Private(pc, A01, A10, &B01)); 1956 if (A11) { 1957 PetscCall(MatCreateVecs(A11, &diagonal, nullptr)); 1958 PetscCall(MatGetDiagonal(A11, diagonal)); 1959 } 1960 PetscCall(MatCreateVecs(P00, &v, nullptr)); 1961 PetscCall(MatSchurComplementGetAinvType(P, &type)); 1962 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]); 1963 if (type == MAT_SCHUR_COMPLEMENT_AINV_LUMP) { 1964 PetscCall(MatGetRowSum(P00, v)); 1965 if (A00 == P00) PetscCall(PetscObjectReference((PetscObject)A00)); 1966 PetscCall(MatDestroy(&P00)); 1967 PetscCall(VecGetArrayRead(v, &array)); 1968 PetscCall(MatCreateAIJ(PetscObjectComm((PetscObject)A00), A00->rmap->n, A00->cmap->n, A00->rmap->N, A00->cmap->N, 1, nullptr, 0, nullptr, &P00)); 1969 PetscCall(MatSetOption(P00, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE)); 1970 for (n = A00->rmap->rstart; n < A00->rmap->rend; ++n) PetscCall(MatSetValue(P00, n, n, array[n - A00->rmap->rstart], INSERT_VALUES)); 1971 PetscCall(MatAssemblyBegin(P00, MAT_FINAL_ASSEMBLY)); 1972 PetscCall(MatAssemblyEnd(P00, MAT_FINAL_ASSEMBLY)); 1973 PetscCall(VecRestoreArrayRead(v, &array)); 1974 PetscCall(MatSchurComplementUpdateSubMatrices(P, A00, P00, A01, A10, A11)); /* replace P00 by diag(sum of each row of P00) */ 1975 PetscCall(MatDestroy(&P00)); 1976 } else PetscCall(MatGetDiagonal(P00, v)); 1977 PetscCall(VecReciprocal(v)); /* inv(diag(P00)) */ 1978 PetscCall(VecSqrtAbs(v)); /* sqrt(inv(diag(P00))) */ 1979 PetscCall(MatDuplicate(A01, MAT_COPY_VALUES, &B)); 1980 PetscCall(MatDiagonalScale(B, v, nullptr)); 1981 if (B01) PetscCall(MatDiagonalScale(B01, v, nullptr)); 1982 PetscCall(VecDestroy(&v)); 1983 PetscCall(MatCreateNormalHermitian(B, &N)); 1984 PetscCall(PCHPDDMSetAuxiliaryMatNormal_Private(pc, B, N, &P, pcpre, &diagonal, B01)); 1985 PetscCall(PetscObjectTypeCompare((PetscObject)data->aux, MATSEQAIJ, &flg)); 1986 if (!flg) { 1987 PetscCall(MatDestroy(&P)); 1988 P = N; 1989 PetscCall(PetscObjectReference((PetscObject)P)); 1990 } 1991 if (diagonal) { 1992 PetscCall(MatDiagonalSet(P, diagonal, ADD_VALUES)); 1993 PetscCall(PCSetOperators(pc, P, P)); /* replace P by A01^T inv(diag(P00)) A01 - diag(P11) */ 1994 PetscCall(VecDestroy(&diagonal)); 1995 } else PetscCall(PCSetOperators(pc, B01 ? P : N, P)); /* replace P by A01^T inv(diag(P00)) A01 */ 1996 pc->ops->postsolve = PCPostSolve_SchurPreLeastSquares; /* PCFIELDSPLIT expect a KSP for (P11 - A10 inv(diag(P00)) A01) */ 1997 PetscCall(MatDestroy(&N)); /* but a PC for (A10 inv(diag(P00)) A10 - P11) is setup instead */ 1998 PetscCall(MatDestroy(&P)); /* so the sign of the solution must be flipped */ 1999 PetscCall(MatDestroy(&B)); 2000 } else 2001 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 : ""); 2002 for (std::vector<Vec>::iterator it = initial.begin(); it != initial.end(); ++it) PetscCall(VecDestroy(&*it)); 2003 PetscFunctionReturn(PETSC_SUCCESS); 2004 } else { 2005 PetscCall(PetscOptionsGetString(nullptr, pcpre, "-pc_hpddm_levels_1_st_pc_type", type, sizeof(type), nullptr)); 2006 PetscCall(PetscStrcmp(type, PCMAT, &algebraic)); 2007 PetscCheck(!algebraic || !block, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-pc_hpddm_levels_1_st_pc_type mat and -pc_hpddm_block_splitting"); 2008 if (overlap != -1) { 2009 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"); 2010 PetscCheck(overlap >= 1, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_WRONG, "-pc_hpddm_harmonic_overlap %" PetscInt_FMT " < 1", overlap); 2011 } 2012 if (block || overlap != -1) algebraic = PETSC_TRUE; 2013 if (algebraic) { 2014 PetscCall(ISCreateStride(PETSC_COMM_SELF, P->rmap->n, P->rmap->rstart, 1, &data->is)); 2015 PetscCall(MatIncreaseOverlap(P, 1, &data->is, 1)); 2016 PetscCall(ISSort(data->is)); 2017 } else 2018 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 : "")); 2019 } 2020 } 2021 } 2022 } 2023 #if PetscDefined(USE_DEBUG) 2024 if (data->is) PetscCall(ISDuplicate(data->is, &dis)); 2025 if (data->aux) PetscCall(MatDuplicate(data->aux, MAT_COPY_VALUES, &daux)); 2026 #endif 2027 if (data->is || (ismatis && data->N > 1)) { 2028 if (ismatis) { 2029 std::initializer_list<std::string> list = {MATSEQBAIJ, MATSEQSBAIJ}; 2030 PetscCall(MatISGetLocalMat(P, &N)); 2031 std::initializer_list<std::string>::const_iterator it = std::find(list.begin(), list.end(), ((PetscObject)N)->type_name); 2032 PetscCall(MatISRestoreLocalMat(P, &N)); 2033 switch (std::distance(list.begin(), it)) { 2034 case 0: 2035 PetscCall(MatConvert(P, MATMPIBAIJ, MAT_INITIAL_MATRIX, &C)); 2036 break; 2037 case 1: 2038 /* MatCreateSubMatrices() does not work with MATSBAIJ and unsorted ISes, so convert to MPIBAIJ */ 2039 PetscCall(MatConvert(P, MATMPIBAIJ, MAT_INITIAL_MATRIX, &C)); 2040 PetscCall(MatSetOption(C, MAT_SYMMETRIC, PETSC_TRUE)); 2041 break; 2042 default: 2043 PetscCall(MatConvert(P, MATMPIAIJ, MAT_INITIAL_MATRIX, &C)); 2044 } 2045 PetscCall(MatISGetLocalToGlobalMapping(P, &l2g, nullptr)); 2046 PetscCall(PetscObjectReference((PetscObject)P)); 2047 PetscCall(KSPSetOperators(data->levels[0]->ksp, A, C)); 2048 std::swap(C, P); 2049 PetscCall(ISLocalToGlobalMappingGetSize(l2g, &n)); 2050 PetscCall(ISCreateStride(PETSC_COMM_SELF, n, 0, 1, &loc)); 2051 PetscCall(ISLocalToGlobalMappingApplyIS(l2g, loc, &is[0])); 2052 PetscCall(ISDestroy(&loc)); 2053 /* the auxiliary Mat is _not_ the local Neumann matrix */ 2054 /* it is the local Neumann matrix augmented (with zeros) through MatIncreaseOverlap() */ 2055 data->Neumann = PETSC_BOOL3_FALSE; 2056 structure = SAME_NONZERO_PATTERN; 2057 } else { 2058 is[0] = data->is; 2059 if (algebraic || ctx) subdomains = PETSC_TRUE; 2060 PetscCall(PetscOptionsGetBool(nullptr, pcpre, "-pc_hpddm_define_subdomains", &subdomains, nullptr)); 2061 if (ctx) PetscCheck(subdomains, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-%spc_hpddm_schur_precondition geneo and -%spc_hpddm_define_subdomains false", pcpre, pcpre); 2062 if (PetscBool3ToBool(data->Neumann)) { 2063 PetscCheck(!block, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-pc_hpddm_block_splitting and -pc_hpddm_has_neumann"); 2064 PetscCheck(overlap == -1, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-pc_hpddm_harmonic_overlap %" PetscInt_FMT " and -pc_hpddm_has_neumann", overlap); 2065 PetscCheck(!algebraic, PetscObjectComm((PetscObject)P), PETSC_ERR_ARG_INCOMP, "-pc_hpddm_levels_1_st_pc_type mat and -pc_hpddm_has_neumann"); 2066 } 2067 if (PetscBool3ToBool(data->Neumann) || block) structure = SAME_NONZERO_PATTERN; 2068 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)data->is), P->rmap->n, P->rmap->rstart, 1, &loc)); 2069 } 2070 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_levels_1_", pcpre ? pcpre : "")); 2071 PetscCall(PetscOptionsGetEnum(nullptr, prefix, "-st_matstructure", MatStructures, (PetscEnum *)&structure, &flg)); /* if not user-provided, force its value when possible */ 2072 if (!flg && structure == SAME_NONZERO_PATTERN) { /* cannot call STSetMatStructure() yet, insert the appropriate option in the database, parsed by STSetFromOptions() */ 2073 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "-%spc_hpddm_levels_1_st_matstructure", pcpre ? pcpre : "")); 2074 PetscCall(PetscOptionsSetValue(nullptr, prefix, MatStructures[structure])); 2075 } 2076 flg = PETSC_FALSE; 2077 if (data->share) { 2078 data->share = PETSC_FALSE; /* will be reset to PETSC_TRUE if none of the conditions below are true */ 2079 if (!subdomains) PetscCall(PetscInfo(pc, "Cannot share subdomain KSP between SLEPc and PETSc since -%spc_hpddm_define_subdomains is not true\n", pcpre ? pcpre : "")); 2080 else if (data->deflation) PetscCall(PetscInfo(pc, "Nothing to share since PCHPDDMSetDeflationMat() has been called\n")); 2081 else if (ismatis) PetscCall(PetscInfo(pc, "Cannot share subdomain KSP between SLEPc and PETSc with a Pmat of type MATIS\n")); 2082 else if (!algebraic && structure != SAME_NONZERO_PATTERN) 2083 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])); 2084 else data->share = PETSC_TRUE; 2085 } 2086 if (!ismatis) { 2087 if (data->share || (!PetscBool3ToBool(data->Neumann) && subdomains)) PetscCall(ISDuplicate(is[0], &unsorted)); 2088 else unsorted = is[0]; 2089 } 2090 if ((ctx || data->N > 1) && (data->aux || ismatis || algebraic)) { 2091 PetscCheck(loadedSym, PETSC_COMM_SELF, PETSC_ERR_PLIB, "HPDDM library not loaded, cannot use more than one level"); 2092 PetscCall(MatSetOption(P, MAT_SUBMAT_SINGLEIS, PETSC_TRUE)); 2093 if (ismatis) { 2094 /* needed by HPDDM (currently) so that the partition of unity is 0 on subdomain interfaces */ 2095 PetscCall(MatIncreaseOverlap(P, 1, is, 1)); 2096 PetscCall(ISDestroy(&data->is)); 2097 data->is = is[0]; 2098 } else { 2099 if (PetscDefined(USE_DEBUG)) PetscCall(PCHPDDMCheckInclusion_Private(pc, data->is, loc, PETSC_TRUE)); 2100 if (!ctx && overlap == -1) PetscCall(PetscObjectComposeFunction((PetscObject)pc->pmat, "PCHPDDMAlgebraicAuxiliaryMat_Private_C", PCHPDDMAlgebraicAuxiliaryMat_Private)); 2101 if (!PetscBool3ToBool(data->Neumann) && (!algebraic || overlap != -1)) { 2102 PetscCall(PetscObjectTypeCompare((PetscObject)P, MATMPISBAIJ, &flg)); 2103 if (flg) { 2104 /* maybe better to ISSort(is[0]), MatCreateSubMatrices(), and then MatPermute() */ 2105 /* but there is no MatPermute_SeqSBAIJ(), so as before, just use MATMPIBAIJ */ 2106 PetscCall(MatConvert(P, MATMPIBAIJ, MAT_INITIAL_MATRIX, &uaux)); 2107 flg = PETSC_FALSE; 2108 } 2109 } 2110 } 2111 if (algebraic && overlap == -1) { 2112 PetscUseMethod(pc->pmat, "PCHPDDMAlgebraicAuxiliaryMat_Private_C", (Mat, IS *, Mat *[], PetscBool), (P, is, &sub, block)); 2113 if (block) { 2114 PetscCall(PetscObjectQuery((PetscObject)sub[0], "_PCHPDDM_Neumann_Mat", (PetscObject *)&data->aux)); 2115 PetscCall(PetscObjectCompose((PetscObject)sub[0], "_PCHPDDM_Neumann_Mat", nullptr)); 2116 } 2117 } else if (!uaux || overlap != -1) { 2118 if (!ctx) { 2119 if (PetscBool3ToBool(data->Neumann)) sub = &data->aux; 2120 else { 2121 PetscBool flg; 2122 if (overlap != -1) { 2123 Harmonic h; 2124 Mat A0, *a; /* with an SVD: [ A_00 A_01 ] */ 2125 IS ov[2], rows, cols, stride; /* [ A_10 A_11 A_12 ] */ 2126 const PetscInt *i[2], bs = P->cmap->bs; /* with a GEVP: [ A_00 A_01 ] */ 2127 PetscInt n[2]; /* [ A_10 A_11 A_12 ] */ 2128 std::vector<PetscInt> v[2]; /* [ A_21 A_22 ] */ 2129 2130 PetscCall(ISDuplicate(data->is, ov)); 2131 if (overlap > 1) PetscCall(MatIncreaseOverlap(P, 1, ov, overlap - 1)); 2132 PetscCall(ISDuplicate(ov[0], ov + 1)); 2133 PetscCall(MatIncreaseOverlap(P, 1, ov + 1, 1)); 2134 PetscCall(PetscNew(&h)); 2135 h->ksp = nullptr; 2136 PetscCall(PetscCalloc1(2, &h->A)); 2137 PetscCall(PetscOptionsHasName(nullptr, prefix, "-svd_nsv", &flg)); 2138 if (!flg) PetscCall(PetscOptionsHasName(nullptr, prefix, "-svd_relative_threshold", &flg)); 2139 PetscCall(ISSort(ov[0])); 2140 if (!flg) PetscCall(ISSort(ov[1])); 2141 PetscCall(PetscCalloc1(5, &h->is)); 2142 PetscCall(MatCreateSubMatrices(uaux ? uaux : P, 1, ov + !flg, ov + 1, MAT_INITIAL_MATRIX, &a)); /* submatrix from above, either square (!flg) or rectangular (flg) */ 2143 for (PetscInt j = 0; j < 2; ++j) { 2144 PetscCall(ISGetIndices(ov[j], i + j)); 2145 PetscCall(ISGetLocalSize(ov[j], n + j)); 2146 } 2147 v[1].reserve((n[1] - n[0]) / bs); 2148 for (PetscInt j = 0; j < n[1]; j += bs) { /* indices of the (2,2) block */ 2149 PetscInt location; 2150 PetscCall(ISLocate(ov[0], i[1][j], &location)); 2151 if (location < 0) v[1].emplace_back(j / bs); 2152 } 2153 if (!flg) { 2154 h->A[1] = a[0]; 2155 PetscCall(PetscObjectReference((PetscObject)h->A[1])); 2156 v[0].reserve((n[0] - P->rmap->n) / bs); 2157 for (PetscInt j = 0; j < n[1]; j += bs) { /* row indices of the (1,2) block */ 2158 PetscInt location; 2159 PetscCall(ISLocate(loc, i[1][j], &location)); 2160 if (location < 0) { 2161 PetscCall(ISLocate(ov[0], i[1][j], &location)); 2162 if (location >= 0) v[0].emplace_back(j / bs); 2163 } 2164 } 2165 PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, v[0].size(), v[0].data(), PETSC_USE_POINTER, &rows)); 2166 PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, v[1].size(), v[1].data(), PETSC_COPY_VALUES, h->is + 4)); 2167 PetscCall(MatCreateSubMatrix(a[0], rows, h->is[4], MAT_INITIAL_MATRIX, h->A)); /* A_12 submatrix from above */ 2168 PetscCall(ISDestroy(&rows)); 2169 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 */ 2170 PetscCall(ISEmbed(ov[0], ov[1], PETSC_TRUE, &rows)); 2171 PetscCall(MatCreateSubMatrix(a[0], rows, cols = rows, MAT_INITIAL_MATRIX, &A0)); /* [ A_00 A_01 ; A_10 A_11 ] submatrix from above */ 2172 PetscCall(ISDestroy(&rows)); 2173 v[0].clear(); 2174 PetscCall(ISEmbed(loc, ov[1], PETSC_TRUE, h->is + 3)); 2175 PetscCall(ISEmbed(data->is, ov[1], PETSC_TRUE, h->is + 2)); 2176 } 2177 v[0].reserve((n[0] - P->rmap->n) / bs); 2178 for (PetscInt j = 0; j < n[0]; j += bs) { 2179 PetscInt location; 2180 PetscCall(ISLocate(loc, i[0][j], &location)); 2181 if (location < 0) v[0].emplace_back(j / bs); 2182 } 2183 PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, v[0].size(), v[0].data(), PETSC_USE_POINTER, &rows)); 2184 for (PetscInt j = 0; j < 2; ++j) PetscCall(ISRestoreIndices(ov[j], i + j)); 2185 if (flg) { 2186 IS is; 2187 PetscCall(ISCreateStride(PETSC_COMM_SELF, a[0]->rmap->n, 0, 1, &is)); 2188 PetscCall(ISEmbed(ov[0], ov[1], PETSC_TRUE, &cols)); 2189 PetscCall(MatCreateSubMatrix(a[0], is, cols, MAT_INITIAL_MATRIX, &A0)); /* [ A_00 A_01 ; A_10 A_11 ] submatrix from above */ 2190 PetscCall(ISDestroy(&cols)); 2191 PetscCall(ISDestroy(&is)); 2192 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 */ 2193 PetscCall(ISEmbed(loc, data->is, PETSC_TRUE, h->is + 2)); 2194 PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, v[1].size(), v[1].data(), PETSC_USE_POINTER, &cols)); 2195 PetscCall(MatCreateSubMatrix(a[0], rows, cols, MAT_INITIAL_MATRIX, h->A)); /* A_12 submatrix from above */ 2196 PetscCall(ISDestroy(&cols)); 2197 } 2198 PetscCall(ISCreateStride(PETSC_COMM_SELF, A0->rmap->n, 0, 1, &stride)); 2199 PetscCall(ISEmbed(rows, stride, PETSC_TRUE, h->is)); 2200 PetscCall(ISDestroy(&stride)); 2201 PetscCall(ISDestroy(&rows)); 2202 PetscCall(ISEmbed(loc, ov[0], PETSC_TRUE, h->is + 1)); 2203 if (subdomains) { 2204 if (!data->levels[0]->pc) { 2205 PetscCall(PCCreate(PetscObjectComm((PetscObject)pc), &data->levels[0]->pc)); 2206 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_levels_1_", pcpre ? pcpre : "")); 2207 PetscCall(PCSetOptionsPrefix(data->levels[0]->pc, prefix)); 2208 PetscCall(PCSetOperators(data->levels[0]->pc, A, P)); 2209 } 2210 PetscCall(PCSetType(data->levels[0]->pc, PCASM)); 2211 if (!data->levels[0]->pc->setupcalled) PetscCall(PCASMSetLocalSubdomains(data->levels[0]->pc, 1, ov + !flg, &loc)); 2212 PetscCall(PCHPDDMCommunicationAvoidingPCASM_Private(data->levels[0]->pc, flg ? A0 : a[0], PETSC_TRUE)); 2213 if (!flg) ++overlap; 2214 if (data->share) { 2215 PetscInt n = -1; 2216 PetscTryMethod(data->levels[0]->pc, "PCASMGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (data->levels[0]->pc, &n, nullptr, &ksp)); 2217 PetscCheck(n == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of subdomain solver %" PetscInt_FMT " != 1", n); 2218 if (flg) { 2219 h->ksp = ksp[0]; 2220 PetscCall(PetscObjectReference((PetscObject)h->ksp)); 2221 } 2222 } 2223 } 2224 if (!h->ksp) { 2225 PetscBool share = data->share; 2226 PetscCall(KSPCreate(PETSC_COMM_SELF, &h->ksp)); 2227 PetscCall(KSPSetType(h->ksp, KSPPREONLY)); 2228 PetscCall(KSPSetOperators(h->ksp, A0, A0)); 2229 do { 2230 if (!data->share) { 2231 share = PETSC_FALSE; 2232 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_levels_1_%s", pcpre ? pcpre : "", flg ? "svd_" : "eps_")); 2233 PetscCall(KSPSetOptionsPrefix(h->ksp, prefix)); 2234 PetscCall(KSPSetFromOptions(h->ksp)); 2235 } else { 2236 MatSolverType type; 2237 PetscCall(KSPGetPC(ksp[0], &pc)); 2238 PetscCall(PetscObjectTypeCompareAny((PetscObject)pc, &data->share, PCLU, PCCHOLESKY, "")); 2239 if (data->share) { 2240 PetscCall(PCFactorGetMatSolverType(pc, &type)); 2241 if (!type) { 2242 if (PetscDefined(HAVE_MUMPS)) PetscCall(PCFactorSetMatSolverType(pc, MATSOLVERMUMPS)); 2243 else if (PetscDefined(HAVE_MKL_PARDISO)) PetscCall(PCFactorSetMatSolverType(pc, MATSOLVERMKL_PARDISO)); 2244 else data->share = PETSC_FALSE; 2245 if (data->share) PetscCall(PCSetFromOptions(pc)); 2246 } else { 2247 PetscCall(PetscStrcmp(type, MATSOLVERMUMPS, &data->share)); 2248 if (!data->share) PetscCall(PetscStrcmp(type, MATSOLVERMKL_PARDISO, &data->share)); 2249 } 2250 if (data->share) { 2251 std::tuple<KSP, IS, Vec[2]> *p; 2252 PetscCall(PCFactorGetMatrix(pc, &A)); 2253 PetscCall(MatFactorSetSchurIS(A, h->is[4])); 2254 PetscCall(KSPSetUp(ksp[0])); 2255 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_levels_1_eps_shell_", pcpre ? pcpre : "")); 2256 PetscCall(KSPSetOptionsPrefix(h->ksp, prefix)); 2257 PetscCall(KSPSetFromOptions(h->ksp)); 2258 PetscCall(KSPGetPC(h->ksp, &pc)); 2259 PetscCall(PCSetType(pc, PCSHELL)); 2260 PetscCall(PetscNew(&p)); 2261 std::get<0>(*p) = ksp[0]; 2262 PetscCall(ISEmbed(ov[0], ov[1], PETSC_TRUE, &std::get<1>(*p))); 2263 PetscCall(MatCreateVecs(A, std::get<2>(*p), std::get<2>(*p) + 1)); 2264 PetscCall(PCShellSetContext(pc, p)); 2265 PetscCall(PCShellSetApply(pc, PCApply_Schur)); 2266 PetscCall(PCShellSetApplyTranspose(pc, PCApply_Schur<Vec, true>)); 2267 PetscCall(PCShellSetMatApply(pc, PCApply_Schur<Mat>)); 2268 PetscCall(PCShellSetDestroy(pc, PCDestroy_Schur)); 2269 } 2270 } 2271 if (!data->share) PetscCall(PetscInfo(pc, "Cannot share subdomain KSP between SLEPc and PETSc since neither MUMPS nor MKL PARDISO is used\n")); 2272 } 2273 } 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 */ 2274 } 2275 PetscCall(ISDestroy(ov)); 2276 PetscCall(ISDestroy(ov + 1)); 2277 if (overlap == 1 && subdomains && flg) { 2278 *subA = A0; 2279 sub = subA; 2280 if (uaux) PetscCall(MatDestroy(&uaux)); 2281 } else PetscCall(MatDestroy(&A0)); 2282 PetscCall(MatCreateShell(PETSC_COMM_SELF, P->rmap->n, n[1] - n[0], P->rmap->n, n[1] - n[0], h, &data->aux)); 2283 PetscCall(KSPSetErrorIfNotConverged(h->ksp, PETSC_TRUE)); /* bail out as early as possible to avoid (apparently) unrelated error messages */ 2284 PetscCall(MatCreateVecs(h->ksp->pc->pmat, &h->v, nullptr)); 2285 PetscCall(MatShellSetOperation(data->aux, MATOP_MULT, (void (*)(void))MatMult_Harmonic)); 2286 PetscCall(MatShellSetOperation(data->aux, MATOP_MULT_TRANSPOSE, (void (*)(void))MatMultTranspose_Harmonic)); 2287 PetscCall(MatShellSetMatProductOperation(data->aux, MATPRODUCT_AB, nullptr, MatProduct_AB_Harmonic, nullptr, MATDENSE, MATDENSE)); 2288 PetscCall(MatShellSetOperation(data->aux, MATOP_DESTROY, (void (*)(void))MatDestroy_Harmonic)); 2289 PetscCall(MatDestroySubMatrices(1, &a)); 2290 } 2291 if (overlap != 1 || !subdomains) { 2292 PetscCall(MatCreateSubMatrices(uaux ? uaux : P, 1, is, is, MAT_INITIAL_MATRIX, &sub)); 2293 if (ismatis) { 2294 PetscCall(MatISGetLocalMat(C, &N)); 2295 PetscCall(PetscObjectTypeCompare((PetscObject)N, MATSEQSBAIJ, &flg)); 2296 if (flg) PetscCall(MatConvert(sub[0], MATSEQSBAIJ, MAT_INPLACE_MATRIX, sub)); 2297 PetscCall(MatISRestoreLocalMat(C, &N)); 2298 } 2299 } 2300 if (uaux) { 2301 PetscCall(MatDestroy(&uaux)); 2302 PetscCall(MatConvert(sub[0], MATSEQSBAIJ, MAT_INPLACE_MATRIX, sub)); 2303 } 2304 } 2305 } 2306 } else if (!ctx) { 2307 PetscCall(MatCreateSubMatrices(uaux, 1, is, is, MAT_INITIAL_MATRIX, &sub)); 2308 PetscCall(MatDestroy(&uaux)); 2309 PetscCall(MatConvert(sub[0], MATSEQSBAIJ, MAT_INPLACE_MATRIX, sub)); 2310 } 2311 if (data->N > 1) { 2312 /* Vec holding the partition of unity */ 2313 if (!data->levels[0]->D) { 2314 PetscCall(ISGetLocalSize(data->is, &n)); 2315 PetscCall(VecCreateMPI(PETSC_COMM_SELF, n, PETSC_DETERMINE, &data->levels[0]->D)); 2316 } 2317 if (data->share && overlap == -1) { 2318 Mat D; 2319 IS perm = nullptr; 2320 PetscInt size = -1; 2321 2322 if (!data->levels[0]->pc) { 2323 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_hpddm_levels_1_", pcpre ? pcpre : "")); 2324 PetscCall(PCCreate(PetscObjectComm((PetscObject)pc), &data->levels[0]->pc)); 2325 PetscCall(PCSetOptionsPrefix(data->levels[0]->pc, prefix)); 2326 PetscCall(PCSetOperators(data->levels[0]->pc, A, P)); 2327 } 2328 PetscCall(PCSetType(data->levels[0]->pc, PCASM)); 2329 if (!ctx) { 2330 if (!data->levels[0]->pc->setupcalled) { 2331 IS sorted; /* PCASM will sort the input IS, duplicate it to return an unmodified (PCHPDDM) input IS */ 2332 PetscCall(ISDuplicate(is[0], &sorted)); 2333 PetscCall(PCASMSetLocalSubdomains(data->levels[0]->pc, 1, &sorted, &loc)); 2334 PetscCall(PetscObjectDereference((PetscObject)sorted)); 2335 } 2336 PetscCall(PCSetFromOptions(data->levels[0]->pc)); 2337 if (block) { 2338 PetscCall(PCHPDDMPermute_Private(unsorted, data->is, &uis, sub[0], &C, &perm)); 2339 PetscCall(PCHPDDMCommunicationAvoidingPCASM_Private(data->levels[0]->pc, C, algebraic)); 2340 } else PetscCall(PCSetUp(data->levels[0]->pc)); 2341 PetscTryMethod(data->levels[0]->pc, "PCASMGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (data->levels[0]->pc, &size, nullptr, &ksp)); 2342 if (size != 1) { 2343 data->share = PETSC_FALSE; 2344 PetscCheck(size == -1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of subdomain solver %" PetscInt_FMT " != 1", size); 2345 PetscCall(PetscInfo(pc, "Cannot share subdomain KSP between SLEPc and PETSc since PCASMGetSubKSP() not found in fine-level PC\n")); 2346 PetscCall(ISDestroy(&unsorted)); 2347 unsorted = is[0]; 2348 } else { 2349 const char *matpre; 2350 PetscBool cmp[4]; 2351 2352 if (!block && !ctx) PetscCall(PCHPDDMPermute_Private(unsorted, data->is, &uis, PetscBool3ToBool(data->Neumann) ? sub[0] : data->aux, &C, &perm)); 2353 if (perm) { /* unsorted input IS */ 2354 if (!PetscBool3ToBool(data->Neumann) && !block) { 2355 PetscCall(MatPermute(sub[0], perm, perm, &D)); /* permute since PCASM will call ISSort() */ 2356 PetscCall(MatHeaderReplace(sub[0], &D)); 2357 } 2358 if (data->B) { /* see PCHPDDMSetRHSMat() */ 2359 PetscCall(MatPermute(data->B, perm, perm, &D)); 2360 PetscCall(MatHeaderReplace(data->B, &D)); 2361 } 2362 PetscCall(ISDestroy(&perm)); 2363 } 2364 PetscCall(KSPGetOperators(ksp[0], subA, subA + 1)); 2365 PetscCall(PetscObjectReference((PetscObject)subA[0])); 2366 PetscCall(MatDuplicate(subA[1], MAT_SHARE_NONZERO_PATTERN, &D)); 2367 PetscCall(MatGetOptionsPrefix(subA[1], &matpre)); 2368 PetscCall(MatSetOptionsPrefix(D, matpre)); 2369 PetscCall(PetscObjectTypeCompare((PetscObject)D, MATNORMAL, cmp)); 2370 PetscCall(PetscObjectTypeCompare((PetscObject)C, MATNORMAL, cmp + 1)); 2371 if (!cmp[0]) PetscCall(PetscObjectTypeCompare((PetscObject)D, MATNORMALHERMITIAN, cmp + 2)); 2372 else cmp[2] = PETSC_FALSE; 2373 if (!cmp[1]) PetscCall(PetscObjectTypeCompare((PetscObject)C, MATNORMALHERMITIAN, cmp + 3)); 2374 else cmp[3] = PETSC_FALSE; 2375 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); 2376 if (!cmp[0] && !cmp[2]) { 2377 if (!block) PetscCall(MatAXPY(D, 1.0, C, SUBSET_NONZERO_PATTERN)); 2378 else { 2379 PetscCall(MatMissingDiagonal(D, cmp, nullptr)); 2380 if (cmp[0]) structure = DIFFERENT_NONZERO_PATTERN; /* data->aux has no missing diagonal entry */ 2381 PetscCall(MatAXPY(D, 1.0, data->aux, structure)); 2382 } 2383 } else { 2384 Mat mat[2]; 2385 2386 if (cmp[0]) { 2387 PetscCall(MatNormalGetMat(D, mat)); 2388 PetscCall(MatNormalGetMat(C, mat + 1)); 2389 } else { 2390 PetscCall(MatNormalHermitianGetMat(D, mat)); 2391 PetscCall(MatNormalHermitianGetMat(C, mat + 1)); 2392 } 2393 PetscCall(MatAXPY(mat[0], 1.0, mat[1], SUBSET_NONZERO_PATTERN)); 2394 } 2395 PetscCall(MatPropagateSymmetryOptions(C, D)); 2396 PetscCall(MatDestroy(&C)); 2397 C = D; 2398 /* swap pointers so that variables stay consistent throughout PCSetUp() */ 2399 std::swap(C, data->aux); 2400 std::swap(uis, data->is); 2401 swap = PETSC_TRUE; 2402 } 2403 } 2404 } 2405 } 2406 if (ctx) { 2407 PC_HPDDM *data_00 = (PC_HPDDM *)std::get<0>(*ctx)[0]->data; 2408 PC s; 2409 Mat A00, P00, A01 = nullptr, A10, A11, N, b[4]; 2410 IS sorted, is[2], *is_00; 2411 MatSolverType type; 2412 std::pair<PC, Vec[2]> *p; 2413 2414 n = -1; 2415 PetscTryMethod(data_00->levels[0]->pc, "PCASMGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (data_00->levels[0]->pc, &n, nullptr, &ksp)); 2416 PetscCheck(n == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of subdomain solver %" PetscInt_FMT " != 1", n); 2417 PetscCall(KSPGetOperators(ksp[0], subA, subA + 1)); 2418 PetscCall(ISGetLocalSize(data_00->is, &n)); 2419 if (n != subA[0]->rmap->n || n != subA[0]->cmap->n) { 2420 PetscCall(PCASMGetLocalSubdomains(data_00->levels[0]->pc, &n, &is_00, nullptr)); 2421 PetscCall(ISGetLocalSize(*is_00, &n)); 2422 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); 2423 } else is_00 = &data_00->is; 2424 PetscCall(PCHPDDMPermute_Private(unsorted, data->is, &uis, data->aux, &C, nullptr)); /* permute since PCASM works with a sorted IS */ 2425 std::swap(C, data->aux); 2426 std::swap(uis, data->is); 2427 swap = PETSC_TRUE; 2428 PetscCall(MatSchurComplementGetSubMatrices(P, &A00, &P00, std::get<1>(*ctx), &A10, &A11)); 2429 std::get<1>(*ctx)[1] = A10; 2430 PetscCall(PetscObjectTypeCompare((PetscObject)A10, MATTRANSPOSEVIRTUAL, &flg)); 2431 if (flg) PetscCall(MatTransposeGetMat(A10, &A01)); 2432 else { 2433 PetscBool flg; 2434 2435 PetscCall(PetscObjectTypeCompare((PetscObject)A10, MATHERMITIANTRANSPOSEVIRTUAL, &flg)); 2436 if (flg) PetscCall(MatHermitianTransposeGetMat(A10, &A01)); 2437 } 2438 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 */ 2439 PetscCall(ISSort(sorted)); /* this is to avoid changing users inputs, but it requires a new call to ISSort() here */ 2440 if (!A01) { 2441 PetscCall(MatSetOption(A10, MAT_SUBMAT_SINGLEIS, PETSC_TRUE)); 2442 PetscCall(MatCreateSubMatrices(A10, 1, &data->is, &sorted, MAT_INITIAL_MATRIX, &sub)); 2443 b[2] = sub[0]; 2444 PetscCall(PetscObjectReference((PetscObject)sub[0])); 2445 PetscCall(MatDestroySubMatrices(1, &sub)); 2446 PetscCall(PetscObjectTypeCompare((PetscObject)std::get<1>(*ctx)[0], MATTRANSPOSEVIRTUAL, &flg)); 2447 A10 = nullptr; 2448 if (flg) PetscCall(MatTransposeGetMat(std::get<1>(*ctx)[0], &A10)); 2449 else { 2450 PetscBool flg; 2451 2452 PetscCall(PetscObjectTypeCompare((PetscObject)std::get<1>(*ctx)[0], MATHERMITIANTRANSPOSEVIRTUAL, &flg)); 2453 if (flg) PetscCall(MatHermitianTransposeGetMat(std::get<1>(*ctx)[0], &A10)); 2454 } 2455 if (!A10) PetscCall(MatCreateSubMatrices(std::get<1>(*ctx)[0], 1, &sorted, &data->is, MAT_INITIAL_MATRIX, &sub)); 2456 else { 2457 if (flg) PetscCall(MatCreateTranspose(b[2], b + 1)); 2458 else PetscCall(MatCreateHermitianTranspose(b[2], b + 1)); 2459 } 2460 } else { 2461 PetscCall(MatSetOption(A01, MAT_SUBMAT_SINGLEIS, PETSC_TRUE)); 2462 PetscCall(MatCreateSubMatrices(A01, 1, &sorted, &data->is, MAT_INITIAL_MATRIX, &sub)); 2463 if (flg) PetscCall(MatCreateTranspose(*sub, b + 2)); 2464 else PetscCall(MatCreateHermitianTranspose(*sub, b + 2)); 2465 } 2466 if (A01 || !A10) { 2467 b[1] = sub[0]; 2468 PetscCall(PetscObjectReference((PetscObject)sub[0])); 2469 } 2470 PetscCall(MatDestroySubMatrices(1, &sub)); 2471 PetscCall(ISDestroy(&sorted)); 2472 b[3] = data->aux; 2473 PetscCall(MatCreateSchurComplement(subA[0], subA[1], b[1], b[2], b[3], &S)); 2474 PetscCall(MatSchurComplementSetKSP(S, ksp[0])); 2475 if (data->N != 1) { 2476 PetscCall(PCASMSetType(data->levels[0]->pc, PC_ASM_NONE)); /* "Neumann--Neumann" preconditioning with overlap and a Boolean partition of unity */ 2477 PetscCall(PCASMSetLocalSubdomains(data->levels[0]->pc, 1, &data->is, &loc)); 2478 PetscCall(PCSetFromOptions(data->levels[0]->pc)); /* action of eq. (15) of https://hal.science/hal-02343808v6/document (with a sign flip) */ 2479 s = data->levels[0]->pc; 2480 } else { 2481 is[0] = data->is; 2482 PetscCall(PetscObjectReference((PetscObject)is[0])); 2483 PetscCall(PetscObjectReference((PetscObject)b[3])); 2484 PetscCall(PCSetType(pc, PCASM)); /* change the type of the current PC */ 2485 data = nullptr; /* destroyed in the previous PCSetType(), so reset to NULL to avoid any faulty use */ 2486 PetscCall(PCAppendOptionsPrefix(pc, "pc_hpddm_coarse_")); /* same prefix as when using PCHPDDM with a single level */ 2487 PetscCall(PCASMSetLocalSubdomains(pc, 1, is, &loc)); 2488 PetscCall(ISDestroy(is)); 2489 PetscCall(ISDestroy(&loc)); 2490 s = pc; 2491 } 2492 PetscCall(PCHPDDMCommunicationAvoidingPCASM_Private(s, S, PETSC_TRUE)); /* the subdomain Mat is already known and the input IS of PCASMSetLocalSubdomains() is already sorted */ 2493 PetscTryMethod(s, "PCASMGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (s, &n, nullptr, &ksp)); 2494 PetscCheck(n == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of subdomain solver %" PetscInt_FMT " != 1", n); 2495 PetscCall(KSPGetPC(ksp[0], &inner)); 2496 PetscCall(PCSetType(inner, PCSHELL)); /* compute the action of the inverse of the local Schur complement with a PCSHELL */ 2497 b[0] = subA[0]; 2498 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 */ 2499 if (!data) PetscCall(PetscObjectDereference((PetscObject)b[3])); 2500 PetscCall(PetscObjectDereference((PetscObject)b[1])); 2501 PetscCall(PetscObjectDereference((PetscObject)b[2])); 2502 PetscCall(PCCreate(PETSC_COMM_SELF, &s)); 2503 PetscCall(PCSetOptionsPrefix(s, ((PetscObject)inner)->prefix)); 2504 PetscCall(PCSetOptionsPrefix(inner, nullptr)); 2505 PetscCall(KSPSetSkipPCSetFromOptions(ksp[0], PETSC_TRUE)); 2506 PetscCall(PCSetType(s, PCLU)); 2507 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 */ 2508 PetscCall(PCSetFromOptions(s)); 2509 PetscCall(PCFactorGetMatSolverType(s, &type)); 2510 PetscCall(PetscStrcmp(type, MATSOLVERMUMPS, &flg)); 2511 PetscCall(MatGetLocalSize(A11, &n, nullptr)); 2512 if (flg || n == 0) { 2513 PetscCall(PCSetOperators(s, N, N)); 2514 if (n) { 2515 PetscCall(PCFactorGetMatrix(s, b)); 2516 PetscCall(MatSetOptionsPrefix(*b, ((PetscObject)s)->prefix)); 2517 n = -1; 2518 PetscCall(PetscOptionsGetInt(nullptr, ((PetscObject)s)->prefix, "-mat_mumps_icntl_26", &n, nullptr)); 2519 if (n == 1) { /* allocates a square MatDense of size is[1]->map->n, so one */ 2520 PetscCall(MatNestGetISs(N, is, nullptr)); /* needs to be able to deactivate this path when dealing */ 2521 PetscCall(MatFactorSetSchurIS(*b, is[1])); /* with a large constraint space in order to avoid OOM */ 2522 } 2523 } else PetscCall(PCSetType(s, PCNONE)); /* empty local Schur complement (e.g., centralized on another process) */ 2524 } else { 2525 PetscCall(MatConvert(N, MATAIJ, MAT_INITIAL_MATRIX, b)); 2526 PetscCall(PCSetOperators(s, N, *b)); 2527 PetscCall(PetscObjectDereference((PetscObject)*b)); 2528 PetscCall(PetscObjectTypeCompareAny((PetscObject)s, &flg, PCLU, PCCHOLESKY, PCILU, PCICC, PCQR, "")); 2529 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 */ 2530 } 2531 PetscCall(PetscNew(&p)); 2532 p->first = s; 2533 if (n != 0) PetscCall(MatCreateVecs(*b, p->second, p->second + 1)); 2534 else p->second[0] = p->second[1] = nullptr; 2535 PetscCall(PCShellSetContext(inner, p)); 2536 PetscCall(PCShellSetApply(inner, PCApply_Nest)); 2537 PetscCall(PCShellSetView(inner, PCView_Nest)); 2538 PetscCall(PCShellSetDestroy(inner, PCDestroy_Nest)); 2539 PetscCall(PetscObjectDereference((PetscObject)N)); 2540 if (!data) { 2541 PetscCall(MatDestroy(&S)); 2542 PetscCall(ISDestroy(&unsorted)); 2543 PetscCall(MatDestroy(&C)); 2544 PetscCall(ISDestroy(&uis)); 2545 PetscCall(PetscFree(ctx)); 2546 #if PetscDefined(USE_DEBUG) 2547 PetscCall(ISDestroy(&dis)); 2548 PetscCall(MatDestroy(&daux)); 2549 #endif 2550 PetscFunctionReturn(PETSC_SUCCESS); 2551 } 2552 } 2553 if (!data->levels[0]->scatter) { 2554 PetscCall(MatCreateVecs(P, &xin, nullptr)); 2555 if (ismatis) PetscCall(MatDestroy(&P)); 2556 PetscCall(VecScatterCreate(xin, data->is, data->levels[0]->D, nullptr, &data->levels[0]->scatter)); 2557 PetscCall(VecDestroy(&xin)); 2558 } 2559 if (data->levels[0]->P) { 2560 /* if the pattern is the same and PCSetUp() has previously succeeded, reuse HPDDM buffers and connectivity */ 2561 PetscCall(HPDDM::Schwarz<PetscScalar>::destroy(data->levels[0], !pc->setupcalled || pc->flag == DIFFERENT_NONZERO_PATTERN ? PETSC_TRUE : PETSC_FALSE)); 2562 } 2563 if (!data->levels[0]->P) data->levels[0]->P = new HPDDM::Schwarz<PetscScalar>(); 2564 if (data->log_separate) PetscCall(PetscLogEventBegin(PC_HPDDM_SetUp[0], data->levels[0]->ksp, nullptr, nullptr, nullptr)); 2565 else PetscCall(PetscLogEventBegin(PC_HPDDM_Strc, data->levels[0]->ksp, nullptr, nullptr, nullptr)); 2566 /* HPDDM internal data structure */ 2567 PetscCall(data->levels[0]->P->structure(loc, data->is, !ctx ? sub[0] : nullptr, ismatis ? C : data->aux, data->levels)); 2568 if (!data->log_separate) PetscCall(PetscLogEventEnd(PC_HPDDM_Strc, data->levels[0]->ksp, nullptr, nullptr, nullptr)); 2569 /* matrix pencil of the generalized eigenvalue problem on the overlap (GenEO) */ 2570 if (!ctx) { 2571 if (data->deflation || overlap != -1) weighted = data->aux; 2572 else if (!data->B) { 2573 PetscBool cmp; 2574 2575 PetscCall(MatDuplicate(sub[0], MAT_COPY_VALUES, &weighted)); 2576 PetscCall(PetscObjectTypeCompareAny((PetscObject)weighted, &cmp, MATNORMAL, MATNORMALHERMITIAN, "")); 2577 if (cmp) flg = PETSC_FALSE; 2578 PetscCall(MatDiagonalScale(weighted, data->levels[0]->D, data->levels[0]->D)); 2579 /* neither MatDuplicate() nor MatDiagonalScale() handles the symmetry options, so propagate the options explicitly */ 2580 /* only useful for -mat_type baij -pc_hpddm_levels_1_st_pc_type cholesky (no problem with MATAIJ or MATSBAIJ) */ 2581 PetscCall(MatPropagateSymmetryOptions(sub[0], weighted)); 2582 } else weighted = data->B; 2583 } else weighted = nullptr; 2584 /* SLEPc is used inside the loaded symbol */ 2585 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)); 2586 if (!ctx && data->share && overlap == -1) { 2587 Mat st[2]; 2588 PetscCall(KSPGetOperators(ksp[0], st, st + 1)); 2589 PetscCall(MatCopy(subA[0], st[0], structure)); 2590 if (subA[1] != subA[0] || st[1] != st[0]) PetscCall(MatCopy(subA[1], st[1], SAME_NONZERO_PATTERN)); 2591 PetscCall(PetscObjectDereference((PetscObject)subA[0])); 2592 } 2593 if (data->log_separate) PetscCall(PetscLogEventEnd(PC_HPDDM_SetUp[0], data->levels[0]->ksp, nullptr, nullptr, nullptr)); 2594 if (ismatis) PetscCall(MatISGetLocalMat(C, &N)); 2595 else N = data->aux; 2596 if (!ctx) P = sub[0]; 2597 else P = S; 2598 /* going through the grid hierarchy */ 2599 for (n = 1; n < data->N; ++n) { 2600 if (data->log_separate) PetscCall(PetscLogEventBegin(PC_HPDDM_SetUp[n], data->levels[n]->ksp, nullptr, nullptr, nullptr)); 2601 /* method composed in the loaded symbol since there, SLEPc is used as well */ 2602 PetscTryMethod(data->levels[0]->ksp, "PCHPDDMSetUp_Private_C", (Mat *, Mat *, PetscInt, PetscInt *const, PC_HPDDM_Level **const), (&P, &N, n, &data->N, data->levels)); 2603 if (data->log_separate) PetscCall(PetscLogEventEnd(PC_HPDDM_SetUp[n], data->levels[n]->ksp, nullptr, nullptr, nullptr)); 2604 } 2605 /* reset to NULL to avoid any faulty use */ 2606 PetscCall(PetscObjectComposeFunction((PetscObject)data->levels[0]->ksp, "PCHPDDMSetUp_Private_C", nullptr)); 2607 if (!ismatis) PetscCall(PetscObjectComposeFunction((PetscObject)pc->pmat, "PCHPDDMAlgebraicAuxiliaryMat_C", nullptr)); 2608 else PetscCall(PetscObjectDereference((PetscObject)C)); /* matching PetscObjectReference() above */ 2609 for (n = 0; n < data->N - 1; ++n) 2610 if (data->levels[n]->P) { 2611 /* HPDDM internal work buffers */ 2612 data->levels[n]->P->setBuffer(); 2613 data->levels[n]->P->super::start(); 2614 } 2615 if (ismatis || !subdomains) PetscCall(PCHPDDMDestroySubMatrices_Private(PetscBool3ToBool(data->Neumann), PetscBool(algebraic && !block && overlap == -1), sub)); 2616 if (ismatis) data->is = nullptr; 2617 for (n = 0; n < data->N - 1 + (reused > 0); ++n) { 2618 if (data->levels[n]->P) { 2619 PC spc; 2620 2621 /* force the PC to be PCSHELL to do the coarse grid corrections */ 2622 PetscCall(KSPSetSkipPCSetFromOptions(data->levels[n]->ksp, PETSC_TRUE)); 2623 PetscCall(KSPGetPC(data->levels[n]->ksp, &spc)); 2624 PetscCall(PCSetType(spc, PCSHELL)); 2625 PetscCall(PCShellSetContext(spc, data->levels[n])); 2626 PetscCall(PCShellSetSetUp(spc, PCSetUp_HPDDMShell)); 2627 PetscCall(PCShellSetApply(spc, PCApply_HPDDMShell)); 2628 PetscCall(PCShellSetMatApply(spc, PCMatApply_HPDDMShell)); 2629 if (ctx && n == 0) { 2630 Mat Amat, Pmat; 2631 PetscInt m, M; 2632 std::tuple<Mat, PetscSF, Vec[2]> *ctx; 2633 2634 PetscCall(KSPGetOperators(data->levels[n]->ksp, nullptr, &Pmat)); 2635 PetscCall(MatGetLocalSize(Pmat, &m, nullptr)); 2636 PetscCall(MatGetSize(Pmat, &M, nullptr)); 2637 PetscCall(PetscNew(&ctx)); 2638 std::get<0>(*ctx) = S; 2639 std::get<1>(*ctx) = data->levels[n]->scatter; 2640 PetscCall(MatCreateShell(PetscObjectComm((PetscObject)data->levels[n]->ksp), m, m, M, M, ctx, &Amat)); 2641 PetscCall(MatShellSetOperation(Amat, MATOP_MULT, (void (*)(void))MatMult_Schur<false>)); 2642 PetscCall(MatShellSetOperation(Amat, MATOP_MULT_TRANSPOSE, (void (*)(void))MatMult_Schur<true>)); 2643 PetscCall(MatShellSetOperation(Amat, MATOP_DESTROY, (void (*)(void))MatDestroy_Schur)); 2644 PetscCall(MatCreateVecs(S, std::get<2>(*ctx), std::get<2>(*ctx) + 1)); 2645 PetscCall(KSPSetOperators(data->levels[n]->ksp, Amat, Pmat)); 2646 PetscCall(PetscObjectDereference((PetscObject)Amat)); 2647 } 2648 PetscCall(PCShellSetDestroy(spc, PCDestroy_HPDDMShell)); 2649 if (!data->levels[n]->pc) PetscCall(PCCreate(PetscObjectComm((PetscObject)data->levels[n]->ksp), &data->levels[n]->pc)); 2650 if (n < reused) { 2651 PetscCall(PCSetReusePreconditioner(spc, PETSC_TRUE)); 2652 PetscCall(PCSetReusePreconditioner(data->levels[n]->pc, PETSC_TRUE)); 2653 } 2654 PetscCall(PCSetUp(spc)); 2655 } 2656 } 2657 if (ctx) PetscCall(MatDestroy(&S)); 2658 if (overlap == -1) PetscCall(PetscObjectComposeFunction((PetscObject)pc->pmat, "PCHPDDMAlgebraicAuxiliaryMat_Private_C", nullptr)); 2659 } else flg = reused ? PETSC_FALSE : PETSC_TRUE; 2660 if (!ismatis && subdomains) { 2661 if (flg) PetscCall(KSPGetPC(data->levels[0]->ksp, &inner)); 2662 else inner = data->levels[0]->pc; 2663 if (inner) { 2664 if (!inner->setupcalled) PetscCall(PCSetType(inner, PCASM)); 2665 PetscCall(PCSetFromOptions(inner)); 2666 PetscCall(PetscStrcmp(((PetscObject)inner)->type_name, PCASM, &flg)); 2667 if (flg) { 2668 if (!inner->setupcalled) { /* evaluates to PETSC_FALSE when -pc_hpddm_block_splitting */ 2669 IS sorted; /* PCASM will sort the input IS, duplicate it to return an unmodified (PCHPDDM) input IS */ 2670 2671 PetscCall(ISDuplicate(is[0], &sorted)); 2672 PetscCall(PCASMSetLocalSubdomains(inner, 1, &sorted, &loc)); 2673 PetscCall(PetscObjectDereference((PetscObject)sorted)); 2674 } 2675 if (!PetscBool3ToBool(data->Neumann) && data->N > 1) { /* subdomain matrices are already created for the eigenproblem, reuse them for the fine-level PC */ 2676 PetscCall(PCHPDDMPermute_Private(*is, nullptr, nullptr, sub[0], &P, nullptr)); 2677 PetscCall(PCHPDDMCommunicationAvoidingPCASM_Private(inner, P, algebraic)); 2678 PetscCall(PetscObjectDereference((PetscObject)P)); 2679 } 2680 } 2681 } 2682 if (data->N > 1) { 2683 if (overlap != 1) PetscCall(PCHPDDMDestroySubMatrices_Private(PetscBool3ToBool(data->Neumann), PetscBool(algebraic && !block && overlap == -1), sub)); 2684 if (overlap == 1) PetscCall(MatDestroy(subA)); 2685 } 2686 } 2687 PetscCall(ISDestroy(&loc)); 2688 } else data->N = 1 + reused; /* enforce this value to 1 + reused if there is no way to build another level */ 2689 if (requested != data->N + reused) { 2690 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, 2691 data->N, pcpre ? pcpre : "")); 2692 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)); 2693 /* cannot use PCDestroy_HPDDMShell() because PCSHELL not set for unassembled levels */ 2694 for (n = data->N - 1; n < requested - 1; ++n) { 2695 if (data->levels[n]->P) { 2696 PetscCall(HPDDM::Schwarz<PetscScalar>::destroy(data->levels[n], PETSC_TRUE)); 2697 PetscCall(VecDestroyVecs(1, &data->levels[n]->v[0])); 2698 PetscCall(VecDestroyVecs(2, &data->levels[n]->v[1])); 2699 PetscCall(MatDestroy(data->levels[n]->V)); 2700 PetscCall(MatDestroy(data->levels[n]->V + 1)); 2701 PetscCall(MatDestroy(data->levels[n]->V + 2)); 2702 PetscCall(VecDestroy(&data->levels[n]->D)); 2703 PetscCall(PetscSFDestroy(&data->levels[n]->scatter)); 2704 } 2705 } 2706 if (reused) { 2707 for (n = reused; n < PETSC_PCHPDDM_MAXLEVELS && data->levels[n]; ++n) { 2708 PetscCall(KSPDestroy(&data->levels[n]->ksp)); 2709 PetscCall(PCDestroy(&data->levels[n]->pc)); 2710 } 2711 } 2712 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, 2713 data->N, reused, data->N, pcpre ? pcpre : "", pcpre ? pcpre : "", data->N); 2714 } 2715 /* these solvers are created after PCSetFromOptions() is called */ 2716 if (pc->setfromoptionscalled) { 2717 for (n = 0; n < data->N; ++n) { 2718 if (data->levels[n]->ksp) PetscCall(KSPSetFromOptions(data->levels[n]->ksp)); 2719 if (data->levels[n]->pc) PetscCall(PCSetFromOptions(data->levels[n]->pc)); 2720 } 2721 pc->setfromoptionscalled = 0; 2722 } 2723 data->N += reused; 2724 if (data->share && swap) { 2725 /* swap back pointers so that variables follow the user-provided numbering */ 2726 std::swap(C, data->aux); 2727 std::swap(uis, data->is); 2728 PetscCall(MatDestroy(&C)); 2729 PetscCall(ISDestroy(&uis)); 2730 } 2731 if (algebraic) PetscCall(MatDestroy(&data->aux)); 2732 if (unsorted && unsorted != is[0]) { 2733 PetscCall(ISCopy(unsorted, data->is)); 2734 PetscCall(ISDestroy(&unsorted)); 2735 } 2736 #if PetscDefined(USE_DEBUG) 2737 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); 2738 if (data->is) { 2739 PetscCall(ISEqualUnsorted(data->is, dis, &flg)); 2740 PetscCall(ISDestroy(&dis)); 2741 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Input IS and output IS are not equal"); 2742 } 2743 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); 2744 if (data->aux) { 2745 PetscCall(MatMultEqual(data->aux, daux, 20, &flg)); 2746 PetscCall(MatDestroy(&daux)); 2747 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Input Mat and output Mat are not equal"); 2748 } 2749 #endif 2750 PetscFunctionReturn(PETSC_SUCCESS); 2751 } 2752 2753 /*@ 2754 PCHPDDMSetCoarseCorrectionType - Sets the coarse correction type. 2755 2756 Collective 2757 2758 Input Parameters: 2759 + pc - preconditioner context 2760 - type - `PC_HPDDM_COARSE_CORRECTION_DEFLATED`, `PC_HPDDM_COARSE_CORRECTION_ADDITIVE`, `PC_HPDDM_COARSE_CORRECTION_BALANCED`, or `PC_HPDDM_COARSE_CORRECTION_NONE` 2761 2762 Options Database Key: 2763 . -pc_hpddm_coarse_correction <deflated, additive, balanced, none> - type of coarse correction to apply 2764 2765 Level: intermediate 2766 2767 .seealso: [](ch_ksp), `PCHPDDMGetCoarseCorrectionType()`, `PCHPDDM`, `PCHPDDMCoarseCorrectionType` 2768 @*/ 2769 PetscErrorCode PCHPDDMSetCoarseCorrectionType(PC pc, PCHPDDMCoarseCorrectionType type) 2770 { 2771 PetscFunctionBegin; 2772 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2773 PetscValidLogicalCollectiveEnum(pc, type, 2); 2774 PetscTryMethod(pc, "PCHPDDMSetCoarseCorrectionType_C", (PC, PCHPDDMCoarseCorrectionType), (pc, type)); 2775 PetscFunctionReturn(PETSC_SUCCESS); 2776 } 2777 2778 /*@ 2779 PCHPDDMGetCoarseCorrectionType - Gets the coarse correction type. 2780 2781 Input Parameter: 2782 . pc - preconditioner context 2783 2784 Output Parameter: 2785 . type - `PC_HPDDM_COARSE_CORRECTION_DEFLATED`, `PC_HPDDM_COARSE_CORRECTION_ADDITIVE`, `PC_HPDDM_COARSE_CORRECTION_BALANCED`, or `PC_HPDDM_COARSE_CORRECTION_NONE` 2786 2787 Level: intermediate 2788 2789 .seealso: [](ch_ksp), `PCHPDDMSetCoarseCorrectionType()`, `PCHPDDM`, `PCHPDDMCoarseCorrectionType` 2790 @*/ 2791 PetscErrorCode PCHPDDMGetCoarseCorrectionType(PC pc, PCHPDDMCoarseCorrectionType *type) 2792 { 2793 PetscFunctionBegin; 2794 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2795 if (type) { 2796 PetscAssertPointer(type, 2); 2797 PetscUseMethod(pc, "PCHPDDMGetCoarseCorrectionType_C", (PC, PCHPDDMCoarseCorrectionType *), (pc, type)); 2798 } 2799 PetscFunctionReturn(PETSC_SUCCESS); 2800 } 2801 2802 static PetscErrorCode PCHPDDMSetCoarseCorrectionType_HPDDM(PC pc, PCHPDDMCoarseCorrectionType type) 2803 { 2804 PC_HPDDM *data = (PC_HPDDM *)pc->data; 2805 2806 PetscFunctionBegin; 2807 data->correction = type; 2808 PetscFunctionReturn(PETSC_SUCCESS); 2809 } 2810 2811 static PetscErrorCode PCHPDDMGetCoarseCorrectionType_HPDDM(PC pc, PCHPDDMCoarseCorrectionType *type) 2812 { 2813 PC_HPDDM *data = (PC_HPDDM *)pc->data; 2814 2815 PetscFunctionBegin; 2816 *type = data->correction; 2817 PetscFunctionReturn(PETSC_SUCCESS); 2818 } 2819 2820 /*@ 2821 PCHPDDMSetSTShareSubKSP - Sets whether the `KSP` in SLEPc `ST` and the fine-level subdomain solver should be shared. 2822 2823 Input Parameters: 2824 + pc - preconditioner context 2825 - share - whether the `KSP` should be shared or not 2826 2827 Note: 2828 This is not the same as `PCSetReusePreconditioner()`. Given certain conditions (visible using -info), a symbolic factorization can be skipped 2829 when using a subdomain `PCType` such as `PCLU` or `PCCHOLESKY`. 2830 2831 Level: advanced 2832 2833 .seealso: [](ch_ksp), `PCHPDDM`, `PCHPDDMGetSTShareSubKSP()` 2834 @*/ 2835 PetscErrorCode PCHPDDMSetSTShareSubKSP(PC pc, PetscBool share) 2836 { 2837 PetscFunctionBegin; 2838 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2839 PetscTryMethod(pc, "PCHPDDMSetSTShareSubKSP_C", (PC, PetscBool), (pc, share)); 2840 PetscFunctionReturn(PETSC_SUCCESS); 2841 } 2842 2843 /*@ 2844 PCHPDDMGetSTShareSubKSP - Gets whether the `KSP` in SLEPc `ST` and the fine-level subdomain solver is shared. 2845 2846 Input Parameter: 2847 . pc - preconditioner context 2848 2849 Output Parameter: 2850 . share - whether the `KSP` is shared or not 2851 2852 Note: 2853 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 2854 when using a subdomain `PCType` such as `PCLU` or `PCCHOLESKY`. 2855 2856 Level: advanced 2857 2858 .seealso: [](ch_ksp), `PCHPDDM`, `PCHPDDMSetSTShareSubKSP()` 2859 @*/ 2860 PetscErrorCode PCHPDDMGetSTShareSubKSP(PC pc, PetscBool *share) 2861 { 2862 PetscFunctionBegin; 2863 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2864 if (share) { 2865 PetscAssertPointer(share, 2); 2866 PetscUseMethod(pc, "PCHPDDMGetSTShareSubKSP_C", (PC, PetscBool *), (pc, share)); 2867 } 2868 PetscFunctionReturn(PETSC_SUCCESS); 2869 } 2870 2871 static PetscErrorCode PCHPDDMSetSTShareSubKSP_HPDDM(PC pc, PetscBool share) 2872 { 2873 PC_HPDDM *data = (PC_HPDDM *)pc->data; 2874 2875 PetscFunctionBegin; 2876 data->share = share; 2877 PetscFunctionReturn(PETSC_SUCCESS); 2878 } 2879 2880 static PetscErrorCode PCHPDDMGetSTShareSubKSP_HPDDM(PC pc, PetscBool *share) 2881 { 2882 PC_HPDDM *data = (PC_HPDDM *)pc->data; 2883 2884 PetscFunctionBegin; 2885 *share = data->share; 2886 PetscFunctionReturn(PETSC_SUCCESS); 2887 } 2888 2889 /*@ 2890 PCHPDDMSetDeflationMat - Sets the deflation space used to assemble a coarser operator. 2891 2892 Input Parameters: 2893 + pc - preconditioner context 2894 . is - index set of the local deflation matrix 2895 - U - deflation sequential matrix stored as a `MATSEQDENSE` 2896 2897 Level: advanced 2898 2899 .seealso: [](ch_ksp), `PCHPDDM`, `PCDeflationSetSpace()`, `PCMGSetRestriction()` 2900 @*/ 2901 PetscErrorCode PCHPDDMSetDeflationMat(PC pc, IS is, Mat U) 2902 { 2903 PetscFunctionBegin; 2904 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2905 PetscValidHeaderSpecific(is, IS_CLASSID, 2); 2906 PetscValidHeaderSpecific(U, MAT_CLASSID, 3); 2907 PetscTryMethod(pc, "PCHPDDMSetDeflationMat_C", (PC, IS, Mat), (pc, is, U)); 2908 PetscFunctionReturn(PETSC_SUCCESS); 2909 } 2910 2911 static PetscErrorCode PCHPDDMSetDeflationMat_HPDDM(PC pc, IS is, Mat U) 2912 { 2913 PetscFunctionBegin; 2914 PetscCall(PCHPDDMSetAuxiliaryMat_Private(pc, is, U, PETSC_TRUE)); 2915 PetscFunctionReturn(PETSC_SUCCESS); 2916 } 2917 2918 PetscErrorCode HPDDMLoadDL_Private(PetscBool *found) 2919 { 2920 PetscBool flg; 2921 char lib[PETSC_MAX_PATH_LEN], dlib[PETSC_MAX_PATH_LEN], dir[PETSC_MAX_PATH_LEN]; 2922 2923 PetscFunctionBegin; 2924 PetscAssertPointer(found, 1); 2925 PetscCall(PetscStrncpy(dir, "${PETSC_LIB_DIR}", sizeof(dir))); 2926 PetscCall(PetscOptionsGetString(nullptr, nullptr, "-hpddm_dir", dir, sizeof(dir), nullptr)); 2927 PetscCall(PetscSNPrintf(lib, sizeof(lib), "%s/libhpddm_petsc", dir)); 2928 PetscCall(PetscDLLibraryRetrieve(PETSC_COMM_SELF, lib, dlib, 1024, found)); 2929 #if defined(SLEPC_LIB_DIR) /* this variable is passed during SLEPc ./configure when PETSc has not been configured */ 2930 if (!*found) { /* with --download-hpddm since slepcconf.h is not yet built (and thus can't be included) */ 2931 PetscCall(PetscStrncpy(dir, HPDDM_STR(SLEPC_LIB_DIR), sizeof(dir))); 2932 PetscCall(PetscSNPrintf(lib, sizeof(lib), "%s/libhpddm_petsc", dir)); 2933 PetscCall(PetscDLLibraryRetrieve(PETSC_COMM_SELF, lib, dlib, 1024, found)); 2934 } 2935 #endif 2936 if (!*found) { /* probable options for this to evaluate to PETSC_TRUE: system inconsistency (libhpddm_petsc moved by user?) or PETSc configured without --download-slepc */ 2937 PetscCall(PetscOptionsGetenv(PETSC_COMM_SELF, "SLEPC_DIR", dir, sizeof(dir), &flg)); 2938 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 */ 2939 PetscCall(PetscSNPrintf(lib, sizeof(lib), "%s/lib/libslepc", dir)); 2940 PetscCall(PetscDLLibraryRetrieve(PETSC_COMM_SELF, lib, dlib, 1024, found)); 2941 PetscCheck(*found, PETSC_COMM_SELF, PETSC_ERR_PLIB, "%s not found but SLEPC_DIR=%s", lib, dir); 2942 PetscCall(PetscDLLibraryAppend(PETSC_COMM_SELF, &PetscDLLibrariesLoaded, dlib)); 2943 PetscCall(PetscSNPrintf(lib, sizeof(lib), "%s/lib/libhpddm_petsc", dir)); /* libhpddm_petsc is always in the same directory as libslepc */ 2944 PetscCall(PetscDLLibraryRetrieve(PETSC_COMM_SELF, lib, dlib, 1024, found)); 2945 } 2946 } 2947 PetscCheck(*found, PETSC_COMM_SELF, PETSC_ERR_PLIB, "%s not found", lib); 2948 PetscCall(PetscDLLibraryAppend(PETSC_COMM_SELF, &PetscDLLibrariesLoaded, dlib)); 2949 PetscFunctionReturn(PETSC_SUCCESS); 2950 } 2951 2952 /*MC 2953 PCHPDDM - Interface with the HPDDM library. 2954 2955 This `PC` may be used to build multilevel spectral domain decomposition methods based on the GenEO framework {cite}`spillane2011robust` {cite}`al2021multilevel`. 2956 It may be viewed as an alternative to spectral 2957 AMGe or `PCBDDC` with adaptive selection of constraints. The interface is explained in details in {cite}`jolivetromanzampini2020` 2958 2959 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`). 2960 2961 For multilevel preconditioning, when using an assembled or hierarchical Pmat, one must provide an auxiliary local `Mat` (unassembled local operator for GenEO) using 2962 `PCHPDDMSetAuxiliaryMat()`. Calling this routine is not needed when using a `MATIS` Pmat, assembly is done internally using `MatConvert()`. 2963 2964 Options Database Keys: 2965 + -pc_hpddm_define_subdomains <true, default=false> - on the finest level, calls `PCASMSetLocalSubdomains()` with the `IS` supplied in `PCHPDDMSetAuxiliaryMat()` 2966 (not relevant with an unassembled Pmat) 2967 . -pc_hpddm_has_neumann <true, default=false> - on the finest level, informs the `PC` that the local Neumann matrix is supplied in `PCHPDDMSetAuxiliaryMat()` 2968 - -pc_hpddm_coarse_correction <type, default=deflated> - determines the `PCHPDDMCoarseCorrectionType` when calling `PCApply()` 2969 2970 Options for subdomain solvers, subdomain eigensolvers (for computing deflation vectors), and the coarse solver can be set using the following options database prefixes. 2971 .vb 2972 -pc_hpddm_levels_%d_pc_ 2973 -pc_hpddm_levels_%d_ksp_ 2974 -pc_hpddm_levels_%d_eps_ 2975 -pc_hpddm_levels_%d_p 2976 -pc_hpddm_levels_%d_mat_type 2977 -pc_hpddm_coarse_ 2978 -pc_hpddm_coarse_p 2979 -pc_hpddm_coarse_mat_type 2980 -pc_hpddm_coarse_mat_filter 2981 .ve 2982 2983 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 2984 -pc_hpddm_coarse_p 2 -pc_hpddm_coarse_mat_type baij will use 10 deflation vectors per subdomain on the fine "level 1", 2985 aggregate the fine subdomains into 4 "level 2" subdomains, then use 10 deflation vectors per subdomain on "level 2", 2986 and assemble the coarse matrix (of dimension 4 x 10 = 40) on two processes as a `MATBAIJ` (default is `MATSBAIJ`). 2987 2988 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. 2989 2990 Level: intermediate 2991 2992 Notes: 2993 This preconditioner requires that PETSc is built with SLEPc (``--download-slepc``). 2994 2995 By default, the underlying concurrent eigenproblems 2996 are solved using SLEPc shift-and-invert spectral transformation. This is usually what gives the best performance for GenEO, cf. 2997 {cite}`spillane2011robust` {cite}`jolivet2013scalabledd`. As 2998 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 2999 -pc_hpddm_levels_1_st_type sinvert. There are furthermore three options related to the (subdomain-wise local) eigensolver that are not described in 3000 SLEPc documentation since they are specific to `PCHPDDM`. 3001 .vb 3002 -pc_hpddm_levels_1_st_share_sub_ksp 3003 -pc_hpddm_levels_%d_eps_threshold 3004 -pc_hpddm_levels_1_eps_use_inertia 3005 .ve 3006 3007 The first option from the list only applies to the fine-level eigensolver, see `PCHPDDMSetSTShareSubKSP()`. The second option from the list is 3008 used to filter eigenmodes retrieved after convergence of `EPSSolve()` at "level N" such that eigenvectors used to define a "level N+1" coarse 3009 correction are associated to eigenvalues whose magnitude are lower or equal than -pc_hpddm_levels_N_eps_threshold. When using an `EPS` which cannot 3010 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 3011 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 3012 to supply -pc_hpddm_levels_1_eps_nev. This last option also only applies to the fine-level (N = 1) eigensolver. 3013 3014 See also {cite}`dolean2015introduction`, {cite}`al2022robust`, {cite}`al2022robustpd`, and {cite}`nataf2022recent` 3015 3016 .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCHPDDMSetAuxiliaryMat()`, `MATIS`, `PCBDDC`, `PCDEFLATION`, `PCTELESCOPE`, `PCASM`, 3017 `PCHPDDMSetCoarseCorrectionType()`, `PCHPDDMHasNeumannMat()`, `PCHPDDMSetRHSMat()`, `PCHPDDMSetDeflationMat()`, `PCHPDDMSetSTShareSubKSP()`, 3018 `PCHPDDMGetSTShareSubKSP()`, `PCHPDDMGetCoarseCorrectionType()`, `PCHPDDMGetComplexities()` 3019 M*/ 3020 PETSC_EXTERN PetscErrorCode PCCreate_HPDDM(PC pc) 3021 { 3022 PC_HPDDM *data; 3023 PetscBool found; 3024 3025 PetscFunctionBegin; 3026 if (!loadedSym) { 3027 PetscCall(HPDDMLoadDL_Private(&found)); 3028 if (found) PetscCall(PetscDLLibrarySym(PETSC_COMM_SELF, &PetscDLLibrariesLoaded, nullptr, "PCHPDDM_Internal", (void **)&loadedSym)); 3029 } 3030 PetscCheck(loadedSym, PETSC_COMM_SELF, PETSC_ERR_PLIB, "PCHPDDM_Internal symbol not found in loaded libhpddm_petsc"); 3031 PetscCall(PetscNew(&data)); 3032 pc->data = data; 3033 data->Neumann = PETSC_BOOL3_UNKNOWN; 3034 pc->ops->reset = PCReset_HPDDM; 3035 pc->ops->destroy = PCDestroy_HPDDM; 3036 pc->ops->setfromoptions = PCSetFromOptions_HPDDM; 3037 pc->ops->setup = PCSetUp_HPDDM; 3038 pc->ops->apply = PCApply_HPDDM; 3039 pc->ops->matapply = PCMatApply_HPDDM; 3040 pc->ops->view = PCView_HPDDM; 3041 pc->ops->presolve = PCPreSolve_HPDDM; 3042 3043 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetAuxiliaryMat_C", PCHPDDMSetAuxiliaryMat_HPDDM)); 3044 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMHasNeumannMat_C", PCHPDDMHasNeumannMat_HPDDM)); 3045 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetRHSMat_C", PCHPDDMSetRHSMat_HPDDM)); 3046 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetCoarseCorrectionType_C", PCHPDDMSetCoarseCorrectionType_HPDDM)); 3047 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMGetCoarseCorrectionType_C", PCHPDDMGetCoarseCorrectionType_HPDDM)); 3048 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetSTShareSubKSP_C", PCHPDDMSetSTShareSubKSP_HPDDM)); 3049 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMGetSTShareSubKSP_C", PCHPDDMGetSTShareSubKSP_HPDDM)); 3050 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHPDDMSetDeflationMat_C", PCHPDDMSetDeflationMat_HPDDM)); 3051 PetscFunctionReturn(PETSC_SUCCESS); 3052 } 3053 3054 /*@C 3055 PCHPDDMInitializePackage - This function initializes everything in the `PCHPDDM` package. It is called from `PCInitializePackage()`. 3056 3057 Level: developer 3058 3059 .seealso: [](ch_ksp), `PetscInitialize()` 3060 @*/ 3061 PetscErrorCode PCHPDDMInitializePackage(void) 3062 { 3063 char ename[32]; 3064 3065 PetscFunctionBegin; 3066 if (PCHPDDMPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS); 3067 PCHPDDMPackageInitialized = PETSC_TRUE; 3068 PetscCall(PetscRegisterFinalize(PCHPDDMFinalizePackage)); 3069 /* general events registered once during package initialization */ 3070 /* some of these events are not triggered in libpetsc, */ 3071 /* but rather directly in libhpddm_petsc, */ 3072 /* which is in charge of performing the following operations */ 3073 3074 /* domain decomposition structure from Pmat sparsity pattern */ 3075 PetscCall(PetscLogEventRegister("PCHPDDMStrc", PC_CLASSID, &PC_HPDDM_Strc)); 3076 /* Galerkin product, redistribution, and setup (not triggered in libpetsc) */ 3077 PetscCall(PetscLogEventRegister("PCHPDDMPtAP", PC_CLASSID, &PC_HPDDM_PtAP)); 3078 /* Galerkin product with summation, redistribution, and setup (not triggered in libpetsc) */ 3079 PetscCall(PetscLogEventRegister("PCHPDDMPtBP", PC_CLASSID, &PC_HPDDM_PtBP)); 3080 /* next level construction using PtAP and PtBP (not triggered in libpetsc) */ 3081 PetscCall(PetscLogEventRegister("PCHPDDMNext", PC_CLASSID, &PC_HPDDM_Next)); 3082 static_assert(PETSC_PCHPDDM_MAXLEVELS <= 9, "PETSC_PCHPDDM_MAXLEVELS value is too high"); 3083 for (PetscInt i = 1; i < PETSC_PCHPDDM_MAXLEVELS; ++i) { 3084 PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCHPDDMSetUp L%1" PetscInt_FMT, i)); 3085 /* events during a PCSetUp() at level #i _except_ the assembly */ 3086 /* of the Galerkin operator of the coarser level #(i + 1) */ 3087 PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_HPDDM_SetUp[i - 1])); 3088 PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCHPDDMSolve L%1" PetscInt_FMT, i)); 3089 /* events during a PCApply() at level #i _except_ */ 3090 /* the KSPSolve() of the coarser level #(i + 1) */ 3091 PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_HPDDM_Solve[i - 1])); 3092 } 3093 PetscFunctionReturn(PETSC_SUCCESS); 3094 } 3095 3096 /*@C 3097 PCHPDDMFinalizePackage - This function frees everything from the `PCHPDDM` package. It is called from `PetscFinalize()`. 3098 3099 Level: developer 3100 3101 .seealso: [](ch_ksp), `PetscFinalize()` 3102 @*/ 3103 PetscErrorCode PCHPDDMFinalizePackage(void) 3104 { 3105 PetscFunctionBegin; 3106 PCHPDDMPackageInitialized = PETSC_FALSE; 3107 PetscFunctionReturn(PETSC_SUCCESS); 3108 } 3109 3110 static PetscErrorCode MatMult_Harmonic(Mat A, Vec x, Vec y) 3111 { 3112 Harmonic h; /* [ A_00 A_01 ], furthermore, A_00 = [ A_loc,loc A_loc,ovl ], thus, A_01 = [ ] */ 3113 /* [ A_10 A_11 A_12 ] [ A_ovl,loc A_ovl,ovl ] [ A_ovl,1 ] */ 3114 Vec sub; /* y = A x = R_loc R_0 [ A_00 A_01 ]^-1 R_loc = [ I_loc ] */ 3115 /* [ A_10 A_11 ] R_1^T A_12 x [ ] */ 3116 PetscFunctionBegin; 3117 PetscCall(MatShellGetContext(A, &h)); 3118 PetscCall(VecSet(h->v, 0.0)); 3119 PetscCall(VecGetSubVector(h->v, h->is[0], &sub)); 3120 PetscCall(MatMult(h->A[0], x, sub)); 3121 PetscCall(VecRestoreSubVector(h->v, h->is[0], &sub)); 3122 PetscCall(KSPSolve(h->ksp, h->v, h->v)); 3123 PetscCall(VecISCopy(h->v, h->is[1], SCATTER_REVERSE, y)); 3124 PetscFunctionReturn(PETSC_SUCCESS); 3125 } 3126 3127 static PetscErrorCode MatMultTranspose_Harmonic(Mat A, Vec y, Vec x) 3128 { 3129 Harmonic h; /* x = A^T y = [ A_00 A_01 ]^-T R_0^T R_loc^T y */ 3130 Vec sub; /* A_12^T R_1 [ A_10 A_11 ] */ 3131 3132 PetscFunctionBegin; 3133 PetscCall(MatShellGetContext(A, &h)); 3134 PetscCall(VecSet(h->v, 0.0)); 3135 PetscCall(VecISCopy(h->v, h->is[1], SCATTER_FORWARD, y)); 3136 PetscCall(KSPSolveTranspose(h->ksp, h->v, h->v)); 3137 PetscCall(VecGetSubVector(h->v, h->is[0], &sub)); 3138 PetscCall(MatMultTranspose(h->A[0], sub, x)); 3139 PetscCall(VecRestoreSubVector(h->v, h->is[0], &sub)); 3140 PetscFunctionReturn(PETSC_SUCCESS); 3141 } 3142 3143 static PetscErrorCode MatProduct_AB_Harmonic(Mat S, Mat X, Mat Y, void *) 3144 { 3145 Harmonic h; 3146 Mat A, B; 3147 Vec a, b; 3148 3149 PetscFunctionBegin; 3150 PetscCall(MatShellGetContext(S, &h)); 3151 PetscCall(MatMatMult(h->A[0], X, MAT_INITIAL_MATRIX, PETSC_CURRENT, &A)); 3152 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, h->ksp->pc->mat->rmap->n, A->cmap->n, nullptr, &B)); 3153 for (PetscInt i = 0; i < A->cmap->n; ++i) { 3154 PetscCall(MatDenseGetColumnVecRead(A, i, &a)); 3155 PetscCall(MatDenseGetColumnVecWrite(B, i, &b)); 3156 PetscCall(VecISCopy(b, h->is[0], SCATTER_FORWARD, a)); 3157 PetscCall(MatDenseRestoreColumnVecWrite(B, i, &b)); 3158 PetscCall(MatDenseRestoreColumnVecRead(A, i, &a)); 3159 } 3160 PetscCall(MatDestroy(&A)); 3161 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, h->ksp->pc->mat->rmap->n, B->cmap->n, nullptr, &A)); 3162 PetscCall(KSPMatSolve(h->ksp, B, A)); 3163 PetscCall(MatDestroy(&B)); 3164 for (PetscInt i = 0; i < A->cmap->n; ++i) { 3165 PetscCall(MatDenseGetColumnVecRead(A, i, &a)); 3166 PetscCall(MatDenseGetColumnVecWrite(Y, i, &b)); 3167 PetscCall(VecISCopy(a, h->is[1], SCATTER_REVERSE, b)); 3168 PetscCall(MatDenseRestoreColumnVecWrite(Y, i, &b)); 3169 PetscCall(MatDenseRestoreColumnVecRead(A, i, &a)); 3170 } 3171 PetscCall(MatDestroy(&A)); 3172 PetscFunctionReturn(PETSC_SUCCESS); 3173 } 3174 3175 static PetscErrorCode MatDestroy_Harmonic(Mat A) 3176 { 3177 Harmonic h; 3178 3179 PetscFunctionBegin; 3180 PetscCall(MatShellGetContext(A, &h)); 3181 for (PetscInt i = 0; i < 5; ++i) PetscCall(ISDestroy(h->is + i)); 3182 PetscCall(PetscFree(h->is)); 3183 PetscCall(VecDestroy(&h->v)); 3184 for (PetscInt i = 0; i < 2; ++i) PetscCall(MatDestroy(h->A + i)); 3185 PetscCall(PetscFree(h->A)); 3186 PetscCall(KSPDestroy(&h->ksp)); 3187 PetscCall(PetscFree(h)); 3188 PetscFunctionReturn(PETSC_SUCCESS); 3189 } 3190