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