1 #include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/ 2 #include <petsc/private/kspimpl.h> /* This is needed to provide the appropriate PETSC_EXTERN for KSP_Solve_FS ....*/ 3 #include <petsc/private/matimpl.h> /* MatScatterDense_Private() for PCMatApply() */ 4 #include <petscdm.h> 5 #include <petscdevice.h> 6 #if PetscDefined(HAVE_CUDA) 7 #include <petscdevice_cuda.h> 8 #endif 9 #if PetscDefined(HAVE_HIP) 10 #include <petscdevice_hip.h> 11 #endif 12 13 const char *const PCFieldSplitSchurPreTypes[] = {"SELF", "SELFP", "A11", "USER", "FULL", "PCFieldSplitSchurPreType", "PC_FIELDSPLIT_SCHUR_PRE_", NULL}; 14 const char *const PCFieldSplitSchurFactTypes[] = {"DIAG", "LOWER", "UPPER", "FULL", "PCFieldSplitSchurFactType", "PC_FIELDSPLIT_SCHUR_FACT_", NULL}; 15 16 PetscLogEvent KSP_Solve_FS_0, KSP_Solve_FS_1, KSP_Solve_FS_S, KSP_Solve_FS_U, KSP_Solve_FS_L, KSP_Solve_FS_2, KSP_Solve_FS_3, KSP_Solve_FS_4; 17 18 typedef struct _PC_FieldSplitLink *PC_FieldSplitLink; 19 struct _PC_FieldSplitLink { 20 KSP ksp; 21 Vec x, y, z; 22 Mat X, Y, Z; 23 char *splitname; 24 PetscInt nfields; 25 PetscInt *fields, *fields_col; 26 VecScatter sctx; 27 IS is, is_col; 28 PC_FieldSplitLink next, previous; 29 PetscLogEvent event; 30 31 /* Used only when setting coordinates with PCSetCoordinates */ 32 PetscInt dim; 33 PetscInt ndofs; 34 PetscReal *coords; 35 }; 36 37 typedef struct { 38 PCCompositeType type; 39 PetscBool defaultsplit; /* Flag for a system with a set of 'k' scalar fields with the same layout (and bs = k) */ 40 PetscBool splitdefined; /* Flag is set after the splits have been defined, to prevent more splits from being added */ 41 PetscInt bs; /* Block size for IS and Mat structures */ 42 PetscInt nsplits; /* Number of field divisions defined */ 43 Vec *x, *y, w1, w2; 44 Mat *mat; /* The diagonal block for each split */ 45 Mat *pmat; /* The preconditioning diagonal block for each split */ 46 Mat *Afield; /* The rows of the matrix associated with each split */ 47 PetscBool issetup; 48 49 /* Only used when Schur complement preconditioning is used */ 50 Mat B; /* The (0,1) block */ 51 Mat C; /* The (1,0) block */ 52 Mat schur; /* The Schur complement S = A11 - A10 A00^{-1} A01, the KSP here, kspinner, is H_1 in [El08] */ 53 Mat schurp; /* Assembled approximation to S built by MatSchurComplement to be used as a matrix for constructing the preconditioner when solving with S */ 54 Mat schur_user; /* User-provided matrix for constructing the preconditioner for the Schur complement */ 55 PCFieldSplitSchurPreType schurpre; /* Determines which matrix is used for the Schur complement */ 56 PCFieldSplitSchurFactType schurfactorization; 57 KSP kspschur; /* The solver for S */ 58 KSP kspupper; /* The solver for A in the upper diagonal part of the factorization (H_2 in [El08]) */ 59 PetscScalar schurscale; /* Scaling factor for the Schur complement solution with DIAG factorization */ 60 61 /* Only used when Golub-Kahan bidiagonalization preconditioning is used */ 62 Mat H; /* The modified matrix H = A00 + nu*A01*A01' */ 63 PetscReal gkbtol; /* Stopping tolerance for lower bound estimate */ 64 PetscInt gkbdelay; /* The delay window for the stopping criterion */ 65 PetscReal gkbnu; /* Parameter for augmented Lagrangian H = A + nu*A01*A01' */ 66 PetscInt gkbmaxit; /* Maximum number of iterations for outer loop */ 67 PetscBool gkbmonitor; /* Monitor for gkb iterations and the lower bound error */ 68 PetscViewer gkbviewer; /* Viewer context for gkbmonitor */ 69 Vec u, v, d, Hu; /* Work vectors for the GKB algorithm */ 70 PetscScalar *vecz; /* Contains intermediate values, eg for lower bound */ 71 72 PC_FieldSplitLink head; 73 PetscBool isrestrict; /* indicates PCFieldSplitRestrictIS() has been last called on this object, hack */ 74 PetscBool suboptionsset; /* Indicates that the KSPSetFromOptions() has been called on the sub-KSPs */ 75 PetscBool dm_splits; /* Whether to use DMCreateFieldDecomposition() whenever possible */ 76 PetscBool diag_use_amat; /* Whether to extract diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */ 77 PetscBool offdiag_use_amat; /* Whether to extract off-diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */ 78 PetscBool detect; /* Whether to form 2-way split by finding zero diagonal entries */ 79 PetscBool coordinates_set; /* Whether PCSetCoordinates has been called */ 80 } PC_FieldSplit; 81 82 /* 83 Note: 84 there is no particular reason that pmat, x, and y are stored as arrays in PC_FieldSplit instead of 85 inside PC_FieldSplitLink, just historical. If you want to be able to add new fields after already using the 86 PC you could change this. 87 */ 88 89 /* This helper is so that setting a user-provided matrix is orthogonal to choosing to use it. This way the 90 * application-provided FormJacobian can provide this matrix without interfering with the user's (command-line) choices. */ 91 static Mat FieldSplitSchurPre(PC_FieldSplit *jac) 92 { 93 switch (jac->schurpre) { 94 case PC_FIELDSPLIT_SCHUR_PRE_SELF: 95 return jac->schur; 96 case PC_FIELDSPLIT_SCHUR_PRE_SELFP: 97 return jac->schurp; 98 case PC_FIELDSPLIT_SCHUR_PRE_A11: 99 return jac->pmat[1]; 100 case PC_FIELDSPLIT_SCHUR_PRE_FULL: /* We calculate this and store it in schur_user */ 101 case PC_FIELDSPLIT_SCHUR_PRE_USER: /* Use a user-provided matrix if it is given, otherwise diagonal block */ 102 default: 103 return jac->schur_user ? jac->schur_user : jac->pmat[1]; 104 } 105 } 106 107 #include <petscdraw.h> 108 static PetscErrorCode PCView_FieldSplit(PC pc, PetscViewer viewer) 109 { 110 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 111 PetscBool isascii, isdraw; 112 PetscInt i, j; 113 PC_FieldSplitLink ilink = jac->head; 114 115 PetscFunctionBegin; 116 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); 117 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw)); 118 if (isascii) { 119 if (jac->bs > 0) { 120 PetscCall(PetscViewerASCIIPrintf(viewer, " FieldSplit with %s composition: total splits = %" PetscInt_FMT ", blocksize = %" PetscInt_FMT "\n", PCCompositeTypes[jac->type], jac->nsplits, jac->bs)); 121 } else { 122 PetscCall(PetscViewerASCIIPrintf(viewer, " FieldSplit with %s composition: total splits = %" PetscInt_FMT "\n", PCCompositeTypes[jac->type], jac->nsplits)); 123 } 124 if (pc->useAmat) PetscCall(PetscViewerASCIIPrintf(viewer, " using Amat (not Pmat) as operator for blocks\n")); 125 if (jac->diag_use_amat) PetscCall(PetscViewerASCIIPrintf(viewer, " using Amat (not Pmat) as operator for diagonal blocks\n")); 126 if (jac->offdiag_use_amat) PetscCall(PetscViewerASCIIPrintf(viewer, " using Amat (not Pmat) as operator for off-diagonal blocks\n")); 127 PetscCall(PetscViewerASCIIPrintf(viewer, " Solver info for each split is in the following KSP objects:\n")); 128 for (i = 0; i < jac->nsplits; i++) { 129 if (ilink->fields) { 130 PetscCall(PetscViewerASCIIPrintf(viewer, "Split number %" PetscInt_FMT " Fields ", i)); 131 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE)); 132 for (j = 0; j < ilink->nfields; j++) { 133 if (j > 0) PetscCall(PetscViewerASCIIPrintf(viewer, ",")); 134 PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, ilink->fields[j])); 135 } 136 PetscCall(PetscViewerASCIIPrintf(viewer, "\n")); 137 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE)); 138 } else { 139 PetscCall(PetscViewerASCIIPrintf(viewer, "Split number %" PetscInt_FMT " Defined by IS\n", i)); 140 } 141 PetscCall(KSPView(ilink->ksp, viewer)); 142 ilink = ilink->next; 143 } 144 } 145 146 if (isdraw) { 147 PetscDraw draw; 148 PetscReal x, y, w, wd; 149 150 PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw)); 151 PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y)); 152 w = 2 * PetscMin(1.0 - x, x); 153 wd = w / (jac->nsplits + 1); 154 x = x - wd * (jac->nsplits - 1) / 2.0; 155 for (i = 0; i < jac->nsplits; i++) { 156 PetscCall(PetscDrawPushCurrentPoint(draw, x, y)); 157 PetscCall(KSPView(ilink->ksp, viewer)); 158 PetscCall(PetscDrawPopCurrentPoint(draw)); 159 x += wd; 160 ilink = ilink->next; 161 } 162 } 163 PetscFunctionReturn(PETSC_SUCCESS); 164 } 165 166 static PetscErrorCode PCView_FieldSplit_Schur(PC pc, PetscViewer viewer) 167 { 168 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 169 PetscBool isascii, isdraw; 170 PetscInt i, j; 171 PC_FieldSplitLink ilink = jac->head; 172 MatSchurComplementAinvType atype; 173 174 PetscFunctionBegin; 175 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); 176 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw)); 177 if (isascii) { 178 if (jac->bs > 0) { 179 PetscCall(PetscViewerASCIIPrintf(viewer, " FieldSplit with Schur preconditioner, blocksize = %" PetscInt_FMT ", factorization %s\n", jac->bs, PCFieldSplitSchurFactTypes[jac->schurfactorization])); 180 } else { 181 PetscCall(PetscViewerASCIIPrintf(viewer, " FieldSplit with Schur preconditioner, factorization %s\n", PCFieldSplitSchurFactTypes[jac->schurfactorization])); 182 } 183 if (pc->useAmat) PetscCall(PetscViewerASCIIPrintf(viewer, " using Amat (not Pmat) as operator for blocks\n")); 184 switch (jac->schurpre) { 185 case PC_FIELDSPLIT_SCHUR_PRE_SELF: 186 PetscCall(PetscViewerASCIIPrintf(viewer, " Preconditioner for the Schur complement formed from S itself\n")); 187 break; 188 case PC_FIELDSPLIT_SCHUR_PRE_SELFP: 189 if (jac->schur) { 190 PetscCall(MatSchurComplementGetAinvType(jac->schur, &atype)); 191 PetscCall(PetscViewerASCIIPrintf(viewer, " Preconditioner for the Schur complement formed from Sp, an assembled approximation to S, which uses A00's %sinverse\n", atype == MAT_SCHUR_COMPLEMENT_AINV_DIAG ? "diagonal's " : (atype == MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG ? "block diagonal's " : (atype == MAT_SCHUR_COMPLEMENT_AINV_FULL ? "full " : "lumped diagonal's ")))); 192 } 193 break; 194 case PC_FIELDSPLIT_SCHUR_PRE_A11: 195 PetscCall(PetscViewerASCIIPrintf(viewer, " Preconditioner for the Schur complement formed from A11\n")); 196 break; 197 case PC_FIELDSPLIT_SCHUR_PRE_FULL: 198 PetscCall(PetscViewerASCIIPrintf(viewer, " Preconditioner for the Schur complement formed from the exact Schur complement\n")); 199 break; 200 case PC_FIELDSPLIT_SCHUR_PRE_USER: 201 if (jac->schur_user) { 202 PetscCall(PetscViewerASCIIPrintf(viewer, " Preconditioner for the Schur complement formed from user provided matrix\n")); 203 } else { 204 PetscCall(PetscViewerASCIIPrintf(viewer, " Preconditioner for the Schur complement formed from A11\n")); 205 } 206 break; 207 default: 208 SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Invalid Schur preconditioning type: %d", jac->schurpre); 209 } 210 PetscCall(PetscViewerASCIIPrintf(viewer, " Split info:\n")); 211 PetscCall(PetscViewerASCIIPushTab(viewer)); 212 for (i = 0; i < jac->nsplits; i++) { 213 if (ilink->fields) { 214 PetscCall(PetscViewerASCIIPrintf(viewer, "Split number %" PetscInt_FMT " Fields ", i)); 215 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE)); 216 for (j = 0; j < ilink->nfields; j++) { 217 if (j > 0) PetscCall(PetscViewerASCIIPrintf(viewer, ",")); 218 PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, ilink->fields[j])); 219 } 220 PetscCall(PetscViewerASCIIPrintf(viewer, "\n")); 221 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE)); 222 } else { 223 PetscCall(PetscViewerASCIIPrintf(viewer, "Split number %" PetscInt_FMT " Defined by IS\n", i)); 224 } 225 ilink = ilink->next; 226 } 227 PetscCall(PetscViewerASCIIPrintf(viewer, "KSP solver for A00 block\n")); 228 PetscCall(PetscViewerASCIIPushTab(viewer)); 229 if (jac->head) PetscCall(KSPView(jac->head->ksp, viewer)); 230 else PetscCall(PetscViewerASCIIPrintf(viewer, " not yet available\n")); 231 PetscCall(PetscViewerASCIIPopTab(viewer)); 232 if (jac->head && jac->kspupper != jac->head->ksp) { 233 PetscCall(PetscViewerASCIIPrintf(viewer, "KSP solver for upper A00 in upper triangular factor\n")); 234 PetscCall(PetscViewerASCIIPushTab(viewer)); 235 if (jac->kspupper) PetscCall(KSPView(jac->kspupper, viewer)); 236 else PetscCall(PetscViewerASCIIPrintf(viewer, " not yet available\n")); 237 PetscCall(PetscViewerASCIIPopTab(viewer)); 238 } 239 PetscCall(PetscViewerASCIIPrintf(viewer, "KSP solver for S = A11 - A10 inv(A00) A01\n")); 240 PetscCall(PetscViewerASCIIPushTab(viewer)); 241 if (jac->kspschur) { 242 PetscCall(KSPView(jac->kspschur, viewer)); 243 } else { 244 PetscCall(PetscViewerASCIIPrintf(viewer, " not yet available\n")); 245 } 246 PetscCall(PetscViewerASCIIPopTab(viewer)); 247 PetscCall(PetscViewerASCIIPopTab(viewer)); 248 } else if (isdraw && jac->head) { 249 PetscDraw draw; 250 PetscReal x, y, w, wd, h; 251 PetscInt cnt = 2; 252 char str[32]; 253 254 PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw)); 255 PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y)); 256 if (jac->kspupper != jac->head->ksp) cnt++; 257 w = 2 * PetscMin(1.0 - x, x); 258 wd = w / (cnt + 1); 259 260 PetscCall(PetscSNPrintf(str, 32, "Schur fact. %s", PCFieldSplitSchurFactTypes[jac->schurfactorization])); 261 PetscCall(PetscDrawStringBoxed(draw, x, y, PETSC_DRAW_RED, PETSC_DRAW_BLACK, str, NULL, &h)); 262 y -= h; 263 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_USER && !jac->schur_user) { 264 PetscCall(PetscSNPrintf(str, 32, "Prec. for Schur from %s", PCFieldSplitSchurPreTypes[PC_FIELDSPLIT_SCHUR_PRE_A11])); 265 } else { 266 PetscCall(PetscSNPrintf(str, 32, "Prec. for Schur from %s", PCFieldSplitSchurPreTypes[jac->schurpre])); 267 } 268 PetscCall(PetscDrawStringBoxed(draw, x + wd * (cnt - 1) / 2.0, y, PETSC_DRAW_RED, PETSC_DRAW_BLACK, str, NULL, &h)); 269 y -= h; 270 x = x - wd * (cnt - 1) / 2.0; 271 272 PetscCall(PetscDrawPushCurrentPoint(draw, x, y)); 273 PetscCall(KSPView(jac->head->ksp, viewer)); 274 PetscCall(PetscDrawPopCurrentPoint(draw)); 275 if (jac->kspupper != jac->head->ksp) { 276 x += wd; 277 PetscCall(PetscDrawPushCurrentPoint(draw, x, y)); 278 PetscCall(KSPView(jac->kspupper, viewer)); 279 PetscCall(PetscDrawPopCurrentPoint(draw)); 280 } 281 x += wd; 282 PetscCall(PetscDrawPushCurrentPoint(draw, x, y)); 283 PetscCall(KSPView(jac->kspschur, viewer)); 284 PetscCall(PetscDrawPopCurrentPoint(draw)); 285 } 286 PetscFunctionReturn(PETSC_SUCCESS); 287 } 288 289 static PetscErrorCode PCView_FieldSplit_GKB(PC pc, PetscViewer viewer) 290 { 291 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 292 PetscBool isascii, isdraw; 293 PetscInt i, j; 294 PC_FieldSplitLink ilink = jac->head; 295 296 PetscFunctionBegin; 297 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); 298 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw)); 299 if (isascii) { 300 if (jac->bs > 0) { 301 PetscCall(PetscViewerASCIIPrintf(viewer, " FieldSplit with %s composition: total splits = %" PetscInt_FMT ", blocksize = %" PetscInt_FMT "\n", PCCompositeTypes[jac->type], jac->nsplits, jac->bs)); 302 } else { 303 PetscCall(PetscViewerASCIIPrintf(viewer, " FieldSplit with %s composition: total splits = %" PetscInt_FMT "\n", PCCompositeTypes[jac->type], jac->nsplits)); 304 } 305 if (pc->useAmat) PetscCall(PetscViewerASCIIPrintf(viewer, " using Amat (not Pmat) as operator for blocks\n")); 306 if (jac->diag_use_amat) PetscCall(PetscViewerASCIIPrintf(viewer, " using Amat (not Pmat) as operator for diagonal blocks\n")); 307 if (jac->offdiag_use_amat) PetscCall(PetscViewerASCIIPrintf(viewer, " using Amat (not Pmat) as operator for off-diagonal blocks\n")); 308 309 PetscCall(PetscViewerASCIIPrintf(viewer, " Stopping tolerance=%.1e, delay in error estimate=%" PetscInt_FMT ", maximum iterations=%" PetscInt_FMT "\n", (double)jac->gkbtol, jac->gkbdelay, jac->gkbmaxit)); 310 PetscCall(PetscViewerASCIIPrintf(viewer, " Solver info for H = A00 + nu*A01*A01' matrix:\n")); 311 PetscCall(PetscViewerASCIIPushTab(viewer)); 312 313 if (ilink->fields) { 314 PetscCall(PetscViewerASCIIPrintf(viewer, "Split number 0 Fields ")); 315 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE)); 316 for (j = 0; j < ilink->nfields; j++) { 317 if (j > 0) PetscCall(PetscViewerASCIIPrintf(viewer, ",")); 318 PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, ilink->fields[j])); 319 } 320 PetscCall(PetscViewerASCIIPrintf(viewer, "\n")); 321 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE)); 322 } else { 323 PetscCall(PetscViewerASCIIPrintf(viewer, "Split number 0 Defined by IS\n")); 324 } 325 PetscCall(KSPView(ilink->ksp, viewer)); 326 327 PetscCall(PetscViewerASCIIPopTab(viewer)); 328 } 329 330 if (isdraw) { 331 PetscDraw draw; 332 PetscReal x, y, w, wd; 333 334 PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw)); 335 PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y)); 336 w = 2 * PetscMin(1.0 - x, x); 337 wd = w / (jac->nsplits + 1); 338 x = x - wd * (jac->nsplits - 1) / 2.0; 339 for (i = 0; i < jac->nsplits; i++) { 340 PetscCall(PetscDrawPushCurrentPoint(draw, x, y)); 341 PetscCall(KSPView(ilink->ksp, viewer)); 342 PetscCall(PetscDrawPopCurrentPoint(draw)); 343 x += wd; 344 ilink = ilink->next; 345 } 346 } 347 PetscFunctionReturn(PETSC_SUCCESS); 348 } 349 350 /* Precondition: jac->bs is set to a meaningful value or MATNEST */ 351 static PetscErrorCode PCFieldSplitSetRuntimeSplits_Private(PC pc) 352 { 353 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 354 PetscInt bs, i, nfields, *ifields, nfields_col, *ifields_col; 355 PetscBool flg, flg_col, mnest; 356 char optionname[128], splitname[8], optionname_col[128]; 357 358 PetscFunctionBegin; 359 PetscCall(PetscObjectTypeCompare((PetscObject)pc->mat, MATNEST, &mnest)); 360 if (mnest) { 361 PetscCall(MatNestGetSize(pc->pmat, &bs, NULL)); 362 } else { 363 bs = jac->bs; 364 } 365 PetscCall(PetscMalloc2(bs, &ifields, bs, &ifields_col)); 366 for (i = 0, flg = PETSC_TRUE;; i++) { 367 PetscCall(PetscSNPrintf(splitname, sizeof(splitname), "%" PetscInt_FMT, i)); 368 PetscCall(PetscSNPrintf(optionname, sizeof(optionname), "-pc_fieldsplit_%" PetscInt_FMT "_fields", i)); 369 PetscCall(PetscSNPrintf(optionname_col, sizeof(optionname_col), "-pc_fieldsplit_%" PetscInt_FMT "_fields_col", i)); 370 nfields = bs; 371 nfields_col = bs; 372 PetscCall(PetscOptionsGetIntArray(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, optionname, ifields, &nfields, &flg)); 373 PetscCall(PetscOptionsGetIntArray(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, optionname_col, ifields_col, &nfields_col, &flg_col)); 374 if (!flg) break; 375 else if (flg && !flg_col) { 376 PetscCheck(nfields, PETSC_COMM_SELF, PETSC_ERR_USER, "Cannot list zero fields"); 377 PetscCall(PCFieldSplitSetFields(pc, splitname, nfields, ifields, ifields)); 378 } else { 379 PetscCheck(nfields && nfields_col, PETSC_COMM_SELF, PETSC_ERR_USER, "Cannot list zero fields"); 380 PetscCheck(nfields == nfields_col, PETSC_COMM_SELF, PETSC_ERR_USER, "Number of row and column fields must match"); 381 PetscCall(PCFieldSplitSetFields(pc, splitname, nfields, ifields, ifields_col)); 382 } 383 } 384 if (i > 0) { 385 /* Makes command-line setting of splits take precedence over setting them in code. 386 Otherwise subsequent calls to PCFieldSplitSetIS() or PCFieldSplitSetFields() would 387 create new splits, which would probably not be what the user wanted. */ 388 jac->splitdefined = PETSC_TRUE; 389 } 390 PetscCall(PetscFree2(ifields, ifields_col)); 391 PetscFunctionReturn(PETSC_SUCCESS); 392 } 393 394 static PetscErrorCode PCFieldSplitSetDefaults(PC pc) 395 { 396 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 397 PC_FieldSplitLink ilink = jac->head; 398 PetscBool fieldsplit_default = PETSC_FALSE, coupling = PETSC_FALSE; 399 PetscInt i; 400 401 PetscFunctionBegin; 402 /* 403 Kinda messy, but at least this now uses DMCreateFieldDecomposition(). 404 Should probably be rewritten. 405 */ 406 if (!ilink) { 407 PetscCall(PetscOptionsGetBool(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_fieldsplit_detect_coupling", &coupling, NULL)); 408 if (pc->dm && jac->dm_splits && !jac->detect && !coupling) { 409 PetscInt numFields, f, i, j; 410 char **fieldNames; 411 IS *fields; 412 DM *dms; 413 DM subdm[128]; 414 PetscBool flg; 415 416 PetscCall(DMCreateFieldDecomposition(pc->dm, &numFields, &fieldNames, &fields, &dms)); 417 /* Allow the user to prescribe the splits */ 418 for (i = 0, flg = PETSC_TRUE;; i++) { 419 PetscInt ifields[128]; 420 IS compField; 421 char optionname[128], splitname[8]; 422 PetscInt nfields = numFields; 423 424 PetscCall(PetscSNPrintf(optionname, sizeof(optionname), "-pc_fieldsplit_%" PetscInt_FMT "_fields", i)); 425 PetscCall(PetscOptionsGetIntArray(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, optionname, ifields, &nfields, &flg)); 426 if (!flg) break; 427 PetscCheck(numFields <= 128, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Cannot currently support %" PetscInt_FMT " > 128 fields", numFields); 428 PetscCall(DMCreateSubDM(pc->dm, nfields, ifields, &compField, &subdm[i])); 429 if (nfields == 1) { 430 PetscCall(PCFieldSplitSetIS(pc, fieldNames[ifields[0]], compField)); 431 } else { 432 PetscCall(PetscSNPrintf(splitname, sizeof(splitname), "%" PetscInt_FMT, i)); 433 PetscCall(PCFieldSplitSetIS(pc, splitname, compField)); 434 } 435 PetscCall(ISDestroy(&compField)); 436 for (j = 0; j < nfields; ++j) { 437 f = ifields[j]; 438 PetscCall(PetscFree(fieldNames[f])); 439 PetscCall(ISDestroy(&fields[f])); 440 } 441 } 442 if (i == 0) { 443 for (f = 0; f < numFields; ++f) { 444 PetscCall(PCFieldSplitSetIS(pc, fieldNames[f], fields[f])); 445 PetscCall(PetscFree(fieldNames[f])); 446 PetscCall(ISDestroy(&fields[f])); 447 } 448 } else { 449 for (j = 0; j < numFields; j++) PetscCall(DMDestroy(dms + j)); 450 PetscCall(PetscFree(dms)); 451 PetscCall(PetscMalloc1(i, &dms)); 452 for (j = 0; j < i; ++j) dms[j] = subdm[j]; 453 } 454 PetscCall(PetscFree(fieldNames)); 455 PetscCall(PetscFree(fields)); 456 if (dms) { 457 PetscCall(PetscInfo(pc, "Setting up physics based fieldsplit preconditioner using the embedded DM\n")); 458 for (ilink = jac->head, i = 0; ilink; ilink = ilink->next, ++i) { 459 const char *prefix; 460 PetscCall(PetscObjectGetOptionsPrefix((PetscObject)ilink->ksp, &prefix)); 461 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)dms[i], prefix)); 462 PetscCall(KSPSetDM(ilink->ksp, dms[i])); 463 PetscCall(KSPSetDMActive(ilink->ksp, KSP_DMACTIVE_ALL, PETSC_FALSE)); 464 PetscCall(PetscObjectIncrementTabLevel((PetscObject)dms[i], (PetscObject)ilink->ksp, 0)); 465 PetscCall(DMDestroy(&dms[i])); 466 } 467 PetscCall(PetscFree(dms)); 468 } 469 } else { 470 if (jac->bs <= 0) { 471 if (pc->pmat) PetscCall(MatGetBlockSize(pc->pmat, &jac->bs)); 472 else jac->bs = 1; 473 } 474 475 if (jac->detect) { 476 IS zerodiags, rest; 477 PetscInt nmin, nmax; 478 479 PetscCall(MatGetOwnershipRange(pc->mat, &nmin, &nmax)); 480 if (jac->diag_use_amat) { 481 PetscCall(MatFindZeroDiagonals(pc->mat, &zerodiags)); 482 } else { 483 PetscCall(MatFindZeroDiagonals(pc->pmat, &zerodiags)); 484 } 485 PetscCall(ISComplement(zerodiags, nmin, nmax, &rest)); 486 PetscCall(PCFieldSplitSetIS(pc, "0", rest)); 487 PetscCall(PCFieldSplitSetIS(pc, "1", zerodiags)); 488 PetscCall(ISDestroy(&zerodiags)); 489 PetscCall(ISDestroy(&rest)); 490 } else if (coupling) { 491 IS coupling, rest; 492 PetscInt nmin, nmax; 493 494 PetscCall(MatGetOwnershipRange(pc->mat, &nmin, &nmax)); 495 if (jac->offdiag_use_amat) { 496 PetscCall(MatFindOffBlockDiagonalEntries(pc->mat, &coupling)); 497 } else { 498 PetscCall(MatFindOffBlockDiagonalEntries(pc->pmat, &coupling)); 499 } 500 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)pc->mat), nmax - nmin, nmin, 1, &rest)); 501 PetscCall(ISSetIdentity(rest)); 502 PetscCall(PCFieldSplitSetIS(pc, "0", rest)); 503 PetscCall(PCFieldSplitSetIS(pc, "1", coupling)); 504 PetscCall(ISDestroy(&coupling)); 505 PetscCall(ISDestroy(&rest)); 506 } else { 507 PetscCall(PetscOptionsGetBool(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_fieldsplit_default", &fieldsplit_default, NULL)); 508 if (!fieldsplit_default) { 509 /* Allow user to set fields from command line, if bs was known at the time of PCSetFromOptions_FieldSplit() 510 then it is set there. This is not ideal because we should only have options set in XXSetFromOptions(). */ 511 PetscCall(PCFieldSplitSetRuntimeSplits_Private(pc)); 512 if (jac->splitdefined) PetscCall(PetscInfo(pc, "Splits defined using the options database\n")); 513 } 514 if ((fieldsplit_default || !jac->splitdefined) && !jac->isrestrict) { 515 Mat M = pc->pmat; 516 PetscBool isnest; 517 PetscInt nf; 518 519 PetscCall(PetscInfo(pc, "Using default splitting of fields\n")); 520 PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATNEST, &isnest)); 521 if (!isnest) { 522 M = pc->mat; 523 PetscCall(PetscObjectTypeCompare((PetscObject)pc->mat, MATNEST, &isnest)); 524 } 525 if (!isnest) nf = jac->bs; 526 else PetscCall(MatNestGetSize(M, &nf, NULL)); 527 for (i = 0; i < nf; i++) { 528 char splitname[8]; 529 530 PetscCall(PetscSNPrintf(splitname, sizeof(splitname), "%" PetscInt_FMT, i)); 531 PetscCall(PCFieldSplitSetFields(pc, splitname, 1, &i, &i)); 532 } 533 jac->defaultsplit = PETSC_TRUE; 534 } 535 } 536 } 537 } else if (jac->nsplits == 1) { 538 IS is2; 539 PetscInt nmin, nmax; 540 541 PetscCheck(ilink->is, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Must provide at least two sets of fields to PCFieldSplit()"); 542 PetscCall(MatGetOwnershipRange(pc->mat, &nmin, &nmax)); 543 PetscCall(ISComplement(ilink->is, nmin, nmax, &is2)); 544 PetscCall(PCFieldSplitSetIS(pc, "1", is2)); 545 PetscCall(ISDestroy(&is2)); 546 } 547 548 PetscCheck(jac->nsplits >= 2, PetscObjectComm((PetscObject)pc), PETSC_ERR_PLIB, "Unhandled case, must have at least two fields, not %" PetscInt_FMT, jac->nsplits); 549 PetscFunctionReturn(PETSC_SUCCESS); 550 } 551 552 static PetscErrorCode MatGolubKahanComputeExplicitOperator(Mat A, Mat B, Mat C, Mat *H, PetscReal gkbnu) 553 { 554 Mat BT, T; 555 PetscReal nrmT, nrmB; 556 557 PetscFunctionBegin; 558 PetscCall(MatHermitianTranspose(C, MAT_INITIAL_MATRIX, &T)); /* Test if augmented matrix is symmetric */ 559 PetscCall(MatAXPY(T, -1.0, B, DIFFERENT_NONZERO_PATTERN)); 560 PetscCall(MatNorm(T, NORM_1, &nrmT)); 561 PetscCall(MatNorm(B, NORM_1, &nrmB)); 562 PetscCheck(nrmB <= 0 || nrmT / nrmB < PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Matrix is not symmetric/Hermitian, GKB is not applicable."); 563 564 /* Compute augmented Lagrangian matrix H = A00 + nu*A01*A01'. This corresponds to */ 565 /* setting N := 1/nu*I in [Ar13]. */ 566 PetscCall(MatHermitianTranspose(B, MAT_INITIAL_MATRIX, &BT)); 567 PetscCall(MatMatMult(B, BT, MAT_INITIAL_MATRIX, PETSC_CURRENT, H)); /* H = A01*A01' */ 568 PetscCall(MatAYPX(*H, gkbnu, A, DIFFERENT_NONZERO_PATTERN)); /* H = A00 + nu*A01*A01' */ 569 570 PetscCall(MatDestroy(&BT)); 571 PetscCall(MatDestroy(&T)); 572 PetscFunctionReturn(PETSC_SUCCESS); 573 } 574 575 PETSC_EXTERN PetscErrorCode PetscOptionsFindPairPrefix_Private(PetscOptions, const char pre[], const char name[], const char *option[], const char *value[], PetscBool *flg); 576 577 static PetscErrorCode PCSetUp_FieldSplit(PC pc) 578 { 579 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 580 PC_FieldSplitLink ilink; 581 PetscInt i, nsplit; 582 PetscBool matnest = PETSC_FALSE; 583 584 PetscFunctionBegin; 585 pc->failedreason = PC_NOERROR; 586 PetscCall(PCFieldSplitSetDefaults(pc)); 587 nsplit = jac->nsplits; 588 ilink = jac->head; 589 if (pc->pmat) PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATNEST, &matnest)); 590 591 /* get the matrices for each split */ 592 if (!jac->issetup) { 593 PetscInt rstart, rend, nslots, bs; 594 595 jac->issetup = PETSC_TRUE; 596 597 /* This is done here instead of in PCFieldSplitSetFields() because may not have matrix at that point */ 598 if (jac->defaultsplit || !ilink->is) { 599 if (jac->bs <= 0) jac->bs = nsplit; 600 } 601 602 /* MatCreateSubMatrix() for [S]BAIJ matrices can only work if the indices include entire blocks of the matrix */ 603 PetscCall(MatGetBlockSize(pc->pmat, &bs)); 604 if (bs > 1 && (jac->bs <= bs || jac->bs % bs)) { 605 PetscBool blk; 606 607 PetscCall(PetscObjectTypeCompareAny((PetscObject)pc->pmat, &blk, MATBAIJ, MATSBAIJ, MATSEQBAIJ, MATSEQSBAIJ, MATMPIBAIJ, MATMPISBAIJ, NULL)); 608 PetscCheck(!blk, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Cannot use MATBAIJ with PCFIELDSPLIT and currently set matrix and PC blocksizes"); 609 } 610 611 if (!matnest) { /* use the matrix blocksize and stride IS to determine the index sets that define the submatrices */ 612 bs = jac->bs; 613 PetscCall(MatGetOwnershipRange(pc->pmat, &rstart, &rend)); 614 nslots = (rend - rstart) / bs; 615 for (i = 0; i < nsplit; i++) { 616 if (jac->defaultsplit) { 617 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)pc), nslots, rstart + i, nsplit, &ilink->is)); 618 PetscCall(PetscObjectReference((PetscObject)ilink->is)); 619 ilink->is_col = ilink->is; 620 } else if (!ilink->is) { 621 PetscBool same_fields = PETSC_TRUE; 622 623 for (PetscInt k = 0; k < ilink->nfields; k++) { 624 if (ilink->fields[k] != ilink->fields_col[k]) same_fields = PETSC_FALSE; 625 } 626 627 if (ilink->nfields > 1) { 628 PetscInt *ii, *jj, j, k, nfields = ilink->nfields, *fields = ilink->fields, *fields_col = ilink->fields_col; 629 630 PetscCall(PetscMalloc1(ilink->nfields * nslots, &ii)); 631 if (!same_fields) PetscCall(PetscMalloc1(ilink->nfields * nslots, &jj)); 632 for (j = 0; j < nslots; j++) { 633 for (k = 0; k < nfields; k++) { 634 ii[nfields * j + k] = rstart + bs * j + fields[k]; 635 if (!same_fields) jj[nfields * j + k] = rstart + bs * j + fields_col[k]; 636 } 637 } 638 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pc), nslots * nfields, ii, PETSC_OWN_POINTER, &ilink->is)); 639 if (!same_fields) PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pc), nslots * nfields, jj, PETSC_OWN_POINTER, &ilink->is_col)); 640 else { 641 PetscCall(PetscObjectReference((PetscObject)ilink->is)); 642 ilink->is_col = ilink->is; 643 } 644 PetscCall(ISSetBlockSize(ilink->is, nfields)); 645 PetscCall(ISSetBlockSize(ilink->is_col, nfields)); 646 } else { 647 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)pc), nslots, rstart + ilink->fields[0], bs, &ilink->is)); 648 if (!same_fields) PetscCall(ISCreateStride(PetscObjectComm((PetscObject)pc), nslots, rstart + ilink->fields_col[0], bs, &ilink->is_col)); 649 else { 650 PetscCall(PetscObjectReference((PetscObject)ilink->is)); 651 ilink->is_col = ilink->is; 652 } 653 } 654 } 655 ilink = ilink->next; 656 } 657 } else { /* use the IS that define the MATNEST to determine the index sets that define the submatrices */ 658 IS *rowis, *colis, *ises = NULL; 659 PetscInt mis, nis; 660 661 PetscCall(MatNestGetSize(pc->pmat, &mis, &nis)); 662 PetscCall(PetscMalloc2(mis, &rowis, nis, &colis)); 663 PetscCall(MatNestGetISs(pc->pmat, rowis, colis)); 664 if (!jac->defaultsplit) PetscCall(PetscMalloc1(mis, &ises)); 665 666 for (i = 0; i < nsplit; i++) { 667 if (jac->defaultsplit) { 668 PetscCall(ISDuplicate(rowis[i], &ilink->is)); 669 PetscCall(PetscObjectReference((PetscObject)ilink->is)); 670 ilink->is_col = ilink->is; 671 } else if (!ilink->is) { 672 if (ilink->nfields > 1) { 673 for (PetscInt j = 0; j < ilink->nfields; j++) ises[j] = rowis[ilink->fields[j]]; 674 PetscCall(ISConcatenate(PetscObjectComm((PetscObject)pc), ilink->nfields, ises, &ilink->is)); 675 } else { 676 PetscCall(ISDuplicate(rowis[ilink->fields[0]], &ilink->is)); 677 } 678 PetscCall(PetscObjectReference((PetscObject)ilink->is)); 679 ilink->is_col = ilink->is; 680 } 681 ilink = ilink->next; 682 } 683 PetscCall(PetscFree2(rowis, colis)); 684 PetscCall(PetscFree(ises)); 685 } 686 } 687 688 ilink = jac->head; 689 if (!jac->pmat) { 690 Vec xtmp; 691 692 PetscCall(MatCreateVecs(pc->pmat, &xtmp, NULL)); 693 PetscCall(PetscMalloc1(nsplit, &jac->pmat)); 694 PetscCall(PetscMalloc2(nsplit, &jac->x, nsplit, &jac->y)); 695 for (i = 0; i < nsplit; i++) { 696 MatNullSpace sp; 697 698 /* Check for matrix attached to IS */ 699 PetscCall(PetscObjectQuery((PetscObject)ilink->is, "pmat", (PetscObject *)&jac->pmat[i])); 700 if (jac->pmat[i]) { 701 PetscCall(PetscObjectReference((PetscObject)jac->pmat[i])); 702 if (jac->type == PC_COMPOSITE_SCHUR) { 703 jac->schur_user = jac->pmat[i]; 704 705 PetscCall(PetscObjectReference((PetscObject)jac->schur_user)); 706 } 707 } else { 708 const char *prefix; 709 PetscCall(MatCreateSubMatrix(pc->pmat, ilink->is, ilink->is_col, MAT_INITIAL_MATRIX, &jac->pmat[i])); 710 PetscCall(MatGetOptionsPrefix(jac->pmat[i], &prefix)); 711 if (!prefix) { 712 PetscCall(KSPGetOptionsPrefix(ilink->ksp, &prefix)); 713 PetscCall(MatSetOptionsPrefix(jac->pmat[i], prefix)); 714 } 715 PetscCall(MatSetFromOptions(jac->pmat[i])); 716 PetscCall(MatViewFromOptions(jac->pmat[i], NULL, "-mat_view")); 717 } 718 /* create work vectors for each split */ 719 PetscCall(MatCreateVecs(jac->pmat[i], &jac->x[i], &jac->y[i])); 720 ilink->x = jac->x[i]; 721 ilink->y = jac->y[i]; 722 ilink->z = NULL; 723 /* compute scatter contexts needed by multiplicative versions and non-default splits */ 724 PetscCall(VecScatterCreate(xtmp, ilink->is, jac->x[i], NULL, &ilink->sctx)); 725 PetscCall(PetscObjectQuery((PetscObject)ilink->is, "nearnullspace", (PetscObject *)&sp)); 726 if (sp) PetscCall(MatSetNearNullSpace(jac->pmat[i], sp)); 727 ilink = ilink->next; 728 } 729 PetscCall(VecDestroy(&xtmp)); 730 } else { 731 MatReuse scall; 732 MatNullSpace *nullsp = NULL; 733 734 if (pc->flag == DIFFERENT_NONZERO_PATTERN) { 735 PetscCall(MatGetNullSpaces(nsplit, jac->pmat, &nullsp)); 736 for (i = 0; i < nsplit; i++) PetscCall(MatDestroy(&jac->pmat[i])); 737 scall = MAT_INITIAL_MATRIX; 738 } else scall = MAT_REUSE_MATRIX; 739 740 for (i = 0; i < nsplit; i++) { 741 Mat pmat; 742 743 /* Check for matrix attached to IS */ 744 PetscCall(PetscObjectQuery((PetscObject)ilink->is, "pmat", (PetscObject *)&pmat)); 745 if (!pmat) PetscCall(MatCreateSubMatrix(pc->pmat, ilink->is, ilink->is_col, scall, &jac->pmat[i])); 746 ilink = ilink->next; 747 } 748 if (nullsp) PetscCall(MatRestoreNullSpaces(nsplit, jac->pmat, &nullsp)); 749 } 750 if (jac->diag_use_amat) { 751 ilink = jac->head; 752 if (!jac->mat) { 753 PetscCall(PetscMalloc1(nsplit, &jac->mat)); 754 for (i = 0; i < nsplit; i++) { 755 PetscCall(MatCreateSubMatrix(pc->mat, ilink->is, ilink->is_col, MAT_INITIAL_MATRIX, &jac->mat[i])); 756 ilink = ilink->next; 757 } 758 } else { 759 MatReuse scall; 760 MatNullSpace *nullsp = NULL; 761 762 if (pc->flag == DIFFERENT_NONZERO_PATTERN) { 763 PetscCall(MatGetNullSpaces(nsplit, jac->mat, &nullsp)); 764 for (i = 0; i < nsplit; i++) PetscCall(MatDestroy(&jac->mat[i])); 765 scall = MAT_INITIAL_MATRIX; 766 } else scall = MAT_REUSE_MATRIX; 767 768 for (i = 0; i < nsplit; i++) { 769 PetscCall(MatCreateSubMatrix(pc->mat, ilink->is, ilink->is_col, scall, &jac->mat[i])); 770 ilink = ilink->next; 771 } 772 if (nullsp) PetscCall(MatRestoreNullSpaces(nsplit, jac->mat, &nullsp)); 773 } 774 } else { 775 jac->mat = jac->pmat; 776 } 777 778 /* Check for null space attached to IS */ 779 ilink = jac->head; 780 for (i = 0; i < nsplit; i++) { 781 MatNullSpace sp; 782 783 PetscCall(PetscObjectQuery((PetscObject)ilink->is, "nullspace", (PetscObject *)&sp)); 784 if (sp) PetscCall(MatSetNullSpace(jac->mat[i], sp)); 785 ilink = ilink->next; 786 } 787 788 if (jac->type != PC_COMPOSITE_ADDITIVE && jac->type != PC_COMPOSITE_SCHUR && jac->type != PC_COMPOSITE_GKB) { 789 /* extract the rows of the matrix associated with each field: used for efficient computation of residual inside algorithm */ 790 /* FIXME: Can/should we reuse jac->mat whenever (jac->diag_use_amat) is true? */ 791 ilink = jac->head; 792 if (nsplit == 2 && jac->type == PC_COMPOSITE_MULTIPLICATIVE) { 793 /* special case need where Afield[0] is not needed and only certain columns of Afield[1] are needed since update is only on those rows of the solution */ 794 if (!jac->Afield) { 795 PetscCall(PetscCalloc1(nsplit, &jac->Afield)); 796 if (jac->offdiag_use_amat) { 797 PetscCall(MatCreateSubMatrix(pc->mat, ilink->next->is, ilink->is, MAT_INITIAL_MATRIX, &jac->Afield[1])); 798 } else { 799 PetscCall(MatCreateSubMatrix(pc->pmat, ilink->next->is, ilink->is, MAT_INITIAL_MATRIX, &jac->Afield[1])); 800 } 801 } else { 802 MatReuse scall; 803 804 if (pc->flag == DIFFERENT_NONZERO_PATTERN) { 805 PetscCall(MatDestroy(&jac->Afield[1])); 806 scall = MAT_INITIAL_MATRIX; 807 } else scall = MAT_REUSE_MATRIX; 808 809 if (jac->offdiag_use_amat) { 810 PetscCall(MatCreateSubMatrix(pc->mat, ilink->next->is, ilink->is, scall, &jac->Afield[1])); 811 } else { 812 PetscCall(MatCreateSubMatrix(pc->pmat, ilink->next->is, ilink->is, scall, &jac->Afield[1])); 813 } 814 } 815 } else { 816 if (!jac->Afield) { 817 PetscCall(PetscMalloc1(nsplit, &jac->Afield)); 818 for (i = 0; i < nsplit; i++) { 819 if (jac->offdiag_use_amat) { 820 PetscCall(MatCreateSubMatrix(pc->mat, ilink->is, NULL, MAT_INITIAL_MATRIX, &jac->Afield[i])); 821 } else { 822 PetscCall(MatCreateSubMatrix(pc->pmat, ilink->is, NULL, MAT_INITIAL_MATRIX, &jac->Afield[i])); 823 } 824 ilink = ilink->next; 825 } 826 } else { 827 MatReuse scall; 828 if (pc->flag == DIFFERENT_NONZERO_PATTERN) { 829 for (i = 0; i < nsplit; i++) PetscCall(MatDestroy(&jac->Afield[i])); 830 scall = MAT_INITIAL_MATRIX; 831 } else scall = MAT_REUSE_MATRIX; 832 833 for (i = 0; i < nsplit; i++) { 834 if (jac->offdiag_use_amat) { 835 PetscCall(MatCreateSubMatrix(pc->mat, ilink->is, NULL, scall, &jac->Afield[i])); 836 } else { 837 PetscCall(MatCreateSubMatrix(pc->pmat, ilink->is, NULL, scall, &jac->Afield[i])); 838 } 839 ilink = ilink->next; 840 } 841 } 842 } 843 } 844 845 if (jac->type == PC_COMPOSITE_SCHUR) { 846 PetscBool isset, isspd = PETSC_FALSE, issym = PETSC_FALSE, flg; 847 char lscname[256]; 848 PetscObject LSC_L; 849 850 PetscCheck(nsplit == 2, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "To use Schur complement preconditioner you must have exactly 2 fields"); 851 852 /* If pc->mat is SPD, don't scale by -1 the Schur complement */ 853 PetscCall(MatIsSPDKnown(pc->pmat, &isset, &isspd)); 854 if (jac->schurscale == (PetscScalar)-1.0) jac->schurscale = (isset && isspd) ? 1.0 : -1.0; 855 PetscCall(MatIsSymmetricKnown(pc->pmat, &isset, &issym)); 856 857 PetscCall(PetscObjectTypeCompareAny(jac->offdiag_use_amat ? (PetscObject)pc->mat : (PetscObject)pc->pmat, &flg, MATSEQSBAIJ, MATMPISBAIJ, "")); 858 859 if (jac->schur) { 860 KSP kspA = jac->head->ksp, kspInner = NULL, kspUpper = jac->kspupper; 861 MatReuse scall; 862 863 if (pc->flag == DIFFERENT_NONZERO_PATTERN) { 864 scall = MAT_INITIAL_MATRIX; 865 PetscCall(MatDestroy(&jac->B)); 866 PetscCall(MatDestroy(&jac->C)); 867 } else scall = MAT_REUSE_MATRIX; 868 869 PetscCall(MatSchurComplementGetKSP(jac->schur, &kspInner)); 870 ilink = jac->head; 871 PetscCall(MatCreateSubMatrix(jac->offdiag_use_amat ? pc->mat : pc->pmat, ilink->is, ilink->next->is, scall, &jac->B)); 872 if (!flg) PetscCall(MatCreateSubMatrix(jac->offdiag_use_amat ? pc->mat : pc->pmat, ilink->next->is, ilink->is, scall, &jac->C)); 873 else { 874 PetscCall(MatIsHermitianKnown(jac->offdiag_use_amat ? pc->mat : pc->pmat, &isset, &flg)); 875 if (isset && flg) PetscCall(MatCreateHermitianTranspose(jac->B, &jac->C)); 876 else PetscCall(MatCreateTranspose(jac->B, &jac->C)); 877 } 878 ilink = ilink->next; 879 PetscCall(MatSchurComplementUpdateSubMatrices(jac->schur, jac->mat[0], jac->pmat[0], jac->B, jac->C, jac->mat[1])); 880 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) { 881 PetscCall(MatDestroy(&jac->schurp)); 882 PetscCall(MatSchurComplementGetPmat(jac->schur, MAT_INITIAL_MATRIX, &jac->schurp)); 883 } else if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_FULL && jac->kspupper != jac->head->ksp) { 884 PetscCall(MatDestroy(&jac->schur_user)); 885 PetscCall(MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user)); 886 } 887 if (kspA != kspInner) PetscCall(KSPSetOperators(kspA, jac->mat[0], jac->pmat[0])); 888 if (kspUpper != kspA) PetscCall(KSPSetOperators(kspUpper, jac->mat[0], jac->pmat[0])); 889 PetscCall(KSPSetOperators(jac->kspschur, jac->schur, FieldSplitSchurPre(jac))); 890 } else { 891 const char *Dprefix; 892 char schurprefix[256], schurmatprefix[256]; 893 char schurtestoption[256]; 894 MatNullSpace sp; 895 KSP kspt; 896 897 /* extract the A01 and A10 matrices */ 898 ilink = jac->head; 899 PetscCall(MatCreateSubMatrix(jac->offdiag_use_amat ? pc->mat : pc->pmat, ilink->is, ilink->next->is, MAT_INITIAL_MATRIX, &jac->B)); 900 if (!flg) PetscCall(MatCreateSubMatrix(jac->offdiag_use_amat ? pc->mat : pc->pmat, ilink->next->is, ilink->is, MAT_INITIAL_MATRIX, &jac->C)); 901 else { 902 PetscCall(MatIsHermitianKnown(jac->offdiag_use_amat ? pc->mat : pc->pmat, &isset, &flg)); 903 if (isset && flg) PetscCall(MatCreateHermitianTranspose(jac->B, &jac->C)); 904 else PetscCall(MatCreateTranspose(jac->B, &jac->C)); 905 } 906 ilink = ilink->next; 907 /* Use mat[0] (diagonal block of Amat) preconditioned by pmat[0] to define Schur complement */ 908 PetscCall(MatCreate(((PetscObject)jac->mat[0])->comm, &jac->schur)); 909 PetscCall(MatSetType(jac->schur, MATSCHURCOMPLEMENT)); 910 PetscCall(MatSchurComplementSetSubMatrices(jac->schur, jac->mat[0], jac->pmat[0], jac->B, jac->C, jac->mat[1])); 911 PetscCall(PetscSNPrintf(schurmatprefix, sizeof(schurmatprefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname)); 912 PetscCall(MatSetOptionsPrefix(jac->schur, schurmatprefix)); 913 PetscCall(MatSchurComplementGetKSP(jac->schur, &kspt)); 914 PetscCall(KSPSetOptionsPrefix(kspt, schurmatprefix)); 915 916 /* Note: this is not true in general */ 917 PetscCall(MatGetNullSpace(jac->mat[1], &sp)); 918 if (sp) PetscCall(MatSetNullSpace(jac->schur, sp)); 919 920 PetscCall(PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_inner_", ilink->splitname)); 921 PetscCall(PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, schurtestoption, NULL, NULL, &flg)); 922 if (flg) { 923 DM dmInner; 924 KSP kspInner; 925 PC pcInner; 926 927 PetscCall(MatSchurComplementGetKSP(jac->schur, &kspInner)); 928 PetscCall(KSPReset(kspInner)); 929 PetscCall(KSPSetOperators(kspInner, jac->mat[0], jac->pmat[0])); 930 PetscCall(PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_inner_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname)); 931 /* Indent this deeper to emphasize the "inner" nature of this solver. */ 932 PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspInner, (PetscObject)pc, 2)); 933 PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspInner->pc, (PetscObject)pc, 2)); 934 PetscCall(KSPSetOptionsPrefix(kspInner, schurprefix)); 935 936 /* Set DM for new solver */ 937 PetscCall(KSPGetDM(jac->head->ksp, &dmInner)); 938 PetscCall(KSPSetDM(kspInner, dmInner)); 939 PetscCall(KSPSetDMActive(kspInner, KSP_DMACTIVE_ALL, PETSC_FALSE)); 940 941 /* Defaults to PCKSP as preconditioner */ 942 PetscCall(KSPGetPC(kspInner, &pcInner)); 943 PetscCall(PCSetType(pcInner, PCKSP)); 944 PetscCall(PCKSPSetKSP(pcInner, jac->head->ksp)); 945 } else { 946 /* Use the outer solver for the inner solve, but revert the KSPPREONLY from PCFieldSplitSetFields_FieldSplit or 947 * PCFieldSplitSetIS_FieldSplit. We don't want KSPPREONLY because it makes the Schur complement inexact, 948 * preventing Schur complement reduction to be an accurate solve. Usually when an iterative solver is used for 949 * S = D - C A_inner^{-1} B, we expect S to be defined using an accurate definition of A_inner^{-1}, so we make 950 * GMRES the default. Note that it is also common to use PREONLY for S, in which case S may not be used 951 * directly, and the user is responsible for setting an inexact method for fieldsplit's A^{-1}. */ 952 PetscCall(KSPSetType(jac->head->ksp, KSPGMRES)); 953 PetscCall(MatSchurComplementSetKSP(jac->schur, jac->head->ksp)); 954 } 955 PetscCall(KSPSetOperators(jac->head->ksp, jac->mat[0], jac->pmat[0])); 956 PetscCall(KSPSetFromOptions(jac->head->ksp)); 957 PetscCall(MatSetFromOptions(jac->schur)); 958 959 PetscCall(PetscObjectTypeCompare((PetscObject)jac->schur, MATSCHURCOMPLEMENT, &flg)); 960 if (flg) { /* Need to do this otherwise PCSetUp_KSP will overwrite the amat of jac->head->ksp */ 961 KSP kspInner; 962 PC pcInner; 963 964 PetscCall(MatSchurComplementGetKSP(jac->schur, &kspInner)); 965 PetscCall(KSPGetPC(kspInner, &pcInner)); 966 PetscCall(PetscObjectTypeCompare((PetscObject)pcInner, PCKSP, &flg)); 967 if (flg) { 968 KSP ksp; 969 970 PetscCall(PCKSPGetKSP(pcInner, &ksp)); 971 if (ksp == jac->head->ksp) PetscCall(PCSetUseAmat(pcInner, PETSC_TRUE)); 972 } 973 } 974 PetscCall(PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_upper_", ilink->splitname)); 975 PetscCall(PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, schurtestoption, NULL, NULL, &flg)); 976 if (flg) { 977 DM dmInner; 978 979 PetscCall(PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_upper_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname)); 980 PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspupper)); 981 PetscCall(KSPSetNestLevel(jac->kspupper, pc->kspnestlevel)); 982 PetscCall(KSPSetErrorIfNotConverged(jac->kspupper, pc->erroriffailure)); 983 PetscCall(KSPSetOptionsPrefix(jac->kspupper, schurprefix)); 984 PetscCall(PetscObjectIncrementTabLevel((PetscObject)jac->kspupper, (PetscObject)pc, 1)); 985 PetscCall(PetscObjectIncrementTabLevel((PetscObject)jac->kspupper->pc, (PetscObject)pc, 1)); 986 PetscCall(KSPGetDM(jac->head->ksp, &dmInner)); 987 PetscCall(KSPSetDM(jac->kspupper, dmInner)); 988 PetscCall(KSPSetDMActive(jac->kspupper, KSP_DMACTIVE_ALL, PETSC_FALSE)); 989 PetscCall(KSPSetFromOptions(jac->kspupper)); 990 PetscCall(KSPSetOperators(jac->kspupper, jac->mat[0], jac->pmat[0])); 991 PetscCall(VecDuplicate(jac->head->x, &jac->head->z)); 992 } else { 993 jac->kspupper = jac->head->ksp; 994 PetscCall(PetscObjectReference((PetscObject)jac->head->ksp)); 995 } 996 997 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) PetscCall(MatSchurComplementGetPmat(jac->schur, MAT_INITIAL_MATRIX, &jac->schurp)); 998 PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspschur)); 999 PetscCall(KSPSetNestLevel(jac->kspschur, pc->kspnestlevel)); 1000 PetscCall(KSPSetErrorIfNotConverged(jac->kspschur, pc->erroriffailure)); 1001 PetscCall(PetscObjectIncrementTabLevel((PetscObject)jac->kspschur, (PetscObject)pc, 1)); 1002 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) { 1003 PC pcschur; 1004 PetscCall(KSPGetPC(jac->kspschur, &pcschur)); 1005 PetscCall(PCSetType(pcschur, PCNONE)); 1006 /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */ 1007 } else if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_FULL) { 1008 if (jac->schurfactorization != PC_FIELDSPLIT_SCHUR_FACT_FULL || jac->kspupper != jac->head->ksp) PetscCall(MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user)); 1009 } 1010 PetscCall(KSPSetOperators(jac->kspschur, jac->schur, FieldSplitSchurPre(jac))); 1011 PetscCall(KSPGetOptionsPrefix(jac->head->next->ksp, &Dprefix)); 1012 PetscCall(KSPSetOptionsPrefix(jac->kspschur, Dprefix)); 1013 /* propagate DM */ 1014 { 1015 DM sdm; 1016 PetscCall(KSPGetDM(jac->head->next->ksp, &sdm)); 1017 if (sdm) { 1018 PetscCall(KSPSetDM(jac->kspschur, sdm)); 1019 PetscCall(KSPSetDMActive(jac->kspschur, KSP_DMACTIVE_ALL, PETSC_FALSE)); 1020 } 1021 } 1022 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 1023 /* need to call this every time, since the jac->kspschur is freshly created, otherwise its options never get set */ 1024 PetscCall(KSPSetFromOptions(jac->kspschur)); 1025 } 1026 PetscCall(MatAssemblyBegin(jac->schur, MAT_FINAL_ASSEMBLY)); 1027 PetscCall(MatAssemblyEnd(jac->schur, MAT_FINAL_ASSEMBLY)); 1028 if (issym) PetscCall(MatSetOption(jac->schur, MAT_SYMMETRIC, PETSC_TRUE)); 1029 if (isspd) PetscCall(MatSetOption(jac->schur, MAT_SPD, PETSC_TRUE)); 1030 1031 /* HACK: special support to forward L and Lp matrices that might be used by PCLSC */ 1032 PetscCall(PetscSNPrintf(lscname, sizeof(lscname), "%s_LSC_L", ilink->splitname)); 1033 PetscCall(PetscObjectQuery((PetscObject)pc->mat, lscname, &LSC_L)); 1034 if (!LSC_L) PetscCall(PetscObjectQuery((PetscObject)pc->pmat, lscname, &LSC_L)); 1035 if (LSC_L) PetscCall(PetscObjectCompose((PetscObject)jac->schur, "LSC_L", LSC_L)); 1036 PetscCall(PetscSNPrintf(lscname, sizeof(lscname), "%s_LSC_Lp", ilink->splitname)); 1037 PetscCall(PetscObjectQuery((PetscObject)pc->pmat, lscname, &LSC_L)); 1038 if (!LSC_L) PetscCall(PetscObjectQuery((PetscObject)pc->mat, lscname, &LSC_L)); 1039 if (LSC_L) PetscCall(PetscObjectCompose((PetscObject)jac->schur, "LSC_Lp", LSC_L)); 1040 } else if (jac->type == PC_COMPOSITE_GKB) { 1041 PetscCheck(nsplit == 2, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "To use GKB preconditioner you must have exactly 2 fields"); 1042 ilink = jac->head; 1043 PetscCall(MatCreateSubMatrix(jac->offdiag_use_amat ? pc->mat : pc->pmat, ilink->is, ilink->next->is, MAT_INITIAL_MATRIX, &jac->B)); 1044 /* Create work vectors for GKB algorithm */ 1045 PetscCall(VecDuplicate(ilink->x, &jac->u)); 1046 PetscCall(VecDuplicate(ilink->x, &jac->Hu)); 1047 PetscCall(VecDuplicate(ilink->x, &jac->w2)); 1048 PetscCall(MatCreateSubMatrix(jac->offdiag_use_amat ? pc->mat : pc->pmat, ilink->next->is, ilink->is, MAT_INITIAL_MATRIX, &jac->C)); 1049 ilink = ilink->next; 1050 /* Create work vectors for GKB algorithm */ 1051 PetscCall(VecDuplicate(ilink->x, &jac->v)); 1052 PetscCall(VecDuplicate(ilink->x, &jac->d)); 1053 PetscCall(VecDuplicate(ilink->x, &jac->w1)); 1054 PetscCall(MatGolubKahanComputeExplicitOperator(jac->mat[0], jac->B, jac->C, &jac->H, jac->gkbnu)); 1055 PetscCall(PetscCalloc1(jac->gkbdelay, &jac->vecz)); 1056 1057 ilink = jac->head; 1058 PetscCall(KSPSetOperators(ilink->ksp, jac->H, jac->H)); 1059 if (!jac->suboptionsset) PetscCall(KSPSetFromOptions(ilink->ksp)); 1060 /* Create gkb_monitor context */ 1061 if (jac->gkbmonitor) { 1062 PetscInt tablevel; 1063 PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &jac->gkbviewer)); 1064 PetscCall(PetscViewerSetType(jac->gkbviewer, PETSCVIEWERASCII)); 1065 PetscCall(PetscObjectGetTabLevel((PetscObject)ilink->ksp, &tablevel)); 1066 PetscCall(PetscViewerASCIISetTab(jac->gkbviewer, tablevel)); 1067 PetscCall(PetscObjectIncrementTabLevel((PetscObject)ilink->ksp, (PetscObject)ilink->ksp, 1)); 1068 } 1069 } else { 1070 /* set up the individual splits' PCs */ 1071 i = 0; 1072 ilink = jac->head; 1073 while (ilink) { 1074 PetscCall(KSPSetOperators(ilink->ksp, jac->mat[i], jac->pmat[i])); 1075 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 1076 if (!jac->suboptionsset) PetscCall(KSPSetFromOptions(ilink->ksp)); 1077 i++; 1078 ilink = ilink->next; 1079 } 1080 } 1081 1082 /* Set coordinates to the sub PC objects whenever these are set */ 1083 if (jac->coordinates_set) { 1084 PC pc_coords; 1085 if (jac->type == PC_COMPOSITE_SCHUR) { 1086 // Head is first block. 1087 PetscCall(KSPGetPC(jac->head->ksp, &pc_coords)); 1088 PetscCall(PCSetCoordinates(pc_coords, jac->head->dim, jac->head->ndofs, jac->head->coords)); 1089 // Second one is Schur block, but its KSP object is in kspschur. 1090 PetscCall(KSPGetPC(jac->kspschur, &pc_coords)); 1091 PetscCall(PCSetCoordinates(pc_coords, jac->head->next->dim, jac->head->next->ndofs, jac->head->next->coords)); 1092 } else if (jac->type == PC_COMPOSITE_GKB) { 1093 PetscCall(PetscInfo(pc, "Warning: Setting coordinates does nothing for the GKB Fieldpslit preconditioner\n")); 1094 } else { 1095 ilink = jac->head; 1096 while (ilink) { 1097 PetscCall(KSPGetPC(ilink->ksp, &pc_coords)); 1098 PetscCall(PCSetCoordinates(pc_coords, ilink->dim, ilink->ndofs, ilink->coords)); 1099 ilink = ilink->next; 1100 } 1101 } 1102 } 1103 1104 jac->suboptionsset = PETSC_TRUE; 1105 PetscFunctionReturn(PETSC_SUCCESS); 1106 } 1107 1108 static PetscErrorCode PCSetUpOnBlocks_FieldSplit_Schur(PC pc) 1109 { 1110 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 1111 PC_FieldSplitLink ilinkA = jac->head; 1112 KSP kspA = ilinkA->ksp, kspUpper = jac->kspupper; 1113 1114 PetscFunctionBegin; 1115 if (jac->schurfactorization == PC_FIELDSPLIT_SCHUR_FACT_FULL && kspUpper != kspA) { 1116 PetscCall(KSPSetUp(kspUpper)); 1117 PetscCall(KSPSetUpOnBlocks(kspUpper)); 1118 } 1119 PetscCall(KSPSetUp(kspA)); 1120 PetscCall(KSPSetUpOnBlocks(kspA)); 1121 if (jac->schurpre != PC_FIELDSPLIT_SCHUR_PRE_FULL) { 1122 PetscCall(KSPSetUp(jac->kspschur)); 1123 PetscCall(KSPSetUpOnBlocks(jac->kspschur)); 1124 } else if (kspUpper == kspA && jac->schurfactorization == PC_FIELDSPLIT_SCHUR_FACT_FULL) { 1125 Mat A; 1126 PetscInt m, M, N; 1127 VecType vtype; 1128 PetscMemType mtype; 1129 PetscScalar *array; 1130 1131 PetscCall(MatGetSize(jac->B, &M, &N)); 1132 PetscCall(MatGetLocalSize(jac->B, &m, NULL)); 1133 PetscCall(MatGetVecType(jac->B, &vtype)); 1134 PetscCall(VecGetArrayAndMemType(ilinkA->x, &array, &mtype)); 1135 PetscCall(VecRestoreArrayAndMemType(ilinkA->x, &array)); 1136 PetscCall(PetscObjectQuery((PetscObject)jac->schur, "AinvB", (PetscObject *)&A)); 1137 if (A) { 1138 PetscInt P; 1139 1140 PetscCall(MatGetSize(A, NULL, &P)); 1141 if (P < N + 1) { // need to recreate AinvB, otherwise, the Schur complement won't be updated 1142 PetscCall(PetscObjectCompose((PetscObject)jac->schur, "AinvB", NULL)); 1143 A = NULL; 1144 } 1145 } 1146 if (!A) { 1147 if (PetscMemTypeHost(mtype) || (!PetscDefined(HAVE_CUDA) && !PetscDefined(HAVE_HIP))) PetscCall(PetscMalloc1(m * (N + 1), &array)); 1148 #if PetscDefined(HAVE_CUDA) 1149 else if (PetscMemTypeCUDA(mtype)) PetscCallCUDA(cudaMalloc((void **)&array, sizeof(PetscScalar) * m * (N + 1))); 1150 #endif 1151 #if PetscDefined(HAVE_HIP) 1152 else if (PetscMemTypeHIP(mtype)) PetscCallHIP(hipMalloc((void **)&array, sizeof(PetscScalar) * m * (N + 1))); 1153 #endif 1154 PetscCall(MatCreateDenseFromVecType(PetscObjectComm((PetscObject)jac->schur), vtype, m, PETSC_DECIDE, M, N + 1, -1, array, &A)); // number of columns of the Schur complement plus one 1155 PetscCall(PetscObjectCompose((PetscObject)jac->schur, "AinvB", (PetscObject)A)); 1156 PetscCall(MatDestroy(&A)); 1157 } 1158 } 1159 PetscFunctionReturn(PETSC_SUCCESS); 1160 } 1161 1162 static PetscErrorCode PCSetUpOnBlocks_FieldSplit(PC pc) 1163 { 1164 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 1165 PC_FieldSplitLink ilink = jac->head; 1166 1167 PetscFunctionBegin; 1168 while (ilink) { 1169 PetscCall(KSPSetUp(ilink->ksp)); 1170 PetscCall(KSPSetUpOnBlocks(ilink->ksp)); 1171 ilink = ilink->next; 1172 } 1173 PetscFunctionReturn(PETSC_SUCCESS); 1174 } 1175 1176 static PetscErrorCode PCSetUpOnBlocks_FieldSplit_GKB(PC pc) 1177 { 1178 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 1179 PC_FieldSplitLink ilinkA = jac->head; 1180 KSP ksp = ilinkA->ksp; 1181 1182 PetscFunctionBegin; 1183 PetscCall(KSPSetUp(ksp)); 1184 PetscCall(KSPSetUpOnBlocks(ksp)); 1185 PetscFunctionReturn(PETSC_SUCCESS); 1186 } 1187 1188 static PetscErrorCode PCApply_FieldSplit_Schur(PC pc, Vec x, Vec y) 1189 { 1190 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 1191 PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next; 1192 KSP kspA = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper; 1193 Mat AinvB = NULL; 1194 PetscInt N, P; 1195 1196 PetscFunctionBegin; 1197 switch (jac->schurfactorization) { 1198 case PC_FIELDSPLIT_SCHUR_FACT_DIAG: 1199 /* [A00 0; 0 -S], positive definite, suitable for MINRES */ 1200 PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD)); 1201 PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD)); 1202 PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD)); 1203 PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL)); 1204 PetscCall(KSPSolve(kspA, ilinkA->x, ilinkA->y)); 1205 PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y)); 1206 PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL)); 1207 PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1208 PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD)); 1209 PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL)); 1210 PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1)); 1211 PetscCall(KSPSolve(jac->kspschur, ilinkD->x, ilinkD->y)); 1212 PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y)); 1213 PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1)); 1214 PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL)); 1215 PetscCall(VecScale(ilinkD->y, jac->schurscale)); 1216 PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1217 PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1218 PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1219 break; 1220 case PC_FIELDSPLIT_SCHUR_FACT_LOWER: 1221 /* [A00 0; A10 S], suitable for left preconditioning */ 1222 PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD)); 1223 PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD)); 1224 PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL)); 1225 PetscCall(KSPSolve(kspA, ilinkA->x, ilinkA->y)); 1226 PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y)); 1227 PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL)); 1228 PetscCall(MatMult(jac->C, ilinkA->y, ilinkD->x)); 1229 PetscCall(VecScale(ilinkD->x, -1.)); 1230 PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD)); 1231 PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1232 PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD)); 1233 PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL)); 1234 PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1)); 1235 PetscCall(KSPSolve(jac->kspschur, ilinkD->x, ilinkD->y)); 1236 PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1)); 1237 PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y)); 1238 PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL)); 1239 PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1240 PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1241 PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1242 break; 1243 case PC_FIELDSPLIT_SCHUR_FACT_UPPER: 1244 /* [A00 A01; 0 S], suitable for right preconditioning */ 1245 PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD)); 1246 PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD)); 1247 PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL)); 1248 PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1)); 1249 PetscCall(KSPSolve(jac->kspschur, ilinkD->x, ilinkD->y)); 1250 PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1)); 1251 PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y)); 1252 PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL)); 1253 PetscCall(MatMult(jac->B, ilinkD->y, ilinkA->x)); 1254 PetscCall(VecScale(ilinkA->x, -1.)); 1255 PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, ADD_VALUES, SCATTER_FORWARD)); 1256 PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1257 PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, ADD_VALUES, SCATTER_FORWARD)); 1258 PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL)); 1259 PetscCall(KSPSolve(kspA, ilinkA->x, ilinkA->y)); 1260 PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y)); 1261 PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL)); 1262 PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1263 PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1264 PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1265 break; 1266 case PC_FIELDSPLIT_SCHUR_FACT_FULL: 1267 /* [1 0; A10 A00^{-1} 1] [A00 0; 0 S] [1 A00^{-1}A01; 0 1] */ 1268 PetscCall(MatGetSize(jac->B, NULL, &P)); 1269 N = P; 1270 PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD)); 1271 PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD)); 1272 PetscCall(PetscLogEventBegin(KSP_Solve_FS_L, kspLower, ilinkA->x, ilinkA->y, NULL)); 1273 if (kspUpper == kspA) { 1274 PetscCall(PetscObjectQuery((PetscObject)jac->schur, "AinvB", (PetscObject *)&AinvB)); 1275 if (AinvB) { 1276 PetscCall(MatGetSize(AinvB, NULL, &N)); 1277 if (N > P) { // first time PCApply_FieldSplit_Schur() is called 1278 PetscMemType mtype; 1279 Vec c = NULL; 1280 PetscScalar *array; 1281 PetscInt m, M; 1282 1283 PetscCall(MatGetSize(jac->B, &M, NULL)); 1284 PetscCall(MatGetLocalSize(jac->B, &m, NULL)); 1285 PetscCall(MatDenseGetArrayAndMemType(AinvB, &array, &mtype)); 1286 if (PetscMemTypeHost(mtype) || (!PetscDefined(HAVE_CUDA) && !PetscDefined(HAVE_HIP))) PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)jac->schur), 1, m, M, array + m * P, &c)); 1287 #if PetscDefined(HAVE_CUDA) 1288 else if (PetscMemTypeCUDA(mtype)) PetscCall(VecCreateMPICUDAWithArray(PetscObjectComm((PetscObject)jac->schur), 1, m, M, array + m * P, &c)); 1289 #endif 1290 #if PetscDefined(HAVE_HIP) 1291 else if (PetscMemTypeHIP(mtype)) PetscCall(VecCreateMPIHIPWithArray(PetscObjectComm((PetscObject)jac->schur), 1, m, M, array + m * P, &c)); 1292 #endif 1293 PetscCall(MatDenseRestoreArrayAndMemType(AinvB, &array)); 1294 PetscCall(VecCopy(ilinkA->x, c)); 1295 PetscCall(MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user)); 1296 PetscCall(KSPSetOperators(jac->kspschur, jac->schur, jac->schur_user)); 1297 PetscCall(VecCopy(c, ilinkA->y)); // retrieve the solution as the last column of the composed Mat 1298 PetscCall(VecDestroy(&c)); 1299 } 1300 } 1301 } 1302 if (N == P) PetscCall(KSPSolve(kspLower, ilinkA->x, ilinkA->y)); 1303 PetscCall(KSPCheckSolve(kspLower, pc, ilinkA->y)); 1304 PetscCall(PetscLogEventEnd(KSP_Solve_FS_L, kspLower, ilinkA->x, ilinkA->y, NULL)); 1305 PetscCall(MatMult(jac->C, ilinkA->y, ilinkD->x)); 1306 PetscCall(VecScale(ilinkD->x, -1.0)); 1307 PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD)); 1308 PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD)); 1309 1310 PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL)); 1311 PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1)); 1312 PetscCall(KSPSolve(jac->kspschur, ilinkD->x, ilinkD->y)); 1313 PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1)); 1314 PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y)); 1315 PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL)); 1316 PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1317 1318 if (kspUpper == kspA) { 1319 if (!AinvB) { 1320 PetscCall(MatMult(jac->B, ilinkD->y, ilinkA->y)); 1321 PetscCall(VecAXPY(ilinkA->x, -1.0, ilinkA->y)); 1322 PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL)); 1323 PetscCall(KSPSolve(kspA, ilinkA->x, ilinkA->y)); 1324 PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y)); 1325 PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL)); 1326 } else PetscCall(MatMultAdd(AinvB, ilinkD->y, ilinkA->y, ilinkA->y)); 1327 } else { 1328 PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL)); 1329 PetscCall(KSPSolve(kspA, ilinkA->x, ilinkA->y)); 1330 PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y)); 1331 PetscCall(MatMult(jac->B, ilinkD->y, ilinkA->x)); 1332 PetscCall(PetscLogEventBegin(KSP_Solve_FS_U, kspUpper, ilinkA->x, ilinkA->z, NULL)); 1333 PetscCall(KSPSolve(kspUpper, ilinkA->x, ilinkA->z)); 1334 PetscCall(KSPCheckSolve(kspUpper, pc, ilinkA->z)); 1335 PetscCall(PetscLogEventEnd(KSP_Solve_FS_U, kspUpper, ilinkA->x, ilinkA->z, NULL)); 1336 PetscCall(VecAXPY(ilinkA->y, -1.0, ilinkA->z)); 1337 } 1338 PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1339 PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1340 PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1341 } 1342 PetscFunctionReturn(PETSC_SUCCESS); 1343 } 1344 1345 /* 1346 PCFieldSplitCreateWorkMats_Private - Allocate per-field dense work matrices for multi-RHS 1347 1348 Input Parameters: 1349 + pc - the PC context 1350 - X - matrix to copy column-layout from 1351 1352 Notes: 1353 If matrices already exist with correct column count, they are reused. 1354 If column count changed, old matrices are destroyed and new ones created. 1355 */ 1356 static PetscErrorCode PCFieldSplitCreateWorkMats_Private(PC pc, Mat X) 1357 { 1358 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 1359 PC_FieldSplitLink ilink = jac->head; 1360 PetscInt mx, Mx, my, My, N; 1361 1362 PetscFunctionBegin; 1363 while (ilink) { 1364 /* check if reallocation needed (previous allocation with wrong column count) */ 1365 if (ilink->X) { 1366 PetscCall(MatGetSize(ilink->X, NULL, &N)); 1367 if (N != X->cmap->N) { 1368 PetscCall(MatDestroy(&ilink->X)); 1369 PetscCall(MatDestroy(&ilink->Y)); 1370 PetscCall(MatDestroy(&ilink->Z)); 1371 } 1372 } 1373 /* create if needed */ 1374 if (!ilink->X) { 1375 VecType xtype, ytype; 1376 1377 PetscCall(VecGetType(ilink->x, &xtype)); 1378 PetscCall(VecGetType(ilink->y, &ytype)); 1379 PetscCall(VecGetLocalSize(ilink->x, &mx)); 1380 PetscCall(VecGetSize(ilink->x, &Mx)); 1381 PetscCall(VecGetLocalSize(ilink->y, &my)); 1382 PetscCall(VecGetSize(ilink->y, &My)); 1383 /* use default lda */ 1384 PetscCall(MatCreateDenseFromVecType(PetscObjectComm((PetscObject)pc), xtype, mx, X->cmap->n, Mx, X->cmap->N, -1, NULL, &ilink->X)); 1385 PetscCall(MatCreateDenseFromVecType(PetscObjectComm((PetscObject)pc), ytype, my, X->cmap->n, My, X->cmap->N, -1, NULL, &ilink->Y)); 1386 } 1387 ilink = ilink->next; 1388 } 1389 PetscFunctionReturn(PETSC_SUCCESS); 1390 } 1391 1392 static PetscErrorCode PCMatApply_FieldSplit_Schur(PC pc, Mat X, Mat Y) 1393 { 1394 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 1395 PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next; 1396 KSP kspA = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper; 1397 Mat AinvB = NULL; 1398 PetscInt N, P; 1399 1400 PetscFunctionBegin; 1401 /* create working matrices with the correct number of columns */ 1402 PetscCall(PCFieldSplitCreateWorkMats_Private(pc, X)); 1403 switch (jac->schurfactorization) { 1404 case PC_FIELDSPLIT_SCHUR_FACT_DIAG: 1405 /* [A00 0; 0 -S], positive definite, suitable for MINRES */ 1406 PetscCall(MatDenseScatter_Private(ilinkA->sctx, X, ilinkA->X, INSERT_VALUES, SCATTER_FORWARD)); 1407 PetscCall(MatDenseScatter_Private(ilinkD->sctx, X, ilinkD->X, INSERT_VALUES, SCATTER_FORWARD)); 1408 PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->X, ilinkA->Y, NULL)); 1409 PetscCall(KSPMatSolve(kspA, ilinkA->X, ilinkA->Y)); 1410 PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->X, ilinkA->Y, NULL)); 1411 PetscCall(MatDenseScatter_Private(ilinkA->sctx, ilinkA->Y, Y, INSERT_VALUES, SCATTER_REVERSE)); 1412 PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->X, ilinkD->Y, NULL)); 1413 PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1)); 1414 PetscCall(KSPMatSolve(jac->kspschur, ilinkD->X, ilinkD->Y)); 1415 PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1)); 1416 PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->X, ilinkD->Y, NULL)); 1417 PetscCall(MatScale(ilinkD->Y, jac->schurscale)); 1418 PetscCall(MatDenseScatter_Private(ilinkD->sctx, ilinkD->Y, Y, INSERT_VALUES, SCATTER_REVERSE)); 1419 break; 1420 case PC_FIELDSPLIT_SCHUR_FACT_LOWER: 1421 /* [A00 0; A10 S], suitable for left preconditioning */ 1422 PetscCall(MatDenseScatter_Private(ilinkA->sctx, X, ilinkA->X, INSERT_VALUES, SCATTER_FORWARD)); 1423 PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->X, ilinkA->Y, NULL)); 1424 PetscCall(KSPMatSolve(kspA, ilinkA->X, ilinkA->Y)); 1425 PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->X, ilinkA->Y, NULL)); 1426 PetscCall(MatMatMult(jac->C, ilinkA->Y, MAT_REUSE_MATRIX, PETSC_DETERMINE, &ilinkD->X)); 1427 PetscCall(MatScale(ilinkD->X, -1.0)); 1428 PetscCall(MatDenseScatter_Private(ilinkD->sctx, X, ilinkD->X, ADD_VALUES, SCATTER_FORWARD)); 1429 PetscCall(MatDenseScatter_Private(ilinkA->sctx, ilinkA->Y, Y, INSERT_VALUES, SCATTER_REVERSE)); 1430 PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->X, ilinkD->Y, NULL)); 1431 PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1)); 1432 PetscCall(KSPMatSolve(jac->kspschur, ilinkD->X, ilinkD->Y)); 1433 PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1)); 1434 PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->X, ilinkD->Y, NULL)); 1435 PetscCall(MatDenseScatter_Private(ilinkD->sctx, ilinkD->Y, Y, INSERT_VALUES, SCATTER_REVERSE)); 1436 break; 1437 case PC_FIELDSPLIT_SCHUR_FACT_UPPER: 1438 /* [A00 A01; 0 S], suitable for right preconditioning */ 1439 PetscCall(MatDenseScatter_Private(ilinkD->sctx, X, ilinkD->X, INSERT_VALUES, SCATTER_FORWARD)); 1440 PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->X, ilinkD->Y, NULL)); 1441 PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1)); 1442 PetscCall(KSPMatSolve(jac->kspschur, ilinkD->X, ilinkD->Y)); 1443 PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1)); 1444 PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->X, ilinkD->Y, NULL)); 1445 PetscCall(MatMatMult(jac->B, ilinkD->Y, MAT_REUSE_MATRIX, PETSC_DETERMINE, &ilinkA->X)); 1446 PetscCall(MatScale(ilinkA->X, -1.0)); 1447 PetscCall(MatDenseScatter_Private(ilinkA->sctx, X, ilinkA->X, ADD_VALUES, SCATTER_FORWARD)); 1448 PetscCall(MatDenseScatter_Private(ilinkD->sctx, ilinkD->Y, Y, INSERT_VALUES, SCATTER_REVERSE)); 1449 PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->X, ilinkA->Y, NULL)); 1450 PetscCall(KSPMatSolve(kspA, ilinkA->X, ilinkA->Y)); 1451 PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->X, ilinkA->Y, NULL)); 1452 PetscCall(MatDenseScatter_Private(ilinkA->sctx, ilinkA->Y, Y, INSERT_VALUES, SCATTER_REVERSE)); 1453 break; 1454 case PC_FIELDSPLIT_SCHUR_FACT_FULL: 1455 /* [1 0; A10 A00^{-1} 1] [A00 0; 0 S] [1 A00^{-1}A01; 0 1] */ 1456 PetscCall(MatGetSize(jac->B, NULL, &P)); 1457 N = P; 1458 PetscCall(MatDenseScatter_Private(ilinkA->sctx, X, ilinkA->X, INSERT_VALUES, SCATTER_FORWARD)); 1459 PetscCall(PetscLogEventBegin(KSP_Solve_FS_L, kspLower, ilinkA->X, ilinkA->Y, NULL)); 1460 if (kspUpper == kspA) { 1461 PetscCall(PetscObjectQuery((PetscObject)jac->schur, "AinvB", (PetscObject *)&AinvB)); 1462 if (AinvB) { 1463 PetscCall(MatGetSize(AinvB, NULL, &N)); 1464 if (N > P) { // first time PCApply_FieldSplit_Schur() is called 1465 PetscMemType mtype; 1466 Mat C = NULL; 1467 PetscScalar *array; 1468 PetscInt m, M, q, Q, p; 1469 1470 PetscCall(MatGetSize(jac->B, &M, NULL)); 1471 PetscCall(MatGetLocalSize(jac->B, &m, NULL)); 1472 PetscCall(MatGetSize(X, NULL, &Q)); 1473 PetscCall(MatGetLocalSize(X, NULL, &q)); 1474 PetscCall(MatDenseGetArrayAndMemType(AinvB, &array, &mtype)); 1475 if (N != P + Q) { 1476 Mat replace; 1477 1478 PetscCall(MatGetLocalSize(jac->B, NULL, &p)); 1479 if (PetscMemTypeHost(mtype) || (!PetscDefined(HAVE_CUDA) && !PetscDefined(HAVE_HIP))) { 1480 PetscCall(PetscFree(array)); 1481 PetscCall(PetscMalloc1(m * (P + Q), &array)); 1482 PetscCall(MatCreateDense(PetscObjectComm((PetscObject)jac->schur), m, PETSC_DECIDE, M, P + Q, array, &replace)); 1483 } 1484 #if PetscDefined(HAVE_CUDA) 1485 else if (PetscMemTypeCUDA(mtype)) { 1486 PetscCallCUDA(cudaFree(array)); 1487 PetscCallCUDA(cudaMalloc((void **)&array, sizeof(PetscScalar) * m * (P + Q))); 1488 PetscCall(MatCreateDenseCUDA(PetscObjectComm((PetscObject)jac->schur), m, PETSC_DECIDE, M, P + Q, array, &replace)); 1489 } 1490 #endif 1491 #if PetscDefined(HAVE_HIP) 1492 else if (PetscMemTypeHIP(mtype)) { 1493 PetscCallHIP(hipFree(array)); 1494 PetscCallHIP(hipMalloc((void **)&array, sizeof(PetscScalar) * m * (P + Q))); 1495 PetscCall(MatCreateDenseHIP(PetscObjectComm((PetscObject)jac->schur), m, PETSC_DECIDE, M, P + Q, array, &replace)); 1496 } 1497 #endif 1498 PetscCall(MatHeaderReplace(AinvB, &replace)); 1499 } 1500 if (PetscMemTypeHost(mtype) || (!PetscDefined(HAVE_CUDA) && !PetscDefined(HAVE_HIP))) PetscCall(MatCreateDense(PetscObjectComm((PetscObject)jac->schur), m, q, M, Q, array + m * P, &C)); 1501 #if PetscDefined(HAVE_CUDA) 1502 else if (PetscMemTypeCUDA(mtype)) PetscCall(MatCreateDenseCUDA(PetscObjectComm((PetscObject)jac->schur), m, q, M, Q, array + m * P, &C)); 1503 #endif 1504 #if PetscDefined(HAVE_HIP) 1505 else if (PetscMemTypeHIP(mtype)) PetscCall(MatCreateDenseHIP(PetscObjectComm((PetscObject)jac->schur), m, q, M, Q, array + m * P, &C)); 1506 #endif 1507 PetscCall(MatDenseRestoreArrayAndMemType(AinvB, &array)); 1508 PetscCall(MatCopy(ilinkA->X, C, SAME_NONZERO_PATTERN)); 1509 PetscCall(MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user)); 1510 PetscCall(KSPSetOperators(jac->kspschur, jac->schur, jac->schur_user)); 1511 PetscCall(MatCopy(C, ilinkA->Y, SAME_NONZERO_PATTERN)); // retrieve solutions as last columns of the composed Mat 1512 PetscCall(MatDestroy(&C)); 1513 } 1514 } 1515 } 1516 if (N == P) PetscCall(KSPMatSolve(kspLower, ilinkA->X, ilinkA->Y)); 1517 PetscCall(PetscLogEventEnd(KSP_Solve_FS_L, kspLower, ilinkA->X, ilinkA->Y, NULL)); 1518 PetscCall(MatMatMult(jac->C, ilinkA->Y, MAT_REUSE_MATRIX, PETSC_DETERMINE, &ilinkD->X)); 1519 PetscCall(MatScale(ilinkD->X, -1.0)); 1520 PetscCall(MatDenseScatter_Private(ilinkD->sctx, X, ilinkD->X, ADD_VALUES, SCATTER_FORWARD)); 1521 1522 PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->X, ilinkD->Y, NULL)); 1523 PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1)); 1524 PetscCall(KSPMatSolve(jac->kspschur, ilinkD->X, ilinkD->Y)); 1525 PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1)); 1526 PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->X, ilinkD->Y, NULL)); 1527 PetscCall(MatDenseScatter_Private(ilinkD->sctx, ilinkD->Y, Y, INSERT_VALUES, SCATTER_REVERSE)); 1528 1529 if (kspUpper == kspA) { 1530 if (!AinvB) { 1531 PetscCall(MatMatMult(jac->B, ilinkD->Y, MAT_REUSE_MATRIX, PETSC_DETERMINE, &ilinkA->Y)); 1532 PetscCall(MatAXPY(ilinkA->X, -1.0, ilinkA->Y, SAME_NONZERO_PATTERN)); 1533 PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->X, ilinkA->Y, NULL)); 1534 PetscCall(KSPMatSolve(kspA, ilinkA->X, ilinkA->Y)); 1535 PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->X, ilinkA->Y, NULL)); 1536 } else { 1537 PetscCall(MatMatMult(AinvB, ilinkD->Y, MAT_REUSE_MATRIX, PETSC_DETERMINE, &ilinkA->X)); 1538 PetscCall(MatAXPY(ilinkA->Y, 1.0, ilinkA->X, SAME_NONZERO_PATTERN)); 1539 } 1540 } else { 1541 PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->X, ilinkA->Y, NULL)); 1542 PetscCall(KSPMatSolve(kspA, ilinkA->X, ilinkA->Y)); 1543 PetscCall(MatMatMult(jac->B, ilinkD->Y, MAT_REUSE_MATRIX, PETSC_DETERMINE, &ilinkA->X)); 1544 if (!ilinkA->Z) PetscCall(MatDuplicate(ilinkA->X, MAT_DO_NOT_COPY_VALUES, &ilinkA->Z)); 1545 PetscCall(PetscLogEventBegin(KSP_Solve_FS_U, kspUpper, ilinkA->X, ilinkA->Z, NULL)); 1546 PetscCall(KSPMatSolve(kspUpper, ilinkA->X, ilinkA->Z)); 1547 PetscCall(PetscLogEventEnd(KSP_Solve_FS_U, kspUpper, ilinkA->X, ilinkA->Z, NULL)); 1548 PetscCall(MatAXPY(ilinkA->Y, -1.0, ilinkA->Z, SAME_NONZERO_PATTERN)); 1549 } 1550 PetscCall(MatDenseScatter_Private(ilinkA->sctx, ilinkA->Y, Y, INSERT_VALUES, SCATTER_REVERSE)); 1551 } 1552 PetscFunctionReturn(PETSC_SUCCESS); 1553 } 1554 1555 static PetscErrorCode PCApplyTranspose_FieldSplit_Schur(PC pc, Vec x, Vec y) 1556 { 1557 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 1558 PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next; 1559 KSP kspA = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper; 1560 1561 PetscFunctionBegin; 1562 switch (jac->schurfactorization) { 1563 case PC_FIELDSPLIT_SCHUR_FACT_DIAG: 1564 /* [A00 0; 0 -S], positive definite, suitable for MINRES */ 1565 PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD)); 1566 PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD)); 1567 PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD)); 1568 PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL)); 1569 PetscCall(KSPSolveTranspose(kspA, ilinkA->x, ilinkA->y)); 1570 PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y)); 1571 PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL)); 1572 PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1573 PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD)); 1574 PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL)); 1575 PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1)); 1576 PetscCall(KSPSolveTranspose(jac->kspschur, ilinkD->x, ilinkD->y)); 1577 PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1)); 1578 PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y)); 1579 PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL)); 1580 PetscCall(VecScale(ilinkD->y, jac->schurscale)); 1581 PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1582 PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1583 PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1584 break; 1585 case PC_FIELDSPLIT_SCHUR_FACT_UPPER: 1586 PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD)); 1587 PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD)); 1588 PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL)); 1589 PetscCall(KSPSolveTranspose(kspA, ilinkA->x, ilinkA->y)); 1590 PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y)); 1591 PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL)); 1592 PetscCall(MatMultTranspose(jac->B, ilinkA->y, ilinkD->x)); 1593 PetscCall(VecScale(ilinkD->x, -1.)); 1594 PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD)); 1595 PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1596 PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD)); 1597 PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL)); 1598 PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1)); 1599 PetscCall(KSPSolveTranspose(jac->kspschur, ilinkD->x, ilinkD->y)); 1600 PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1)); 1601 PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y)); 1602 PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL)); 1603 PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1604 PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1605 PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1606 break; 1607 case PC_FIELDSPLIT_SCHUR_FACT_LOWER: 1608 PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD)); 1609 PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD)); 1610 PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL)); 1611 PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1)); 1612 PetscCall(KSPSolveTranspose(jac->kspschur, ilinkD->x, ilinkD->y)); 1613 PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1)); 1614 PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y)); 1615 PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL)); 1616 PetscCall(MatMultTranspose(jac->C, ilinkD->y, ilinkA->x)); 1617 PetscCall(VecScale(ilinkA->x, -1.)); 1618 PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, ADD_VALUES, SCATTER_FORWARD)); 1619 PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1620 PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, ADD_VALUES, SCATTER_FORWARD)); 1621 PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL)); 1622 PetscCall(KSPSolveTranspose(kspA, ilinkA->x, ilinkA->y)); 1623 PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y)); 1624 PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL)); 1625 PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1626 PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1627 PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1628 break; 1629 case PC_FIELDSPLIT_SCHUR_FACT_FULL: 1630 PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD)); 1631 PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD)); 1632 PetscCall(PetscLogEventBegin(KSP_Solve_FS_U, kspUpper, ilinkA->x, ilinkA->y, NULL)); 1633 PetscCall(KSPSolveTranspose(kspUpper, ilinkA->x, ilinkA->y)); 1634 PetscCall(KSPCheckSolve(kspUpper, pc, ilinkA->y)); 1635 PetscCall(PetscLogEventEnd(KSP_Solve_FS_U, kspUpper, ilinkA->x, ilinkA->y, NULL)); 1636 PetscCall(MatMultTranspose(jac->B, ilinkA->y, ilinkD->x)); 1637 PetscCall(VecScale(ilinkD->x, -1.0)); 1638 PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD)); 1639 PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD)); 1640 1641 PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL)); 1642 PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1)); 1643 PetscCall(KSPSolveTranspose(jac->kspschur, ilinkD->x, ilinkD->y)); 1644 PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1)); 1645 PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y)); 1646 PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL)); 1647 PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1648 1649 if (kspLower == kspA) { 1650 PetscCall(MatMultTranspose(jac->C, ilinkD->y, ilinkA->y)); 1651 PetscCall(VecAXPY(ilinkA->x, -1.0, ilinkA->y)); 1652 PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL)); 1653 PetscCall(KSPSolveTranspose(kspA, ilinkA->x, ilinkA->y)); 1654 PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y)); 1655 PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL)); 1656 } else { 1657 PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL)); 1658 PetscCall(KSPSolveTranspose(kspA, ilinkA->x, ilinkA->y)); 1659 PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y)); 1660 PetscCall(MatMultTranspose(jac->C, ilinkD->y, ilinkA->x)); 1661 PetscCall(PetscLogEventBegin(KSP_Solve_FS_L, kspLower, ilinkA->x, ilinkA->z, NULL)); 1662 PetscCall(KSPSolveTranspose(kspLower, ilinkA->x, ilinkA->z)); 1663 PetscCall(KSPCheckSolve(kspLower, pc, ilinkA->z)); 1664 PetscCall(PetscLogEventEnd(KSP_Solve_FS_L, kspLower, ilinkA->x, ilinkA->z, NULL)); 1665 PetscCall(VecAXPY(ilinkA->y, -1.0, ilinkA->z)); 1666 } 1667 PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1668 PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1669 PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1670 } 1671 PetscFunctionReturn(PETSC_SUCCESS); 1672 } 1673 1674 #define FieldSplitSplitSolveAdd(ilink, xx, yy) \ 1675 ((PetscErrorCode)(VecScatterBegin(ilink->sctx, xx, ilink->x, INSERT_VALUES, SCATTER_FORWARD) || VecScatterEnd(ilink->sctx, xx, ilink->x, INSERT_VALUES, SCATTER_FORWARD) || PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL) || \ 1676 KSPSolve(ilink->ksp, ilink->x, ilink->y) || KSPCheckSolve(ilink->ksp, pc, ilink->y) || PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL) || VecScatterBegin(ilink->sctx, ilink->y, yy, ADD_VALUES, SCATTER_REVERSE) || \ 1677 VecScatterEnd(ilink->sctx, ilink->y, yy, ADD_VALUES, SCATTER_REVERSE))) 1678 1679 static PetscErrorCode PCApply_FieldSplit(PC pc, Vec x, Vec y) 1680 { 1681 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 1682 PC_FieldSplitLink ilink = jac->head; 1683 PetscInt cnt, bs; 1684 1685 PetscFunctionBegin; 1686 if (jac->type == PC_COMPOSITE_ADDITIVE) { 1687 PetscBool matnest; 1688 1689 PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATNEST, &matnest)); 1690 if (jac->defaultsplit && !matnest) { 1691 PetscCall(VecGetBlockSize(x, &bs)); 1692 PetscCheck(jac->bs <= 0 || bs == jac->bs, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Blocksize of x vector %" PetscInt_FMT " does not match fieldsplit blocksize %" PetscInt_FMT, bs, jac->bs); 1693 PetscCall(VecGetBlockSize(y, &bs)); 1694 PetscCheck(jac->bs <= 0 || bs == jac->bs, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Blocksize of y vector %" PetscInt_FMT " does not match fieldsplit blocksize %" PetscInt_FMT, bs, jac->bs); 1695 PetscCall(VecStrideGatherAll(x, jac->x, INSERT_VALUES)); 1696 while (ilink) { 1697 PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL)); 1698 PetscCall(KSPSolve(ilink->ksp, ilink->x, ilink->y)); 1699 PetscCall(KSPCheckSolve(ilink->ksp, pc, ilink->y)); 1700 PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL)); 1701 ilink = ilink->next; 1702 } 1703 PetscCall(VecStrideScatterAll(jac->y, y, INSERT_VALUES)); 1704 } else { 1705 PetscCall(VecSet(y, 0.0)); 1706 while (ilink) { 1707 PetscCall(FieldSplitSplitSolveAdd(ilink, x, y)); 1708 ilink = ilink->next; 1709 } 1710 } 1711 } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE && jac->nsplits == 2) { 1712 PetscCall(VecSet(y, 0.0)); 1713 /* solve on first block for first block variables */ 1714 PetscCall(VecScatterBegin(ilink->sctx, x, ilink->x, INSERT_VALUES, SCATTER_FORWARD)); 1715 PetscCall(VecScatterEnd(ilink->sctx, x, ilink->x, INSERT_VALUES, SCATTER_FORWARD)); 1716 PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL)); 1717 PetscCall(KSPSolve(ilink->ksp, ilink->x, ilink->y)); 1718 PetscCall(KSPCheckSolve(ilink->ksp, pc, ilink->y)); 1719 PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL)); 1720 PetscCall(VecScatterBegin(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE)); 1721 PetscCall(VecScatterEnd(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE)); 1722 1723 /* compute the residual only onto second block variables using first block variables */ 1724 PetscCall(MatMult(jac->Afield[1], ilink->y, ilink->next->x)); 1725 ilink = ilink->next; 1726 PetscCall(VecScale(ilink->x, -1.0)); 1727 PetscCall(VecScatterBegin(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD)); 1728 PetscCall(VecScatterEnd(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD)); 1729 1730 /* solve on second block variables */ 1731 PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL)); 1732 PetscCall(KSPSolve(ilink->ksp, ilink->x, ilink->y)); 1733 PetscCall(KSPCheckSolve(ilink->ksp, pc, ilink->y)); 1734 PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL)); 1735 PetscCall(VecScatterBegin(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE)); 1736 PetscCall(VecScatterEnd(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE)); 1737 } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 1738 if (!jac->w1) { 1739 PetscCall(VecDuplicate(x, &jac->w1)); 1740 PetscCall(VecDuplicate(x, &jac->w2)); 1741 } 1742 PetscCall(VecSet(y, 0.0)); 1743 PetscCall(FieldSplitSplitSolveAdd(ilink, x, y)); 1744 cnt = 1; 1745 while (ilink->next) { 1746 ilink = ilink->next; 1747 /* compute the residual only over the part of the vector needed */ 1748 PetscCall(MatMult(jac->Afield[cnt++], y, ilink->x)); 1749 PetscCall(VecScale(ilink->x, -1.0)); 1750 PetscCall(VecScatterBegin(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD)); 1751 PetscCall(VecScatterEnd(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD)); 1752 PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL)); 1753 PetscCall(KSPSolve(ilink->ksp, ilink->x, ilink->y)); 1754 PetscCall(KSPCheckSolve(ilink->ksp, pc, ilink->y)); 1755 PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL)); 1756 PetscCall(VecScatterBegin(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE)); 1757 PetscCall(VecScatterEnd(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE)); 1758 } 1759 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 1760 cnt -= 2; 1761 while (ilink->previous) { 1762 ilink = ilink->previous; 1763 /* compute the residual only over the part of the vector needed */ 1764 PetscCall(MatMult(jac->Afield[cnt--], y, ilink->x)); 1765 PetscCall(VecScale(ilink->x, -1.0)); 1766 PetscCall(VecScatterBegin(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD)); 1767 PetscCall(VecScatterEnd(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD)); 1768 PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL)); 1769 PetscCall(KSPSolve(ilink->ksp, ilink->x, ilink->y)); 1770 PetscCall(KSPCheckSolve(ilink->ksp, pc, ilink->y)); 1771 PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL)); 1772 PetscCall(VecScatterBegin(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE)); 1773 PetscCall(VecScatterEnd(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE)); 1774 } 1775 } 1776 } else SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Unsupported or unknown composition %d", (int)jac->type); 1777 PetscFunctionReturn(PETSC_SUCCESS); 1778 } 1779 1780 static PetscErrorCode PCMatApply_FieldSplit(PC pc, Mat X, Mat Y) 1781 { 1782 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 1783 PC_FieldSplitLink ilink = jac->head; 1784 PetscInt cnt; 1785 1786 PetscFunctionBegin; 1787 /* create working matrices with the correct number of columns */ 1788 PetscCall(PCFieldSplitCreateWorkMats_Private(pc, X)); 1789 if (jac->type == PC_COMPOSITE_ADDITIVE) { 1790 PetscCall(MatZeroEntries(Y)); 1791 while (ilink) { 1792 PetscCall(MatDenseScatter_Private(ilink->sctx, X, ilink->X, INSERT_VALUES, SCATTER_FORWARD)); 1793 PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL)); 1794 PetscCall(KSPMatSolve(ilink->ksp, ilink->X, ilink->Y)); 1795 PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL)); 1796 PetscCall(MatDenseScatter_Private(ilink->sctx, ilink->Y, Y, ADD_VALUES, SCATTER_REVERSE)); 1797 ilink = ilink->next; 1798 } 1799 } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE && jac->nsplits == 2) { 1800 PetscCall(MatZeroEntries(Y)); 1801 PetscCall(MatDenseScatter_Private(ilink->sctx, X, ilink->X, INSERT_VALUES, SCATTER_FORWARD)); 1802 PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL)); 1803 PetscCall(KSPMatSolve(ilink->ksp, ilink->X, ilink->Y)); 1804 PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL)); 1805 PetscCall(MatDenseScatter_Private(ilink->sctx, ilink->Y, Y, ADD_VALUES, SCATTER_REVERSE)); 1806 1807 /* compute the residual only onto second block variables using first block variables */ 1808 PetscCall(MatMatMult(jac->Afield[1], ilink->Y, MAT_REUSE_MATRIX, PETSC_DETERMINE, &ilink->next->X)); 1809 ilink = ilink->next; 1810 PetscCall(MatScale(ilink->X, -1.0)); 1811 PetscCall(MatDenseScatter_Private(ilink->sctx, X, ilink->X, ADD_VALUES, SCATTER_FORWARD)); 1812 1813 /* solve on second block variables */ 1814 PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL)); 1815 PetscCall(KSPMatSolve(ilink->ksp, ilink->X, ilink->Y)); 1816 PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL)); 1817 PetscCall(MatDenseScatter_Private(ilink->sctx, ilink->Y, Y, ADD_VALUES, SCATTER_REVERSE)); 1818 } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 1819 /* general multiplicative with any number of splits */ 1820 PetscCall(MatZeroEntries(Y)); 1821 /* first split */ 1822 PetscCall(MatDenseScatter_Private(ilink->sctx, X, ilink->X, INSERT_VALUES, SCATTER_FORWARD)); 1823 PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL)); 1824 PetscCall(KSPMatSolve(ilink->ksp, ilink->X, ilink->Y)); 1825 PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL)); 1826 PetscCall(MatDenseScatter_Private(ilink->sctx, ilink->Y, Y, ADD_VALUES, SCATTER_REVERSE)); 1827 cnt = 1; 1828 /* forward sweep */ 1829 while (ilink->next) { 1830 ilink = ilink->next; 1831 /* compute the residual only over the part of the vector needed */ 1832 PetscCall(MatMatMult(jac->Afield[cnt++], Y, MAT_REUSE_MATRIX, PETSC_DETERMINE, &ilink->X)); 1833 PetscCall(MatScale(ilink->X, -1.0)); 1834 PetscCall(MatDenseScatter_Private(ilink->sctx, X, ilink->X, ADD_VALUES, SCATTER_FORWARD)); 1835 PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL)); 1836 PetscCall(KSPMatSolve(ilink->ksp, ilink->X, ilink->Y)); 1837 PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL)); 1838 PetscCall(MatDenseScatter_Private(ilink->sctx, ilink->Y, Y, ADD_VALUES, SCATTER_REVERSE)); 1839 } 1840 /* backward sweep for symmetric multiplicative */ 1841 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 1842 cnt -= 2; 1843 while (ilink->previous) { 1844 ilink = ilink->previous; 1845 /* compute the residual only over the part of the vector needed */ 1846 PetscCall(MatMatMult(jac->Afield[cnt--], Y, MAT_REUSE_MATRIX, PETSC_DETERMINE, &ilink->X)); 1847 PetscCall(MatScale(ilink->X, -1.0)); 1848 PetscCall(MatDenseScatter_Private(ilink->sctx, X, ilink->X, ADD_VALUES, SCATTER_FORWARD)); 1849 PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL)); 1850 PetscCall(KSPMatSolve(ilink->ksp, ilink->X, ilink->Y)); 1851 PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL)); 1852 PetscCall(MatDenseScatter_Private(ilink->sctx, ilink->Y, Y, ADD_VALUES, SCATTER_REVERSE)); 1853 } 1854 } 1855 } else SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "PCMatApply() not implemented for this fieldsplit type"); 1856 PetscFunctionReturn(PETSC_SUCCESS); 1857 } 1858 1859 static PetscErrorCode PCApply_FieldSplit_GKB(PC pc, Vec x, Vec y) 1860 { 1861 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 1862 PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next; 1863 KSP ksp = ilinkA->ksp; 1864 Vec u, v, Hu, d, work1, work2; 1865 PetscScalar alpha, z, nrmz2, *vecz; 1866 PetscReal lowbnd, nu, beta; 1867 PetscInt j, iterGKB; 1868 1869 PetscFunctionBegin; 1870 PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD)); 1871 PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD)); 1872 PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD)); 1873 PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD)); 1874 1875 u = jac->u; 1876 v = jac->v; 1877 Hu = jac->Hu; 1878 d = jac->d; 1879 work1 = jac->w1; 1880 work2 = jac->w2; 1881 vecz = jac->vecz; 1882 1883 /* Change RHS to comply with matrix regularization H = A + nu*B*B' */ 1884 /* Add q = q + nu*B*b */ 1885 if (jac->gkbnu) { 1886 nu = jac->gkbnu; 1887 PetscCall(VecScale(ilinkD->x, jac->gkbnu)); 1888 PetscCall(MatMultAdd(jac->B, ilinkD->x, ilinkA->x, ilinkA->x)); /* q = q + nu*B*b */ 1889 } else { 1890 /* Situation when no augmented Lagrangian is used. Then we set inner */ 1891 /* matrix N = I in [Ar13], and thus nu = 1. */ 1892 nu = 1; 1893 } 1894 1895 /* Transform rhs from [q,tilde{b}] to [0,b] */ 1896 PetscCall(PetscLogEventBegin(ilinkA->event, ksp, ilinkA->x, ilinkA->y, NULL)); 1897 PetscCall(KSPSolve(ksp, ilinkA->x, ilinkA->y)); 1898 PetscCall(KSPCheckSolve(ksp, pc, ilinkA->y)); 1899 PetscCall(PetscLogEventEnd(ilinkA->event, ksp, ilinkA->x, ilinkA->y, NULL)); 1900 PetscCall(MatMultHermitianTranspose(jac->B, ilinkA->y, work1)); 1901 PetscCall(VecAXPBY(work1, 1.0 / nu, -1.0, ilinkD->x)); /* c = b - B'*x */ 1902 1903 /* First step of algorithm */ 1904 PetscCall(VecNorm(work1, NORM_2, &beta)); /* beta = sqrt(nu*c'*c)*/ 1905 KSPCheckDot(ksp, beta); 1906 beta = PetscSqrtReal(nu) * beta; 1907 PetscCall(VecAXPBY(v, nu / beta, 0.0, work1)); /* v = nu/beta *c */ 1908 PetscCall(MatMult(jac->B, v, work2)); /* u = H^{-1}*B*v */ 1909 PetscCall(PetscLogEventBegin(ilinkA->event, ksp, work2, u, NULL)); 1910 PetscCall(KSPSolve(ksp, work2, u)); 1911 PetscCall(KSPCheckSolve(ksp, pc, u)); 1912 PetscCall(PetscLogEventEnd(ilinkA->event, ksp, work2, u, NULL)); 1913 PetscCall(MatMult(jac->H, u, Hu)); /* alpha = u'*H*u */ 1914 PetscCall(VecDot(Hu, u, &alpha)); 1915 KSPCheckDot(ksp, alpha); 1916 PetscCheck(PetscRealPart(alpha) > 0.0, PETSC_COMM_SELF, PETSC_ERR_NOT_CONVERGED, "GKB preconditioner diverged, H is not positive definite"); 1917 alpha = PetscSqrtReal(PetscAbsScalar(alpha)); 1918 PetscCall(VecScale(u, 1.0 / alpha)); 1919 PetscCall(VecAXPBY(d, 1.0 / alpha, 0.0, v)); /* v = nu/beta *c */ 1920 1921 z = beta / alpha; 1922 vecz[1] = z; 1923 1924 /* Computation of first iterate x(1) and p(1) */ 1925 PetscCall(VecAXPY(ilinkA->y, z, u)); 1926 PetscCall(VecCopy(d, ilinkD->y)); 1927 PetscCall(VecScale(ilinkD->y, -z)); 1928 1929 iterGKB = 1; 1930 lowbnd = 2 * jac->gkbtol; 1931 if (jac->gkbmonitor) PetscCall(PetscViewerASCIIPrintf(jac->gkbviewer, "%3" PetscInt_FMT " GKB Lower bound estimate %14.12e\n", iterGKB, (double)lowbnd)); 1932 1933 while (iterGKB < jac->gkbmaxit && lowbnd > jac->gkbtol) { 1934 iterGKB += 1; 1935 PetscCall(MatMultHermitianTranspose(jac->B, u, work1)); /* v <- nu*(B'*u-alpha/nu*v) */ 1936 PetscCall(VecAXPBY(v, nu, -alpha, work1)); 1937 PetscCall(VecNorm(v, NORM_2, &beta)); /* beta = sqrt(nu)*v'*v */ 1938 beta = beta / PetscSqrtReal(nu); 1939 PetscCall(VecScale(v, 1.0 / beta)); 1940 PetscCall(MatMult(jac->B, v, work2)); /* u <- H^{-1}*(B*v-beta*H*u) */ 1941 PetscCall(MatMult(jac->H, u, Hu)); 1942 PetscCall(VecAXPY(work2, -beta, Hu)); 1943 PetscCall(PetscLogEventBegin(ilinkA->event, ksp, work2, u, NULL)); 1944 PetscCall(KSPSolve(ksp, work2, u)); 1945 PetscCall(KSPCheckSolve(ksp, pc, u)); 1946 PetscCall(PetscLogEventEnd(ilinkA->event, ksp, work2, u, NULL)); 1947 PetscCall(MatMult(jac->H, u, Hu)); /* alpha = u'*H*u */ 1948 PetscCall(VecDot(Hu, u, &alpha)); 1949 KSPCheckDot(ksp, alpha); 1950 PetscCheck(PetscRealPart(alpha) > 0.0, PETSC_COMM_SELF, PETSC_ERR_NOT_CONVERGED, "GKB preconditioner diverged, H is not positive definite"); 1951 alpha = PetscSqrtReal(PetscAbsScalar(alpha)); 1952 PetscCall(VecScale(u, 1.0 / alpha)); 1953 1954 z = -beta / alpha * z; /* z <- beta/alpha*z */ 1955 vecz[0] = z; 1956 1957 /* Computation of new iterate x(i+1) and p(i+1) */ 1958 PetscCall(VecAXPBY(d, 1.0 / alpha, -beta / alpha, v)); /* d = (v-beta*d)/alpha */ 1959 PetscCall(VecAXPY(ilinkA->y, z, u)); /* r = r + z*u */ 1960 PetscCall(VecAXPY(ilinkD->y, -z, d)); /* p = p - z*d */ 1961 PetscCall(MatMult(jac->H, ilinkA->y, Hu)); /* ||u||_H = u'*H*u */ 1962 PetscCall(VecDot(Hu, ilinkA->y, &nrmz2)); 1963 1964 /* Compute Lower Bound estimate */ 1965 if (iterGKB > jac->gkbdelay) { 1966 lowbnd = 0.0; 1967 for (j = 0; j < jac->gkbdelay; j++) lowbnd += PetscAbsScalar(vecz[j] * vecz[j]); 1968 lowbnd = PetscSqrtReal(lowbnd / PetscAbsScalar(nrmz2)); 1969 } 1970 1971 for (j = 0; j < jac->gkbdelay - 1; j++) vecz[jac->gkbdelay - j - 1] = vecz[jac->gkbdelay - j - 2]; 1972 if (jac->gkbmonitor) PetscCall(PetscViewerASCIIPrintf(jac->gkbviewer, "%3" PetscInt_FMT " GKB Lower bound estimate %14.12e\n", iterGKB, (double)lowbnd)); 1973 } 1974 1975 PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1976 PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1977 PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1978 PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1979 PetscFunctionReturn(PETSC_SUCCESS); 1980 } 1981 1982 #define FieldSplitSplitSolveAddTranspose(ilink, xx, yy) \ 1983 ((PetscErrorCode)(VecScatterBegin(ilink->sctx, xx, ilink->y, INSERT_VALUES, SCATTER_FORWARD) || VecScatterEnd(ilink->sctx, xx, ilink->y, INSERT_VALUES, SCATTER_FORWARD) || PetscLogEventBegin(ilink->event, ilink->ksp, ilink->y, ilink->x, NULL) || \ 1984 KSPSolveTranspose(ilink->ksp, ilink->y, ilink->x) || KSPCheckSolve(ilink->ksp, pc, ilink->x) || PetscLogEventEnd(ilink->event, ilink->ksp, ilink->y, ilink->x, NULL) || VecScatterBegin(ilink->sctx, ilink->x, yy, ADD_VALUES, SCATTER_REVERSE) || \ 1985 VecScatterEnd(ilink->sctx, ilink->x, yy, ADD_VALUES, SCATTER_REVERSE))) 1986 1987 static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc, Vec x, Vec y) 1988 { 1989 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 1990 PC_FieldSplitLink ilink = jac->head; 1991 PetscInt bs; 1992 1993 PetscFunctionBegin; 1994 if (jac->type == PC_COMPOSITE_ADDITIVE) { 1995 PetscBool matnest; 1996 1997 PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATNEST, &matnest)); 1998 if (jac->defaultsplit && !matnest) { 1999 PetscCall(VecGetBlockSize(x, &bs)); 2000 PetscCheck(jac->bs <= 0 || bs == jac->bs, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Blocksize of x vector %" PetscInt_FMT " does not match fieldsplit blocksize %" PetscInt_FMT, bs, jac->bs); 2001 PetscCall(VecGetBlockSize(y, &bs)); 2002 PetscCheck(jac->bs <= 0 || bs == jac->bs, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Blocksize of y vector %" PetscInt_FMT " does not match fieldsplit blocksize %" PetscInt_FMT, bs, jac->bs); 2003 PetscCall(VecStrideGatherAll(x, jac->x, INSERT_VALUES)); 2004 while (ilink) { 2005 PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL)); 2006 PetscCall(KSPSolveTranspose(ilink->ksp, ilink->x, ilink->y)); 2007 PetscCall(KSPCheckSolve(ilink->ksp, pc, ilink->y)); 2008 PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL)); 2009 ilink = ilink->next; 2010 } 2011 PetscCall(VecStrideScatterAll(jac->y, y, INSERT_VALUES)); 2012 } else { 2013 PetscCall(VecSet(y, 0.0)); 2014 while (ilink) { 2015 PetscCall(FieldSplitSplitSolveAddTranspose(ilink, x, y)); 2016 ilink = ilink->next; 2017 } 2018 } 2019 } else { 2020 if (!jac->w1) { 2021 PetscCall(VecDuplicate(x, &jac->w1)); 2022 PetscCall(VecDuplicate(x, &jac->w2)); 2023 } 2024 PetscCall(VecSet(y, 0.0)); 2025 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 2026 PetscCall(FieldSplitSplitSolveAddTranspose(ilink, x, y)); 2027 while (ilink->next) { 2028 ilink = ilink->next; 2029 PetscCall(MatMultTranspose(pc->mat, y, jac->w1)); 2030 PetscCall(VecWAXPY(jac->w2, -1.0, jac->w1, x)); 2031 PetscCall(FieldSplitSplitSolveAddTranspose(ilink, jac->w2, y)); 2032 } 2033 while (ilink->previous) { 2034 ilink = ilink->previous; 2035 PetscCall(MatMultTranspose(pc->mat, y, jac->w1)); 2036 PetscCall(VecWAXPY(jac->w2, -1.0, jac->w1, x)); 2037 PetscCall(FieldSplitSplitSolveAddTranspose(ilink, jac->w2, y)); 2038 } 2039 } else { 2040 while (ilink->next) { /* get to last entry in linked list */ 2041 ilink = ilink->next; 2042 } 2043 PetscCall(FieldSplitSplitSolveAddTranspose(ilink, x, y)); 2044 while (ilink->previous) { 2045 ilink = ilink->previous; 2046 PetscCall(MatMultTranspose(pc->mat, y, jac->w1)); 2047 PetscCall(VecWAXPY(jac->w2, -1.0, jac->w1, x)); 2048 PetscCall(FieldSplitSplitSolveAddTranspose(ilink, jac->w2, y)); 2049 } 2050 } 2051 } 2052 PetscFunctionReturn(PETSC_SUCCESS); 2053 } 2054 2055 static PetscErrorCode PCReset_FieldSplit(PC pc) 2056 { 2057 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 2058 PC_FieldSplitLink ilink = jac->head, next; 2059 2060 PetscFunctionBegin; 2061 while (ilink) { 2062 PetscCall(KSPDestroy(&ilink->ksp)); 2063 PetscCall(VecDestroy(&ilink->x)); 2064 PetscCall(VecDestroy(&ilink->y)); 2065 PetscCall(VecDestroy(&ilink->z)); 2066 PetscCall(MatDestroy(&ilink->X)); 2067 PetscCall(MatDestroy(&ilink->Y)); 2068 PetscCall(MatDestroy(&ilink->Z)); 2069 PetscCall(VecScatterDestroy(&ilink->sctx)); 2070 PetscCall(ISDestroy(&ilink->is)); 2071 PetscCall(ISDestroy(&ilink->is_col)); 2072 PetscCall(PetscFree(ilink->splitname)); 2073 PetscCall(PetscFree(ilink->fields)); 2074 PetscCall(PetscFree(ilink->fields_col)); 2075 next = ilink->next; 2076 PetscCall(PetscFree(ilink)); 2077 ilink = next; 2078 } 2079 jac->head = NULL; 2080 PetscCall(PetscFree2(jac->x, jac->y)); 2081 if (jac->mat && jac->mat != jac->pmat) { 2082 PetscCall(MatDestroyMatrices(jac->nsplits, &jac->mat)); 2083 } else if (jac->mat) { 2084 jac->mat = NULL; 2085 } 2086 if (jac->pmat) PetscCall(MatDestroyMatrices(jac->nsplits, &jac->pmat)); 2087 if (jac->Afield) PetscCall(MatDestroyMatrices(jac->nsplits, &jac->Afield)); 2088 jac->nsplits = 0; 2089 PetscCall(VecDestroy(&jac->w1)); 2090 PetscCall(VecDestroy(&jac->w2)); 2091 if (jac->schur) PetscCall(PetscObjectCompose((PetscObject)jac->schur, "AinvB", NULL)); 2092 PetscCall(MatDestroy(&jac->schur)); 2093 PetscCall(MatDestroy(&jac->schurp)); 2094 PetscCall(MatDestroy(&jac->schur_user)); 2095 PetscCall(KSPDestroy(&jac->kspschur)); 2096 PetscCall(KSPDestroy(&jac->kspupper)); 2097 PetscCall(MatDestroy(&jac->B)); 2098 PetscCall(MatDestroy(&jac->C)); 2099 PetscCall(MatDestroy(&jac->H)); 2100 PetscCall(VecDestroy(&jac->u)); 2101 PetscCall(VecDestroy(&jac->v)); 2102 PetscCall(VecDestroy(&jac->Hu)); 2103 PetscCall(VecDestroy(&jac->d)); 2104 PetscCall(PetscFree(jac->vecz)); 2105 PetscCall(PetscViewerDestroy(&jac->gkbviewer)); 2106 jac->isrestrict = PETSC_FALSE; 2107 PetscFunctionReturn(PETSC_SUCCESS); 2108 } 2109 2110 static PetscErrorCode PCDestroy_FieldSplit(PC pc) 2111 { 2112 PetscFunctionBegin; 2113 PetscCall(PCReset_FieldSplit(pc)); 2114 PetscCall(PetscFree(pc->data)); 2115 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", NULL)); 2116 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetFields_C", NULL)); 2117 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetIS_C", NULL)); 2118 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetType_C", NULL)); 2119 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetBlockSize_C", NULL)); 2120 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitRestrictIS_C", NULL)); 2121 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSchurGetSubKSP_C", NULL)); 2122 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", NULL)); 2123 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBTol_C", NULL)); 2124 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBMaxit_C", NULL)); 2125 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBNu_C", NULL)); 2126 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBDelay_C", NULL)); 2127 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurPre_C", NULL)); 2128 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSchurPre_C", NULL)); 2129 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurFactType_C", NULL)); 2130 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurScale_C", NULL)); 2131 PetscFunctionReturn(PETSC_SUCCESS); 2132 } 2133 2134 static PetscErrorCode PCSetFromOptions_FieldSplit(PC pc, PetscOptionItems PetscOptionsObject) 2135 { 2136 PetscInt bs; 2137 PetscBool flg; 2138 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 2139 PCCompositeType ctype; 2140 2141 PetscFunctionBegin; 2142 PetscOptionsHeadBegin(PetscOptionsObject, "FieldSplit options"); 2143 PetscCall(PetscOptionsBool("-pc_fieldsplit_dm_splits", "Whether to use DMCreateFieldDecomposition() for splits", "PCFieldSplitSetDMSplits", jac->dm_splits, &jac->dm_splits, NULL)); 2144 PetscCall(PetscOptionsInt("-pc_fieldsplit_block_size", "Blocksize that defines number of fields", "PCFieldSplitSetBlockSize", jac->bs, &bs, &flg)); 2145 if (flg) PetscCall(PCFieldSplitSetBlockSize(pc, bs)); 2146 jac->diag_use_amat = pc->useAmat; 2147 PetscCall(PetscOptionsBool("-pc_fieldsplit_diag_use_amat", "Use Amat (not Pmat) to extract diagonal fieldsplit blocks", "PCFieldSplitSetDiagUseAmat", jac->diag_use_amat, &jac->diag_use_amat, NULL)); 2148 jac->offdiag_use_amat = pc->useAmat; 2149 PetscCall(PetscOptionsBool("-pc_fieldsplit_off_diag_use_amat", "Use Amat (not Pmat) to extract off-diagonal fieldsplit blocks", "PCFieldSplitSetOffDiagUseAmat", jac->offdiag_use_amat, &jac->offdiag_use_amat, NULL)); 2150 PetscCall(PetscOptionsBool("-pc_fieldsplit_detect_saddle_point", "Form 2-way split by detecting zero diagonal entries", "PCFieldSplitSetDetectSaddlePoint", jac->detect, &jac->detect, NULL)); 2151 PetscCall(PCFieldSplitSetDetectSaddlePoint(pc, jac->detect)); /* Sets split type and Schur PC type */ 2152 PetscCall(PetscOptionsEnum("-pc_fieldsplit_type", "Type of composition", "PCFieldSplitSetType", PCCompositeTypes, (PetscEnum)jac->type, (PetscEnum *)&ctype, &flg)); 2153 if (flg) PetscCall(PCFieldSplitSetType(pc, ctype)); 2154 /* Only setup fields once */ 2155 if (jac->bs > 0 && jac->nsplits == 0) { 2156 /* only allow user to set fields from command line. 2157 otherwise user can set them in PCFieldSplitSetDefaults() */ 2158 PetscCall(PCFieldSplitSetRuntimeSplits_Private(pc)); 2159 if (jac->splitdefined) PetscCall(PetscInfo(pc, "Splits defined using the options database\n")); 2160 } 2161 if (jac->type == PC_COMPOSITE_SCHUR) { 2162 PetscCall(PetscOptionsGetEnum(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_fieldsplit_schur_factorization_type", PCFieldSplitSchurFactTypes, (PetscEnum *)&jac->schurfactorization, &flg)); 2163 if (flg) PetscCall(PetscInfo(pc, "Deprecated use of -pc_fieldsplit_schur_factorization_type\n")); 2164 PetscCall(PetscOptionsEnum("-pc_fieldsplit_schur_fact_type", "Which off-diagonal parts of the block factorization to use", "PCFieldSplitSetSchurFactType", PCFieldSplitSchurFactTypes, (PetscEnum)jac->schurfactorization, (PetscEnum *)&jac->schurfactorization, NULL)); 2165 PetscCall(PetscOptionsEnum("-pc_fieldsplit_schur_precondition", "How to build preconditioner for Schur complement", "PCFieldSplitSetSchurPre", PCFieldSplitSchurPreTypes, (PetscEnum)jac->schurpre, (PetscEnum *)&jac->schurpre, NULL)); 2166 PetscCall(PetscOptionsScalar("-pc_fieldsplit_schur_scale", "Scale Schur complement", "PCFieldSplitSetSchurScale", jac->schurscale, &jac->schurscale, NULL)); 2167 } else if (jac->type == PC_COMPOSITE_GKB) { 2168 PetscCall(PetscOptionsReal("-pc_fieldsplit_gkb_tol", "The tolerance for the lower bound stopping criterion", "PCFieldSplitSetGKBTol", jac->gkbtol, &jac->gkbtol, NULL)); 2169 PetscCall(PetscOptionsInt("-pc_fieldsplit_gkb_delay", "The delay value for lower bound criterion", "PCFieldSplitSetGKBDelay", jac->gkbdelay, &jac->gkbdelay, NULL)); 2170 PetscCall(PetscOptionsBoundedReal("-pc_fieldsplit_gkb_nu", "Parameter in augmented Lagrangian approach", "PCFieldSplitSetGKBNu", jac->gkbnu, &jac->gkbnu, NULL, 0.0)); 2171 PetscCall(PetscOptionsInt("-pc_fieldsplit_gkb_maxit", "Maximum allowed number of iterations", "PCFieldSplitSetGKBMaxit", jac->gkbmaxit, &jac->gkbmaxit, NULL)); 2172 PetscCall(PetscOptionsBool("-pc_fieldsplit_gkb_monitor", "Prints number of GKB iterations and error", "PCFieldSplitGKB", jac->gkbmonitor, &jac->gkbmonitor, NULL)); 2173 } 2174 /* 2175 In the initial call to this routine the sub-solver data structures do not exist so we cannot call KSPSetFromOptions() on them yet. 2176 But after the initial setup of ALL the layers of sub-solvers is completed we do want to call KSPSetFromOptions() on the sub-solvers every time it 2177 is called on the outer solver in case changes were made in the options database 2178 2179 But even after PCSetUp_FieldSplit() is called all the options inside the inner levels of sub-solvers may still not have been set thus we only call the KSPSetFromOptions() 2180 if we know that the entire stack of sub-solvers below this have been complete instantiated, we check this by seeing if any solver iterations are complete. 2181 Without this extra check test p2p1fetidp_olof_full and others fail with incorrect matrix types. 2182 2183 There could be a negative side effect of calling the KSPSetFromOptions() below. 2184 2185 If one captured the PetscObjectState of the options database one could skip these calls if the database has not changed from the previous call 2186 */ 2187 if (jac->issetup) { 2188 PC_FieldSplitLink ilink = jac->head; 2189 if (jac->type == PC_COMPOSITE_SCHUR) { 2190 if (jac->kspupper && jac->kspupper->totalits > 0) PetscCall(KSPSetFromOptions(jac->kspupper)); 2191 if (jac->kspschur && jac->kspschur->totalits > 0) PetscCall(KSPSetFromOptions(jac->kspschur)); 2192 } 2193 while (ilink) { 2194 if (ilink->ksp->totalits > 0) PetscCall(KSPSetFromOptions(ilink->ksp)); 2195 ilink = ilink->next; 2196 } 2197 } 2198 PetscOptionsHeadEnd(); 2199 PetscFunctionReturn(PETSC_SUCCESS); 2200 } 2201 2202 static PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc, const char splitname[], PetscInt n, const PetscInt *fields, const PetscInt *fields_col) 2203 { 2204 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 2205 PC_FieldSplitLink ilink, next = jac->head; 2206 char prefix[128]; 2207 PetscInt i; 2208 PetscLogEvent nse; 2209 2210 PetscFunctionBegin; 2211 if (jac->splitdefined) { 2212 PetscCall(PetscInfo(pc, "Ignoring new split \"%s\" because the splits have already been defined\n", splitname)); 2213 PetscFunctionReturn(PETSC_SUCCESS); 2214 } 2215 for (i = 0; i < n; i++) PetscCheck(fields[i] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative field %" PetscInt_FMT " requested", fields[i]); 2216 PetscCall(PetscNew(&ilink)); 2217 if (splitname) { 2218 PetscCall(PetscStrallocpy(splitname, &ilink->splitname)); 2219 } else { 2220 PetscCall(PetscMalloc1(3, &ilink->splitname)); 2221 PetscCall(PetscSNPrintf(ilink->splitname, 2, "%" PetscInt_FMT, jac->nsplits)); 2222 } 2223 PetscCall(PetscMPIIntCast(jac->nsplits, &nse)); 2224 ilink->event = jac->nsplits < 5 ? KSP_Solve_FS_0 + nse : KSP_Solve_FS_0 + 4; /* Splits greater than 4 logged in 4th split */ 2225 PetscCall(PetscMalloc1(n, &ilink->fields)); 2226 PetscCall(PetscArraycpy(ilink->fields, fields, n)); 2227 PetscCall(PetscMalloc1(n, &ilink->fields_col)); 2228 PetscCall(PetscArraycpy(ilink->fields_col, fields_col, n)); 2229 2230 ilink->nfields = n; 2231 ilink->next = NULL; 2232 PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &ilink->ksp)); 2233 PetscCall(KSPSetNestLevel(ilink->ksp, pc->kspnestlevel)); 2234 PetscCall(KSPSetErrorIfNotConverged(ilink->ksp, pc->erroriffailure)); 2235 PetscCall(PetscObjectIncrementTabLevel((PetscObject)ilink->ksp, (PetscObject)pc, 1)); 2236 PetscCall(KSPSetType(ilink->ksp, KSPPREONLY)); 2237 2238 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname)); 2239 PetscCall(KSPSetOptionsPrefix(ilink->ksp, prefix)); 2240 2241 if (!next) { 2242 jac->head = ilink; 2243 ilink->previous = NULL; 2244 } else { 2245 while (next->next) next = next->next; 2246 next->next = ilink; 2247 ilink->previous = next; 2248 } 2249 jac->nsplits++; 2250 PetscFunctionReturn(PETSC_SUCCESS); 2251 } 2252 2253 static PetscErrorCode PCFieldSplitSchurGetSubKSP_FieldSplit(PC pc, PetscInt *n, KSP **subksp) 2254 { 2255 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 2256 2257 PetscFunctionBegin; 2258 *subksp = NULL; 2259 if (n) *n = 0; 2260 if (jac->type == PC_COMPOSITE_SCHUR) { 2261 PetscInt nn; 2262 2263 PetscCheck(jac->schur, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must call KSPSetUp() or PCSetUp() before calling PCFieldSplitSchurGetSubKSP()"); 2264 PetscCheck(jac->nsplits == 2, PetscObjectComm((PetscObject)pc), PETSC_ERR_PLIB, "Unexpected number of splits %" PetscInt_FMT " != 2", jac->nsplits); 2265 nn = jac->nsplits + (jac->kspupper != jac->head->ksp ? 1 : 0); 2266 PetscCall(PetscMalloc1(nn, subksp)); 2267 (*subksp)[0] = jac->head->ksp; 2268 (*subksp)[1] = jac->kspschur; 2269 if (jac->kspupper != jac->head->ksp) (*subksp)[2] = jac->kspupper; 2270 if (n) *n = nn; 2271 } 2272 PetscFunctionReturn(PETSC_SUCCESS); 2273 } 2274 2275 static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc, PetscInt *n, KSP **subksp) 2276 { 2277 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 2278 2279 PetscFunctionBegin; 2280 PetscCheck(jac->schur, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must call KSPSetUp() or PCSetUp() before calling PCFieldSplitGetSubKSP()"); 2281 PetscCall(PetscMalloc1(jac->nsplits, subksp)); 2282 PetscCall(MatSchurComplementGetKSP(jac->schur, *subksp)); 2283 2284 (*subksp)[1] = jac->kspschur; 2285 if (n) *n = jac->nsplits; 2286 PetscFunctionReturn(PETSC_SUCCESS); 2287 } 2288 2289 static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc, PetscInt *n, KSP **subksp) 2290 { 2291 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 2292 PetscInt cnt = 0; 2293 PC_FieldSplitLink ilink = jac->head; 2294 2295 PetscFunctionBegin; 2296 PetscCall(PetscMalloc1(jac->nsplits, subksp)); 2297 while (ilink) { 2298 (*subksp)[cnt++] = ilink->ksp; 2299 ilink = ilink->next; 2300 } 2301 PetscCheck(cnt == jac->nsplits, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Corrupt PCFIELDSPLIT object: number of splits in linked list %" PetscInt_FMT " does not match number in object %" PetscInt_FMT, cnt, jac->nsplits); 2302 if (n) *n = jac->nsplits; 2303 PetscFunctionReturn(PETSC_SUCCESS); 2304 } 2305 2306 /*@ 2307 PCFieldSplitRestrictIS - Restricts the fieldsplit `IS`s to be within a given `IS`. 2308 2309 Input Parameters: 2310 + pc - the preconditioner context 2311 - isy - the index set that defines the indices to which the fieldsplit is to be restricted 2312 2313 Level: advanced 2314 2315 Developer Notes: 2316 It seems the resulting `IS`s will not cover the entire space, so 2317 how can they define a convergent preconditioner? Needs explaining. 2318 2319 .seealso: [](sec_block_matrices), `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()` 2320 @*/ 2321 PetscErrorCode PCFieldSplitRestrictIS(PC pc, IS isy) 2322 { 2323 PetscFunctionBegin; 2324 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2325 PetscValidHeaderSpecific(isy, IS_CLASSID, 2); 2326 PetscTryMethod(pc, "PCFieldSplitRestrictIS_C", (PC, IS), (pc, isy)); 2327 PetscFunctionReturn(PETSC_SUCCESS); 2328 } 2329 2330 static PetscErrorCode PCFieldSplitRestrictIS_FieldSplit(PC pc, IS isy) 2331 { 2332 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 2333 PC_FieldSplitLink ilink = jac->head, next; 2334 PetscInt localsize, size, sizez, i; 2335 const PetscInt *ind, *indz; 2336 PetscInt *indc, *indcz; 2337 PetscBool flg; 2338 2339 PetscFunctionBegin; 2340 PetscCall(ISGetLocalSize(isy, &localsize)); 2341 PetscCallMPI(MPI_Scan(&localsize, &size, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)isy))); 2342 size -= localsize; 2343 while (ilink) { 2344 IS isrl, isr; 2345 PC subpc; 2346 PetscCall(ISEmbed(ilink->is, isy, PETSC_TRUE, &isrl)); 2347 PetscCall(ISGetLocalSize(isrl, &localsize)); 2348 PetscCall(PetscMalloc1(localsize, &indc)); 2349 PetscCall(ISGetIndices(isrl, &ind)); 2350 PetscCall(PetscArraycpy(indc, ind, localsize)); 2351 PetscCall(ISRestoreIndices(isrl, &ind)); 2352 PetscCall(ISDestroy(&isrl)); 2353 for (i = 0; i < localsize; i++) *(indc + i) += size; 2354 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)isy), localsize, indc, PETSC_OWN_POINTER, &isr)); 2355 PetscCall(PetscObjectReference((PetscObject)isr)); 2356 PetscCall(ISDestroy(&ilink->is)); 2357 ilink->is = isr; 2358 PetscCall(PetscObjectReference((PetscObject)isr)); 2359 PetscCall(ISDestroy(&ilink->is_col)); 2360 ilink->is_col = isr; 2361 PetscCall(ISDestroy(&isr)); 2362 PetscCall(KSPGetPC(ilink->ksp, &subpc)); 2363 PetscCall(PetscObjectTypeCompare((PetscObject)subpc, PCFIELDSPLIT, &flg)); 2364 if (flg) { 2365 IS iszl, isz; 2366 MPI_Comm comm; 2367 PetscCall(ISGetLocalSize(ilink->is, &localsize)); 2368 comm = PetscObjectComm((PetscObject)ilink->is); 2369 PetscCall(ISEmbed(isy, ilink->is, PETSC_TRUE, &iszl)); 2370 PetscCallMPI(MPI_Scan(&localsize, &sizez, 1, MPIU_INT, MPI_SUM, comm)); 2371 sizez -= localsize; 2372 PetscCall(ISGetLocalSize(iszl, &localsize)); 2373 PetscCall(PetscMalloc1(localsize, &indcz)); 2374 PetscCall(ISGetIndices(iszl, &indz)); 2375 PetscCall(PetscArraycpy(indcz, indz, localsize)); 2376 PetscCall(ISRestoreIndices(iszl, &indz)); 2377 PetscCall(ISDestroy(&iszl)); 2378 for (i = 0; i < localsize; i++) *(indcz + i) += sizez; 2379 PetscCall(ISCreateGeneral(comm, localsize, indcz, PETSC_OWN_POINTER, &isz)); 2380 PetscCall(PCFieldSplitRestrictIS(subpc, isz)); 2381 PetscCall(ISDestroy(&isz)); 2382 } 2383 next = ilink->next; 2384 ilink = next; 2385 } 2386 jac->isrestrict = PETSC_TRUE; 2387 PetscFunctionReturn(PETSC_SUCCESS); 2388 } 2389 2390 static PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc, const char splitname[], IS is) 2391 { 2392 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 2393 PC_FieldSplitLink ilink, next = jac->head; 2394 char prefix[128]; 2395 PetscLogEvent nse; 2396 2397 PetscFunctionBegin; 2398 if (jac->splitdefined) { 2399 PetscCall(PetscInfo(pc, "Ignoring new split \"%s\" because the splits have already been defined\n", splitname)); 2400 PetscFunctionReturn(PETSC_SUCCESS); 2401 } 2402 PetscCall(PetscNew(&ilink)); 2403 if (splitname) { 2404 PetscCall(PetscStrallocpy(splitname, &ilink->splitname)); 2405 } else { 2406 PetscCall(PetscMalloc1(8, &ilink->splitname)); 2407 PetscCall(PetscSNPrintf(ilink->splitname, 7, "%" PetscInt_FMT, jac->nsplits)); 2408 } 2409 PetscCall(PetscMPIIntCast(jac->nsplits, &nse)); 2410 ilink->event = jac->nsplits < 5 ? KSP_Solve_FS_0 + nse : KSP_Solve_FS_0 + 4; /* Splits greater than 4 logged in 4th split */ 2411 PetscCall(PetscObjectReference((PetscObject)is)); 2412 PetscCall(ISDestroy(&ilink->is)); 2413 ilink->is = is; 2414 PetscCall(PetscObjectReference((PetscObject)is)); 2415 PetscCall(ISDestroy(&ilink->is_col)); 2416 ilink->is_col = is; 2417 ilink->next = NULL; 2418 PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &ilink->ksp)); 2419 PetscCall(KSPSetNestLevel(ilink->ksp, pc->kspnestlevel)); 2420 PetscCall(KSPSetErrorIfNotConverged(ilink->ksp, pc->erroriffailure)); 2421 PetscCall(PetscObjectIncrementTabLevel((PetscObject)ilink->ksp, (PetscObject)pc, 1)); 2422 PetscCall(KSPSetType(ilink->ksp, KSPPREONLY)); 2423 2424 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname)); 2425 PetscCall(KSPSetOptionsPrefix(ilink->ksp, prefix)); 2426 2427 if (!next) { 2428 jac->head = ilink; 2429 ilink->previous = NULL; 2430 } else { 2431 while (next->next) next = next->next; 2432 next->next = ilink; 2433 ilink->previous = next; 2434 } 2435 jac->nsplits++; 2436 PetscFunctionReturn(PETSC_SUCCESS); 2437 } 2438 2439 /*@ 2440 PCFieldSplitSetFields - Sets the fields that define one particular split in `PCFIELDSPLIT` 2441 2442 Logically Collective 2443 2444 Input Parameters: 2445 + pc - the preconditioner context 2446 . splitname - name of this split, if `NULL` the number of the split is used 2447 . n - the number of fields in this split 2448 . fields - the fields in this split 2449 - fields_col - generally the same as `fields`, if it does not match `fields` then the submatrix that is solved for this set of fields comes from an off-diagonal block 2450 of the matrix and `fields_col` provides the column indices for that block 2451 2452 Options Database Key: 2453 . -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the `%d`'th split 2454 2455 Level: intermediate 2456 2457 Notes: 2458 Use `PCFieldSplitSetIS()` to set a general set of indices as a split. 2459 2460 If the matrix used to construct the preconditioner is `MATNEST` then field i refers to the `is_row[i]` `IS` passed to `MatCreateNest()`. 2461 2462 If the matrix used to construct the preconditioner is not `MATNEST` then 2463 `PCFieldSplitSetFields()` is for defining fields as strided blocks (based on the block size provided to the matrix with `MatSetBlockSize()` or 2464 to the `PC` with `PCFieldSplitSetBlockSize()`). For example, if the block 2465 size is three then one can define a split as 0, or 1 or 2 or 0,1 or 0,2 or 1,2 which mean 2466 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x23x56x8.. x12x45x78x.... 2467 where the numbered entries indicate what is in the split. 2468 2469 This function is called once per split (it creates a new split each time). Solve options 2470 for this split will be available under the prefix `-fieldsplit_SPLITNAME_`. 2471 2472 `PCFieldSplitSetIS()` does not support having a `fields_col` different from `fields` 2473 2474 Developer Notes: 2475 This routine does not actually create the `IS` representing the split, that is delayed 2476 until `PCSetUp_FieldSplit()`, because information about the vector/matrix layouts may not be 2477 available when this routine is called. 2478 2479 .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetBlockSize()`, `PCFieldSplitSetIS()`, `PCFieldSplitRestrictIS()`, 2480 `MatSetBlockSize()`, `MatCreateNest()` 2481 @*/ 2482 PetscErrorCode PCFieldSplitSetFields(PC pc, const char splitname[], PetscInt n, const PetscInt fields[], const PetscInt fields_col[]) 2483 { 2484 PetscFunctionBegin; 2485 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2486 PetscAssertPointer(splitname, 2); 2487 PetscCheck(n >= 1, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Provided number of fields %" PetscInt_FMT " in split \"%s\" not positive", n, splitname); 2488 PetscAssertPointer(fields, 4); 2489 PetscTryMethod(pc, "PCFieldSplitSetFields_C", (PC, const char[], PetscInt, const PetscInt *, const PetscInt *), (pc, splitname, n, fields, fields_col)); 2490 PetscFunctionReturn(PETSC_SUCCESS); 2491 } 2492 2493 /*@ 2494 PCFieldSplitSetDiagUseAmat - set flag indicating whether to extract diagonal blocks from Amat (rather than Pmat) to build 2495 the sub-matrices associated with each split. Where `KSPSetOperators`(ksp,Amat,Pmat) was used to supply the operators. 2496 2497 Logically Collective 2498 2499 Input Parameters: 2500 + pc - the preconditioner object 2501 - flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from 2502 2503 Options Database Key: 2504 . -pc_fieldsplit_diag_use_amat - use the Amat to provide the diagonal blocks 2505 2506 Level: intermediate 2507 2508 .seealso: [](sec_block_matrices), `PC`, `PCSetOperators()`, `KSPSetOperators()`, `PCFieldSplitGetDiagUseAmat()`, `PCFieldSplitSetOffDiagUseAmat()`, `PCFIELDSPLIT` 2509 @*/ 2510 PetscErrorCode PCFieldSplitSetDiagUseAmat(PC pc, PetscBool flg) 2511 { 2512 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 2513 PetscBool isfs; 2514 2515 PetscFunctionBegin; 2516 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2517 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs)); 2518 PetscCheck(isfs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "PC not of type %s", PCFIELDSPLIT); 2519 jac->diag_use_amat = flg; 2520 PetscFunctionReturn(PETSC_SUCCESS); 2521 } 2522 2523 /*@ 2524 PCFieldSplitGetDiagUseAmat - get the flag indicating whether to extract diagonal blocks from Amat (rather than Pmat) to build 2525 the sub-matrices associated with each split. Where `KSPSetOperators`(ksp,Amat,Pmat) was used to supply the operators. 2526 2527 Logically Collective 2528 2529 Input Parameter: 2530 . pc - the preconditioner object 2531 2532 Output Parameter: 2533 . flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from 2534 2535 Level: intermediate 2536 2537 .seealso: [](sec_block_matrices), `PC`, `PCSetOperators()`, `KSPSetOperators()`, `PCFieldSplitSetDiagUseAmat()`, `PCFieldSplitGetOffDiagUseAmat()`, `PCFIELDSPLIT` 2538 @*/ 2539 PetscErrorCode PCFieldSplitGetDiagUseAmat(PC pc, PetscBool *flg) 2540 { 2541 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 2542 PetscBool isfs; 2543 2544 PetscFunctionBegin; 2545 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2546 PetscAssertPointer(flg, 2); 2547 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs)); 2548 PetscCheck(isfs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "PC not of type %s", PCFIELDSPLIT); 2549 *flg = jac->diag_use_amat; 2550 PetscFunctionReturn(PETSC_SUCCESS); 2551 } 2552 2553 /*@ 2554 PCFieldSplitSetOffDiagUseAmat - set flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat) to build 2555 the sub-matrices associated with each split. Where `KSPSetOperators`(ksp,Amat,Pmat) was used to supply the operators. 2556 2557 Logically Collective 2558 2559 Input Parameters: 2560 + pc - the preconditioner object 2561 - flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from 2562 2563 Options Database Key: 2564 . -pc_fieldsplit_off_diag_use_amat <bool> - use the Amat to extract the off-diagonal blocks 2565 2566 Level: intermediate 2567 2568 .seealso: [](sec_block_matrices), `PC`, `PCSetOperators()`, `KSPSetOperators()`, `PCFieldSplitGetOffDiagUseAmat()`, `PCFieldSplitSetDiagUseAmat()`, `PCFIELDSPLIT` 2569 @*/ 2570 PetscErrorCode PCFieldSplitSetOffDiagUseAmat(PC pc, PetscBool flg) 2571 { 2572 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 2573 PetscBool isfs; 2574 2575 PetscFunctionBegin; 2576 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2577 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs)); 2578 PetscCheck(isfs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "PC not of type %s", PCFIELDSPLIT); 2579 jac->offdiag_use_amat = flg; 2580 PetscFunctionReturn(PETSC_SUCCESS); 2581 } 2582 2583 /*@ 2584 PCFieldSplitGetOffDiagUseAmat - get the flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat) to build 2585 the sub-matrices associated with each split. Where `KSPSetOperators`(ksp,Amat,Pmat) was used to supply the operators. 2586 2587 Logically Collective 2588 2589 Input Parameter: 2590 . pc - the preconditioner object 2591 2592 Output Parameter: 2593 . flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from 2594 2595 Level: intermediate 2596 2597 .seealso: [](sec_block_matrices), `PC`, `PCSetOperators()`, `KSPSetOperators()`, `PCFieldSplitSetOffDiagUseAmat()`, `PCFieldSplitGetDiagUseAmat()`, `PCFIELDSPLIT` 2598 @*/ 2599 PetscErrorCode PCFieldSplitGetOffDiagUseAmat(PC pc, PetscBool *flg) 2600 { 2601 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 2602 PetscBool isfs; 2603 2604 PetscFunctionBegin; 2605 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2606 PetscAssertPointer(flg, 2); 2607 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs)); 2608 PetscCheck(isfs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "PC not of type %s", PCFIELDSPLIT); 2609 *flg = jac->offdiag_use_amat; 2610 PetscFunctionReturn(PETSC_SUCCESS); 2611 } 2612 2613 /*@ 2614 PCFieldSplitSetIS - Sets the exact elements for a split in a `PCFIELDSPLIT` 2615 2616 Logically Collective 2617 2618 Input Parameters: 2619 + pc - the preconditioner context 2620 . splitname - name of this split, if `NULL` the number of the split is used 2621 - is - the index set that defines the elements in this split 2622 2623 Level: intermediate 2624 2625 Notes: 2626 Use `PCFieldSplitSetFields()`, for splits defined by strided `IS` based on the matrix block size or the `is_rows[]` passed into `MATNEST` 2627 2628 This function is called once per split (it creates a new split each time). Solve options 2629 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 2630 2631 .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetBlockSize()`, `PCFieldSplitSetFields()` 2632 @*/ 2633 PetscErrorCode PCFieldSplitSetIS(PC pc, const char splitname[], IS is) 2634 { 2635 PetscFunctionBegin; 2636 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2637 if (splitname) PetscAssertPointer(splitname, 2); 2638 PetscValidHeaderSpecific(is, IS_CLASSID, 3); 2639 PetscTryMethod(pc, "PCFieldSplitSetIS_C", (PC, const char[], IS), (pc, splitname, is)); 2640 PetscFunctionReturn(PETSC_SUCCESS); 2641 } 2642 2643 /*@ 2644 PCFieldSplitGetIS - Retrieves the elements for a split as an `IS` 2645 2646 Logically Collective 2647 2648 Input Parameters: 2649 + pc - the preconditioner context 2650 - splitname - name of this split 2651 2652 Output Parameter: 2653 . is - the index set that defines the elements in this split, or `NULL` if the split is not found 2654 2655 Level: intermediate 2656 2657 .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetIS()`, `PCFieldSplitGetISByIndex()` 2658 @*/ 2659 PetscErrorCode PCFieldSplitGetIS(PC pc, const char splitname[], IS *is) 2660 { 2661 PetscFunctionBegin; 2662 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2663 PetscAssertPointer(splitname, 2); 2664 PetscAssertPointer(is, 3); 2665 { 2666 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 2667 PC_FieldSplitLink ilink = jac->head; 2668 PetscBool found; 2669 2670 *is = NULL; 2671 while (ilink) { 2672 PetscCall(PetscStrcmp(ilink->splitname, splitname, &found)); 2673 if (found) { 2674 *is = ilink->is; 2675 break; 2676 } 2677 ilink = ilink->next; 2678 } 2679 } 2680 PetscFunctionReturn(PETSC_SUCCESS); 2681 } 2682 2683 /*@ 2684 PCFieldSplitGetISByIndex - Retrieves the elements for a given split as an `IS` 2685 2686 Logically Collective 2687 2688 Input Parameters: 2689 + pc - the preconditioner context 2690 - index - index of this split 2691 2692 Output Parameter: 2693 . is - the index set that defines the elements in this split 2694 2695 Level: intermediate 2696 2697 .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitGetIS()`, `PCFieldSplitSetIS()`, 2698 2699 @*/ 2700 PetscErrorCode PCFieldSplitGetISByIndex(PC pc, PetscInt index, IS *is) 2701 { 2702 PetscFunctionBegin; 2703 PetscCheck(index >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative field %" PetscInt_FMT " requested", index); 2704 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2705 PetscAssertPointer(is, 3); 2706 { 2707 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 2708 PC_FieldSplitLink ilink = jac->head; 2709 PetscInt i = 0; 2710 PetscCheck(index < jac->nsplits, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field %" PetscInt_FMT " requested but only %" PetscInt_FMT " exist", index, jac->nsplits); 2711 2712 while (i < index) { 2713 ilink = ilink->next; 2714 ++i; 2715 } 2716 PetscCall(PCFieldSplitGetIS(pc, ilink->splitname, is)); 2717 } 2718 PetscFunctionReturn(PETSC_SUCCESS); 2719 } 2720 2721 /*@ 2722 PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the 2723 fieldsplit preconditioner when calling `PCFieldSplitSetFields()`. If not set the matrix block size is used. 2724 2725 Logically Collective 2726 2727 Input Parameters: 2728 + pc - the preconditioner context 2729 - bs - the block size 2730 2731 Level: intermediate 2732 2733 Note: 2734 If the matrix is a `MATNEST` then the `is_rows[]` passed to `MatCreateNest()` determines the fields. 2735 2736 .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()` 2737 @*/ 2738 PetscErrorCode PCFieldSplitSetBlockSize(PC pc, PetscInt bs) 2739 { 2740 PetscFunctionBegin; 2741 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2742 PetscValidLogicalCollectiveInt(pc, bs, 2); 2743 PetscTryMethod(pc, "PCFieldSplitSetBlockSize_C", (PC, PetscInt), (pc, bs)); 2744 PetscFunctionReturn(PETSC_SUCCESS); 2745 } 2746 2747 /*@C 2748 PCFieldSplitGetSubKSP - Gets the `KSP` contexts for all splits 2749 2750 Collective 2751 2752 Input Parameter: 2753 . pc - the preconditioner context 2754 2755 Output Parameters: 2756 + n - the number of splits 2757 - subksp - the array of `KSP` contexts 2758 2759 Level: advanced 2760 2761 Notes: 2762 After `PCFieldSplitGetSubKSP()` the array of `KSP`s is to be freed by the user with `PetscFree()` 2763 (not the `KSP`, just the array that contains them). 2764 2765 You must call `PCSetUp()` before calling `PCFieldSplitGetSubKSP()`. 2766 2767 If the fieldsplit is of type `PC_COMPOSITE_SCHUR`, it returns the `KSP` object used inside the 2768 Schur complement and the `KSP` object used to iterate over the Schur complement. 2769 To access all the `KSP` objects used in `PC_COMPOSITE_SCHUR`, use `PCFieldSplitSchurGetSubKSP()`. 2770 2771 If the fieldsplit is of type `PC_COMPOSITE_GKB`, it returns the `KSP` object used to solve the 2772 inner linear system defined by the matrix H in each loop. 2773 2774 Fortran Note: 2775 Call `PCFieldSplitRestoreSubKSP()` when the array of `KSP` is no longer needed 2776 2777 Developer Notes: 2778 There should be a `PCFieldSplitRestoreSubKSP()` instead of requiring the user to call `PetscFree()` 2779 2780 .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()`, `PCFieldSplitSchurGetSubKSP()` 2781 @*/ 2782 PetscErrorCode PCFieldSplitGetSubKSP(PC pc, PetscInt *n, KSP *subksp[]) 2783 { 2784 PetscFunctionBegin; 2785 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2786 if (n) PetscAssertPointer(n, 2); 2787 PetscUseMethod(pc, "PCFieldSplitGetSubKSP_C", (PC, PetscInt *, KSP **), (pc, n, subksp)); 2788 PetscFunctionReturn(PETSC_SUCCESS); 2789 } 2790 2791 /*@C 2792 PCFieldSplitSchurGetSubKSP - Gets the `KSP` contexts used inside the Schur complement based `PCFIELDSPLIT` 2793 2794 Collective 2795 2796 Input Parameter: 2797 . pc - the preconditioner context 2798 2799 Output Parameters: 2800 + n - the number of splits 2801 - subksp - the array of `KSP` contexts 2802 2803 Level: advanced 2804 2805 Notes: 2806 After `PCFieldSplitSchurGetSubKSP()` the array of `KSP`s is to be freed by the user with `PetscFree()` 2807 (not the `KSP` just the array that contains them). 2808 2809 You must call `PCSetUp()` before calling `PCFieldSplitSchurGetSubKSP()`. 2810 2811 If the fieldsplit type is of type `PC_COMPOSITE_SCHUR`, it returns (in order) 2812 + 1 - the `KSP` used for the (1,1) block 2813 . 2 - the `KSP` used for the Schur complement (not the one used for the interior Schur solver) 2814 - 3 - the `KSP` used for the (1,1) block in the upper triangular factor (if different from that of the (1,1) block). 2815 2816 It returns a null array if the fieldsplit is not of type `PC_COMPOSITE_SCHUR`; in this case, you should use `PCFieldSplitGetSubKSP()`. 2817 2818 Fortran Note: 2819 Call `PCFieldSplitSchurRestoreSubKSP()` when the array of `KSP` is no longer needed 2820 2821 Developer Notes: 2822 There should be a `PCFieldSplitRestoreSubKSP()` instead of requiring the user to call `PetscFree()` 2823 2824 Should the functionality of `PCFieldSplitSchurGetSubKSP()` and `PCFieldSplitGetSubKSP()` be merged? 2825 2826 .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()`, `PCFieldSplitGetSubKSP()` 2827 @*/ 2828 PetscErrorCode PCFieldSplitSchurGetSubKSP(PC pc, PetscInt *n, KSP *subksp[]) 2829 { 2830 PetscFunctionBegin; 2831 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2832 if (n) PetscAssertPointer(n, 2); 2833 PetscUseMethod(pc, "PCFieldSplitSchurGetSubKSP_C", (PC, PetscInt *, KSP **), (pc, n, subksp)); 2834 PetscFunctionReturn(PETSC_SUCCESS); 2835 } 2836 2837 /*@ 2838 PCFieldSplitSetSchurPre - Indicates from what operator the preconditioner is constructed for the Schur complement. 2839 The default is the A11 matrix. 2840 2841 Collective 2842 2843 Input Parameters: 2844 + pc - the preconditioner context 2845 . ptype - which matrix to use for preconditioning the Schur complement: `PC_FIELDSPLIT_SCHUR_PRE_A11` (default), 2846 `PC_FIELDSPLIT_SCHUR_PRE_SELF`, `PC_FIELDSPLIT_SCHUR_PRE_USER`, 2847 `PC_FIELDSPLIT_SCHUR_PRE_SELFP`, and `PC_FIELDSPLIT_SCHUR_PRE_FULL` 2848 - pre - matrix to use for preconditioning, or `NULL` 2849 2850 Options Database Keys: 2851 + -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is `a11`. See notes for meaning of various arguments 2852 - -fieldsplit_1_pc_type <pctype> - the preconditioner algorithm that is used to construct the preconditioner from the operator 2853 2854 Level: intermediate 2855 2856 Notes: 2857 If ptype is 2858 + a11 - the preconditioner for the Schur complement is generated from the block diagonal part of the preconditioner 2859 matrix associated with the Schur complement (i.e. A11), not the Schur complement matrix 2860 . self - the preconditioner for the Schur complement is generated from the symbolic representation of the Schur complement matrix: 2861 The only preconditioners that currently work with this symbolic representation matrix object are `PCLSC` and `PCHPDDM` 2862 . user - the preconditioner for the Schur complement is generated from the user provided matrix (pre argument 2863 to this function). 2864 . selfp - the preconditioning for the Schur complement is generated from an explicitly-assembled approximation $ Sp = A11 - A10 inv(diag(A00)) A01 $ 2865 This is only a good preconditioner when diag(A00) is a good preconditioner for A00. Optionally, A00 can be 2866 lumped before extracting the diagonal using the additional option `-fieldsplit_1_mat_schur_complement_ainv_type lump` 2867 - full - the preconditioner for the Schur complement is generated from the exact Schur complement matrix representation 2868 computed internally by `PCFIELDSPLIT` (this is expensive) 2869 useful mostly as a test that the Schur complement approach can work for your problem 2870 2871 When solving a saddle point problem, where the A11 block is identically zero, using `a11` as the ptype only makes sense 2872 with the additional option `-fieldsplit_1_pc_type none`. Usually for saddle point problems one would use a `ptype` of `self` and 2873 `-fieldsplit_1_pc_type lsc` which uses the least squares commutator to compute a preconditioner for the Schur complement. 2874 2875 Developer Note: 2876 The name of this function and the option `-pc_fieldsplit_schur_precondition` are inconsistent; precondition should be used everywhere. 2877 2878 .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSchurPre()`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurPreType`, 2879 `MatSchurComplementSetAinvType()`, `PCLSC`, `PCFieldSplitSetSchurFactType()` 2880 @*/ 2881 PetscErrorCode PCFieldSplitSetSchurPre(PC pc, PCFieldSplitSchurPreType ptype, Mat pre) 2882 { 2883 PetscFunctionBegin; 2884 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2885 PetscTryMethod(pc, "PCFieldSplitSetSchurPre_C", (PC, PCFieldSplitSchurPreType, Mat), (pc, ptype, pre)); 2886 PetscFunctionReturn(PETSC_SUCCESS); 2887 } 2888 2889 PetscErrorCode PCFieldSplitSchurPrecondition(PC pc, PCFieldSplitSchurPreType ptype, Mat pre) 2890 { 2891 return PCFieldSplitSetSchurPre(pc, ptype, pre); 2892 } /* Deprecated name */ 2893 2894 /*@ 2895 PCFieldSplitGetSchurPre - For Schur complement fieldsplit, determine how the Schur complement will be 2896 preconditioned. See `PCFieldSplitSetSchurPre()` for details. 2897 2898 Logically Collective 2899 2900 Input Parameter: 2901 . pc - the preconditioner context 2902 2903 Output Parameters: 2904 + ptype - which matrix to use for preconditioning the Schur complement: `PC_FIELDSPLIT_SCHUR_PRE_A11`, `PC_FIELDSPLIT_SCHUR_PRE_SELF`, `PC_FIELDSPLIT_SCHUR_PRE_USER` 2905 - pre - matrix to use for preconditioning (with `PC_FIELDSPLIT_SCHUR_PRE_USER`), or `NULL` 2906 2907 Level: intermediate 2908 2909 .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitSetSchurPre()`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurPreType`, `PCLSC` 2910 @*/ 2911 PetscErrorCode PCFieldSplitGetSchurPre(PC pc, PCFieldSplitSchurPreType *ptype, Mat *pre) 2912 { 2913 PetscFunctionBegin; 2914 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2915 PetscUseMethod(pc, "PCFieldSplitGetSchurPre_C", (PC, PCFieldSplitSchurPreType *, Mat *), (pc, ptype, pre)); 2916 PetscFunctionReturn(PETSC_SUCCESS); 2917 } 2918 2919 /*@ 2920 PCFieldSplitSchurGetS - extract the `MATSCHURCOMPLEMENT` object used by this `PCFIELDSPLIT` in case it needs to be configured separately 2921 2922 Not Collective 2923 2924 Input Parameter: 2925 . pc - the preconditioner context 2926 2927 Output Parameter: 2928 . S - the Schur complement matrix 2929 2930 Level: advanced 2931 2932 Note: 2933 This matrix should not be destroyed using `MatDestroy()`; rather, use `PCFieldSplitSchurRestoreS()`. 2934 2935 .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSchurPreType`, `PCFieldSplitSetSchurPre()`, `MATSCHURCOMPLEMENT`, `PCFieldSplitSchurRestoreS()`, 2936 `MatCreateSchurComplement()`, `MatSchurComplementGetKSP()`, `MatSchurComplementComputeExplicitOperator()`, `MatGetSchurComplement()` 2937 @*/ 2938 PetscErrorCode PCFieldSplitSchurGetS(PC pc, Mat *S) 2939 { 2940 const char *t; 2941 PetscBool isfs; 2942 PC_FieldSplit *jac; 2943 2944 PetscFunctionBegin; 2945 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2946 PetscCall(PetscObjectGetType((PetscObject)pc, &t)); 2947 PetscCall(PetscStrcmp(t, PCFIELDSPLIT, &isfs)); 2948 PetscCheck(isfs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Expected PC of type PCFIELDSPLIT, got %s instead", t); 2949 jac = (PC_FieldSplit *)pc->data; 2950 PetscCheck(jac->type == PC_COMPOSITE_SCHUR, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Expected PCFIELDSPLIT of type SCHUR, got %d instead", jac->type); 2951 if (S) *S = jac->schur; 2952 PetscFunctionReturn(PETSC_SUCCESS); 2953 } 2954 2955 /*@ 2956 PCFieldSplitSchurRestoreS - returns the `MATSCHURCOMPLEMENT` matrix used by this `PC` 2957 2958 Not Collective 2959 2960 Input Parameters: 2961 + pc - the preconditioner context 2962 - S - the Schur complement matrix 2963 2964 Level: advanced 2965 2966 .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSchurPreType`, `PCFieldSplitSetSchurPre()`, `MatSchurComplement`, `PCFieldSplitSchurGetS()` 2967 @*/ 2968 PetscErrorCode PCFieldSplitSchurRestoreS(PC pc, Mat *S) 2969 { 2970 const char *t; 2971 PetscBool isfs; 2972 PC_FieldSplit *jac; 2973 2974 PetscFunctionBegin; 2975 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2976 PetscCall(PetscObjectGetType((PetscObject)pc, &t)); 2977 PetscCall(PetscStrcmp(t, PCFIELDSPLIT, &isfs)); 2978 PetscCheck(isfs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Expected PC of type PCFIELDSPLIT, got %s instead", t); 2979 jac = (PC_FieldSplit *)pc->data; 2980 PetscCheck(jac->type == PC_COMPOSITE_SCHUR, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Expected PCFIELDSPLIT of type SCHUR, got %d instead", jac->type); 2981 PetscCheck(S && (*S == jac->schur), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "MatSchurComplement restored is not the same as gotten"); 2982 PetscFunctionReturn(PETSC_SUCCESS); 2983 } 2984 2985 static PetscErrorCode PCFieldSplitSetSchurPre_FieldSplit(PC pc, PCFieldSplitSchurPreType ptype, Mat pre) 2986 { 2987 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 2988 2989 PetscFunctionBegin; 2990 jac->schurpre = ptype; 2991 if (ptype == PC_FIELDSPLIT_SCHUR_PRE_USER && pre) { 2992 PetscCall(MatDestroy(&jac->schur_user)); 2993 jac->schur_user = pre; 2994 PetscCall(PetscObjectReference((PetscObject)jac->schur_user)); 2995 } 2996 PetscFunctionReturn(PETSC_SUCCESS); 2997 } 2998 2999 static PetscErrorCode PCFieldSplitGetSchurPre_FieldSplit(PC pc, PCFieldSplitSchurPreType *ptype, Mat *pre) 3000 { 3001 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 3002 3003 PetscFunctionBegin; 3004 if (ptype) *ptype = jac->schurpre; 3005 if (pre) *pre = jac->schur_user; 3006 PetscFunctionReturn(PETSC_SUCCESS); 3007 } 3008 3009 /*@ 3010 PCFieldSplitSetSchurFactType - sets which blocks of the approximate block factorization to retain in the preconditioner {cite}`murphy2000note` and {cite}`ipsen2001note` 3011 3012 Collective 3013 3014 Input Parameters: 3015 + pc - the preconditioner context 3016 - ftype - which blocks of factorization to retain, `PC_FIELDSPLIT_SCHUR_FACT_FULL` is default 3017 3018 Options Database Key: 3019 . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> - default is `full` 3020 3021 Level: intermediate 3022 3023 Notes: 3024 The `full` factorization is 3025 3026 ```{math} 3027 \left(\begin{array}{cc} A & B \\ 3028 C & E \\ 3029 \end{array}\right) = 3030 \left(\begin{array}{cc} I & 0 \\ 3031 C A^{-1} & I \\ 3032 \end{array}\right) 3033 \left(\begin{array}{cc} A & 0 \\ 3034 0 & S \\ 3035 \end{array}\right) 3036 \left(\begin{array}{cc} I & A^{-1}B \\ 3037 0 & I \\ 3038 \end{array}\right) = L D U, 3039 ``` 3040 3041 where $ S = E - C A^{-1} B $. In practice, the full factorization is applied via block triangular solves with the grouping $L(DU)$. `upper` uses $DU$, `lower` uses $LD$, 3042 and `diag` is the diagonal part with the sign of $S$ flipped (because this makes the preconditioner positive definite for many formulations, 3043 thus allowing the use of `KSPMINRES)`. Sign flipping of $S$ can be turned off with `PCFieldSplitSetSchurScale()`. 3044 3045 If $A$ and $S$ are solved exactly 3046 + 1 - `full` factorization is a direct solver. 3047 . 2 - The preconditioned operator with `lower` or `upper` has all eigenvalues equal to 1 and minimal polynomial of degree 2, so `KSPGMRES` converges in 2 iterations. 3048 - 3 - With `diag`, the preconditioned operator has three distinct nonzero eigenvalues and minimal polynomial of degree at most 4, so `KSPGMRES` converges in at most 4 iterations. 3049 3050 If the iteration count is very low, consider using `KSPFGMRES` or `KSPGCR` which can use one less preconditioner 3051 application in this case. Note that the preconditioned operator may be highly non-normal, so such fast convergence may not be observed in practice. 3052 3053 For symmetric problems in which $A$ is positive definite and $S$ is negative definite, `diag` can be used with `KSPMINRES`. 3054 3055 A flexible method like `KSPFGMRES` or `KSPGCR`, [](sec_flexibleksp), must be used if the fieldsplit preconditioner is nonlinear (e.g., a few iterations of a Krylov method is used to solve with $A$ or $S$). 3056 3057 .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurPreType`, `PCFieldSplitSetSchurScale()`, 3058 [](sec_flexibleksp), `PCFieldSplitSetSchurPre()` 3059 @*/ 3060 PetscErrorCode PCFieldSplitSetSchurFactType(PC pc, PCFieldSplitSchurFactType ftype) 3061 { 3062 PetscFunctionBegin; 3063 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 3064 PetscTryMethod(pc, "PCFieldSplitSetSchurFactType_C", (PC, PCFieldSplitSchurFactType), (pc, ftype)); 3065 PetscFunctionReturn(PETSC_SUCCESS); 3066 } 3067 3068 static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc, PCFieldSplitSchurFactType ftype) 3069 { 3070 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 3071 3072 PetscFunctionBegin; 3073 jac->schurfactorization = ftype; 3074 PetscFunctionReturn(PETSC_SUCCESS); 3075 } 3076 3077 /*@ 3078 PCFieldSplitSetSchurScale - Controls the sign flip of S for `PC_FIELDSPLIT_SCHUR_FACT_DIAG`. 3079 3080 Collective 3081 3082 Input Parameters: 3083 + pc - the preconditioner context 3084 - scale - scaling factor for the Schur complement 3085 3086 Options Database Key: 3087 . -pc_fieldsplit_schur_scale <scale> - default is -1.0 3088 3089 Level: intermediate 3090 3091 .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurFactType`, `PCFieldSplitSetSchurFactType()` 3092 @*/ 3093 PetscErrorCode PCFieldSplitSetSchurScale(PC pc, PetscScalar scale) 3094 { 3095 PetscFunctionBegin; 3096 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 3097 PetscValidLogicalCollectiveScalar(pc, scale, 2); 3098 PetscTryMethod(pc, "PCFieldSplitSetSchurScale_C", (PC, PetscScalar), (pc, scale)); 3099 PetscFunctionReturn(PETSC_SUCCESS); 3100 } 3101 3102 static PetscErrorCode PCFieldSplitSetSchurScale_FieldSplit(PC pc, PetscScalar scale) 3103 { 3104 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 3105 3106 PetscFunctionBegin; 3107 jac->schurscale = scale; 3108 PetscFunctionReturn(PETSC_SUCCESS); 3109 } 3110 3111 /*@C 3112 PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement 3113 3114 Collective 3115 3116 Input Parameter: 3117 . pc - the preconditioner context 3118 3119 Output Parameters: 3120 + A00 - the (0,0) block 3121 . A01 - the (0,1) block 3122 . A10 - the (1,0) block 3123 - A11 - the (1,1) block 3124 3125 Level: advanced 3126 3127 Note: 3128 Use `NULL` for any unneeded output arguments 3129 3130 .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `MatSchurComplementGetSubMatrices()`, `MatSchurComplementSetSubMatrices()` 3131 @*/ 3132 PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc, Mat *A00, Mat *A01, Mat *A10, Mat *A11) 3133 { 3134 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 3135 3136 PetscFunctionBegin; 3137 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 3138 PetscCheck(jac->type == PC_COMPOSITE_SCHUR, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach."); 3139 if (A00) *A00 = jac->pmat[0]; 3140 if (A01) *A01 = jac->B; 3141 if (A10) *A10 = jac->C; 3142 if (A11) *A11 = jac->pmat[1]; 3143 PetscFunctionReturn(PETSC_SUCCESS); 3144 } 3145 3146 /*@ 3147 PCFieldSplitSetGKBTol - Sets the solver tolerance for the generalized Golub-Kahan bidiagonalization preconditioner {cite}`arioli2013` in `PCFIELDSPLIT` 3148 3149 Collective 3150 3151 Input Parameters: 3152 + pc - the preconditioner context 3153 - tolerance - the solver tolerance 3154 3155 Options Database Key: 3156 . -pc_fieldsplit_gkb_tol <tolerance> - default is 1e-5 3157 3158 Level: intermediate 3159 3160 Note: 3161 The generalized GKB algorithm {cite}`arioli2013` uses a lower bound estimate of the error in energy norm as stopping criterion. 3162 It stops once the lower bound estimate undershoots the required solver tolerance. Although the actual error might be bigger than 3163 this estimate, the stopping criterion is satisfactory in practical cases. 3164 3165 .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetGKBDelay()`, `PCFieldSplitSetGKBNu()`, `PCFieldSplitSetGKBMaxit()` 3166 @*/ 3167 PetscErrorCode PCFieldSplitSetGKBTol(PC pc, PetscReal tolerance) 3168 { 3169 PetscFunctionBegin; 3170 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 3171 PetscValidLogicalCollectiveReal(pc, tolerance, 2); 3172 PetscTryMethod(pc, "PCFieldSplitSetGKBTol_C", (PC, PetscReal), (pc, tolerance)); 3173 PetscFunctionReturn(PETSC_SUCCESS); 3174 } 3175 3176 static PetscErrorCode PCFieldSplitSetGKBTol_FieldSplit(PC pc, PetscReal tolerance) 3177 { 3178 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 3179 3180 PetscFunctionBegin; 3181 jac->gkbtol = tolerance; 3182 PetscFunctionReturn(PETSC_SUCCESS); 3183 } 3184 3185 /*@ 3186 PCFieldSplitSetGKBMaxit - Sets the maximum number of iterations for the generalized Golub-Kahan bidiagonalization preconditioner {cite}`arioli2013` in `PCFIELDSPLIT` 3187 3188 Collective 3189 3190 Input Parameters: 3191 + pc - the preconditioner context 3192 - maxit - the maximum number of iterations 3193 3194 Options Database Key: 3195 . -pc_fieldsplit_gkb_maxit <maxit> - default is 100 3196 3197 Level: intermediate 3198 3199 .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetGKBDelay()`, `PCFieldSplitSetGKBTol()`, `PCFieldSplitSetGKBNu()` 3200 @*/ 3201 PetscErrorCode PCFieldSplitSetGKBMaxit(PC pc, PetscInt maxit) 3202 { 3203 PetscFunctionBegin; 3204 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 3205 PetscValidLogicalCollectiveInt(pc, maxit, 2); 3206 PetscTryMethod(pc, "PCFieldSplitSetGKBMaxit_C", (PC, PetscInt), (pc, maxit)); 3207 PetscFunctionReturn(PETSC_SUCCESS); 3208 } 3209 3210 static PetscErrorCode PCFieldSplitSetGKBMaxit_FieldSplit(PC pc, PetscInt maxit) 3211 { 3212 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 3213 3214 PetscFunctionBegin; 3215 jac->gkbmaxit = maxit; 3216 PetscFunctionReturn(PETSC_SUCCESS); 3217 } 3218 3219 /*@ 3220 PCFieldSplitSetGKBDelay - Sets the delay in the lower bound error estimate in the generalized Golub-Kahan bidiagonalization {cite}`arioli2013` in `PCFIELDSPLIT` 3221 preconditioner. 3222 3223 Collective 3224 3225 Input Parameters: 3226 + pc - the preconditioner context 3227 - delay - the delay window in the lower bound estimate 3228 3229 Options Database Key: 3230 . -pc_fieldsplit_gkb_delay <delay> - default is 5 3231 3232 Level: intermediate 3233 3234 Notes: 3235 The algorithm uses a lower bound estimate of the error in energy norm as stopping criterion. The lower bound of the error $ ||u-u^k||_H $ 3236 is expressed as a truncated sum. The error at iteration k can only be measured at iteration (k + `delay`), and thus the algorithm needs 3237 at least (`delay` + 1) iterations to stop. 3238 3239 For more details on the generalized Golub-Kahan bidiagonalization method and its lower bound stopping criterion, please refer to {cite}`arioli2013` 3240 3241 .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetGKBNu()`, `PCFieldSplitSetGKBTol()`, `PCFieldSplitSetGKBMaxit()` 3242 @*/ 3243 PetscErrorCode PCFieldSplitSetGKBDelay(PC pc, PetscInt delay) 3244 { 3245 PetscFunctionBegin; 3246 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 3247 PetscValidLogicalCollectiveInt(pc, delay, 2); 3248 PetscTryMethod(pc, "PCFieldSplitSetGKBDelay_C", (PC, PetscInt), (pc, delay)); 3249 PetscFunctionReturn(PETSC_SUCCESS); 3250 } 3251 3252 static PetscErrorCode PCFieldSplitSetGKBDelay_FieldSplit(PC pc, PetscInt delay) 3253 { 3254 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 3255 3256 PetscFunctionBegin; 3257 jac->gkbdelay = delay; 3258 PetscFunctionReturn(PETSC_SUCCESS); 3259 } 3260 3261 /*@ 3262 PCFieldSplitSetGKBNu - Sets the scalar value nu >= 0 in the transformation H = A00 + nu*A01*A01' of the (1,1) block in the 3263 Golub-Kahan bidiagonalization preconditioner {cite}`arioli2013` in `PCFIELDSPLIT` 3264 3265 Collective 3266 3267 Input Parameters: 3268 + pc - the preconditioner context 3269 - nu - the shift parameter 3270 3271 Options Database Key: 3272 . -pc_fieldsplit_gkb_nu <nu> - default is 1 3273 3274 Level: intermediate 3275 3276 Notes: 3277 This shift is in general done to obtain better convergence properties for the outer loop of the algorithm. This is often achieved by choosing `nu` sufficiently large. However, 3278 if `nu` is chosen too large, the matrix H might be badly conditioned and the solution of the linear system $Hx = b$ in the inner loop becomes difficult. It is therefore 3279 necessary to find a good balance in between the convergence of the inner and outer loop. 3280 3281 For `nu` = 0, no shift is done. In this case A00 has to be positive definite. The matrix N in {cite}`arioli2013` is then chosen as identity. 3282 3283 .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetGKBDelay()`, `PCFieldSplitSetGKBTol()`, `PCFieldSplitSetGKBMaxit()` 3284 @*/ 3285 PetscErrorCode PCFieldSplitSetGKBNu(PC pc, PetscReal nu) 3286 { 3287 PetscFunctionBegin; 3288 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 3289 PetscValidLogicalCollectiveReal(pc, nu, 2); 3290 PetscTryMethod(pc, "PCFieldSplitSetGKBNu_C", (PC, PetscReal), (pc, nu)); 3291 PetscFunctionReturn(PETSC_SUCCESS); 3292 } 3293 3294 static PetscErrorCode PCFieldSplitSetGKBNu_FieldSplit(PC pc, PetscReal nu) 3295 { 3296 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 3297 3298 PetscFunctionBegin; 3299 jac->gkbnu = nu; 3300 PetscFunctionReturn(PETSC_SUCCESS); 3301 } 3302 3303 static PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc, PCCompositeType type) 3304 { 3305 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 3306 3307 PetscFunctionBegin; 3308 jac->type = type; 3309 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", NULL)); 3310 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurPre_C", NULL)); 3311 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSchurPre_C", NULL)); 3312 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurFactType_C", NULL)); 3313 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurScale_C", NULL)); 3314 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBTol_C", NULL)); 3315 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBMaxit_C", NULL)); 3316 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBNu_C", NULL)); 3317 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBDelay_C", NULL)); 3318 3319 if (type == PC_COMPOSITE_SCHUR) { 3320 pc->ops->apply = PCApply_FieldSplit_Schur; 3321 pc->ops->applytranspose = PCApplyTranspose_FieldSplit_Schur; 3322 pc->ops->matapply = PCMatApply_FieldSplit_Schur; 3323 pc->ops->view = PCView_FieldSplit_Schur; 3324 pc->ops->setuponblocks = PCSetUpOnBlocks_FieldSplit_Schur; 3325 3326 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", PCFieldSplitGetSubKSP_FieldSplit_Schur)); 3327 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurPre_C", PCFieldSplitSetSchurPre_FieldSplit)); 3328 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSchurPre_C", PCFieldSplitGetSchurPre_FieldSplit)); 3329 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurFactType_C", PCFieldSplitSetSchurFactType_FieldSplit)); 3330 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurScale_C", PCFieldSplitSetSchurScale_FieldSplit)); 3331 } else if (type == PC_COMPOSITE_GKB) { 3332 pc->ops->apply = PCApply_FieldSplit_GKB; 3333 pc->ops->applytranspose = NULL; 3334 pc->ops->matapply = NULL; 3335 pc->ops->view = PCView_FieldSplit_GKB; 3336 pc->ops->setuponblocks = PCSetUpOnBlocks_FieldSplit_GKB; 3337 3338 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", PCFieldSplitGetSubKSP_FieldSplit)); 3339 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBTol_C", PCFieldSplitSetGKBTol_FieldSplit)); 3340 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBMaxit_C", PCFieldSplitSetGKBMaxit_FieldSplit)); 3341 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBNu_C", PCFieldSplitSetGKBNu_FieldSplit)); 3342 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBDelay_C", PCFieldSplitSetGKBDelay_FieldSplit)); 3343 } else { 3344 pc->ops->apply = PCApply_FieldSplit; 3345 pc->ops->applytranspose = PCApplyTranspose_FieldSplit; 3346 pc->ops->matapply = PCMatApply_FieldSplit; 3347 pc->ops->view = PCView_FieldSplit; 3348 pc->ops->setuponblocks = PCSetUpOnBlocks_FieldSplit; 3349 3350 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", PCFieldSplitGetSubKSP_FieldSplit)); 3351 } 3352 PetscFunctionReturn(PETSC_SUCCESS); 3353 } 3354 3355 static PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc, PetscInt bs) 3356 { 3357 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 3358 3359 PetscFunctionBegin; 3360 PetscCheck(bs >= 1, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Blocksize must be positive, you gave %" PetscInt_FMT, bs); 3361 PetscCheck(jac->bs <= 0 || jac->bs == bs, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Cannot change fieldsplit blocksize from %" PetscInt_FMT " to %" PetscInt_FMT " after it has been set", jac->bs, bs); 3362 jac->bs = bs; 3363 PetscFunctionReturn(PETSC_SUCCESS); 3364 } 3365 3366 static PetscErrorCode PCSetCoordinates_FieldSplit(PC pc, PetscInt dim, PetscInt nloc, PetscReal coords[]) 3367 { 3368 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 3369 PC_FieldSplitLink ilink_current = jac->head; 3370 IS is_owned; 3371 3372 PetscFunctionBegin; 3373 jac->coordinates_set = PETSC_TRUE; // Internal flag 3374 PetscCall(MatGetOwnershipIS(pc->mat, &is_owned, NULL)); 3375 3376 while (ilink_current) { 3377 // For each IS, embed it to get local coords indces 3378 IS is_coords; 3379 PetscInt ndofs_block; 3380 const PetscInt *block_dofs_enumeration; // Numbering of the dofs relevant to the current block 3381 3382 // Setting drop to true for safety. It should make no difference. 3383 PetscCall(ISEmbed(ilink_current->is, is_owned, PETSC_TRUE, &is_coords)); 3384 PetscCall(ISGetLocalSize(is_coords, &ndofs_block)); 3385 PetscCall(ISGetIndices(is_coords, &block_dofs_enumeration)); 3386 3387 // Allocate coordinates vector and set it directly 3388 PetscCall(PetscMalloc1(ndofs_block * dim, &ilink_current->coords)); 3389 for (PetscInt dof = 0; dof < ndofs_block; ++dof) { 3390 for (PetscInt d = 0; d < dim; ++d) (ilink_current->coords)[dim * dof + d] = coords[dim * block_dofs_enumeration[dof] + d]; 3391 } 3392 ilink_current->dim = dim; 3393 ilink_current->ndofs = ndofs_block; 3394 PetscCall(ISRestoreIndices(is_coords, &block_dofs_enumeration)); 3395 PetscCall(ISDestroy(&is_coords)); 3396 ilink_current = ilink_current->next; 3397 } 3398 PetscCall(ISDestroy(&is_owned)); 3399 PetscFunctionReturn(PETSC_SUCCESS); 3400 } 3401 3402 /*@ 3403 PCFieldSplitSetType - Sets the type, `PCCompositeType`, of a `PCFIELDSPLIT` 3404 3405 Collective 3406 3407 Input Parameters: 3408 + pc - the preconditioner context 3409 - type - `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE` (default), `PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`, `PC_COMPOSITE_SCHUR`, 3410 `PC_COMPOSITE_GKB` 3411 3412 Options Database Key: 3413 . -pc_fieldsplit_type <one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type 3414 3415 Level: intermediate 3416 3417 .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCCompositeType`, `PCCompositeGetType()`, `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE`, 3418 `PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`, `PC_COMPOSITE_SCHUR`, `PCFieldSplitSetSchurFactType()` 3419 @*/ 3420 PetscErrorCode PCFieldSplitSetType(PC pc, PCCompositeType type) 3421 { 3422 PetscFunctionBegin; 3423 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 3424 PetscTryMethod(pc, "PCFieldSplitSetType_C", (PC, PCCompositeType), (pc, type)); 3425 PetscFunctionReturn(PETSC_SUCCESS); 3426 } 3427 3428 /*@ 3429 PCFieldSplitGetType - Gets the type, `PCCompositeType`, of a `PCFIELDSPLIT` 3430 3431 Not collective 3432 3433 Input Parameter: 3434 . pc - the preconditioner context 3435 3436 Output Parameter: 3437 . type - `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE` (default), `PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`, `PC_COMPOSITE_SCHUR` 3438 3439 Level: intermediate 3440 3441 .seealso: [](sec_block_matrices), `PC`, `PCCompositeSetType()`, `PCFIELDSPLIT`, `PCCompositeType`, `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE`, 3442 `PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`, `PC_COMPOSITE_SCHUR` 3443 @*/ 3444 PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type) 3445 { 3446 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 3447 3448 PetscFunctionBegin; 3449 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 3450 PetscAssertPointer(type, 2); 3451 *type = jac->type; 3452 PetscFunctionReturn(PETSC_SUCCESS); 3453 } 3454 3455 /*@ 3456 PCFieldSplitSetDMSplits - Flags whether `DMCreateFieldDecomposition()` should be used to define the splits in a `PCFIELDSPLIT`, whenever possible. 3457 3458 Logically Collective 3459 3460 Input Parameters: 3461 + pc - the preconditioner context 3462 - flg - boolean indicating whether to use field splits defined by the `DM` 3463 3464 Options Database Key: 3465 . -pc_fieldsplit_dm_splits <bool> - use the field splits defined by the `DM` 3466 3467 Level: intermediate 3468 3469 Developer Note: 3470 The name should be `PCFieldSplitSetUseDMSplits()`, similar change to options database 3471 3472 .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitGetDMSplits()`, `DMCreateFieldDecomposition()`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()` 3473 @*/ 3474 PetscErrorCode PCFieldSplitSetDMSplits(PC pc, PetscBool flg) 3475 { 3476 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 3477 PetscBool isfs; 3478 3479 PetscFunctionBegin; 3480 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 3481 PetscValidLogicalCollectiveBool(pc, flg, 2); 3482 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs)); 3483 if (isfs) jac->dm_splits = flg; 3484 PetscFunctionReturn(PETSC_SUCCESS); 3485 } 3486 3487 /*@ 3488 PCFieldSplitGetDMSplits - Returns flag indicating whether `DMCreateFieldDecomposition()` should be used to define the splits in a `PCFIELDSPLIT`, whenever possible. 3489 3490 Logically Collective 3491 3492 Input Parameter: 3493 . pc - the preconditioner context 3494 3495 Output Parameter: 3496 . flg - boolean indicating whether to use field splits defined by the `DM` 3497 3498 Level: intermediate 3499 3500 Developer Note: 3501 The name should be `PCFieldSplitGetUseDMSplits()` 3502 3503 .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetDMSplits()`, `DMCreateFieldDecomposition()`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()` 3504 @*/ 3505 PetscErrorCode PCFieldSplitGetDMSplits(PC pc, PetscBool *flg) 3506 { 3507 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 3508 PetscBool isfs; 3509 3510 PetscFunctionBegin; 3511 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 3512 PetscAssertPointer(flg, 2); 3513 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs)); 3514 if (isfs) { 3515 if (flg) *flg = jac->dm_splits; 3516 } 3517 PetscFunctionReturn(PETSC_SUCCESS); 3518 } 3519 3520 /*@ 3521 PCFieldSplitGetDetectSaddlePoint - Returns flag indicating whether `PCFIELDSPLIT` will attempt to automatically determine fields based on zero diagonal entries. 3522 3523 Logically Collective 3524 3525 Input Parameter: 3526 . pc - the preconditioner context 3527 3528 Output Parameter: 3529 . flg - boolean indicating whether to detect fields or not 3530 3531 Level: intermediate 3532 3533 .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetDetectSaddlePoint()` 3534 @*/ 3535 PetscErrorCode PCFieldSplitGetDetectSaddlePoint(PC pc, PetscBool *flg) 3536 { 3537 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 3538 3539 PetscFunctionBegin; 3540 *flg = jac->detect; 3541 PetscFunctionReturn(PETSC_SUCCESS); 3542 } 3543 3544 /*@ 3545 PCFieldSplitSetDetectSaddlePoint - Sets flag indicating whether `PCFIELDSPLIT` will attempt to automatically determine fields based on zero diagonal entries. 3546 3547 Logically Collective 3548 3549 Input Parameter: 3550 . pc - the preconditioner context 3551 3552 Output Parameter: 3553 . flg - boolean indicating whether to detect fields or not 3554 3555 Options Database Key: 3556 . -pc_fieldsplit_detect_saddle_point <bool> - detect and use the saddle point 3557 3558 Level: intermediate 3559 3560 Note: 3561 Also sets the split type to `PC_COMPOSITE_SCHUR` (see `PCFieldSplitSetType()`) and the Schur preconditioner type to `PC_FIELDSPLIT_SCHUR_PRE_SELF` (see `PCFieldSplitSetSchurPre()`). 3562 3563 .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitGetDetectSaddlePoint()`, `PCFieldSplitSetType()`, `PCFieldSplitSetSchurPre()`, `PC_FIELDSPLIT_SCHUR_PRE_SELF` 3564 @*/ 3565 PetscErrorCode PCFieldSplitSetDetectSaddlePoint(PC pc, PetscBool flg) 3566 { 3567 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 3568 3569 PetscFunctionBegin; 3570 jac->detect = flg; 3571 if (jac->detect) { 3572 PetscCall(PCFieldSplitSetType(pc, PC_COMPOSITE_SCHUR)); 3573 PetscCall(PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_SELF, NULL)); 3574 } 3575 PetscFunctionReturn(PETSC_SUCCESS); 3576 } 3577 3578 /*MC 3579 PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual 3580 collections of variables (that may overlap) called fields or splits. Each field often represents a different continuum variable 3581 represented on a grid, such as velocity, pressure, or temperature. 3582 In the literature these are sometimes called block preconditioners; but should not be confused with `PCBJACOBI`. 3583 See [the users manual section on "Solving Block Matrices"](sec_block_matrices) for more details. 3584 3585 Options Database Keys: 3586 + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the `%d`'th split 3587 . -pc_fieldsplit_default - automatically add any fields to additional splits that have not 3588 been supplied explicitly by `-pc_fieldsplit_%d_fields` 3589 . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields) 3590 when the matrix is not of `MatType` `MATNEST` 3591 . -pc_fieldsplit_type <additive,multiplicative,symmetric_multiplicative,schur,gkb> - type of relaxation or factorization splitting 3592 . -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is `a11`; see `PCFieldSplitSetSchurPre()` 3593 . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> - set factorization type when using `-pc_fieldsplit_type schur`; 3594 see `PCFieldSplitSetSchurFactType()` 3595 . -pc_fieldsplit_dm_splits <true,false> (default is true) - Whether to use `DMCreateFieldDecomposition()` for splits 3596 - -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero diagonal and uses Schur complement with no preconditioner as the solver 3597 3598 Options prefixes for inner solvers when using the Schur complement preconditioner are `-fieldsplit_0_` and `-fieldsplit_1_` . 3599 The options prefix for the inner solver when using the Golub-Kahan biadiagonalization preconditioner is `-fieldsplit_0_` 3600 For all other solvers they are `-fieldsplit_%d_` for the `%d`'th field; use `-fieldsplit_` for all fields. 3601 3602 To set options on the solvers for all blocks, prepend `-fieldsplit_` to all the `PC` 3603 options database keys. For example, `-fieldsplit_pc_type ilu` `-fieldsplit_pc_factor_levels 1`. 3604 3605 To set the options on the solvers separate for each block call `PCFieldSplitGetSubKSP()` 3606 and set the options directly on the resulting `KSP` object 3607 3608 Level: intermediate 3609 3610 Notes: 3611 Use `PCFieldSplitSetFields()` to set splits defined by "strided" entries or with a `MATNEST` and `PCFieldSplitSetIS()` 3612 to define a split by an arbitrary collection of entries. 3613 3614 If no splits are set, the default is used. If a `DM` is associated with the `PC` and it supports 3615 `DMCreateFieldDecomposition()`, then that is used for the default. Otherwise if the matrix is not `MATNEST`, the splits are defined by entries strided by bs, 3616 beginning at 0 then 1, etc to bs-1. The block size can be set with `PCFieldSplitSetBlockSize()`, 3617 if this is not called the block size defaults to the blocksize of the second matrix passed 3618 to `KSPSetOperators()`/`PCSetOperators()`. 3619 3620 For the Schur complement preconditioner if 3621 ```{math} 3622 J = \left[\begin{array}{cc} A_{00} & A_{01} \\ A_{10} & A_{11} \end{array}\right] 3623 ``` 3624 3625 the preconditioner using `full` factorization is logically 3626 ```{math} 3627 \left[\begin{array}{cc} I & -\text{ksp}(A_{00}) A_{01} \\ 0 & I \end{array}\right] \left[\begin{array}{cc} \text{ksp}(A_{00}) & 0 \\ 0 & \text{ksp}(S) \end{array}\right] \left[\begin{array}{cc} I & 0 \\ -A_{10} \text{ksp}(A_{00}) & I \end{array}\right] 3628 ``` 3629 where the action of $\text{ksp}(A_{00})$ is applied using the `KSP` solver with prefix `-fieldsplit_0_`. $S$ is the Schur complement 3630 ```{math} 3631 S = A_{11} - A_{10} \text{ksp}(A_{00}) A_{01} 3632 ``` 3633 which is usually dense and not stored explicitly. The action of $\text{ksp}(S)$ is computed using the `KSP` solver with prefix `-fieldsplit_splitname_` (where `splitname` 3634 was given in providing the SECOND split or 1 if not given). Accordingly, if using `PCFieldSplitGetSubKSP()`, the array of sub-`KSP` contexts will hold two `KSP`s: at its 3635 0th index, the `KSP` associated with `-fieldsplit_0_`, and at its 1st index, the `KSP` corresponding to `-fieldsplit_1_`. 3636 By default, $A_{11}$ is used to construct a preconditioner for $S$, use `PCFieldSplitSetSchurPre()` for all the possible ways to construct the preconditioner for $S$. 3637 3638 The factorization type is set using `-pc_fieldsplit_schur_fact_type <diag, lower, upper, full>`. `full` is shown above, 3639 `diag` gives 3640 ```{math} 3641 \left[\begin{array}{cc} \text{ksp}(A_{00}) & 0 \\ 0 & -\text{ksp}(S) \end{array}\right] 3642 ``` 3643 Note that, slightly counter intuitively, there is a negative in front of the $\text{ksp}(S)$ so that the preconditioner is positive definite. For SPD matrices $J$, the sign flip 3644 can be turned off with `PCFieldSplitSetSchurScale()` or by command line `-pc_fieldsplit_schur_scale 1.0`. The `lower` factorization is the inverse of 3645 ```{math} 3646 \left[\begin{array}{cc} A_{00} & 0 \\ A_{10} & S \end{array}\right] 3647 ``` 3648 where the inverses of $A_{00}$ and $S$ are applied using `KSP`s. The upper factorization is the inverse of 3649 ```{math} 3650 \left[\begin{array}{cc} A_{00} & A_{01} \\ 0 & S \end{array}\right] 3651 ``` 3652 where again the inverses of $A_{00}$ and $S$ are applied using `KSP`s. 3653 3654 If only one set of indices (one `IS`) is provided with `PCFieldSplitSetIS()` then the complement of that `IS` 3655 is used automatically for a second submatrix. 3656 3657 The fieldsplit preconditioner cannot currently be used with the `MATBAIJ` or `MATSBAIJ` data formats if the blocksize is larger than 1. 3658 Generally it should be used with the `MATAIJ` or `MATNEST` `MatType` 3659 3660 The forms of these preconditioners are closely related, if not identical, to forms derived as "Distributive Iterations", see, 3661 for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling {cite}`wesseling2009`. 3662 One can also use `PCFIELDSPLIT` inside a smoother resulting in "Distributive Smoothers". 3663 3664 See "A taxonomy and comparison of parallel block multi-level preconditioners for the incompressible Navier-Stokes equations" {cite}`elman2008tcp`. 3665 3666 The Constrained Pressure Preconditioner (CPR) can be implemented using `PCCOMPOSITE` with `PCGALERKIN`. CPR first solves an $R A P$ subsystem, updates the 3667 residual on all variables (`PCCompositeSetType(pc,PC_COMPOSITE_MULTIPLICATIVE)`), and then applies a simple ILU like preconditioner on all the variables. 3668 3669 The generalized Golub-Kahan bidiagonalization preconditioner (GKB) can be applied to symmetric $2 \times 2$ block matrices of the shape 3670 ```{math} 3671 \left[\begin{array}{cc} A_{00} & A_{01} \\ A_{01}' & 0 \end{array}\right] 3672 ``` 3673 with $A_{00}$ positive semi-definite. The implementation follows {cite}`arioli2013`. Therein, we choose $N := 1/\nu * I$ and the $(1,1)$-block of the matrix is modified to $H = _{A00} + \nu*A_{01}*A_{01}'$. 3674 A linear system $Hx = b$ has to be solved in each iteration of the GKB algorithm. This solver is chosen with the option prefix `-fieldsplit_0_`. 3675 3676 Some `PCFIELDSPLIT` variants are called physics-based preconditioners, since the preconditioner takes into account the underlying physics of the 3677 problem. But this nomenclature is not well-defined. 3678 3679 Developer Note: 3680 The Schur complement functionality of `PCFIELDSPLIT` should likely be factored into its own `PC` thus simplifying the implementation of the preconditioners and their 3681 user API. 3682 3683 .seealso: [](sec_block_matrices), `PC`, `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCLSC`, 3684 `PCFieldSplitGetSubKSP()`, `PCFieldSplitSchurGetSubKSP()`, `PCFieldSplitSetFields()`, 3685 `PCFieldSplitSetType()`, `PCFieldSplitSetIS()`, `PCFieldSplitSetSchurPre()`, `PCFieldSplitSetSchurFactType()`, 3686 `MatSchurComplementSetAinvType()`, `PCFieldSplitSetSchurScale()`, `PCFieldSplitSetDetectSaddlePoint()` 3687 M*/ 3688 3689 PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc) 3690 { 3691 PC_FieldSplit *jac; 3692 3693 PetscFunctionBegin; 3694 PetscCall(PetscNew(&jac)); 3695 3696 jac->bs = -1; 3697 jac->type = PC_COMPOSITE_MULTIPLICATIVE; 3698 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */ 3699 jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL; 3700 jac->schurscale = -1.0; 3701 jac->dm_splits = PETSC_TRUE; 3702 jac->gkbtol = 1e-5; 3703 jac->gkbdelay = 5; 3704 jac->gkbnu = 1; 3705 jac->gkbmaxit = 100; 3706 3707 pc->data = (void *)jac; 3708 3709 pc->ops->setup = PCSetUp_FieldSplit; 3710 pc->ops->reset = PCReset_FieldSplit; 3711 pc->ops->destroy = PCDestroy_FieldSplit; 3712 pc->ops->setfromoptions = PCSetFromOptions_FieldSplit; 3713 pc->ops->applyrichardson = NULL; 3714 3715 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSchurGetSubKSP_C", PCFieldSplitSchurGetSubKSP_FieldSplit)); 3716 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetFields_C", PCFieldSplitSetFields_FieldSplit)); 3717 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetIS_C", PCFieldSplitSetIS_FieldSplit)); 3718 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetType_C", PCFieldSplitSetType_FieldSplit)); 3719 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetBlockSize_C", PCFieldSplitSetBlockSize_FieldSplit)); 3720 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitRestrictIS_C", PCFieldSplitRestrictIS_FieldSplit)); 3721 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", PCSetCoordinates_FieldSplit)); 3722 3723 /* Initialize function pointers */ 3724 PetscCall(PCFieldSplitSetType(pc, jac->type)); 3725 PetscFunctionReturn(PETSC_SUCCESS); 3726 } 3727