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,"PCSetCoordinates_C",NULL)); 1605 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",NULL)); 1606 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",NULL)); 1607 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",NULL)); 1608 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",NULL)); 1609 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitRestrictIS_C",NULL)); 1610 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSchurGetSubKSP_C",NULL)); 1611 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",NULL)); 1612 1613 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBTol_C",NULL)); 1614 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBMaxit_C",NULL)); 1615 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBNu_C",NULL)); 1616 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBDelay_C",NULL)); 1617 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",NULL)); 1618 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",NULL)); 1619 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",NULL)); 1620 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",NULL)); 1621 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurScale_C",NULL)); 1622 PetscFunctionReturn(0); 1623 } 1624 1625 static PetscErrorCode PCSetFromOptions_FieldSplit(PetscOptionItems *PetscOptionsObject,PC pc) 1626 { 1627 PetscInt bs; 1628 PetscBool flg; 1629 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1630 PCCompositeType ctype; 1631 1632 PetscFunctionBegin; 1633 PetscOptionsHeadBegin(PetscOptionsObject,"FieldSplit options"); 1634 PetscCall(PetscOptionsBool("-pc_fieldsplit_dm_splits","Whether to use DMCreateFieldDecomposition() for splits","PCFieldSplitSetDMSplits",jac->dm_splits,&jac->dm_splits,NULL)); 1635 PetscCall(PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg)); 1636 if (flg) { 1637 PetscCall(PCFieldSplitSetBlockSize(pc,bs)); 1638 } 1639 jac->diag_use_amat = pc->useAmat; 1640 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)); 1641 jac->offdiag_use_amat = pc->useAmat; 1642 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)); 1643 PetscCall(PetscOptionsBool("-pc_fieldsplit_detect_saddle_point","Form 2-way split by detecting zero diagonal entries", "PCFieldSplitSetDetectSaddlePoint",jac->detect,&jac->detect,NULL)); 1644 PetscCall(PCFieldSplitSetDetectSaddlePoint(pc,jac->detect)); /* Sets split type and Schur PC type */ 1645 PetscCall(PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg)); 1646 if (flg) { 1647 PetscCall(PCFieldSplitSetType(pc,ctype)); 1648 } 1649 /* Only setup fields once */ 1650 if ((jac->bs > 0) && (jac->nsplits == 0)) { 1651 /* only allow user to set fields from command line if bs is already known. 1652 otherwise user can set them in PCFieldSplitSetDefaults() */ 1653 PetscCall(PCFieldSplitSetRuntimeSplits_Private(pc)); 1654 if (jac->splitdefined) PetscCall(PetscInfo(pc,"Splits defined using the options database\n")); 1655 } 1656 if (jac->type == PC_COMPOSITE_SCHUR) { 1657 PetscCall(PetscOptionsGetEnum(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_schur_factorization_type",PCFieldSplitSchurFactTypes,(PetscEnum*)&jac->schurfactorization,&flg)); 1658 if (flg) PetscCall(PetscInfo(pc,"Deprecated use of -pc_fieldsplit_schur_factorization_type\n")); 1659 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)); 1660 PetscCall(PetscOptionsEnum("-pc_fieldsplit_schur_precondition","How to build preconditioner for Schur complement","PCFieldSplitSetSchurPre",PCFieldSplitSchurPreTypes,(PetscEnum)jac->schurpre,(PetscEnum*)&jac->schurpre,NULL)); 1661 PetscCall(PetscOptionsScalar("-pc_fieldsplit_schur_scale","Scale Schur complement","PCFieldSplitSetSchurScale",jac->schurscale,&jac->schurscale,NULL)); 1662 } else if (jac->type == PC_COMPOSITE_GKB) { 1663 PetscCall(PetscOptionsReal("-pc_fieldsplit_gkb_tol","The tolerance for the lower bound stopping criterion","PCFieldSplitGKBTol",jac->gkbtol,&jac->gkbtol,NULL)); 1664 PetscCall(PetscOptionsInt("-pc_fieldsplit_gkb_delay","The delay value for lower bound criterion","PCFieldSplitGKBDelay",jac->gkbdelay,&jac->gkbdelay,NULL)); 1665 PetscCall(PetscOptionsReal("-pc_fieldsplit_gkb_nu","Parameter in augmented Lagrangian approach","PCFieldSplitGKBNu",jac->gkbnu,&jac->gkbnu,NULL)); 1666 PetscCheck(jac->gkbnu >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nu cannot be less than 0: value %g",(double)jac->gkbnu); 1667 PetscCall(PetscOptionsInt("-pc_fieldsplit_gkb_maxit","Maximum allowed number of iterations","PCFieldSplitGKBMaxit",jac->gkbmaxit,&jac->gkbmaxit,NULL)); 1668 PetscCall(PetscOptionsBool("-pc_fieldsplit_gkb_monitor","Prints number of GKB iterations and error","PCFieldSplitGKB",jac->gkbmonitor,&jac->gkbmonitor,NULL)); 1669 } 1670 PetscOptionsHeadEnd(); 1671 PetscFunctionReturn(0); 1672 } 1673 1674 /*------------------------------------------------------------------------------------*/ 1675 1676 static PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 1677 { 1678 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1679 PC_FieldSplitLink ilink,next = jac->head; 1680 char prefix[128]; 1681 PetscInt i; 1682 1683 PetscFunctionBegin; 1684 if (jac->splitdefined) { 1685 PetscCall(PetscInfo(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname)); 1686 PetscFunctionReturn(0); 1687 } 1688 for (i=0; i<n; i++) { 1689 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); 1690 PetscCheck(fields[i] >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %" PetscInt_FMT " requested",fields[i]); 1691 } 1692 PetscCall(PetscNew(&ilink)); 1693 if (splitname) { 1694 PetscCall(PetscStrallocpy(splitname,&ilink->splitname)); 1695 } else { 1696 PetscCall(PetscMalloc1(3,&ilink->splitname)); 1697 PetscCall(PetscSNPrintf(ilink->splitname,2,"%" PetscInt_FMT,jac->nsplits)); 1698 } 1699 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 */ 1700 PetscCall(PetscMalloc1(n,&ilink->fields)); 1701 PetscCall(PetscArraycpy(ilink->fields,fields,n)); 1702 PetscCall(PetscMalloc1(n,&ilink->fields_col)); 1703 PetscCall(PetscArraycpy(ilink->fields_col,fields_col,n)); 1704 1705 ilink->nfields = n; 1706 ilink->next = NULL; 1707 PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp)); 1708 PetscCall(KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure)); 1709 PetscCall(PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1)); 1710 PetscCall(KSPSetType(ilink->ksp,KSPPREONLY)); 1711 PetscCall(PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp)); 1712 1713 PetscCall(PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname)); 1714 PetscCall(KSPSetOptionsPrefix(ilink->ksp,prefix)); 1715 1716 if (!next) { 1717 jac->head = ilink; 1718 ilink->previous = NULL; 1719 } else { 1720 while (next->next) { 1721 next = next->next; 1722 } 1723 next->next = ilink; 1724 ilink->previous = next; 1725 } 1726 jac->nsplits++; 1727 PetscFunctionReturn(0); 1728 } 1729 1730 static PetscErrorCode PCFieldSplitSchurGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp) 1731 { 1732 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1733 1734 PetscFunctionBegin; 1735 *subksp = NULL; 1736 if (n) *n = 0; 1737 if (jac->type == PC_COMPOSITE_SCHUR) { 1738 PetscInt nn; 1739 1740 PetscCheck(jac->schur,PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call KSPSetUp() or PCSetUp() before calling PCFieldSplitSchurGetSubKSP()"); 1741 PetscCheck(jac->nsplits == 2,PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Unexpected number of splits %" PetscInt_FMT " != 2",jac->nsplits); 1742 nn = jac->nsplits + (jac->kspupper != jac->head->ksp ? 1 : 0); 1743 PetscCall(PetscMalloc1(nn,subksp)); 1744 (*subksp)[0] = jac->head->ksp; 1745 (*subksp)[1] = jac->kspschur; 1746 if (jac->kspupper != jac->head->ksp) (*subksp)[2] = jac->kspupper; 1747 if (n) *n = nn; 1748 } 1749 PetscFunctionReturn(0); 1750 } 1751 1752 static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp) 1753 { 1754 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1755 1756 PetscFunctionBegin; 1757 PetscCheck(jac->schur,PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call KSPSetUp() or PCSetUp() before calling PCFieldSplitGetSubKSP()"); 1758 PetscCall(PetscMalloc1(jac->nsplits,subksp)); 1759 PetscCall(MatSchurComplementGetKSP(jac->schur,*subksp)); 1760 1761 (*subksp)[1] = jac->kspschur; 1762 if (n) *n = jac->nsplits; 1763 PetscFunctionReturn(0); 1764 } 1765 1766 static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp) 1767 { 1768 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1769 PetscInt cnt = 0; 1770 PC_FieldSplitLink ilink = jac->head; 1771 1772 PetscFunctionBegin; 1773 PetscCall(PetscMalloc1(jac->nsplits,subksp)); 1774 while (ilink) { 1775 (*subksp)[cnt++] = ilink->ksp; 1776 ilink = ilink->next; 1777 } 1778 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); 1779 if (n) *n = jac->nsplits; 1780 PetscFunctionReturn(0); 1781 } 1782 1783 /*@C 1784 PCFieldSplitRestrictIS - Restricts the fieldsplit ISs to be within a given IS. 1785 1786 Input Parameters: 1787 + pc - the preconditioner context 1788 - is - the index set that defines the indices to which the fieldsplit is to be restricted 1789 1790 Level: advanced 1791 1792 @*/ 1793 PetscErrorCode PCFieldSplitRestrictIS(PC pc,IS isy) 1794 { 1795 PetscFunctionBegin; 1796 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1797 PetscValidHeaderSpecific(isy,IS_CLASSID,2); 1798 PetscTryMethod(pc,"PCFieldSplitRestrictIS_C",(PC,IS),(pc,isy)); 1799 PetscFunctionReturn(0); 1800 } 1801 1802 static PetscErrorCode PCFieldSplitRestrictIS_FieldSplit(PC pc, IS isy) 1803 { 1804 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1805 PC_FieldSplitLink ilink = jac->head, next; 1806 PetscInt localsize,size,sizez,i; 1807 const PetscInt *ind, *indz; 1808 PetscInt *indc, *indcz; 1809 PetscBool flg; 1810 1811 PetscFunctionBegin; 1812 PetscCall(ISGetLocalSize(isy,&localsize)); 1813 PetscCallMPI(MPI_Scan(&localsize,&size,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)isy))); 1814 size -= localsize; 1815 while (ilink) { 1816 IS isrl,isr; 1817 PC subpc; 1818 PetscCall(ISEmbed(ilink->is, isy, PETSC_TRUE, &isrl)); 1819 PetscCall(ISGetLocalSize(isrl,&localsize)); 1820 PetscCall(PetscMalloc1(localsize,&indc)); 1821 PetscCall(ISGetIndices(isrl,&ind)); 1822 PetscCall(PetscArraycpy(indc,ind,localsize)); 1823 PetscCall(ISRestoreIndices(isrl,&ind)); 1824 PetscCall(ISDestroy(&isrl)); 1825 for (i=0; i<localsize; i++) *(indc+i) += size; 1826 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)isy),localsize,indc,PETSC_OWN_POINTER,&isr)); 1827 PetscCall(PetscObjectReference((PetscObject)isr)); 1828 PetscCall(ISDestroy(&ilink->is)); 1829 ilink->is = isr; 1830 PetscCall(PetscObjectReference((PetscObject)isr)); 1831 PetscCall(ISDestroy(&ilink->is_col)); 1832 ilink->is_col = isr; 1833 PetscCall(ISDestroy(&isr)); 1834 PetscCall(KSPGetPC(ilink->ksp, &subpc)); 1835 PetscCall(PetscObjectTypeCompare((PetscObject)subpc,PCFIELDSPLIT,&flg)); 1836 if (flg) { 1837 IS iszl,isz; 1838 MPI_Comm comm; 1839 PetscCall(ISGetLocalSize(ilink->is,&localsize)); 1840 comm = PetscObjectComm((PetscObject)ilink->is); 1841 PetscCall(ISEmbed(isy, ilink->is, PETSC_TRUE, &iszl)); 1842 PetscCallMPI(MPI_Scan(&localsize,&sizez,1,MPIU_INT,MPI_SUM,comm)); 1843 sizez -= localsize; 1844 PetscCall(ISGetLocalSize(iszl,&localsize)); 1845 PetscCall(PetscMalloc1(localsize,&indcz)); 1846 PetscCall(ISGetIndices(iszl,&indz)); 1847 PetscCall(PetscArraycpy(indcz,indz,localsize)); 1848 PetscCall(ISRestoreIndices(iszl,&indz)); 1849 PetscCall(ISDestroy(&iszl)); 1850 for (i=0; i<localsize; i++) *(indcz+i) += sizez; 1851 PetscCall(ISCreateGeneral(comm,localsize,indcz,PETSC_OWN_POINTER,&isz)); 1852 PetscCall(PCFieldSplitRestrictIS(subpc,isz)); 1853 PetscCall(ISDestroy(&isz)); 1854 } 1855 next = ilink->next; 1856 ilink = next; 1857 } 1858 jac->isrestrict = PETSC_TRUE; 1859 PetscFunctionReturn(0); 1860 } 1861 1862 static PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is) 1863 { 1864 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1865 PC_FieldSplitLink ilink, next = jac->head; 1866 char prefix[128]; 1867 1868 PetscFunctionBegin; 1869 if (jac->splitdefined) { 1870 PetscCall(PetscInfo(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname)); 1871 PetscFunctionReturn(0); 1872 } 1873 PetscCall(PetscNew(&ilink)); 1874 if (splitname) { 1875 PetscCall(PetscStrallocpy(splitname,&ilink->splitname)); 1876 } else { 1877 PetscCall(PetscMalloc1(8,&ilink->splitname)); 1878 PetscCall(PetscSNPrintf(ilink->splitname,7,"%" PetscInt_FMT,jac->nsplits)); 1879 } 1880 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 */ 1881 PetscCall(PetscObjectReference((PetscObject)is)); 1882 PetscCall(ISDestroy(&ilink->is)); 1883 ilink->is = is; 1884 PetscCall(PetscObjectReference((PetscObject)is)); 1885 PetscCall(ISDestroy(&ilink->is_col)); 1886 ilink->is_col = is; 1887 ilink->next = NULL; 1888 PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp)); 1889 PetscCall(KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure)); 1890 PetscCall(PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1)); 1891 PetscCall(KSPSetType(ilink->ksp,KSPPREONLY)); 1892 PetscCall(PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp)); 1893 1894 PetscCall(PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname)); 1895 PetscCall(KSPSetOptionsPrefix(ilink->ksp,prefix)); 1896 1897 if (!next) { 1898 jac->head = ilink; 1899 ilink->previous = NULL; 1900 } else { 1901 while (next->next) { 1902 next = next->next; 1903 } 1904 next->next = ilink; 1905 ilink->previous = next; 1906 } 1907 jac->nsplits++; 1908 PetscFunctionReturn(0); 1909 } 1910 1911 /*@C 1912 PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner 1913 1914 Logically Collective on PC 1915 1916 Input Parameters: 1917 + pc - the preconditioner context 1918 . splitname - name of this split, if NULL the number of the split is used 1919 . n - the number of fields in this split 1920 - fields - the fields in this split 1921 1922 Level: intermediate 1923 1924 Notes: 1925 Use PCFieldSplitSetIS() to set a completely general set of indices as a field. 1926 1927 The PCFieldSplitSetFields() is for defining fields as strided blocks. For example, if the block 1928 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 1929 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x.... 1930 where the numbered entries indicate what is in the field. 1931 1932 This function is called once per split (it creates a new split each time). Solve options 1933 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1934 1935 Developer Note: This routine does not actually create the IS representing the split, that is delayed 1936 until PCSetUp_FieldSplit(), because information about the vector/matrix layouts may not be 1937 available when this routine is called. 1938 1939 .seealso: `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetBlockSize()`, `PCFieldSplitSetIS()` 1940 1941 @*/ 1942 PetscErrorCode PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 1943 { 1944 PetscFunctionBegin; 1945 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1946 PetscValidCharPointer(splitname,2); 1947 PetscCheck(n >= 1,PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Provided number of fields %" PetscInt_FMT " in split \"%s\" not positive",n,splitname); 1948 PetscValidIntPointer(fields,4); 1949 PetscTryMethod(pc,"PCFieldSplitSetFields_C",(PC,const char[],PetscInt,const PetscInt*,const PetscInt*),(pc,splitname,n,fields,fields_col)); 1950 PetscFunctionReturn(0); 1951 } 1952 1953 /*@ 1954 PCFieldSplitSetDiagUseAmat - set flag indicating whether to extract diagonal blocks from Amat (rather than Pmat) 1955 1956 Logically Collective on PC 1957 1958 Input Parameters: 1959 + pc - the preconditioner object 1960 - flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from 1961 1962 Options Database: 1963 . -pc_fieldsplit_diag_use_amat - use the Amat to provide the diagonal blocks 1964 1965 Level: intermediate 1966 1967 .seealso: `PCFieldSplitGetDiagUseAmat()`, `PCFieldSplitSetOffDiagUseAmat()`, `PCFIELDSPLIT` 1968 1969 @*/ 1970 PetscErrorCode PCFieldSplitSetDiagUseAmat(PC pc,PetscBool flg) 1971 { 1972 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1973 PetscBool isfs; 1974 1975 PetscFunctionBegin; 1976 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1977 PetscCall(PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs)); 1978 PetscCheck(isfs,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 1979 jac->diag_use_amat = flg; 1980 PetscFunctionReturn(0); 1981 } 1982 1983 /*@ 1984 PCFieldSplitGetDiagUseAmat - get the flag indicating whether to extract diagonal blocks from Amat (rather than Pmat) 1985 1986 Logically Collective on PC 1987 1988 Input Parameters: 1989 . pc - the preconditioner object 1990 1991 Output Parameters: 1992 . flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from 1993 1994 Level: intermediate 1995 1996 .seealso: `PCFieldSplitSetDiagUseAmat()`, `PCFieldSplitGetOffDiagUseAmat()`, `PCFIELDSPLIT` 1997 1998 @*/ 1999 PetscErrorCode PCFieldSplitGetDiagUseAmat(PC pc,PetscBool *flg) 2000 { 2001 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2002 PetscBool isfs; 2003 2004 PetscFunctionBegin; 2005 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2006 PetscValidBoolPointer(flg,2); 2007 PetscCall(PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs)); 2008 PetscCheck(isfs,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 2009 *flg = jac->diag_use_amat; 2010 PetscFunctionReturn(0); 2011 } 2012 2013 /*@ 2014 PCFieldSplitSetOffDiagUseAmat - set flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat) 2015 2016 Logically Collective on PC 2017 2018 Input Parameters: 2019 + pc - the preconditioner object 2020 - flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from 2021 2022 Options Database: 2023 . -pc_fieldsplit_off_diag_use_amat <bool> - use the Amat to extract the off-diagonal blocks 2024 2025 Level: intermediate 2026 2027 .seealso: `PCFieldSplitGetOffDiagUseAmat()`, `PCFieldSplitSetDiagUseAmat()`, `PCFIELDSPLIT` 2028 2029 @*/ 2030 PetscErrorCode PCFieldSplitSetOffDiagUseAmat(PC pc,PetscBool flg) 2031 { 2032 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2033 PetscBool isfs; 2034 2035 PetscFunctionBegin; 2036 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2037 PetscCall(PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs)); 2038 PetscCheck(isfs,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 2039 jac->offdiag_use_amat = flg; 2040 PetscFunctionReturn(0); 2041 } 2042 2043 /*@ 2044 PCFieldSplitGetOffDiagUseAmat - get the flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat) 2045 2046 Logically Collective on PC 2047 2048 Input Parameters: 2049 . pc - the preconditioner object 2050 2051 Output Parameters: 2052 . flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from 2053 2054 Level: intermediate 2055 2056 .seealso: `PCFieldSplitSetOffDiagUseAmat()`, `PCFieldSplitGetDiagUseAmat()`, `PCFIELDSPLIT` 2057 2058 @*/ 2059 PetscErrorCode PCFieldSplitGetOffDiagUseAmat(PC pc,PetscBool *flg) 2060 { 2061 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2062 PetscBool isfs; 2063 2064 PetscFunctionBegin; 2065 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2066 PetscValidBoolPointer(flg,2); 2067 PetscCall(PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs)); 2068 PetscCheck(isfs,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 2069 *flg = jac->offdiag_use_amat; 2070 PetscFunctionReturn(0); 2071 } 2072 2073 /*@C 2074 PCFieldSplitSetIS - Sets the exact elements for field 2075 2076 Logically Collective on PC 2077 2078 Input Parameters: 2079 + pc - the preconditioner context 2080 . splitname - name of this split, if NULL the number of the split is used 2081 - is - the index set that defines the vector elements in this field 2082 2083 Notes: 2084 Use PCFieldSplitSetFields(), for fields defined by strided types. 2085 2086 This function is called once per split (it creates a new split each time). Solve options 2087 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 2088 2089 Level: intermediate 2090 2091 .seealso: `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetBlockSize()` 2092 2093 @*/ 2094 PetscErrorCode PCFieldSplitSetIS(PC pc,const char splitname[],IS is) 2095 { 2096 PetscFunctionBegin; 2097 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2098 if (splitname) PetscValidCharPointer(splitname,2); 2099 PetscValidHeaderSpecific(is,IS_CLASSID,3); 2100 PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is)); 2101 PetscFunctionReturn(0); 2102 } 2103 2104 /*@C 2105 PCFieldSplitGetIS - Retrieves the elements for a field as an IS 2106 2107 Logically Collective on PC 2108 2109 Input Parameters: 2110 + pc - the preconditioner context 2111 - splitname - name of this split 2112 2113 Output Parameter: 2114 - is - the index set that defines the vector elements in this field, or NULL if the field is not found 2115 2116 Level: intermediate 2117 2118 .seealso: `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetIS()` 2119 2120 @*/ 2121 PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is) 2122 { 2123 PetscFunctionBegin; 2124 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2125 PetscValidCharPointer(splitname,2); 2126 PetscValidPointer(is,3); 2127 { 2128 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 2129 PC_FieldSplitLink ilink = jac->head; 2130 PetscBool found; 2131 2132 *is = NULL; 2133 while (ilink) { 2134 PetscCall(PetscStrcmp(ilink->splitname, splitname, &found)); 2135 if (found) { 2136 *is = ilink->is; 2137 break; 2138 } 2139 ilink = ilink->next; 2140 } 2141 } 2142 PetscFunctionReturn(0); 2143 } 2144 2145 /*@C 2146 PCFieldSplitGetISByIndex - Retrieves the elements for a given index field as an IS 2147 2148 Logically Collective on PC 2149 2150 Input Parameters: 2151 + pc - the preconditioner context 2152 - index - index of this split 2153 2154 Output Parameter: 2155 - is - the index set that defines the vector elements in this field 2156 2157 Level: intermediate 2158 2159 .seealso: `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitGetIS()`, `PCFieldSplitSetIS()` 2160 2161 @*/ 2162 PetscErrorCode PCFieldSplitGetISByIndex(PC pc,PetscInt index,IS *is) 2163 { 2164 PetscFunctionBegin; 2165 PetscCheck(index >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %" PetscInt_FMT " requested",index); 2166 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2167 PetscValidPointer(is,3); 2168 { 2169 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 2170 PC_FieldSplitLink ilink = jac->head; 2171 PetscInt i = 0; 2172 PetscCheck(index < jac->nsplits,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %" PetscInt_FMT " requested but only %" PetscInt_FMT " exist",index,jac->nsplits); 2173 2174 while (i < index) { 2175 ilink = ilink->next; 2176 ++i; 2177 } 2178 PetscCall(PCFieldSplitGetIS(pc,ilink->splitname,is)); 2179 } 2180 PetscFunctionReturn(0); 2181 } 2182 2183 /*@ 2184 PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the 2185 fieldsplit preconditioner. If not set the matrix block size is used. 2186 2187 Logically Collective on PC 2188 2189 Input Parameters: 2190 + pc - the preconditioner context 2191 - bs - the block size 2192 2193 Level: intermediate 2194 2195 .seealso: `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()` 2196 2197 @*/ 2198 PetscErrorCode PCFieldSplitSetBlockSize(PC pc,PetscInt bs) 2199 { 2200 PetscFunctionBegin; 2201 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2202 PetscValidLogicalCollectiveInt(pc,bs,2); 2203 PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs)); 2204 PetscFunctionReturn(0); 2205 } 2206 2207 /*@C 2208 PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits 2209 2210 Collective on KSP 2211 2212 Input Parameter: 2213 . pc - the preconditioner context 2214 2215 Output Parameters: 2216 + n - the number of splits 2217 - subksp - the array of KSP contexts 2218 2219 Note: 2220 After PCFieldSplitGetSubKSP() the array of KSPs is to be freed by the user with PetscFree() 2221 (not the KSP just the array that contains them). 2222 2223 You must call PCSetUp() before calling PCFieldSplitGetSubKSP(). 2224 2225 If the fieldsplit is of type PC_COMPOSITE_SCHUR, it returns the KSP object used inside the 2226 Schur complement and the KSP object used to iterate over the Schur complement. 2227 To access all the KSP objects used in PC_COMPOSITE_SCHUR, use PCFieldSplitSchurGetSubKSP(). 2228 2229 If the fieldsplit is of type PC_COMPOSITE_GKB, it returns the KSP object used to solve the 2230 inner linear system defined by the matrix H in each loop. 2231 2232 Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs. 2233 You can call PCFieldSplitGetSubKSP(pc,n,PETSC_NULL_KSP,ierr) to determine how large the 2234 KSP array must be. 2235 2236 Level: advanced 2237 2238 .seealso: `PCFIELDSPLIT` 2239 @*/ 2240 PetscErrorCode PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[]) 2241 { 2242 PetscFunctionBegin; 2243 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2244 if (n) PetscValidIntPointer(n,2); 2245 PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp)); 2246 PetscFunctionReturn(0); 2247 } 2248 2249 /*@C 2250 PCFieldSplitSchurGetSubKSP - Gets the KSP contexts used inside the Schur complement based PCFIELDSPLIT 2251 2252 Collective on KSP 2253 2254 Input Parameter: 2255 . pc - the preconditioner context 2256 2257 Output Parameters: 2258 + n - the number of splits 2259 - subksp - the array of KSP contexts 2260 2261 Note: 2262 After PCFieldSplitSchurGetSubKSP() the array of KSPs is to be freed by the user with PetscFree() 2263 (not the KSP just the array that contains them). 2264 2265 You must call PCSetUp() before calling PCFieldSplitSchurGetSubKSP(). 2266 2267 If the fieldsplit type is of type PC_COMPOSITE_SCHUR, it returns (in order) 2268 - the KSP used for the (1,1) block 2269 - the KSP used for the Schur complement (not the one used for the interior Schur solver) 2270 - the KSP used for the (1,1) block in the upper triangular factor (if different from that of the (1,1) block). 2271 2272 It returns a null array if the fieldsplit is not of type PC_COMPOSITE_SCHUR; in this case, you should use PCFieldSplitGetSubKSP(). 2273 2274 Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs. 2275 You can call PCFieldSplitSchurGetSubKSP(pc,n,PETSC_NULL_KSP,ierr) to determine how large the 2276 KSP array must be. 2277 2278 Level: advanced 2279 2280 .seealso: `PCFIELDSPLIT` 2281 @*/ 2282 PetscErrorCode PCFieldSplitSchurGetSubKSP(PC pc,PetscInt *n,KSP *subksp[]) 2283 { 2284 PetscFunctionBegin; 2285 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2286 if (n) PetscValidIntPointer(n,2); 2287 PetscUseMethod(pc,"PCFieldSplitSchurGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp)); 2288 PetscFunctionReturn(0); 2289 } 2290 2291 /*@ 2292 PCFieldSplitSetSchurPre - Indicates from what operator the preconditioner is constructucted for the Schur complement. 2293 The default is the A11 matrix. 2294 2295 Collective on PC 2296 2297 Input Parameters: 2298 + pc - the preconditioner context 2299 . 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 2300 PC_FIELDSPLIT_SCHUR_PRE_SELFP, and PC_FIELDSPLIT_SCHUR_PRE_FULL 2301 - userpre - matrix to use for preconditioning, or NULL 2302 2303 Options Database: 2304 + -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is a11. See notes for meaning of various arguments 2305 - -fieldsplit_1_pc_type <pctype> - the preconditioner algorithm that is used to construct the preconditioner from the operator 2306 2307 Notes: 2308 $ If ptype is 2309 $ a11 - the preconditioner for the Schur complement is generated from the block diagonal part of the preconditioner 2310 $ matrix associated with the Schur complement (i.e. A11), not the Schur complement matrix 2311 $ self - the preconditioner for the Schur complement is generated from the symbolic representation of the Schur complement matrix: 2312 $ The only preconditioner that currently works with this symbolic respresentation matrix object is the PCLSC 2313 $ preconditioner 2314 $ user - the preconditioner for the Schur complement is generated from the user provided matrix (pre argument 2315 $ to this function). 2316 $ selfp - the preconditioning for the Schur complement is generated from an explicitly-assembled approximation Sp = A11 - A10 inv(diag(A00)) A01 2317 $ This is only a good preconditioner when diag(A00) is a good preconditioner for A00. Optionally, A00 can be 2318 $ lumped before extracting the diagonal using the additional option -fieldsplit_1_mat_schur_complement_ainv_type lump 2319 $ full - the preconditioner for the Schur complement is generated from the exact Schur complement matrix representation computed internally by PCFIELDSPLIT (this is expensive) 2320 $ useful mostly as a test that the Schur complement approach can work for your problem 2321 2322 When solving a saddle point problem, where the A11 block is identically zero, using a11 as the ptype only makes sense 2323 with the additional option -fieldsplit_1_pc_type none. Usually for saddle point problems one would use a ptype of self and 2324 -fieldsplit_1_pc_type lsc which uses the least squares commutator to compute a preconditioner for the Schur complement. 2325 2326 Level: intermediate 2327 2328 .seealso: `PCFieldSplitGetSchurPre()`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurPreType`, 2329 `MatSchurComplementSetAinvType()`, `PCLSC` 2330 2331 @*/ 2332 PetscErrorCode PCFieldSplitSetSchurPre(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 2333 { 2334 PetscFunctionBegin; 2335 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2336 PetscTryMethod(pc,"PCFieldSplitSetSchurPre_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre)); 2337 PetscFunctionReturn(0); 2338 } 2339 2340 PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) {return PCFieldSplitSetSchurPre(pc,ptype,pre);} /* Deprecated name */ 2341 2342 /*@ 2343 PCFieldSplitGetSchurPre - For Schur complement fieldsplit, determine how the Schur complement will be 2344 preconditioned. See PCFieldSplitSetSchurPre() for details. 2345 2346 Logically Collective on PC 2347 2348 Input Parameter: 2349 . pc - the preconditioner context 2350 2351 Output Parameters: 2352 + ptype - which matrix to use for preconditioning the Schur complement: PC_FIELDSPLIT_SCHUR_PRE_A11, PC_FIELDSPLIT_SCHUR_PRE_SELF, PC_FIELDSPLIT_PRE_USER 2353 - userpre - matrix to use for preconditioning (with PC_FIELDSPLIT_PRE_USER), or NULL 2354 2355 Level: intermediate 2356 2357 .seealso: `PCFieldSplitSetSchurPre()`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurPreType`, `PCLSC` 2358 2359 @*/ 2360 PetscErrorCode PCFieldSplitGetSchurPre(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre) 2361 { 2362 PetscFunctionBegin; 2363 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2364 PetscUseMethod(pc,"PCFieldSplitGetSchurPre_C",(PC,PCFieldSplitSchurPreType*,Mat*),(pc,ptype,pre)); 2365 PetscFunctionReturn(0); 2366 } 2367 2368 /*@ 2369 PCFieldSplitSchurGetS - extract the MatSchurComplement object used by this PC in case it needs to be configured separately 2370 2371 Not collective 2372 2373 Input Parameter: 2374 . pc - the preconditioner context 2375 2376 Output Parameter: 2377 . S - the Schur complement matrix 2378 2379 Notes: 2380 This matrix should not be destroyed using MatDestroy(); rather, use PCFieldSplitSchurRestoreS(). 2381 2382 Level: advanced 2383 2384 .seealso: `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSchurPreType`, `PCFieldSplitSetSchurPre()`, `MatSchurComplement`, `PCFieldSplitSchurRestoreS()` 2385 2386 @*/ 2387 PetscErrorCode PCFieldSplitSchurGetS(PC pc,Mat *S) 2388 { 2389 const char* t; 2390 PetscBool isfs; 2391 PC_FieldSplit *jac; 2392 2393 PetscFunctionBegin; 2394 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2395 PetscCall(PetscObjectGetType((PetscObject)pc,&t)); 2396 PetscCall(PetscStrcmp(t,PCFIELDSPLIT,&isfs)); 2397 PetscCheck(isfs,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t); 2398 jac = (PC_FieldSplit*)pc->data; 2399 PetscCheck(jac->type == PC_COMPOSITE_SCHUR,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %d instead",jac->type); 2400 if (S) *S = jac->schur; 2401 PetscFunctionReturn(0); 2402 } 2403 2404 /*@ 2405 PCFieldSplitSchurRestoreS - restores the MatSchurComplement object used by this PC 2406 2407 Not collective 2408 2409 Input Parameters: 2410 + pc - the preconditioner context 2411 - S - the Schur complement matrix 2412 2413 Level: advanced 2414 2415 .seealso: `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSchurPreType`, `PCFieldSplitSetSchurPre()`, `MatSchurComplement`, `PCFieldSplitSchurGetS()` 2416 2417 @*/ 2418 PetscErrorCode PCFieldSplitSchurRestoreS(PC pc,Mat *S) 2419 { 2420 const char* t; 2421 PetscBool isfs; 2422 PC_FieldSplit *jac; 2423 2424 PetscFunctionBegin; 2425 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2426 PetscCall(PetscObjectGetType((PetscObject)pc,&t)); 2427 PetscCall(PetscStrcmp(t,PCFIELDSPLIT,&isfs)); 2428 PetscCheck(isfs,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t); 2429 jac = (PC_FieldSplit*)pc->data; 2430 PetscCheck(jac->type == PC_COMPOSITE_SCHUR,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %d instead",jac->type); 2431 PetscCheck(S && (*S == jac->schur),PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"MatSchurComplement restored is not the same as gotten"); 2432 PetscFunctionReturn(0); 2433 } 2434 2435 static PetscErrorCode PCFieldSplitSetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 2436 { 2437 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2438 2439 PetscFunctionBegin; 2440 jac->schurpre = ptype; 2441 if (ptype == PC_FIELDSPLIT_SCHUR_PRE_USER && pre) { 2442 PetscCall(MatDestroy(&jac->schur_user)); 2443 jac->schur_user = pre; 2444 PetscCall(PetscObjectReference((PetscObject)jac->schur_user)); 2445 } 2446 PetscFunctionReturn(0); 2447 } 2448 2449 static PetscErrorCode PCFieldSplitGetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre) 2450 { 2451 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2452 2453 PetscFunctionBegin; 2454 *ptype = jac->schurpre; 2455 *pre = jac->schur_user; 2456 PetscFunctionReturn(0); 2457 } 2458 2459 /*@ 2460 PCFieldSplitSetSchurFactType - sets which blocks of the approximate block factorization to retain in the preconditioner 2461 2462 Collective on PC 2463 2464 Input Parameters: 2465 + pc - the preconditioner context 2466 - ftype - which blocks of factorization to retain, PC_FIELDSPLIT_SCHUR_FACT_FULL is default 2467 2468 Options Database: 2469 . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> - default is full 2470 2471 Level: intermediate 2472 2473 Notes: 2474 The FULL factorization is 2475 2476 .vb 2477 (A B) = (1 0) (A 0) (1 Ainv*B) = L D U 2478 (C E) (C*Ainv 1) (0 S) (0 1) 2479 .vb 2480 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, 2481 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(). 2482 2483 If A and S are solved exactly 2484 .vb 2485 *) FULL factorization is a direct solver. 2486 *) 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. 2487 *) 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. 2488 .ve 2489 2490 If the iteration count is very low, consider using KSPFGMRES or KSPGCR which can use one less preconditioner 2491 application in this case. Note that the preconditioned operator may be highly non-normal, so such fast convergence may not be observed in practice. 2492 2493 For symmetric problems in which A is positive definite and S is negative definite, DIAG can be used with KSPMINRES. 2494 2495 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). 2496 2497 References: 2498 + * - Murphy, Golub, and Wathen, A note on preconditioning indefinite linear systems, SIAM J. Sci. Comput., 21 (2000). 2499 - * - Ipsen, A note on preconditioning nonsymmetric matrices, SIAM J. Sci. Comput., 23 (2001). 2500 2501 .seealso: `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurPreType`, `PCFieldSplitSetSchurScale()` 2502 @*/ 2503 PetscErrorCode PCFieldSplitSetSchurFactType(PC pc,PCFieldSplitSchurFactType ftype) 2504 { 2505 PetscFunctionBegin; 2506 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2507 PetscTryMethod(pc,"PCFieldSplitSetSchurFactType_C",(PC,PCFieldSplitSchurFactType),(pc,ftype)); 2508 PetscFunctionReturn(0); 2509 } 2510 2511 static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc,PCFieldSplitSchurFactType ftype) 2512 { 2513 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2514 2515 PetscFunctionBegin; 2516 jac->schurfactorization = ftype; 2517 PetscFunctionReturn(0); 2518 } 2519 2520 /*@ 2521 PCFieldSplitSetSchurScale - Controls the sign flip of S for PC_FIELDSPLIT_SCHUR_FACT_DIAG. 2522 2523 Collective on PC 2524 2525 Input Parameters: 2526 + pc - the preconditioner context 2527 - scale - scaling factor for the Schur complement 2528 2529 Options Database: 2530 . -pc_fieldsplit_schur_scale - default is -1.0 2531 2532 Level: intermediate 2533 2534 .seealso: `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurFactType`, `PCFieldSplitSetSchurScale()` 2535 @*/ 2536 PetscErrorCode PCFieldSplitSetSchurScale(PC pc,PetscScalar scale) 2537 { 2538 PetscFunctionBegin; 2539 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2540 PetscValidLogicalCollectiveScalar(pc,scale,2); 2541 PetscTryMethod(pc,"PCFieldSplitSetSchurScale_C",(PC,PetscScalar),(pc,scale)); 2542 PetscFunctionReturn(0); 2543 } 2544 2545 static PetscErrorCode PCFieldSplitSetSchurScale_FieldSplit(PC pc,PetscScalar scale) 2546 { 2547 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2548 2549 PetscFunctionBegin; 2550 jac->schurscale = scale; 2551 PetscFunctionReturn(0); 2552 } 2553 2554 /*@C 2555 PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement 2556 2557 Collective on KSP 2558 2559 Input Parameter: 2560 . pc - the preconditioner context 2561 2562 Output Parameters: 2563 + A00 - the (0,0) block 2564 . A01 - the (0,1) block 2565 . A10 - the (1,0) block 2566 - A11 - the (1,1) block 2567 2568 Level: advanced 2569 2570 .seealso: `PCFIELDSPLIT` 2571 @*/ 2572 PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11) 2573 { 2574 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 2575 2576 PetscFunctionBegin; 2577 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2578 PetscCheck(jac->type == PC_COMPOSITE_SCHUR,PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach."); 2579 if (A00) *A00 = jac->pmat[0]; 2580 if (A01) *A01 = jac->B; 2581 if (A10) *A10 = jac->C; 2582 if (A11) *A11 = jac->pmat[1]; 2583 PetscFunctionReturn(0); 2584 } 2585 2586 /*@ 2587 PCFieldSplitSetGKBTol - Sets the solver tolerance for the generalized Golub-Kahan bidiagonalization preconditioner. 2588 2589 Collective on PC 2590 2591 Notes: 2592 The generalized GKB algorithm uses a lower bound estimate of the error in energy norm as stopping criterion. 2593 It stops once the lower bound estimate undershoots the required solver tolerance. Although the actual error might be bigger than 2594 this estimate, the stopping criterion is satisfactory in practical cases [A13]. 2595 2596 [Ar13] Generalized Golub-Kahan bidiagonalization and stopping criteria, SIAM J. Matrix Anal. Appl., Vol. 34, No. 2, pp. 571-592, 2013. 2597 2598 Input Parameters: 2599 + pc - the preconditioner context 2600 - tolerance - the solver tolerance 2601 2602 Options Database: 2603 . -pc_fieldsplit_gkb_tol - default is 1e-5 2604 2605 Level: intermediate 2606 2607 .seealso: `PCFIELDSPLIT`, `PCFieldSplitSetGKBDelay()`, `PCFieldSplitSetGKBNu()`, `PCFieldSplitSetGKBMaxit()` 2608 @*/ 2609 PetscErrorCode PCFieldSplitSetGKBTol(PC pc,PetscReal tolerance) 2610 { 2611 PetscFunctionBegin; 2612 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2613 PetscValidLogicalCollectiveReal(pc,tolerance,2); 2614 PetscTryMethod(pc,"PCFieldSplitSetGKBTol_C",(PC,PetscReal),(pc,tolerance)); 2615 PetscFunctionReturn(0); 2616 } 2617 2618 static PetscErrorCode PCFieldSplitSetGKBTol_FieldSplit(PC pc,PetscReal tolerance) 2619 { 2620 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2621 2622 PetscFunctionBegin; 2623 jac->gkbtol = tolerance; 2624 PetscFunctionReturn(0); 2625 } 2626 2627 /*@ 2628 PCFieldSplitSetGKBMaxit - Sets the maximum number of iterations for the generalized Golub-Kahan bidiagonalization 2629 preconditioner. 2630 2631 Collective on PC 2632 2633 Input Parameters: 2634 + pc - the preconditioner context 2635 - maxit - the maximum number of iterations 2636 2637 Options Database: 2638 . -pc_fieldsplit_gkb_maxit - default is 100 2639 2640 Level: intermediate 2641 2642 .seealso: `PCFIELDSPLIT`, `PCFieldSplitSetGKBDelay()`, `PCFieldSplitSetGKBTol()`, `PCFieldSplitSetGKBNu()` 2643 @*/ 2644 PetscErrorCode PCFieldSplitSetGKBMaxit(PC pc,PetscInt maxit) 2645 { 2646 PetscFunctionBegin; 2647 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2648 PetscValidLogicalCollectiveInt(pc,maxit,2); 2649 PetscTryMethod(pc,"PCFieldSplitSetGKBMaxit_C",(PC,PetscInt),(pc,maxit)); 2650 PetscFunctionReturn(0); 2651 } 2652 2653 static PetscErrorCode PCFieldSplitSetGKBMaxit_FieldSplit(PC pc,PetscInt maxit) 2654 { 2655 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2656 2657 PetscFunctionBegin; 2658 jac->gkbmaxit = maxit; 2659 PetscFunctionReturn(0); 2660 } 2661 2662 /*@ 2663 PCFieldSplitSetGKBDelay - Sets the delay in the lower bound error estimate in the generalized Golub-Kahan bidiagonalization 2664 preconditioner. 2665 2666 Collective on PC 2667 2668 Notes: 2669 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 2670 is expressed as a truncated sum. The error at iteration k can only be measured at iteration (k + delay), and thus the algorithm needs 2671 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 2672 2673 [Ar13] Generalized Golub-Kahan bidiagonalization and stopping criteria, SIAM J. Matrix Anal. Appl., Vol. 34, No. 2, pp. 571-592, 2013. 2674 2675 Input Parameters: 2676 + pc - the preconditioner context 2677 - delay - the delay window in the lower bound estimate 2678 2679 Options Database: 2680 . -pc_fieldsplit_gkb_delay - default is 5 2681 2682 Level: intermediate 2683 2684 .seealso: `PCFIELDSPLIT`, `PCFieldSplitSetGKBNu()`, `PCFieldSplitSetGKBTol()`, `PCFieldSplitSetGKBMaxit()` 2685 @*/ 2686 PetscErrorCode PCFieldSplitSetGKBDelay(PC pc,PetscInt delay) 2687 { 2688 PetscFunctionBegin; 2689 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2690 PetscValidLogicalCollectiveInt(pc,delay,2); 2691 PetscTryMethod(pc,"PCFieldSplitSetGKBDelay_C",(PC,PetscInt),(pc,delay)); 2692 PetscFunctionReturn(0); 2693 } 2694 2695 static PetscErrorCode PCFieldSplitSetGKBDelay_FieldSplit(PC pc,PetscInt delay) 2696 { 2697 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2698 2699 PetscFunctionBegin; 2700 jac->gkbdelay = delay; 2701 PetscFunctionReturn(0); 2702 } 2703 2704 /*@ 2705 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. 2706 2707 Collective on PC 2708 2709 Notes: 2710 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, 2711 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 2712 necessary to find a good balance in between the convergence of the inner and outer loop. 2713 2714 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. 2715 2716 [Ar13] Generalized Golub-Kahan bidiagonalization and stopping criteria, SIAM J. Matrix Anal. Appl., Vol. 34, No. 2, pp. 571-592, 2013. 2717 2718 Input Parameters: 2719 + pc - the preconditioner context 2720 - nu - the shift parameter 2721 2722 Options Database: 2723 . -pc_fieldsplit_gkb_nu - default is 1 2724 2725 Level: intermediate 2726 2727 .seealso: `PCFIELDSPLIT`, `PCFieldSplitSetGKBDelay()`, `PCFieldSplitSetGKBTol()`, `PCFieldSplitSetGKBMaxit()` 2728 @*/ 2729 PetscErrorCode PCFieldSplitSetGKBNu(PC pc,PetscReal nu) 2730 { 2731 PetscFunctionBegin; 2732 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2733 PetscValidLogicalCollectiveReal(pc,nu,2); 2734 PetscTryMethod(pc,"PCFieldSplitSetGKBNu_C",(PC,PetscReal),(pc,nu)); 2735 PetscFunctionReturn(0); 2736 } 2737 2738 static PetscErrorCode PCFieldSplitSetGKBNu_FieldSplit(PC pc,PetscReal nu) 2739 { 2740 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2741 2742 PetscFunctionBegin; 2743 jac->gkbnu = nu; 2744 PetscFunctionReturn(0); 2745 } 2746 2747 static PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type) 2748 { 2749 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2750 2751 PetscFunctionBegin; 2752 jac->type = type; 2753 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",NULL)); 2754 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",NULL)); 2755 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",NULL)); 2756 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",NULL)); 2757 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurScale_C",NULL)); 2758 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBTol_C",NULL)); 2759 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBMaxit_C",NULL)); 2760 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBNu_C",NULL)); 2761 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBDelay_C",NULL)); 2762 2763 if (type == PC_COMPOSITE_SCHUR) { 2764 pc->ops->apply = PCApply_FieldSplit_Schur; 2765 pc->ops->view = PCView_FieldSplit_Schur; 2766 2767 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit_Schur)); 2768 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",PCFieldSplitSetSchurPre_FieldSplit)); 2769 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",PCFieldSplitGetSchurPre_FieldSplit)); 2770 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",PCFieldSplitSetSchurFactType_FieldSplit)); 2771 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurScale_C",PCFieldSplitSetSchurScale_FieldSplit)); 2772 } else if (type == PC_COMPOSITE_GKB) { 2773 pc->ops->apply = PCApply_FieldSplit_GKB; 2774 pc->ops->view = PCView_FieldSplit_GKB; 2775 2776 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit)); 2777 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBTol_C",PCFieldSplitSetGKBTol_FieldSplit)); 2778 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBMaxit_C",PCFieldSplitSetGKBMaxit_FieldSplit)); 2779 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBNu_C",PCFieldSplitSetGKBNu_FieldSplit)); 2780 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBDelay_C",PCFieldSplitSetGKBDelay_FieldSplit)); 2781 } else { 2782 pc->ops->apply = PCApply_FieldSplit; 2783 pc->ops->view = PCView_FieldSplit; 2784 2785 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit)); 2786 } 2787 PetscFunctionReturn(0); 2788 } 2789 2790 static PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs) 2791 { 2792 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2793 2794 PetscFunctionBegin; 2795 PetscCheck(bs >= 1,PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %" PetscInt_FMT,bs); 2796 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); 2797 jac->bs = bs; 2798 PetscFunctionReturn(0); 2799 } 2800 2801 static PetscErrorCode PCSetCoordinates_FieldSplit(PC pc, PetscInt dim, PetscInt nloc, PetscReal coords[]) 2802 { 2803 PC_FieldSplit * jac = (PC_FieldSplit*)pc->data; 2804 PC_FieldSplitLink ilink_current = jac->head; 2805 PetscInt ii; 2806 IS is_owned; 2807 2808 PetscFunctionBegin; 2809 jac->coordinates_set = PETSC_TRUE; // Internal flag 2810 PetscCall(MatGetOwnershipIS(pc->mat,&is_owned,PETSC_NULL)); 2811 2812 ii=0; 2813 while (ilink_current) { 2814 // For each IS, embed it to get local coords indces 2815 IS is_coords; 2816 PetscInt ndofs_block; 2817 const PetscInt *block_dofs_enumeration; // Numbering of the dofs relevant to the current block 2818 2819 // Setting drop to true for safety. It should make no difference. 2820 PetscCall(ISEmbed(ilink_current->is, is_owned, PETSC_TRUE, &is_coords)); 2821 PetscCall(ISGetLocalSize(is_coords, &ndofs_block)); 2822 PetscCall(ISGetIndices(is_coords, &block_dofs_enumeration)); 2823 2824 // Allocate coordinates vector and set it directly 2825 PetscCall(PetscMalloc1(ndofs_block * dim, &(ilink_current->coords))); 2826 for (PetscInt dof=0;dof<ndofs_block;++dof) { 2827 for (PetscInt d=0;d<dim;++d) { 2828 (ilink_current->coords)[dim*dof + d] = coords[dim * block_dofs_enumeration[dof] + d]; 2829 } 2830 } 2831 ilink_current->dim = dim; 2832 ilink_current->ndofs = ndofs_block; 2833 PetscCall(ISRestoreIndices(is_coords, &block_dofs_enumeration)); 2834 PetscCall(ISDestroy(&is_coords)); 2835 ilink_current = ilink_current->next; 2836 ++ii; 2837 } 2838 PetscCall(ISDestroy(&is_owned)); 2839 PetscFunctionReturn(0); 2840 } 2841 2842 /*@ 2843 PCFieldSplitSetType - Sets the type of fieldsplit preconditioner. 2844 2845 Collective on PC 2846 2847 Input Parameters: 2848 + pc - the preconditioner context 2849 - type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 2850 2851 Options Database Key: 2852 . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type 2853 2854 Level: Intermediate 2855 2856 .seealso: `PCCompositeSetType()` 2857 2858 @*/ 2859 PetscErrorCode PCFieldSplitSetType(PC pc,PCCompositeType type) 2860 { 2861 PetscFunctionBegin; 2862 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2863 PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type)); 2864 PetscFunctionReturn(0); 2865 } 2866 2867 /*@ 2868 PCFieldSplitGetType - Gets the type of fieldsplit preconditioner. 2869 2870 Not collective 2871 2872 Input Parameter: 2873 . pc - the preconditioner context 2874 2875 Output Parameter: 2876 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 2877 2878 Level: Intermediate 2879 2880 .seealso: `PCCompositeSetType()` 2881 @*/ 2882 PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type) 2883 { 2884 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 2885 2886 PetscFunctionBegin; 2887 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2888 PetscValidIntPointer(type,2); 2889 *type = jac->type; 2890 PetscFunctionReturn(0); 2891 } 2892 2893 /*@ 2894 PCFieldSplitSetDMSplits - Flags whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible. 2895 2896 Logically Collective 2897 2898 Input Parameters: 2899 + pc - the preconditioner context 2900 - flg - boolean indicating whether to use field splits defined by the DM 2901 2902 Options Database Key: 2903 . -pc_fieldsplit_dm_splits <bool> - use the field splits defined by the DM 2904 2905 Level: Intermediate 2906 2907 .seealso: `PCFieldSplitGetDMSplits()` 2908 2909 @*/ 2910 PetscErrorCode PCFieldSplitSetDMSplits(PC pc,PetscBool flg) 2911 { 2912 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2913 PetscBool isfs; 2914 2915 PetscFunctionBegin; 2916 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2917 PetscValidLogicalCollectiveBool(pc,flg,2); 2918 PetscCall(PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs)); 2919 if (isfs) { 2920 jac->dm_splits = flg; 2921 } 2922 PetscFunctionReturn(0); 2923 } 2924 2925 /*@ 2926 PCFieldSplitGetDMSplits - Returns flag indicating whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible. 2927 2928 Logically Collective 2929 2930 Input Parameter: 2931 . pc - the preconditioner context 2932 2933 Output Parameter: 2934 . flg - boolean indicating whether to use field splits defined by the DM 2935 2936 Level: Intermediate 2937 2938 .seealso: `PCFieldSplitSetDMSplits()` 2939 2940 @*/ 2941 PetscErrorCode PCFieldSplitGetDMSplits(PC pc,PetscBool* flg) 2942 { 2943 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2944 PetscBool isfs; 2945 2946 PetscFunctionBegin; 2947 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2948 PetscValidBoolPointer(flg,2); 2949 PetscCall(PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs)); 2950 if (isfs) { 2951 if (flg) *flg = jac->dm_splits; 2952 } 2953 PetscFunctionReturn(0); 2954 } 2955 2956 /*@ 2957 PCFieldSplitGetDetectSaddlePoint - Returns flag indicating whether PCFieldSplit will attempt to automatically determine fields based on zero diagonal entries. 2958 2959 Logically Collective 2960 2961 Input Parameter: 2962 . pc - the preconditioner context 2963 2964 Output Parameter: 2965 . flg - boolean indicating whether to detect fields or not 2966 2967 Level: Intermediate 2968 2969 .seealso: `PCFIELDSPLIT`, `PCFieldSplitSetDetectSaddlePoint()` 2970 2971 @*/ 2972 PetscErrorCode PCFieldSplitGetDetectSaddlePoint(PC pc,PetscBool *flg) 2973 { 2974 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2975 2976 PetscFunctionBegin; 2977 *flg = jac->detect; 2978 PetscFunctionReturn(0); 2979 } 2980 2981 /*@ 2982 PCFieldSplitSetDetectSaddlePoint - Sets flag indicating whether PCFieldSplit will attempt to automatically determine fields based on zero diagonal entries. 2983 2984 Logically Collective 2985 2986 Input Parameter: 2987 . pc - the preconditioner context 2988 2989 Output Parameter: 2990 . flg - boolean indicating whether to detect fields or not 2991 2992 Options Database Key: 2993 . -pc_fieldsplit_detect_saddle_point <bool> - detect and use the saddle point 2994 2995 Notes: 2996 Also sets the split type to PC_COMPOSITE_SCHUR (see PCFieldSplitSetType()) and the Schur preconditioner type to PC_FIELDSPLIT_SCHUR_PRE_SELF (see PCFieldSplitSetSchurPre()). 2997 2998 Level: Intermediate 2999 3000 .seealso: `PCFIELDSPLIT`, `PCFieldSplitSetDetectSaddlePoint()`, `PCFieldSplitSetType()`, `PCFieldSplitSetSchurPre()` 3001 3002 @*/ 3003 PetscErrorCode PCFieldSplitSetDetectSaddlePoint(PC pc,PetscBool flg) 3004 { 3005 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 3006 3007 PetscFunctionBegin; 3008 jac->detect = flg; 3009 if (jac->detect) { 3010 PetscCall(PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR)); 3011 PetscCall(PCFieldSplitSetSchurPre(pc,PC_FIELDSPLIT_SCHUR_PRE_SELF,NULL)); 3012 } 3013 PetscFunctionReturn(0); 3014 } 3015 3016 /* -------------------------------------------------------------------------------------*/ 3017 /*MC 3018 PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual 3019 fields or groups of fields. See [the users manual section on "Solving Block Matrices"](sec_block_matrices) for more details. 3020 3021 To set options on the solvers for each block append `-fieldsplit_` to all the PC 3022 options database keys. For example, `-fieldsplit_pc_type ilu -fieldsplit_pc_factor_levels 1` 3023 3024 To set the options on the solvers separate for each block call `PCFieldSplitGetSubKSP()` 3025 and set the options directly on the resulting `KSP` object 3026 3027 Level: intermediate 3028 3029 Options Database Keys: 3030 + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the `%d`'th split 3031 . -pc_fieldsplit_default - automatically add any fields to additional splits that have not 3032 been supplied explicitly by `-pc_fieldsplit_%d_fields` 3033 . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields) 3034 . -pc_fieldsplit_type <additive,multiplicative,symmetric_multiplicative,schur,gkb> - type of relaxation or factorization splitting 3035 . -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is a11; see `PCFieldSplitSetSchurPre()` 3036 . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> - set factorization type when using `-pc_fieldsplit_type schur`; see `PCFieldSplitSetSchurFactType()` 3037 - -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero diagonal and uses Schur complement with no preconditioner as the solver 3038 3039 Options prefixes for inner solvers when using the Schur complement preconditioner are `-fieldsplit_0_` and `-fieldsplit_1_` . 3040 For all other solvers they are `-fieldsplit_%d_` for the `d`th field; use `-fieldsplit_` for all fields. 3041 The options prefix for the inner solver when using the Golub-Kahan biadiagonalization preconditioner is `-fieldsplit_0_` 3042 3043 Notes: 3044 Use `PCFieldSplitSetFields()` to set fields defined by "strided" entries and `PCFieldSplitSetIS()` 3045 to define a field by an arbitrary collection of entries. 3046 3047 If no fields are set the default is used. The fields are defined by entries strided by bs, 3048 beginning at 0 then 1, etc to bs-1. The block size can be set with `PCFieldSplitSetBlockSize()`, 3049 if this is not called the block size defaults to the blocksize of the second matrix passed 3050 to `KSPSetOperators()`/`PCSetOperators()`. 3051 3052 For the Schur complement preconditioner if 3053 3054 ```{math} 3055 J = \left[\begin{array}{cc} A_{00} & A_{01} \\ A_{10} & A_{11} \end{array}\right] 3056 ``` 3057 3058 the preconditioner using `full` factorization is logically 3059 ```{math} 3060 \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] 3061 ``` 3062 where the action of $\text{inv}(A_{00})$ is applied using the KSP solver with prefix `-fieldsplit_0_`. $S$ is the Schur complement 3063 ```{math} 3064 S = A_{11} - A_{10} \text{ksp}(A_{00}) A_{01} 3065 ``` 3066 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 3067 in providing the SECOND split or 1 if not given). For `PCFieldSplitGetSub\text{ksp}()` when field number is 0, 3068 it returns the KSP associated with `-fieldsplit_0_` while field number 1 gives `-fieldsplit_1_` KSP. By default 3069 $A_{11}$ is used to construct a preconditioner for $S$, use `PCFieldSplitSetSchurPre()` for all the possible ways to construct the preconditioner for $S$. 3070 3071 The factorization type is set using `-pc_fieldsplit_schur_fact_type <diag, lower, upper, full>`. `full` is shown above, 3072 `diag` gives 3073 ```{math} 3074 \left[\begin{array}{cc} \text{inv}(A_{00}) & 0 \\ 0 & -\text{ksp}(S) \end{array}\right] 3075 ``` 3076 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 3077 can be turned off with `PCFieldSplitSetSchurScale()` or by command line `-pc_fieldsplit_schur_scale 1.0`. The `lower` factorization is the inverse of 3078 ```{math} 3079 \left[\begin{array}{cc} A_{00} & 0 \\ A_{10} & S \end{array}\right] 3080 ``` 3081 where the inverses of A_{00} and S are applied using KSPs. The upper factorization is the inverse of 3082 ```{math} 3083 \left[\begin{array}{cc} A_{00} & A_{01} \\ 0 & S \end{array}\right] 3084 ``` 3085 where again the inverses of $A_{00}$ and $S$ are applied using `KSP`s. 3086 3087 If only one set of indices (one `IS`) is provided with `PCFieldSplitSetIS()` then the complement of that `IS` 3088 is used automatically for a second block. 3089 3090 The fieldsplit preconditioner cannot currently be used with the BAIJ or SBAIJ data formats if the blocksize is larger than 1. 3091 Generally it should be used with the AIJ format. 3092 3093 The forms of these preconditioners are closely related if not identical to forms derived as "Distributive Iterations", see, 3094 for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling {cite}`Wesseling2009`. Note that one can also use `PCFIELDSPLIT` 3095 inside a smoother resulting in "Distributive Smoothers". 3096 3097 References: 3098 3099 See "A taxonomy and comparison of parallel block multi-level preconditioners for the incompressible Navier-Stokes equations" {cite}`elman2008tcp`. 3100 3101 The Constrained Pressure Preconditioner (CPR) can be implemented using `PCCOMPOSITE` with `PCGALERKIN`. CPR first solves an $R A P$ subsystem, updates the 3102 residual on all variables (`PCCompositeSetType(pc,PC_COMPOSITE_MULTIPLICATIVE)`), and then applies a simple ILU like preconditioner on all the variables. 3103 3104 The generalized Golub-Kahan bidiagonalization preconditioner (GKB) can be applied to symmetric $2 \times 2$ block matrices of the shape 3105 ```{math} 3106 \left[\begin{array}{cc} A_{00} & A_{01} \\ A_{01}' & 0 \end{array}\right] 3107 ``` 3108 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}'$. 3109 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_`. 3110 3111 ```{bibliography} 3112 :filter: docname in docnames 3113 ``` 3114 3115 .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCLSC`, 3116 `PCFieldSplitGetSubKSP()`, `PCFieldSplitSchurGetSubKSP()`, `PCFieldSplitSetFields()`, 3117 `PCFieldSplitSetType()`, `PCFieldSplitSetIS()`, `PCFieldSplitSetSchurPre()`, `PCFieldSplitSetSchurFactType()`, 3118 `MatSchurComplementSetAinvType()`, `PCFieldSplitSetSchurScale()`, `PCFieldSplitSetDetectSaddlePoint()` 3119 M*/ 3120 3121 PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc) 3122 { 3123 PC_FieldSplit *jac; 3124 3125 PetscFunctionBegin; 3126 PetscCall(PetscNewLog(pc,&jac)); 3127 3128 jac->bs = -1; 3129 jac->nsplits = 0; 3130 jac->type = PC_COMPOSITE_MULTIPLICATIVE; 3131 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */ 3132 jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL; 3133 jac->schurscale = -1.0; 3134 jac->dm_splits = PETSC_TRUE; 3135 jac->detect = PETSC_FALSE; 3136 jac->gkbtol = 1e-5; 3137 jac->gkbdelay = 5; 3138 jac->gkbnu = 1; 3139 jac->gkbmaxit = 100; 3140 jac->gkbmonitor = PETSC_FALSE; 3141 jac->coordinates_set = PETSC_FALSE; 3142 3143 pc->data = (void*)jac; 3144 3145 pc->ops->apply = PCApply_FieldSplit; 3146 pc->ops->applytranspose = PCApplyTranspose_FieldSplit; 3147 pc->ops->setup = PCSetUp_FieldSplit; 3148 pc->ops->reset = PCReset_FieldSplit; 3149 pc->ops->destroy = PCDestroy_FieldSplit; 3150 pc->ops->setfromoptions = PCSetFromOptions_FieldSplit; 3151 pc->ops->view = PCView_FieldSplit; 3152 pc->ops->applyrichardson = NULL; 3153 3154 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSchurGetSubKSP_C",PCFieldSplitSchurGetSubKSP_FieldSplit)); 3155 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit)); 3156 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",PCFieldSplitSetFields_FieldSplit)); 3157 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",PCFieldSplitSetIS_FieldSplit)); 3158 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",PCFieldSplitSetType_FieldSplit)); 3159 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",PCFieldSplitSetBlockSize_FieldSplit)); 3160 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitRestrictIS_C",PCFieldSplitRestrictIS_FieldSplit)); 3161 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",PCSetCoordinates_FieldSplit)); 3162 PetscFunctionReturn(0); 3163 } 3164