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