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