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