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