1 #include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/ 2 #include <petsc/private/kspimpl.h> /* This is needed to provide the appropriate PETSC_EXTERN for KSP_Solve_FS ....*/ 3 #include <petscdm.h> 4 5 const char *const PCFieldSplitSchurPreTypes[] = {"SELF","SELFP","A11","USER","FULL","PCFieldSplitSchurPreType","PC_FIELDSPLIT_SCHUR_PRE_",0}; 6 const char *const PCFieldSplitSchurFactTypes[] = {"DIAG","LOWER","UPPER","FULL","PCFieldSplitSchurFactType","PC_FIELDSPLIT_SCHUR_FACT_",0}; 7 8 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; 9 10 typedef struct _PC_FieldSplitLink *PC_FieldSplitLink; 11 struct _PC_FieldSplitLink { 12 KSP ksp; 13 Vec x,y,z; 14 char *splitname; 15 PetscInt nfields; 16 PetscInt *fields,*fields_col; 17 VecScatter sctx; 18 IS is,is_col; 19 PC_FieldSplitLink next,previous; 20 PetscLogEvent event; 21 }; 22 23 typedef struct { 24 PCCompositeType type; 25 PetscBool defaultsplit; /* Flag for a system with a set of 'k' scalar fields with the same layout (and bs = k) */ 26 PetscBool splitdefined; /* Flag is set after the splits have been defined, to prevent more splits from being added */ 27 PetscInt bs; /* Block size for IS and Mat structures */ 28 PetscInt nsplits; /* Number of field divisions defined */ 29 Vec *x,*y,w1,w2; 30 Mat *mat; /* The diagonal block for each split */ 31 Mat *pmat; /* The preconditioning diagonal block for each split */ 32 Mat *Afield; /* The rows of the matrix associated with each split */ 33 PetscBool issetup; 34 35 /* Only used when Schur complement preconditioning is used */ 36 Mat B; /* The (0,1) block */ 37 Mat C; /* The (1,0) block */ 38 Mat schur; /* The Schur complement S = A11 - A10 A00^{-1} A01, the KSP here, kspinner, is H_1 in [El08] */ 39 Mat schurp; /* Assembled approximation to S built by MatSchurComplement to be used as a preconditioning matrix when solving with S */ 40 Mat schur_user; /* User-provided preconditioning matrix for the Schur complement */ 41 PCFieldSplitSchurPreType schurpre; /* Determines which preconditioning matrix is used for the Schur complement */ 42 PCFieldSplitSchurFactType schurfactorization; 43 KSP kspschur; /* The solver for S */ 44 KSP kspupper; /* The solver for A in the upper diagonal part of the factorization (H_2 in [El08]) */ 45 PetscScalar schurscale; /* Scaling factor for the Schur complement solution with DIAG factorization */ 46 47 /* Only used when Golub-Kahan bidiagonalization preconditioning is used */ 48 Mat H; /* The modified matrix H = A00 + nu*A01*A01' */ 49 PetscReal gkbtol; /* Stopping tolerance for lower bound estimate */ 50 PetscInt gkbdelay; /* The delay window for the stopping criterion */ 51 PetscReal gkbnu; /* Parameter for augmented Lagrangian H = A + nu*A01*A01' */ 52 PetscInt gkbmaxit; /* Maximum number of iterations for outer loop */ 53 PetscBool gkbmonitor; /* Monitor for gkb iterations and the lower bound error */ 54 PetscViewer gkbviewer; /* Viewer context for gkbmonitor */ 55 Vec u,v,d,Hu; /* Work vectors for the GKB algorithm */ 56 PetscScalar *vecz; /* Contains intermediate values, eg for lower bound */ 57 58 PC_FieldSplitLink head; 59 PetscBool isrestrict; /* indicates PCFieldSplitRestrictIS() has been last called on this object, hack */ 60 PetscBool suboptionsset; /* Indicates that the KSPSetFromOptions() has been called on the sub-KSPs */ 61 PetscBool dm_splits; /* Whether to use DMCreateFieldDecomposition() whenever possible */ 62 PetscBool diag_use_amat; /* Whether to extract diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */ 63 PetscBool offdiag_use_amat; /* Whether to extract off-diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */ 64 PetscBool detect; /* Whether to form 2-way split by finding zero diagonal entries */ 65 } PC_FieldSplit; 66 67 /* 68 Notes: 69 there is no particular reason that pmat, x, and y are stored as arrays in PC_FieldSplit instead of 70 inside PC_FieldSplitLink, just historical. If you want to be able to add new fields after already using the 71 PC you could change this. 72 */ 73 74 /* This helper is so that setting a user-provided preconditioning matrix is orthogonal to choosing to use it. This way the 75 * application-provided FormJacobian can provide this matrix without interfering with the user's (command-line) choices. */ 76 static Mat FieldSplitSchurPre(PC_FieldSplit *jac) 77 { 78 switch (jac->schurpre) { 79 case PC_FIELDSPLIT_SCHUR_PRE_SELF: return jac->schur; 80 case PC_FIELDSPLIT_SCHUR_PRE_SELFP: return jac->schurp; 81 case PC_FIELDSPLIT_SCHUR_PRE_A11: return jac->pmat[1]; 82 case PC_FIELDSPLIT_SCHUR_PRE_FULL: /* We calculate this and store it in schur_user */ 83 case PC_FIELDSPLIT_SCHUR_PRE_USER: /* Use a user-provided matrix if it is given, otherwise diagonal block */ 84 default: 85 return jac->schur_user ? jac->schur_user : jac->pmat[1]; 86 } 87 } 88 89 90 #include <petscdraw.h> 91 static PetscErrorCode PCView_FieldSplit(PC pc,PetscViewer viewer) 92 { 93 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 94 PetscErrorCode ierr; 95 PetscBool iascii,isdraw; 96 PetscInt i,j; 97 PC_FieldSplitLink ilink = jac->head; 98 99 PetscFunctionBegin; 100 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 101 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 102 if (iascii) { 103 if (jac->bs > 0) { 104 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D, blocksize = %D\n",PCCompositeTypes[jac->type],jac->nsplits,jac->bs);CHKERRQ(ierr); 105 } else { 106 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D\n",PCCompositeTypes[jac->type],jac->nsplits);CHKERRQ(ierr); 107 } 108 if (pc->useAmat) { 109 ierr = PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for blocks\n");CHKERRQ(ierr); 110 } 111 if (jac->diag_use_amat) { 112 ierr = PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for diagonal blocks\n");CHKERRQ(ierr); 113 } 114 if (jac->offdiag_use_amat) { 115 ierr = PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for off-diagonal blocks\n");CHKERRQ(ierr); 116 } 117 ierr = PetscViewerASCIIPrintf(viewer," Solver info for each split is in the following KSP objects:\n");CHKERRQ(ierr); 118 for (i=0; i<jac->nsplits; i++) { 119 if (ilink->fields) { 120 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);CHKERRQ(ierr); 121 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 122 for (j=0; j<ilink->nfields; j++) { 123 if (j > 0) { 124 ierr = PetscViewerASCIIPrintf(viewer,",");CHKERRQ(ierr); 125 } 126 ierr = PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);CHKERRQ(ierr); 127 } 128 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 129 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 130 } else { 131 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i);CHKERRQ(ierr); 132 } 133 ierr = KSPView(ilink->ksp,viewer);CHKERRQ(ierr); 134 ilink = ilink->next; 135 } 136 } 137 138 if (isdraw) { 139 PetscDraw draw; 140 PetscReal x,y,w,wd; 141 142 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 143 ierr = PetscDrawGetCurrentPoint(draw,&x,&y);CHKERRQ(ierr); 144 w = 2*PetscMin(1.0 - x,x); 145 wd = w/(jac->nsplits + 1); 146 x = x - wd*(jac->nsplits-1)/2.0; 147 for (i=0; i<jac->nsplits; i++) { 148 ierr = PetscDrawPushCurrentPoint(draw,x,y);CHKERRQ(ierr); 149 ierr = KSPView(ilink->ksp,viewer);CHKERRQ(ierr); 150 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr); 151 x += wd; 152 ilink = ilink->next; 153 } 154 } 155 PetscFunctionReturn(0); 156 } 157 158 static PetscErrorCode PCView_FieldSplit_Schur(PC pc,PetscViewer viewer) 159 { 160 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 161 PetscErrorCode ierr; 162 PetscBool iascii,isdraw; 163 PetscInt i,j; 164 PC_FieldSplitLink ilink = jac->head; 165 MatSchurComplementAinvType atype; 166 167 PetscFunctionBegin; 168 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 169 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 170 if (iascii) { 171 if (jac->bs > 0) { 172 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit with Schur preconditioner, blocksize = %D, factorization %s\n",jac->bs,PCFieldSplitSchurFactTypes[jac->schurfactorization]);CHKERRQ(ierr); 173 } else { 174 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit with Schur preconditioner, factorization %s\n",PCFieldSplitSchurFactTypes[jac->schurfactorization]);CHKERRQ(ierr); 175 } 176 if (pc->useAmat) { 177 ierr = PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for blocks\n");CHKERRQ(ierr); 178 } 179 switch (jac->schurpre) { 180 case PC_FIELDSPLIT_SCHUR_PRE_SELF: 181 ierr = PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from S itself\n");CHKERRQ(ierr); 182 break; 183 case PC_FIELDSPLIT_SCHUR_PRE_SELFP: 184 ierr = MatSchurComplementGetAinvType(jac->schur,&atype);CHKERRQ(ierr); 185 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; 186 case PC_FIELDSPLIT_SCHUR_PRE_A11: 187 ierr = PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from A11\n");CHKERRQ(ierr); 188 break; 189 case PC_FIELDSPLIT_SCHUR_PRE_FULL: 190 ierr = PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from the exact Schur complement\n");CHKERRQ(ierr); 191 break; 192 case PC_FIELDSPLIT_SCHUR_PRE_USER: 193 if (jac->schur_user) { 194 ierr = PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from user provided matrix\n");CHKERRQ(ierr); 195 } else { 196 ierr = PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from A11\n");CHKERRQ(ierr); 197 } 198 break; 199 default: 200 SETERRQ1(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Invalid Schur preconditioning type: %d", jac->schurpre); 201 } 202 ierr = PetscViewerASCIIPrintf(viewer," Split info:\n");CHKERRQ(ierr); 203 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 204 for (i=0; i<jac->nsplits; i++) { 205 if (ilink->fields) { 206 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);CHKERRQ(ierr); 207 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 208 for (j=0; j<ilink->nfields; j++) { 209 if (j > 0) { 210 ierr = PetscViewerASCIIPrintf(viewer,",");CHKERRQ(ierr); 211 } 212 ierr = PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);CHKERRQ(ierr); 213 } 214 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 215 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 216 } else { 217 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i);CHKERRQ(ierr); 218 } 219 ilink = ilink->next; 220 } 221 ierr = PetscViewerASCIIPrintf(viewer,"KSP solver for A00 block\n");CHKERRQ(ierr); 222 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 223 if (jac->head) { 224 ierr = KSPView(jac->head->ksp,viewer);CHKERRQ(ierr); 225 } else {ierr = PetscViewerASCIIPrintf(viewer," not yet available\n");CHKERRQ(ierr);} 226 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 227 if (jac->head && jac->kspupper != jac->head->ksp) { 228 ierr = PetscViewerASCIIPrintf(viewer,"KSP solver for upper A00 in upper triangular factor \n");CHKERRQ(ierr); 229 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 230 if (jac->kspupper) {ierr = KSPView(jac->kspupper,viewer);CHKERRQ(ierr);} 231 else {ierr = PetscViewerASCIIPrintf(viewer," not yet available\n");CHKERRQ(ierr);} 232 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 233 } 234 ierr = PetscViewerASCIIPrintf(viewer,"KSP solver for S = A11 - A10 inv(A00) A01 \n");CHKERRQ(ierr); 235 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 236 if (jac->kspschur) { 237 ierr = KSPView(jac->kspschur,viewer);CHKERRQ(ierr); 238 } else { 239 ierr = PetscViewerASCIIPrintf(viewer," not yet available\n");CHKERRQ(ierr); 240 } 241 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 242 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 243 } else if (isdraw && jac->head) { 244 PetscDraw draw; 245 PetscReal x,y,w,wd,h; 246 PetscInt cnt = 2; 247 char str[32]; 248 249 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 250 ierr = PetscDrawGetCurrentPoint(draw,&x,&y);CHKERRQ(ierr); 251 if (jac->kspupper != jac->head->ksp) cnt++; 252 w = 2*PetscMin(1.0 - x,x); 253 wd = w/(cnt + 1); 254 255 ierr = PetscSNPrintf(str,32,"Schur fact. %s",PCFieldSplitSchurFactTypes[jac->schurfactorization]);CHKERRQ(ierr); 256 ierr = PetscDrawStringBoxed(draw,x,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);CHKERRQ(ierr); 257 y -= h; 258 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_USER && !jac->schur_user) { 259 ierr = PetscSNPrintf(str,32,"Prec. for Schur from %s",PCFieldSplitSchurPreTypes[PC_FIELDSPLIT_SCHUR_PRE_A11]);CHKERRQ(ierr); 260 } else { 261 ierr = PetscSNPrintf(str,32,"Prec. for Schur from %s",PCFieldSplitSchurPreTypes[jac->schurpre]);CHKERRQ(ierr); 262 } 263 ierr = PetscDrawStringBoxed(draw,x+wd*(cnt-1)/2.0,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);CHKERRQ(ierr); 264 y -= h; 265 x = x - wd*(cnt-1)/2.0; 266 267 ierr = PetscDrawPushCurrentPoint(draw,x,y);CHKERRQ(ierr); 268 ierr = KSPView(jac->head->ksp,viewer);CHKERRQ(ierr); 269 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr); 270 if (jac->kspupper != jac->head->ksp) { 271 x += wd; 272 ierr = PetscDrawPushCurrentPoint(draw,x,y);CHKERRQ(ierr); 273 ierr = KSPView(jac->kspupper,viewer);CHKERRQ(ierr); 274 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr); 275 } 276 x += wd; 277 ierr = PetscDrawPushCurrentPoint(draw,x,y);CHKERRQ(ierr); 278 ierr = KSPView(jac->kspschur,viewer);CHKERRQ(ierr); 279 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr); 280 } 281 PetscFunctionReturn(0); 282 } 283 284 static PetscErrorCode PCView_FieldSplit_GKB(PC pc,PetscViewer viewer) 285 { 286 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 287 PetscErrorCode ierr; 288 PetscBool iascii,isdraw; 289 PetscInt i,j; 290 PC_FieldSplitLink ilink = jac->head; 291 292 PetscFunctionBegin; 293 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 294 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 295 if (iascii) { 296 if (jac->bs > 0) { 297 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D, blocksize = %D\n",PCCompositeTypes[jac->type],jac->nsplits,jac->bs);CHKERRQ(ierr); 298 } else { 299 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D\n",PCCompositeTypes[jac->type],jac->nsplits);CHKERRQ(ierr); 300 } 301 if (pc->useAmat) { 302 ierr = PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for blocks\n");CHKERRQ(ierr); 303 } 304 if (jac->diag_use_amat) { 305 ierr = PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for diagonal blocks\n");CHKERRQ(ierr); 306 } 307 if (jac->offdiag_use_amat) { 308 ierr = PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for off-diagonal blocks\n");CHKERRQ(ierr); 309 } 310 311 ierr = PetscViewerASCIIPrintf(viewer," Stopping tolerance=%.1e, delay in error estimate=%D, maximum iterations=%D\n",jac->gkbtol,jac->gkbdelay,jac->gkbmaxit);CHKERRQ(ierr); 312 ierr = PetscViewerASCIIPrintf(viewer," Solver info for H = A00 + nu*A01*A01' matrix:\n");CHKERRQ(ierr); 313 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 314 315 if (ilink->fields) { 316 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",0);CHKERRQ(ierr); 317 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 318 for (j=0; j<ilink->nfields; j++) { 319 if (j > 0) { 320 ierr = PetscViewerASCIIPrintf(viewer,",");CHKERRQ(ierr); 321 } 322 ierr = PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);CHKERRQ(ierr); 323 } 324 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 325 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 326 } else { 327 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",0);CHKERRQ(ierr); 328 } 329 ierr = KSPView(ilink->ksp,viewer);CHKERRQ(ierr); 330 331 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 332 } 333 334 if (isdraw) { 335 PetscDraw draw; 336 PetscReal x,y,w,wd; 337 338 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 339 ierr = PetscDrawGetCurrentPoint(draw,&x,&y);CHKERRQ(ierr); 340 w = 2*PetscMin(1.0 - x,x); 341 wd = w/(jac->nsplits + 1); 342 x = x - wd*(jac->nsplits-1)/2.0; 343 for (i=0; i<jac->nsplits; i++) { 344 ierr = PetscDrawPushCurrentPoint(draw,x,y);CHKERRQ(ierr); 345 ierr = KSPView(ilink->ksp,viewer);CHKERRQ(ierr); 346 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr); 347 x += wd; 348 ilink = ilink->next; 349 } 350 } 351 PetscFunctionReturn(0); 352 } 353 354 355 /* Precondition: jac->bs is set to a meaningful value */ 356 static PetscErrorCode PCFieldSplitSetRuntimeSplits_Private(PC pc) 357 { 358 PetscErrorCode ierr; 359 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 360 PetscInt i,nfields,*ifields,nfields_col,*ifields_col; 361 PetscBool flg,flg_col; 362 char optionname[128],splitname[8],optionname_col[128]; 363 364 PetscFunctionBegin; 365 ierr = PetscMalloc1(jac->bs,&ifields);CHKERRQ(ierr); 366 ierr = PetscMalloc1(jac->bs,&ifields_col);CHKERRQ(ierr); 367 for (i=0,flg=PETSC_TRUE;; i++) { 368 ierr = PetscSNPrintf(splitname,sizeof(splitname),"%D",i);CHKERRQ(ierr); 369 ierr = PetscSNPrintf(optionname,sizeof(optionname),"-pc_fieldsplit_%D_fields",i);CHKERRQ(ierr); 370 ierr = PetscSNPrintf(optionname_col,sizeof(optionname_col),"-pc_fieldsplit_%D_fields_col",i);CHKERRQ(ierr); 371 nfields = jac->bs; 372 nfields_col = jac->bs; 373 ierr = PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix,optionname,ifields,&nfields,&flg);CHKERRQ(ierr); 374 ierr = PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix,optionname_col,ifields_col,&nfields_col,&flg_col);CHKERRQ(ierr); 375 if (!flg) break; 376 else if (flg && !flg_col) { 377 if (!nfields) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot list zero fields"); 378 ierr = PCFieldSplitSetFields(pc,splitname,nfields,ifields,ifields);CHKERRQ(ierr); 379 } else { 380 if (!nfields || !nfields_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot list zero fields"); 381 if (nfields != nfields_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Number of row and column fields must match"); 382 ierr = PCFieldSplitSetFields(pc,splitname,nfields,ifields,ifields_col);CHKERRQ(ierr); 383 } 384 } 385 if (i > 0) { 386 /* Makes command-line setting of splits take precedence over setting them in code. 387 Otherwise subsequent calls to PCFieldSplitSetIS() or PCFieldSplitSetFields() would 388 create new splits, which would probably not be what the user wanted. */ 389 jac->splitdefined = PETSC_TRUE; 390 } 391 ierr = PetscFree(ifields);CHKERRQ(ierr); 392 ierr = PetscFree(ifields_col);CHKERRQ(ierr); 393 PetscFunctionReturn(0); 394 } 395 396 static PetscErrorCode PCFieldSplitSetDefaults(PC pc) 397 { 398 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 399 PetscErrorCode ierr; 400 PC_FieldSplitLink ilink = jac->head; 401 PetscBool fieldsplit_default = PETSC_FALSE,coupling = PETSC_FALSE; 402 PetscInt i; 403 404 PetscFunctionBegin; 405 /* 406 Kinda messy, but at least this now uses DMCreateFieldDecomposition(). 407 Should probably be rewritten. 408 */ 409 if (!ilink) { 410 ierr = PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_coupling",&coupling,NULL);CHKERRQ(ierr); 411 if (pc->dm && jac->dm_splits && !jac->detect && !coupling) { 412 PetscInt numFields, f, i, j; 413 char **fieldNames; 414 IS *fields; 415 DM *dms; 416 DM subdm[128]; 417 PetscBool flg; 418 419 ierr = DMCreateFieldDecomposition(pc->dm, &numFields, &fieldNames, &fields, &dms);CHKERRQ(ierr); 420 /* Allow the user to prescribe the splits */ 421 for (i = 0, flg = PETSC_TRUE;; i++) { 422 PetscInt ifields[128]; 423 IS compField; 424 char optionname[128], splitname[8]; 425 PetscInt nfields = numFields; 426 427 ierr = PetscSNPrintf(optionname, sizeof(optionname), "-pc_fieldsplit_%D_fields", i);CHKERRQ(ierr); 428 ierr = PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix, optionname, ifields, &nfields, &flg);CHKERRQ(ierr); 429 if (!flg) break; 430 if (numFields > 128) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot currently support %d > 128 fields", numFields); 431 ierr = DMCreateSubDM(pc->dm, nfields, ifields, &compField, &subdm[i]);CHKERRQ(ierr); 432 if (nfields == 1) { 433 ierr = PCFieldSplitSetIS(pc, fieldNames[ifields[0]], compField);CHKERRQ(ierr); 434 } else { 435 ierr = PetscSNPrintf(splitname, sizeof(splitname), "%D", i);CHKERRQ(ierr); 436 ierr = PCFieldSplitSetIS(pc, splitname, compField);CHKERRQ(ierr); 437 } 438 ierr = ISDestroy(&compField);CHKERRQ(ierr); 439 for (j = 0; j < nfields; ++j) { 440 f = ifields[j]; 441 ierr = PetscFree(fieldNames[f]);CHKERRQ(ierr); 442 ierr = ISDestroy(&fields[f]);CHKERRQ(ierr); 443 } 444 } 445 if (i == 0) { 446 for (f = 0; f < numFields; ++f) { 447 ierr = PCFieldSplitSetIS(pc, fieldNames[f], fields[f]);CHKERRQ(ierr); 448 ierr = PetscFree(fieldNames[f]);CHKERRQ(ierr); 449 ierr = ISDestroy(&fields[f]);CHKERRQ(ierr); 450 } 451 } else { 452 for (j=0; j<numFields; j++) { 453 ierr = DMDestroy(dms+j);CHKERRQ(ierr); 454 } 455 ierr = PetscFree(dms);CHKERRQ(ierr); 456 ierr = PetscMalloc1(i, &dms);CHKERRQ(ierr); 457 for (j = 0; j < i; ++j) dms[j] = subdm[j]; 458 } 459 ierr = PetscFree(fieldNames);CHKERRQ(ierr); 460 ierr = PetscFree(fields);CHKERRQ(ierr); 461 if (dms) { 462 ierr = PetscInfo(pc, "Setting up physics based fieldsplit preconditioner using the embedded DM\n");CHKERRQ(ierr); 463 for (ilink = jac->head, i = 0; ilink; ilink = ilink->next, ++i) { 464 const char *prefix; 465 ierr = PetscObjectGetOptionsPrefix((PetscObject)(ilink->ksp),&prefix);CHKERRQ(ierr); 466 ierr = PetscObjectSetOptionsPrefix((PetscObject)(dms[i]), prefix);CHKERRQ(ierr); 467 ierr = KSPSetDM(ilink->ksp, dms[i]);CHKERRQ(ierr); 468 ierr = KSPSetDMActive(ilink->ksp, PETSC_FALSE);CHKERRQ(ierr); 469 { 470 PetscErrorCode (*func)(KSP,Mat,Mat,void*); 471 void *ctx; 472 473 ierr = DMKSPGetComputeOperators(pc->dm, &func, &ctx);CHKERRQ(ierr); 474 ierr = DMKSPSetComputeOperators(dms[i], func, ctx);CHKERRQ(ierr); 475 } 476 ierr = PetscObjectIncrementTabLevel((PetscObject)dms[i],(PetscObject)ilink->ksp,0);CHKERRQ(ierr); 477 ierr = DMDestroy(&dms[i]);CHKERRQ(ierr); 478 } 479 ierr = PetscFree(dms);CHKERRQ(ierr); 480 } 481 } else { 482 if (jac->bs <= 0) { 483 if (pc->pmat) { 484 ierr = MatGetBlockSize(pc->pmat,&jac->bs);CHKERRQ(ierr); 485 } else jac->bs = 1; 486 } 487 488 if (jac->detect) { 489 IS zerodiags,rest; 490 PetscInt nmin,nmax; 491 492 ierr = MatGetOwnershipRange(pc->mat,&nmin,&nmax);CHKERRQ(ierr); 493 ierr = MatFindZeroDiagonals(pc->mat,&zerodiags);CHKERRQ(ierr); 494 ierr = ISComplement(zerodiags,nmin,nmax,&rest);CHKERRQ(ierr); 495 ierr = PCFieldSplitSetIS(pc,"0",rest);CHKERRQ(ierr); 496 ierr = PCFieldSplitSetIS(pc,"1",zerodiags);CHKERRQ(ierr); 497 ierr = ISDestroy(&zerodiags);CHKERRQ(ierr); 498 ierr = ISDestroy(&rest);CHKERRQ(ierr); 499 } else if (coupling) { 500 IS coupling,rest; 501 PetscInt nmin,nmax; 502 503 ierr = MatGetOwnershipRange(pc->mat,&nmin,&nmax);CHKERRQ(ierr); 504 ierr = MatFindOffBlockDiagonalEntries(pc->mat,&coupling);CHKERRQ(ierr); 505 ierr = ISCreateStride(PetscObjectComm((PetscObject)pc->mat),nmax-nmin,nmin,1,&rest);CHKERRQ(ierr); 506 ierr = ISSetIdentity(rest);CHKERRQ(ierr); 507 ierr = PCFieldSplitSetIS(pc,"0",rest);CHKERRQ(ierr); 508 ierr = PCFieldSplitSetIS(pc,"1",coupling);CHKERRQ(ierr); 509 ierr = ISDestroy(&coupling);CHKERRQ(ierr); 510 ierr = ISDestroy(&rest);CHKERRQ(ierr); 511 } else { 512 ierr = PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_default",&fieldsplit_default,NULL);CHKERRQ(ierr); 513 if (!fieldsplit_default) { 514 /* Allow user to set fields from command line, if bs was known at the time of PCSetFromOptions_FieldSplit() 515 then it is set there. This is not ideal because we should only have options set in XXSetFromOptions(). */ 516 ierr = PCFieldSplitSetRuntimeSplits_Private(pc);CHKERRQ(ierr); 517 if (jac->splitdefined) {ierr = PetscInfo(pc,"Splits defined using the options database\n");CHKERRQ(ierr);} 518 } 519 if ((fieldsplit_default || !jac->splitdefined) && !jac->isrestrict) { 520 Mat M = pc->pmat; 521 PetscBool isnest; 522 523 ierr = PetscInfo(pc,"Using default splitting of fields\n");CHKERRQ(ierr); 524 ierr = PetscObjectTypeCompare((PetscObject)pc->pmat,MATNEST,&isnest);CHKERRQ(ierr); 525 if (!isnest) { 526 M = pc->mat; 527 ierr = PetscObjectTypeCompare((PetscObject)pc->mat,MATNEST,&isnest);CHKERRQ(ierr); 528 } 529 if (isnest) { 530 IS *fields; 531 PetscInt nf; 532 533 ierr = MatNestGetSize(M,&nf,NULL);CHKERRQ(ierr); 534 ierr = PetscMalloc1(nf,&fields);CHKERRQ(ierr); 535 ierr = MatNestGetISs(M,fields,NULL);CHKERRQ(ierr); 536 for (i=0;i<nf;i++) { 537 ierr = PCFieldSplitSetIS(pc,NULL,fields[i]);CHKERRQ(ierr); 538 } 539 ierr = PetscFree(fields);CHKERRQ(ierr); 540 } else { 541 for (i=0; i<jac->bs; i++) { 542 char splitname[8]; 543 ierr = PetscSNPrintf(splitname,sizeof(splitname),"%D",i);CHKERRQ(ierr); 544 ierr = PCFieldSplitSetFields(pc,splitname,1,&i,&i);CHKERRQ(ierr); 545 } 546 jac->defaultsplit = PETSC_TRUE; 547 } 548 } 549 } 550 } 551 } else if (jac->nsplits == 1) { 552 if (ilink->is) { 553 IS is2; 554 PetscInt nmin,nmax; 555 556 ierr = MatGetOwnershipRange(pc->mat,&nmin,&nmax);CHKERRQ(ierr); 557 ierr = ISComplement(ilink->is,nmin,nmax,&is2);CHKERRQ(ierr); 558 ierr = PCFieldSplitSetIS(pc,"1",is2);CHKERRQ(ierr); 559 ierr = ISDestroy(&is2);CHKERRQ(ierr); 560 } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Must provide at least two sets of fields to PCFieldSplit()"); 561 } 562 563 if (jac->nsplits < 2) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Unhandled case, must have at least two fields, not %d", jac->nsplits); 564 PetscFunctionReturn(0); 565 } 566 567 static PetscErrorCode MatGolubKahanComputeExplicitOperator(Mat A,Mat B,Mat C,Mat *H,PetscReal gkbnu) 568 { 569 PetscErrorCode ierr; 570 Mat BT,T; 571 PetscReal nrmT,nrmB; 572 573 PetscFunctionBegin; 574 ierr = MatHermitianTranspose(C,MAT_INITIAL_MATRIX,&T);CHKERRQ(ierr); /* Test if augmented matrix is symmetric */ 575 ierr = MatAXPY(T,-1.0,B,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 576 ierr = MatNorm(T,NORM_1,&nrmT);CHKERRQ(ierr); 577 ierr = MatNorm(B,NORM_1,&nrmB);CHKERRQ(ierr); 578 if (nrmB > 0) { 579 if (nrmT/nrmB >= PETSC_SMALL) { 580 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Matrix is not symmetric/hermitian, GKB is not applicable."); 581 } 582 } 583 /* Compute augmented Lagrangian matrix H = A00 + nu*A01*A01'. This corresponds to */ 584 /* setting N := 1/nu*I in [Ar13]. */ 585 ierr = MatHermitianTranspose(B,MAT_INITIAL_MATRIX,&BT);CHKERRQ(ierr); 586 ierr = MatMatMult(B,BT,MAT_INITIAL_MATRIX,PETSC_DEFAULT,H);CHKERRQ(ierr); /* H = A01*A01' */ 587 ierr = MatAYPX(*H,gkbnu,A,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); /* H = A00 + nu*A01*A01' */ 588 589 ierr = MatDestroy(&BT);CHKERRQ(ierr); 590 ierr = MatDestroy(&T);CHKERRQ(ierr); 591 PetscFunctionReturn(0); 592 } 593 594 PETSC_EXTERN PetscErrorCode PetscOptionsFindPairPrefix_Private(PetscOptions,const char pre[], const char name[],const char *value[],PetscBool *flg); 595 596 static PetscErrorCode PCSetUp_FieldSplit(PC pc) 597 { 598 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 599 PetscErrorCode ierr; 600 PC_FieldSplitLink ilink; 601 PetscInt i,nsplit; 602 PetscBool sorted, sorted_col; 603 604 PetscFunctionBegin; 605 pc->failedreason = PC_NOERROR; 606 ierr = PCFieldSplitSetDefaults(pc);CHKERRQ(ierr); 607 nsplit = jac->nsplits; 608 ilink = jac->head; 609 610 /* get the matrices for each split */ 611 if (!jac->issetup) { 612 PetscInt rstart,rend,nslots,bs; 613 614 jac->issetup = PETSC_TRUE; 615 616 /* This is done here instead of in PCFieldSplitSetFields() because may not have matrix at that point */ 617 if (jac->defaultsplit || !ilink->is) { 618 if (jac->bs <= 0) jac->bs = nsplit; 619 } 620 bs = jac->bs; 621 ierr = MatGetOwnershipRange(pc->pmat,&rstart,&rend);CHKERRQ(ierr); 622 nslots = (rend - rstart)/bs; 623 for (i=0; i<nsplit; i++) { 624 if (jac->defaultsplit) { 625 ierr = ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+i,nsplit,&ilink->is);CHKERRQ(ierr); 626 ierr = ISDuplicate(ilink->is,&ilink->is_col);CHKERRQ(ierr); 627 } else if (!ilink->is) { 628 if (ilink->nfields > 1) { 629 PetscInt *ii,*jj,j,k,nfields = ilink->nfields,*fields = ilink->fields,*fields_col = ilink->fields_col; 630 ierr = PetscMalloc1(ilink->nfields*nslots,&ii);CHKERRQ(ierr); 631 ierr = PetscMalloc1(ilink->nfields*nslots,&jj);CHKERRQ(ierr); 632 for (j=0; j<nslots; j++) { 633 for (k=0; k<nfields; k++) { 634 ii[nfields*j + k] = rstart + bs*j + fields[k]; 635 jj[nfields*j + k] = rstart + bs*j + fields_col[k]; 636 } 637 } 638 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,ii,PETSC_OWN_POINTER,&ilink->is);CHKERRQ(ierr); 639 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,jj,PETSC_OWN_POINTER,&ilink->is_col);CHKERRQ(ierr); 640 ierr = ISSetBlockSize(ilink->is, nfields);CHKERRQ(ierr); 641 ierr = ISSetBlockSize(ilink->is_col, nfields);CHKERRQ(ierr); 642 } else { 643 ierr = ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields[0],bs,&ilink->is);CHKERRQ(ierr); 644 ierr = ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields_col[0],bs,&ilink->is_col);CHKERRQ(ierr); 645 } 646 } 647 ierr = ISSorted(ilink->is,&sorted);CHKERRQ(ierr); 648 if (ilink->is_col) { ierr = ISSorted(ilink->is_col,&sorted_col);CHKERRQ(ierr); } 649 if (!sorted || !sorted_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Fields must be sorted when creating split"); 650 ilink = ilink->next; 651 } 652 } 653 654 ilink = jac->head; 655 if (!jac->pmat) { 656 Vec xtmp; 657 658 ierr = MatCreateVecs(pc->pmat,&xtmp,NULL);CHKERRQ(ierr); 659 ierr = PetscMalloc1(nsplit,&jac->pmat);CHKERRQ(ierr); 660 ierr = PetscMalloc2(nsplit,&jac->x,nsplit,&jac->y);CHKERRQ(ierr); 661 for (i=0; i<nsplit; i++) { 662 MatNullSpace sp; 663 664 /* Check for preconditioning matrix attached to IS */ 665 ierr = PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &jac->pmat[i]);CHKERRQ(ierr); 666 if (jac->pmat[i]) { 667 ierr = PetscObjectReference((PetscObject) jac->pmat[i]);CHKERRQ(ierr); 668 if (jac->type == PC_COMPOSITE_SCHUR) { 669 jac->schur_user = jac->pmat[i]; 670 671 ierr = PetscObjectReference((PetscObject) jac->schur_user);CHKERRQ(ierr); 672 } 673 } else { 674 const char *prefix; 675 ierr = MatCreateSubMatrix(pc->pmat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->pmat[i]);CHKERRQ(ierr); 676 ierr = KSPGetOptionsPrefix(ilink->ksp,&prefix);CHKERRQ(ierr); 677 ierr = MatSetOptionsPrefix(jac->pmat[i],prefix);CHKERRQ(ierr); 678 ierr = MatViewFromOptions(jac->pmat[i],NULL,"-mat_view");CHKERRQ(ierr); 679 } 680 /* create work vectors for each split */ 681 ierr = MatCreateVecs(jac->pmat[i],&jac->x[i],&jac->y[i]);CHKERRQ(ierr); 682 ilink->x = jac->x[i]; ilink->y = jac->y[i]; ilink->z = NULL; 683 /* compute scatter contexts needed by multiplicative versions and non-default splits */ 684 ierr = VecScatterCreate(xtmp,ilink->is,jac->x[i],NULL,&ilink->sctx);CHKERRQ(ierr); 685 ierr = PetscObjectQuery((PetscObject) ilink->is, "nearnullspace", (PetscObject*) &sp);CHKERRQ(ierr); 686 if (sp) { 687 ierr = MatSetNearNullSpace(jac->pmat[i], sp);CHKERRQ(ierr); 688 } 689 ilink = ilink->next; 690 } 691 ierr = VecDestroy(&xtmp);CHKERRQ(ierr); 692 } else { 693 MatReuse scall; 694 if (pc->flag == DIFFERENT_NONZERO_PATTERN) { 695 for (i=0; i<nsplit; i++) { 696 ierr = MatDestroy(&jac->pmat[i]);CHKERRQ(ierr); 697 } 698 scall = MAT_INITIAL_MATRIX; 699 } else scall = MAT_REUSE_MATRIX; 700 701 for (i=0; i<nsplit; i++) { 702 Mat pmat; 703 704 /* Check for preconditioning matrix attached to IS */ 705 ierr = PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &pmat);CHKERRQ(ierr); 706 if (!pmat) { 707 ierr = MatCreateSubMatrix(pc->pmat,ilink->is,ilink->is_col,scall,&jac->pmat[i]);CHKERRQ(ierr); 708 } 709 ilink = ilink->next; 710 } 711 } 712 if (jac->diag_use_amat) { 713 ilink = jac->head; 714 if (!jac->mat) { 715 ierr = PetscMalloc1(nsplit,&jac->mat);CHKERRQ(ierr); 716 for (i=0; i<nsplit; i++) { 717 ierr = MatCreateSubMatrix(pc->mat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->mat[i]);CHKERRQ(ierr); 718 ilink = ilink->next; 719 } 720 } else { 721 MatReuse scall; 722 if (pc->flag == DIFFERENT_NONZERO_PATTERN) { 723 for (i=0; i<nsplit; i++) { 724 ierr = MatDestroy(&jac->mat[i]);CHKERRQ(ierr); 725 } 726 scall = MAT_INITIAL_MATRIX; 727 } else scall = MAT_REUSE_MATRIX; 728 729 for (i=0; i<nsplit; i++) { 730 ierr = MatCreateSubMatrix(pc->mat,ilink->is,ilink->is_col,scall,&jac->mat[i]);CHKERRQ(ierr); 731 ilink = ilink->next; 732 } 733 } 734 } else { 735 jac->mat = jac->pmat; 736 } 737 738 /* Check for null space attached to IS */ 739 ilink = jac->head; 740 for (i=0; i<nsplit; i++) { 741 MatNullSpace sp; 742 743 ierr = PetscObjectQuery((PetscObject) ilink->is, "nullspace", (PetscObject*) &sp);CHKERRQ(ierr); 744 if (sp) { 745 ierr = MatSetNullSpace(jac->mat[i], sp);CHKERRQ(ierr); 746 } 747 ilink = ilink->next; 748 } 749 750 if (jac->type != PC_COMPOSITE_ADDITIVE && jac->type != PC_COMPOSITE_SCHUR && jac->type != PC_COMPOSITE_GKB) { 751 /* extract the rows of the matrix associated with each field: used for efficient computation of residual inside algorithm */ 752 /* FIXME: Can/should we reuse jac->mat whenever (jac->diag_use_amat) is true? */ 753 ilink = jac->head; 754 if (nsplit == 2 && jac->type == PC_COMPOSITE_MULTIPLICATIVE) { 755 /* special case need where Afield[0] is not needed and only certain columns of Afield[1] are needed since update is only on those rows of the solution */ 756 if (!jac->Afield) { 757 ierr = PetscCalloc1(nsplit,&jac->Afield);CHKERRQ(ierr); 758 if (jac->offdiag_use_amat) { 759 ierr = MatCreateSubMatrix(pc->mat,ilink->next->is,ilink->is,MAT_INITIAL_MATRIX,&jac->Afield[1]);CHKERRQ(ierr); 760 } else { 761 ierr = MatCreateSubMatrix(pc->pmat,ilink->next->is,ilink->is,MAT_INITIAL_MATRIX,&jac->Afield[1]);CHKERRQ(ierr); 762 } 763 } else { 764 MatReuse scall; 765 766 if (pc->flag == DIFFERENT_NONZERO_PATTERN) { 767 ierr = MatDestroy(&jac->Afield[1]);CHKERRQ(ierr); 768 scall = MAT_INITIAL_MATRIX; 769 } else scall = MAT_REUSE_MATRIX; 770 771 if (jac->offdiag_use_amat) { 772 ierr = MatCreateSubMatrix(pc->mat,ilink->next->is,ilink->is,scall,&jac->Afield[1]);CHKERRQ(ierr); 773 } else { 774 ierr = MatCreateSubMatrix(pc->pmat,ilink->next->is,ilink->is,scall,&jac->Afield[1]);CHKERRQ(ierr); 775 } 776 } 777 } else { 778 if (!jac->Afield) { 779 ierr = PetscMalloc1(nsplit,&jac->Afield);CHKERRQ(ierr); 780 for (i=0; i<nsplit; i++) { 781 if (jac->offdiag_use_amat) { 782 ierr = MatCreateSubMatrix(pc->mat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);CHKERRQ(ierr); 783 } else { 784 ierr = MatCreateSubMatrix(pc->pmat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);CHKERRQ(ierr); 785 } 786 ilink = ilink->next; 787 } 788 } else { 789 MatReuse scall; 790 if (pc->flag == DIFFERENT_NONZERO_PATTERN) { 791 for (i=0; i<nsplit; i++) { 792 ierr = MatDestroy(&jac->Afield[i]);CHKERRQ(ierr); 793 } 794 scall = MAT_INITIAL_MATRIX; 795 } else scall = MAT_REUSE_MATRIX; 796 797 for (i=0; i<nsplit; i++) { 798 if (jac->offdiag_use_amat) { 799 ierr = MatCreateSubMatrix(pc->mat,ilink->is,NULL,scall,&jac->Afield[i]);CHKERRQ(ierr); 800 } else { 801 ierr = MatCreateSubMatrix(pc->pmat,ilink->is,NULL,scall,&jac->Afield[i]);CHKERRQ(ierr); 802 } 803 ilink = ilink->next; 804 } 805 } 806 } 807 } 808 809 if (jac->type == PC_COMPOSITE_SCHUR) { 810 IS ccis; 811 PetscBool isspd; 812 PetscInt rstart,rend; 813 char lscname[256]; 814 PetscObject LSC_L; 815 816 if (nsplit != 2) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"To use Schur complement preconditioner you must have exactly 2 fields"); 817 818 /* If pc->mat is SPD, don't scale by -1 the Schur complement */ 819 if (jac->schurscale == (PetscScalar)-1.0) { 820 ierr = MatGetOption(pc->pmat,MAT_SPD,&isspd);CHKERRQ(ierr); 821 jac->schurscale = (isspd == PETSC_TRUE) ? 1.0 : -1.0; 822 } 823 824 /* When extracting off-diagonal submatrices, we take complements from this range */ 825 ierr = MatGetOwnershipRangeColumn(pc->mat,&rstart,&rend);CHKERRQ(ierr); 826 827 if (jac->schur) { 828 KSP kspA = jac->head->ksp, kspInner = NULL, kspUpper = jac->kspupper; 829 MatReuse scall; 830 831 if (pc->flag == DIFFERENT_NONZERO_PATTERN) { 832 scall = MAT_INITIAL_MATRIX; 833 ierr = MatDestroy(&jac->B);CHKERRQ(ierr); 834 ierr = MatDestroy(&jac->C);CHKERRQ(ierr); 835 } else scall = MAT_REUSE_MATRIX; 836 837 ierr = MatSchurComplementGetKSP(jac->schur, &kspInner);CHKERRQ(ierr); 838 ilink = jac->head; 839 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 840 if (jac->offdiag_use_amat) { 841 ierr = MatCreateSubMatrix(pc->mat,ilink->is,ccis,scall,&jac->B);CHKERRQ(ierr); 842 } else { 843 ierr = MatCreateSubMatrix(pc->pmat,ilink->is,ccis,scall,&jac->B);CHKERRQ(ierr); 844 } 845 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 846 ilink = ilink->next; 847 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 848 if (jac->offdiag_use_amat) { 849 ierr = MatCreateSubMatrix(pc->mat,ilink->is,ccis,scall,&jac->C);CHKERRQ(ierr); 850 } else { 851 ierr = MatCreateSubMatrix(pc->pmat,ilink->is,ccis,scall,&jac->C);CHKERRQ(ierr); 852 } 853 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 854 ierr = MatSchurComplementUpdateSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);CHKERRQ(ierr); 855 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) { 856 ierr = MatDestroy(&jac->schurp);CHKERRQ(ierr); 857 ierr = MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);CHKERRQ(ierr); 858 } 859 if (kspA != kspInner) { 860 ierr = KSPSetOperators(kspA,jac->mat[0],jac->pmat[0]);CHKERRQ(ierr); 861 } 862 if (kspUpper != kspA) { 863 ierr = KSPSetOperators(kspUpper,jac->mat[0],jac->pmat[0]);CHKERRQ(ierr); 864 } 865 ierr = KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac));CHKERRQ(ierr); 866 } else { 867 const char *Dprefix; 868 char schurprefix[256], schurmatprefix[256]; 869 char schurtestoption[256]; 870 MatNullSpace sp; 871 PetscBool flg; 872 KSP kspt; 873 874 /* extract the A01 and A10 matrices */ 875 ilink = jac->head; 876 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 877 if (jac->offdiag_use_amat) { 878 ierr = MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);CHKERRQ(ierr); 879 } else { 880 ierr = MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);CHKERRQ(ierr); 881 } 882 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 883 ilink = ilink->next; 884 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 885 if (jac->offdiag_use_amat) { 886 ierr = MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);CHKERRQ(ierr); 887 } else { 888 ierr = MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);CHKERRQ(ierr); 889 } 890 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 891 892 /* Use mat[0] (diagonal block of Amat) preconditioned by pmat[0] to define Schur complement */ 893 ierr = MatCreate(((PetscObject)jac->mat[0])->comm,&jac->schur);CHKERRQ(ierr); 894 ierr = MatSetType(jac->schur,MATSCHURCOMPLEMENT);CHKERRQ(ierr); 895 ierr = MatSchurComplementSetSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);CHKERRQ(ierr); 896 ierr = PetscSNPrintf(schurmatprefix, sizeof(schurmatprefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);CHKERRQ(ierr); 897 ierr = MatSetOptionsPrefix(jac->schur,schurmatprefix);CHKERRQ(ierr); 898 ierr = MatSchurComplementGetKSP(jac->schur,&kspt);CHKERRQ(ierr); 899 ierr = KSPSetOptionsPrefix(kspt,schurmatprefix);CHKERRQ(ierr); 900 901 /* Note: this is not true in general */ 902 ierr = MatGetNullSpace(jac->mat[1], &sp);CHKERRQ(ierr); 903 if (sp) { 904 ierr = MatSetNullSpace(jac->schur, sp);CHKERRQ(ierr); 905 } 906 907 ierr = PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_inner_", ilink->splitname);CHKERRQ(ierr); 908 ierr = PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options,((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);CHKERRQ(ierr); 909 if (flg) { 910 DM dmInner; 911 KSP kspInner; 912 PC pcInner; 913 914 ierr = MatSchurComplementGetKSP(jac->schur, &kspInner);CHKERRQ(ierr); 915 ierr = KSPReset(kspInner);CHKERRQ(ierr); 916 ierr = KSPSetOperators(kspInner,jac->mat[0],jac->pmat[0]);CHKERRQ(ierr); 917 ierr = PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_inner_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);CHKERRQ(ierr); 918 /* Indent this deeper to emphasize the "inner" nature of this solver. */ 919 ierr = PetscObjectIncrementTabLevel((PetscObject)kspInner, (PetscObject) pc, 2);CHKERRQ(ierr); 920 ierr = PetscObjectIncrementTabLevel((PetscObject)kspInner->pc, (PetscObject) pc, 2);CHKERRQ(ierr); 921 ierr = KSPSetOptionsPrefix(kspInner, schurprefix);CHKERRQ(ierr); 922 923 /* Set DM for new solver */ 924 ierr = KSPGetDM(jac->head->ksp, &dmInner);CHKERRQ(ierr); 925 ierr = KSPSetDM(kspInner, dmInner);CHKERRQ(ierr); 926 ierr = KSPSetDMActive(kspInner, PETSC_FALSE);CHKERRQ(ierr); 927 928 /* Defaults to PCKSP as preconditioner */ 929 ierr = KSPGetPC(kspInner, &pcInner);CHKERRQ(ierr); 930 ierr = PCSetType(pcInner, PCKSP);CHKERRQ(ierr); 931 ierr = PCKSPSetKSP(pcInner, jac->head->ksp);CHKERRQ(ierr); 932 } else { 933 /* Use the outer solver for the inner solve, but revert the KSPPREONLY from PCFieldSplitSetFields_FieldSplit or 934 * PCFieldSplitSetIS_FieldSplit. We don't want KSPPREONLY because it makes the Schur complement inexact, 935 * preventing Schur complement reduction to be an accurate solve. Usually when an iterative solver is used for 936 * S = D - C A_inner^{-1} B, we expect S to be defined using an accurate definition of A_inner^{-1}, so we make 937 * GMRES the default. Note that it is also common to use PREONLY for S, in which case S may not be used 938 * directly, and the user is responsible for setting an inexact method for fieldsplit's A^{-1}. */ 939 ierr = KSPSetType(jac->head->ksp,KSPGMRES);CHKERRQ(ierr); 940 ierr = MatSchurComplementSetKSP(jac->schur,jac->head->ksp);CHKERRQ(ierr); 941 } 942 ierr = KSPSetOperators(jac->head->ksp,jac->mat[0],jac->pmat[0]);CHKERRQ(ierr); 943 ierr = KSPSetFromOptions(jac->head->ksp);CHKERRQ(ierr); 944 ierr = MatSetFromOptions(jac->schur);CHKERRQ(ierr); 945 946 ierr = PetscObjectTypeCompare((PetscObject)jac->schur, MATSCHURCOMPLEMENT, &flg);CHKERRQ(ierr); 947 if (flg) { /* Need to do this otherwise PCSetUp_KSP will overwrite the amat of jac->head->ksp */ 948 KSP kspInner; 949 PC pcInner; 950 951 ierr = MatSchurComplementGetKSP(jac->schur, &kspInner);CHKERRQ(ierr); 952 ierr = KSPGetPC(kspInner, &pcInner);CHKERRQ(ierr); 953 ierr = PetscObjectTypeCompare((PetscObject)pcInner, PCKSP, &flg);CHKERRQ(ierr); 954 if (flg) { 955 KSP ksp; 956 957 ierr = PCKSPGetKSP(pcInner, &ksp);CHKERRQ(ierr); 958 if (ksp == jac->head->ksp) { 959 ierr = PCSetUseAmat(pcInner, PETSC_TRUE);CHKERRQ(ierr); 960 } 961 } 962 } 963 ierr = PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_upper_", ilink->splitname);CHKERRQ(ierr); 964 ierr = PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options,((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);CHKERRQ(ierr); 965 if (flg) { 966 DM dmInner; 967 968 ierr = PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_upper_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);CHKERRQ(ierr); 969 ierr = KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspupper);CHKERRQ(ierr); 970 ierr = KSPSetErrorIfNotConverged(jac->kspupper,pc->erroriffailure);CHKERRQ(ierr); 971 ierr = KSPSetOptionsPrefix(jac->kspupper, schurprefix);CHKERRQ(ierr); 972 ierr = PetscObjectIncrementTabLevel((PetscObject)jac->kspupper, (PetscObject) pc, 1);CHKERRQ(ierr); 973 ierr = PetscObjectIncrementTabLevel((PetscObject)jac->kspupper->pc, (PetscObject) pc, 1);CHKERRQ(ierr); 974 ierr = KSPGetDM(jac->head->ksp, &dmInner);CHKERRQ(ierr); 975 ierr = KSPSetDM(jac->kspupper, dmInner);CHKERRQ(ierr); 976 ierr = KSPSetDMActive(jac->kspupper, PETSC_FALSE);CHKERRQ(ierr); 977 ierr = KSPSetFromOptions(jac->kspupper);CHKERRQ(ierr); 978 ierr = KSPSetOperators(jac->kspupper,jac->mat[0],jac->pmat[0]);CHKERRQ(ierr); 979 ierr = VecDuplicate(jac->head->x, &jac->head->z);CHKERRQ(ierr); 980 } else { 981 jac->kspupper = jac->head->ksp; 982 ierr = PetscObjectReference((PetscObject) jac->head->ksp);CHKERRQ(ierr); 983 } 984 985 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) { 986 ierr = MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);CHKERRQ(ierr); 987 } 988 ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&jac->kspschur);CHKERRQ(ierr); 989 ierr = KSPSetErrorIfNotConverged(jac->kspschur,pc->erroriffailure);CHKERRQ(ierr); 990 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->kspschur);CHKERRQ(ierr); 991 ierr = PetscObjectIncrementTabLevel((PetscObject)jac->kspschur,(PetscObject)pc,1);CHKERRQ(ierr); 992 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) { 993 PC pcschur; 994 ierr = KSPGetPC(jac->kspschur,&pcschur);CHKERRQ(ierr); 995 ierr = PCSetType(pcschur,PCNONE);CHKERRQ(ierr); 996 /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */ 997 } else if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_FULL) { 998 ierr = MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user);CHKERRQ(ierr); 999 } 1000 ierr = KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac));CHKERRQ(ierr); 1001 ierr = KSPGetOptionsPrefix(jac->head->next->ksp, &Dprefix);CHKERRQ(ierr); 1002 ierr = KSPSetOptionsPrefix(jac->kspschur, Dprefix);CHKERRQ(ierr); 1003 /* propagate DM */ 1004 { 1005 DM sdm; 1006 ierr = KSPGetDM(jac->head->next->ksp, &sdm);CHKERRQ(ierr); 1007 if (sdm) { 1008 ierr = KSPSetDM(jac->kspschur, sdm);CHKERRQ(ierr); 1009 ierr = KSPSetDMActive(jac->kspschur, PETSC_FALSE);CHKERRQ(ierr); 1010 } 1011 } 1012 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 1013 /* need to call this every time, since the jac->kspschur is freshly created, otherwise its options never get set */ 1014 ierr = KSPSetFromOptions(jac->kspschur);CHKERRQ(ierr); 1015 } 1016 ierr = MatAssemblyBegin(jac->schur,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1017 ierr = MatAssemblyEnd(jac->schur,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1018 1019 /* HACK: special support to forward L and Lp matrices that might be used by PCLSC */ 1020 ierr = PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_L",ilink->splitname);CHKERRQ(ierr); 1021 ierr = PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr); 1022 if (!LSC_L) {ierr = PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr);} 1023 if (LSC_L) {ierr = PetscObjectCompose((PetscObject)jac->schur,"LSC_L",(PetscObject)LSC_L);CHKERRQ(ierr);} 1024 ierr = PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_Lp",ilink->splitname);CHKERRQ(ierr); 1025 ierr = PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr); 1026 if (!LSC_L) {ierr = PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr);} 1027 if (LSC_L) {ierr = PetscObjectCompose((PetscObject)jac->schur,"LSC_Lp",(PetscObject)LSC_L);CHKERRQ(ierr);} 1028 } else if (jac->type == PC_COMPOSITE_GKB) { 1029 IS ccis; 1030 PetscInt rstart,rend; 1031 1032 if (nsplit != 2) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"To use GKB preconditioner you must have exactly 2 fields"); 1033 1034 ilink = jac->head; 1035 1036 /* When extracting off-diagonal submatrices, we take complements from this range */ 1037 ierr = MatGetOwnershipRangeColumn(pc->mat,&rstart,&rend);CHKERRQ(ierr); 1038 1039 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 1040 if (jac->offdiag_use_amat) { 1041 ierr = MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);CHKERRQ(ierr); 1042 } else { 1043 ierr = MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);CHKERRQ(ierr); 1044 } 1045 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 1046 /* Create work vectors for GKB algorithm */ 1047 ierr = VecDuplicate(ilink->x,&jac->u);CHKERRQ(ierr); 1048 ierr = VecDuplicate(ilink->x,&jac->Hu);CHKERRQ(ierr); 1049 ierr = VecDuplicate(ilink->x,&jac->w2);CHKERRQ(ierr); 1050 ilink = ilink->next; 1051 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 1052 if (jac->offdiag_use_amat) { 1053 ierr = MatCreateSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);CHKERRQ(ierr); 1054 } else { 1055 ierr = MatCreateSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);CHKERRQ(ierr); 1056 } 1057 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 1058 /* Create work vectors for GKB algorithm */ 1059 ierr = VecDuplicate(ilink->x,&jac->v);CHKERRQ(ierr); 1060 ierr = VecDuplicate(ilink->x,&jac->d);CHKERRQ(ierr); 1061 ierr = VecDuplicate(ilink->x,&jac->w1);CHKERRQ(ierr); 1062 ierr = MatGolubKahanComputeExplicitOperator(jac->mat[0],jac->B,jac->C,&jac->H,jac->gkbnu);CHKERRQ(ierr); 1063 ierr = PetscCalloc1(jac->gkbdelay,&jac->vecz);CHKERRQ(ierr); 1064 1065 ilink = jac->head; 1066 ierr = KSPSetOperators(ilink->ksp,jac->H,jac->H);CHKERRQ(ierr); 1067 if (!jac->suboptionsset) {ierr = KSPSetFromOptions(ilink->ksp);CHKERRQ(ierr);} 1068 /* Create gkb_monitor context */ 1069 if (jac->gkbmonitor) { 1070 PetscInt tablevel; 1071 ierr = PetscViewerCreate(PETSC_COMM_WORLD,&jac->gkbviewer);CHKERRQ(ierr); 1072 ierr = PetscViewerSetType(jac->gkbviewer,PETSCVIEWERASCII);CHKERRQ(ierr); 1073 ierr = PetscObjectGetTabLevel((PetscObject)ilink->ksp,&tablevel);CHKERRQ(ierr); 1074 ierr = PetscViewerASCIISetTab(jac->gkbviewer,tablevel);CHKERRQ(ierr); 1075 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)ilink->ksp,1);CHKERRQ(ierr); 1076 } 1077 } else { 1078 /* set up the individual splits' PCs */ 1079 i = 0; 1080 ilink = jac->head; 1081 while (ilink) { 1082 ierr = KSPSetOperators(ilink->ksp,jac->mat[i],jac->pmat[i]);CHKERRQ(ierr); 1083 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 1084 if (!jac->suboptionsset) {ierr = KSPSetFromOptions(ilink->ksp);CHKERRQ(ierr);} 1085 i++; 1086 ilink = ilink->next; 1087 } 1088 } 1089 1090 jac->suboptionsset = PETSC_TRUE; 1091 PetscFunctionReturn(0); 1092 } 1093 1094 #define FieldSplitSplitSolveAdd(ilink,xx,yy) \ 1095 (VecScatterBegin(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 1096 VecScatterEnd(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 1097 PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL) ||\ 1098 KSPSolve(ilink->ksp,ilink->x,ilink->y) || \ 1099 KSPCheckSolve(ilink->ksp,pc,ilink->y) || \ 1100 PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL) ||\ 1101 VecScatterBegin(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE) || \ 1102 VecScatterEnd(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE)) 1103 1104 static PetscErrorCode PCApply_FieldSplit_Schur(PC pc,Vec x,Vec y) 1105 { 1106 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1107 PetscErrorCode ierr; 1108 PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next; 1109 KSP kspA = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper; 1110 1111 PetscFunctionBegin; 1112 switch (jac->schurfactorization) { 1113 case PC_FIELDSPLIT_SCHUR_FACT_DIAG: 1114 /* [A00 0; 0 -S], positive definite, suitable for MINRES */ 1115 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1116 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1117 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1118 ierr = PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 1119 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 1120 ierr = KSPCheckSolve(kspA,pc,ilinkA->y);CHKERRQ(ierr); 1121 ierr = PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 1122 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1123 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1124 ierr = PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);CHKERRQ(ierr); 1125 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 1126 ierr = KSPCheckSolve(jac->kspschur,pc,ilinkD->y);CHKERRQ(ierr); 1127 ierr = PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);CHKERRQ(ierr); 1128 ierr = VecScale(ilinkD->y,jac->schurscale);CHKERRQ(ierr); 1129 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1130 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1131 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1132 break; 1133 case PC_FIELDSPLIT_SCHUR_FACT_LOWER: 1134 /* [A00 0; A10 S], suitable for left preconditioning */ 1135 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1136 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1137 ierr = PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 1138 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 1139 ierr = KSPCheckSolve(kspA,pc,ilinkA->y);CHKERRQ(ierr); 1140 ierr = PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 1141 ierr = MatMult(jac->C,ilinkA->y,ilinkD->x);CHKERRQ(ierr); 1142 ierr = VecScale(ilinkD->x,-1.);CHKERRQ(ierr); 1143 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1144 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1145 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1146 ierr = PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);CHKERRQ(ierr); 1147 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 1148 ierr = KSPCheckSolve(jac->kspschur,pc,ilinkD->y);CHKERRQ(ierr); 1149 ierr = PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);CHKERRQ(ierr); 1150 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1151 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1152 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1153 break; 1154 case PC_FIELDSPLIT_SCHUR_FACT_UPPER: 1155 /* [A00 A01; 0 S], suitable for right preconditioning */ 1156 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1157 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1158 ierr = PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);CHKERRQ(ierr); 1159 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 1160 ierr = KSPCheckSolve(jac->kspschur,pc,ilinkD->y);CHKERRQ(ierr); 1161 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); 1162 ierr = VecScale(ilinkA->x,-1.);CHKERRQ(ierr); 1163 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1164 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1165 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1166 ierr = PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 1167 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 1168 ierr = KSPCheckSolve(kspA,pc,ilinkA->y);CHKERRQ(ierr); 1169 ierr = PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 1170 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1171 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1172 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1173 break; 1174 case PC_FIELDSPLIT_SCHUR_FACT_FULL: 1175 /* [1 0; A10 A00^{-1} 1] [A00 0; 0 S] [1 A00^{-1}A01; 0 1] */ 1176 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1177 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1178 ierr = PetscLogEventBegin(KSP_Solve_FS_L,kspLower,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 1179 ierr = KSPSolve(kspLower,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 1180 ierr = KSPCheckSolve(kspLower,pc,ilinkA->y);CHKERRQ(ierr); 1181 ierr = PetscLogEventEnd(KSP_Solve_FS_L,kspLower,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 1182 ierr = MatMult(jac->C,ilinkA->y,ilinkD->x);CHKERRQ(ierr); 1183 ierr = VecScale(ilinkD->x,-1.0);CHKERRQ(ierr); 1184 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1185 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1186 1187 ierr = PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);CHKERRQ(ierr); 1188 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 1189 ierr = KSPCheckSolve(jac->kspschur,pc,ilinkD->y);CHKERRQ(ierr); 1190 ierr = PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);CHKERRQ(ierr); 1191 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1192 1193 if (kspUpper == kspA) { 1194 ierr = MatMult(jac->B,ilinkD->y,ilinkA->y);CHKERRQ(ierr); 1195 ierr = VecAXPY(ilinkA->x,-1.0,ilinkA->y);CHKERRQ(ierr); 1196 ierr = PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 1197 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 1198 ierr = KSPCheckSolve(kspA,pc,ilinkA->y);CHKERRQ(ierr); 1199 ierr = PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 1200 } else { 1201 ierr = PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 1202 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 1203 ierr = KSPCheckSolve(kspA,pc,ilinkA->y);CHKERRQ(ierr); 1204 ierr = MatMult(jac->B,ilinkD->y,ilinkA->x);CHKERRQ(ierr); 1205 ierr = PetscLogEventBegin(KSP_Solve_FS_U,kspUpper,ilinkA->x,ilinkA->z,NULL);CHKERRQ(ierr); 1206 ierr = KSPSolve(kspUpper,ilinkA->x,ilinkA->z);CHKERRQ(ierr); 1207 ierr = KSPCheckSolve(kspUpper,pc,ilinkA->z);CHKERRQ(ierr); 1208 ierr = PetscLogEventEnd(KSP_Solve_FS_U,kspUpper,ilinkA->x,ilinkA->z,NULL);CHKERRQ(ierr); 1209 ierr = VecAXPY(ilinkA->y,-1.0,ilinkA->z);CHKERRQ(ierr); 1210 } 1211 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1212 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1213 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1214 } 1215 PetscFunctionReturn(0); 1216 } 1217 1218 static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y) 1219 { 1220 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1221 PetscErrorCode ierr; 1222 PC_FieldSplitLink ilink = jac->head; 1223 PetscInt cnt,bs; 1224 1225 PetscFunctionBegin; 1226 if (jac->type == PC_COMPOSITE_ADDITIVE) { 1227 if (jac->defaultsplit) { 1228 ierr = VecGetBlockSize(x,&bs);CHKERRQ(ierr); 1229 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); 1230 ierr = VecGetBlockSize(y,&bs);CHKERRQ(ierr); 1231 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); 1232 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 1233 while (ilink) { 1234 ierr = PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1235 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 1236 ierr = KSPCheckSolve(ilink->ksp,pc,ilink->y);CHKERRQ(ierr); 1237 ierr = PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1238 ilink = ilink->next; 1239 } 1240 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 1241 } else { 1242 ierr = VecSet(y,0.0);CHKERRQ(ierr); 1243 while (ilink) { 1244 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 1245 ilink = ilink->next; 1246 } 1247 } 1248 } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE && jac->nsplits == 2) { 1249 ierr = VecSet(y,0.0);CHKERRQ(ierr); 1250 /* solve on first block for first block variables */ 1251 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1252 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1253 ierr = PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1254 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 1255 ierr = KSPCheckSolve(ilink->ksp,pc,ilink->y);CHKERRQ(ierr); 1256 ierr = PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1257 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1258 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1259 1260 /* compute the residual only onto second block variables using first block variables */ 1261 ierr = MatMult(jac->Afield[1],ilink->y,ilink->next->x);CHKERRQ(ierr); 1262 ilink = ilink->next; 1263 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 1264 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1265 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1266 1267 /* solve on second block variables */ 1268 ierr = PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1269 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 1270 ierr = KSPCheckSolve(ilink->ksp,pc,ilink->y);CHKERRQ(ierr); 1271 ierr = PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1272 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1273 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1274 } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 1275 if (!jac->w1) { 1276 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 1277 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 1278 } 1279 ierr = VecSet(y,0.0);CHKERRQ(ierr); 1280 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 1281 cnt = 1; 1282 while (ilink->next) { 1283 ilink = ilink->next; 1284 /* compute the residual only over the part of the vector needed */ 1285 ierr = MatMult(jac->Afield[cnt++],y,ilink->x);CHKERRQ(ierr); 1286 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 1287 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1288 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1289 ierr = PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1290 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 1291 ierr = KSPCheckSolve(ilink->ksp,pc,ilink->y);CHKERRQ(ierr); 1292 ierr = PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1293 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1294 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1295 } 1296 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 1297 cnt -= 2; 1298 while (ilink->previous) { 1299 ilink = ilink->previous; 1300 /* compute the residual only over the part of the vector needed */ 1301 ierr = MatMult(jac->Afield[cnt--],y,ilink->x);CHKERRQ(ierr); 1302 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 1303 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1304 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1305 ierr = PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1306 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 1307 ierr = KSPCheckSolve(ilink->ksp,pc,ilink->y);CHKERRQ(ierr); 1308 ierr = PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1309 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1310 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1311 } 1312 } 1313 } else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Unsupported or unknown composition",(int) jac->type); 1314 PetscFunctionReturn(0); 1315 } 1316 1317 1318 static PetscErrorCode PCApply_FieldSplit_GKB(PC pc,Vec x,Vec y) 1319 { 1320 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1321 PetscErrorCode ierr; 1322 PC_FieldSplitLink ilinkA = jac->head,ilinkD = ilinkA->next; 1323 KSP ksp = ilinkA->ksp; 1324 Vec u,v,Hu,d,work1,work2; 1325 PetscScalar alpha,z,nrmz2,*vecz; 1326 PetscReal lowbnd,nu,beta; 1327 PetscInt j,iterGKB; 1328 1329 PetscFunctionBegin; 1330 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1331 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1332 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1333 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1334 1335 u = jac->u; 1336 v = jac->v; 1337 Hu = jac->Hu; 1338 d = jac->d; 1339 work1 = jac->w1; 1340 work2 = jac->w2; 1341 vecz = jac->vecz; 1342 1343 /* Change RHS to comply with matrix regularization H = A + nu*B*B' */ 1344 /* Add q = q + nu*B*b */ 1345 if (jac->gkbnu) { 1346 nu = jac->gkbnu; 1347 ierr = VecScale(ilinkD->x,jac->gkbnu);CHKERRQ(ierr); 1348 ierr = MatMultAdd(jac->B,ilinkD->x,ilinkA->x,ilinkA->x);CHKERRQ(ierr); /* q = q + nu*B*b */ 1349 } else { 1350 /* Situation when no augmented Lagrangian is used. Then we set inner */ 1351 /* matrix N = I in [Ar13], and thus nu = 1. */ 1352 nu = 1; 1353 } 1354 1355 /* Transform rhs from [q,tilde{b}] to [0,b] */ 1356 ierr = PetscLogEventBegin(ilinkA->event,ksp,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 1357 ierr = KSPSolve(ksp,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 1358 ierr = KSPCheckSolve(ksp,pc,ilinkA->y);CHKERRQ(ierr); 1359 ierr = PetscLogEventEnd(ilinkA->event,ksp,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 1360 ierr = MatMultHermitianTranspose(jac->B,ilinkA->y,work1);CHKERRQ(ierr); 1361 ierr = VecAXPBY(work1,1.0/nu,-1.0,ilinkD->x);CHKERRQ(ierr); /* c = b - B'*x */ 1362 1363 /* First step of algorithm */ 1364 ierr = VecNorm(work1,NORM_2,&beta);CHKERRQ(ierr); /* beta = sqrt(nu*c'*c)*/ 1365 KSPCheckDot(ksp,beta); 1366 beta = PetscSqrtScalar(nu)*beta; 1367 ierr = VecAXPBY(v,nu/beta,0.0,work1);CHKERRQ(ierr); /* v = nu/beta *c */ 1368 ierr = MatMult(jac->B,v,work2);CHKERRQ(ierr); /* u = H^{-1}*B*v */ 1369 ierr = PetscLogEventBegin(ilinkA->event,ksp,work2,u,NULL);CHKERRQ(ierr); 1370 ierr = KSPSolve(ksp,work2,u);CHKERRQ(ierr); 1371 ierr = KSPCheckSolve(ksp,pc,u);CHKERRQ(ierr); 1372 ierr = PetscLogEventEnd(ilinkA->event,ksp,work2,u,NULL);CHKERRQ(ierr); 1373 ierr = MatMult(jac->H,u,Hu);CHKERRQ(ierr); /* alpha = u'*H*u */ 1374 ierr = VecDot(Hu,u,&alpha);CHKERRQ(ierr); 1375 KSPCheckDot(ksp,alpha); 1376 if (PetscRealPart(alpha) <= 0.0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_NOT_CONVERGED,"GKB preconditioner diverged, H is not positive definite"); 1377 alpha = PetscSqrtScalar(PetscAbsScalar(alpha)); 1378 ierr = VecScale(u,1.0/alpha);CHKERRQ(ierr); 1379 ierr = VecAXPBY(d,1.0/alpha,0.0,v);CHKERRQ(ierr); /* v = nu/beta *c */ 1380 1381 z = beta/alpha; 1382 vecz[1] = z; 1383 1384 /* Computation of first iterate x(1) and p(1) */ 1385 ierr = VecAXPY(ilinkA->y,z,u);CHKERRQ(ierr); 1386 ierr = VecCopy(d,ilinkD->y);CHKERRQ(ierr); 1387 ierr = VecScale(ilinkD->y,-z);CHKERRQ(ierr); 1388 1389 iterGKB = 1; lowbnd = 2*jac->gkbtol; 1390 if (jac->gkbmonitor) { 1391 ierr = PetscViewerASCIIPrintf(jac->gkbviewer,"%3D GKB Lower bound estimate %14.12e\n",iterGKB,lowbnd);CHKERRQ(ierr); 1392 } 1393 1394 while (iterGKB < jac->gkbmaxit && lowbnd > jac->gkbtol) { 1395 iterGKB += 1; 1396 ierr = MatMultHermitianTranspose(jac->B,u,work1);CHKERRQ(ierr); /* v <- nu*(B'*u-alpha/nu*v) */ 1397 ierr = VecAXPBY(v,nu,-alpha,work1);CHKERRQ(ierr); 1398 ierr = VecNorm(v,NORM_2,&beta);CHKERRQ(ierr); /* beta = sqrt(nu)*v'*v */ 1399 beta = beta/PetscSqrtScalar(nu); 1400 ierr = VecScale(v,1.0/beta);CHKERRQ(ierr); 1401 ierr = MatMult(jac->B,v,work2);CHKERRQ(ierr); /* u <- H^{-1}*(B*v-beta*H*u) */ 1402 ierr = MatMult(jac->H,u,Hu);CHKERRQ(ierr); 1403 ierr = VecAXPY(work2,-beta,Hu);CHKERRQ(ierr); 1404 ierr = PetscLogEventBegin(ilinkA->event,ksp,work2,u,NULL);CHKERRQ(ierr); 1405 ierr = KSPSolve(ksp,work2,u);CHKERRQ(ierr); 1406 ierr = KSPCheckSolve(ksp,pc,u);CHKERRQ(ierr); 1407 ierr = PetscLogEventEnd(ilinkA->event,ksp,work2,u,NULL);CHKERRQ(ierr); 1408 ierr = MatMult(jac->H,u,Hu);CHKERRQ(ierr); /* alpha = u'*H*u */ 1409 ierr = VecDot(Hu,u,&alpha);CHKERRQ(ierr); 1410 KSPCheckDot(ksp,alpha); 1411 if (PetscRealPart(alpha) <= 0.0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_NOT_CONVERGED,"GKB preconditioner diverged, H is not positive definite"); 1412 alpha = PetscSqrtScalar(PetscAbsScalar(alpha)); 1413 ierr = VecScale(u,1.0/alpha);CHKERRQ(ierr); 1414 1415 z = -beta/alpha*z; /* z <- beta/alpha*z */ 1416 vecz[0] = z; 1417 1418 /* Computation of new iterate x(i+1) and p(i+1) */ 1419 ierr = VecAXPBY(d,1.0/alpha,-beta/alpha,v);CHKERRQ(ierr); /* d = (v-beta*d)/alpha */ 1420 ierr = VecAXPY(ilinkA->y,z,u);CHKERRQ(ierr); /* r = r + z*u */ 1421 ierr = VecAXPY(ilinkD->y,-z,d);CHKERRQ(ierr); /* p = p - z*d */ 1422 ierr = MatMult(jac->H,ilinkA->y,Hu);CHKERRQ(ierr); /* ||u||_H = u'*H*u */ 1423 ierr = VecDot(Hu,ilinkA->y,&nrmz2);CHKERRQ(ierr); 1424 1425 /* Compute Lower Bound estimate */ 1426 if (iterGKB > jac->gkbdelay) { 1427 lowbnd = 0.0; 1428 for (j=0; j<jac->gkbdelay; j++) { 1429 lowbnd += PetscAbsScalar(vecz[j]*vecz[j]); 1430 } 1431 lowbnd = PetscSqrtScalar(lowbnd/PetscAbsScalar(nrmz2)); 1432 } 1433 1434 for (j=0; j<jac->gkbdelay-1; j++) { 1435 vecz[jac->gkbdelay-j-1] = vecz[jac->gkbdelay-j-2]; 1436 } 1437 if (jac->gkbmonitor) { 1438 ierr = PetscViewerASCIIPrintf(jac->gkbviewer,"%3D GKB Lower bound estimate %14.12e\n",iterGKB,lowbnd);CHKERRQ(ierr); 1439 } 1440 } 1441 1442 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1443 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1444 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1445 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1446 1447 PetscFunctionReturn(0); 1448 } 1449 1450 1451 #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \ 1452 (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 1453 VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 1454 PetscLogEventBegin(ilink->event,ilink->ksp,ilink->y,ilink->x,NULL) || \ 1455 KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) || \ 1456 KSPCheckSolve(ilink->ksp,pc,ilink->x) || \ 1457 PetscLogEventEnd(ilink->event,ilink->ksp,ilink->y,ilink->x,NULL) || \ 1458 VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \ 1459 VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE)) 1460 1461 static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y) 1462 { 1463 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1464 PetscErrorCode ierr; 1465 PC_FieldSplitLink ilink = jac->head; 1466 PetscInt bs; 1467 1468 PetscFunctionBegin; 1469 if (jac->type == PC_COMPOSITE_ADDITIVE) { 1470 if (jac->defaultsplit) { 1471 ierr = VecGetBlockSize(x,&bs);CHKERRQ(ierr); 1472 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); 1473 ierr = VecGetBlockSize(y,&bs);CHKERRQ(ierr); 1474 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); 1475 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 1476 while (ilink) { 1477 ierr = PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1478 ierr = KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 1479 ierr = KSPCheckSolve(ilink->ksp,pc,ilink->y);CHKERRQ(ierr); 1480 ierr = PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1481 ilink = ilink->next; 1482 } 1483 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 1484 } else { 1485 ierr = VecSet(y,0.0);CHKERRQ(ierr); 1486 while (ilink) { 1487 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 1488 ilink = ilink->next; 1489 } 1490 } 1491 } else { 1492 if (!jac->w1) { 1493 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 1494 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 1495 } 1496 ierr = VecSet(y,0.0);CHKERRQ(ierr); 1497 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 1498 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 1499 while (ilink->next) { 1500 ilink = ilink->next; 1501 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 1502 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 1503 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 1504 } 1505 while (ilink->previous) { 1506 ilink = ilink->previous; 1507 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 1508 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 1509 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 1510 } 1511 } else { 1512 while (ilink->next) { /* get to last entry in linked list */ 1513 ilink = ilink->next; 1514 } 1515 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 1516 while (ilink->previous) { 1517 ilink = ilink->previous; 1518 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 1519 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 1520 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 1521 } 1522 } 1523 } 1524 PetscFunctionReturn(0); 1525 } 1526 1527 static PetscErrorCode PCReset_FieldSplit(PC pc) 1528 { 1529 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1530 PetscErrorCode ierr; 1531 PC_FieldSplitLink ilink = jac->head,next; 1532 1533 PetscFunctionBegin; 1534 while (ilink) { 1535 ierr = KSPDestroy(&ilink->ksp);CHKERRQ(ierr); 1536 ierr = VecDestroy(&ilink->x);CHKERRQ(ierr); 1537 ierr = VecDestroy(&ilink->y);CHKERRQ(ierr); 1538 ierr = VecDestroy(&ilink->z);CHKERRQ(ierr); 1539 ierr = VecScatterDestroy(&ilink->sctx);CHKERRQ(ierr); 1540 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 1541 ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr); 1542 ierr = PetscFree(ilink->splitname);CHKERRQ(ierr); 1543 ierr = PetscFree(ilink->fields);CHKERRQ(ierr); 1544 ierr = PetscFree(ilink->fields_col);CHKERRQ(ierr); 1545 next = ilink->next; 1546 ierr = PetscFree(ilink);CHKERRQ(ierr); 1547 ilink = next; 1548 } 1549 jac->head = NULL; 1550 ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr); 1551 if (jac->mat && jac->mat != jac->pmat) { 1552 ierr = MatDestroyMatrices(jac->nsplits,&jac->mat);CHKERRQ(ierr); 1553 } else if (jac->mat) { 1554 jac->mat = NULL; 1555 } 1556 if (jac->pmat) {ierr = MatDestroyMatrices(jac->nsplits,&jac->pmat);CHKERRQ(ierr);} 1557 if (jac->Afield) {ierr = MatDestroyMatrices(jac->nsplits,&jac->Afield);CHKERRQ(ierr);} 1558 jac->nsplits = 0; 1559 ierr = VecDestroy(&jac->w1);CHKERRQ(ierr); 1560 ierr = VecDestroy(&jac->w2);CHKERRQ(ierr); 1561 ierr = MatDestroy(&jac->schur);CHKERRQ(ierr); 1562 ierr = MatDestroy(&jac->schurp);CHKERRQ(ierr); 1563 ierr = MatDestroy(&jac->schur_user);CHKERRQ(ierr); 1564 ierr = KSPDestroy(&jac->kspschur);CHKERRQ(ierr); 1565 ierr = KSPDestroy(&jac->kspupper);CHKERRQ(ierr); 1566 ierr = MatDestroy(&jac->B);CHKERRQ(ierr); 1567 ierr = MatDestroy(&jac->C);CHKERRQ(ierr); 1568 ierr = MatDestroy(&jac->H);CHKERRQ(ierr); 1569 ierr = VecDestroy(&jac->u);CHKERRQ(ierr); 1570 ierr = VecDestroy(&jac->v);CHKERRQ(ierr); 1571 ierr = VecDestroy(&jac->Hu);CHKERRQ(ierr); 1572 ierr = VecDestroy(&jac->d);CHKERRQ(ierr); 1573 ierr = PetscFree(jac->vecz);CHKERRQ(ierr); 1574 ierr = PetscViewerDestroy(&jac->gkbviewer);CHKERRQ(ierr); 1575 jac->isrestrict = PETSC_FALSE; 1576 PetscFunctionReturn(0); 1577 } 1578 1579 static PetscErrorCode PCDestroy_FieldSplit(PC pc) 1580 { 1581 PetscErrorCode ierr; 1582 1583 PetscFunctionBegin; 1584 ierr = PCReset_FieldSplit(pc);CHKERRQ(ierr); 1585 ierr = PetscFree(pc->data);CHKERRQ(ierr); 1586 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSchurGetSubKSP_C",NULL);CHKERRQ(ierr); 1587 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",NULL);CHKERRQ(ierr); 1588 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",NULL);CHKERRQ(ierr); 1589 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",NULL);CHKERRQ(ierr); 1590 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",NULL);CHKERRQ(ierr); 1591 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",NULL);CHKERRQ(ierr); 1592 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",NULL);CHKERRQ(ierr); 1593 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",NULL);CHKERRQ(ierr); 1594 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",NULL);CHKERRQ(ierr); 1595 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitRestrictIS_C",NULL);CHKERRQ(ierr); 1596 PetscFunctionReturn(0); 1597 } 1598 1599 static PetscErrorCode PCSetFromOptions_FieldSplit(PetscOptionItems *PetscOptionsObject,PC pc) 1600 { 1601 PetscErrorCode ierr; 1602 PetscInt bs; 1603 PetscBool flg; 1604 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1605 PCCompositeType ctype; 1606 1607 PetscFunctionBegin; 1608 ierr = PetscOptionsHead(PetscOptionsObject,"FieldSplit options");CHKERRQ(ierr); 1609 ierr = PetscOptionsBool("-pc_fieldsplit_dm_splits","Whether to use DMCreateFieldDecomposition() for splits","PCFieldSplitSetDMSplits",jac->dm_splits,&jac->dm_splits,NULL);CHKERRQ(ierr); 1610 ierr = PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);CHKERRQ(ierr); 1611 if (flg) { 1612 ierr = PCFieldSplitSetBlockSize(pc,bs);CHKERRQ(ierr); 1613 } 1614 jac->diag_use_amat = pc->useAmat; 1615 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); 1616 jac->offdiag_use_amat = pc->useAmat; 1617 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); 1618 ierr = PetscOptionsBool("-pc_fieldsplit_detect_saddle_point","Form 2-way split by detecting zero diagonal entries", "PCFieldSplitSetDetectSaddlePoint",jac->detect,&jac->detect,NULL);CHKERRQ(ierr); 1619 ierr = PCFieldSplitSetDetectSaddlePoint(pc,jac->detect);CHKERRQ(ierr); /* Sets split type and Schur PC type */ 1620 ierr = PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg);CHKERRQ(ierr); 1621 if (flg) { 1622 ierr = PCFieldSplitSetType(pc,ctype);CHKERRQ(ierr); 1623 } 1624 /* Only setup fields once */ 1625 if ((jac->bs > 0) && (jac->nsplits == 0)) { 1626 /* only allow user to set fields from command line if bs is already known. 1627 otherwise user can set them in PCFieldSplitSetDefaults() */ 1628 ierr = PCFieldSplitSetRuntimeSplits_Private(pc);CHKERRQ(ierr); 1629 if (jac->splitdefined) {ierr = PetscInfo(pc,"Splits defined using the options database\n");CHKERRQ(ierr);} 1630 } 1631 if (jac->type == PC_COMPOSITE_SCHUR) { 1632 ierr = PetscOptionsGetEnum(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_schur_factorization_type",PCFieldSplitSchurFactTypes,(PetscEnum*)&jac->schurfactorization,&flg);CHKERRQ(ierr); 1633 if (flg) {ierr = PetscInfo(pc,"Deprecated use of -pc_fieldsplit_schur_factorization_type\n");CHKERRQ(ierr);} 1634 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); 1635 ierr = PetscOptionsEnum("-pc_fieldsplit_schur_precondition","How to build preconditioner for Schur complement","PCFieldSplitSetSchurPre",PCFieldSplitSchurPreTypes,(PetscEnum)jac->schurpre,(PetscEnum*)&jac->schurpre,NULL);CHKERRQ(ierr); 1636 ierr = PetscOptionsScalar("-pc_fieldsplit_schur_scale","Scale Schur complement","PCFieldSplitSetSchurScale",jac->schurscale,&jac->schurscale,NULL);CHKERRQ(ierr); 1637 } else if (jac->type == PC_COMPOSITE_GKB) { 1638 ierr = PetscOptionsReal("-pc_fieldsplit_gkb_tol","The tolerance for the lower bound stopping criterion","PCFieldSplitGKBTol",jac->gkbtol,&jac->gkbtol,NULL);CHKERRQ(ierr); 1639 ierr = PetscOptionsInt("-pc_fieldsplit_gkb_delay","The delay value for lower bound criterion","PCFieldSplitGKBDelay",jac->gkbdelay,&jac->gkbdelay,NULL);CHKERRQ(ierr); 1640 ierr = PetscOptionsReal("-pc_fieldsplit_gkb_nu","Parameter in augmented Lagrangian approach","PCFieldSplitGKBNu",jac->gkbnu,&jac->gkbnu,NULL);CHKERRQ(ierr); 1641 if (jac->gkbnu < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nu cannot be less than 0: value %f",jac->gkbnu); 1642 ierr = PetscOptionsInt("-pc_fieldsplit_gkb_maxit","Maximum allowed number of iterations","PCFieldSplitGKBMaxit",jac->gkbmaxit,&jac->gkbmaxit,NULL);CHKERRQ(ierr); 1643 ierr = PetscOptionsBool("-pc_fieldsplit_gkb_monitor","Prints number of GKB iterations and error","PCFieldSplitGKB",jac->gkbmonitor,&jac->gkbmonitor,NULL);CHKERRQ(ierr); 1644 } 1645 ierr = PetscOptionsTail();CHKERRQ(ierr); 1646 PetscFunctionReturn(0); 1647 } 1648 1649 /*------------------------------------------------------------------------------------*/ 1650 1651 static PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 1652 { 1653 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1654 PetscErrorCode ierr; 1655 PC_FieldSplitLink ilink,next = jac->head; 1656 char prefix[128]; 1657 PetscInt i; 1658 1659 PetscFunctionBegin; 1660 if (jac->splitdefined) { 1661 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 1662 PetscFunctionReturn(0); 1663 } 1664 for (i=0; i<n; i++) { 1665 if (fields[i] >= jac->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs); 1666 if (fields[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]); 1667 } 1668 ierr = PetscNew(&ilink);CHKERRQ(ierr); 1669 if (splitname) { 1670 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 1671 } else { 1672 ierr = PetscMalloc1(3,&ilink->splitname);CHKERRQ(ierr); 1673 ierr = PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);CHKERRQ(ierr); 1674 } 1675 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 */ 1676 ierr = PetscMalloc1(n,&ilink->fields);CHKERRQ(ierr); 1677 ierr = PetscArraycpy(ilink->fields,fields,n);CHKERRQ(ierr); 1678 ierr = PetscMalloc1(n,&ilink->fields_col);CHKERRQ(ierr); 1679 ierr = PetscArraycpy(ilink->fields_col,fields_col,n);CHKERRQ(ierr); 1680 1681 ilink->nfields = n; 1682 ilink->next = NULL; 1683 ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);CHKERRQ(ierr); 1684 ierr = KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure);CHKERRQ(ierr); 1685 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 1686 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 1687 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 1688 1689 ierr = PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);CHKERRQ(ierr); 1690 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 1691 1692 if (!next) { 1693 jac->head = ilink; 1694 ilink->previous = NULL; 1695 } else { 1696 while (next->next) { 1697 next = next->next; 1698 } 1699 next->next = ilink; 1700 ilink->previous = next; 1701 } 1702 jac->nsplits++; 1703 PetscFunctionReturn(0); 1704 } 1705 1706 static PetscErrorCode PCFieldSplitSchurGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp) 1707 { 1708 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1709 PetscErrorCode ierr; 1710 1711 PetscFunctionBegin; 1712 *subksp = NULL; 1713 if (n) *n = 0; 1714 if (jac->type == PC_COMPOSITE_SCHUR) { 1715 PetscInt nn; 1716 1717 if (!jac->schur) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call KSPSetUp() or PCSetUp() before calling PCFieldSplitSchurGetSubKSP()"); 1718 if (jac->nsplits != 2) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Unexpected number of splits %D != 2",jac->nsplits); 1719 nn = jac->nsplits + (jac->kspupper != jac->head->ksp ? 1 : 0); 1720 ierr = PetscMalloc1(nn,subksp);CHKERRQ(ierr); 1721 (*subksp)[0] = jac->head->ksp; 1722 (*subksp)[1] = jac->kspschur; 1723 if (jac->kspupper != jac->head->ksp) (*subksp)[2] = jac->kspupper; 1724 if (n) *n = nn; 1725 } 1726 PetscFunctionReturn(0); 1727 } 1728 1729 static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp) 1730 { 1731 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1732 PetscErrorCode ierr; 1733 1734 PetscFunctionBegin; 1735 if (!jac->schur) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call KSPSetUp() or PCSetUp() before calling PCFieldSplitGetSubKSP()"); 1736 ierr = PetscMalloc1(jac->nsplits,subksp);CHKERRQ(ierr); 1737 ierr = MatSchurComplementGetKSP(jac->schur,*subksp);CHKERRQ(ierr); 1738 1739 (*subksp)[1] = jac->kspschur; 1740 if (n) *n = jac->nsplits; 1741 PetscFunctionReturn(0); 1742 } 1743 1744 static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp) 1745 { 1746 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1747 PetscErrorCode ierr; 1748 PetscInt cnt = 0; 1749 PC_FieldSplitLink ilink = jac->head; 1750 1751 PetscFunctionBegin; 1752 ierr = PetscMalloc1(jac->nsplits,subksp);CHKERRQ(ierr); 1753 while (ilink) { 1754 (*subksp)[cnt++] = ilink->ksp; 1755 ilink = ilink->next; 1756 } 1757 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); 1758 if (n) *n = jac->nsplits; 1759 PetscFunctionReturn(0); 1760 } 1761 1762 /*@C 1763 PCFieldSplitRestrictIS - Restricts the fieldsplit ISs to be within a given IS. 1764 1765 Input Parameters: 1766 + pc - the preconditioner context 1767 - is - the index set that defines the indices to which the fieldsplit is to be restricted 1768 1769 Level: advanced 1770 1771 @*/ 1772 PetscErrorCode PCFieldSplitRestrictIS(PC pc,IS isy) 1773 { 1774 PetscErrorCode ierr; 1775 1776 PetscFunctionBegin; 1777 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1778 PetscValidHeaderSpecific(isy,IS_CLASSID,2); 1779 ierr = PetscTryMethod(pc,"PCFieldSplitRestrictIS_C",(PC,IS),(pc,isy));CHKERRQ(ierr); 1780 PetscFunctionReturn(0); 1781 } 1782 1783 1784 static PetscErrorCode PCFieldSplitRestrictIS_FieldSplit(PC pc, IS isy) 1785 { 1786 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1787 PetscErrorCode ierr; 1788 PC_FieldSplitLink ilink = jac->head, next; 1789 PetscInt localsize,size,sizez,i; 1790 const PetscInt *ind, *indz; 1791 PetscInt *indc, *indcz; 1792 PetscBool flg; 1793 1794 PetscFunctionBegin; 1795 ierr = ISGetLocalSize(isy,&localsize);CHKERRQ(ierr); 1796 ierr = MPI_Scan(&localsize,&size,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)isy));CHKERRQ(ierr); 1797 size -= localsize; 1798 while(ilink) { 1799 IS isrl,isr; 1800 PC subpc; 1801 ierr = ISEmbed(ilink->is, isy, PETSC_TRUE, &isrl);CHKERRQ(ierr); 1802 ierr = ISGetLocalSize(isrl,&localsize);CHKERRQ(ierr); 1803 ierr = PetscMalloc1(localsize,&indc);CHKERRQ(ierr); 1804 ierr = ISGetIndices(isrl,&ind);CHKERRQ(ierr); 1805 ierr = PetscArraycpy(indc,ind,localsize);CHKERRQ(ierr); 1806 ierr = ISRestoreIndices(isrl,&ind);CHKERRQ(ierr); 1807 ierr = ISDestroy(&isrl);CHKERRQ(ierr); 1808 for (i=0; i<localsize; i++) *(indc+i) += size; 1809 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)isy),localsize,indc,PETSC_OWN_POINTER,&isr);CHKERRQ(ierr); 1810 ierr = PetscObjectReference((PetscObject)isr);CHKERRQ(ierr); 1811 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 1812 ilink->is = isr; 1813 ierr = PetscObjectReference((PetscObject)isr);CHKERRQ(ierr); 1814 ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr); 1815 ilink->is_col = isr; 1816 ierr = ISDestroy(&isr);CHKERRQ(ierr); 1817 ierr = KSPGetPC(ilink->ksp, &subpc);CHKERRQ(ierr); 1818 ierr = PetscObjectTypeCompare((PetscObject)subpc,PCFIELDSPLIT,&flg);CHKERRQ(ierr); 1819 if(flg) { 1820 IS iszl,isz; 1821 MPI_Comm comm; 1822 ierr = ISGetLocalSize(ilink->is,&localsize);CHKERRQ(ierr); 1823 comm = PetscObjectComm((PetscObject)ilink->is); 1824 ierr = ISEmbed(isy, ilink->is, PETSC_TRUE, &iszl);CHKERRQ(ierr); 1825 ierr = MPI_Scan(&localsize,&sizez,1,MPIU_INT,MPI_SUM,comm);CHKERRQ(ierr); 1826 sizez -= localsize; 1827 ierr = ISGetLocalSize(iszl,&localsize);CHKERRQ(ierr); 1828 ierr = PetscMalloc1(localsize,&indcz);CHKERRQ(ierr); 1829 ierr = ISGetIndices(iszl,&indz);CHKERRQ(ierr); 1830 ierr = PetscArraycpy(indcz,indz,localsize);CHKERRQ(ierr); 1831 ierr = ISRestoreIndices(iszl,&indz);CHKERRQ(ierr); 1832 ierr = ISDestroy(&iszl);CHKERRQ(ierr); 1833 for (i=0; i<localsize; i++) *(indcz+i) += sizez; 1834 ierr = ISCreateGeneral(comm,localsize,indcz,PETSC_OWN_POINTER,&isz);CHKERRQ(ierr); 1835 ierr = PCFieldSplitRestrictIS(subpc,isz);CHKERRQ(ierr); 1836 ierr = ISDestroy(&isz);CHKERRQ(ierr); 1837 } 1838 next = ilink->next; 1839 ilink = next; 1840 } 1841 jac->isrestrict = PETSC_TRUE; 1842 PetscFunctionReturn(0); 1843 } 1844 1845 static PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is) 1846 { 1847 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1848 PetscErrorCode ierr; 1849 PC_FieldSplitLink ilink, next = jac->head; 1850 char prefix[128]; 1851 1852 PetscFunctionBegin; 1853 if (jac->splitdefined) { 1854 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 1855 PetscFunctionReturn(0); 1856 } 1857 ierr = PetscNew(&ilink);CHKERRQ(ierr); 1858 if (splitname) { 1859 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 1860 } else { 1861 ierr = PetscMalloc1(8,&ilink->splitname);CHKERRQ(ierr); 1862 ierr = PetscSNPrintf(ilink->splitname,7,"%D",jac->nsplits);CHKERRQ(ierr); 1863 } 1864 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 */ 1865 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 1866 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 1867 ilink->is = is; 1868 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 1869 ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr); 1870 ilink->is_col = is; 1871 ilink->next = NULL; 1872 ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);CHKERRQ(ierr); 1873 ierr = KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure);CHKERRQ(ierr); 1874 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 1875 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 1876 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 1877 1878 ierr = PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);CHKERRQ(ierr); 1879 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 1880 1881 if (!next) { 1882 jac->head = ilink; 1883 ilink->previous = NULL; 1884 } else { 1885 while (next->next) { 1886 next = next->next; 1887 } 1888 next->next = ilink; 1889 ilink->previous = next; 1890 } 1891 jac->nsplits++; 1892 PetscFunctionReturn(0); 1893 } 1894 1895 /*@C 1896 PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner 1897 1898 Logically Collective on PC 1899 1900 Input Parameters: 1901 + pc - the preconditioner context 1902 . splitname - name of this split, if NULL the number of the split is used 1903 . n - the number of fields in this split 1904 - fields - the fields in this split 1905 1906 Level: intermediate 1907 1908 Notes: 1909 Use PCFieldSplitSetIS() to set a completely general set of indices as a field. 1910 1911 The PCFieldSplitSetFields() is for defining fields as strided blocks. For example, if the block 1912 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 1913 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x.... 1914 where the numbered entries indicate what is in the field. 1915 1916 This function is called once per split (it creates a new split each time). Solve options 1917 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1918 1919 Developer Note: This routine does not actually create the IS representing the split, that is delayed 1920 until PCSetUp_FieldSplit(), because information about the vector/matrix layouts may not be 1921 available when this routine is called. 1922 1923 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize(), PCFieldSplitSetIS() 1924 1925 @*/ 1926 PetscErrorCode PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 1927 { 1928 PetscErrorCode ierr; 1929 1930 PetscFunctionBegin; 1931 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1932 PetscValidCharPointer(splitname,2); 1933 if (n < 1) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Provided number of fields %D in split \"%s\" not positive",n,splitname); 1934 PetscValidIntPointer(fields,3); 1935 ierr = PetscTryMethod(pc,"PCFieldSplitSetFields_C",(PC,const char[],PetscInt,const PetscInt*,const PetscInt*),(pc,splitname,n,fields,fields_col));CHKERRQ(ierr); 1936 PetscFunctionReturn(0); 1937 } 1938 1939 /*@ 1940 PCFieldSplitSetDiagUseAmat - set flag indicating whether to extract diagonal blocks from Amat (rather than Pmat) 1941 1942 Logically Collective on PC 1943 1944 Input Parameters: 1945 + pc - the preconditioner object 1946 - flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from 1947 1948 Options Database: 1949 . -pc_fieldsplit_diag_use_amat 1950 1951 Level: intermediate 1952 1953 .seealso: PCFieldSplitGetDiagUseAmat(), PCFieldSplitSetOffDiagUseAmat(), PCFIELDSPLIT 1954 1955 @*/ 1956 PetscErrorCode PCFieldSplitSetDiagUseAmat(PC pc,PetscBool flg) 1957 { 1958 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1959 PetscBool isfs; 1960 PetscErrorCode ierr; 1961 1962 PetscFunctionBegin; 1963 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1964 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1965 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 1966 jac->diag_use_amat = flg; 1967 PetscFunctionReturn(0); 1968 } 1969 1970 /*@ 1971 PCFieldSplitGetDiagUseAmat - get the flag indicating whether to extract diagonal blocks from Amat (rather than Pmat) 1972 1973 Logically Collective on PC 1974 1975 Input Parameters: 1976 . pc - the preconditioner object 1977 1978 Output Parameters: 1979 . flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from 1980 1981 1982 Level: intermediate 1983 1984 .seealso: PCFieldSplitSetDiagUseAmat(), PCFieldSplitGetOffDiagUseAmat(), PCFIELDSPLIT 1985 1986 @*/ 1987 PetscErrorCode PCFieldSplitGetDiagUseAmat(PC pc,PetscBool *flg) 1988 { 1989 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1990 PetscBool isfs; 1991 PetscErrorCode ierr; 1992 1993 PetscFunctionBegin; 1994 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1995 PetscValidBoolPointer(flg,2); 1996 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1997 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 1998 *flg = jac->diag_use_amat; 1999 PetscFunctionReturn(0); 2000 } 2001 2002 /*@ 2003 PCFieldSplitSetOffDiagUseAmat - set flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat) 2004 2005 Logically Collective on PC 2006 2007 Input Parameters: 2008 + pc - the preconditioner object 2009 - flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from 2010 2011 Options Database: 2012 . -pc_fieldsplit_off_diag_use_amat 2013 2014 Level: intermediate 2015 2016 .seealso: PCFieldSplitGetOffDiagUseAmat(), PCFieldSplitSetDiagUseAmat(), PCFIELDSPLIT 2017 2018 @*/ 2019 PetscErrorCode PCFieldSplitSetOffDiagUseAmat(PC pc,PetscBool flg) 2020 { 2021 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2022 PetscBool isfs; 2023 PetscErrorCode ierr; 2024 2025 PetscFunctionBegin; 2026 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2027 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 2028 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 2029 jac->offdiag_use_amat = flg; 2030 PetscFunctionReturn(0); 2031 } 2032 2033 /*@ 2034 PCFieldSplitGetOffDiagUseAmat - get the flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat) 2035 2036 Logically Collective on PC 2037 2038 Input Parameters: 2039 . pc - the preconditioner object 2040 2041 Output Parameters: 2042 . flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from 2043 2044 2045 Level: intermediate 2046 2047 .seealso: PCFieldSplitSetOffDiagUseAmat(), PCFieldSplitGetDiagUseAmat(), PCFIELDSPLIT 2048 2049 @*/ 2050 PetscErrorCode PCFieldSplitGetOffDiagUseAmat(PC pc,PetscBool *flg) 2051 { 2052 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2053 PetscBool isfs; 2054 PetscErrorCode ierr; 2055 2056 PetscFunctionBegin; 2057 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2058 PetscValidBoolPointer(flg,2); 2059 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 2060 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 2061 *flg = jac->offdiag_use_amat; 2062 PetscFunctionReturn(0); 2063 } 2064 2065 2066 2067 /*@C 2068 PCFieldSplitSetIS - Sets the exact elements for field 2069 2070 Logically Collective on PC 2071 2072 Input Parameters: 2073 + pc - the preconditioner context 2074 . splitname - name of this split, if NULL the number of the split is used 2075 - is - the index set that defines the vector elements in this field 2076 2077 2078 Notes: 2079 Use PCFieldSplitSetFields(), for fields defined by strided types. 2080 2081 This function is called once per split (it creates a new split each time). Solve options 2082 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 2083 2084 Level: intermediate 2085 2086 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize() 2087 2088 @*/ 2089 PetscErrorCode PCFieldSplitSetIS(PC pc,const char splitname[],IS is) 2090 { 2091 PetscErrorCode ierr; 2092 2093 PetscFunctionBegin; 2094 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2095 if (splitname) PetscValidCharPointer(splitname,2); 2096 PetscValidHeaderSpecific(is,IS_CLASSID,3); 2097 ierr = PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));CHKERRQ(ierr); 2098 PetscFunctionReturn(0); 2099 } 2100 2101 /*@C 2102 PCFieldSplitGetIS - Retrieves the elements for a field as an IS 2103 2104 Logically Collective on PC 2105 2106 Input Parameters: 2107 + pc - the preconditioner context 2108 - splitname - name of this split 2109 2110 Output Parameter: 2111 - is - the index set that defines the vector elements in this field, or NULL if the field is not found 2112 2113 Level: intermediate 2114 2115 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetIS() 2116 2117 @*/ 2118 PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is) 2119 { 2120 PetscErrorCode ierr; 2121 2122 PetscFunctionBegin; 2123 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2124 PetscValidCharPointer(splitname,2); 2125 PetscValidPointer(is,3); 2126 { 2127 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 2128 PC_FieldSplitLink ilink = jac->head; 2129 PetscBool found; 2130 2131 *is = NULL; 2132 while (ilink) { 2133 ierr = PetscStrcmp(ilink->splitname, splitname, &found);CHKERRQ(ierr); 2134 if (found) { 2135 *is = ilink->is; 2136 break; 2137 } 2138 ilink = ilink->next; 2139 } 2140 } 2141 PetscFunctionReturn(0); 2142 } 2143 2144 /*@C 2145 PCFieldSplitGetISByIndex - Retrieves the elements for a given index field as an IS 2146 2147 Logically Collective on PC 2148 2149 Input Parameters: 2150 + pc - the preconditioner context 2151 - index - index of this split 2152 2153 Output Parameter: 2154 - is - the index set that defines the vector elements in this field 2155 2156 Level: intermediate 2157 2158 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitGetIS(), PCFieldSplitSetIS() 2159 2160 @*/ 2161 PetscErrorCode PCFieldSplitGetISByIndex(PC pc,PetscInt index,IS *is) 2162 { 2163 PetscErrorCode ierr; 2164 2165 PetscFunctionBegin; 2166 if (index < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",index); 2167 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2168 PetscValidPointer(is,3); 2169 { 2170 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 2171 PC_FieldSplitLink ilink = jac->head; 2172 PetscInt i = 0; 2173 if (index >= jac->nsplits) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",index,jac->nsplits); 2174 2175 while (i < index) { 2176 ilink = ilink->next; 2177 ++i; 2178 } 2179 ierr = PCFieldSplitGetIS(pc,ilink->splitname,is);CHKERRQ(ierr); 2180 } 2181 PetscFunctionReturn(0); 2182 } 2183 2184 /*@ 2185 PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the 2186 fieldsplit preconditioner. If not set the matrix block size is used. 2187 2188 Logically Collective on PC 2189 2190 Input Parameters: 2191 + pc - the preconditioner context 2192 - bs - the block size 2193 2194 Level: intermediate 2195 2196 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields() 2197 2198 @*/ 2199 PetscErrorCode PCFieldSplitSetBlockSize(PC pc,PetscInt bs) 2200 { 2201 PetscErrorCode ierr; 2202 2203 PetscFunctionBegin; 2204 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2205 PetscValidLogicalCollectiveInt(pc,bs,2); 2206 ierr = PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));CHKERRQ(ierr); 2207 PetscFunctionReturn(0); 2208 } 2209 2210 /*@C 2211 PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits 2212 2213 Collective on KSP 2214 2215 Input Parameter: 2216 . pc - the preconditioner context 2217 2218 Output Parameters: 2219 + n - the number of splits 2220 - subksp - the array of KSP contexts 2221 2222 Note: 2223 After PCFieldSplitGetSubKSP() the array of KSPs is to be freed by the user with PetscFree() 2224 (not the KSP just the array that contains them). 2225 2226 You must call PCSetUp() before calling PCFieldSplitGetSubKSP(). 2227 2228 If the fieldsplit is of type PC_COMPOSITE_SCHUR, it returns the KSP object used inside the 2229 Schur complement and the KSP object used to iterate over the Schur complement. 2230 To access all the KSP objects used in PC_COMPOSITE_SCHUR, use PCFieldSplitSchurGetSubKSP(). 2231 2232 If the fieldsplit is of type PC_COMPOSITE_GKB, it returns the KSP object used to solve the 2233 inner linear system defined by the matrix H in each loop. 2234 2235 Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs. 2236 You can call PCFieldSplitGetSubKSP(pc,n,PETSC_NULL_KSP,ierr) to determine how large the 2237 KSP array must be. 2238 2239 2240 Level: advanced 2241 2242 .seealso: PCFIELDSPLIT 2243 @*/ 2244 PetscErrorCode PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[]) 2245 { 2246 PetscErrorCode ierr; 2247 2248 PetscFunctionBegin; 2249 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2250 if (n) PetscValidIntPointer(n,2); 2251 ierr = PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));CHKERRQ(ierr); 2252 PetscFunctionReturn(0); 2253 } 2254 2255 /*@C 2256 PCFieldSplitSchurGetSubKSP - Gets the KSP contexts used inside the Schur complement based PCFIELDSPLIT 2257 2258 Collective on KSP 2259 2260 Input Parameter: 2261 . pc - the preconditioner context 2262 2263 Output Parameters: 2264 + n - the number of splits 2265 - subksp - the array of KSP contexts 2266 2267 Note: 2268 After PCFieldSplitSchurGetSubKSP() the array of KSPs is to be freed by the user with PetscFree() 2269 (not the KSP just the array that contains them). 2270 2271 You must call PCSetUp() before calling PCFieldSplitSchurGetSubKSP(). 2272 2273 If the fieldsplit type is of type PC_COMPOSITE_SCHUR, it returns (in order) 2274 - the KSP used for the (1,1) block 2275 - the KSP used for the Schur complement (not the one used for the interior Schur solver) 2276 - the KSP used for the (1,1) block in the upper triangular factor (if different from that of the (1,1) block). 2277 2278 It returns a null array if the fieldsplit is not of type PC_COMPOSITE_SCHUR; in this case, you should use PCFieldSplitGetSubKSP(). 2279 2280 Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs. 2281 You can call PCFieldSplitSchurGetSubKSP(pc,n,PETSC_NULL_KSP,ierr) to determine how large the 2282 KSP array must be. 2283 2284 Level: advanced 2285 2286 .seealso: PCFIELDSPLIT 2287 @*/ 2288 PetscErrorCode PCFieldSplitSchurGetSubKSP(PC pc,PetscInt *n,KSP *subksp[]) 2289 { 2290 PetscErrorCode ierr; 2291 2292 PetscFunctionBegin; 2293 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2294 if (n) PetscValidIntPointer(n,2); 2295 ierr = PetscUseMethod(pc,"PCFieldSplitSchurGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));CHKERRQ(ierr); 2296 PetscFunctionReturn(0); 2297 } 2298 2299 /*@ 2300 PCFieldSplitSetSchurPre - Indicates from what operator the preconditioner is constructucted for the Schur complement. 2301 The default is the A11 matrix. 2302 2303 Collective on PC 2304 2305 Input Parameters: 2306 + pc - the preconditioner context 2307 . 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 2308 PC_FIELDSPLIT_SCHUR_PRE_SELFP, and PC_FIELDSPLIT_SCHUR_PRE_FULL 2309 - userpre - matrix to use for preconditioning, or NULL 2310 2311 Options Database: 2312 + -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is a11. See notes for meaning of various arguments 2313 - -fieldsplit_1_pc_type <pctype> - the preconditioner algorithm that is used to construct the preconditioner from the operator 2314 2315 Notes: 2316 $ If ptype is 2317 $ a11 - 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 - the preconditioner for the Schur complement is generated from the user provided matrix (pre argument 2323 $ to this function). 2324 $ selfp - 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 - 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