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