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