1 #include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/ 2 #include <petsc/private/kspimpl.h> /* This is needed to provide the appropriate PETSC_EXTERN for KSP_Solve_FS ....*/ 3 #include <petscdm.h> 4 5 const char *const PCFieldSplitSchurPreTypes[] = {"SELF","SELFP","A11","USER","FULL","PCFieldSplitSchurPreType","PC_FIELDSPLIT_SCHUR_PRE_",NULL}; 6 const char *const PCFieldSplitSchurFactTypes[] = {"DIAG","LOWER","UPPER","FULL","PCFieldSplitSchurFactType","PC_FIELDSPLIT_SCHUR_FACT_",NULL}; 7 8 PetscLogEvent KSP_Solve_FS_0,KSP_Solve_FS_1,KSP_Solve_FS_S,KSP_Solve_FS_U,KSP_Solve_FS_L,KSP_Solve_FS_2,KSP_Solve_FS_3,KSP_Solve_FS_4; 9 10 typedef struct _PC_FieldSplitLink *PC_FieldSplitLink; 11 struct _PC_FieldSplitLink { 12 KSP ksp; 13 Vec x,y,z; 14 char *splitname; 15 PetscInt nfields; 16 PetscInt *fields,*fields_col; 17 VecScatter sctx; 18 IS is,is_col; 19 PC_FieldSplitLink next,previous; 20 PetscLogEvent event; 21 22 /* Used only when setting coordinates with PCSetCoordinates */ 23 PetscInt dim; 24 PetscInt ndofs; 25 PetscReal *coords; 26 }; 27 28 typedef struct { 29 PCCompositeType type; 30 PetscBool defaultsplit; /* Flag for a system with a set of 'k' scalar fields with the same layout (and bs = k) */ 31 PetscBool splitdefined; /* Flag is set after the splits have been defined, to prevent more splits from being added */ 32 PetscInt bs; /* Block size for IS and Mat structures */ 33 PetscInt nsplits; /* Number of field divisions defined */ 34 Vec *x,*y,w1,w2; 35 Mat *mat; /* The diagonal block for each split */ 36 Mat *pmat; /* The preconditioning diagonal block for each split */ 37 Mat *Afield; /* The rows of the matrix associated with each split */ 38 PetscBool issetup; 39 40 /* Only used when Schur complement preconditioning is used */ 41 Mat B; /* The (0,1) block */ 42 Mat C; /* The (1,0) block */ 43 Mat schur; /* The Schur complement S = A11 - A10 A00^{-1} A01, the KSP here, kspinner, is H_1 in [El08] */ 44 Mat schurp; /* Assembled approximation to S built by MatSchurComplement to be used as a preconditioning matrix when solving with S */ 45 Mat schur_user; /* User-provided preconditioning matrix for the Schur complement */ 46 PCFieldSplitSchurPreType schurpre; /* Determines which preconditioning matrix is used for the Schur complement */ 47 PCFieldSplitSchurFactType schurfactorization; 48 KSP kspschur; /* The solver for S */ 49 KSP kspupper; /* The solver for A in the upper diagonal part of the factorization (H_2 in [El08]) */ 50 PetscScalar schurscale; /* Scaling factor for the Schur complement solution with DIAG factorization */ 51 52 /* Only used when Golub-Kahan bidiagonalization preconditioning is used */ 53 Mat H; /* The modified matrix H = A00 + nu*A01*A01' */ 54 PetscReal gkbtol; /* Stopping tolerance for lower bound estimate */ 55 PetscInt gkbdelay; /* The delay window for the stopping criterion */ 56 PetscReal gkbnu; /* Parameter for augmented Lagrangian H = A + nu*A01*A01' */ 57 PetscInt gkbmaxit; /* Maximum number of iterations for outer loop */ 58 PetscBool gkbmonitor; /* Monitor for gkb iterations and the lower bound error */ 59 PetscViewer gkbviewer; /* Viewer context for gkbmonitor */ 60 Vec u,v,d,Hu; /* Work vectors for the GKB algorithm */ 61 PetscScalar *vecz; /* Contains intermediate values, eg for lower bound */ 62 63 PC_FieldSplitLink head; 64 PetscBool isrestrict; /* indicates PCFieldSplitRestrictIS() has been last called on this object, hack */ 65 PetscBool suboptionsset; /* Indicates that the KSPSetFromOptions() has been called on the sub-KSPs */ 66 PetscBool dm_splits; /* Whether to use DMCreateFieldDecomposition() whenever possible */ 67 PetscBool diag_use_amat; /* Whether to extract diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */ 68 PetscBool offdiag_use_amat; /* Whether to extract off-diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */ 69 PetscBool detect; /* Whether to form 2-way split by finding zero diagonal entries */ 70 PetscBool coordinates_set; /* Whether PCSetCoordinates has been called */ 71 } PC_FieldSplit; 72 73 /* 74 Notes: 75 there is no particular reason that pmat, x, and y are stored as arrays in PC_FieldSplit instead of 76 inside PC_FieldSplitLink, just historical. If you want to be able to add new fields after already using the 77 PC you could change this. 78 */ 79 80 /* This helper is so that setting a user-provided preconditioning matrix is orthogonal to choosing to use it. This way the 81 * application-provided FormJacobian can provide this matrix without interfering with the user's (command-line) choices. */ 82 static Mat FieldSplitSchurPre(PC_FieldSplit *jac) 83 { 84 switch (jac->schurpre) { 85 case PC_FIELDSPLIT_SCHUR_PRE_SELF: return jac->schur; 86 case PC_FIELDSPLIT_SCHUR_PRE_SELFP: return jac->schurp; 87 case PC_FIELDSPLIT_SCHUR_PRE_A11: return jac->pmat[1]; 88 case PC_FIELDSPLIT_SCHUR_PRE_FULL: /* We calculate this and store it in schur_user */ 89 case PC_FIELDSPLIT_SCHUR_PRE_USER: /* Use a user-provided matrix if it is given, otherwise diagonal block */ 90 default: 91 return jac->schur_user ? jac->schur_user : jac->pmat[1]; 92 } 93 } 94 95 #include <petscdraw.h> 96 static PetscErrorCode PCView_FieldSplit(PC pc,PetscViewer viewer) 97 { 98 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 99 PetscBool iascii,isdraw; 100 PetscInt i,j; 101 PC_FieldSplitLink ilink = jac->head; 102 103 PetscFunctionBegin; 104 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 105 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw)); 106 if (iascii) { 107 if (jac->bs > 0) { 108 PetscCall(PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %" PetscInt_FMT ", blocksize = %" PetscInt_FMT "\n",PCCompositeTypes[jac->type],jac->nsplits,jac->bs)); 109 } else { 110 PetscCall(PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %" PetscInt_FMT "\n",PCCompositeTypes[jac->type],jac->nsplits)); 111 } 112 if (pc->useAmat) { 113 PetscCall(PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for blocks\n")); 114 } 115 if (jac->diag_use_amat) { 116 PetscCall(PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for diagonal blocks\n")); 117 } 118 if (jac->offdiag_use_amat) { 119 PetscCall(PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for off-diagonal blocks\n")); 120 } 121 PetscCall(PetscViewerASCIIPrintf(viewer," Solver info for each split is in the following KSP objects:\n")); 122 for (i=0; i<jac->nsplits; i++) { 123 if (ilink->fields) { 124 PetscCall(PetscViewerASCIIPrintf(viewer,"Split number %" PetscInt_FMT " Fields ",i)); 125 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 126 for (j=0; j<ilink->nfields; j++) { 127 if (j > 0) { 128 PetscCall(PetscViewerASCIIPrintf(viewer,",")); 129 } 130 PetscCall(PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT,ilink->fields[j])); 131 } 132 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 133 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 134 } else { 135 PetscCall(PetscViewerASCIIPrintf(viewer,"Split number %" PetscInt_FMT " Defined by IS\n",i)); 136 } 137 PetscCall(KSPView(ilink->ksp,viewer)); 138 ilink = ilink->next; 139 } 140 } 141 142 if (isdraw) { 143 PetscDraw draw; 144 PetscReal x,y,w,wd; 145 146 PetscCall(PetscViewerDrawGetDraw(viewer,0,&draw)); 147 PetscCall(PetscDrawGetCurrentPoint(draw,&x,&y)); 148 w = 2*PetscMin(1.0 - x,x); 149 wd = w/(jac->nsplits + 1); 150 x = x - wd*(jac->nsplits-1)/2.0; 151 for (i=0; i<jac->nsplits; i++) { 152 PetscCall(PetscDrawPushCurrentPoint(draw,x,y)); 153 PetscCall(KSPView(ilink->ksp,viewer)); 154 PetscCall(PetscDrawPopCurrentPoint(draw)); 155 x += wd; 156 ilink = ilink->next; 157 } 158 } 159 PetscFunctionReturn(0); 160 } 161 162 static PetscErrorCode PCView_FieldSplit_Schur(PC pc,PetscViewer viewer) 163 { 164 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 165 PetscBool iascii,isdraw; 166 PetscInt i,j; 167 PC_FieldSplitLink ilink = jac->head; 168 MatSchurComplementAinvType atype; 169 170 PetscFunctionBegin; 171 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 172 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw)); 173 if (iascii) { 174 if (jac->bs > 0) { 175 PetscCall(PetscViewerASCIIPrintf(viewer," FieldSplit with Schur preconditioner, blocksize = %" PetscInt_FMT ", factorization %s\n",jac->bs,PCFieldSplitSchurFactTypes[jac->schurfactorization])); 176 } else { 177 PetscCall(PetscViewerASCIIPrintf(viewer," FieldSplit with Schur preconditioner, factorization %s\n",PCFieldSplitSchurFactTypes[jac->schurfactorization])); 178 } 179 if (pc->useAmat) { 180 PetscCall(PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for blocks\n")); 181 } 182 switch (jac->schurpre) { 183 case PC_FIELDSPLIT_SCHUR_PRE_SELF: 184 PetscCall(PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from S itself\n")); 185 break; 186 case PC_FIELDSPLIT_SCHUR_PRE_SELFP: 187 PetscCall(MatSchurComplementGetAinvType(jac->schur,&atype)); 188 PetscCall(PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from Sp, an assembled approximation to S, which uses A00's %sdiagonal's inverse\n",atype == MAT_SCHUR_COMPLEMENT_AINV_DIAG ? "" : (atype == MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG ? "block " : "lumped ")));break; 189 case PC_FIELDSPLIT_SCHUR_PRE_A11: 190 PetscCall(PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from A11\n")); 191 break; 192 case PC_FIELDSPLIT_SCHUR_PRE_FULL: 193 PetscCall(PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from the exact Schur complement\n")); 194 break; 195 case PC_FIELDSPLIT_SCHUR_PRE_USER: 196 if (jac->schur_user) { 197 PetscCall(PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from user provided matrix\n")); 198 } else { 199 PetscCall(PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from A11\n")); 200 } 201 break; 202 default: 203 SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Invalid Schur preconditioning type: %d", jac->schurpre); 204 } 205 PetscCall(PetscViewerASCIIPrintf(viewer," Split info:\n")); 206 PetscCall(PetscViewerASCIIPushTab(viewer)); 207 for (i=0; i<jac->nsplits; i++) { 208 if (ilink->fields) { 209 PetscCall(PetscViewerASCIIPrintf(viewer,"Split number %" PetscInt_FMT " Fields ",i)); 210 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 211 for (j=0; j<ilink->nfields; j++) { 212 if (j > 0) { 213 PetscCall(PetscViewerASCIIPrintf(viewer,",")); 214 } 215 PetscCall(PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT,ilink->fields[j])); 216 } 217 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 218 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 219 } else { 220 PetscCall(PetscViewerASCIIPrintf(viewer,"Split number %" PetscInt_FMT " Defined by IS\n",i)); 221 } 222 ilink = ilink->next; 223 } 224 PetscCall(PetscViewerASCIIPrintf(viewer,"KSP solver for A00 block\n")); 225 PetscCall(PetscViewerASCIIPushTab(viewer)); 226 if (jac->head) { 227 PetscCall(KSPView(jac->head->ksp,viewer)); 228 } else PetscCall(PetscViewerASCIIPrintf(viewer," not yet available\n")); 229 PetscCall(PetscViewerASCIIPopTab(viewer)); 230 if (jac->head && jac->kspupper != jac->head->ksp) { 231 PetscCall(PetscViewerASCIIPrintf(viewer,"KSP solver for upper A00 in upper triangular factor \n")); 232 PetscCall(PetscViewerASCIIPushTab(viewer)); 233 if (jac->kspupper) PetscCall(KSPView(jac->kspupper,viewer)); 234 else PetscCall(PetscViewerASCIIPrintf(viewer," not yet available\n")); 235 PetscCall(PetscViewerASCIIPopTab(viewer)); 236 } 237 PetscCall(PetscViewerASCIIPrintf(viewer,"KSP solver for S = A11 - A10 inv(A00) A01 \n")); 238 PetscCall(PetscViewerASCIIPushTab(viewer)); 239 if (jac->kspschur) { 240 PetscCall(KSPView(jac->kspschur,viewer)); 241 } else { 242 PetscCall(PetscViewerASCIIPrintf(viewer," not yet available\n")); 243 } 244 PetscCall(PetscViewerASCIIPopTab(viewer)); 245 PetscCall(PetscViewerASCIIPopTab(viewer)); 246 } else if (isdraw && jac->head) { 247 PetscDraw draw; 248 PetscReal x,y,w,wd,h; 249 PetscInt cnt = 2; 250 char str[32]; 251 252 PetscCall(PetscViewerDrawGetDraw(viewer,0,&draw)); 253 PetscCall(PetscDrawGetCurrentPoint(draw,&x,&y)); 254 if (jac->kspupper != jac->head->ksp) cnt++; 255 w = 2*PetscMin(1.0 - x,x); 256 wd = w/(cnt + 1); 257 258 PetscCall(PetscSNPrintf(str,32,"Schur fact. %s",PCFieldSplitSchurFactTypes[jac->schurfactorization])); 259 PetscCall(PetscDrawStringBoxed(draw,x,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h)); 260 y -= h; 261 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_USER && !jac->schur_user) { 262 PetscCall(PetscSNPrintf(str,32,"Prec. for Schur from %s",PCFieldSplitSchurPreTypes[PC_FIELDSPLIT_SCHUR_PRE_A11])); 263 } else { 264 PetscCall(PetscSNPrintf(str,32,"Prec. for Schur from %s",PCFieldSplitSchurPreTypes[jac->schurpre])); 265 } 266 PetscCall(PetscDrawStringBoxed(draw,x+wd*(cnt-1)/2.0,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h)); 267 y -= h; 268 x = x - wd*(cnt-1)/2.0; 269 270 PetscCall(PetscDrawPushCurrentPoint(draw,x,y)); 271 PetscCall(KSPView(jac->head->ksp,viewer)); 272 PetscCall(PetscDrawPopCurrentPoint(draw)); 273 if (jac->kspupper != jac->head->ksp) { 274 x += wd; 275 PetscCall(PetscDrawPushCurrentPoint(draw,x,y)); 276 PetscCall(KSPView(jac->kspupper,viewer)); 277 PetscCall(PetscDrawPopCurrentPoint(draw)); 278 } 279 x += wd; 280 PetscCall(PetscDrawPushCurrentPoint(draw,x,y)); 281 PetscCall(KSPView(jac->kspschur,viewer)); 282 PetscCall(PetscDrawPopCurrentPoint(draw)); 283 } 284 PetscFunctionReturn(0); 285 } 286 287 static PetscErrorCode PCView_FieldSplit_GKB(PC pc,PetscViewer viewer) 288 { 289 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 290 PetscBool iascii,isdraw; 291 PetscInt i,j; 292 PC_FieldSplitLink ilink = jac->head; 293 294 PetscFunctionBegin; 295 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 296 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw)); 297 if (iascii) { 298 if (jac->bs > 0) { 299 PetscCall(PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %" PetscInt_FMT ", blocksize = %" PetscInt_FMT "\n",PCCompositeTypes[jac->type],jac->nsplits,jac->bs)); 300 } else { 301 PetscCall(PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %" PetscInt_FMT "\n",PCCompositeTypes[jac->type],jac->nsplits)); 302 } 303 if (pc->useAmat) { 304 PetscCall(PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for blocks\n")); 305 } 306 if (jac->diag_use_amat) { 307 PetscCall(PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for diagonal blocks\n")); 308 } 309 if (jac->offdiag_use_amat) { 310 PetscCall(PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for off-diagonal blocks\n")); 311 } 312 313 PetscCall(PetscViewerASCIIPrintf(viewer," Stopping tolerance=%.1e, delay in error estimate=%" PetscInt_FMT ", maximum iterations=%" PetscInt_FMT "\n",(double)jac->gkbtol,jac->gkbdelay,jac->gkbmaxit)); 314 PetscCall(PetscViewerASCIIPrintf(viewer," Solver info for H = A00 + nu*A01*A01' matrix:\n")); 315 PetscCall(PetscViewerASCIIPushTab(viewer)); 316 317 if (ilink->fields) { 318 PetscCall(PetscViewerASCIIPrintf(viewer,"Split number 0 Fields ")); 319 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 320 for (j=0; j<ilink->nfields; j++) { 321 if (j > 0) { 322 PetscCall(PetscViewerASCIIPrintf(viewer,",")); 323 } 324 PetscCall(PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT,ilink->fields[j])); 325 } 326 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 327 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 328 } else { 329 PetscCall(PetscViewerASCIIPrintf(viewer,"Split number 0 Defined by IS\n")); 330 } 331 PetscCall(KSPView(ilink->ksp,viewer)); 332 333 PetscCall(PetscViewerASCIIPopTab(viewer)); 334 } 335 336 if (isdraw) { 337 PetscDraw draw; 338 PetscReal x,y,w,wd; 339 340 PetscCall(PetscViewerDrawGetDraw(viewer,0,&draw)); 341 PetscCall(PetscDrawGetCurrentPoint(draw,&x,&y)); 342 w = 2*PetscMin(1.0 - x,x); 343 wd = w/(jac->nsplits + 1); 344 x = x - wd*(jac->nsplits-1)/2.0; 345 for (i=0; i<jac->nsplits; i++) { 346 PetscCall(PetscDrawPushCurrentPoint(draw,x,y)); 347 PetscCall(KSPView(ilink->ksp,viewer)); 348 PetscCall(PetscDrawPopCurrentPoint(draw)); 349 x += wd; 350 ilink = ilink->next; 351 } 352 } 353 PetscFunctionReturn(0); 354 } 355 356 /* Precondition: jac->bs is set to a meaningful value */ 357 static PetscErrorCode PCFieldSplitSetRuntimeSplits_Private(PC pc) 358 { 359 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 360 PetscInt i,nfields,*ifields,nfields_col,*ifields_col; 361 PetscBool flg,flg_col; 362 char optionname[128],splitname[8],optionname_col[128]; 363 364 PetscFunctionBegin; 365 PetscCall(PetscMalloc1(jac->bs,&ifields)); 366 PetscCall(PetscMalloc1(jac->bs,&ifields_col)); 367 for (i=0,flg=PETSC_TRUE;; i++) { 368 PetscCall(PetscSNPrintf(splitname,sizeof(splitname),"%" PetscInt_FMT,i)); 369 PetscCall(PetscSNPrintf(optionname,sizeof(optionname),"-pc_fieldsplit_%" PetscInt_FMT "_fields",i)); 370 PetscCall(PetscSNPrintf(optionname_col,sizeof(optionname_col),"-pc_fieldsplit_%" PetscInt_FMT "_fields_col",i)); 371 nfields = jac->bs; 372 nfields_col = jac->bs; 373 PetscCall(PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix,optionname,ifields,&nfields,&flg)); 374 PetscCall(PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix,optionname_col,ifields_col,&nfields_col,&flg_col)); 375 if (!flg) break; 376 else if (flg && !flg_col) { 377 PetscCheck(nfields,PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot list zero fields"); 378 PetscCall(PCFieldSplitSetFields(pc,splitname,nfields,ifields,ifields)); 379 } else { 380 PetscCheck(nfields && nfields_col,PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot list zero fields"); 381 PetscCheck(nfields == nfields_col,PETSC_COMM_SELF,PETSC_ERR_USER,"Number of row and column fields must match"); 382 PetscCall(PCFieldSplitSetFields(pc,splitname,nfields,ifields,ifields_col)); 383 } 384 } 385 if (i > 0) { 386 /* Makes command-line setting of splits take precedence over setting them in code. 387 Otherwise subsequent calls to PCFieldSplitSetIS() or PCFieldSplitSetFields() would 388 create new splits, which would probably not be what the user wanted. */ 389 jac->splitdefined = PETSC_TRUE; 390 } 391 PetscCall(PetscFree(ifields)); 392 PetscCall(PetscFree(ifields_col)); 393 PetscFunctionReturn(0); 394 } 395 396 static PetscErrorCode PCFieldSplitSetDefaults(PC pc) 397 { 398 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 399 PC_FieldSplitLink ilink = jac->head; 400 PetscBool fieldsplit_default = PETSC_FALSE,coupling = PETSC_FALSE; 401 PetscInt i; 402 403 PetscFunctionBegin; 404 /* 405 Kinda messy, but at least this now uses DMCreateFieldDecomposition(). 406 Should probably be rewritten. 407 */ 408 if (!ilink) { 409 PetscCall(PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_coupling",&coupling,NULL)); 410 if (pc->dm && jac->dm_splits && !jac->detect && !coupling) { 411 PetscInt numFields, f, i, j; 412 char **fieldNames; 413 IS *fields; 414 DM *dms; 415 DM subdm[128]; 416 PetscBool flg; 417 418 PetscCall(DMCreateFieldDecomposition(pc->dm, &numFields, &fieldNames, &fields, &dms)); 419 /* Allow the user to prescribe the splits */ 420 for (i = 0, flg = PETSC_TRUE;; i++) { 421 PetscInt ifields[128]; 422 IS compField; 423 char optionname[128], splitname[8]; 424 PetscInt nfields = numFields; 425 426 PetscCall(PetscSNPrintf(optionname, sizeof(optionname), "-pc_fieldsplit_%" PetscInt_FMT "_fields", i)); 427 PetscCall(PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix, optionname, ifields, &nfields, &flg)); 428 if (!flg) break; 429 PetscCheck(numFields <= 128,PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot currently support %" PetscInt_FMT " > 128 fields", numFields); 430 PetscCall(DMCreateSubDM(pc->dm, nfields, ifields, &compField, &subdm[i])); 431 if (nfields == 1) { 432 PetscCall(PCFieldSplitSetIS(pc, fieldNames[ifields[0]], compField)); 433 } else { 434 PetscCall(PetscSNPrintf(splitname, sizeof(splitname), "%" PetscInt_FMT, i)); 435 PetscCall(PCFieldSplitSetIS(pc, splitname, compField)); 436 } 437 PetscCall(ISDestroy(&compField)); 438 for (j = 0; j < nfields; ++j) { 439 f = ifields[j]; 440 PetscCall(PetscFree(fieldNames[f])); 441 PetscCall(ISDestroy(&fields[f])); 442 } 443 } 444 if (i == 0) { 445 for (f = 0; f < numFields; ++f) { 446 PetscCall(PCFieldSplitSetIS(pc, fieldNames[f], fields[f])); 447 PetscCall(PetscFree(fieldNames[f])); 448 PetscCall(ISDestroy(&fields[f])); 449 } 450 } else { 451 for (j=0; j<numFields; j++) { 452 PetscCall(DMDestroy(dms+j)); 453 } 454 PetscCall(PetscFree(dms)); 455 PetscCall(PetscMalloc1(i, &dms)); 456 for (j = 0; j < i; ++j) dms[j] = subdm[j]; 457 } 458 PetscCall(PetscFree(fieldNames)); 459 PetscCall(PetscFree(fields)); 460 if (dms) { 461 PetscCall(PetscInfo(pc, "Setting up physics based fieldsplit preconditioner using the embedded DM\n")); 462 for (ilink = jac->head, i = 0; ilink; ilink = ilink->next, ++i) { 463 const char *prefix; 464 PetscCall(PetscObjectGetOptionsPrefix((PetscObject)(ilink->ksp),&prefix)); 465 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)(dms[i]), prefix)); 466 PetscCall(KSPSetDM(ilink->ksp, dms[i])); 467 PetscCall(KSPSetDMActive(ilink->ksp, PETSC_FALSE)); 468 { 469 PetscErrorCode (*func)(KSP,Mat,Mat,void*); 470 void *ctx; 471 472 PetscCall(DMKSPGetComputeOperators(pc->dm, &func, &ctx)); 473 PetscCall(DMKSPSetComputeOperators(dms[i], func, ctx)); 474 } 475 PetscCall(PetscObjectIncrementTabLevel((PetscObject)dms[i],(PetscObject)ilink->ksp,0)); 476 PetscCall(DMDestroy(&dms[i])); 477 } 478 PetscCall(PetscFree(dms)); 479 } 480 } else { 481 if (jac->bs <= 0) { 482 if (pc->pmat) { 483 PetscCall(MatGetBlockSize(pc->pmat,&jac->bs)); 484 } else jac->bs = 1; 485 } 486 487 if (jac->detect) { 488 IS zerodiags,rest; 489 PetscInt nmin,nmax; 490 491 PetscCall(MatGetOwnershipRange(pc->mat,&nmin,&nmax)); 492 if (jac->diag_use_amat) { 493 PetscCall(MatFindZeroDiagonals(pc->mat,&zerodiags)); 494 } else { 495 PetscCall(MatFindZeroDiagonals(pc->pmat,&zerodiags)); 496 } 497 PetscCall(ISComplement(zerodiags,nmin,nmax,&rest)); 498 PetscCall(PCFieldSplitSetIS(pc,"0",rest)); 499 PetscCall(PCFieldSplitSetIS(pc,"1",zerodiags)); 500 PetscCall(ISDestroy(&zerodiags)); 501 PetscCall(ISDestroy(&rest)); 502 } else if (coupling) { 503 IS coupling,rest; 504 PetscInt nmin,nmax; 505 506 PetscCall(MatGetOwnershipRange(pc->mat,&nmin,&nmax)); 507 if (jac->offdiag_use_amat) { 508 PetscCall(MatFindOffBlockDiagonalEntries(pc->mat,&coupling)); 509 } else { 510 PetscCall(MatFindOffBlockDiagonalEntries(pc->pmat,&coupling)); 511 } 512 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)pc->mat),nmax-nmin,nmin,1,&rest)); 513 PetscCall(ISSetIdentity(rest)); 514 PetscCall(PCFieldSplitSetIS(pc,"0",rest)); 515 PetscCall(PCFieldSplitSetIS(pc,"1",coupling)); 516 PetscCall(ISDestroy(&coupling)); 517 PetscCall(ISDestroy(&rest)); 518 } else { 519 PetscCall(PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_default",&fieldsplit_default,NULL)); 520 if (!fieldsplit_default) { 521 /* Allow user to set fields from command line, if bs was known at the time of PCSetFromOptions_FieldSplit() 522 then it is set there. This is not ideal because we should only have options set in XXSetFromOptions(). */ 523 PetscCall(PCFieldSplitSetRuntimeSplits_Private(pc)); 524 if (jac->splitdefined) PetscCall(PetscInfo(pc,"Splits defined using the options database\n")); 525 } 526 if ((fieldsplit_default || !jac->splitdefined) && !jac->isrestrict) { 527 Mat M = pc->pmat; 528 PetscBool isnest; 529 530 PetscCall(PetscInfo(pc,"Using default splitting of fields\n")); 531 PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat,MATNEST,&isnest)); 532 if (!isnest) { 533 M = pc->mat; 534 PetscCall(PetscObjectTypeCompare((PetscObject)pc->mat,MATNEST,&isnest)); 535 } 536 if (isnest) { 537 IS *fields; 538 PetscInt nf; 539 540 PetscCall(MatNestGetSize(M,&nf,NULL)); 541 PetscCall(PetscMalloc1(nf,&fields)); 542 PetscCall(MatNestGetISs(M,fields,NULL)); 543 for (i=0;i<nf;i++) { 544 PetscCall(PCFieldSplitSetIS(pc,NULL,fields[i])); 545 } 546 PetscCall(PetscFree(fields)); 547 } else { 548 for (i=0; i<jac->bs; i++) { 549 char splitname[8]; 550 PetscCall(PetscSNPrintf(splitname,sizeof(splitname),"%" PetscInt_FMT,i)); 551 PetscCall(PCFieldSplitSetFields(pc,splitname,1,&i,&i)); 552 } 553 jac->defaultsplit = PETSC_TRUE; 554 } 555 } 556 } 557 } 558 } else if (jac->nsplits == 1) { 559 if (ilink->is) { 560 IS is2; 561 PetscInt nmin,nmax; 562 563 PetscCall(MatGetOwnershipRange(pc->mat,&nmin,&nmax)); 564 PetscCall(ISComplement(ilink->is,nmin,nmax,&is2)); 565 PetscCall(PCFieldSplitSetIS(pc,"1",is2)); 566 PetscCall(ISDestroy(&is2)); 567 } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Must provide at least two sets of fields to PCFieldSplit()"); 568 } 569 570 PetscCheck(jac->nsplits >= 2,PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Unhandled case, must have at least two fields, not %" PetscInt_FMT, jac->nsplits); 571 PetscFunctionReturn(0); 572 } 573 574 static PetscErrorCode MatGolubKahanComputeExplicitOperator(Mat A,Mat B,Mat C,Mat *H,PetscReal gkbnu) 575 { 576 Mat BT,T; 577 PetscReal nrmT,nrmB; 578 579 PetscFunctionBegin; 580 PetscCall(MatHermitianTranspose(C,MAT_INITIAL_MATRIX,&T)); /* Test if augmented matrix is symmetric */ 581 PetscCall(MatAXPY(T,-1.0,B,DIFFERENT_NONZERO_PATTERN)); 582 PetscCall(MatNorm(T,NORM_1,&nrmT)); 583 PetscCall(MatNorm(B,NORM_1,&nrmB)); 584 if (nrmB > 0) { 585 if (nrmT/nrmB >= PETSC_SMALL) { 586 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Matrix is not symmetric/hermitian, GKB is not applicable."); 587 } 588 } 589 /* Compute augmented Lagrangian matrix H = A00 + nu*A01*A01'. This corresponds to */ 590 /* setting N := 1/nu*I in [Ar13]. */ 591 PetscCall(MatHermitianTranspose(B,MAT_INITIAL_MATRIX,&BT)); 592 PetscCall(MatMatMult(B,BT,MAT_INITIAL_MATRIX,PETSC_DEFAULT,H)); /* H = A01*A01' */ 593 PetscCall(MatAYPX(*H,gkbnu,A,DIFFERENT_NONZERO_PATTERN)); /* H = A00 + nu*A01*A01' */ 594 595 PetscCall(MatDestroy(&BT)); 596 PetscCall(MatDestroy(&T)); 597 PetscFunctionReturn(0); 598 } 599 600 PETSC_EXTERN PetscErrorCode PetscOptionsFindPairPrefix_Private(PetscOptions,const char pre[], const char name[],const char *value[],PetscBool *flg); 601 602 static PetscErrorCode PCSetUp_FieldSplit(PC pc) 603 { 604 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 605 PC_FieldSplitLink ilink; 606 PetscInt i,nsplit; 607 PetscBool sorted, sorted_col; 608 609 PetscFunctionBegin; 610 pc->failedreason = PC_NOERROR; 611 PetscCall(PCFieldSplitSetDefaults(pc)); 612 nsplit = jac->nsplits; 613 ilink = jac->head; 614 615 /* get the matrices for each split */ 616 if (!jac->issetup) { 617 PetscInt rstart,rend,nslots,bs; 618 619 jac->issetup = PETSC_TRUE; 620 621 /* This is done here instead of in PCFieldSplitSetFields() because may not have matrix at that point */ 622 if (jac->defaultsplit || !ilink->is) { 623 if (jac->bs <= 0) jac->bs = nsplit; 624 } 625 bs = jac->bs; 626 PetscCall(MatGetOwnershipRange(pc->pmat,&rstart,&rend)); 627 nslots = (rend - rstart)/bs; 628 for (i=0; i<nsplit; i++) { 629 if (jac->defaultsplit) { 630 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+i,nsplit,&ilink->is)); 631 PetscCall(ISDuplicate(ilink->is,&ilink->is_col)); 632 } else if (!ilink->is) { 633 if (ilink->nfields > 1) { 634 PetscInt *ii,*jj,j,k,nfields = ilink->nfields,*fields = ilink->fields,*fields_col = ilink->fields_col; 635 PetscCall(PetscMalloc1(ilink->nfields*nslots,&ii)); 636 PetscCall(PetscMalloc1(ilink->nfields*nslots,&jj)); 637 for (j=0; j<nslots; j++) { 638 for (k=0; k<nfields; k++) { 639 ii[nfields*j + k] = rstart + bs*j + fields[k]; 640 jj[nfields*j + k] = rstart + bs*j + fields_col[k]; 641 } 642 } 643 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,ii,PETSC_OWN_POINTER,&ilink->is)); 644 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,jj,PETSC_OWN_POINTER,&ilink->is_col)); 645 PetscCall(ISSetBlockSize(ilink->is, nfields)); 646 PetscCall(ISSetBlockSize(ilink->is_col, nfields)); 647 } else { 648 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields[0],bs,&ilink->is)); 649 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields_col[0],bs,&ilink->is_col)); 650 } 651 } 652 PetscCall(ISSorted(ilink->is,&sorted)); 653 if (ilink->is_col) PetscCall(ISSorted(ilink->is_col,&sorted_col)); 654 PetscCheck(sorted && sorted_col,PETSC_COMM_SELF,PETSC_ERR_USER,"Fields must be sorted when creating split"); 655 ilink = ilink->next; 656 } 657 } 658 659 ilink = jac->head; 660 if (!jac->pmat) { 661 Vec xtmp; 662 663 PetscCall(MatCreateVecs(pc->pmat,&xtmp,NULL)); 664 PetscCall(PetscMalloc1(nsplit,&jac->pmat)); 665 PetscCall(PetscMalloc2(nsplit,&jac->x,nsplit,&jac->y)); 666 for (i=0; i<nsplit; i++) { 667 MatNullSpace sp; 668 669 /* Check for preconditioning matrix attached to IS */ 670 PetscCall(PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &jac->pmat[i])); 671 if (jac->pmat[i]) { 672 PetscCall(PetscObjectReference((PetscObject) jac->pmat[i])); 673 if (jac->type == PC_COMPOSITE_SCHUR) { 674 jac->schur_user = jac->pmat[i]; 675 676 PetscCall(PetscObjectReference((PetscObject) jac->schur_user)); 677 } 678 } else { 679 const char *prefix; 680 PetscCall(MatCreateSubMatrix(pc->pmat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->pmat[i])); 681 PetscCall(KSPGetOptionsPrefix(ilink->ksp,&prefix)); 682 PetscCall(MatSetOptionsPrefix(jac->pmat[i],prefix)); 683 PetscCall(MatViewFromOptions(jac->pmat[i],NULL,"-mat_view")); 684 } 685 /* create work vectors for each split */ 686 PetscCall(MatCreateVecs(jac->pmat[i],&jac->x[i],&jac->y[i])); 687 ilink->x = jac->x[i]; ilink->y = jac->y[i]; ilink->z = NULL; 688 /* compute scatter contexts needed by multiplicative versions and non-default splits */ 689 PetscCall(VecScatterCreate(xtmp,ilink->is,jac->x[i],NULL,&ilink->sctx)); 690 PetscCall(PetscObjectQuery((PetscObject) ilink->is, "nearnullspace", (PetscObject*) &sp)); 691 if (sp) { 692 PetscCall(MatSetNearNullSpace(jac->pmat[i], sp)); 693 } 694 ilink = ilink->next; 695 } 696 PetscCall(VecDestroy(&xtmp)); 697 } else { 698 MatReuse scall; 699 if (pc->flag == DIFFERENT_NONZERO_PATTERN) { 700 for (i=0; i<nsplit; i++) { 701 PetscCall(MatDestroy(&jac->pmat[i])); 702 } 703 scall = MAT_INITIAL_MATRIX; 704 } else scall = MAT_REUSE_MATRIX; 705 706 for (i=0; i<nsplit; i++) { 707 Mat pmat; 708 709 /* Check for preconditioning matrix attached to IS */ 710 PetscCall(PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &pmat)); 711 if (!pmat) { 712 PetscCall(MatCreateSubMatrix(pc->pmat,ilink->is,ilink->is_col,scall,&jac->pmat[i])); 713 } 714 ilink = ilink->next; 715 } 716 } 717 if (jac->diag_use_amat) { 718 ilink = jac->head; 719 if (!jac->mat) { 720 PetscCall(PetscMalloc1(nsplit,&jac->mat)); 721 for (i=0; i<nsplit; i++) { 722 PetscCall(MatCreateSubMatrix(pc->mat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->mat[i])); 723 ilink = ilink->next; 724 } 725 } else { 726 MatReuse scall; 727 if (pc->flag == DIFFERENT_NONZERO_PATTERN) { 728 for (i=0; i<nsplit; i++) { 729 PetscCall(MatDestroy(&jac->mat[i])); 730 } 731 scall = MAT_INITIAL_MATRIX; 732 } else scall = MAT_REUSE_MATRIX; 733 734 for (i=0; i<nsplit; i++) { 735 PetscCall(MatCreateSubMatrix(pc->mat,ilink->is,ilink->is_col,scall,&jac->mat[i])); 736 ilink = ilink->next; 737 } 738 } 739 } else { 740 jac->mat = jac->pmat; 741 } 742 743 /* Check for null space attached to IS */ 744 ilink = jac->head; 745 for (i=0; i<nsplit; i++) { 746 MatNullSpace sp; 747 748 PetscCall(PetscObjectQuery((PetscObject) ilink->is, "nullspace", (PetscObject*) &sp)); 749 if (sp) { 750 PetscCall(MatSetNullSpace(jac->mat[i], sp)); 751 } 752 ilink = ilink->next; 753 } 754 755 if (jac->type != PC_COMPOSITE_ADDITIVE && jac->type != PC_COMPOSITE_SCHUR && jac->type != PC_COMPOSITE_GKB) { 756 /* extract the rows of the matrix associated with each field: used for efficient computation of residual inside algorithm */ 757 /* FIXME: Can/should we reuse jac->mat whenever (jac->diag_use_amat) is true? */ 758 ilink = jac->head; 759 if (nsplit == 2 && jac->type == PC_COMPOSITE_MULTIPLICATIVE) { 760 /* special case need where Afield[0] is not needed and only certain columns of Afield[1] are needed since update is only on those rows of the solution */ 761 if (!jac->Afield) { 762 PetscCall(PetscCalloc1(nsplit,&jac->Afield)); 763 if (jac->offdiag_use_amat) { 764 PetscCall(MatCreateSubMatrix(pc->mat,ilink->next->is,ilink->is,MAT_INITIAL_MATRIX,&jac->Afield[1])); 765 } else { 766 PetscCall(MatCreateSubMatrix(pc->pmat,ilink->next->is,ilink->is,MAT_INITIAL_MATRIX,&jac->Afield[1])); 767 } 768 } else { 769 MatReuse scall; 770 771 if (pc->flag == DIFFERENT_NONZERO_PATTERN) { 772 PetscCall(MatDestroy(&jac->Afield[1])); 773 scall = MAT_INITIAL_MATRIX; 774 } else scall = MAT_REUSE_MATRIX; 775 776 if (jac->offdiag_use_amat) { 777 PetscCall(MatCreateSubMatrix(pc->mat,ilink->next->is,ilink->is,scall,&jac->Afield[1])); 778 } else { 779 PetscCall(MatCreateSubMatrix(pc->pmat,ilink->next->is,ilink->is,scall,&jac->Afield[1])); 780 } 781 } 782 } else { 783 if (!jac->Afield) { 784 PetscCall(PetscMalloc1(nsplit,&jac->Afield)); 785 for (i=0; i<nsplit; i++) { 786 if (jac->offdiag_use_amat) { 787 PetscCall(MatCreateSubMatrix(pc->mat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i])); 788 } else { 789 PetscCall(MatCreateSubMatrix(pc->pmat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i])); 790 } 791 ilink = ilink->next; 792 } 793 } else { 794 MatReuse scall; 795 if (pc->flag == DIFFERENT_NONZERO_PATTERN) { 796 for (i=0; i<nsplit; i++) { 797 PetscCall(MatDestroy(&jac->Afield[i])); 798 } 799 scall = MAT_INITIAL_MATRIX; 800 } else scall = MAT_REUSE_MATRIX; 801 802 for (i=0; i<nsplit; i++) { 803 if (jac->offdiag_use_amat) { 804 PetscCall(MatCreateSubMatrix(pc->mat,ilink->is,NULL,scall,&jac->Afield[i])); 805 } else { 806 PetscCall(MatCreateSubMatrix(pc->pmat,ilink->is,NULL,scall,&jac->Afield[i])); 807 } 808 ilink = ilink->next; 809 } 810 } 811 } 812 } 813 814 if (jac->type == PC_COMPOSITE_SCHUR) { 815 IS ccis; 816 PetscBool isspd; 817 PetscInt rstart,rend; 818 char lscname[256]; 819 PetscObject LSC_L; 820 821 PetscCheck(nsplit == 2,PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"To use Schur complement preconditioner you must have exactly 2 fields"); 822 823 /* If pc->mat is SPD, don't scale by -1 the Schur complement */ 824 if (jac->schurscale == (PetscScalar)-1.0) { 825 PetscCall(MatGetOption(pc->pmat,MAT_SPD,&isspd)); 826 jac->schurscale = (isspd == PETSC_TRUE) ? 1.0 : -1.0; 827 } 828 829 /* When extracting off-diagonal submatrices, we take complements from this range */ 830 PetscCall(MatGetOwnershipRangeColumn(pc->mat,&rstart,&rend)); 831 832 if (jac->schur) { 833 KSP kspA = jac->head->ksp, kspInner = NULL, kspUpper = jac->kspupper; 834 MatReuse scall; 835 836 if (pc->flag == DIFFERENT_NONZERO_PATTERN) { 837 scall = MAT_INITIAL_MATRIX; 838 PetscCall(MatDestroy(&jac->B)); 839 PetscCall(MatDestroy(&jac->C)); 840 } else scall = MAT_REUSE_MATRIX; 841 842 PetscCall(MatSchurComplementGetKSP(jac->schur, &kspInner)); 843 ilink = jac->head; 844 PetscCall(ISComplement(ilink->is_col,rstart,rend,&ccis)); 845 if (jac->offdiag_use_amat) { 846 PetscCall(MatCreateSubMatrix(pc->mat,ilink->is,ccis,scall,&jac->B)); 847 } else { 848 PetscCall(MatCreateSubMatrix(pc->pmat,ilink->is,ccis,scall,&jac->B)); 849 } 850 PetscCall(ISDestroy(&ccis)); 851 ilink = ilink->next; 852 PetscCall(ISComplement(ilink->is_col,rstart,rend,&ccis)); 853 if (jac->offdiag_use_amat) { 854 PetscCall(MatCreateSubMatrix(pc->mat,ilink->is,ccis,scall,&jac->C)); 855 } else { 856 PetscCall(MatCreateSubMatrix(pc->pmat,ilink->is,ccis,scall,&jac->C)); 857 } 858 PetscCall(ISDestroy(&ccis)); 859 PetscCall(MatSchurComplementUpdateSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1])); 860 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) { 861 PetscCall(MatDestroy(&jac->schurp)); 862 PetscCall(MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp)); 863 } 864 if (kspA != kspInner) { 865 PetscCall(KSPSetOperators(kspA,jac->mat[0],jac->pmat[0])); 866 } 867 if (kspUpper != kspA) { 868 PetscCall(KSPSetOperators(kspUpper,jac->mat[0],jac->pmat[0])); 869 } 870 PetscCall(KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac))); 871 } else { 872 const char *Dprefix; 873 char schurprefix[256], schurmatprefix[256]; 874 char schurtestoption[256]; 875 MatNullSpace sp; 876 PetscBool flg; 877 KSP kspt; 878 879 /* extract the A01 and A10 matrices */ 880 ilink = jac->head; 881 PetscCall(ISComplement(ilink->is_col,rstart,rend,&ccis)); 882 if (jac->offdiag_use_amat) { 883 PetscCall(MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B)); 884 } else { 885 PetscCall(MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B)); 886 } 887 PetscCall(ISDestroy(&ccis)); 888 ilink = ilink->next; 889 PetscCall(ISComplement(ilink->is_col,rstart,rend,&ccis)); 890 if (jac->offdiag_use_amat) { 891 PetscCall(MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C)); 892 } else { 893 PetscCall(MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C)); 894 } 895 PetscCall(ISDestroy(&ccis)); 896 897 /* Use mat[0] (diagonal block of Amat) preconditioned by pmat[0] to define Schur complement */ 898 PetscCall(MatCreate(((PetscObject)jac->mat[0])->comm,&jac->schur)); 899 PetscCall(MatSetType(jac->schur,MATSCHURCOMPLEMENT)); 900 PetscCall(MatSchurComplementSetSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1])); 901 PetscCall(PetscSNPrintf(schurmatprefix, sizeof(schurmatprefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname)); 902 PetscCall(MatSetOptionsPrefix(jac->schur,schurmatprefix)); 903 PetscCall(MatSchurComplementGetKSP(jac->schur,&kspt)); 904 PetscCall(KSPSetOptionsPrefix(kspt,schurmatprefix)); 905 906 /* Note: this is not true in general */ 907 PetscCall(MatGetNullSpace(jac->mat[1], &sp)); 908 if (sp) 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 PetscOptionsHeadEnd(); 1665 PetscFunctionReturn(0); 1666 } 1667 1668 /*------------------------------------------------------------------------------------*/ 1669 1670 static PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 1671 { 1672 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1673 PC_FieldSplitLink ilink,next = jac->head; 1674 char prefix[128]; 1675 PetscInt i; 1676 1677 PetscFunctionBegin; 1678 if (jac->splitdefined) { 1679 PetscCall(PetscInfo(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname)); 1680 PetscFunctionReturn(0); 1681 } 1682 for (i=0; i<n; i++) { 1683 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); 1684 PetscCheck(fields[i] >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %" PetscInt_FMT " requested",fields[i]); 1685 } 1686 PetscCall(PetscNew(&ilink)); 1687 if (splitname) { 1688 PetscCall(PetscStrallocpy(splitname,&ilink->splitname)); 1689 } else { 1690 PetscCall(PetscMalloc1(3,&ilink->splitname)); 1691 PetscCall(PetscSNPrintf(ilink->splitname,2,"%" PetscInt_FMT,jac->nsplits)); 1692 } 1693 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 */ 1694 PetscCall(PetscMalloc1(n,&ilink->fields)); 1695 PetscCall(PetscArraycpy(ilink->fields,fields,n)); 1696 PetscCall(PetscMalloc1(n,&ilink->fields_col)); 1697 PetscCall(PetscArraycpy(ilink->fields_col,fields_col,n)); 1698 1699 ilink->nfields = n; 1700 ilink->next = NULL; 1701 PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp)); 1702 PetscCall(KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure)); 1703 PetscCall(PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1)); 1704 PetscCall(KSPSetType(ilink->ksp,KSPPREONLY)); 1705 PetscCall(PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp)); 1706 1707 PetscCall(PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname)); 1708 PetscCall(KSPSetOptionsPrefix(ilink->ksp,prefix)); 1709 1710 if (!next) { 1711 jac->head = ilink; 1712 ilink->previous = NULL; 1713 } else { 1714 while (next->next) { 1715 next = next->next; 1716 } 1717 next->next = ilink; 1718 ilink->previous = next; 1719 } 1720 jac->nsplits++; 1721 PetscFunctionReturn(0); 1722 } 1723 1724 static PetscErrorCode PCFieldSplitSchurGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp) 1725 { 1726 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1727 1728 PetscFunctionBegin; 1729 *subksp = NULL; 1730 if (n) *n = 0; 1731 if (jac->type == PC_COMPOSITE_SCHUR) { 1732 PetscInt nn; 1733 1734 PetscCheck(jac->schur,PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call KSPSetUp() or PCSetUp() before calling PCFieldSplitSchurGetSubKSP()"); 1735 PetscCheck(jac->nsplits == 2,PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Unexpected number of splits %" PetscInt_FMT " != 2",jac->nsplits); 1736 nn = jac->nsplits + (jac->kspupper != jac->head->ksp ? 1 : 0); 1737 PetscCall(PetscMalloc1(nn,subksp)); 1738 (*subksp)[0] = jac->head->ksp; 1739 (*subksp)[1] = jac->kspschur; 1740 if (jac->kspupper != jac->head->ksp) (*subksp)[2] = jac->kspupper; 1741 if (n) *n = nn; 1742 } 1743 PetscFunctionReturn(0); 1744 } 1745 1746 static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp) 1747 { 1748 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1749 1750 PetscFunctionBegin; 1751 PetscCheck(jac->schur,PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call KSPSetUp() or PCSetUp() before calling PCFieldSplitGetSubKSP()"); 1752 PetscCall(PetscMalloc1(jac->nsplits,subksp)); 1753 PetscCall(MatSchurComplementGetKSP(jac->schur,*subksp)); 1754 1755 (*subksp)[1] = jac->kspschur; 1756 if (n) *n = jac->nsplits; 1757 PetscFunctionReturn(0); 1758 } 1759 1760 static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp) 1761 { 1762 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1763 PetscInt cnt = 0; 1764 PC_FieldSplitLink ilink = jac->head; 1765 1766 PetscFunctionBegin; 1767 PetscCall(PetscMalloc1(jac->nsplits,subksp)); 1768 while (ilink) { 1769 (*subksp)[cnt++] = ilink->ksp; 1770 ilink = ilink->next; 1771 } 1772 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); 1773 if (n) *n = jac->nsplits; 1774 PetscFunctionReturn(0); 1775 } 1776 1777 /*@C 1778 PCFieldSplitRestrictIS - Restricts the fieldsplit ISs to be within a given IS. 1779 1780 Input Parameters: 1781 + pc - the preconditioner context 1782 - is - the index set that defines the indices to which the fieldsplit is to be restricted 1783 1784 Level: advanced 1785 1786 @*/ 1787 PetscErrorCode PCFieldSplitRestrictIS(PC pc,IS isy) 1788 { 1789 PetscFunctionBegin; 1790 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1791 PetscValidHeaderSpecific(isy,IS_CLASSID,2); 1792 PetscTryMethod(pc,"PCFieldSplitRestrictIS_C",(PC,IS),(pc,isy)); 1793 PetscFunctionReturn(0); 1794 } 1795 1796 static PetscErrorCode PCFieldSplitRestrictIS_FieldSplit(PC pc, IS isy) 1797 { 1798 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1799 PC_FieldSplitLink ilink = jac->head, next; 1800 PetscInt localsize,size,sizez,i; 1801 const PetscInt *ind, *indz; 1802 PetscInt *indc, *indcz; 1803 PetscBool flg; 1804 1805 PetscFunctionBegin; 1806 PetscCall(ISGetLocalSize(isy,&localsize)); 1807 PetscCallMPI(MPI_Scan(&localsize,&size,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)isy))); 1808 size -= localsize; 1809 while (ilink) { 1810 IS isrl,isr; 1811 PC subpc; 1812 PetscCall(ISEmbed(ilink->is, isy, PETSC_TRUE, &isrl)); 1813 PetscCall(ISGetLocalSize(isrl,&localsize)); 1814 PetscCall(PetscMalloc1(localsize,&indc)); 1815 PetscCall(ISGetIndices(isrl,&ind)); 1816 PetscCall(PetscArraycpy(indc,ind,localsize)); 1817 PetscCall(ISRestoreIndices(isrl,&ind)); 1818 PetscCall(ISDestroy(&isrl)); 1819 for (i=0; i<localsize; i++) *(indc+i) += size; 1820 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)isy),localsize,indc,PETSC_OWN_POINTER,&isr)); 1821 PetscCall(PetscObjectReference((PetscObject)isr)); 1822 PetscCall(ISDestroy(&ilink->is)); 1823 ilink->is = isr; 1824 PetscCall(PetscObjectReference((PetscObject)isr)); 1825 PetscCall(ISDestroy(&ilink->is_col)); 1826 ilink->is_col = isr; 1827 PetscCall(ISDestroy(&isr)); 1828 PetscCall(KSPGetPC(ilink->ksp, &subpc)); 1829 PetscCall(PetscObjectTypeCompare((PetscObject)subpc,PCFIELDSPLIT,&flg)); 1830 if (flg) { 1831 IS iszl,isz; 1832 MPI_Comm comm; 1833 PetscCall(ISGetLocalSize(ilink->is,&localsize)); 1834 comm = PetscObjectComm((PetscObject)ilink->is); 1835 PetscCall(ISEmbed(isy, ilink->is, PETSC_TRUE, &iszl)); 1836 PetscCallMPI(MPI_Scan(&localsize,&sizez,1,MPIU_INT,MPI_SUM,comm)); 1837 sizez -= localsize; 1838 PetscCall(ISGetLocalSize(iszl,&localsize)); 1839 PetscCall(PetscMalloc1(localsize,&indcz)); 1840 PetscCall(ISGetIndices(iszl,&indz)); 1841 PetscCall(PetscArraycpy(indcz,indz,localsize)); 1842 PetscCall(ISRestoreIndices(iszl,&indz)); 1843 PetscCall(ISDestroy(&iszl)); 1844 for (i=0; i<localsize; i++) *(indcz+i) += sizez; 1845 PetscCall(ISCreateGeneral(comm,localsize,indcz,PETSC_OWN_POINTER,&isz)); 1846 PetscCall(PCFieldSplitRestrictIS(subpc,isz)); 1847 PetscCall(ISDestroy(&isz)); 1848 } 1849 next = ilink->next; 1850 ilink = next; 1851 } 1852 jac->isrestrict = PETSC_TRUE; 1853 PetscFunctionReturn(0); 1854 } 1855 1856 static PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is) 1857 { 1858 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1859 PC_FieldSplitLink ilink, next = jac->head; 1860 char prefix[128]; 1861 1862 PetscFunctionBegin; 1863 if (jac->splitdefined) { 1864 PetscCall(PetscInfo(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname)); 1865 PetscFunctionReturn(0); 1866 } 1867 PetscCall(PetscNew(&ilink)); 1868 if (splitname) { 1869 PetscCall(PetscStrallocpy(splitname,&ilink->splitname)); 1870 } else { 1871 PetscCall(PetscMalloc1(8,&ilink->splitname)); 1872 PetscCall(PetscSNPrintf(ilink->splitname,7,"%" PetscInt_FMT,jac->nsplits)); 1873 } 1874 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 */ 1875 PetscCall(PetscObjectReference((PetscObject)is)); 1876 PetscCall(ISDestroy(&ilink->is)); 1877 ilink->is = is; 1878 PetscCall(PetscObjectReference((PetscObject)is)); 1879 PetscCall(ISDestroy(&ilink->is_col)); 1880 ilink->is_col = is; 1881 ilink->next = NULL; 1882 PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp)); 1883 PetscCall(KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure)); 1884 PetscCall(PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1)); 1885 PetscCall(KSPSetType(ilink->ksp,KSPPREONLY)); 1886 PetscCall(PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp)); 1887 1888 PetscCall(PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname)); 1889 PetscCall(KSPSetOptionsPrefix(ilink->ksp,prefix)); 1890 1891 if (!next) { 1892 jac->head = ilink; 1893 ilink->previous = NULL; 1894 } else { 1895 while (next->next) { 1896 next = next->next; 1897 } 1898 next->next = ilink; 1899 ilink->previous = next; 1900 } 1901 jac->nsplits++; 1902 PetscFunctionReturn(0); 1903 } 1904 1905 /*@C 1906 PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner 1907 1908 Logically Collective on PC 1909 1910 Input Parameters: 1911 + pc - the preconditioner context 1912 . splitname - name of this split, if NULL the number of the split is used 1913 . n - the number of fields in this split 1914 - fields - the fields in this split 1915 1916 Level: intermediate 1917 1918 Notes: 1919 Use PCFieldSplitSetIS() to set a completely general set of indices as a field. 1920 1921 The PCFieldSplitSetFields() is for defining fields as strided blocks. For example, if the block 1922 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 1923 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x.... 1924 where the numbered entries indicate what is in the field. 1925 1926 This function is called once per split (it creates a new split each time). Solve options 1927 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1928 1929 Developer Note: This routine does not actually create the IS representing the split, that is delayed 1930 until PCSetUp_FieldSplit(), because information about the vector/matrix layouts may not be 1931 available when this routine is called. 1932 1933 .seealso: `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetBlockSize()`, `PCFieldSplitSetIS()` 1934 1935 @*/ 1936 PetscErrorCode PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 1937 { 1938 PetscFunctionBegin; 1939 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1940 PetscValidCharPointer(splitname,2); 1941 PetscCheck(n >= 1,PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Provided number of fields %" PetscInt_FMT " in split \"%s\" not positive",n,splitname); 1942 PetscValidIntPointer(fields,4); 1943 PetscTryMethod(pc,"PCFieldSplitSetFields_C",(PC,const char[],PetscInt,const PetscInt*,const PetscInt*),(pc,splitname,n,fields,fields_col)); 1944 PetscFunctionReturn(0); 1945 } 1946 1947 /*@ 1948 PCFieldSplitSetDiagUseAmat - set flag indicating whether to extract diagonal blocks from Amat (rather than Pmat) 1949 1950 Logically Collective on PC 1951 1952 Input Parameters: 1953 + pc - the preconditioner object 1954 - flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from 1955 1956 Options Database: 1957 . -pc_fieldsplit_diag_use_amat - use the Amat to provide the diagonal blocks 1958 1959 Level: intermediate 1960 1961 .seealso: `PCFieldSplitGetDiagUseAmat()`, `PCFieldSplitSetOffDiagUseAmat()`, `PCFIELDSPLIT` 1962 1963 @*/ 1964 PetscErrorCode PCFieldSplitSetDiagUseAmat(PC pc,PetscBool flg) 1965 { 1966 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1967 PetscBool isfs; 1968 1969 PetscFunctionBegin; 1970 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1971 PetscCall(PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs)); 1972 PetscCheck(isfs,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 1973 jac->diag_use_amat = flg; 1974 PetscFunctionReturn(0); 1975 } 1976 1977 /*@ 1978 PCFieldSplitGetDiagUseAmat - get the flag indicating whether to extract diagonal blocks from Amat (rather than Pmat) 1979 1980 Logically Collective on PC 1981 1982 Input Parameters: 1983 . pc - the preconditioner object 1984 1985 Output Parameters: 1986 . flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from 1987 1988 Level: intermediate 1989 1990 .seealso: `PCFieldSplitSetDiagUseAmat()`, `PCFieldSplitGetOffDiagUseAmat()`, `PCFIELDSPLIT` 1991 1992 @*/ 1993 PetscErrorCode PCFieldSplitGetDiagUseAmat(PC pc,PetscBool *flg) 1994 { 1995 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1996 PetscBool isfs; 1997 1998 PetscFunctionBegin; 1999 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2000 PetscValidBoolPointer(flg,2); 2001 PetscCall(PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs)); 2002 PetscCheck(isfs,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 2003 *flg = jac->diag_use_amat; 2004 PetscFunctionReturn(0); 2005 } 2006 2007 /*@ 2008 PCFieldSplitSetOffDiagUseAmat - set flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat) 2009 2010 Logically Collective on PC 2011 2012 Input Parameters: 2013 + pc - the preconditioner object 2014 - flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from 2015 2016 Options Database: 2017 . -pc_fieldsplit_off_diag_use_amat <bool> - use the Amat to extract the off-diagonal blocks 2018 2019 Level: intermediate 2020 2021 .seealso: `PCFieldSplitGetOffDiagUseAmat()`, `PCFieldSplitSetDiagUseAmat()`, `PCFIELDSPLIT` 2022 2023 @*/ 2024 PetscErrorCode PCFieldSplitSetOffDiagUseAmat(PC pc,PetscBool flg) 2025 { 2026 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2027 PetscBool isfs; 2028 2029 PetscFunctionBegin; 2030 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2031 PetscCall(PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs)); 2032 PetscCheck(isfs,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 2033 jac->offdiag_use_amat = flg; 2034 PetscFunctionReturn(0); 2035 } 2036 2037 /*@ 2038 PCFieldSplitGetOffDiagUseAmat - get the flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat) 2039 2040 Logically Collective on PC 2041 2042 Input Parameters: 2043 . pc - the preconditioner object 2044 2045 Output Parameters: 2046 . flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from 2047 2048 Level: intermediate 2049 2050 .seealso: `PCFieldSplitSetOffDiagUseAmat()`, `PCFieldSplitGetDiagUseAmat()`, `PCFIELDSPLIT` 2051 2052 @*/ 2053 PetscErrorCode PCFieldSplitGetOffDiagUseAmat(PC pc,PetscBool *flg) 2054 { 2055 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2056 PetscBool isfs; 2057 2058 PetscFunctionBegin; 2059 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2060 PetscValidBoolPointer(flg,2); 2061 PetscCall(PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs)); 2062 PetscCheck(isfs,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 2063 *flg = jac->offdiag_use_amat; 2064 PetscFunctionReturn(0); 2065 } 2066 2067 /*@C 2068 PCFieldSplitSetIS - Sets the exact elements for field 2069 2070 Logically Collective on PC 2071 2072 Input Parameters: 2073 + pc - the preconditioner context 2074 . splitname - name of this split, if NULL the number of the split is used 2075 - is - the index set that defines the vector elements in this field 2076 2077 Notes: 2078 Use PCFieldSplitSetFields(), for fields defined by strided types. 2079 2080 This function is called once per split (it creates a new split each time). Solve options 2081 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 2082 2083 Level: intermediate 2084 2085 .seealso: `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetBlockSize()` 2086 2087 @*/ 2088 PetscErrorCode PCFieldSplitSetIS(PC pc,const char splitname[],IS is) 2089 { 2090 PetscFunctionBegin; 2091 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2092 if (splitname) PetscValidCharPointer(splitname,2); 2093 PetscValidHeaderSpecific(is,IS_CLASSID,3); 2094 PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is)); 2095 PetscFunctionReturn(0); 2096 } 2097 2098 /*@C 2099 PCFieldSplitGetIS - Retrieves the elements for a field as an IS 2100 2101 Logically Collective on PC 2102 2103 Input Parameters: 2104 + pc - the preconditioner context 2105 - splitname - name of this split 2106 2107 Output Parameter: 2108 - is - the index set that defines the vector elements in this field, or NULL if the field is not found 2109 2110 Level: intermediate 2111 2112 .seealso: `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetIS()` 2113 2114 @*/ 2115 PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is) 2116 { 2117 PetscFunctionBegin; 2118 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2119 PetscValidCharPointer(splitname,2); 2120 PetscValidPointer(is,3); 2121 { 2122 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 2123 PC_FieldSplitLink ilink = jac->head; 2124 PetscBool found; 2125 2126 *is = NULL; 2127 while (ilink) { 2128 PetscCall(PetscStrcmp(ilink->splitname, splitname, &found)); 2129 if (found) { 2130 *is = ilink->is; 2131 break; 2132 } 2133 ilink = ilink->next; 2134 } 2135 } 2136 PetscFunctionReturn(0); 2137 } 2138 2139 /*@C 2140 PCFieldSplitGetISByIndex - Retrieves the elements for a given index field as an IS 2141 2142 Logically Collective on PC 2143 2144 Input Parameters: 2145 + pc - the preconditioner context 2146 - index - index of this split 2147 2148 Output Parameter: 2149 - is - the index set that defines the vector elements in this field 2150 2151 Level: intermediate 2152 2153 .seealso: `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitGetIS()`, `PCFieldSplitSetIS()` 2154 2155 @*/ 2156 PetscErrorCode PCFieldSplitGetISByIndex(PC pc,PetscInt index,IS *is) 2157 { 2158 PetscFunctionBegin; 2159 PetscCheck(index >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %" PetscInt_FMT " requested",index); 2160 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2161 PetscValidPointer(is,3); 2162 { 2163 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 2164 PC_FieldSplitLink ilink = jac->head; 2165 PetscInt i = 0; 2166 PetscCheck(index < jac->nsplits,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %" PetscInt_FMT " requested but only %" PetscInt_FMT " exist",index,jac->nsplits); 2167 2168 while (i < index) { 2169 ilink = ilink->next; 2170 ++i; 2171 } 2172 PetscCall(PCFieldSplitGetIS(pc,ilink->splitname,is)); 2173 } 2174 PetscFunctionReturn(0); 2175 } 2176 2177 /*@ 2178 PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the 2179 fieldsplit preconditioner. If not set the matrix block size is used. 2180 2181 Logically Collective on PC 2182 2183 Input Parameters: 2184 + pc - the preconditioner context 2185 - bs - the block size 2186 2187 Level: intermediate 2188 2189 .seealso: `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()` 2190 2191 @*/ 2192 PetscErrorCode PCFieldSplitSetBlockSize(PC pc,PetscInt bs) 2193 { 2194 PetscFunctionBegin; 2195 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2196 PetscValidLogicalCollectiveInt(pc,bs,2); 2197 PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs)); 2198 PetscFunctionReturn(0); 2199 } 2200 2201 /*@C 2202 PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits 2203 2204 Collective on KSP 2205 2206 Input Parameter: 2207 . pc - the preconditioner context 2208 2209 Output Parameters: 2210 + n - the number of splits 2211 - subksp - the array of KSP contexts 2212 2213 Note: 2214 After PCFieldSplitGetSubKSP() the array of KSPs is to be freed by the user with PetscFree() 2215 (not the KSP just the array that contains them). 2216 2217 You must call PCSetUp() before calling PCFieldSplitGetSubKSP(). 2218 2219 If the fieldsplit is of type PC_COMPOSITE_SCHUR, it returns the KSP object used inside the 2220 Schur complement and the KSP object used to iterate over the Schur complement. 2221 To access all the KSP objects used in PC_COMPOSITE_SCHUR, use PCFieldSplitSchurGetSubKSP(). 2222 2223 If the fieldsplit is of type PC_COMPOSITE_GKB, it returns the KSP object used to solve the 2224 inner linear system defined by the matrix H in each loop. 2225 2226 Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs. 2227 You can call PCFieldSplitGetSubKSP(pc,n,PETSC_NULL_KSP,ierr) to determine how large the 2228 KSP array must be. 2229 2230 Level: advanced 2231 2232 .seealso: `PCFIELDSPLIT` 2233 @*/ 2234 PetscErrorCode PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[]) 2235 { 2236 PetscFunctionBegin; 2237 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2238 if (n) PetscValidIntPointer(n,2); 2239 PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp)); 2240 PetscFunctionReturn(0); 2241 } 2242 2243 /*@C 2244 PCFieldSplitSchurGetSubKSP - Gets the KSP contexts used inside the Schur complement based PCFIELDSPLIT 2245 2246 Collective on KSP 2247 2248 Input Parameter: 2249 . pc - the preconditioner context 2250 2251 Output Parameters: 2252 + n - the number of splits 2253 - subksp - the array of KSP contexts 2254 2255 Note: 2256 After PCFieldSplitSchurGetSubKSP() the array of KSPs is to be freed by the user with PetscFree() 2257 (not the KSP just the array that contains them). 2258 2259 You must call PCSetUp() before calling PCFieldSplitSchurGetSubKSP(). 2260 2261 If the fieldsplit type is of type PC_COMPOSITE_SCHUR, it returns (in order) 2262 - the KSP used for the (1,1) block 2263 - the KSP used for the Schur complement (not the one used for the interior Schur solver) 2264 - the KSP used for the (1,1) block in the upper triangular factor (if different from that of the (1,1) block). 2265 2266 It returns a null array if the fieldsplit is not of type PC_COMPOSITE_SCHUR; in this case, you should use PCFieldSplitGetSubKSP(). 2267 2268 Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs. 2269 You can call PCFieldSplitSchurGetSubKSP(pc,n,PETSC_NULL_KSP,ierr) to determine how large the 2270 KSP array must be. 2271 2272 Level: advanced 2273 2274 .seealso: `PCFIELDSPLIT` 2275 @*/ 2276 PetscErrorCode PCFieldSplitSchurGetSubKSP(PC pc,PetscInt *n,KSP *subksp[]) 2277 { 2278 PetscFunctionBegin; 2279 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2280 if (n) PetscValidIntPointer(n,2); 2281 PetscUseMethod(pc,"PCFieldSplitSchurGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp)); 2282 PetscFunctionReturn(0); 2283 } 2284 2285 /*@ 2286 PCFieldSplitSetSchurPre - Indicates from what operator the preconditioner is constructucted for the Schur complement. 2287 The default is the A11 matrix. 2288 2289 Collective on PC 2290 2291 Input Parameters: 2292 + pc - the preconditioner context 2293 . 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 2294 PC_FIELDSPLIT_SCHUR_PRE_SELFP, and PC_FIELDSPLIT_SCHUR_PRE_FULL 2295 - userpre - matrix to use for preconditioning, or NULL 2296 2297 Options Database: 2298 + -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is a11. See notes for meaning of various arguments 2299 - -fieldsplit_1_pc_type <pctype> - the preconditioner algorithm that is used to construct the preconditioner from the operator 2300 2301 Notes: 2302 $ If ptype is 2303 $ a11 - the preconditioner for the Schur complement is generated from the block diagonal part of the preconditioner 2304 $ matrix associated with the Schur complement (i.e. A11), not the Schur complement matrix 2305 $ self - the preconditioner for the Schur complement is generated from the symbolic representation of the Schur complement matrix: 2306 $ The only preconditioner that currently works with this symbolic respresentation matrix object is the PCLSC 2307 $ preconditioner 2308 $ user - the preconditioner for the Schur complement is generated from the user provided matrix (pre argument 2309 $ to this function). 2310 $ selfp - the preconditioning for the Schur complement is generated from an explicitly-assembled approximation Sp = A11 - A10 inv(diag(A00)) A01 2311 $ This is only a good preconditioner when diag(A00) is a good preconditioner for A00. Optionally, A00 can be 2312 $ lumped before extracting the diagonal using the additional option -fieldsplit_1_mat_schur_complement_ainv_type lump 2313 $ full - the preconditioner for the Schur complement is generated from the exact Schur complement matrix representation computed internally by PCFIELDSPLIT (this is expensive) 2314 $ useful mostly as a test that the Schur complement approach can work for your problem 2315 2316 When solving a saddle point problem, where the A11 block is identically zero, using a11 as the ptype only makes sense 2317 with the additional option -fieldsplit_1_pc_type none. Usually for saddle point problems one would use a ptype of self and 2318 -fieldsplit_1_pc_type lsc which uses the least squares commutator to compute a preconditioner for the Schur complement. 2319 2320 Level: intermediate 2321 2322 .seealso: `PCFieldSplitGetSchurPre()`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurPreType`, 2323 `MatSchurComplementSetAinvType()`, `PCLSC` 2324 2325 @*/ 2326 PetscErrorCode PCFieldSplitSetSchurPre(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 2327 { 2328 PetscFunctionBegin; 2329 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2330 PetscTryMethod(pc,"PCFieldSplitSetSchurPre_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre)); 2331 PetscFunctionReturn(0); 2332 } 2333 2334 PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) {return PCFieldSplitSetSchurPre(pc,ptype,pre);} /* Deprecated name */ 2335 2336 /*@ 2337 PCFieldSplitGetSchurPre - For Schur complement fieldsplit, determine how the Schur complement will be 2338 preconditioned. See PCFieldSplitSetSchurPre() for details. 2339 2340 Logically Collective on PC 2341 2342 Input Parameter: 2343 . pc - the preconditioner context 2344 2345 Output Parameters: 2346 + ptype - which matrix to use for preconditioning the Schur complement: PC_FIELDSPLIT_SCHUR_PRE_A11, PC_FIELDSPLIT_SCHUR_PRE_SELF, PC_FIELDSPLIT_PRE_USER 2347 - userpre - matrix to use for preconditioning (with PC_FIELDSPLIT_PRE_USER), or NULL 2348 2349 Level: intermediate 2350 2351 .seealso: `PCFieldSplitSetSchurPre()`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurPreType`, `PCLSC` 2352 2353 @*/ 2354 PetscErrorCode PCFieldSplitGetSchurPre(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre) 2355 { 2356 PetscFunctionBegin; 2357 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2358 PetscUseMethod(pc,"PCFieldSplitGetSchurPre_C",(PC,PCFieldSplitSchurPreType*,Mat*),(pc,ptype,pre)); 2359 PetscFunctionReturn(0); 2360 } 2361 2362 /*@ 2363 PCFieldSplitSchurGetS - extract the MatSchurComplement object used by this PC in case it needs to be configured separately 2364 2365 Not collective 2366 2367 Input Parameter: 2368 . pc - the preconditioner context 2369 2370 Output Parameter: 2371 . S - the Schur complement matrix 2372 2373 Notes: 2374 This matrix should not be destroyed using MatDestroy(); rather, use PCFieldSplitSchurRestoreS(). 2375 2376 Level: advanced 2377 2378 .seealso: `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSchurPreType`, `PCFieldSplitSetSchurPre()`, `MatSchurComplement`, `PCFieldSplitSchurRestoreS()` 2379 2380 @*/ 2381 PetscErrorCode PCFieldSplitSchurGetS(PC pc,Mat *S) 2382 { 2383 const char* t; 2384 PetscBool isfs; 2385 PC_FieldSplit *jac; 2386 2387 PetscFunctionBegin; 2388 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2389 PetscCall(PetscObjectGetType((PetscObject)pc,&t)); 2390 PetscCall(PetscStrcmp(t,PCFIELDSPLIT,&isfs)); 2391 PetscCheck(isfs,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t); 2392 jac = (PC_FieldSplit*)pc->data; 2393 PetscCheck(jac->type == PC_COMPOSITE_SCHUR,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %d instead",jac->type); 2394 if (S) *S = jac->schur; 2395 PetscFunctionReturn(0); 2396 } 2397 2398 /*@ 2399 PCFieldSplitSchurRestoreS - restores the MatSchurComplement object used by this PC 2400 2401 Not collective 2402 2403 Input Parameters: 2404 + pc - the preconditioner context 2405 - S - the Schur complement matrix 2406 2407 Level: advanced 2408 2409 .seealso: `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSchurPreType`, `PCFieldSplitSetSchurPre()`, `MatSchurComplement`, `PCFieldSplitSchurGetS()` 2410 2411 @*/ 2412 PetscErrorCode PCFieldSplitSchurRestoreS(PC pc,Mat *S) 2413 { 2414 const char* t; 2415 PetscBool isfs; 2416 PC_FieldSplit *jac; 2417 2418 PetscFunctionBegin; 2419 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2420 PetscCall(PetscObjectGetType((PetscObject)pc,&t)); 2421 PetscCall(PetscStrcmp(t,PCFIELDSPLIT,&isfs)); 2422 PetscCheck(isfs,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t); 2423 jac = (PC_FieldSplit*)pc->data; 2424 PetscCheck(jac->type == PC_COMPOSITE_SCHUR,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %d instead",jac->type); 2425 PetscCheck(S && (*S == jac->schur),PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"MatSchurComplement restored is not the same as gotten"); 2426 PetscFunctionReturn(0); 2427 } 2428 2429 static PetscErrorCode PCFieldSplitSetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 2430 { 2431 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2432 2433 PetscFunctionBegin; 2434 jac->schurpre = ptype; 2435 if (ptype == PC_FIELDSPLIT_SCHUR_PRE_USER && pre) { 2436 PetscCall(MatDestroy(&jac->schur_user)); 2437 jac->schur_user = pre; 2438 PetscCall(PetscObjectReference((PetscObject)jac->schur_user)); 2439 } 2440 PetscFunctionReturn(0); 2441 } 2442 2443 static PetscErrorCode PCFieldSplitGetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre) 2444 { 2445 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2446 2447 PetscFunctionBegin; 2448 *ptype = jac->schurpre; 2449 *pre = jac->schur_user; 2450 PetscFunctionReturn(0); 2451 } 2452 2453 /*@ 2454 PCFieldSplitSetSchurFactType - sets which blocks of the approximate block factorization to retain in the preconditioner 2455 2456 Collective on PC 2457 2458 Input Parameters: 2459 + pc - the preconditioner context 2460 - ftype - which blocks of factorization to retain, PC_FIELDSPLIT_SCHUR_FACT_FULL is default 2461 2462 Options Database: 2463 . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> - default is full 2464 2465 Level: intermediate 2466 2467 Notes: 2468 The FULL factorization is 2469 2470 .vb 2471 (A B) = (1 0) (A 0) (1 Ainv*B) = L D U 2472 (C E) (C*Ainv 1) (0 S) (0 1) 2473 .vb 2474 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, 2475 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(). 2476 2477 If A and S are solved exactly 2478 .vb 2479 *) FULL factorization is a direct solver. 2480 *) 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. 2481 *) 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. 2482 .ve 2483 2484 If the iteration count is very low, consider using KSPFGMRES or KSPGCR which can use one less preconditioner 2485 application in this case. Note that the preconditioned operator may be highly non-normal, so such fast convergence may not be observed in practice. 2486 2487 For symmetric problems in which A is positive definite and S is negative definite, DIAG can be used with KSPMINRES. 2488 2489 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). 2490 2491 References: 2492 + * - Murphy, Golub, and Wathen, A note on preconditioning indefinite linear systems, SIAM J. Sci. Comput., 21 (2000). 2493 - * - Ipsen, A note on preconditioning nonsymmetric matrices, SIAM J. Sci. Comput., 23 (2001). 2494 2495 .seealso: `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurPreType`, `PCFieldSplitSetSchurScale()` 2496 @*/ 2497 PetscErrorCode PCFieldSplitSetSchurFactType(PC pc,PCFieldSplitSchurFactType ftype) 2498 { 2499 PetscFunctionBegin; 2500 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2501 PetscTryMethod(pc,"PCFieldSplitSetSchurFactType_C",(PC,PCFieldSplitSchurFactType),(pc,ftype)); 2502 PetscFunctionReturn(0); 2503 } 2504 2505 static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc,PCFieldSplitSchurFactType ftype) 2506 { 2507 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2508 2509 PetscFunctionBegin; 2510 jac->schurfactorization = ftype; 2511 PetscFunctionReturn(0); 2512 } 2513 2514 /*@ 2515 PCFieldSplitSetSchurScale - Controls the sign flip of S for PC_FIELDSPLIT_SCHUR_FACT_DIAG. 2516 2517 Collective on PC 2518 2519 Input Parameters: 2520 + pc - the preconditioner context 2521 - scale - scaling factor for the Schur complement 2522 2523 Options Database: 2524 . -pc_fieldsplit_schur_scale - default is -1.0 2525 2526 Level: intermediate 2527 2528 .seealso: `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurFactType`, `PCFieldSplitSetSchurScale()` 2529 @*/ 2530 PetscErrorCode PCFieldSplitSetSchurScale(PC pc,PetscScalar scale) 2531 { 2532 PetscFunctionBegin; 2533 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2534 PetscValidLogicalCollectiveScalar(pc,scale,2); 2535 PetscTryMethod(pc,"PCFieldSplitSetSchurScale_C",(PC,PetscScalar),(pc,scale)); 2536 PetscFunctionReturn(0); 2537 } 2538 2539 static PetscErrorCode PCFieldSplitSetSchurScale_FieldSplit(PC pc,PetscScalar scale) 2540 { 2541 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2542 2543 PetscFunctionBegin; 2544 jac->schurscale = scale; 2545 PetscFunctionReturn(0); 2546 } 2547 2548 /*@C 2549 PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement 2550 2551 Collective on KSP 2552 2553 Input Parameter: 2554 . pc - the preconditioner context 2555 2556 Output Parameters: 2557 + A00 - the (0,0) block 2558 . A01 - the (0,1) block 2559 . A10 - the (1,0) block 2560 - A11 - the (1,1) block 2561 2562 Level: advanced 2563 2564 .seealso: `PCFIELDSPLIT` 2565 @*/ 2566 PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11) 2567 { 2568 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 2569 2570 PetscFunctionBegin; 2571 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2572 PetscCheck(jac->type == PC_COMPOSITE_SCHUR,PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach."); 2573 if (A00) *A00 = jac->pmat[0]; 2574 if (A01) *A01 = jac->B; 2575 if (A10) *A10 = jac->C; 2576 if (A11) *A11 = jac->pmat[1]; 2577 PetscFunctionReturn(0); 2578 } 2579 2580 /*@ 2581 PCFieldSplitSetGKBTol - Sets the solver tolerance for the generalized Golub-Kahan bidiagonalization preconditioner. 2582 2583 Collective on PC 2584 2585 Notes: 2586 The generalized GKB algorithm uses a lower bound estimate of the error in energy norm as stopping criterion. 2587 It stops once the lower bound estimate undershoots the required solver tolerance. Although the actual error might be bigger than 2588 this estimate, the stopping criterion is satisfactory in practical cases [A13]. 2589 2590 [Ar13] Generalized Golub-Kahan bidiagonalization and stopping criteria, SIAM J. Matrix Anal. Appl., Vol. 34, No. 2, pp. 571-592, 2013. 2591 2592 Input Parameters: 2593 + pc - the preconditioner context 2594 - tolerance - the solver tolerance 2595 2596 Options Database: 2597 . -pc_fieldsplit_gkb_tol - default is 1e-5 2598 2599 Level: intermediate 2600 2601 .seealso: `PCFIELDSPLIT`, `PCFieldSplitSetGKBDelay()`, `PCFieldSplitSetGKBNu()`, `PCFieldSplitSetGKBMaxit()` 2602 @*/ 2603 PetscErrorCode PCFieldSplitSetGKBTol(PC pc,PetscReal tolerance) 2604 { 2605 PetscFunctionBegin; 2606 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2607 PetscValidLogicalCollectiveReal(pc,tolerance,2); 2608 PetscTryMethod(pc,"PCFieldSplitSetGKBTol_C",(PC,PetscReal),(pc,tolerance)); 2609 PetscFunctionReturn(0); 2610 } 2611 2612 static PetscErrorCode PCFieldSplitSetGKBTol_FieldSplit(PC pc,PetscReal tolerance) 2613 { 2614 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2615 2616 PetscFunctionBegin; 2617 jac->gkbtol = tolerance; 2618 PetscFunctionReturn(0); 2619 } 2620 2621 /*@ 2622 PCFieldSplitSetGKBMaxit - Sets the maximum number of iterations for the generalized Golub-Kahan bidiagonalization 2623 preconditioner. 2624 2625 Collective on PC 2626 2627 Input Parameters: 2628 + pc - the preconditioner context 2629 - maxit - the maximum number of iterations 2630 2631 Options Database: 2632 . -pc_fieldsplit_gkb_maxit - default is 100 2633 2634 Level: intermediate 2635 2636 .seealso: `PCFIELDSPLIT`, `PCFieldSplitSetGKBDelay()`, `PCFieldSplitSetGKBTol()`, `PCFieldSplitSetGKBNu()` 2637 @*/ 2638 PetscErrorCode PCFieldSplitSetGKBMaxit(PC pc,PetscInt maxit) 2639 { 2640 PetscFunctionBegin; 2641 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2642 PetscValidLogicalCollectiveInt(pc,maxit,2); 2643 PetscTryMethod(pc,"PCFieldSplitSetGKBMaxit_C",(PC,PetscInt),(pc,maxit)); 2644 PetscFunctionReturn(0); 2645 } 2646 2647 static PetscErrorCode PCFieldSplitSetGKBMaxit_FieldSplit(PC pc,PetscInt maxit) 2648 { 2649 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2650 2651 PetscFunctionBegin; 2652 jac->gkbmaxit = maxit; 2653 PetscFunctionReturn(0); 2654 } 2655 2656 /*@ 2657 PCFieldSplitSetGKBDelay - Sets the delay in the lower bound error estimate in the generalized Golub-Kahan bidiagonalization 2658 preconditioner. 2659 2660 Collective on PC 2661 2662 Notes: 2663 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 2664 is expressed as a truncated sum. The error at iteration k can only be measured at iteration (k + delay), and thus the algorithm needs 2665 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 2666 2667 [Ar13] Generalized Golub-Kahan bidiagonalization and stopping criteria, SIAM J. Matrix Anal. Appl., Vol. 34, No. 2, pp. 571-592, 2013. 2668 2669 Input Parameters: 2670 + pc - the preconditioner context 2671 - delay - the delay window in the lower bound estimate 2672 2673 Options Database: 2674 . -pc_fieldsplit_gkb_delay - default is 5 2675 2676 Level: intermediate 2677 2678 .seealso: `PCFIELDSPLIT`, `PCFieldSplitSetGKBNu()`, `PCFieldSplitSetGKBTol()`, `PCFieldSplitSetGKBMaxit()` 2679 @*/ 2680 PetscErrorCode PCFieldSplitSetGKBDelay(PC pc,PetscInt delay) 2681 { 2682 PetscFunctionBegin; 2683 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2684 PetscValidLogicalCollectiveInt(pc,delay,2); 2685 PetscTryMethod(pc,"PCFieldSplitSetGKBDelay_C",(PC,PetscInt),(pc,delay)); 2686 PetscFunctionReturn(0); 2687 } 2688 2689 static PetscErrorCode PCFieldSplitSetGKBDelay_FieldSplit(PC pc,PetscInt delay) 2690 { 2691 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2692 2693 PetscFunctionBegin; 2694 jac->gkbdelay = delay; 2695 PetscFunctionReturn(0); 2696 } 2697 2698 /*@ 2699 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. 2700 2701 Collective on PC 2702 2703 Notes: 2704 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, 2705 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 2706 necessary to find a good balance in between the convergence of the inner and outer loop. 2707 2708 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. 2709 2710 [Ar13] Generalized Golub-Kahan bidiagonalization and stopping criteria, SIAM J. Matrix Anal. Appl., Vol. 34, No. 2, pp. 571-592, 2013. 2711 2712 Input Parameters: 2713 + pc - the preconditioner context 2714 - nu - the shift parameter 2715 2716 Options Database: 2717 . -pc_fieldsplit_gkb_nu - default is 1 2718 2719 Level: intermediate 2720 2721 .seealso: `PCFIELDSPLIT`, `PCFieldSplitSetGKBDelay()`, `PCFieldSplitSetGKBTol()`, `PCFieldSplitSetGKBMaxit()` 2722 @*/ 2723 PetscErrorCode PCFieldSplitSetGKBNu(PC pc,PetscReal nu) 2724 { 2725 PetscFunctionBegin; 2726 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2727 PetscValidLogicalCollectiveReal(pc,nu,2); 2728 PetscTryMethod(pc,"PCFieldSplitSetGKBNu_C",(PC,PetscReal),(pc,nu)); 2729 PetscFunctionReturn(0); 2730 } 2731 2732 static PetscErrorCode PCFieldSplitSetGKBNu_FieldSplit(PC pc,PetscReal nu) 2733 { 2734 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2735 2736 PetscFunctionBegin; 2737 jac->gkbnu = nu; 2738 PetscFunctionReturn(0); 2739 } 2740 2741 static PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type) 2742 { 2743 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2744 2745 PetscFunctionBegin; 2746 jac->type = type; 2747 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",NULL)); 2748 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",NULL)); 2749 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",NULL)); 2750 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",NULL)); 2751 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurScale_C",NULL)); 2752 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBTol_C",NULL)); 2753 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBMaxit_C",NULL)); 2754 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBNu_C",NULL)); 2755 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBDelay_C",NULL)); 2756 2757 if (type == PC_COMPOSITE_SCHUR) { 2758 pc->ops->apply = PCApply_FieldSplit_Schur; 2759 pc->ops->view = PCView_FieldSplit_Schur; 2760 2761 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit_Schur)); 2762 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",PCFieldSplitSetSchurPre_FieldSplit)); 2763 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",PCFieldSplitGetSchurPre_FieldSplit)); 2764 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",PCFieldSplitSetSchurFactType_FieldSplit)); 2765 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurScale_C",PCFieldSplitSetSchurScale_FieldSplit)); 2766 } else if (type == PC_COMPOSITE_GKB) { 2767 pc->ops->apply = PCApply_FieldSplit_GKB; 2768 pc->ops->view = PCView_FieldSplit_GKB; 2769 2770 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit)); 2771 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBTol_C",PCFieldSplitSetGKBTol_FieldSplit)); 2772 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBMaxit_C",PCFieldSplitSetGKBMaxit_FieldSplit)); 2773 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBNu_C",PCFieldSplitSetGKBNu_FieldSplit)); 2774 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetGKBDelay_C",PCFieldSplitSetGKBDelay_FieldSplit)); 2775 } else { 2776 pc->ops->apply = PCApply_FieldSplit; 2777 pc->ops->view = PCView_FieldSplit; 2778 2779 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit)); 2780 } 2781 PetscFunctionReturn(0); 2782 } 2783 2784 static PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs) 2785 { 2786 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2787 2788 PetscFunctionBegin; 2789 PetscCheck(bs >= 1,PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %" PetscInt_FMT,bs); 2790 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); 2791 jac->bs = bs; 2792 PetscFunctionReturn(0); 2793 } 2794 2795 static PetscErrorCode PCSetCoordinates_FieldSplit(PC pc, PetscInt dim, PetscInt nloc, PetscReal coords[]) 2796 { 2797 PC_FieldSplit * jac = (PC_FieldSplit*)pc->data; 2798 PC_FieldSplitLink ilink_current = jac->head; 2799 PetscInt ii; 2800 IS is_owned; 2801 2802 PetscFunctionBegin; 2803 jac->coordinates_set = PETSC_TRUE; // Internal flag 2804 PetscCall(MatGetOwnershipIS(pc->mat,&is_owned,PETSC_NULL)); 2805 2806 ii=0; 2807 while (ilink_current) { 2808 // For each IS, embed it to get local coords indces 2809 IS is_coords; 2810 PetscInt ndofs_block; 2811 const PetscInt *block_dofs_enumeration; // Numbering of the dofs relevant to the current block 2812 2813 // Setting drop to true for safety. It should make no difference. 2814 PetscCall(ISEmbed(ilink_current->is, is_owned, PETSC_TRUE, &is_coords)); 2815 PetscCall(ISGetLocalSize(is_coords, &ndofs_block)); 2816 PetscCall(ISGetIndices(is_coords, &block_dofs_enumeration)); 2817 2818 // Allocate coordinates vector and set it directly 2819 PetscCall(PetscMalloc1(ndofs_block * dim, &(ilink_current->coords))); 2820 for (PetscInt dof=0;dof<ndofs_block;++dof) { 2821 for (PetscInt d=0;d<dim;++d) { 2822 (ilink_current->coords)[dim*dof + d] = coords[dim * block_dofs_enumeration[dof] + d]; 2823 } 2824 } 2825 ilink_current->dim = dim; 2826 ilink_current->ndofs = ndofs_block; 2827 PetscCall(ISRestoreIndices(is_coords, &block_dofs_enumeration)); 2828 PetscCall(ISDestroy(&is_coords)); 2829 ilink_current = ilink_current->next; 2830 ++ii; 2831 } 2832 PetscCall(ISDestroy(&is_owned)); 2833 PetscFunctionReturn(0); 2834 } 2835 2836 /*@ 2837 PCFieldSplitSetType - Sets the type of fieldsplit preconditioner. 2838 2839 Collective on PC 2840 2841 Input Parameters: 2842 + pc - the preconditioner context 2843 - type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 2844 2845 Options Database Key: 2846 . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type 2847 2848 Level: Intermediate 2849 2850 .seealso: `PCCompositeSetType()` 2851 2852 @*/ 2853 PetscErrorCode PCFieldSplitSetType(PC pc,PCCompositeType type) 2854 { 2855 PetscFunctionBegin; 2856 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2857 PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type)); 2858 PetscFunctionReturn(0); 2859 } 2860 2861 /*@ 2862 PCFieldSplitGetType - Gets the type of fieldsplit preconditioner. 2863 2864 Not collective 2865 2866 Input Parameter: 2867 . pc - the preconditioner context 2868 2869 Output Parameter: 2870 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 2871 2872 Level: Intermediate 2873 2874 .seealso: `PCCompositeSetType()` 2875 @*/ 2876 PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type) 2877 { 2878 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 2879 2880 PetscFunctionBegin; 2881 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2882 PetscValidIntPointer(type,2); 2883 *type = jac->type; 2884 PetscFunctionReturn(0); 2885 } 2886 2887 /*@ 2888 PCFieldSplitSetDMSplits - Flags whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible. 2889 2890 Logically Collective 2891 2892 Input Parameters: 2893 + pc - the preconditioner context 2894 - flg - boolean indicating whether to use field splits defined by the DM 2895 2896 Options Database Key: 2897 . -pc_fieldsplit_dm_splits <bool> - use the field splits defined by the DM 2898 2899 Level: Intermediate 2900 2901 .seealso: `PCFieldSplitGetDMSplits()` 2902 2903 @*/ 2904 PetscErrorCode PCFieldSplitSetDMSplits(PC pc,PetscBool flg) 2905 { 2906 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2907 PetscBool isfs; 2908 2909 PetscFunctionBegin; 2910 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2911 PetscValidLogicalCollectiveBool(pc,flg,2); 2912 PetscCall(PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs)); 2913 if (isfs) { 2914 jac->dm_splits = flg; 2915 } 2916 PetscFunctionReturn(0); 2917 } 2918 2919 /*@ 2920 PCFieldSplitGetDMSplits - Returns flag indicating whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible. 2921 2922 Logically Collective 2923 2924 Input Parameter: 2925 . pc - the preconditioner context 2926 2927 Output Parameter: 2928 . flg - boolean indicating whether to use field splits defined by the DM 2929 2930 Level: Intermediate 2931 2932 .seealso: `PCFieldSplitSetDMSplits()` 2933 2934 @*/ 2935 PetscErrorCode PCFieldSplitGetDMSplits(PC pc,PetscBool* flg) 2936 { 2937 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2938 PetscBool isfs; 2939 2940 PetscFunctionBegin; 2941 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2942 PetscValidBoolPointer(flg,2); 2943 PetscCall(PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs)); 2944 if (isfs) { 2945 if (flg) *flg = jac->dm_splits; 2946 } 2947 PetscFunctionReturn(0); 2948 } 2949 2950 /*@ 2951 PCFieldSplitGetDetectSaddlePoint - Returns flag indicating whether PCFieldSplit will attempt to automatically determine fields based on zero diagonal entries. 2952 2953 Logically Collective 2954 2955 Input Parameter: 2956 . pc - the preconditioner context 2957 2958 Output Parameter: 2959 . flg - boolean indicating whether to detect fields or not 2960 2961 Level: Intermediate 2962 2963 .seealso: `PCFIELDSPLIT`, `PCFieldSplitSetDetectSaddlePoint()` 2964 2965 @*/ 2966 PetscErrorCode PCFieldSplitGetDetectSaddlePoint(PC pc,PetscBool *flg) 2967 { 2968 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2969 2970 PetscFunctionBegin; 2971 *flg = jac->detect; 2972 PetscFunctionReturn(0); 2973 } 2974 2975 /*@ 2976 PCFieldSplitSetDetectSaddlePoint - Sets flag indicating whether PCFieldSplit will attempt to automatically determine fields based on zero diagonal entries. 2977 2978 Logically Collective 2979 2980 Input Parameter: 2981 . pc - the preconditioner context 2982 2983 Output Parameter: 2984 . flg - boolean indicating whether to detect fields or not 2985 2986 Options Database Key: 2987 . -pc_fieldsplit_detect_saddle_point <bool> - detect and use the saddle point 2988 2989 Notes: 2990 Also sets the split type to PC_COMPOSITE_SCHUR (see PCFieldSplitSetType()) and the Schur preconditioner type to PC_FIELDSPLIT_SCHUR_PRE_SELF (see PCFieldSplitSetSchurPre()). 2991 2992 Level: Intermediate 2993 2994 .seealso: `PCFIELDSPLIT`, `PCFieldSplitSetDetectSaddlePoint()`, `PCFieldSplitSetType()`, `PCFieldSplitSetSchurPre()` 2995 2996 @*/ 2997 PetscErrorCode PCFieldSplitSetDetectSaddlePoint(PC pc,PetscBool flg) 2998 { 2999 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 3000 3001 PetscFunctionBegin; 3002 jac->detect = flg; 3003 if (jac->detect) { 3004 PetscCall(PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR)); 3005 PetscCall(PCFieldSplitSetSchurPre(pc,PC_FIELDSPLIT_SCHUR_PRE_SELF,NULL)); 3006 } 3007 PetscFunctionReturn(0); 3008 } 3009 3010 /* -------------------------------------------------------------------------------------*/ 3011 /*MC 3012 PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual 3013 fields or groups of fields. See [the users manual section on "Solving Block Matrices"](sec_block_matrices) for more details. 3014 3015 To set options on the solvers for each block append `-fieldsplit_` to all the PC 3016 options database keys. For example, `-fieldsplit_pc_type ilu -fieldsplit_pc_factor_levels 1` 3017 3018 To set the options on the solvers separate for each block call `PCFieldSplitGetSubKSP()` 3019 and set the options directly on the resulting `KSP` object 3020 3021 Level: intermediate 3022 3023 Options Database Keys: 3024 + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the `%d`'th split 3025 . -pc_fieldsplit_default - automatically add any fields to additional splits that have not 3026 been supplied explicitly by `-pc_fieldsplit_%d_fields` 3027 . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields) 3028 . -pc_fieldsplit_type <additive,multiplicative,symmetric_multiplicative,schur,gkb> - type of relaxation or factorization splitting 3029 . -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is a11; see `PCFieldSplitSetSchurPre()` 3030 . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> - set factorization type when using `-pc_fieldsplit_type schur`; see `PCFieldSplitSetSchurFactType()` 3031 - -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero diagonal and uses Schur complement with no preconditioner as the solver 3032 3033 Options prefixes for inner solvers when using the Schur complement preconditioner are `-fieldsplit_0_` and `-fieldsplit_1_` . 3034 For all other solvers they are `-fieldsplit_%d_` for the `d`th field; use `-fieldsplit_` for all fields. 3035 The options prefix for the inner solver when using the Golub-Kahan biadiagonalization preconditioner is `-fieldsplit_0_` 3036 3037 Notes: 3038 Use `PCFieldSplitSetFields()` to set fields defined by "strided" entries and `PCFieldSplitSetIS()` 3039 to define a field by an arbitrary collection of entries. 3040 3041 If no fields are set the default is used. The fields are defined by entries strided by bs, 3042 beginning at 0 then 1, etc to bs-1. The block size can be set with `PCFieldSplitSetBlockSize()`, 3043 if this is not called the block size defaults to the blocksize of the second matrix passed 3044 to `KSPSetOperators()`/`PCSetOperators()`. 3045 3046 For the Schur complement preconditioner if 3047 3048 ```{math} 3049 J = \left[\begin{array}{cc} A_{00} & A_{01} \\ A_{10} & A_{11} \end{array}\right] 3050 ``` 3051 3052 the preconditioner using `full` factorization is logically 3053 ```{math} 3054 \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] 3055 ``` 3056 where the action of $\text{inv}(A_{00})$ is applied using the KSP solver with prefix `-fieldsplit_0_`. $S$ is the Schur complement 3057 ```{math} 3058 S = A_{11} - A_{10} \text{ksp}(A_{00}) A_{01} 3059 ``` 3060 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 3061 in providing the SECOND split or 1 if not given). For `PCFieldSplitGetSub\text{ksp}()` when field number is 0, 3062 it returns the KSP associated with `-fieldsplit_0_` while field number 1 gives `-fieldsplit_1_` KSP. By default 3063 $A_{11}$ is used to construct a preconditioner for $S$, use `PCFieldSplitSetSchurPre()` for all the possible ways to construct the preconditioner for $S$. 3064 3065 The factorization type is set using `-pc_fieldsplit_schur_fact_type <diag, lower, upper, full>`. `full` is shown above, 3066 `diag` gives 3067 ```{math} 3068 \left[\begin{array}{cc} \text{inv}(A_{00}) & 0 \\ 0 & -\text{ksp}(S) \end{array}\right] 3069 ``` 3070 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 3071 can be turned off with `PCFieldSplitSetSchurScale()` or by command line `-pc_fieldsplit_schur_scale 1.0`. The `lower` factorization is the inverse of 3072 ```{math} 3073 \left[\begin{array}{cc} A_{00} & 0 \\ A_{10} & S \end{array}\right] 3074 ``` 3075 where the inverses of A_{00} and S are applied using KSPs. The upper factorization is the inverse of 3076 ```{math} 3077 \left[\begin{array}{cc} A_{00} & A_{01} \\ 0 & S \end{array}\right] 3078 ``` 3079 where again the inverses of $A_{00}$ and $S$ are applied using `KSP`s. 3080 3081 If only one set of indices (one `IS`) is provided with `PCFieldSplitSetIS()` then the complement of that `IS` 3082 is used automatically for a second block. 3083 3084 The fieldsplit preconditioner cannot currently be used with the BAIJ or SBAIJ data formats if the blocksize is larger than 1. 3085 Generally it should be used with the AIJ format. 3086 3087 The forms of these preconditioners are closely related if not identical to forms derived as "Distributive Iterations", see, 3088 for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling {cite}`Wesseling2009`. Note that one can also use `PCFIELDSPLIT` 3089 inside a smoother resulting in "Distributive Smoothers". 3090 3091 References: 3092 3093 See "A taxonomy and comparison of parallel block multi-level preconditioners for the incompressible Navier-Stokes equations" {cite}`elman2008tcp`. 3094 3095 The Constrained Pressure Preconditioner (CPR) can be implemented using `PCCOMPOSITE` with `PCGALERKIN`. CPR first solves an $R A P$ subsystem, updates the 3096 residual on all variables (`PCCompositeSetType(pc,PC_COMPOSITE_MULTIPLICATIVE)`), and then applies a simple ILU like preconditioner on all the variables. 3097 3098 The generalized Golub-Kahan bidiagonalization preconditioner (GKB) can be applied to symmetric $2 \times 2$ block matrices of the shape 3099 ```{math} 3100 \left[\begin{array}{cc} A_{00} & A_{01} \\ A_{01}' & 0 \end{array}\right] 3101 ``` 3102 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}'$. 3103 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_`. 3104 3105 ```{bibliography} 3106 :filter: docname in docnames 3107 ``` 3108 3109 .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCLSC`, 3110 `PCFieldSplitGetSubKSP()`, `PCFieldSplitSchurGetSubKSP()`, `PCFieldSplitSetFields()`, 3111 `PCFieldSplitSetType()`, `PCFieldSplitSetIS()`, `PCFieldSplitSetSchurPre()`, `PCFieldSplitSetSchurFactType()`, 3112 `MatSchurComplementSetAinvType()`, `PCFieldSplitSetSchurScale()`, `PCFieldSplitSetDetectSaddlePoint()` 3113 M*/ 3114 3115 PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc) 3116 { 3117 PC_FieldSplit *jac; 3118 3119 PetscFunctionBegin; 3120 PetscCall(PetscNewLog(pc,&jac)); 3121 3122 jac->bs = -1; 3123 jac->nsplits = 0; 3124 jac->type = PC_COMPOSITE_MULTIPLICATIVE; 3125 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */ 3126 jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL; 3127 jac->schurscale = -1.0; 3128 jac->dm_splits = PETSC_TRUE; 3129 jac->detect = PETSC_FALSE; 3130 jac->gkbtol = 1e-5; 3131 jac->gkbdelay = 5; 3132 jac->gkbnu = 1; 3133 jac->gkbmaxit = 100; 3134 jac->gkbmonitor = PETSC_FALSE; 3135 jac->coordinates_set = PETSC_FALSE; 3136 3137 pc->data = (void*)jac; 3138 3139 pc->ops->apply = PCApply_FieldSplit; 3140 pc->ops->applytranspose = PCApplyTranspose_FieldSplit; 3141 pc->ops->setup = PCSetUp_FieldSplit; 3142 pc->ops->reset = PCReset_FieldSplit; 3143 pc->ops->destroy = PCDestroy_FieldSplit; 3144 pc->ops->setfromoptions = PCSetFromOptions_FieldSplit; 3145 pc->ops->view = PCView_FieldSplit; 3146 pc->ops->applyrichardson = NULL; 3147 3148 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSchurGetSubKSP_C",PCFieldSplitSchurGetSubKSP_FieldSplit)); 3149 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit)); 3150 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",PCFieldSplitSetFields_FieldSplit)); 3151 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",PCFieldSplitSetIS_FieldSplit)); 3152 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",PCFieldSplitSetType_FieldSplit)); 3153 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",PCFieldSplitSetBlockSize_FieldSplit)); 3154 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitRestrictIS_C",PCFieldSplitRestrictIS_FieldSplit)); 3155 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",PCSetCoordinates_FieldSplit)); 3156 PetscFunctionReturn(0); 3157 } 3158