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