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