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