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