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