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