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