1 2 #include <petsc-private/pcimpl.h> /*I "petscpc.h" I*/ 3 #include <petscdm.h> 4 5 /* 6 There is a nice discussion of block preconditioners in 7 8 [El08] A taxonomy and comparison of parallel block multi-level preconditioners for the incompressible Navier–Stokes equations 9 Howard Elman, V.E. Howle, John Shadid, Robert Shuttleworth, Ray Tuminaro, Journal of Computational Physics 227 (2008) 1790--1808 10 http://chess.cs.umd.edu/~elman/papers/tax.pdf 11 */ 12 13 const char *const PCFieldSplitSchurPreTypes[] = {"SELF","A11","USER","PCFieldSplitSchurPreType","PC_FIELDSPLIT_SCHUR_PRE_",0}; 14 const char *const PCFieldSplitSchurFactTypes[] = {"DIAG","LOWER","UPPER","FULL","PCFieldSplitSchurFactType","PC_FIELDSPLIT_SCHUR_FACT_",0}; 15 16 typedef struct _PC_FieldSplitLink *PC_FieldSplitLink; 17 struct _PC_FieldSplitLink { 18 KSP ksp; 19 Vec x,y,z; 20 char *splitname; 21 PetscInt nfields; 22 PetscInt *fields,*fields_col; 23 VecScatter sctx; 24 IS is,is_col; 25 PC_FieldSplitLink next,previous; 26 }; 27 28 typedef struct { 29 PCCompositeType type; 30 PetscBool defaultsplit; /* Flag for a system with a set of 'k' scalar fields with the same layout (and bs = k) */ 31 PetscBool splitdefined; /* Flag is set after the splits have been defined, to prevent more splits from being added */ 32 PetscInt bs; /* Block size for IS and Mat structures */ 33 PetscInt nsplits; /* Number of field divisions defined */ 34 Vec *x,*y,w1,w2; 35 Mat *mat; /* The diagonal block for each split */ 36 Mat *pmat; /* The preconditioning diagonal block for each split */ 37 Mat *Afield; /* The rows of the matrix associated with each split */ 38 PetscBool issetup; 39 40 /* Only used when Schur complement preconditioning is used */ 41 Mat B; /* The (0,1) block */ 42 Mat C; /* The (1,0) block */ 43 Mat schur; /* The Schur complement S = A11 - A10 A00^{-1} A01, the KSP here, kspinner, is H_1 in [El08] */ 44 Mat schur_user; /* User-provided preconditioning matrix for the Schur complement */ 45 PCFieldSplitSchurPreType schurpre; /* Determines which preconditioning matrix is used for the Schur complement */ 46 PCFieldSplitSchurFactType schurfactorization; 47 KSP kspschur; /* The solver for S */ 48 KSP kspupper; /* The solver for A in the upper diagonal part of the factorization (H_2 in [El08]) */ 49 PC_FieldSplitLink head; 50 PetscBool reset; /* indicates PCReset() has been last called on this object, hack */ 51 PetscBool suboptionsset; /* Indicates that the KSPSetFromOptions() has been called on the sub-KSPs */ 52 PetscBool dm_splits; /* Whether to use DMCreateFieldDecomposition() whenever possible */ 53 } PC_FieldSplit; 54 55 /* 56 Notes: there is no particular reason that pmat, x, and y are stored as arrays in PC_FieldSplit instead of 57 inside PC_FieldSplitLink, just historical. If you want to be able to add new fields after already using the 58 PC you could change this. 59 */ 60 61 /* This helper is so that setting a user-provided preconditioning matrix is orthogonal to choosing to use it. This way the 62 * application-provided FormJacobian can provide this matrix without interfering with the user's (command-line) choices. */ 63 static Mat FieldSplitSchurPre(PC_FieldSplit *jac) 64 { 65 switch (jac->schurpre) { 66 case PC_FIELDSPLIT_SCHUR_PRE_SELF: return jac->schur; 67 case PC_FIELDSPLIT_SCHUR_PRE_A11: return jac->pmat[1]; 68 case PC_FIELDSPLIT_SCHUR_PRE_USER: /* Use a user-provided matrix if it is given, otherwise diagonal block */ 69 default: 70 return jac->schur_user ? jac->schur_user : jac->pmat[1]; 71 } 72 } 73 74 75 #include <petscdraw.h> 76 #undef __FUNCT__ 77 #define __FUNCT__ "PCView_FieldSplit" 78 static PetscErrorCode PCView_FieldSplit(PC pc,PetscViewer viewer) 79 { 80 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 81 PetscErrorCode ierr; 82 PetscBool iascii,isdraw; 83 PetscInt i,j; 84 PC_FieldSplitLink ilink = jac->head; 85 86 PetscFunctionBegin; 87 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 88 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 89 if (iascii) { 90 if (jac->bs > 0) { 91 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D, blocksize = %D\n",PCCompositeTypes[jac->type],jac->nsplits,jac->bs);CHKERRQ(ierr); 92 } else { 93 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D\n",PCCompositeTypes[jac->type],jac->nsplits);CHKERRQ(ierr); 94 } 95 if (pc->useAmat) { 96 ierr = PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for blocks\n");CHKERRQ(ierr); 97 } 98 ierr = PetscViewerASCIIPrintf(viewer," Solver info for each split is in the following KSP objects:\n");CHKERRQ(ierr); 99 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 100 for (i=0; i<jac->nsplits; i++) { 101 if (ilink->fields) { 102 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);CHKERRQ(ierr); 103 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 104 for (j=0; j<ilink->nfields; j++) { 105 if (j > 0) { 106 ierr = PetscViewerASCIIPrintf(viewer,",");CHKERRQ(ierr); 107 } 108 ierr = PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);CHKERRQ(ierr); 109 } 110 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 111 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 112 } else { 113 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i);CHKERRQ(ierr); 114 } 115 ierr = KSPView(ilink->ksp,viewer);CHKERRQ(ierr); 116 ilink = ilink->next; 117 } 118 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 119 } 120 121 if (isdraw) { 122 PetscDraw draw; 123 PetscReal x,y,w,wd; 124 125 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 126 ierr = PetscDrawGetCurrentPoint(draw,&x,&y);CHKERRQ(ierr); 127 w = 2*PetscMin(1.0 - x,x); 128 wd = w/(jac->nsplits + 1); 129 x = x - wd*(jac->nsplits-1)/2.0; 130 for (i=0; i<jac->nsplits; i++) { 131 ierr = PetscDrawPushCurrentPoint(draw,x,y);CHKERRQ(ierr); 132 ierr = KSPView(ilink->ksp,viewer);CHKERRQ(ierr); 133 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr); 134 x += wd; 135 ilink = ilink->next; 136 } 137 } 138 PetscFunctionReturn(0); 139 } 140 141 #undef __FUNCT__ 142 #define __FUNCT__ "PCView_FieldSplit_Schur" 143 static PetscErrorCode PCView_FieldSplit_Schur(PC pc,PetscViewer viewer) 144 { 145 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 146 PetscErrorCode ierr; 147 PetscBool iascii,isdraw; 148 PetscInt i,j; 149 PC_FieldSplitLink ilink = jac->head; 150 151 PetscFunctionBegin; 152 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 153 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 154 if (iascii) { 155 if (jac->bs > 0) { 156 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit with Schur preconditioner, blocksize = %D, factorization %s\n",jac->bs,PCFieldSplitSchurFactTypes[jac->schurfactorization]);CHKERRQ(ierr); 157 } else { 158 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit with Schur preconditioner, factorization %s\n",PCFieldSplitSchurFactTypes[jac->schurfactorization]);CHKERRQ(ierr); 159 } 160 if (pc->useAmat) { 161 ierr = PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for blocks\n");CHKERRQ(ierr); 162 } 163 switch (jac->schurpre) { 164 case PC_FIELDSPLIT_SCHUR_PRE_SELF: 165 ierr = PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from S itself\n");CHKERRQ(ierr);break; 166 case PC_FIELDSPLIT_SCHUR_PRE_A11: 167 ierr = PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from A11\n");CHKERRQ(ierr);break; 168 case PC_FIELDSPLIT_SCHUR_PRE_USER: 169 if (jac->schur_user) { 170 ierr = PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from user provided matrix\n");CHKERRQ(ierr); 171 } else { 172 ierr = PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from A11\n");CHKERRQ(ierr); 173 } 174 break; 175 default: 176 SETERRQ1(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Invalid Schur preconditioning type: %d", jac->schurpre); 177 } 178 ierr = PetscViewerASCIIPrintf(viewer," Split info:\n");CHKERRQ(ierr); 179 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 180 for (i=0; i<jac->nsplits; i++) { 181 if (ilink->fields) { 182 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);CHKERRQ(ierr); 183 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 184 for (j=0; j<ilink->nfields; j++) { 185 if (j > 0) { 186 ierr = PetscViewerASCIIPrintf(viewer,",");CHKERRQ(ierr); 187 } 188 ierr = PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);CHKERRQ(ierr); 189 } 190 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 191 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 192 } else { 193 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i);CHKERRQ(ierr); 194 } 195 ilink = ilink->next; 196 } 197 ierr = PetscViewerASCIIPrintf(viewer,"KSP solver for A00 block \n");CHKERRQ(ierr); 198 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 199 ierr = KSPView(jac->head->ksp,viewer);CHKERRQ(ierr); 200 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 201 if (jac->kspupper != jac->head->ksp) { 202 ierr = PetscViewerASCIIPrintf(viewer,"KSP solver for upper A00 in upper triangular factor \n");CHKERRQ(ierr); 203 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 204 ierr = KSPView(jac->kspupper,viewer);CHKERRQ(ierr); 205 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 206 } 207 ierr = PetscViewerASCIIPrintf(viewer,"KSP solver for S = A11 - A10 inv(A00) A01 \n");CHKERRQ(ierr); 208 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 209 if (jac->kspschur) { 210 ierr = KSPView(jac->kspschur,viewer);CHKERRQ(ierr); 211 } else { 212 ierr = PetscViewerASCIIPrintf(viewer," not yet available\n");CHKERRQ(ierr); 213 } 214 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 215 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 216 } else if (isdraw) { 217 PetscDraw draw; 218 PetscReal x,y,w,wd,h; 219 PetscInt cnt = 2; 220 char str[32]; 221 222 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 223 ierr = PetscDrawGetCurrentPoint(draw,&x,&y);CHKERRQ(ierr); 224 if (jac->kspupper != jac->head->ksp) cnt++; 225 w = 2*PetscMin(1.0 - x,x); 226 wd = w/(cnt + 1); 227 228 ierr = PetscSNPrintf(str,32,"Schur fact. %s",PCFieldSplitSchurFactTypes[jac->schurfactorization]);CHKERRQ(ierr); 229 ierr = PetscDrawBoxedString(draw,x,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);CHKERRQ(ierr); 230 y -= h; 231 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_USER && !jac->schur_user) { 232 ierr = PetscSNPrintf(str,32,"Prec. for Schur from %s",PCFieldSplitSchurPreTypes[PC_FIELDSPLIT_SCHUR_PRE_A11]);CHKERRQ(ierr); 233 } else { 234 ierr = PetscSNPrintf(str,32,"Prec. for Schur from %s",PCFieldSplitSchurPreTypes[jac->schurpre]);CHKERRQ(ierr); 235 } 236 ierr = PetscDrawBoxedString(draw,x+wd*(cnt-1)/2.0,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);CHKERRQ(ierr); 237 y -= h; 238 x = x - wd*(cnt-1)/2.0; 239 240 ierr = PetscDrawPushCurrentPoint(draw,x,y);CHKERRQ(ierr); 241 ierr = KSPView(jac->head->ksp,viewer);CHKERRQ(ierr); 242 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr); 243 if (jac->kspupper != jac->head->ksp) { 244 x += wd; 245 ierr = PetscDrawPushCurrentPoint(draw,x,y);CHKERRQ(ierr); 246 ierr = KSPView(jac->kspupper,viewer);CHKERRQ(ierr); 247 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr); 248 } 249 x += wd; 250 ierr = PetscDrawPushCurrentPoint(draw,x,y);CHKERRQ(ierr); 251 ierr = KSPView(jac->kspschur,viewer);CHKERRQ(ierr); 252 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr); 253 } 254 PetscFunctionReturn(0); 255 } 256 257 #undef __FUNCT__ 258 #define __FUNCT__ "PCFieldSplitSetRuntimeSplits_Private" 259 /* Precondition: jac->bs is set to a meaningful value */ 260 static PetscErrorCode PCFieldSplitSetRuntimeSplits_Private(PC pc) 261 { 262 PetscErrorCode ierr; 263 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 264 PetscInt i,nfields,*ifields,nfields_col,*ifields_col; 265 PetscBool flg,flg_col; 266 char optionname[128],splitname[8],optionname_col[128]; 267 268 PetscFunctionBegin; 269 ierr = PetscMalloc(jac->bs*sizeof(PetscInt),&ifields);CHKERRQ(ierr); 270 ierr = PetscMalloc(jac->bs*sizeof(PetscInt),&ifields_col);CHKERRQ(ierr); 271 for (i=0,flg=PETSC_TRUE;; i++) { 272 ierr = PetscSNPrintf(splitname,sizeof(splitname),"%D",i);CHKERRQ(ierr); 273 ierr = PetscSNPrintf(optionname,sizeof(optionname),"-pc_fieldsplit_%D_fields",i);CHKERRQ(ierr); 274 ierr = PetscSNPrintf(optionname_col,sizeof(optionname_col),"-pc_fieldsplit_%D_fields_col",i);CHKERRQ(ierr); 275 nfields = jac->bs; 276 nfields_col = jac->bs; 277 ierr = PetscOptionsGetIntArray(((PetscObject)pc)->prefix,optionname,ifields,&nfields,&flg);CHKERRQ(ierr); 278 ierr = PetscOptionsGetIntArray(((PetscObject)pc)->prefix,optionname_col,ifields_col,&nfields_col,&flg_col);CHKERRQ(ierr); 279 if (!flg) break; 280 else if (flg && !flg_col) { 281 if (!nfields) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot list zero fields"); 282 ierr = PCFieldSplitSetFields(pc,splitname,nfields,ifields,ifields);CHKERRQ(ierr); 283 } else { 284 if (!nfields || !nfields_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot list zero fields"); 285 if (nfields != nfields_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Number of row and column fields must match"); 286 ierr = PCFieldSplitSetFields(pc,splitname,nfields,ifields,ifields_col);CHKERRQ(ierr); 287 } 288 } 289 if (i > 0) { 290 /* Makes command-line setting of splits take precedence over setting them in code. 291 Otherwise subsequent calls to PCFieldSplitSetIS() or PCFieldSplitSetFields() would 292 create new splits, which would probably not be what the user wanted. */ 293 jac->splitdefined = PETSC_TRUE; 294 } 295 ierr = PetscFree(ifields);CHKERRQ(ierr); 296 ierr = PetscFree(ifields_col);CHKERRQ(ierr); 297 PetscFunctionReturn(0); 298 } 299 300 #undef __FUNCT__ 301 #define __FUNCT__ "PCFieldSplitSetDefaults" 302 static PetscErrorCode PCFieldSplitSetDefaults(PC pc) 303 { 304 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 305 PetscErrorCode ierr; 306 PC_FieldSplitLink ilink = jac->head; 307 PetscBool fieldsplit_default = PETSC_FALSE,stokes = PETSC_FALSE; 308 PetscInt i; 309 310 PetscFunctionBegin; 311 /* 312 Kinda messy, but at least this now uses DMCreateFieldDecomposition() even with jac->reset. 313 Should probably be rewritten. 314 */ 315 if (!ilink || jac->reset) { 316 ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,NULL);CHKERRQ(ierr); 317 if (pc->dm && jac->dm_splits && !stokes) { 318 PetscInt numFields, f, i, j; 319 char **fieldNames; 320 IS *fields; 321 DM *dms; 322 DM subdm[128]; 323 PetscBool flg; 324 325 ierr = DMCreateFieldDecomposition(pc->dm, &numFields, &fieldNames, &fields, &dms);CHKERRQ(ierr); 326 /* Allow the user to prescribe the splits */ 327 for (i = 0, flg = PETSC_TRUE;; i++) { 328 PetscInt ifields[128]; 329 IS compField; 330 char optionname[128], splitname[8]; 331 PetscInt nfields = numFields; 332 333 ierr = PetscSNPrintf(optionname, sizeof(optionname), "-pc_fieldsplit_%D_fields", i);CHKERRQ(ierr); 334 ierr = PetscOptionsGetIntArray(((PetscObject) pc)->prefix, optionname, ifields, &nfields, &flg);CHKERRQ(ierr); 335 if (!flg) break; 336 if (numFields > 128) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot currently support %d > 128 fields", numFields); 337 ierr = DMCreateSubDM(pc->dm, nfields, ifields, &compField, &subdm[i]);CHKERRQ(ierr); 338 if (nfields == 1) { 339 ierr = PCFieldSplitSetIS(pc, fieldNames[ifields[0]], compField);CHKERRQ(ierr); 340 /* ierr = PetscPrintf(PetscObjectComm((PetscObject)pc), "%s Field Indices:", fieldNames[ifields[0]]);CHKERRQ(ierr); 341 ierr = ISView(compField, NULL);CHKERRQ(ierr); */ 342 } else { 343 ierr = PetscSNPrintf(splitname, sizeof(splitname), "%D", i);CHKERRQ(ierr); 344 ierr = PCFieldSplitSetIS(pc, splitname, compField);CHKERRQ(ierr); 345 /* ierr = PetscPrintf(PetscObjectComm((PetscObject)pc), "%s Field Indices:", splitname);CHKERRQ(ierr); 346 ierr = ISView(compField, NULL);CHKERRQ(ierr); */ 347 } 348 ierr = ISDestroy(&compField);CHKERRQ(ierr); 349 for (j = 0; j < nfields; ++j) { 350 f = ifields[j]; 351 ierr = PetscFree(fieldNames[f]);CHKERRQ(ierr); 352 ierr = ISDestroy(&fields[f]);CHKERRQ(ierr); 353 } 354 } 355 if (i == 0) { 356 for (f = 0; f < numFields; ++f) { 357 ierr = PCFieldSplitSetIS(pc, fieldNames[f], fields[f]);CHKERRQ(ierr); 358 ierr = PetscFree(fieldNames[f]);CHKERRQ(ierr); 359 ierr = ISDestroy(&fields[f]);CHKERRQ(ierr); 360 } 361 } else { 362 for (j=0; j<numFields; j++) { 363 ierr = DMDestroy(dms+j);CHKERRQ(ierr); 364 } 365 ierr = PetscFree(dms);CHKERRQ(ierr); 366 ierr = PetscMalloc(i * sizeof(DM), &dms);CHKERRQ(ierr); 367 for (j = 0; j < i; ++j) dms[j] = subdm[j]; 368 } 369 ierr = PetscFree(fieldNames);CHKERRQ(ierr); 370 ierr = PetscFree(fields);CHKERRQ(ierr); 371 if (dms) { 372 ierr = PetscInfo(pc, "Setting up physics based fieldsplit preconditioner using the embedded DM\n");CHKERRQ(ierr); 373 for (ilink = jac->head, i = 0; ilink; ilink = ilink->next, ++i) { 374 const char *prefix; 375 ierr = PetscObjectGetOptionsPrefix((PetscObject)(ilink->ksp),&prefix);CHKERRQ(ierr); 376 ierr = PetscObjectSetOptionsPrefix((PetscObject)(dms[i]), prefix);CHKERRQ(ierr); 377 ierr = KSPSetDM(ilink->ksp, dms[i]);CHKERRQ(ierr); 378 ierr = KSPSetDMActive(ilink->ksp, PETSC_FALSE);CHKERRQ(ierr); 379 ierr = PetscObjectIncrementTabLevel((PetscObject)dms[i],(PetscObject)ilink->ksp,0);CHKERRQ(ierr); 380 ierr = DMDestroy(&dms[i]);CHKERRQ(ierr); 381 } 382 ierr = PetscFree(dms);CHKERRQ(ierr); 383 } 384 } else { 385 if (jac->bs <= 0) { 386 if (pc->pmat) { 387 ierr = MatGetBlockSize(pc->pmat,&jac->bs);CHKERRQ(ierr); 388 } else jac->bs = 1; 389 } 390 391 if (stokes) { 392 IS zerodiags,rest; 393 PetscInt nmin,nmax; 394 395 ierr = MatGetOwnershipRange(pc->mat,&nmin,&nmax);CHKERRQ(ierr); 396 ierr = MatFindZeroDiagonals(pc->mat,&zerodiags);CHKERRQ(ierr); 397 ierr = ISComplement(zerodiags,nmin,nmax,&rest);CHKERRQ(ierr); 398 if (jac->reset) { 399 jac->head->is = rest; 400 jac->head->next->is = zerodiags; 401 } else { 402 ierr = PCFieldSplitSetIS(pc,"0",rest);CHKERRQ(ierr); 403 ierr = PCFieldSplitSetIS(pc,"1",zerodiags);CHKERRQ(ierr); 404 } 405 ierr = ISDestroy(&zerodiags);CHKERRQ(ierr); 406 ierr = ISDestroy(&rest);CHKERRQ(ierr); 407 } else { 408 if (jac->reset) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cases not yet handled when PCReset() was used"); 409 ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_default",&fieldsplit_default,NULL);CHKERRQ(ierr); 410 if (!fieldsplit_default) { 411 /* Allow user to set fields from command line, if bs was known at the time of PCSetFromOptions_FieldSplit() 412 then it is set there. This is not ideal because we should only have options set in XXSetFromOptions(). */ 413 ierr = PCFieldSplitSetRuntimeSplits_Private(pc);CHKERRQ(ierr); 414 if (jac->splitdefined) {ierr = PetscInfo(pc,"Splits defined using the options database\n");CHKERRQ(ierr);} 415 } 416 if (fieldsplit_default || !jac->splitdefined) { 417 ierr = PetscInfo(pc,"Using default splitting of fields\n");CHKERRQ(ierr); 418 for (i=0; i<jac->bs; i++) { 419 char splitname[8]; 420 ierr = PetscSNPrintf(splitname,sizeof(splitname),"%D",i);CHKERRQ(ierr); 421 ierr = PCFieldSplitSetFields(pc,splitname,1,&i,&i);CHKERRQ(ierr); 422 } 423 jac->defaultsplit = PETSC_TRUE; 424 } 425 } 426 } 427 } else if (jac->nsplits == 1) { 428 if (ilink->is) { 429 IS is2; 430 PetscInt nmin,nmax; 431 432 ierr = MatGetOwnershipRange(pc->mat,&nmin,&nmax);CHKERRQ(ierr); 433 ierr = ISComplement(ilink->is,nmin,nmax,&is2);CHKERRQ(ierr); 434 ierr = PCFieldSplitSetIS(pc,"1",is2);CHKERRQ(ierr); 435 ierr = ISDestroy(&is2);CHKERRQ(ierr); 436 } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Must provide at least two sets of fields to PCFieldSplit()"); 437 } 438 439 440 if (jac->nsplits < 2) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Unhandled case, must have at least two fields, not %d", jac->nsplits); 441 PetscFunctionReturn(0); 442 } 443 444 PETSC_EXTERN PetscErrorCode PetscOptionsFindPairPrefix_Private(const char pre[], const char name[], char *value[], PetscBool *flg); 445 446 #undef __FUNCT__ 447 #define __FUNCT__ "PCSetUp_FieldSplit" 448 static PetscErrorCode PCSetUp_FieldSplit(PC pc) 449 { 450 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 451 PetscErrorCode ierr; 452 PC_FieldSplitLink ilink; 453 PetscInt i,nsplit; 454 MatStructure flag = pc->flag; 455 PetscBool sorted, sorted_col; 456 457 PetscFunctionBegin; 458 ierr = PCFieldSplitSetDefaults(pc);CHKERRQ(ierr); 459 nsplit = jac->nsplits; 460 ilink = jac->head; 461 462 /* get the matrices for each split */ 463 if (!jac->issetup) { 464 PetscInt rstart,rend,nslots,bs; 465 466 jac->issetup = PETSC_TRUE; 467 468 /* This is done here instead of in PCFieldSplitSetFields() because may not have matrix at that point */ 469 if (jac->defaultsplit || !ilink->is) { 470 if (jac->bs <= 0) jac->bs = nsplit; 471 } 472 bs = jac->bs; 473 ierr = MatGetOwnershipRange(pc->pmat,&rstart,&rend);CHKERRQ(ierr); 474 nslots = (rend - rstart)/bs; 475 for (i=0; i<nsplit; i++) { 476 if (jac->defaultsplit) { 477 ierr = ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+i,nsplit,&ilink->is);CHKERRQ(ierr); 478 ierr = ISDuplicate(ilink->is,&ilink->is_col);CHKERRQ(ierr); 479 } else if (!ilink->is) { 480 if (ilink->nfields > 1) { 481 PetscInt *ii,*jj,j,k,nfields = ilink->nfields,*fields = ilink->fields,*fields_col = ilink->fields_col; 482 ierr = PetscMalloc(ilink->nfields*nslots*sizeof(PetscInt),&ii);CHKERRQ(ierr); 483 ierr = PetscMalloc(ilink->nfields*nslots*sizeof(PetscInt),&jj);CHKERRQ(ierr); 484 for (j=0; j<nslots; j++) { 485 for (k=0; k<nfields; k++) { 486 ii[nfields*j + k] = rstart + bs*j + fields[k]; 487 jj[nfields*j + k] = rstart + bs*j + fields_col[k]; 488 } 489 } 490 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,ii,PETSC_OWN_POINTER,&ilink->is);CHKERRQ(ierr); 491 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,jj,PETSC_OWN_POINTER,&ilink->is_col);CHKERRQ(ierr); 492 } else { 493 ierr = ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields[0],bs,&ilink->is);CHKERRQ(ierr); 494 ierr = ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields_col[0],bs,&ilink->is_col);CHKERRQ(ierr); 495 } 496 } 497 ierr = ISSorted(ilink->is,&sorted);CHKERRQ(ierr); 498 if (ilink->is_col) { ierr = ISSorted(ilink->is_col,&sorted_col);CHKERRQ(ierr); } 499 if (!sorted || !sorted_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Fields must be sorted when creating split"); 500 ilink = ilink->next; 501 } 502 } 503 504 ilink = jac->head; 505 if (!jac->pmat) { 506 Vec xtmp; 507 508 ierr = MatGetVecs(pc->pmat,&xtmp,NULL);CHKERRQ(ierr); 509 ierr = PetscMalloc(nsplit*sizeof(Mat),&jac->pmat);CHKERRQ(ierr); 510 ierr = PetscMalloc2(nsplit,Vec,&jac->x,nsplit,Vec,&jac->y);CHKERRQ(ierr); 511 for (i=0; i<nsplit; i++) { 512 MatNullSpace sp; 513 514 /* Check for preconditioning matrix attached to IS */ 515 ierr = PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &jac->pmat[i]);CHKERRQ(ierr); 516 if (jac->pmat[i]) { 517 ierr = PetscObjectReference((PetscObject) jac->pmat[i]);CHKERRQ(ierr); 518 if (jac->type == PC_COMPOSITE_SCHUR) { 519 jac->schur_user = jac->pmat[i]; 520 521 ierr = PetscObjectReference((PetscObject) jac->schur_user);CHKERRQ(ierr); 522 } 523 } else { 524 ierr = MatGetSubMatrix(pc->pmat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->pmat[i]);CHKERRQ(ierr); 525 } 526 /* create work vectors for each split */ 527 ierr = MatGetVecs(jac->pmat[i],&jac->x[i],&jac->y[i]);CHKERRQ(ierr); 528 ilink->x = jac->x[i]; ilink->y = jac->y[i]; ilink->z = NULL; 529 /* compute scatter contexts needed by multiplicative versions and non-default splits */ 530 ierr = VecScatterCreate(xtmp,ilink->is,jac->x[i],NULL,&ilink->sctx);CHKERRQ(ierr); 531 /* Check for null space attached to IS */ 532 ierr = PetscObjectQuery((PetscObject) ilink->is, "nullspace", (PetscObject*) &sp);CHKERRQ(ierr); 533 if (sp) { 534 ierr = MatSetNullSpace(jac->pmat[i], sp);CHKERRQ(ierr); 535 } 536 ierr = PetscObjectQuery((PetscObject) ilink->is, "nearnullspace", (PetscObject*) &sp);CHKERRQ(ierr); 537 if (sp) { 538 ierr = MatSetNearNullSpace(jac->pmat[i], sp);CHKERRQ(ierr); 539 } 540 ilink = ilink->next; 541 } 542 ierr = VecDestroy(&xtmp);CHKERRQ(ierr); 543 } else { 544 for (i=0; i<nsplit; i++) { 545 Mat pmat; 546 547 /* Check for preconditioning matrix attached to IS */ 548 ierr = PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &pmat);CHKERRQ(ierr); 549 if (!pmat) { 550 ierr = MatGetSubMatrix(pc->pmat,ilink->is,ilink->is_col,MAT_REUSE_MATRIX,&jac->pmat[i]);CHKERRQ(ierr); 551 } 552 ilink = ilink->next; 553 } 554 } 555 if (pc->useAmat) { 556 ilink = jac->head; 557 if (!jac->mat) { 558 ierr = PetscMalloc(nsplit*sizeof(Mat),&jac->mat);CHKERRQ(ierr); 559 for (i=0; i<nsplit; i++) { 560 ierr = MatGetSubMatrix(pc->mat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->mat[i]);CHKERRQ(ierr); 561 ilink = ilink->next; 562 } 563 } else { 564 for (i=0; i<nsplit; i++) { 565 if (jac->mat[i]) {ierr = MatGetSubMatrix(pc->mat,ilink->is,ilink->is_col,MAT_REUSE_MATRIX,&jac->mat[i]);CHKERRQ(ierr);} 566 ilink = ilink->next; 567 } 568 } 569 } else { 570 jac->mat = jac->pmat; 571 } 572 573 if (jac->type != PC_COMPOSITE_ADDITIVE && jac->type != PC_COMPOSITE_SCHUR) { 574 /* extract the rows of the matrix associated with each field: used for efficient computation of residual inside algorithm */ 575 ilink = jac->head; 576 if (!jac->Afield) { 577 ierr = PetscMalloc(nsplit*sizeof(Mat),&jac->Afield);CHKERRQ(ierr); 578 for (i=0; i<nsplit; i++) { 579 ierr = MatGetSubMatrix(pc->mat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);CHKERRQ(ierr); 580 ilink = ilink->next; 581 } 582 } else { 583 for (i=0; i<nsplit; i++) { 584 ierr = MatGetSubMatrix(pc->mat,ilink->is,NULL,MAT_REUSE_MATRIX,&jac->Afield[i]);CHKERRQ(ierr); 585 ilink = ilink->next; 586 } 587 } 588 } 589 590 if (jac->type == PC_COMPOSITE_SCHUR) { 591 IS ccis; 592 PetscInt rstart,rend; 593 char lscname[256]; 594 PetscObject LSC_L; 595 596 if (nsplit != 2) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"To use Schur complement preconditioner you must have exactly 2 fields"); 597 598 /* When extracting off-diagonal submatrices, we take complements from this range */ 599 ierr = MatGetOwnershipRangeColumn(pc->mat,&rstart,&rend);CHKERRQ(ierr); 600 601 /* need to handle case when one is resetting up the preconditioner */ 602 if (jac->schur) { 603 KSP kspA = jac->head->ksp, kspInner = NULL, kspUpper = jac->kspupper; 604 605 ierr = MatSchurComplementGetKSP(jac->schur, &kspInner);CHKERRQ(ierr); 606 ilink = jac->head; 607 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 608 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->B);CHKERRQ(ierr); 609 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 610 ilink = ilink->next; 611 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 612 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->C);CHKERRQ(ierr); 613 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 614 ierr = MatSchurComplementUpdate(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1],pc->flag);CHKERRQ(ierr); 615 if (kspA != kspInner) { 616 ierr = KSPSetOperators(kspA,jac->mat[0],jac->pmat[0],pc->flag);CHKERRQ(ierr); 617 } 618 if (kspUpper != kspA) { 619 ierr = KSPSetOperators(kspUpper,jac->mat[0],jac->pmat[0],pc->flag);CHKERRQ(ierr); 620 } 621 ierr = KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac),pc->flag);CHKERRQ(ierr); 622 } else { 623 const char *Dprefix; 624 char schurprefix[256]; 625 char schurtestoption[256]; 626 MatNullSpace sp; 627 PetscBool flg; 628 629 /* extract the A01 and A10 matrices */ 630 ilink = jac->head; 631 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 632 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);CHKERRQ(ierr); 633 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 634 ilink = ilink->next; 635 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 636 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);CHKERRQ(ierr); 637 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 638 639 /* Use mat[0] (diagonal block of Amat) preconditioned by pmat[0] to define Schur complement */ 640 ierr = MatCreate(((PetscObject)jac->mat[0])->comm,&jac->schur);CHKERRQ(ierr); 641 ierr = MatSetType(jac->schur,MATSCHURCOMPLEMENT);CHKERRQ(ierr); 642 ierr = MatSchurComplementSet(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);CHKERRQ(ierr); 643 644 ierr = MatGetNullSpace(jac->pmat[1], &sp);CHKERRQ(ierr); 645 if (sp) { 646 ierr = MatSetNullSpace(jac->schur, sp);CHKERRQ(ierr); 647 } 648 649 ierr = PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_inner_", ilink->splitname);CHKERRQ(ierr); 650 ierr = PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);CHKERRQ(ierr); 651 if (flg) { 652 DM dmInner; 653 KSP kspInner; 654 655 ierr = MatSchurComplementGetKSP(jac->schur, &kspInner);CHKERRQ(ierr); 656 ierr = PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_inner_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);CHKERRQ(ierr); 657 /* Indent this deeper to emphasize the "inner" nature of this solver. */ 658 ierr = PetscObjectIncrementTabLevel((PetscObject)kspInner, (PetscObject) pc, 2);CHKERRQ(ierr); 659 ierr = KSPSetOptionsPrefix(kspInner, schurprefix);CHKERRQ(ierr); 660 661 /* Set DM for new solver */ 662 ierr = KSPGetDM(jac->head->ksp, &dmInner);CHKERRQ(ierr); 663 ierr = KSPSetDM(kspInner, dmInner);CHKERRQ(ierr); 664 ierr = KSPSetDMActive(kspInner, PETSC_FALSE);CHKERRQ(ierr); 665 } else { 666 /* Use the outer solver for the inner solve, but revert the KSPPREONLY from PCFieldSplitSetFields_FieldSplit or 667 * PCFieldSplitSetIS_FieldSplit. We don't want KSPPREONLY because it makes the Schur complement inexact, 668 * preventing Schur complement reduction to be an accurate solve. Usually when an iterative solver is used for 669 * S = D - C A_inner^{-1} B, we expect S to be defined using an accurate definition of A_inner^{-1}, so we make 670 * GMRES the default. Note that it is also common to use PREONLY for S, in which case S may not be used 671 * directly, and the user is responsible for setting an inexact method for fieldsplit's A^{-1}. */ 672 ierr = KSPSetType(jac->head->ksp,KSPGMRES);CHKERRQ(ierr); 673 ierr = MatSchurComplementSetKSP(jac->schur,jac->head->ksp);CHKERRQ(ierr); 674 } 675 ierr = KSPSetOperators(jac->head->ksp,jac->mat[0],jac->pmat[0],flag);CHKERRQ(ierr); 676 ierr = KSPSetFromOptions(jac->head->ksp);CHKERRQ(ierr); 677 ierr = MatSetFromOptions(jac->schur);CHKERRQ(ierr); 678 679 ierr = PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_upper_", ilink->splitname);CHKERRQ(ierr); 680 ierr = PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);CHKERRQ(ierr); 681 if (flg) { 682 DM dmInner; 683 684 ierr = PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_upper_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);CHKERRQ(ierr); 685 ierr = KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspupper);CHKERRQ(ierr); 686 ierr = KSPSetOptionsPrefix(jac->kspupper, schurprefix);CHKERRQ(ierr); 687 ierr = KSPGetDM(jac->head->ksp, &dmInner);CHKERRQ(ierr); 688 ierr = KSPSetDM(jac->kspupper, dmInner);CHKERRQ(ierr); 689 ierr = KSPSetDMActive(jac->kspupper, PETSC_FALSE);CHKERRQ(ierr); 690 ierr = KSPSetFromOptions(jac->kspupper);CHKERRQ(ierr); 691 ierr = KSPSetOperators(jac->kspupper,jac->mat[0],jac->pmat[0],flag);CHKERRQ(ierr); 692 ierr = VecDuplicate(jac->head->x, &jac->head->z);CHKERRQ(ierr); 693 } else { 694 jac->kspupper = jac->head->ksp; 695 ierr = PetscObjectReference((PetscObject) jac->head->ksp);CHKERRQ(ierr); 696 } 697 698 ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&jac->kspschur);CHKERRQ(ierr); 699 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->kspschur);CHKERRQ(ierr); 700 ierr = PetscObjectIncrementTabLevel((PetscObject)jac->kspschur,(PetscObject)pc,1);CHKERRQ(ierr); 701 ierr = KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac),DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 702 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) { 703 PC pcschur; 704 ierr = KSPGetPC(jac->kspschur,&pcschur);CHKERRQ(ierr); 705 ierr = PCSetType(pcschur,PCNONE);CHKERRQ(ierr); 706 /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */ 707 } 708 ierr = KSPGetOptionsPrefix(jac->head->next->ksp, &Dprefix);CHKERRQ(ierr); 709 ierr = KSPSetOptionsPrefix(jac->kspschur, Dprefix);CHKERRQ(ierr); 710 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 711 /* need to call this every time, since the jac->kspschur is freshly created, otherwise its options never get set */ 712 ierr = KSPSetFromOptions(jac->kspschur);CHKERRQ(ierr); 713 } 714 715 /* HACK: special support to forward L and Lp matrices that might be used by PCLSC */ 716 ierr = PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_L",ilink->splitname);CHKERRQ(ierr); 717 ierr = PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr); 718 if (!LSC_L) {ierr = PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr);} 719 if (LSC_L) {ierr = PetscObjectCompose((PetscObject)jac->schur,"LSC_L",(PetscObject)LSC_L);CHKERRQ(ierr);} 720 ierr = PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_Lp",ilink->splitname);CHKERRQ(ierr); 721 ierr = PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr); 722 if (!LSC_L) {ierr = PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr);} 723 if (LSC_L) {ierr = PetscObjectCompose((PetscObject)jac->schur,"LSC_Lp",(PetscObject)LSC_L);CHKERRQ(ierr);} 724 } else { 725 /* set up the individual splits' PCs */ 726 i = 0; 727 ilink = jac->head; 728 while (ilink) { 729 ierr = KSPSetOperators(ilink->ksp,jac->mat[i],jac->pmat[i],flag);CHKERRQ(ierr); 730 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 731 if (!jac->suboptionsset) {ierr = KSPSetFromOptions(ilink->ksp);CHKERRQ(ierr);} 732 i++; 733 ilink = ilink->next; 734 } 735 } 736 737 jac->suboptionsset = PETSC_TRUE; 738 PetscFunctionReturn(0); 739 } 740 741 #define FieldSplitSplitSolveAdd(ilink,xx,yy) \ 742 (VecScatterBegin(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 743 VecScatterEnd(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 744 KSPSolve(ilink->ksp,ilink->x,ilink->y) || \ 745 VecScatterBegin(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE) || \ 746 VecScatterEnd(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE)) 747 748 #undef __FUNCT__ 749 #define __FUNCT__ "PCApply_FieldSplit_Schur" 750 static PetscErrorCode PCApply_FieldSplit_Schur(PC pc,Vec x,Vec y) 751 { 752 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 753 PetscErrorCode ierr; 754 PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next; 755 KSP kspA = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper; 756 757 PetscFunctionBegin; 758 switch (jac->schurfactorization) { 759 case PC_FIELDSPLIT_SCHUR_FACT_DIAG: 760 /* [A00 0; 0 -S], positive definite, suitable for MINRES */ 761 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 762 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 763 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 764 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 765 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 766 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 767 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 768 ierr = VecScale(ilinkD->y,-1.);CHKERRQ(ierr); 769 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 770 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 771 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 772 break; 773 case PC_FIELDSPLIT_SCHUR_FACT_LOWER: 774 /* [A00 0; A10 S], suitable for left preconditioning */ 775 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 776 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 777 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 778 ierr = MatMult(jac->C,ilinkA->y,ilinkD->x);CHKERRQ(ierr); 779 ierr = VecScale(ilinkD->x,-1.);CHKERRQ(ierr); 780 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 781 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 782 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 783 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 784 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 785 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 786 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 787 break; 788 case PC_FIELDSPLIT_SCHUR_FACT_UPPER: 789 /* [A00 A01; 0 S], suitable for right preconditioning */ 790 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 791 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 792 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 793 ierr = MatMult(jac->B,ilinkD->y,ilinkA->x);CHKERRQ(ierr); 794 ierr = VecScale(ilinkA->x,-1.);CHKERRQ(ierr); 795 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 796 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 797 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 798 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 799 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 800 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 801 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 802 break; 803 case PC_FIELDSPLIT_SCHUR_FACT_FULL: 804 /* [1 0; A10 A00^{-1} 1] [A00 0; 0 S] [1 A00^{-1}A01; 0 1], an exact solve if applied exactly, needs one extra solve with A */ 805 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 806 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 807 ierr = KSPSolve(kspLower,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 808 ierr = MatMult(jac->C,ilinkA->y,ilinkD->x);CHKERRQ(ierr); 809 ierr = VecScale(ilinkD->x,-1.0);CHKERRQ(ierr); 810 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 811 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 812 813 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 814 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 815 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 816 817 if (kspUpper == kspA) { 818 ierr = MatMult(jac->B,ilinkD->y,ilinkA->y);CHKERRQ(ierr); 819 ierr = VecAXPY(ilinkA->x,-1.0,ilinkA->y);CHKERRQ(ierr); 820 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 821 } else { 822 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 823 ierr = MatMult(jac->B,ilinkD->y,ilinkA->x);CHKERRQ(ierr); 824 ierr = KSPSolve(kspUpper,ilinkA->x,ilinkA->z);CHKERRQ(ierr); 825 ierr = VecAXPY(ilinkA->y,-1.0,ilinkA->z);CHKERRQ(ierr); 826 } 827 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 828 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 829 } 830 PetscFunctionReturn(0); 831 } 832 833 #undef __FUNCT__ 834 #define __FUNCT__ "PCApply_FieldSplit" 835 static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y) 836 { 837 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 838 PetscErrorCode ierr; 839 PC_FieldSplitLink ilink = jac->head; 840 PetscInt cnt,bs; 841 842 PetscFunctionBegin; 843 if (jac->type == PC_COMPOSITE_ADDITIVE) { 844 if (jac->defaultsplit) { 845 ierr = VecGetBlockSize(x,&bs);CHKERRQ(ierr); 846 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); 847 ierr = VecGetBlockSize(y,&bs);CHKERRQ(ierr); 848 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); 849 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 850 while (ilink) { 851 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 852 ilink = ilink->next; 853 } 854 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 855 } else { 856 ierr = VecSet(y,0.0);CHKERRQ(ierr); 857 while (ilink) { 858 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 859 ilink = ilink->next; 860 } 861 } 862 } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 863 if (!jac->w1) { 864 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 865 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 866 } 867 ierr = VecSet(y,0.0);CHKERRQ(ierr); 868 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 869 cnt = 1; 870 while (ilink->next) { 871 ilink = ilink->next; 872 /* compute the residual only over the part of the vector needed */ 873 ierr = MatMult(jac->Afield[cnt++],y,ilink->x);CHKERRQ(ierr); 874 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 875 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 876 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 877 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 878 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 879 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 880 } 881 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 882 cnt -= 2; 883 while (ilink->previous) { 884 ilink = ilink->previous; 885 /* compute the residual only over the part of the vector needed */ 886 ierr = MatMult(jac->Afield[cnt--],y,ilink->x);CHKERRQ(ierr); 887 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 888 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 889 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 890 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 891 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 892 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 893 } 894 } 895 } else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Unsupported or unknown composition",(int) jac->type); 896 PetscFunctionReturn(0); 897 } 898 899 #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \ 900 (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 901 VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 902 KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) || \ 903 VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \ 904 VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE)) 905 906 #undef __FUNCT__ 907 #define __FUNCT__ "PCApplyTranspose_FieldSplit" 908 static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y) 909 { 910 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 911 PetscErrorCode ierr; 912 PC_FieldSplitLink ilink = jac->head; 913 PetscInt bs; 914 915 PetscFunctionBegin; 916 if (jac->type == PC_COMPOSITE_ADDITIVE) { 917 if (jac->defaultsplit) { 918 ierr = VecGetBlockSize(x,&bs);CHKERRQ(ierr); 919 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); 920 ierr = VecGetBlockSize(y,&bs);CHKERRQ(ierr); 921 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); 922 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 923 while (ilink) { 924 ierr = KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 925 ilink = ilink->next; 926 } 927 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 928 } else { 929 ierr = VecSet(y,0.0);CHKERRQ(ierr); 930 while (ilink) { 931 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 932 ilink = ilink->next; 933 } 934 } 935 } else { 936 if (!jac->w1) { 937 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 938 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 939 } 940 ierr = VecSet(y,0.0);CHKERRQ(ierr); 941 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 942 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 943 while (ilink->next) { 944 ilink = ilink->next; 945 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 946 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 947 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 948 } 949 while (ilink->previous) { 950 ilink = ilink->previous; 951 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 952 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 953 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 954 } 955 } else { 956 while (ilink->next) { /* get to last entry in linked list */ 957 ilink = ilink->next; 958 } 959 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 960 while (ilink->previous) { 961 ilink = ilink->previous; 962 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 963 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 964 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 965 } 966 } 967 } 968 PetscFunctionReturn(0); 969 } 970 971 #undef __FUNCT__ 972 #define __FUNCT__ "PCReset_FieldSplit" 973 static PetscErrorCode PCReset_FieldSplit(PC pc) 974 { 975 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 976 PetscErrorCode ierr; 977 PC_FieldSplitLink ilink = jac->head,next; 978 979 PetscFunctionBegin; 980 while (ilink) { 981 ierr = KSPReset(ilink->ksp);CHKERRQ(ierr); 982 ierr = VecDestroy(&ilink->x);CHKERRQ(ierr); 983 ierr = VecDestroy(&ilink->y);CHKERRQ(ierr); 984 ierr = VecDestroy(&ilink->z);CHKERRQ(ierr); 985 ierr = VecScatterDestroy(&ilink->sctx);CHKERRQ(ierr); 986 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 987 ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr); 988 next = ilink->next; 989 ilink = next; 990 } 991 ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr); 992 if (jac->mat && jac->mat != jac->pmat) { 993 ierr = MatDestroyMatrices(jac->nsplits,&jac->mat);CHKERRQ(ierr); 994 } else if (jac->mat) { 995 jac->mat = NULL; 996 } 997 if (jac->pmat) {ierr = MatDestroyMatrices(jac->nsplits,&jac->pmat);CHKERRQ(ierr);} 998 if (jac->Afield) {ierr = MatDestroyMatrices(jac->nsplits,&jac->Afield);CHKERRQ(ierr);} 999 ierr = VecDestroy(&jac->w1);CHKERRQ(ierr); 1000 ierr = VecDestroy(&jac->w2);CHKERRQ(ierr); 1001 ierr = MatDestroy(&jac->schur);CHKERRQ(ierr); 1002 ierr = MatDestroy(&jac->schur_user);CHKERRQ(ierr); 1003 ierr = KSPDestroy(&jac->kspschur);CHKERRQ(ierr); 1004 ierr = KSPDestroy(&jac->kspupper);CHKERRQ(ierr); 1005 ierr = MatDestroy(&jac->B);CHKERRQ(ierr); 1006 ierr = MatDestroy(&jac->C);CHKERRQ(ierr); 1007 jac->reset = PETSC_TRUE; 1008 PetscFunctionReturn(0); 1009 } 1010 1011 #undef __FUNCT__ 1012 #define __FUNCT__ "PCDestroy_FieldSplit" 1013 static PetscErrorCode PCDestroy_FieldSplit(PC pc) 1014 { 1015 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1016 PetscErrorCode ierr; 1017 PC_FieldSplitLink ilink = jac->head,next; 1018 1019 PetscFunctionBegin; 1020 ierr = PCReset_FieldSplit(pc);CHKERRQ(ierr); 1021 while (ilink) { 1022 ierr = KSPDestroy(&ilink->ksp);CHKERRQ(ierr); 1023 next = ilink->next; 1024 ierr = PetscFree(ilink->splitname);CHKERRQ(ierr); 1025 ierr = PetscFree(ilink->fields);CHKERRQ(ierr); 1026 ierr = PetscFree(ilink->fields_col);CHKERRQ(ierr); 1027 ierr = PetscFree(ilink);CHKERRQ(ierr); 1028 ilink = next; 1029 } 1030 ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr); 1031 ierr = PetscFree(pc->data);CHKERRQ(ierr); 1032 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",NULL);CHKERRQ(ierr); 1033 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",NULL);CHKERRQ(ierr); 1034 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",NULL);CHKERRQ(ierr); 1035 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",NULL);CHKERRQ(ierr); 1036 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",NULL);CHKERRQ(ierr); 1037 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSchurPrecondition_C",NULL);CHKERRQ(ierr); 1038 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",NULL);CHKERRQ(ierr); 1039 PetscFunctionReturn(0); 1040 } 1041 1042 #undef __FUNCT__ 1043 #define __FUNCT__ "PCSetFromOptions_FieldSplit" 1044 static PetscErrorCode PCSetFromOptions_FieldSplit(PC pc) 1045 { 1046 PetscErrorCode ierr; 1047 PetscInt bs; 1048 PetscBool flg,stokes = PETSC_FALSE; 1049 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1050 PCCompositeType ctype; 1051 1052 PetscFunctionBegin; 1053 ierr = PetscOptionsHead("FieldSplit options");CHKERRQ(ierr); 1054 ierr = PetscOptionsBool("-pc_fieldsplit_dm_splits","Whether to use DMCreateFieldDecomposition() for splits","PCFieldSplitSetDMSplits",jac->dm_splits,&jac->dm_splits,NULL);CHKERRQ(ierr); 1055 ierr = PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);CHKERRQ(ierr); 1056 if (flg) { 1057 ierr = PCFieldSplitSetBlockSize(pc,bs);CHKERRQ(ierr); 1058 } 1059 1060 ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,NULL);CHKERRQ(ierr); 1061 if (stokes) { 1062 ierr = PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);CHKERRQ(ierr); 1063 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_SELF; 1064 } 1065 1066 ierr = PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg);CHKERRQ(ierr); 1067 if (flg) { 1068 ierr = PCFieldSplitSetType(pc,ctype);CHKERRQ(ierr); 1069 } 1070 /* Only setup fields once */ 1071 if ((jac->bs > 0) && (jac->nsplits == 0)) { 1072 /* only allow user to set fields from command line if bs is already known. 1073 otherwise user can set them in PCFieldSplitSetDefaults() */ 1074 ierr = PCFieldSplitSetRuntimeSplits_Private(pc);CHKERRQ(ierr); 1075 if (jac->splitdefined) {ierr = PetscInfo(pc,"Splits defined using the options database\n");CHKERRQ(ierr);} 1076 } 1077 if (jac->type == PC_COMPOSITE_SCHUR) { 1078 ierr = PetscOptionsGetEnum(((PetscObject)pc)->prefix,"-pc_fieldsplit_schur_factorization_type",PCFieldSplitSchurFactTypes,(PetscEnum*)&jac->schurfactorization,&flg);CHKERRQ(ierr); 1079 if (flg) {ierr = PetscInfo(pc,"Deprecated use of -pc_fieldsplit_schur_factorization_type\n");CHKERRQ(ierr);} 1080 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); 1081 ierr = PetscOptionsEnum("-pc_fieldsplit_schur_precondition","How to build preconditioner for Schur complement","PCFieldSplitSchurPrecondition",PCFieldSplitSchurPreTypes,(PetscEnum)jac->schurpre,(PetscEnum*)&jac->schurpre,NULL);CHKERRQ(ierr); 1082 } 1083 ierr = PetscOptionsTail();CHKERRQ(ierr); 1084 PetscFunctionReturn(0); 1085 } 1086 1087 /*------------------------------------------------------------------------------------*/ 1088 1089 #undef __FUNCT__ 1090 #define __FUNCT__ "PCFieldSplitSetFields_FieldSplit" 1091 static PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 1092 { 1093 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1094 PetscErrorCode ierr; 1095 PC_FieldSplitLink ilink,next = jac->head; 1096 char prefix[128]; 1097 PetscInt i; 1098 1099 PetscFunctionBegin; 1100 if (jac->splitdefined) { 1101 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 1102 PetscFunctionReturn(0); 1103 } 1104 for (i=0; i<n; i++) { 1105 if (fields[i] >= jac->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs); 1106 if (fields[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]); 1107 } 1108 ierr = PetscNew(struct _PC_FieldSplitLink,&ilink);CHKERRQ(ierr); 1109 if (splitname) { 1110 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 1111 } else { 1112 ierr = PetscMalloc(3*sizeof(char),&ilink->splitname);CHKERRQ(ierr); 1113 ierr = PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);CHKERRQ(ierr); 1114 } 1115 ierr = PetscMalloc(n*sizeof(PetscInt),&ilink->fields);CHKERRQ(ierr); 1116 ierr = PetscMemcpy(ilink->fields,fields,n*sizeof(PetscInt));CHKERRQ(ierr); 1117 ierr = PetscMalloc(n*sizeof(PetscInt),&ilink->fields_col);CHKERRQ(ierr); 1118 ierr = PetscMemcpy(ilink->fields_col,fields_col,n*sizeof(PetscInt));CHKERRQ(ierr); 1119 1120 ilink->nfields = n; 1121 ilink->next = NULL; 1122 ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);CHKERRQ(ierr); 1123 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 1124 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 1125 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 1126 1127 ierr = PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);CHKERRQ(ierr); 1128 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 1129 1130 if (!next) { 1131 jac->head = ilink; 1132 ilink->previous = NULL; 1133 } else { 1134 while (next->next) { 1135 next = next->next; 1136 } 1137 next->next = ilink; 1138 ilink->previous = next; 1139 } 1140 jac->nsplits++; 1141 PetscFunctionReturn(0); 1142 } 1143 1144 #undef __FUNCT__ 1145 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit_Schur" 1146 static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp) 1147 { 1148 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1149 PetscErrorCode ierr; 1150 1151 PetscFunctionBegin; 1152 ierr = PetscMalloc(jac->nsplits*sizeof(KSP),subksp);CHKERRQ(ierr); 1153 ierr = MatSchurComplementGetKSP(jac->schur,*subksp);CHKERRQ(ierr); 1154 1155 (*subksp)[1] = jac->kspschur; 1156 if (n) *n = jac->nsplits; 1157 PetscFunctionReturn(0); 1158 } 1159 1160 #undef __FUNCT__ 1161 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit" 1162 static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp) 1163 { 1164 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1165 PetscErrorCode ierr; 1166 PetscInt cnt = 0; 1167 PC_FieldSplitLink ilink = jac->head; 1168 1169 PetscFunctionBegin; 1170 ierr = PetscMalloc(jac->nsplits*sizeof(KSP),subksp);CHKERRQ(ierr); 1171 while (ilink) { 1172 (*subksp)[cnt++] = ilink->ksp; 1173 ilink = ilink->next; 1174 } 1175 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); 1176 if (n) *n = jac->nsplits; 1177 PetscFunctionReturn(0); 1178 } 1179 1180 #undef __FUNCT__ 1181 #define __FUNCT__ "PCFieldSplitSetIS_FieldSplit" 1182 static PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is) 1183 { 1184 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1185 PetscErrorCode ierr; 1186 PC_FieldSplitLink ilink, next = jac->head; 1187 char prefix[128]; 1188 1189 PetscFunctionBegin; 1190 if (jac->splitdefined) { 1191 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 1192 PetscFunctionReturn(0); 1193 } 1194 ierr = PetscNew(struct _PC_FieldSplitLink,&ilink);CHKERRQ(ierr); 1195 if (splitname) { 1196 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 1197 } else { 1198 ierr = PetscMalloc(8*sizeof(char),&ilink->splitname);CHKERRQ(ierr); 1199 ierr = PetscSNPrintf(ilink->splitname,7,"%D",jac->nsplits);CHKERRQ(ierr); 1200 } 1201 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 1202 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 1203 ilink->is = is; 1204 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 1205 ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr); 1206 ilink->is_col = is; 1207 ilink->next = NULL; 1208 ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);CHKERRQ(ierr); 1209 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 1210 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 1211 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 1212 1213 ierr = PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);CHKERRQ(ierr); 1214 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 1215 1216 if (!next) { 1217 jac->head = ilink; 1218 ilink->previous = NULL; 1219 } else { 1220 while (next->next) { 1221 next = next->next; 1222 } 1223 next->next = ilink; 1224 ilink->previous = next; 1225 } 1226 jac->nsplits++; 1227 PetscFunctionReturn(0); 1228 } 1229 1230 #undef __FUNCT__ 1231 #define __FUNCT__ "PCFieldSplitSetFields" 1232 /*@ 1233 PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner 1234 1235 Logically Collective on PC 1236 1237 Input Parameters: 1238 + pc - the preconditioner context 1239 . splitname - name of this split, if NULL the number of the split is used 1240 . n - the number of fields in this split 1241 - fields - the fields in this split 1242 1243 Level: intermediate 1244 1245 Notes: Use PCFieldSplitSetIS() to set a completely general set of indices as a field. 1246 1247 The PCFieldSplitSetFields() is for defining fields as strided blocks. For example, if the block 1248 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 1249 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x.... 1250 where the numbered entries indicate what is in the field. 1251 1252 This function is called once per split (it creates a new split each time). Solve options 1253 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1254 1255 Developer Note: This routine does not actually create the IS representing the split, that is delayed 1256 until PCSetUp_FieldSplit(), because information about the vector/matrix layouts may not be 1257 available when this routine is called. 1258 1259 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize(), PCFieldSplitSetIS() 1260 1261 @*/ 1262 PetscErrorCode PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 1263 { 1264 PetscErrorCode ierr; 1265 1266 PetscFunctionBegin; 1267 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1268 PetscValidCharPointer(splitname,2); 1269 if (n < 1) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Provided number of fields %D in split \"%s\" not positive",n,splitname); 1270 PetscValidIntPointer(fields,3); 1271 ierr = PetscTryMethod(pc,"PCFieldSplitSetFields_C",(PC,const char[],PetscInt,const PetscInt*,const PetscInt*),(pc,splitname,n,fields,fields_col));CHKERRQ(ierr); 1272 PetscFunctionReturn(0); 1273 } 1274 1275 #undef __FUNCT__ 1276 #define __FUNCT__ "PCFieldSplitSetIS" 1277 /*@ 1278 PCFieldSplitSetIS - Sets the exact elements for field 1279 1280 Logically Collective on PC 1281 1282 Input Parameters: 1283 + pc - the preconditioner context 1284 . splitname - name of this split, if NULL the number of the split is used 1285 - is - the index set that defines the vector elements in this field 1286 1287 1288 Notes: 1289 Use PCFieldSplitSetFields(), for fields defined by strided types. 1290 1291 This function is called once per split (it creates a new split each time). Solve options 1292 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1293 1294 Level: intermediate 1295 1296 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize() 1297 1298 @*/ 1299 PetscErrorCode PCFieldSplitSetIS(PC pc,const char splitname[],IS is) 1300 { 1301 PetscErrorCode ierr; 1302 1303 PetscFunctionBegin; 1304 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1305 if (splitname) PetscValidCharPointer(splitname,2); 1306 PetscValidHeaderSpecific(is,IS_CLASSID,3); 1307 ierr = PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));CHKERRQ(ierr); 1308 PetscFunctionReturn(0); 1309 } 1310 1311 #undef __FUNCT__ 1312 #define __FUNCT__ "PCFieldSplitGetIS" 1313 /*@ 1314 PCFieldSplitGetIS - Retrieves the elements for a field as an IS 1315 1316 Logically Collective on PC 1317 1318 Input Parameters: 1319 + pc - the preconditioner context 1320 - splitname - name of this split 1321 1322 Output Parameter: 1323 - is - the index set that defines the vector elements in this field, or NULL if the field is not found 1324 1325 Level: intermediate 1326 1327 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetIS() 1328 1329 @*/ 1330 PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is) 1331 { 1332 PetscErrorCode ierr; 1333 1334 PetscFunctionBegin; 1335 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1336 PetscValidCharPointer(splitname,2); 1337 PetscValidPointer(is,3); 1338 { 1339 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 1340 PC_FieldSplitLink ilink = jac->head; 1341 PetscBool found; 1342 1343 *is = NULL; 1344 while (ilink) { 1345 ierr = PetscStrcmp(ilink->splitname, splitname, &found);CHKERRQ(ierr); 1346 if (found) { 1347 *is = ilink->is; 1348 break; 1349 } 1350 ilink = ilink->next; 1351 } 1352 } 1353 PetscFunctionReturn(0); 1354 } 1355 1356 #undef __FUNCT__ 1357 #define __FUNCT__ "PCFieldSplitSetBlockSize" 1358 /*@ 1359 PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the 1360 fieldsplit preconditioner. If not set the matrix block size is used. 1361 1362 Logically Collective on PC 1363 1364 Input Parameters: 1365 + pc - the preconditioner context 1366 - bs - the block size 1367 1368 Level: intermediate 1369 1370 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields() 1371 1372 @*/ 1373 PetscErrorCode PCFieldSplitSetBlockSize(PC pc,PetscInt bs) 1374 { 1375 PetscErrorCode ierr; 1376 1377 PetscFunctionBegin; 1378 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1379 PetscValidLogicalCollectiveInt(pc,bs,2); 1380 ierr = PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));CHKERRQ(ierr); 1381 PetscFunctionReturn(0); 1382 } 1383 1384 #undef __FUNCT__ 1385 #define __FUNCT__ "PCFieldSplitGetSubKSP" 1386 /*@C 1387 PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits 1388 1389 Collective on KSP 1390 1391 Input Parameter: 1392 . pc - the preconditioner context 1393 1394 Output Parameters: 1395 + n - the number of splits 1396 - pc - the array of KSP contexts 1397 1398 Note: 1399 After PCFieldSplitGetSubKSP() the array of KSPs IS to be freed by the user 1400 (not the KSP just the array that contains them). 1401 1402 You must call KSPSetUp() before calling PCFieldSplitGetSubKSP(). 1403 1404 Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs. 1405 You can call PCFieldSplitGetSubKSP(pc,n,NULL_OBJECT,ierr) to determine how large the 1406 KSP array must be. 1407 1408 1409 Level: advanced 1410 1411 .seealso: PCFIELDSPLIT 1412 @*/ 1413 PetscErrorCode PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[]) 1414 { 1415 PetscErrorCode ierr; 1416 1417 PetscFunctionBegin; 1418 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1419 if (n) PetscValidIntPointer(n,2); 1420 ierr = PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));CHKERRQ(ierr); 1421 PetscFunctionReturn(0); 1422 } 1423 1424 #undef __FUNCT__ 1425 #define __FUNCT__ "PCFieldSplitSchurPrecondition" 1426 /*@ 1427 PCFieldSplitSchurPrecondition - Indicates if the Schur complement is preconditioned by a preconditioner constructed by the 1428 A11 matrix. Otherwise no preconditioner is used. 1429 1430 Collective on PC 1431 1432 Input Parameters: 1433 + pc - the preconditioner context 1434 . ptype - which matrix to use for preconditioning the Schur complement, PC_FIELDSPLIT_SCHUR_PRE_A11 (diag) is default 1435 - userpre - matrix to use for preconditioning, or NULL 1436 1437 Options Database: 1438 . -pc_fieldsplit_schur_precondition <self,user,a11> default is a11 1439 1440 Notes: 1441 $ If ptype is 1442 $ user then the preconditioner for the Schur complement is generated by the provided matrix (pre argument 1443 $ to this function). 1444 $ a11 then the preconditioner for the Schur complement is generated by the block diagonal part of the original 1445 $ matrix associated with the Schur complement (i.e. A11) 1446 $ self the preconditioner for the Schur complement is generated from the Schur complement matrix itself: 1447 $ The only preconditioner that currently works directly with the Schur complement matrix object is the PCLSC 1448 $ preconditioner 1449 1450 When solving a saddle point problem, where the A11 block is identically zero, using a11 as the ptype only makes sense 1451 with the additional option -fieldsplit_1_pc_type none. Usually for saddle point problems one would use a ptype of self and 1452 -fieldsplit_1_pc_type lsc which uses the least squares commutator compute a preconditioner for the Schur complement. 1453 1454 Developer Notes: This is a terrible name, gives no good indication of what the function does and should also have Set in 1455 the name since it sets a proceedure to use. 1456 1457 Level: intermediate 1458 1459 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCLSC 1460 1461 @*/ 1462 PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1463 { 1464 PetscErrorCode ierr; 1465 1466 PetscFunctionBegin; 1467 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1468 ierr = PetscTryMethod(pc,"PCFieldSplitSchurPrecondition_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));CHKERRQ(ierr); 1469 PetscFunctionReturn(0); 1470 } 1471 1472 #undef __FUNCT__ 1473 #define __FUNCT__ "PCFieldSplitSchurPrecondition_FieldSplit" 1474 static PetscErrorCode PCFieldSplitSchurPrecondition_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1475 { 1476 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1477 PetscErrorCode ierr; 1478 1479 PetscFunctionBegin; 1480 jac->schurpre = ptype; 1481 if (pre) { 1482 ierr = MatDestroy(&jac->schur_user);CHKERRQ(ierr); 1483 jac->schur_user = pre; 1484 ierr = PetscObjectReference((PetscObject)jac->schur_user);CHKERRQ(ierr); 1485 } 1486 PetscFunctionReturn(0); 1487 } 1488 1489 #undef __FUNCT__ 1490 #define __FUNCT__ "PCFieldSplitSetSchurFactType" 1491 /*@ 1492 PCFieldSplitSetSchurFactType - sets which blocks of the approximate block factorization to retain 1493 1494 Collective on PC 1495 1496 Input Parameters: 1497 + pc - the preconditioner context 1498 - ftype - which blocks of factorization to retain, PC_FIELDSPLIT_SCHUR_FACT_FULL is default 1499 1500 Options Database: 1501 . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> default is full 1502 1503 1504 Level: intermediate 1505 1506 Notes: 1507 The FULL factorization is 1508 1509 $ (A B) = (1 0) (A 0) (1 Ainv*B) 1510 $ (C D) (C*Ainv 1) (0 S) (0 1 ) 1511 1512 where S = D - 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, 1513 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). 1514 1515 If applied exactly, FULL factorization is a direct solver. The preconditioned operator with LOWER or UPPER has all eigenvalues equal to 1 and minimal polynomial 1516 of degree 2, so KSPGMRES converges in 2 iterations. If the iteration count is very low, consider using KSPFGMRES or KSPGCR which can use one less preconditioner 1517 application in this case. Note that the preconditioned operator may be highly non-normal, so such fast convergence may not be observed in practice. With DIAG, 1518 the preconditioned operator has three distinct nonzero eigenvalues and minimal polynomial of degree at most 4, so KSPGMRES converges in at most 4 iterations. 1519 1520 For symmetric problems in which A is positive definite and S is negative definite, DIAG can be used with KSPMINRES. Note that a flexible method like KSPFGMRES 1521 or KSPGCR must be used if the fieldsplit preconditioner is nonlinear (e.g. a few iterations of a Krylov method is used inside a split). 1522 1523 References: 1524 Murphy, Golub, and Wathen, A note on preconditioning indefinite linear systems, SIAM J. Sci. Comput., 21 (2000) pp. 1969-1972. 1525 1526 Ipsen, A note on preconditioning nonsymmetric matrices, SIAM J. Sci. Comput., 23 (2001), pp. 1050-1051. 1527 1528 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType 1529 @*/ 1530 PetscErrorCode PCFieldSplitSetSchurFactType(PC pc,PCFieldSplitSchurFactType ftype) 1531 { 1532 PetscErrorCode ierr; 1533 1534 PetscFunctionBegin; 1535 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1536 ierr = PetscTryMethod(pc,"PCFieldSplitSetSchurFactType_C",(PC,PCFieldSplitSchurFactType),(pc,ftype));CHKERRQ(ierr); 1537 PetscFunctionReturn(0); 1538 } 1539 1540 #undef __FUNCT__ 1541 #define __FUNCT__ "PCFieldSplitSetSchurFactType_FieldSplit" 1542 static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc,PCFieldSplitSchurFactType ftype) 1543 { 1544 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1545 1546 PetscFunctionBegin; 1547 jac->schurfactorization = ftype; 1548 PetscFunctionReturn(0); 1549 } 1550 1551 #undef __FUNCT__ 1552 #define __FUNCT__ "PCFieldSplitGetSchurBlocks" 1553 /*@C 1554 PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement 1555 1556 Collective on KSP 1557 1558 Input Parameter: 1559 . pc - the preconditioner context 1560 1561 Output Parameters: 1562 + A00 - the (0,0) block 1563 . A01 - the (0,1) block 1564 . A10 - the (1,0) block 1565 - A11 - the (1,1) block 1566 1567 Level: advanced 1568 1569 .seealso: PCFIELDSPLIT 1570 @*/ 1571 PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11) 1572 { 1573 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 1574 1575 PetscFunctionBegin; 1576 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1577 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach."); 1578 if (A00) *A00 = jac->pmat[0]; 1579 if (A01) *A01 = jac->B; 1580 if (A10) *A10 = jac->C; 1581 if (A11) *A11 = jac->pmat[1]; 1582 PetscFunctionReturn(0); 1583 } 1584 1585 #undef __FUNCT__ 1586 #define __FUNCT__ "PCFieldSplitSetType_FieldSplit" 1587 static PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type) 1588 { 1589 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1590 PetscErrorCode ierr; 1591 1592 PetscFunctionBegin; 1593 jac->type = type; 1594 if (type == PC_COMPOSITE_SCHUR) { 1595 pc->ops->apply = PCApply_FieldSplit_Schur; 1596 pc->ops->view = PCView_FieldSplit_Schur; 1597 1598 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit_Schur);CHKERRQ(ierr); 1599 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSchurPrecondition_C",PCFieldSplitSchurPrecondition_FieldSplit);CHKERRQ(ierr); 1600 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",PCFieldSplitSetSchurFactType_FieldSplit);CHKERRQ(ierr); 1601 1602 } else { 1603 pc->ops->apply = PCApply_FieldSplit; 1604 pc->ops->view = PCView_FieldSplit; 1605 1606 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 1607 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSchurPrecondition_C",0);CHKERRQ(ierr); 1608 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",0);CHKERRQ(ierr); 1609 } 1610 PetscFunctionReturn(0); 1611 } 1612 1613 #undef __FUNCT__ 1614 #define __FUNCT__ "PCFieldSplitSetBlockSize_FieldSplit" 1615 static PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs) 1616 { 1617 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1618 1619 PetscFunctionBegin; 1620 if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs); 1621 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); 1622 jac->bs = bs; 1623 PetscFunctionReturn(0); 1624 } 1625 1626 #undef __FUNCT__ 1627 #define __FUNCT__ "PCFieldSplitSetType" 1628 /*@ 1629 PCFieldSplitSetType - Sets the type of fieldsplit preconditioner. 1630 1631 Collective on PC 1632 1633 Input Parameter: 1634 . pc - the preconditioner context 1635 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 1636 1637 Options Database Key: 1638 . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type 1639 1640 Level: Intermediate 1641 1642 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 1643 1644 .seealso: PCCompositeSetType() 1645 1646 @*/ 1647 PetscErrorCode PCFieldSplitSetType(PC pc,PCCompositeType type) 1648 { 1649 PetscErrorCode ierr; 1650 1651 PetscFunctionBegin; 1652 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1653 ierr = PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));CHKERRQ(ierr); 1654 PetscFunctionReturn(0); 1655 } 1656 1657 #undef __FUNCT__ 1658 #define __FUNCT__ "PCFieldSplitGetType" 1659 /*@ 1660 PCFieldSplitGetType - Gets the type of fieldsplit preconditioner. 1661 1662 Not collective 1663 1664 Input Parameter: 1665 . pc - the preconditioner context 1666 1667 Output Parameter: 1668 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 1669 1670 Level: Intermediate 1671 1672 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 1673 .seealso: PCCompositeSetType() 1674 @*/ 1675 PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type) 1676 { 1677 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 1678 1679 PetscFunctionBegin; 1680 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1681 PetscValidIntPointer(type,2); 1682 *type = jac->type; 1683 PetscFunctionReturn(0); 1684 } 1685 1686 #undef __FUNCT__ 1687 #define __FUNCT__ "PCFieldSplitSetDMSplits" 1688 /*@ 1689 PCFieldSplitSetDMSplits - Flags whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible. 1690 1691 Logically Collective 1692 1693 Input Parameters: 1694 + pc - the preconditioner context 1695 - flg - boolean indicating whether to use field splits defined by the DM 1696 1697 Options Database Key: 1698 . -pc_fieldsplit_dm_splits 1699 1700 Level: Intermediate 1701 1702 .keywords: PC, DM, composite preconditioner, additive, multiplicative 1703 1704 .seealso: PCFieldSplitGetDMSplits() 1705 1706 @*/ 1707 PetscErrorCode PCFieldSplitSetDMSplits(PC pc,PetscBool flg) 1708 { 1709 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1710 PetscBool isfs; 1711 PetscErrorCode ierr; 1712 1713 PetscFunctionBegin; 1714 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1715 PetscValidLogicalCollectiveBool(pc,flg,2); 1716 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1717 if (isfs) { 1718 jac->dm_splits = flg; 1719 } 1720 PetscFunctionReturn(0); 1721 } 1722 1723 1724 #undef __FUNCT__ 1725 #define __FUNCT__ "PCFieldSplitGetDMSplits" 1726 /*@ 1727 PCFieldSplitGetDMSplits - Returns flag indicating whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible. 1728 1729 Logically Collective 1730 1731 Input Parameter: 1732 . pc - the preconditioner context 1733 1734 Output Parameter: 1735 . flg - boolean indicating whether to use field splits defined by the DM 1736 1737 Level: Intermediate 1738 1739 .keywords: PC, DM, composite preconditioner, additive, multiplicative 1740 1741 .seealso: PCFieldSplitSetDMSplits() 1742 1743 @*/ 1744 PetscErrorCode PCFieldSplitGetDMSplits(PC pc,PetscBool* flg) 1745 { 1746 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1747 PetscBool isfs; 1748 PetscErrorCode ierr; 1749 1750 PetscFunctionBegin; 1751 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1752 PetscValidPointer(flg,2); 1753 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1754 if (isfs) { 1755 if(flg) *flg = jac->dm_splits; 1756 } 1757 PetscFunctionReturn(0); 1758 } 1759 1760 /* -------------------------------------------------------------------------------------*/ 1761 /*MC 1762 PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual 1763 fields or groups of fields. See the users manual section "Solving Block Matrices" for more details. 1764 1765 To set options on the solvers for each block append -fieldsplit_ to all the PC 1766 options database keys. For example, -fieldsplit_pc_type ilu -fieldsplit_pc_factor_levels 1 1767 1768 To set the options on the solvers separate for each block call PCFieldSplitGetSubKSP() 1769 and set the options directly on the resulting KSP object 1770 1771 Level: intermediate 1772 1773 Options Database Keys: 1774 + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split 1775 . -pc_fieldsplit_default - automatically add any fields to additional splits that have not 1776 been supplied explicitly by -pc_fieldsplit_%d_fields 1777 . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields) 1778 . -pc_fieldsplit_type <additive,multiplicative,symmetric_multiplicative,schur> - type of relaxation or factorization splitting 1779 . -pc_fieldsplit_schur_precondition <self,user,a11> - default is a11 1780 . -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero or negative diagonal and uses Schur complement with no preconditioner as the solver 1781 1782 - Options prefix for inner solvers when using Schur complement preconditioner are -fieldsplit_0_ and -fieldsplit_1_ 1783 for all other solvers they are -fieldsplit_%d_ for the dth field, use -fieldsplit_ for all fields 1784 1785 Notes: 1786 Use PCFieldSplitSetFields() to set fields defined by "strided" entries and PCFieldSplitSetIS() 1787 to define a field by an arbitrary collection of entries. 1788 1789 If no fields are set the default is used. The fields are defined by entries strided by bs, 1790 beginning at 0 then 1, etc to bs-1. The block size can be set with PCFieldSplitSetBlockSize(), 1791 if this is not called the block size defaults to the blocksize of the second matrix passed 1792 to KSPSetOperators()/PCSetOperators(). 1793 1794 $ For the Schur complement preconditioner if J = ( A00 A01 ) 1795 $ ( A10 A11 ) 1796 $ the preconditioner using full factorization is 1797 $ ( I -A10 ksp(A00) ) ( inv(A00) 0 ) ( I 0 ) 1798 $ ( 0 I ) ( 0 ksp(S) ) ( -A10 ksp(A00) I ) 1799 where the action of inv(A00) is applied using the KSP solver with prefix -fieldsplit_0_. The action of 1800 ksp(S) is computed using the KSP solver with prefix -fieldsplit_splitname_ (where splitname was given 1801 in providing the SECOND split or 1 if not give). For PCFieldSplitGetKSP() when field number is 0, 1802 it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default 1803 A11 is used to construct a preconditioner for S, use PCFieldSplitSchurPrecondition() to turn on or off this 1804 option. You can use the preconditioner PCLSC to precondition the Schur complement with -fieldsplit_1_pc_type lsc. The 1805 factorization type is set using -pc_fieldsplit_schur_fact_type <diag, lower, upper, full>. The full is shown above, 1806 diag gives 1807 $ ( inv(A00) 0 ) 1808 $ ( 0 -ksp(S) ) 1809 note that slightly counter intuitively there is a negative in front of the ksp(S) so that the preconditioner is positive definite. The lower factorization is the inverse of 1810 $ ( A00 0 ) 1811 $ ( A10 S ) 1812 where the inverses of A00 and S are applied using KSPs. The upper factorization is the inverse of 1813 $ ( A00 A01 ) 1814 $ ( 0 S ) 1815 where again the inverses of A00 and S are applied using KSPs. 1816 1817 If only one set of indices (one IS) is provided with PCFieldSplitSetIS() then the complement of that IS 1818 is used automatically for a second block. 1819 1820 The fieldsplit preconditioner cannot currently be used with the BAIJ or SBAIJ data formats if the blocksize is larger than 1. 1821 Generally it should be used with the AIJ format. 1822 1823 The forms of these preconditioners are closely related if not identical to forms derived as "Distributive Iterations", see, 1824 for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling. Note that one can also use PCFIELDSPLIT 1825 inside a smoother resulting in "Distributive Smoothers". 1826 1827 Concepts: physics based preconditioners, block preconditioners 1828 1829 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners, PCLSC, 1830 PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSchurPrecondition() 1831 M*/ 1832 1833 #undef __FUNCT__ 1834 #define __FUNCT__ "PCCreate_FieldSplit" 1835 PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc) 1836 { 1837 PetscErrorCode ierr; 1838 PC_FieldSplit *jac; 1839 1840 PetscFunctionBegin; 1841 ierr = PetscNewLog(pc,PC_FieldSplit,&jac);CHKERRQ(ierr); 1842 1843 jac->bs = -1; 1844 jac->nsplits = 0; 1845 jac->type = PC_COMPOSITE_MULTIPLICATIVE; 1846 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */ 1847 jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL; 1848 jac->dm_splits = PETSC_TRUE; 1849 1850 pc->data = (void*)jac; 1851 1852 pc->ops->apply = PCApply_FieldSplit; 1853 pc->ops->applytranspose = PCApplyTranspose_FieldSplit; 1854 pc->ops->setup = PCSetUp_FieldSplit; 1855 pc->ops->reset = PCReset_FieldSplit; 1856 pc->ops->destroy = PCDestroy_FieldSplit; 1857 pc->ops->setfromoptions = PCSetFromOptions_FieldSplit; 1858 pc->ops->view = PCView_FieldSplit; 1859 pc->ops->applyrichardson = 0; 1860 1861 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 1862 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",PCFieldSplitSetFields_FieldSplit);CHKERRQ(ierr); 1863 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",PCFieldSplitSetIS_FieldSplit);CHKERRQ(ierr); 1864 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",PCFieldSplitSetType_FieldSplit);CHKERRQ(ierr); 1865 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",PCFieldSplitSetBlockSize_FieldSplit);CHKERRQ(ierr); 1866 PetscFunctionReturn(0); 1867 } 1868 1869 1870 1871