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