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