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