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 isset,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(MatIsSPDKnown(pc->pmat,&isset,&isspd)); 826 jac->schurscale = (isset && isspd) ? 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) PetscCall(MatSetNullSpace(jac->schur, sp)); 909 910 PetscCall(PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_inner_", ilink->splitname)); 911 PetscCall(PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options,((PetscObject)pc)->prefix, schurtestoption, NULL, &flg)); 912 if (flg) { 913 DM dmInner; 914 KSP kspInner; 915 PC pcInner; 916 917 PetscCall(MatSchurComplementGetKSP(jac->schur, &kspInner)); 918 PetscCall(KSPReset(kspInner)); 919 PetscCall(KSPSetOperators(kspInner,jac->mat[0],jac->pmat[0])); 920 PetscCall(PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_inner_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname)); 921 /* Indent this deeper to emphasize the "inner" nature of this solver. */ 922 PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspInner, (PetscObject) pc, 2)); 923 PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspInner->pc, (PetscObject) pc, 2)); 924 PetscCall(KSPSetOptionsPrefix(kspInner, schurprefix)); 925 926 /* Set DM for new solver */ 927 PetscCall(KSPGetDM(jac->head->ksp, &dmInner)); 928 PetscCall(KSPSetDM(kspInner, dmInner)); 929 PetscCall(KSPSetDMActive(kspInner, PETSC_FALSE)); 930 931 /* Defaults to PCKSP as preconditioner */ 932 PetscCall(KSPGetPC(kspInner, &pcInner)); 933 PetscCall(PCSetType(pcInner, PCKSP)); 934 PetscCall(PCKSPSetKSP(pcInner, jac->head->ksp)); 935 } else { 936 /* Use the outer solver for the inner solve, but revert the KSPPREONLY from PCFieldSplitSetFields_FieldSplit or 937 * PCFieldSplitSetIS_FieldSplit. We don't want KSPPREONLY because it makes the Schur complement inexact, 938 * preventing Schur complement reduction to be an accurate solve. Usually when an iterative solver is used for 939 * S = D - C A_inner^{-1} B, we expect S to be defined using an accurate definition of A_inner^{-1}, so we make 940 * GMRES the default. Note that it is also common to use PREONLY for S, in which case S may not be used 941 * directly, and the user is responsible for setting an inexact method for fieldsplit's A^{-1}. */ 942 PetscCall(KSPSetType(jac->head->ksp,KSPGMRES)); 943 PetscCall(MatSchurComplementSetKSP(jac->schur,jac->head->ksp)); 944 } 945 PetscCall(KSPSetOperators(jac->head->ksp,jac->mat[0],jac->pmat[0])); 946 PetscCall(KSPSetFromOptions(jac->head->ksp)); 947 PetscCall(MatSetFromOptions(jac->schur)); 948 949 PetscCall(PetscObjectTypeCompare((PetscObject)jac->schur, MATSCHURCOMPLEMENT, &flg)); 950 if (flg) { /* Need to do this otherwise PCSetUp_KSP will overwrite the amat of jac->head->ksp */ 951 KSP kspInner; 952 PC pcInner; 953 954 PetscCall(MatSchurComplementGetKSP(jac->schur, &kspInner)); 955 PetscCall(KSPGetPC(kspInner, &pcInner)); 956 PetscCall(PetscObjectTypeCompare((PetscObject)pcInner, PCKSP, &flg)); 957 if (flg) { 958 KSP ksp; 959 960 PetscCall(PCKSPGetKSP(pcInner, &ksp)); 961 if (ksp == jac->head->ksp) { 962 PetscCall(PCSetUseAmat(pcInner, PETSC_TRUE)); 963 } 964 } 965 } 966 PetscCall(PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_upper_", ilink->splitname)); 967 PetscCall(PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options,((PetscObject)pc)->prefix, schurtestoption, NULL, &flg)); 968 if (flg) { 969 DM dmInner; 970 971 PetscCall(PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_upper_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname)); 972 PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspupper)); 973 PetscCall(KSPSetErrorIfNotConverged(jac->kspupper,pc->erroriffailure)); 974 PetscCall(KSPSetOptionsPrefix(jac->kspupper, schurprefix)); 975 PetscCall(PetscObjectIncrementTabLevel((PetscObject)jac->kspupper, (PetscObject) pc, 1)); 976 PetscCall(PetscObjectIncrementTabLevel((PetscObject)jac->kspupper->pc, (PetscObject) pc, 1)); 977 PetscCall(KSPGetDM(jac->head->ksp, &dmInner)); 978 PetscCall(KSPSetDM(jac->kspupper, dmInner)); 979 PetscCall(KSPSetDMActive(jac->kspupper, PETSC_FALSE)); 980 PetscCall(KSPSetFromOptions(jac->kspupper)); 981 PetscCall(KSPSetOperators(jac->kspupper,jac->mat[0],jac->pmat[0])); 982 PetscCall(VecDuplicate(jac->head->x, &jac->head->z)); 983 } else { 984 jac->kspupper = jac->head->ksp; 985 PetscCall(PetscObjectReference((PetscObject) jac->head->ksp)); 986 } 987 988 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) { 989 PetscCall(MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp)); 990 } 991 PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc),&jac->kspschur)); 992 PetscCall(KSPSetErrorIfNotConverged(jac->kspschur,pc->erroriffailure)); 993 PetscCall(PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->kspschur)); 994 PetscCall(PetscObjectIncrementTabLevel((PetscObject)jac->kspschur,(PetscObject)pc,1)); 995 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) { 996 PC pcschur; 997 PetscCall(KSPGetPC(jac->kspschur,&pcschur)); 998 PetscCall(PCSetType(pcschur,PCNONE)); 999 /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */ 1000 } else if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_FULL) { 1001 PetscCall(MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user)); 1002 } 1003 PetscCall(KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac))); 1004 PetscCall(KSPGetOptionsPrefix(jac->head->next->ksp, &Dprefix)); 1005 PetscCall(KSPSetOptionsPrefix(jac->kspschur, Dprefix)); 1006 /* propagate DM */ 1007 { 1008 DM sdm; 1009 PetscCall(KSPGetDM(jac->head->next->ksp, &sdm)); 1010 if (sdm) { 1011 PetscCall(KSPSetDM(jac->kspschur, sdm)); 1012 PetscCall(KSPSetDMActive(jac->kspschur, PETSC_FALSE)); 1013 } 1014 } 1015 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 1016 /* need to call this every time, since the jac->kspschur is freshly created, otherwise its options never get set */ 1017 PetscCall(KSPSetFromOptions(jac->kspschur)); 1018 } 1019 PetscCall(MatAssemblyBegin(jac->schur,MAT_FINAL_ASSEMBLY)); 1020 PetscCall(MatAssemblyEnd(jac->schur,MAT_FINAL_ASSEMBLY)); 1021 1022 /* HACK: special support to forward L and Lp matrices that might be used by PCLSC */ 1023 PetscCall(PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_L",ilink->splitname)); 1024 PetscCall(PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L)); 1025 if (!LSC_L) PetscCall(PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L)); 1026 if (LSC_L) PetscCall(PetscObjectCompose((PetscObject)jac->schur,"LSC_L",(PetscObject)LSC_L)); 1027 PetscCall(PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_Lp",ilink->splitname)); 1028 PetscCall(PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L)); 1029 if (!LSC_L) PetscCall(PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L)); 1030 if (LSC_L) PetscCall(PetscObjectCompose((PetscObject)jac->schur,"LSC_Lp",(PetscObject)LSC_L)); 1031 } else if (jac->type == PC_COMPOSITE_GKB) { 1032 IS ccis; 1033 PetscInt rstart,rend; 1034 1035 PetscCheck(nsplit == 2,PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"To use GKB preconditioner you must have exactly 2 fields"); 1036 1037 ilink = jac->head; 1038 1039 /* When extracting off-diagonal submatrices, we take complements from this range */ 1040 PetscCall(MatGetOwnershipRangeColumn(pc->mat,&rstart,&rend)); 1041 1042 PetscCall(ISComplement(ilink->is_col,rstart,rend,&ccis)); 1043 if (jac->offdiag_use_amat) { 1044 PetscCall(MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B)); 1045 } else { 1046 PetscCall(MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B)); 1047 } 1048 PetscCall(ISDestroy(&ccis)); 1049 /* Create work vectors for GKB algorithm */ 1050 PetscCall(VecDuplicate(ilink->x,&jac->u)); 1051 PetscCall(VecDuplicate(ilink->x,&jac->Hu)); 1052 PetscCall(VecDuplicate(ilink->x,&jac->w2)); 1053 ilink = ilink->next; 1054 PetscCall(ISComplement(ilink->is_col,rstart,rend,&ccis)); 1055 if (jac->offdiag_use_amat) { 1056 PetscCall(MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C)); 1057 } else { 1058 PetscCall(MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C)); 1059 } 1060 PetscCall(ISDestroy(&ccis)); 1061 /* Create work vectors for GKB algorithm */ 1062 PetscCall(VecDuplicate(ilink->x,&jac->v)); 1063 PetscCall(VecDuplicate(ilink->x,&jac->d)); 1064 PetscCall(VecDuplicate(ilink->x,&jac->w1)); 1065 PetscCall(MatGolubKahanComputeExplicitOperator(jac->mat[0],jac->B,jac->C,&jac->H,jac->gkbnu)); 1066 PetscCall(PetscCalloc1(jac->gkbdelay,&jac->vecz)); 1067 1068 ilink = jac->head; 1069 PetscCall(KSPSetOperators(ilink->ksp,jac->H,jac->H)); 1070 if (!jac->suboptionsset) PetscCall(KSPSetFromOptions(ilink->ksp)); 1071 /* Create gkb_monitor context */ 1072 if (jac->gkbmonitor) { 1073 PetscInt tablevel; 1074 PetscCall(PetscViewerCreate(PETSC_COMM_WORLD,&jac->gkbviewer)); 1075 PetscCall(PetscViewerSetType(jac->gkbviewer,PETSCVIEWERASCII)); 1076 PetscCall(PetscObjectGetTabLevel((PetscObject)ilink->ksp,&tablevel)); 1077 PetscCall(PetscViewerASCIISetTab(jac->gkbviewer,tablevel)); 1078 PetscCall(PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)ilink->ksp,1)); 1079 } 1080 } else { 1081 /* set up the individual splits' PCs */ 1082 i = 0; 1083 ilink = jac->head; 1084 while (ilink) { 1085 PetscCall(KSPSetOperators(ilink->ksp,jac->mat[i],jac->pmat[i])); 1086 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 1087 if (!jac->suboptionsset) PetscCall(KSPSetFromOptions(ilink->ksp)); 1088 i++; 1089 ilink = ilink->next; 1090 } 1091 } 1092 1093 /* Set coordinates to the sub PC objects whenever these are set */ 1094 if (jac->coordinates_set) { 1095 PC pc_coords; 1096 if (jac->type == PC_COMPOSITE_SCHUR) { 1097 // Head is first block. 1098 PetscCall(KSPGetPC(jac->head->ksp, &pc_coords)); 1099 PetscCall(PCSetCoordinates(pc_coords, jac->head->dim, jac->head->ndofs, jac->head->coords)); 1100 // Second one is Schur block, but its KSP object is in kspschur. 1101 PetscCall(KSPGetPC(jac->kspschur, &pc_coords)); 1102 PetscCall(PCSetCoordinates(pc_coords, jac->head->next->dim, jac->head->next->ndofs, jac->head->next->coords)); 1103 } else if (jac->type == PC_COMPOSITE_GKB) { 1104 PetscCall(PetscInfo(pc, "Warning: Setting coordinates does nothing for the GKB Fieldpslit preconditioner")); 1105 } else { 1106 ilink = jac->head; 1107 while (ilink) { 1108 PetscCall(KSPGetPC(ilink->ksp, &pc_coords)); 1109 PetscCall(PCSetCoordinates(pc_coords, ilink->dim, ilink->ndofs, ilink->coords)); 1110 ilink = ilink->next; 1111 } 1112 } 1113 } 1114 1115 jac->suboptionsset = PETSC_TRUE; 1116 PetscFunctionReturn(0); 1117 } 1118 1119 #define FieldSplitSplitSolveAdd(ilink,xx,yy) \ 1120 (VecScatterBegin(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 1121 VecScatterEnd(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 1122 PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL) ||\ 1123 KSPSolve(ilink->ksp,ilink->x,ilink->y) || \ 1124 KSPCheckSolve(ilink->ksp,pc,ilink->y) || \ 1125 PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL) ||\ 1126 VecScatterBegin(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE) || \ 1127 VecScatterEnd(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE)) 1128 1129 static PetscErrorCode PCApply_FieldSplit_Schur(PC pc,Vec x,Vec y) 1130 { 1131 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1132 PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next; 1133 KSP kspA = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper; 1134 1135 PetscFunctionBegin; 1136 switch (jac->schurfactorization) { 1137 case PC_FIELDSPLIT_SCHUR_FACT_DIAG: 1138 /* [A00 0; 0 -S], positive definite, suitable for MINRES */ 1139 PetscCall(VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD)); 1140 PetscCall(VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD)); 1141 PetscCall(VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD)); 1142 PetscCall(PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL)); 1143 PetscCall(KSPSolve(kspA,ilinkA->x,ilinkA->y)); 1144 PetscCall(KSPCheckSolve(kspA,pc,ilinkA->y)); 1145 PetscCall(PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL)); 1146 PetscCall(VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE)); 1147 PetscCall(VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD)); 1148 PetscCall(PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL)); 1149 PetscCall(KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y)); 1150 PetscCall(KSPCheckSolve(jac->kspschur,pc,ilinkD->y)); 1151 PetscCall(PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL)); 1152 PetscCall(VecScale(ilinkD->y,jac->schurscale)); 1153 PetscCall(VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE)); 1154 PetscCall(VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE)); 1155 PetscCall(VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE)); 1156 break; 1157 case PC_FIELDSPLIT_SCHUR_FACT_LOWER: 1158 /* [A00 0; A10 S], suitable for left preconditioning */ 1159 PetscCall(VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD)); 1160 PetscCall(VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD)); 1161 PetscCall(PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL)); 1162 PetscCall(KSPSolve(kspA,ilinkA->x,ilinkA->y)); 1163 PetscCall(KSPCheckSolve(kspA,pc,ilinkA->y)); 1164 PetscCall(PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL)); 1165 PetscCall(MatMult(jac->C,ilinkA->y,ilinkD->x)); 1166 PetscCall(VecScale(ilinkD->x,-1.)); 1167 PetscCall(VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD)); 1168 PetscCall(VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE)); 1169 PetscCall(VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD)); 1170 PetscCall(PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL)); 1171 PetscCall(KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y)); 1172 PetscCall(KSPCheckSolve(jac->kspschur,pc,ilinkD->y)); 1173 PetscCall(PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL)); 1174 PetscCall(VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE)); 1175 PetscCall(VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE)); 1176 PetscCall(VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE)); 1177 break; 1178 case PC_FIELDSPLIT_SCHUR_FACT_UPPER: 1179 /* [A00 A01; 0 S], suitable for right preconditioning */ 1180 PetscCall(VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD)); 1181 PetscCall(VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD)); 1182 PetscCall(PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL)); 1183 PetscCall(KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y)); 1184 PetscCall(KSPCheckSolve(jac->kspschur,pc,ilinkD->y)); 1185 PetscCall(PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL)); PetscCall(MatMult(jac->B,ilinkD->y,ilinkA->x)); 1186 PetscCall(VecScale(ilinkA->x,-1.)); 1187 PetscCall(VecScatterBegin(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD)); 1188 PetscCall(VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE)); 1189 PetscCall(VecScatterEnd(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD)); 1190 PetscCall(PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL)); 1191 PetscCall(KSPSolve(kspA,ilinkA->x,ilinkA->y)); 1192 PetscCall(KSPCheckSolve(kspA,pc,ilinkA->y)); 1193 PetscCall(PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL)); 1194 PetscCall(VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE)); 1195 PetscCall(VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE)); 1196 PetscCall(VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE)); 1197 break; 1198 case PC_FIELDSPLIT_SCHUR_FACT_FULL: 1199 /* [1 0; A10 A00^{-1} 1] [A00 0; 0 S] [1 A00^{-1}A01; 0 1] */ 1200 PetscCall(VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD)); 1201 PetscCall(VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD)); 1202 PetscCall(PetscLogEventBegin(KSP_Solve_FS_L,kspLower,ilinkA->x,ilinkA->y,NULL)); 1203 PetscCall(KSPSolve(kspLower,ilinkA->x,ilinkA->y)); 1204 PetscCall(KSPCheckSolve(kspLower,pc,ilinkA->y)); 1205 PetscCall(PetscLogEventEnd(KSP_Solve_FS_L,kspLower,ilinkA->x,ilinkA->y,NULL)); 1206 PetscCall(MatMult(jac->C,ilinkA->y,ilinkD->x)); 1207 PetscCall(VecScale(ilinkD->x,-1.0)); 1208 PetscCall(VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD)); 1209 PetscCall(VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD)); 1210 1211 PetscCall(PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL)); 1212 PetscCall(KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y)); 1213 PetscCall(KSPCheckSolve(jac->kspschur,pc,ilinkD->y)); 1214 PetscCall(PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL)); 1215 PetscCall(VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE)); 1216 1217 if (kspUpper == kspA) { 1218 PetscCall(MatMult(jac->B,ilinkD->y,ilinkA->y)); 1219 PetscCall(VecAXPY(ilinkA->x,-1.0,ilinkA->y)); 1220 PetscCall(PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL)); 1221 PetscCall(KSPSolve(kspA,ilinkA->x,ilinkA->y)); 1222 PetscCall(KSPCheckSolve(kspA,pc,ilinkA->y)); 1223 PetscCall(PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL)); 1224 } else { 1225 PetscCall(PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL)); 1226 PetscCall(KSPSolve(kspA,ilinkA->x,ilinkA->y)); 1227 PetscCall(KSPCheckSolve(kspA,pc,ilinkA->y)); 1228 PetscCall(MatMult(jac->B,ilinkD->y,ilinkA->x)); 1229 PetscCall(PetscLogEventBegin(KSP_Solve_FS_U,kspUpper,ilinkA->x,ilinkA->z,NULL)); 1230 PetscCall(KSPSolve(kspUpper,ilinkA->x,ilinkA->z)); 1231 PetscCall(KSPCheckSolve(kspUpper,pc,ilinkA->z)); 1232 PetscCall(PetscLogEventEnd(KSP_Solve_FS_U,kspUpper,ilinkA->x,ilinkA->z,NULL)); 1233 PetscCall(VecAXPY(ilinkA->y,-1.0,ilinkA->z)); 1234 } 1235 PetscCall(VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE)); 1236 PetscCall(VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE)); 1237 PetscCall(VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE)); 1238 } 1239 PetscFunctionReturn(0); 1240 } 1241 1242 static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y) 1243 { 1244 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1245 PC_FieldSplitLink ilink = jac->head; 1246 PetscInt cnt,bs; 1247 1248 PetscFunctionBegin; 1249 if (jac->type == PC_COMPOSITE_ADDITIVE) { 1250 if (jac->defaultsplit) { 1251 PetscCall(VecGetBlockSize(x,&bs)); 1252 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); 1253 PetscCall(VecGetBlockSize(y,&bs)); 1254 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); 1255 PetscCall(VecStrideGatherAll(x,jac->x,INSERT_VALUES)); 1256 while (ilink) { 1257 PetscCall(PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL)); 1258 PetscCall(KSPSolve(ilink->ksp,ilink->x,ilink->y)); 1259 PetscCall(KSPCheckSolve(ilink->ksp,pc,ilink->y)); 1260 PetscCall(PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL)); 1261 ilink = ilink->next; 1262 } 1263 PetscCall(VecStrideScatterAll(jac->y,y,INSERT_VALUES)); 1264 } else { 1265 PetscCall(VecSet(y,0.0)); 1266 while (ilink) { 1267 PetscCall(FieldSplitSplitSolveAdd(ilink,x,y)); 1268 ilink = ilink->next; 1269 } 1270 } 1271 } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE && jac->nsplits == 2) { 1272 PetscCall(VecSet(y,0.0)); 1273 /* solve on first block for first block variables */ 1274 PetscCall(VecScatterBegin(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD)); 1275 PetscCall(VecScatterEnd(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD)); 1276 PetscCall(PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL)); 1277 PetscCall(KSPSolve(ilink->ksp,ilink->x,ilink->y)); 1278 PetscCall(KSPCheckSolve(ilink->ksp,pc,ilink->y)); 1279 PetscCall(PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL)); 1280 PetscCall(VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE)); 1281 PetscCall(VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE)); 1282 1283 /* compute the residual only onto second block variables using first block variables */ 1284 PetscCall(MatMult(jac->Afield[1],ilink->y,ilink->next->x)); 1285 ilink = ilink->next; 1286 PetscCall(VecScale(ilink->x,-1.0)); 1287 PetscCall(VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD)); 1288 PetscCall(VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD)); 1289 1290 /* solve on second block variables */ 1291 PetscCall(PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL)); 1292 PetscCall(KSPSolve(ilink->ksp,ilink->x,ilink->y)); 1293 PetscCall(KSPCheckSolve(ilink->ksp,pc,ilink->y)); 1294 PetscCall(PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL)); 1295 PetscCall(VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE)); 1296 PetscCall(VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE)); 1297 } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 1298 if (!jac->w1) { 1299 PetscCall(VecDuplicate(x,&jac->w1)); 1300 PetscCall(VecDuplicate(x,&jac->w2)); 1301 } 1302 PetscCall(VecSet(y,0.0)); 1303 PetscCall(FieldSplitSplitSolveAdd(ilink,x,y)); 1304 cnt = 1; 1305 while (ilink->next) { 1306 ilink = ilink->next; 1307 /* compute the residual only over the part of the vector needed */ 1308 PetscCall(MatMult(jac->Afield[cnt++],y,ilink->x)); 1309 PetscCall(VecScale(ilink->x,-1.0)); 1310 PetscCall(VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD)); 1311 PetscCall(VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD)); 1312 PetscCall(PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL)); 1313 PetscCall(KSPSolve(ilink->ksp,ilink->x,ilink->y)); 1314 PetscCall(KSPCheckSolve(ilink->ksp,pc,ilink->y)); 1315 PetscCall(PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL)); 1316 PetscCall(VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE)); 1317 PetscCall(VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE)); 1318 } 1319 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 1320 cnt -= 2; 1321 while (ilink->previous) { 1322 ilink = ilink->previous; 1323 /* compute the residual only over the part of the vector needed */ 1324 PetscCall(MatMult(jac->Afield[cnt--],y,ilink->x)); 1325 PetscCall(VecScale(ilink->x,-1.0)); 1326 PetscCall(VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD)); 1327 PetscCall(VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD)); 1328 PetscCall(PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL)); 1329 PetscCall(KSPSolve(ilink->ksp,ilink->x,ilink->y)); 1330 PetscCall(KSPCheckSolve(ilink->ksp,pc,ilink->y)); 1331 PetscCall(PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL)); 1332 PetscCall(VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE)); 1333 PetscCall(VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE)); 1334 } 1335 } 1336 } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Unsupported or unknown composition %d",(int) jac->type); 1337 PetscFunctionReturn(0); 1338 } 1339 1340 static PetscErrorCode PCApply_FieldSplit_GKB(PC pc,Vec x,Vec y) 1341 { 1342 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1343 PC_FieldSplitLink ilinkA = jac->head,ilinkD = ilinkA->next; 1344 KSP ksp = ilinkA->ksp; 1345 Vec u,v,Hu,d,work1,work2; 1346 PetscScalar alpha,z,nrmz2,*vecz; 1347 PetscReal lowbnd,nu,beta; 1348 PetscInt j,iterGKB; 1349 1350 PetscFunctionBegin; 1351 PetscCall(VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD)); 1352 PetscCall(VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD)); 1353 PetscCall(VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD)); 1354 PetscCall(VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD)); 1355 1356 u = jac->u; 1357 v = jac->v; 1358 Hu = jac->Hu; 1359 d = jac->d; 1360 work1 = jac->w1; 1361 work2 = jac->w2; 1362 vecz = jac->vecz; 1363 1364 /* Change RHS to comply with matrix regularization H = A + nu*B*B' */ 1365 /* Add q = q + nu*B*b */ 1366 if (jac->gkbnu) { 1367 nu = jac->gkbnu; 1368 PetscCall(VecScale(ilinkD->x,jac->gkbnu)); 1369 PetscCall(MatMultAdd(jac->B,ilinkD->x,ilinkA->x,ilinkA->x)); /* q = q + nu*B*b */ 1370 } else { 1371 /* Situation when no augmented Lagrangian is used. Then we set inner */ 1372 /* matrix N = I in [Ar13], and thus nu = 1. */ 1373 nu = 1; 1374 } 1375 1376 /* Transform rhs from [q,tilde{b}] to [0,b] */ 1377 PetscCall(PetscLogEventBegin(ilinkA->event,ksp,ilinkA->x,ilinkA->y,NULL)); 1378 PetscCall(KSPSolve(ksp,ilinkA->x,ilinkA->y)); 1379 PetscCall(KSPCheckSolve(ksp,pc,ilinkA->y)); 1380 PetscCall(PetscLogEventEnd(ilinkA->event,ksp,ilinkA->x,ilinkA->y,NULL)); 1381 PetscCall(MatMultHermitianTranspose(jac->B,ilinkA->y,work1)); 1382 PetscCall(VecAXPBY(work1,1.0/nu,-1.0,ilinkD->x)); /* c = b - B'*x */ 1383 1384 /* First step of algorithm */ 1385 PetscCall(VecNorm(work1,NORM_2,&beta)); /* beta = sqrt(nu*c'*c)*/ 1386 KSPCheckDot(ksp,beta); 1387 beta = PetscSqrtReal(nu)*beta; 1388 PetscCall(VecAXPBY(v,nu/beta,0.0,work1)); /* v = nu/beta *c */ 1389 PetscCall(MatMult(jac->B,v,work2)); /* u = H^{-1}*B*v */ 1390 PetscCall(PetscLogEventBegin(ilinkA->event,ksp,work2,u,NULL)); 1391 PetscCall(KSPSolve(ksp,work2,u)); 1392 PetscCall(KSPCheckSolve(ksp,pc,u)); 1393 PetscCall(PetscLogEventEnd(ilinkA->event,ksp,work2,u,NULL)); 1394 PetscCall(MatMult(jac->H,u,Hu)); /* alpha = u'*H*u */ 1395 PetscCall(VecDot(Hu,u,&alpha)); 1396 KSPCheckDot(ksp,alpha); 1397 PetscCheck(PetscRealPart(alpha) > 0.0,PETSC_COMM_SELF,PETSC_ERR_NOT_CONVERGED,"GKB preconditioner diverged, H is not positive definite"); 1398 alpha = PetscSqrtReal(PetscAbsScalar(alpha)); 1399 PetscCall(VecScale(u,1.0/alpha)); 1400 PetscCall(VecAXPBY(d,1.0/alpha,0.0,v)); /* v = nu/beta *c */ 1401 1402 z = beta/alpha; 1403 vecz[1] = z; 1404 1405 /* Computation of first iterate x(1) and p(1) */ 1406 PetscCall(VecAXPY(ilinkA->y,z,u)); 1407 PetscCall(VecCopy(d,ilinkD->y)); 1408 PetscCall(VecScale(ilinkD->y,-z)); 1409 1410 iterGKB = 1; lowbnd = 2*jac->gkbtol; 1411 if (jac->gkbmonitor) { 1412 PetscCall(PetscViewerASCIIPrintf(jac->gkbviewer,"%3" PetscInt_FMT " GKB Lower bound estimate %14.12e\n",iterGKB,(double)lowbnd)); 1413 } 1414 1415 while (iterGKB < jac->gkbmaxit && lowbnd > jac->gkbtol) { 1416 iterGKB += 1; 1417 PetscCall(MatMultHermitianTranspose(jac->B,u,work1)); /* v <- nu*(B'*u-alpha/nu*v) */ 1418 PetscCall(VecAXPBY(v,nu,-alpha,work1)); 1419 PetscCall(VecNorm(v,NORM_2,&beta)); /* beta = sqrt(nu)*v'*v */ 1420 beta = beta/PetscSqrtReal(nu); 1421 PetscCall(VecScale(v,1.0/beta)); 1422 PetscCall(MatMult(jac->B,v,work2)); /* u <- H^{-1}*(B*v-beta*H*u) */ 1423 PetscCall(MatMult(jac->H,u,Hu)); 1424 PetscCall(VecAXPY(work2,-beta,Hu)); 1425 PetscCall(PetscLogEventBegin(ilinkA->event,ksp,work2,u,NULL)); 1426 PetscCall(KSPSolve(ksp,work2,u)); 1427 PetscCall(KSPCheckSolve(ksp,pc,u)); 1428 PetscCall(PetscLogEventEnd(ilinkA->event,ksp,work2,u,NULL)); 1429 PetscCall(MatMult(jac->H,u,Hu)); /* alpha = u'*H*u */ 1430 PetscCall(VecDot(Hu,u,&alpha)); 1431 KSPCheckDot(ksp,alpha); 1432 PetscCheck(PetscRealPart(alpha) > 0.0,PETSC_COMM_SELF,PETSC_ERR_NOT_CONVERGED,"GKB preconditioner diverged, H is not positive definite"); 1433 alpha = PetscSqrtReal(PetscAbsScalar(alpha)); 1434 PetscCall(VecScale(u,1.0/alpha)); 1435 1436 z = -beta/alpha*z; /* z <- beta/alpha*z */ 1437 vecz[0] = z; 1438 1439 /* Computation of new iterate x(i+1) and p(i+1) */ 1440 PetscCall(VecAXPBY(d,1.0/alpha,-beta/alpha,v)); /* d = (v-beta*d)/alpha */ 1441 PetscCall(VecAXPY(ilinkA->y,z,u)); /* r = r + z*u */ 1442 PetscCall(VecAXPY(ilinkD->y,-z,d)); /* p = p - z*d */ 1443 PetscCall(MatMult(jac->H,ilinkA->y,Hu)); /* ||u||_H = u'*H*u */ 1444 PetscCall(VecDot(Hu,ilinkA->y,&nrmz2)); 1445 1446 /* Compute Lower Bound estimate */ 1447 if (iterGKB > jac->gkbdelay) { 1448 lowbnd = 0.0; 1449 for (j=0; j<jac->gkbdelay; j++) { 1450 lowbnd += PetscAbsScalar(vecz[j]*vecz[j]); 1451 } 1452 lowbnd = PetscSqrtReal(lowbnd/PetscAbsScalar(nrmz2)); 1453 } 1454 1455 for (j=0; j<jac->gkbdelay-1; j++) { 1456 vecz[jac->gkbdelay-j-1] = vecz[jac->gkbdelay-j-2]; 1457 } 1458 if (jac->gkbmonitor) { 1459 PetscCall(PetscViewerASCIIPrintf(jac->gkbviewer,"%3" PetscInt_FMT " GKB Lower bound estimate %14.12e\n",iterGKB,(double)lowbnd)); 1460 } 1461 } 1462 1463 PetscCall(VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE)); 1464 PetscCall(VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE)); 1465 PetscCall(VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE)); 1466 PetscCall(VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE)); 1467 1468 PetscFunctionReturn(0); 1469 } 1470 1471 #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \ 1472 (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 1473 VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 1474 PetscLogEventBegin(ilink->event,ilink->ksp,ilink->y,ilink->x,NULL) || \ 1475 KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) || \ 1476 KSPCheckSolve(ilink->ksp,pc,ilink->x) || \ 1477 PetscLogEventEnd(ilink->event,ilink->ksp,ilink->y,ilink->x,NULL) || \ 1478 VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \ 1479 VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE)) 1480 1481 static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y) 1482 { 1483 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1484 PC_FieldSplitLink ilink = jac->head; 1485 PetscInt bs; 1486 1487 PetscFunctionBegin; 1488 if (jac->type == PC_COMPOSITE_ADDITIVE) { 1489 if (jac->defaultsplit) { 1490 PetscCall(VecGetBlockSize(x,&bs)); 1491 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); 1492 PetscCall(VecGetBlockSize(y,&bs)); 1493 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); 1494 PetscCall(VecStrideGatherAll(x,jac->x,INSERT_VALUES)); 1495 while (ilink) { 1496 PetscCall(PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL)); 1497 PetscCall(KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y)); 1498 PetscCall(KSPCheckSolve(ilink->ksp,pc,ilink->y)); 1499 PetscCall(PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL)); 1500 ilink = ilink->next; 1501 } 1502 PetscCall(VecStrideScatterAll(jac->y,y,INSERT_VALUES)); 1503 } else { 1504 PetscCall(VecSet(y,0.0)); 1505 while (ilink) { 1506 PetscCall(FieldSplitSplitSolveAddTranspose(ilink,x,y)); 1507 ilink = ilink->next; 1508 } 1509 } 1510 } else { 1511 if (!jac->w1) { 1512 PetscCall(VecDuplicate(x,&jac->w1)); 1513 PetscCall(VecDuplicate(x,&jac->w2)); 1514 } 1515 PetscCall(VecSet(y,0.0)); 1516 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 1517 PetscCall(FieldSplitSplitSolveAddTranspose(ilink,x,y)); 1518 while (ilink->next) { 1519 ilink = ilink->next; 1520 PetscCall(MatMultTranspose(pc->mat,y,jac->w1)); 1521 PetscCall(VecWAXPY(jac->w2,-1.0,jac->w1,x)); 1522 PetscCall(FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y)); 1523 } 1524 while (ilink->previous) { 1525 ilink = ilink->previous; 1526 PetscCall(MatMultTranspose(pc->mat,y,jac->w1)); 1527 PetscCall(VecWAXPY(jac->w2,-1.0,jac->w1,x)); 1528 PetscCall(FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y)); 1529 } 1530 } else { 1531 while (ilink->next) { /* get to last entry in linked list */ 1532 ilink = ilink->next; 1533 } 1534 PetscCall(FieldSplitSplitSolveAddTranspose(ilink,x,y)); 1535 while (ilink->previous) { 1536 ilink = ilink->previous; 1537 PetscCall(MatMultTranspose(pc->mat,y,jac->w1)); 1538 PetscCall(VecWAXPY(jac->w2,-1.0,jac->w1,x)); 1539 PetscCall(FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y)); 1540 } 1541 } 1542 } 1543 PetscFunctionReturn(0); 1544 } 1545 1546 static PetscErrorCode PCReset_FieldSplit(PC pc) 1547 { 1548 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1549 PC_FieldSplitLink ilink = jac->head,next; 1550 1551 PetscFunctionBegin; 1552 while (ilink) { 1553 PetscCall(KSPDestroy(&ilink->ksp)); 1554 PetscCall(VecDestroy(&ilink->x)); 1555 PetscCall(VecDestroy(&ilink->y)); 1556 PetscCall(VecDestroy(&ilink->z)); 1557 PetscCall(VecScatterDestroy(&ilink->sctx)); 1558 PetscCall(ISDestroy(&ilink->is)); 1559 PetscCall(ISDestroy(&ilink->is_col)); 1560 PetscCall(PetscFree(ilink->splitname)); 1561 PetscCall(PetscFree(ilink->fields)); 1562 PetscCall(PetscFree(ilink->fields_col)); 1563 next = ilink->next; 1564 PetscCall(PetscFree(ilink)); 1565 ilink = next; 1566 } 1567 jac->head = NULL; 1568 PetscCall(PetscFree2(jac->x,jac->y)); 1569 if (jac->mat && jac->mat != jac->pmat) { 1570 PetscCall(MatDestroyMatrices(jac->nsplits,&jac->mat)); 1571 } else if (jac->mat) { 1572 jac->mat = NULL; 1573 } 1574 if (jac->pmat) PetscCall(MatDestroyMatrices(jac->nsplits,&jac->pmat)); 1575 if (jac->Afield) PetscCall(MatDestroyMatrices(jac->nsplits,&jac->Afield)); 1576 jac->nsplits = 0; 1577 PetscCall(VecDestroy(&jac->w1)); 1578 PetscCall(VecDestroy(&jac->w2)); 1579 PetscCall(MatDestroy(&jac->schur)); 1580 PetscCall(MatDestroy(&jac->schurp)); 1581 PetscCall(MatDestroy(&jac->schur_user)); 1582 PetscCall(KSPDestroy(&jac->kspschur)); 1583 PetscCall(KSPDestroy(&jac->kspupper)); 1584 PetscCall(MatDestroy(&jac->B)); 1585 PetscCall(MatDestroy(&jac->C)); 1586 PetscCall(MatDestroy(&jac->H)); 1587 PetscCall(VecDestroy(&jac->u)); 1588 PetscCall(VecDestroy(&jac->v)); 1589 PetscCall(VecDestroy(&jac->Hu)); 1590 PetscCall(VecDestroy(&jac->d)); 1591 PetscCall(PetscFree(jac->vecz)); 1592 PetscCall(PetscViewerDestroy(&jac->gkbviewer)); 1593 jac->isrestrict = PETSC_FALSE; 1594 PetscFunctionReturn(0); 1595 } 1596 1597 static PetscErrorCode PCDestroy_FieldSplit(PC pc) 1598 { 1599 PetscFunctionBegin; 1600 PetscCall(PCReset_FieldSplit(pc)); 1601 PetscCall(PetscFree(pc->data)); 1602 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",NULL)); 1603 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",NULL)); 1604 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",NULL)); 1605 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",NULL)); 1606 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",NULL)); 1607 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitRestrictIS_C",NULL)); 1608 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSchurGetSubKSP_C",NULL)); 1609 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",NULL)); 1610 1611 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBTol_C",NULL)); 1612 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBMaxit_C",NULL)); 1613 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBNu_C",NULL)); 1614 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBDelay_C",NULL)); 1615 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",NULL)); 1616 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",NULL)); 1617 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",NULL)); 1618 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",NULL)); 1619 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurScale_C",NULL)); 1620 PetscFunctionReturn(0); 1621 } 1622 1623 static PetscErrorCode PCSetFromOptions_FieldSplit(PetscOptionItems *PetscOptionsObject,PC pc) 1624 { 1625 PetscInt bs; 1626 PetscBool flg; 1627 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1628 PCCompositeType ctype; 1629 1630 PetscFunctionBegin; 1631 PetscOptionsHeadBegin(PetscOptionsObject,"FieldSplit options"); 1632 PetscCall(PetscOptionsBool("-pc_fieldsplit_dm_splits","Whether to use DMCreateFieldDecomposition() for splits","PCFieldSplitSetDMSplits",jac->dm_splits,&jac->dm_splits,NULL)); 1633 PetscCall(PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg)); 1634 if (flg) PetscCall(PCFieldSplitSetBlockSize(pc,bs)); 1635 jac->diag_use_amat = pc->useAmat; 1636 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)); 1637 jac->offdiag_use_amat = pc->useAmat; 1638 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)); 1639 PetscCall(PetscOptionsBool("-pc_fieldsplit_detect_saddle_point","Form 2-way split by detecting zero diagonal entries", "PCFieldSplitSetDetectSaddlePoint",jac->detect,&jac->detect,NULL)); 1640 PetscCall(PCFieldSplitSetDetectSaddlePoint(pc,jac->detect)); /* Sets split type and Schur PC type */ 1641 PetscCall(PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg)); 1642 if (flg) PetscCall(PCFieldSplitSetType(pc,ctype)); 1643 /* Only setup fields once */ 1644 if ((jac->bs > 0) && (jac->nsplits == 0)) { 1645 /* only allow user to set fields from command line if bs is already known. 1646 otherwise user can set them in PCFieldSplitSetDefaults() */ 1647 PetscCall(PCFieldSplitSetRuntimeSplits_Private(pc)); 1648 if (jac->splitdefined) PetscCall(PetscInfo(pc,"Splits defined using the options database\n")); 1649 } 1650 if (jac->type == PC_COMPOSITE_SCHUR) { 1651 PetscCall(PetscOptionsGetEnum(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_schur_factorization_type",PCFieldSplitSchurFactTypes,(PetscEnum*)&jac->schurfactorization,&flg)); 1652 if (flg) PetscCall(PetscInfo(pc,"Deprecated use of -pc_fieldsplit_schur_factorization_type\n")); 1653 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)); 1654 PetscCall(PetscOptionsEnum("-pc_fieldsplit_schur_precondition","How to build preconditioner for Schur complement","PCFieldSplitSetSchurPre",PCFieldSplitSchurPreTypes,(PetscEnum)jac->schurpre,(PetscEnum*)&jac->schurpre,NULL)); 1655 PetscCall(PetscOptionsScalar("-pc_fieldsplit_schur_scale","Scale Schur complement","PCFieldSplitSetSchurScale",jac->schurscale,&jac->schurscale,NULL)); 1656 } else if (jac->type == PC_COMPOSITE_GKB) { 1657 PetscCall(PetscOptionsReal("-pc_fieldsplit_gkb_tol","The tolerance for the lower bound stopping criterion","PCFieldSplitGKBTol",jac->gkbtol,&jac->gkbtol,NULL)); 1658 PetscCall(PetscOptionsInt("-pc_fieldsplit_gkb_delay","The delay value for lower bound criterion","PCFieldSplitGKBDelay",jac->gkbdelay,&jac->gkbdelay,NULL)); 1659 PetscCall(PetscOptionsReal("-pc_fieldsplit_gkb_nu","Parameter in augmented Lagrangian approach","PCFieldSplitGKBNu",jac->gkbnu,&jac->gkbnu,NULL)); 1660 PetscCheck(jac->gkbnu >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nu cannot be less than 0: value %g",(double)jac->gkbnu); 1661 PetscCall(PetscOptionsInt("-pc_fieldsplit_gkb_maxit","Maximum allowed number of iterations","PCFieldSplitGKBMaxit",jac->gkbmaxit,&jac->gkbmaxit,NULL)); 1662 PetscCall(PetscOptionsBool("-pc_fieldsplit_gkb_monitor","Prints number of GKB iterations and error","PCFieldSplitGKB",jac->gkbmonitor,&jac->gkbmonitor,NULL)); 1663 } 1664 /* 1665 In the initial call to this routine the sub-solver data structures do not exist so we cannot call KSPSetFromOptions() on them yet. 1666 But after the initial setup of ALL the layers of sub-solvers is completed we do want to call KSPSetFromOptions() on the sub-solvers every time it 1667 is called on the outer solver in case changes were made in the options database 1668 1669 But even after PCSetUp_FieldSplit() is called all the options inside the inner levels of sub-solvers may still not have been set thus we only call the KSPSetFromOptions() 1670 if we know that the entire stack of sub-solvers below this have been complete instantiated, we check this by seeing if any solver iterations are complete. 1671 Without this extra check test p2p1fetidp_olof_full and others fail with incorrect matrix types. 1672 1673 There could be a negative side effect of calling the KSPSetFromOptions() below. 1674 1675 If one captured the PetscObjectState of the options database one could skip these calls if the database has not changed from the previous call 1676 */ 1677 if (jac->issetup) { 1678 PC_FieldSplitLink ilink = jac->head; 1679 if (jac->type == PC_COMPOSITE_SCHUR) { 1680 if (jac->kspupper && jac->kspupper->totalits > 0) PetscCall(KSPSetFromOptions(jac->kspupper)); 1681 if (jac->kspschur && jac->kspschur->totalits > 0) PetscCall(KSPSetFromOptions(jac->kspschur)); 1682 } 1683 while (ilink) { 1684 if (ilink->ksp->totalits > 0) PetscCall(KSPSetFromOptions(ilink->ksp)); 1685 ilink = ilink->next; 1686 } 1687 } 1688 PetscOptionsHeadEnd(); 1689 PetscFunctionReturn(0); 1690 } 1691 1692 /*------------------------------------------------------------------------------------*/ 1693 1694 static PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 1695 { 1696 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1697 PC_FieldSplitLink ilink,next = jac->head; 1698 char prefix[128]; 1699 PetscInt i; 1700 1701 PetscFunctionBegin; 1702 if (jac->splitdefined) { 1703 PetscCall(PetscInfo(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname)); 1704 PetscFunctionReturn(0); 1705 } 1706 for (i=0; i<n; i++) { 1707 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); 1708 PetscCheck(fields[i] >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %" PetscInt_FMT " requested",fields[i]); 1709 } 1710 PetscCall(PetscNew(&ilink)); 1711 if (splitname) { 1712 PetscCall(PetscStrallocpy(splitname,&ilink->splitname)); 1713 } else { 1714 PetscCall(PetscMalloc1(3,&ilink->splitname)); 1715 PetscCall(PetscSNPrintf(ilink->splitname,2,"%" PetscInt_FMT,jac->nsplits)); 1716 } 1717 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 */ 1718 PetscCall(PetscMalloc1(n,&ilink->fields)); 1719 PetscCall(PetscArraycpy(ilink->fields,fields,n)); 1720 PetscCall(PetscMalloc1(n,&ilink->fields_col)); 1721 PetscCall(PetscArraycpy(ilink->fields_col,fields_col,n)); 1722 1723 ilink->nfields = n; 1724 ilink->next = NULL; 1725 PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp)); 1726 PetscCall(KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure)); 1727 PetscCall(PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1)); 1728 PetscCall(KSPSetType(ilink->ksp,KSPPREONLY)); 1729 PetscCall(PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp)); 1730 1731 PetscCall(PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname)); 1732 PetscCall(KSPSetOptionsPrefix(ilink->ksp,prefix)); 1733 1734 if (!next) { 1735 jac->head = ilink; 1736 ilink->previous = NULL; 1737 } else { 1738 while (next->next) { 1739 next = next->next; 1740 } 1741 next->next = ilink; 1742 ilink->previous = next; 1743 } 1744 jac->nsplits++; 1745 PetscFunctionReturn(0); 1746 } 1747 1748 static PetscErrorCode PCFieldSplitSchurGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp) 1749 { 1750 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1751 1752 PetscFunctionBegin; 1753 *subksp = NULL; 1754 if (n) *n = 0; 1755 if (jac->type == PC_COMPOSITE_SCHUR) { 1756 PetscInt nn; 1757 1758 PetscCheck(jac->schur,PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call KSPSetUp() or PCSetUp() before calling PCFieldSplitSchurGetSubKSP()"); 1759 PetscCheck(jac->nsplits == 2,PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Unexpected number of splits %" PetscInt_FMT " != 2",jac->nsplits); 1760 nn = jac->nsplits + (jac->kspupper != jac->head->ksp ? 1 : 0); 1761 PetscCall(PetscMalloc1(nn,subksp)); 1762 (*subksp)[0] = jac->head->ksp; 1763 (*subksp)[1] = jac->kspschur; 1764 if (jac->kspupper != jac->head->ksp) (*subksp)[2] = jac->kspupper; 1765 if (n) *n = nn; 1766 } 1767 PetscFunctionReturn(0); 1768 } 1769 1770 static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp) 1771 { 1772 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1773 1774 PetscFunctionBegin; 1775 PetscCheck(jac->schur,PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call KSPSetUp() or PCSetUp() before calling PCFieldSplitGetSubKSP()"); 1776 PetscCall(PetscMalloc1(jac->nsplits,subksp)); 1777 PetscCall(MatSchurComplementGetKSP(jac->schur,*subksp)); 1778 1779 (*subksp)[1] = jac->kspschur; 1780 if (n) *n = jac->nsplits; 1781 PetscFunctionReturn(0); 1782 } 1783 1784 static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp) 1785 { 1786 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1787 PetscInt cnt = 0; 1788 PC_FieldSplitLink ilink = jac->head; 1789 1790 PetscFunctionBegin; 1791 PetscCall(PetscMalloc1(jac->nsplits,subksp)); 1792 while (ilink) { 1793 (*subksp)[cnt++] = ilink->ksp; 1794 ilink = ilink->next; 1795 } 1796 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); 1797 if (n) *n = jac->nsplits; 1798 PetscFunctionReturn(0); 1799 } 1800 1801 /*@C 1802 PCFieldSplitRestrictIS - Restricts the fieldsplit ISs to be within a given IS. 1803 1804 Input Parameters: 1805 + pc - the preconditioner context 1806 - is - the index set that defines the indices to which the fieldsplit is to be restricted 1807 1808 Level: advanced 1809 1810 @*/ 1811 PetscErrorCode PCFieldSplitRestrictIS(PC pc,IS isy) 1812 { 1813 PetscFunctionBegin; 1814 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1815 PetscValidHeaderSpecific(isy,IS_CLASSID,2); 1816 PetscTryMethod(pc,"PCFieldSplitRestrictIS_C",(PC,IS),(pc,isy)); 1817 PetscFunctionReturn(0); 1818 } 1819 1820 static PetscErrorCode PCFieldSplitRestrictIS_FieldSplit(PC pc, IS isy) 1821 { 1822 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1823 PC_FieldSplitLink ilink = jac->head, next; 1824 PetscInt localsize,size,sizez,i; 1825 const PetscInt *ind, *indz; 1826 PetscInt *indc, *indcz; 1827 PetscBool flg; 1828 1829 PetscFunctionBegin; 1830 PetscCall(ISGetLocalSize(isy,&localsize)); 1831 PetscCallMPI(MPI_Scan(&localsize,&size,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)isy))); 1832 size -= localsize; 1833 while (ilink) { 1834 IS isrl,isr; 1835 PC subpc; 1836 PetscCall(ISEmbed(ilink->is, isy, PETSC_TRUE, &isrl)); 1837 PetscCall(ISGetLocalSize(isrl,&localsize)); 1838 PetscCall(PetscMalloc1(localsize,&indc)); 1839 PetscCall(ISGetIndices(isrl,&ind)); 1840 PetscCall(PetscArraycpy(indc,ind,localsize)); 1841 PetscCall(ISRestoreIndices(isrl,&ind)); 1842 PetscCall(ISDestroy(&isrl)); 1843 for (i=0; i<localsize; i++) *(indc+i) += size; 1844 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)isy),localsize,indc,PETSC_OWN_POINTER,&isr)); 1845 PetscCall(PetscObjectReference((PetscObject)isr)); 1846 PetscCall(ISDestroy(&ilink->is)); 1847 ilink->is = isr; 1848 PetscCall(PetscObjectReference((PetscObject)isr)); 1849 PetscCall(ISDestroy(&ilink->is_col)); 1850 ilink->is_col = isr; 1851 PetscCall(ISDestroy(&isr)); 1852 PetscCall(KSPGetPC(ilink->ksp, &subpc)); 1853 PetscCall(PetscObjectTypeCompare((PetscObject)subpc,PCFIELDSPLIT,&flg)); 1854 if (flg) { 1855 IS iszl,isz; 1856 MPI_Comm comm; 1857 PetscCall(ISGetLocalSize(ilink->is,&localsize)); 1858 comm = PetscObjectComm((PetscObject)ilink->is); 1859 PetscCall(ISEmbed(isy, ilink->is, PETSC_TRUE, &iszl)); 1860 PetscCallMPI(MPI_Scan(&localsize,&sizez,1,MPIU_INT,MPI_SUM,comm)); 1861 sizez -= localsize; 1862 PetscCall(ISGetLocalSize(iszl,&localsize)); 1863 PetscCall(PetscMalloc1(localsize,&indcz)); 1864 PetscCall(ISGetIndices(iszl,&indz)); 1865 PetscCall(PetscArraycpy(indcz,indz,localsize)); 1866 PetscCall(ISRestoreIndices(iszl,&indz)); 1867 PetscCall(ISDestroy(&iszl)); 1868 for (i=0; i<localsize; i++) *(indcz+i) += sizez; 1869 PetscCall(ISCreateGeneral(comm,localsize,indcz,PETSC_OWN_POINTER,&isz)); 1870 PetscCall(PCFieldSplitRestrictIS(subpc,isz)); 1871 PetscCall(ISDestroy(&isz)); 1872 } 1873 next = ilink->next; 1874 ilink = next; 1875 } 1876 jac->isrestrict = PETSC_TRUE; 1877 PetscFunctionReturn(0); 1878 } 1879 1880 static PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is) 1881 { 1882 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1883 PC_FieldSplitLink ilink, next = jac->head; 1884 char prefix[128]; 1885 1886 PetscFunctionBegin; 1887 if (jac->splitdefined) { 1888 PetscCall(PetscInfo(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname)); 1889 PetscFunctionReturn(0); 1890 } 1891 PetscCall(PetscNew(&ilink)); 1892 if (splitname) { 1893 PetscCall(PetscStrallocpy(splitname,&ilink->splitname)); 1894 } else { 1895 PetscCall(PetscMalloc1(8,&ilink->splitname)); 1896 PetscCall(PetscSNPrintf(ilink->splitname,7,"%" PetscInt_FMT,jac->nsplits)); 1897 } 1898 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 */ 1899 PetscCall(PetscObjectReference((PetscObject)is)); 1900 PetscCall(ISDestroy(&ilink->is)); 1901 ilink->is = is; 1902 PetscCall(PetscObjectReference((PetscObject)is)); 1903 PetscCall(ISDestroy(&ilink->is_col)); 1904 ilink->is_col = is; 1905 ilink->next = NULL; 1906 PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp)); 1907 PetscCall(KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure)); 1908 PetscCall(PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1)); 1909 PetscCall(KSPSetType(ilink->ksp,KSPPREONLY)); 1910 PetscCall(PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp)); 1911 1912 PetscCall(PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname)); 1913 PetscCall(KSPSetOptionsPrefix(ilink->ksp,prefix)); 1914 1915 if (!next) { 1916 jac->head = ilink; 1917 ilink->previous = NULL; 1918 } else { 1919 while (next->next) { 1920 next = next->next; 1921 } 1922 next->next = ilink; 1923 ilink->previous = next; 1924 } 1925 jac->nsplits++; 1926 PetscFunctionReturn(0); 1927 } 1928 1929 /*@C 1930 PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner 1931 1932 Logically Collective on PC 1933 1934 Input Parameters: 1935 + pc - the preconditioner context 1936 . splitname - name of this split, if NULL the number of the split is used 1937 . n - the number of fields in this split 1938 - fields - the fields in this split 1939 1940 Level: intermediate 1941 1942 Notes: 1943 Use PCFieldSplitSetIS() to set a completely general set of indices as a field. 1944 1945 The PCFieldSplitSetFields() is for defining fields as strided blocks. For example, if the block 1946 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 1947 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x.... 1948 where the numbered entries indicate what is in the field. 1949 1950 This function is called once per split (it creates a new split each time). Solve options 1951 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1952 1953 Developer Note: This routine does not actually create the IS representing the split, that is delayed 1954 until PCSetUp_FieldSplit(), because information about the vector/matrix layouts may not be 1955 available when this routine is called. 1956 1957 .seealso: `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetBlockSize()`, `PCFieldSplitSetIS()` 1958 1959 @*/ 1960 PetscErrorCode PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 1961 { 1962 PetscFunctionBegin; 1963 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1964 PetscValidCharPointer(splitname,2); 1965 PetscCheck(n >= 1,PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Provided number of fields %" PetscInt_FMT " in split \"%s\" not positive",n,splitname); 1966 PetscValidIntPointer(fields,4); 1967 PetscTryMethod(pc,"PCFieldSplitSetFields_C",(PC,const char[],PetscInt,const PetscInt*,const PetscInt*),(pc,splitname,n,fields,fields_col)); 1968 PetscFunctionReturn(0); 1969 } 1970 1971 /*@ 1972 PCFieldSplitSetDiagUseAmat - set flag indicating whether to extract diagonal blocks from Amat (rather than Pmat) 1973 1974 Logically Collective on PC 1975 1976 Input Parameters: 1977 + pc - the preconditioner object 1978 - flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from 1979 1980 Options Database: 1981 . -pc_fieldsplit_diag_use_amat - use the Amat to provide the diagonal blocks 1982 1983 Level: intermediate 1984 1985 .seealso: `PCFieldSplitGetDiagUseAmat()`, `PCFieldSplitSetOffDiagUseAmat()`, `PCFIELDSPLIT` 1986 1987 @*/ 1988 PetscErrorCode PCFieldSplitSetDiagUseAmat(PC pc,PetscBool flg) 1989 { 1990 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1991 PetscBool isfs; 1992 1993 PetscFunctionBegin; 1994 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1995 PetscCall(PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs)); 1996 PetscCheck(isfs,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 1997 jac->diag_use_amat = flg; 1998 PetscFunctionReturn(0); 1999 } 2000 2001 /*@ 2002 PCFieldSplitGetDiagUseAmat - get the flag indicating whether to extract diagonal blocks from Amat (rather than Pmat) 2003 2004 Logically Collective on PC 2005 2006 Input Parameters: 2007 . pc - the preconditioner object 2008 2009 Output Parameters: 2010 . flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from 2011 2012 Level: intermediate 2013 2014 .seealso: `PCFieldSplitSetDiagUseAmat()`, `PCFieldSplitGetOffDiagUseAmat()`, `PCFIELDSPLIT` 2015 2016 @*/ 2017 PetscErrorCode PCFieldSplitGetDiagUseAmat(PC pc,PetscBool *flg) 2018 { 2019 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2020 PetscBool isfs; 2021 2022 PetscFunctionBegin; 2023 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2024 PetscValidBoolPointer(flg,2); 2025 PetscCall(PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs)); 2026 PetscCheck(isfs,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 2027 *flg = jac->diag_use_amat; 2028 PetscFunctionReturn(0); 2029 } 2030 2031 /*@ 2032 PCFieldSplitSetOffDiagUseAmat - set flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat) 2033 2034 Logically Collective on PC 2035 2036 Input Parameters: 2037 + pc - the preconditioner object 2038 - flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from 2039 2040 Options Database: 2041 . -pc_fieldsplit_off_diag_use_amat <bool> - use the Amat to extract the off-diagonal blocks 2042 2043 Level: intermediate 2044 2045 .seealso: `PCFieldSplitGetOffDiagUseAmat()`, `PCFieldSplitSetDiagUseAmat()`, `PCFIELDSPLIT` 2046 2047 @*/ 2048 PetscErrorCode PCFieldSplitSetOffDiagUseAmat(PC pc,PetscBool flg) 2049 { 2050 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2051 PetscBool isfs; 2052 2053 PetscFunctionBegin; 2054 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2055 PetscCall(PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs)); 2056 PetscCheck(isfs,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 2057 jac->offdiag_use_amat = flg; 2058 PetscFunctionReturn(0); 2059 } 2060 2061 /*@ 2062 PCFieldSplitGetOffDiagUseAmat - get the flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat) 2063 2064 Logically Collective on PC 2065 2066 Input Parameters: 2067 . pc - the preconditioner object 2068 2069 Output Parameters: 2070 . flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from 2071 2072 Level: intermediate 2073 2074 .seealso: `PCFieldSplitSetOffDiagUseAmat()`, `PCFieldSplitGetDiagUseAmat()`, `PCFIELDSPLIT` 2075 2076 @*/ 2077 PetscErrorCode PCFieldSplitGetOffDiagUseAmat(PC pc,PetscBool *flg) 2078 { 2079 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2080 PetscBool isfs; 2081 2082 PetscFunctionBegin; 2083 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2084 PetscValidBoolPointer(flg,2); 2085 PetscCall(PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs)); 2086 PetscCheck(isfs,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 2087 *flg = jac->offdiag_use_amat; 2088 PetscFunctionReturn(0); 2089 } 2090 2091 /*@C 2092 PCFieldSplitSetIS - Sets the exact elements for field 2093 2094 Logically Collective on PC 2095 2096 Input Parameters: 2097 + pc - the preconditioner context 2098 . splitname - name of this split, if NULL the number of the split is used 2099 - is - the index set that defines the vector elements in this field 2100 2101 Notes: 2102 Use PCFieldSplitSetFields(), for fields defined by strided types. 2103 2104 This function is called once per split (it creates a new split each time). Solve options 2105 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 2106 2107 Level: intermediate 2108 2109 .seealso: `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetBlockSize()` 2110 2111 @*/ 2112 PetscErrorCode PCFieldSplitSetIS(PC pc,const char splitname[],IS is) 2113 { 2114 PetscFunctionBegin; 2115 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2116 if (splitname) PetscValidCharPointer(splitname,2); 2117 PetscValidHeaderSpecific(is,IS_CLASSID,3); 2118 PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is)); 2119 PetscFunctionReturn(0); 2120 } 2121 2122 /*@C 2123 PCFieldSplitGetIS - Retrieves the elements for a field as an IS 2124 2125 Logically Collective on PC 2126 2127 Input Parameters: 2128 + pc - the preconditioner context 2129 - splitname - name of this split 2130 2131 Output Parameter: 2132 - is - the index set that defines the vector elements in this field, or NULL if the field is not found 2133 2134 Level: intermediate 2135 2136 .seealso: `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetIS()` 2137 2138 @*/ 2139 PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is) 2140 { 2141 PetscFunctionBegin; 2142 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2143 PetscValidCharPointer(splitname,2); 2144 PetscValidPointer(is,3); 2145 { 2146 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 2147 PC_FieldSplitLink ilink = jac->head; 2148 PetscBool found; 2149 2150 *is = NULL; 2151 while (ilink) { 2152 PetscCall(PetscStrcmp(ilink->splitname, splitname, &found)); 2153 if (found) { 2154 *is = ilink->is; 2155 break; 2156 } 2157 ilink = ilink->next; 2158 } 2159 } 2160 PetscFunctionReturn(0); 2161 } 2162 2163 /*@C 2164 PCFieldSplitGetISByIndex - Retrieves the elements for a given index field as an IS 2165 2166 Logically Collective on PC 2167 2168 Input Parameters: 2169 + pc - the preconditioner context 2170 - index - index of this split 2171 2172 Output Parameter: 2173 - is - the index set that defines the vector elements in this field 2174 2175 Level: intermediate 2176 2177 .seealso: `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitGetIS()`, `PCFieldSplitSetIS()` 2178 2179 @*/ 2180 PetscErrorCode PCFieldSplitGetISByIndex(PC pc,PetscInt index,IS *is) 2181 { 2182 PetscFunctionBegin; 2183 PetscCheck(index >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %" PetscInt_FMT " requested",index); 2184 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2185 PetscValidPointer(is,3); 2186 { 2187 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 2188 PC_FieldSplitLink ilink = jac->head; 2189 PetscInt i = 0; 2190 PetscCheck(index < jac->nsplits,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %" PetscInt_FMT " requested but only %" PetscInt_FMT " exist",index,jac->nsplits); 2191 2192 while (i < index) { 2193 ilink = ilink->next; 2194 ++i; 2195 } 2196 PetscCall(PCFieldSplitGetIS(pc,ilink->splitname,is)); 2197 } 2198 PetscFunctionReturn(0); 2199 } 2200 2201 /*@ 2202 PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the 2203 fieldsplit preconditioner. If not set the matrix block size is used. 2204 2205 Logically Collective on PC 2206 2207 Input Parameters: 2208 + pc - the preconditioner context 2209 - bs - the block size 2210 2211 Level: intermediate 2212 2213 .seealso: `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()` 2214 2215 @*/ 2216 PetscErrorCode PCFieldSplitSetBlockSize(PC pc,PetscInt bs) 2217 { 2218 PetscFunctionBegin; 2219 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2220 PetscValidLogicalCollectiveInt(pc,bs,2); 2221 PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs)); 2222 PetscFunctionReturn(0); 2223 } 2224 2225 /*@C 2226 PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits 2227 2228 Collective on KSP 2229 2230 Input Parameter: 2231 . pc - the preconditioner context 2232 2233 Output Parameters: 2234 + n - the number of splits 2235 - subksp - the array of KSP contexts 2236 2237 Note: 2238 After PCFieldSplitGetSubKSP() the array of KSPs is to be freed by the user with PetscFree() 2239 (not the KSP just the array that contains them). 2240 2241 You must call PCSetUp() before calling PCFieldSplitGetSubKSP(). 2242 2243 If the fieldsplit is of type PC_COMPOSITE_SCHUR, it returns the KSP object used inside the 2244 Schur complement and the KSP object used to iterate over the Schur complement. 2245 To access all the KSP objects used in PC_COMPOSITE_SCHUR, use PCFieldSplitSchurGetSubKSP(). 2246 2247 If the fieldsplit is of type PC_COMPOSITE_GKB, it returns the KSP object used to solve the 2248 inner linear system defined by the matrix H in each loop. 2249 2250 Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs. 2251 You can call PCFieldSplitGetSubKSP(pc,n,PETSC_NULL_KSP,ierr) to determine how large the 2252 KSP array must be. 2253 2254 Level: advanced 2255 2256 .seealso: `PCFIELDSPLIT` 2257 @*/ 2258 PetscErrorCode PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[]) 2259 { 2260 PetscFunctionBegin; 2261 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2262 if (n) PetscValidIntPointer(n,2); 2263 PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp)); 2264 PetscFunctionReturn(0); 2265 } 2266 2267 /*@C 2268 PCFieldSplitSchurGetSubKSP - Gets the KSP contexts used inside the Schur complement based PCFIELDSPLIT 2269 2270 Collective on KSP 2271 2272 Input Parameter: 2273 . pc - the preconditioner context 2274 2275 Output Parameters: 2276 + n - the number of splits 2277 - subksp - the array of KSP contexts 2278 2279 Note: 2280 After PCFieldSplitSchurGetSubKSP() the array of KSPs is to be freed by the user with PetscFree() 2281 (not the KSP just the array that contains them). 2282 2283 You must call PCSetUp() before calling PCFieldSplitSchurGetSubKSP(). 2284 2285 If the fieldsplit type is of type PC_COMPOSITE_SCHUR, it returns (in order) 2286 - the KSP used for the (1,1) block 2287 - the KSP used for the Schur complement (not the one used for the interior Schur solver) 2288 - the KSP used for the (1,1) block in the upper triangular factor (if different from that of the (1,1) block). 2289 2290 It returns a null array if the fieldsplit is not of type PC_COMPOSITE_SCHUR; in this case, you should use PCFieldSplitGetSubKSP(). 2291 2292 Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs. 2293 You can call PCFieldSplitSchurGetSubKSP(pc,n,PETSC_NULL_KSP,ierr) to determine how large the 2294 KSP array must be. 2295 2296 Level: advanced 2297 2298 .seealso: `PCFIELDSPLIT` 2299 @*/ 2300 PetscErrorCode PCFieldSplitSchurGetSubKSP(PC pc,PetscInt *n,KSP *subksp[]) 2301 { 2302 PetscFunctionBegin; 2303 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2304 if (n) PetscValidIntPointer(n,2); 2305 PetscUseMethod(pc,"PCFieldSplitSchurGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp)); 2306 PetscFunctionReturn(0); 2307 } 2308 2309 /*@ 2310 PCFieldSplitSetSchurPre - Indicates from what operator the preconditioner is constructucted for the Schur complement. 2311 The default is the A11 matrix. 2312 2313 Collective on PC 2314 2315 Input Parameters: 2316 + pc - the preconditioner context 2317 . 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 2318 PC_FIELDSPLIT_SCHUR_PRE_SELFP, and PC_FIELDSPLIT_SCHUR_PRE_FULL 2319 - userpre - matrix to use for preconditioning, or NULL 2320 2321 Options Database: 2322 + -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is a11. See notes for meaning of various arguments 2323 - -fieldsplit_1_pc_type <pctype> - the preconditioner algorithm that is used to construct the preconditioner from the operator 2324 2325 Notes: 2326 $ If ptype is 2327 $ a11 - the preconditioner for the Schur complement is generated from the block diagonal part of the preconditioner 2328 $ matrix associated with the Schur complement (i.e. A11), not the Schur complement matrix 2329 $ self - the preconditioner for the Schur complement is generated from the symbolic representation of the Schur complement matrix: 2330 $ The only preconditioner that currently works with this symbolic respresentation matrix object is the PCLSC 2331 $ preconditioner 2332 $ user - the preconditioner for the Schur complement is generated from the user provided matrix (pre argument 2333 $ to this function). 2334 $ selfp - the preconditioning for the Schur complement is generated from an explicitly-assembled approximation Sp = A11 - A10 inv(diag(A00)) A01 2335 $ This is only a good preconditioner when diag(A00) is a good preconditioner for A00. Optionally, A00 can be 2336 $ lumped before extracting the diagonal using the additional option -fieldsplit_1_mat_schur_complement_ainv_type lump 2337 $ full - the preconditioner for the Schur complement is generated from the exact Schur complement matrix representation computed internally by PCFIELDSPLIT (this is expensive) 2338 $ useful mostly as a test that the Schur complement approach can work for your problem 2339 2340 When solving a saddle point problem, where the A11 block is identically zero, using a11 as the ptype only makes sense 2341 with the additional option -fieldsplit_1_pc_type none. Usually for saddle point problems one would use a ptype of self and 2342 -fieldsplit_1_pc_type lsc which uses the least squares commutator to compute a preconditioner for the Schur complement. 2343 2344 Level: intermediate 2345 2346 .seealso: `PCFieldSplitGetSchurPre()`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurPreType`, 2347 `MatSchurComplementSetAinvType()`, `PCLSC` 2348 2349 @*/ 2350 PetscErrorCode PCFieldSplitSetSchurPre(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 2351 { 2352 PetscFunctionBegin; 2353 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2354 PetscTryMethod(pc,"PCFieldSplitSetSchurPre_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre)); 2355 PetscFunctionReturn(0); 2356 } 2357 2358 PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) {return PCFieldSplitSetSchurPre(pc,ptype,pre);} /* Deprecated name */ 2359 2360 /*@ 2361 PCFieldSplitGetSchurPre - For Schur complement fieldsplit, determine how the Schur complement will be 2362 preconditioned. See PCFieldSplitSetSchurPre() for details. 2363 2364 Logically Collective on PC 2365 2366 Input Parameter: 2367 . pc - the preconditioner context 2368 2369 Output Parameters: 2370 + ptype - which matrix to use for preconditioning the Schur complement: PC_FIELDSPLIT_SCHUR_PRE_A11, PC_FIELDSPLIT_SCHUR_PRE_SELF, PC_FIELDSPLIT_PRE_USER 2371 - userpre - matrix to use for preconditioning (with PC_FIELDSPLIT_PRE_USER), or NULL 2372 2373 Level: intermediate 2374 2375 .seealso: `PCFieldSplitSetSchurPre()`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurPreType`, `PCLSC` 2376 2377 @*/ 2378 PetscErrorCode PCFieldSplitGetSchurPre(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre) 2379 { 2380 PetscFunctionBegin; 2381 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2382 PetscUseMethod(pc,"PCFieldSplitGetSchurPre_C",(PC,PCFieldSplitSchurPreType*,Mat*),(pc,ptype,pre)); 2383 PetscFunctionReturn(0); 2384 } 2385 2386 /*@ 2387 PCFieldSplitSchurGetS - extract the MatSchurComplement object used by this PC in case it needs to be configured separately 2388 2389 Not collective 2390 2391 Input Parameter: 2392 . pc - the preconditioner context 2393 2394 Output Parameter: 2395 . S - the Schur complement matrix 2396 2397 Notes: 2398 This matrix should not be destroyed using MatDestroy(); rather, use PCFieldSplitSchurRestoreS(). 2399 2400 Level: advanced 2401 2402 .seealso: `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSchurPreType`, `PCFieldSplitSetSchurPre()`, `MatSchurComplement`, `PCFieldSplitSchurRestoreS()` 2403 2404 @*/ 2405 PetscErrorCode PCFieldSplitSchurGetS(PC pc,Mat *S) 2406 { 2407 const char* t; 2408 PetscBool isfs; 2409 PC_FieldSplit *jac; 2410 2411 PetscFunctionBegin; 2412 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2413 PetscCall(PetscObjectGetType((PetscObject)pc,&t)); 2414 PetscCall(PetscStrcmp(t,PCFIELDSPLIT,&isfs)); 2415 PetscCheck(isfs,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t); 2416 jac = (PC_FieldSplit*)pc->data; 2417 PetscCheck(jac->type == PC_COMPOSITE_SCHUR,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %d instead",jac->type); 2418 if (S) *S = jac->schur; 2419 PetscFunctionReturn(0); 2420 } 2421 2422 /*@ 2423 PCFieldSplitSchurRestoreS - restores the MatSchurComplement object used by this PC 2424 2425 Not collective 2426 2427 Input Parameters: 2428 + pc - the preconditioner context 2429 - S - the Schur complement matrix 2430 2431 Level: advanced 2432 2433 .seealso: `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSchurPreType`, `PCFieldSplitSetSchurPre()`, `MatSchurComplement`, `PCFieldSplitSchurGetS()` 2434 2435 @*/ 2436 PetscErrorCode PCFieldSplitSchurRestoreS(PC pc,Mat *S) 2437 { 2438 const char* t; 2439 PetscBool isfs; 2440 PC_FieldSplit *jac; 2441 2442 PetscFunctionBegin; 2443 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2444 PetscCall(PetscObjectGetType((PetscObject)pc,&t)); 2445 PetscCall(PetscStrcmp(t,PCFIELDSPLIT,&isfs)); 2446 PetscCheck(isfs,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t); 2447 jac = (PC_FieldSplit*)pc->data; 2448 PetscCheck(jac->type == PC_COMPOSITE_SCHUR,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %d instead",jac->type); 2449 PetscCheck(S && (*S == jac->schur),PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"MatSchurComplement restored is not the same as gotten"); 2450 PetscFunctionReturn(0); 2451 } 2452 2453 static PetscErrorCode PCFieldSplitSetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 2454 { 2455 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2456 2457 PetscFunctionBegin; 2458 jac->schurpre = ptype; 2459 if (ptype == PC_FIELDSPLIT_SCHUR_PRE_USER && pre) { 2460 PetscCall(MatDestroy(&jac->schur_user)); 2461 jac->schur_user = pre; 2462 PetscCall(PetscObjectReference((PetscObject)jac->schur_user)); 2463 } 2464 PetscFunctionReturn(0); 2465 } 2466 2467 static PetscErrorCode PCFieldSplitGetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre) 2468 { 2469 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2470 2471 PetscFunctionBegin; 2472 *ptype = jac->schurpre; 2473 *pre = jac->schur_user; 2474 PetscFunctionReturn(0); 2475 } 2476 2477 /*@ 2478 PCFieldSplitSetSchurFactType - sets which blocks of the approximate block factorization to retain in the preconditioner 2479 2480 Collective on PC 2481 2482 Input Parameters: 2483 + pc - the preconditioner context 2484 - ftype - which blocks of factorization to retain, PC_FIELDSPLIT_SCHUR_FACT_FULL is default 2485 2486 Options Database: 2487 . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> - default is full 2488 2489 Level: intermediate 2490 2491 Notes: 2492 The FULL factorization is 2493 2494 .vb 2495 (A B) = (1 0) (A 0) (1 Ainv*B) = L D U 2496 (C E) (C*Ainv 1) (0 S) (0 1) 2497 .vb 2498 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, 2499 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(). 2500 2501 If A and S are solved exactly 2502 .vb 2503 *) FULL factorization is a direct solver. 2504 *) 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. 2505 *) 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. 2506 .ve 2507 2508 If the iteration count is very low, consider using KSPFGMRES or KSPGCR which can use one less preconditioner 2509 application in this case. Note that the preconditioned operator may be highly non-normal, so such fast convergence may not be observed in practice. 2510 2511 For symmetric problems in which A is positive definite and S is negative definite, DIAG can be used with KSPMINRES. 2512 2513 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). 2514 2515 References: 2516 + * - Murphy, Golub, and Wathen, A note on preconditioning indefinite linear systems, SIAM J. Sci. Comput., 21 (2000). 2517 - * - Ipsen, A note on preconditioning nonsymmetric matrices, SIAM J. Sci. Comput., 23 (2001). 2518 2519 .seealso: `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurPreType`, `PCFieldSplitSetSchurScale()` 2520 @*/ 2521 PetscErrorCode PCFieldSplitSetSchurFactType(PC pc,PCFieldSplitSchurFactType ftype) 2522 { 2523 PetscFunctionBegin; 2524 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2525 PetscTryMethod(pc,"PCFieldSplitSetSchurFactType_C",(PC,PCFieldSplitSchurFactType),(pc,ftype)); 2526 PetscFunctionReturn(0); 2527 } 2528 2529 static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc,PCFieldSplitSchurFactType ftype) 2530 { 2531 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2532 2533 PetscFunctionBegin; 2534 jac->schurfactorization = ftype; 2535 PetscFunctionReturn(0); 2536 } 2537 2538 /*@ 2539 PCFieldSplitSetSchurScale - Controls the sign flip of S for PC_FIELDSPLIT_SCHUR_FACT_DIAG. 2540 2541 Collective on PC 2542 2543 Input Parameters: 2544 + pc - the preconditioner context 2545 - scale - scaling factor for the Schur complement 2546 2547 Options Database: 2548 . -pc_fieldsplit_schur_scale - default is -1.0 2549 2550 Level: intermediate 2551 2552 .seealso: `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurFactType`, `PCFieldSplitSetSchurScale()` 2553 @*/ 2554 PetscErrorCode PCFieldSplitSetSchurScale(PC pc,PetscScalar scale) 2555 { 2556 PetscFunctionBegin; 2557 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2558 PetscValidLogicalCollectiveScalar(pc,scale,2); 2559 PetscTryMethod(pc,"PCFieldSplitSetSchurScale_C",(PC,PetscScalar),(pc,scale)); 2560 PetscFunctionReturn(0); 2561 } 2562 2563 static PetscErrorCode PCFieldSplitSetSchurScale_FieldSplit(PC pc,PetscScalar scale) 2564 { 2565 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2566 2567 PetscFunctionBegin; 2568 jac->schurscale = scale; 2569 PetscFunctionReturn(0); 2570 } 2571 2572 /*@C 2573 PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement 2574 2575 Collective on KSP 2576 2577 Input Parameter: 2578 . pc - the preconditioner context 2579 2580 Output Parameters: 2581 + A00 - the (0,0) block 2582 . A01 - the (0,1) block 2583 . A10 - the (1,0) block 2584 - A11 - the (1,1) block 2585 2586 Level: advanced 2587 2588 .seealso: `PCFIELDSPLIT` 2589 @*/ 2590 PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11) 2591 { 2592 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 2593 2594 PetscFunctionBegin; 2595 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2596 PetscCheck(jac->type == PC_COMPOSITE_SCHUR,PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach."); 2597 if (A00) *A00 = jac->pmat[0]; 2598 if (A01) *A01 = jac->B; 2599 if (A10) *A10 = jac->C; 2600 if (A11) *A11 = jac->pmat[1]; 2601 PetscFunctionReturn(0); 2602 } 2603 2604 /*@ 2605 PCFieldSplitSetGKBTol - Sets the solver tolerance for the generalized Golub-Kahan bidiagonalization preconditioner. 2606 2607 Collective on PC 2608 2609 Notes: 2610 The generalized GKB algorithm uses a lower bound estimate of the error in energy norm as stopping criterion. 2611 It stops once the lower bound estimate undershoots the required solver tolerance. Although the actual error might be bigger than 2612 this estimate, the stopping criterion is satisfactory in practical cases [A13]. 2613 2614 [Ar13] Generalized Golub-Kahan bidiagonalization and stopping criteria, SIAM J. Matrix Anal. Appl., Vol. 34, No. 2, pp. 571-592, 2013. 2615 2616 Input Parameters: 2617 + pc - the preconditioner context 2618 - tolerance - the solver tolerance 2619 2620 Options Database: 2621 . -pc_fieldsplit_gkb_tol - default is 1e-5 2622 2623 Level: intermediate 2624 2625 .seealso: `PCFIELDSPLIT`, `PCFieldSplitSetGKBDelay()`, `PCFieldSplitSetGKBNu()`, `PCFieldSplitSetGKBMaxit()` 2626 @*/ 2627 PetscErrorCode PCFieldSplitSetGKBTol(PC pc,PetscReal tolerance) 2628 { 2629 PetscFunctionBegin; 2630 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2631 PetscValidLogicalCollectiveReal(pc,tolerance,2); 2632 PetscTryMethod(pc,"PCFieldSplitSetGKBTol_C",(PC,PetscReal),(pc,tolerance)); 2633 PetscFunctionReturn(0); 2634 } 2635 2636 static PetscErrorCode PCFieldSplitSetGKBTol_FieldSplit(PC pc,PetscReal tolerance) 2637 { 2638 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2639 2640 PetscFunctionBegin; 2641 jac->gkbtol = tolerance; 2642 PetscFunctionReturn(0); 2643 } 2644 2645 /*@ 2646 PCFieldSplitSetGKBMaxit - Sets the maximum number of iterations for the generalized Golub-Kahan bidiagonalization 2647 preconditioner. 2648 2649 Collective on PC 2650 2651 Input Parameters: 2652 + pc - the preconditioner context 2653 - maxit - the maximum number of iterations 2654 2655 Options Database: 2656 . -pc_fieldsplit_gkb_maxit - default is 100 2657 2658 Level: intermediate 2659 2660 .seealso: `PCFIELDSPLIT`, `PCFieldSplitSetGKBDelay()`, `PCFieldSplitSetGKBTol()`, `PCFieldSplitSetGKBNu()` 2661 @*/ 2662 PetscErrorCode PCFieldSplitSetGKBMaxit(PC pc,PetscInt maxit) 2663 { 2664 PetscFunctionBegin; 2665 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2666 PetscValidLogicalCollectiveInt(pc,maxit,2); 2667 PetscTryMethod(pc,"PCFieldSplitSetGKBMaxit_C",(PC,PetscInt),(pc,maxit)); 2668 PetscFunctionReturn(0); 2669 } 2670 2671 static PetscErrorCode PCFieldSplitSetGKBMaxit_FieldSplit(PC pc,PetscInt maxit) 2672 { 2673 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2674 2675 PetscFunctionBegin; 2676 jac->gkbmaxit = maxit; 2677 PetscFunctionReturn(0); 2678 } 2679 2680 /*@ 2681 PCFieldSplitSetGKBDelay - Sets the delay in the lower bound error estimate in the generalized Golub-Kahan bidiagonalization 2682 preconditioner. 2683 2684 Collective on PC 2685 2686 Notes: 2687 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 2688 is expressed as a truncated sum. The error at iteration k can only be measured at iteration (k + delay), and thus the algorithm needs 2689 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 2690 2691 [Ar13] Generalized Golub-Kahan bidiagonalization and stopping criteria, SIAM J. Matrix Anal. Appl., Vol. 34, No. 2, pp. 571-592, 2013. 2692 2693 Input Parameters: 2694 + pc - the preconditioner context 2695 - delay - the delay window in the lower bound estimate 2696 2697 Options Database: 2698 . -pc_fieldsplit_gkb_delay - default is 5 2699 2700 Level: intermediate 2701 2702 .seealso: `PCFIELDSPLIT`, `PCFieldSplitSetGKBNu()`, `PCFieldSplitSetGKBTol()`, `PCFieldSplitSetGKBMaxit()` 2703 @*/ 2704 PetscErrorCode PCFieldSplitSetGKBDelay(PC pc,PetscInt delay) 2705 { 2706 PetscFunctionBegin; 2707 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2708 PetscValidLogicalCollectiveInt(pc,delay,2); 2709 PetscTryMethod(pc,"PCFieldSplitSetGKBDelay_C",(PC,PetscInt),(pc,delay)); 2710 PetscFunctionReturn(0); 2711 } 2712 2713 static PetscErrorCode PCFieldSplitSetGKBDelay_FieldSplit(PC pc,PetscInt delay) 2714 { 2715 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2716 2717 PetscFunctionBegin; 2718 jac->gkbdelay = delay; 2719 PetscFunctionReturn(0); 2720 } 2721 2722 /*@ 2723 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. 2724 2725 Collective on PC 2726 2727 Notes: 2728 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, 2729 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 2730 necessary to find a good balance in between the convergence of the inner and outer loop. 2731 2732 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. 2733 2734 [Ar13] Generalized Golub-Kahan bidiagonalization and stopping criteria, SIAM J. Matrix Anal. Appl., Vol. 34, No. 2, pp. 571-592, 2013. 2735 2736 Input Parameters: 2737 + pc - the preconditioner context 2738 - nu - the shift parameter 2739 2740 Options Database: 2741 . -pc_fieldsplit_gkb_nu - default is 1 2742 2743 Level: intermediate 2744 2745 .seealso: `PCFIELDSPLIT`, `PCFieldSplitSetGKBDelay()`, `PCFieldSplitSetGKBTol()`, `PCFieldSplitSetGKBMaxit()` 2746 @*/ 2747 PetscErrorCode PCFieldSplitSetGKBNu(PC pc,PetscReal nu) 2748 { 2749 PetscFunctionBegin; 2750 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2751 PetscValidLogicalCollectiveReal(pc,nu,2); 2752 PetscTryMethod(pc,"PCFieldSplitSetGKBNu_C",(PC,PetscReal),(pc,nu)); 2753 PetscFunctionReturn(0); 2754 } 2755 2756 static PetscErrorCode PCFieldSplitSetGKBNu_FieldSplit(PC pc,PetscReal nu) 2757 { 2758 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2759 2760 PetscFunctionBegin; 2761 jac->gkbnu = nu; 2762 PetscFunctionReturn(0); 2763 } 2764 2765 static PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type) 2766 { 2767 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2768 2769 PetscFunctionBegin; 2770 jac->type = type; 2771 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",NULL)); 2772 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",NULL)); 2773 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",NULL)); 2774 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",NULL)); 2775 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurScale_C",NULL)); 2776 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBTol_C",NULL)); 2777 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBMaxit_C",NULL)); 2778 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBNu_C",NULL)); 2779 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBDelay_C",NULL)); 2780 2781 if (type == PC_COMPOSITE_SCHUR) { 2782 pc->ops->apply = PCApply_FieldSplit_Schur; 2783 pc->ops->view = PCView_FieldSplit_Schur; 2784 2785 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit_Schur)); 2786 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",PCFieldSplitSetSchurPre_FieldSplit)); 2787 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",PCFieldSplitGetSchurPre_FieldSplit)); 2788 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",PCFieldSplitSetSchurFactType_FieldSplit)); 2789 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurScale_C",PCFieldSplitSetSchurScale_FieldSplit)); 2790 } else if (type == PC_COMPOSITE_GKB) { 2791 pc->ops->apply = PCApply_FieldSplit_GKB; 2792 pc->ops->view = PCView_FieldSplit_GKB; 2793 2794 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit)); 2795 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBTol_C",PCFieldSplitSetGKBTol_FieldSplit)); 2796 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBMaxit_C",PCFieldSplitSetGKBMaxit_FieldSplit)); 2797 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBNu_C",PCFieldSplitSetGKBNu_FieldSplit)); 2798 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBDelay_C",PCFieldSplitSetGKBDelay_FieldSplit)); 2799 } else { 2800 pc->ops->apply = PCApply_FieldSplit; 2801 pc->ops->view = PCView_FieldSplit; 2802 2803 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit)); 2804 } 2805 PetscFunctionReturn(0); 2806 } 2807 2808 static PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs) 2809 { 2810 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2811 2812 PetscFunctionBegin; 2813 PetscCheck(bs >= 1,PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %" PetscInt_FMT,bs); 2814 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); 2815 jac->bs = bs; 2816 PetscFunctionReturn(0); 2817 } 2818 2819 static PetscErrorCode PCSetCoordinates_FieldSplit(PC pc, PetscInt dim, PetscInt nloc, PetscReal coords[]) 2820 { 2821 PC_FieldSplit * jac = (PC_FieldSplit*)pc->data; 2822 PC_FieldSplitLink ilink_current = jac->head; 2823 IS is_owned; 2824 2825 PetscFunctionBegin; 2826 jac->coordinates_set = PETSC_TRUE; // Internal flag 2827 PetscCall(MatGetOwnershipIS(pc->mat,&is_owned,PETSC_NULL)); 2828 2829 while (ilink_current) { 2830 // For each IS, embed it to get local coords indces 2831 IS is_coords; 2832 PetscInt ndofs_block; 2833 const PetscInt *block_dofs_enumeration; // Numbering of the dofs relevant to the current block 2834 2835 // Setting drop to true for safety. It should make no difference. 2836 PetscCall(ISEmbed(ilink_current->is, is_owned, PETSC_TRUE, &is_coords)); 2837 PetscCall(ISGetLocalSize(is_coords, &ndofs_block)); 2838 PetscCall(ISGetIndices(is_coords, &block_dofs_enumeration)); 2839 2840 // Allocate coordinates vector and set it directly 2841 PetscCall(PetscMalloc1(ndofs_block * dim, &(ilink_current->coords))); 2842 for (PetscInt dof=0;dof<ndofs_block;++dof) { 2843 for (PetscInt d=0;d<dim;++d) { 2844 (ilink_current->coords)[dim*dof + d] = coords[dim * block_dofs_enumeration[dof] + d]; 2845 } 2846 } 2847 ilink_current->dim = dim; 2848 ilink_current->ndofs = ndofs_block; 2849 PetscCall(ISRestoreIndices(is_coords, &block_dofs_enumeration)); 2850 PetscCall(ISDestroy(&is_coords)); 2851 ilink_current = ilink_current->next; 2852 } 2853 PetscCall(ISDestroy(&is_owned)); 2854 PetscFunctionReturn(0); 2855 } 2856 2857 /*@ 2858 PCFieldSplitSetType - Sets the type of fieldsplit preconditioner. 2859 2860 Collective on PC 2861 2862 Input Parameters: 2863 + pc - the preconditioner context 2864 - type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 2865 2866 Options Database Key: 2867 . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type 2868 2869 Level: Intermediate 2870 2871 .seealso: `PCCompositeSetType()` 2872 2873 @*/ 2874 PetscErrorCode PCFieldSplitSetType(PC pc,PCCompositeType type) 2875 { 2876 PetscFunctionBegin; 2877 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2878 PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type)); 2879 PetscFunctionReturn(0); 2880 } 2881 2882 /*@ 2883 PCFieldSplitGetType - Gets the type of fieldsplit preconditioner. 2884 2885 Not collective 2886 2887 Input Parameter: 2888 . pc - the preconditioner context 2889 2890 Output Parameter: 2891 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 2892 2893 Level: Intermediate 2894 2895 .seealso: `PCCompositeSetType()` 2896 @*/ 2897 PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type) 2898 { 2899 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 2900 2901 PetscFunctionBegin; 2902 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2903 PetscValidIntPointer(type,2); 2904 *type = jac->type; 2905 PetscFunctionReturn(0); 2906 } 2907 2908 /*@ 2909 PCFieldSplitSetDMSplits - Flags whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible. 2910 2911 Logically Collective 2912 2913 Input Parameters: 2914 + pc - the preconditioner context 2915 - flg - boolean indicating whether to use field splits defined by the DM 2916 2917 Options Database Key: 2918 . -pc_fieldsplit_dm_splits <bool> - use the field splits defined by the DM 2919 2920 Level: Intermediate 2921 2922 .seealso: `PCFieldSplitGetDMSplits()` 2923 2924 @*/ 2925 PetscErrorCode PCFieldSplitSetDMSplits(PC pc,PetscBool flg) 2926 { 2927 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2928 PetscBool isfs; 2929 2930 PetscFunctionBegin; 2931 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2932 PetscValidLogicalCollectiveBool(pc,flg,2); 2933 PetscCall(PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs)); 2934 if (isfs) { 2935 jac->dm_splits = flg; 2936 } 2937 PetscFunctionReturn(0); 2938 } 2939 2940 /*@ 2941 PCFieldSplitGetDMSplits - Returns flag indicating whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible. 2942 2943 Logically Collective 2944 2945 Input Parameter: 2946 . pc - the preconditioner context 2947 2948 Output Parameter: 2949 . flg - boolean indicating whether to use field splits defined by the DM 2950 2951 Level: Intermediate 2952 2953 .seealso: `PCFieldSplitSetDMSplits()` 2954 2955 @*/ 2956 PetscErrorCode PCFieldSplitGetDMSplits(PC pc,PetscBool* flg) 2957 { 2958 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2959 PetscBool isfs; 2960 2961 PetscFunctionBegin; 2962 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2963 PetscValidBoolPointer(flg,2); 2964 PetscCall(PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs)); 2965 if (isfs) { 2966 if (flg) *flg = jac->dm_splits; 2967 } 2968 PetscFunctionReturn(0); 2969 } 2970 2971 /*@ 2972 PCFieldSplitGetDetectSaddlePoint - Returns flag indicating whether PCFieldSplit will attempt to automatically determine fields based on zero diagonal entries. 2973 2974 Logically Collective 2975 2976 Input Parameter: 2977 . pc - the preconditioner context 2978 2979 Output Parameter: 2980 . flg - boolean indicating whether to detect fields or not 2981 2982 Level: Intermediate 2983 2984 .seealso: `PCFIELDSPLIT`, `PCFieldSplitSetDetectSaddlePoint()` 2985 2986 @*/ 2987 PetscErrorCode PCFieldSplitGetDetectSaddlePoint(PC pc,PetscBool *flg) 2988 { 2989 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2990 2991 PetscFunctionBegin; 2992 *flg = jac->detect; 2993 PetscFunctionReturn(0); 2994 } 2995 2996 /*@ 2997 PCFieldSplitSetDetectSaddlePoint - Sets flag indicating whether PCFieldSplit will attempt to automatically determine fields based on zero diagonal entries. 2998 2999 Logically Collective 3000 3001 Input Parameter: 3002 . pc - the preconditioner context 3003 3004 Output Parameter: 3005 . flg - boolean indicating whether to detect fields or not 3006 3007 Options Database Key: 3008 . -pc_fieldsplit_detect_saddle_point <bool> - detect and use the saddle point 3009 3010 Notes: 3011 Also sets the split type to PC_COMPOSITE_SCHUR (see PCFieldSplitSetType()) and the Schur preconditioner type to PC_FIELDSPLIT_SCHUR_PRE_SELF (see PCFieldSplitSetSchurPre()). 3012 3013 Level: Intermediate 3014 3015 .seealso: `PCFIELDSPLIT`, `PCFieldSplitSetDetectSaddlePoint()`, `PCFieldSplitSetType()`, `PCFieldSplitSetSchurPre()` 3016 3017 @*/ 3018 PetscErrorCode PCFieldSplitSetDetectSaddlePoint(PC pc,PetscBool flg) 3019 { 3020 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 3021 3022 PetscFunctionBegin; 3023 jac->detect = flg; 3024 if (jac->detect) { 3025 PetscCall(PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR)); 3026 PetscCall(PCFieldSplitSetSchurPre(pc,PC_FIELDSPLIT_SCHUR_PRE_SELF,NULL)); 3027 } 3028 PetscFunctionReturn(0); 3029 } 3030 3031 /* -------------------------------------------------------------------------------------*/ 3032 /*MC 3033 PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual 3034 fields or groups of fields. See [the users manual section on "Solving Block Matrices"](sec_block_matrices) for more details. 3035 3036 To set options on the solvers for each block append `-fieldsplit_` to all the PC 3037 options database keys. For example, `-fieldsplit_pc_type ilu -fieldsplit_pc_factor_levels 1` 3038 3039 To set the options on the solvers separate for each block call `PCFieldSplitGetSubKSP()` 3040 and set the options directly on the resulting `KSP` object 3041 3042 Level: intermediate 3043 3044 Options Database Keys: 3045 + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the `%d`'th split 3046 . -pc_fieldsplit_default - automatically add any fields to additional splits that have not 3047 been supplied explicitly by `-pc_fieldsplit_%d_fields` 3048 . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields) 3049 . -pc_fieldsplit_type <additive,multiplicative,symmetric_multiplicative,schur,gkb> - type of relaxation or factorization splitting 3050 . -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is a11; see `PCFieldSplitSetSchurPre()` 3051 . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> - set factorization type when using `-pc_fieldsplit_type schur`; see `PCFieldSplitSetSchurFactType()` 3052 - -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero diagonal and uses Schur complement with no preconditioner as the solver 3053 3054 Options prefixes for inner solvers when using the Schur complement preconditioner are `-fieldsplit_0_` and `-fieldsplit_1_` . 3055 For all other solvers they are `-fieldsplit_%d_` for the `d`th field; use `-fieldsplit_` for all fields. 3056 The options prefix for the inner solver when using the Golub-Kahan biadiagonalization preconditioner is `-fieldsplit_0_` 3057 3058 Notes: 3059 Use `PCFieldSplitSetFields()` to set fields defined by "strided" entries and `PCFieldSplitSetIS()` 3060 to define a field by an arbitrary collection of entries. 3061 3062 If no fields are set the default is used. The fields are defined by entries strided by bs, 3063 beginning at 0 then 1, etc to bs-1. The block size can be set with `PCFieldSplitSetBlockSize()`, 3064 if this is not called the block size defaults to the blocksize of the second matrix passed 3065 to `KSPSetOperators()`/`PCSetOperators()`. 3066 3067 For the Schur complement preconditioner if 3068 3069 ```{math} 3070 J = \left[\begin{array}{cc} A_{00} & A_{01} \\ A_{10} & A_{11} \end{array}\right] 3071 ``` 3072 3073 the preconditioner using `full` factorization is logically 3074 ```{math} 3075 \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] 3076 ``` 3077 where the action of $\text{inv}(A_{00})$ is applied using the KSP solver with prefix `-fieldsplit_0_`. $S$ is the Schur complement 3078 ```{math} 3079 S = A_{11} - A_{10} \text{ksp}(A_{00}) A_{01} 3080 ``` 3081 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 3082 in providing the SECOND split or 1 if not given). For `PCFieldSplitGetSub\text{ksp}()` when field number is 0, 3083 it returns the KSP associated with `-fieldsplit_0_` while field number 1 gives `-fieldsplit_1_` KSP. By default 3084 $A_{11}$ is used to construct a preconditioner for $S$, use `PCFieldSplitSetSchurPre()` for all the possible ways to construct the preconditioner for $S$. 3085 3086 The factorization type is set using `-pc_fieldsplit_schur_fact_type <diag, lower, upper, full>`. `full` is shown above, 3087 `diag` gives 3088 ```{math} 3089 \left[\begin{array}{cc} \text{inv}(A_{00}) & 0 \\ 0 & -\text{ksp}(S) \end{array}\right] 3090 ``` 3091 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 3092 can be turned off with `PCFieldSplitSetSchurScale()` or by command line `-pc_fieldsplit_schur_scale 1.0`. The `lower` factorization is the inverse of 3093 ```{math} 3094 \left[\begin{array}{cc} A_{00} & 0 \\ A_{10} & S \end{array}\right] 3095 ``` 3096 where the inverses of A_{00} and S are applied using KSPs. The upper factorization is the inverse of 3097 ```{math} 3098 \left[\begin{array}{cc} A_{00} & A_{01} \\ 0 & S \end{array}\right] 3099 ``` 3100 where again the inverses of $A_{00}$ and $S$ are applied using `KSP`s. 3101 3102 If only one set of indices (one `IS`) is provided with `PCFieldSplitSetIS()` then the complement of that `IS` 3103 is used automatically for a second block. 3104 3105 The fieldsplit preconditioner cannot currently be used with the BAIJ or SBAIJ data formats if the blocksize is larger than 1. 3106 Generally it should be used with the AIJ format. 3107 3108 The forms of these preconditioners are closely related if not identical to forms derived as "Distributive Iterations", see, 3109 for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling {cite}`Wesseling2009`. Note that one can also use `PCFIELDSPLIT` 3110 inside a smoother resulting in "Distributive Smoothers". 3111 3112 References: 3113 3114 See "A taxonomy and comparison of parallel block multi-level preconditioners for the incompressible Navier-Stokes equations" {cite}`elman2008tcp`. 3115 3116 The Constrained Pressure Preconditioner (CPR) can be implemented using `PCCOMPOSITE` with `PCGALERKIN`. CPR first solves an $R A P$ subsystem, updates the 3117 residual on all variables (`PCCompositeSetType(pc,PC_COMPOSITE_MULTIPLICATIVE)`), and then applies a simple ILU like preconditioner on all the variables. 3118 3119 The generalized Golub-Kahan bidiagonalization preconditioner (GKB) can be applied to symmetric $2 \times 2$ block matrices of the shape 3120 ```{math} 3121 \left[\begin{array}{cc} A_{00} & A_{01} \\ A_{01}' & 0 \end{array}\right] 3122 ``` 3123 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}'$. 3124 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_`. 3125 3126 ```{bibliography} 3127 :filter: docname in docnames 3128 ``` 3129 3130 .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCLSC`, 3131 `PCFieldSplitGetSubKSP()`, `PCFieldSplitSchurGetSubKSP()`, `PCFieldSplitSetFields()`, 3132 `PCFieldSplitSetType()`, `PCFieldSplitSetIS()`, `PCFieldSplitSetSchurPre()`, `PCFieldSplitSetSchurFactType()`, 3133 `MatSchurComplementSetAinvType()`, `PCFieldSplitSetSchurScale()`, `PCFieldSplitSetDetectSaddlePoint()` 3134 M*/ 3135 3136 PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc) 3137 { 3138 PC_FieldSplit *jac; 3139 3140 PetscFunctionBegin; 3141 PetscCall(PetscNewLog(pc,&jac)); 3142 3143 jac->bs = -1; 3144 jac->nsplits = 0; 3145 jac->type = PC_COMPOSITE_MULTIPLICATIVE; 3146 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */ 3147 jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL; 3148 jac->schurscale = -1.0; 3149 jac->dm_splits = PETSC_TRUE; 3150 jac->detect = PETSC_FALSE; 3151 jac->gkbtol = 1e-5; 3152 jac->gkbdelay = 5; 3153 jac->gkbnu = 1; 3154 jac->gkbmaxit = 100; 3155 jac->gkbmonitor = PETSC_FALSE; 3156 jac->coordinates_set = PETSC_FALSE; 3157 3158 pc->data = (void*)jac; 3159 3160 pc->ops->apply = PCApply_FieldSplit; 3161 pc->ops->applytranspose = PCApplyTranspose_FieldSplit; 3162 pc->ops->setup = PCSetUp_FieldSplit; 3163 pc->ops->reset = PCReset_FieldSplit; 3164 pc->ops->destroy = PCDestroy_FieldSplit; 3165 pc->ops->setfromoptions = PCSetFromOptions_FieldSplit; 3166 pc->ops->view = PCView_FieldSplit; 3167 pc->ops->applyrichardson = NULL; 3168 3169 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSchurGetSubKSP_C",PCFieldSplitSchurGetSubKSP_FieldSplit)); 3170 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit)); 3171 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",PCFieldSplitSetFields_FieldSplit)); 3172 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",PCFieldSplitSetIS_FieldSplit)); 3173 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",PCFieldSplitSetType_FieldSplit)); 3174 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",PCFieldSplitSetBlockSize_FieldSplit)); 3175 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitRestrictIS_C",PCFieldSplitRestrictIS_FieldSplit)); 3176 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",PCSetCoordinates_FieldSplit)); 3177 PetscFunctionReturn(0); 3178 } 3179