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