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