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