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 /* propogate DM */ 711 { 712 DM sdm; 713 ierr = KSPGetDM(jac->head->next->ksp, &sdm);CHKERRQ(ierr); 714 if (sdm) { 715 ierr = KSPSetDM(jac->kspschur, sdm);CHKERRQ(ierr); 716 ierr = KSPSetDMActive(jac->kspschur, PETSC_FALSE);CHKERRQ(ierr); 717 } 718 } 719 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 720 /* need to call this every time, since the jac->kspschur is freshly created, otherwise its options never get set */ 721 ierr = KSPSetFromOptions(jac->kspschur);CHKERRQ(ierr); 722 } 723 724 /* HACK: special support to forward L and Lp matrices that might be used by PCLSC */ 725 ierr = PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_L",ilink->splitname);CHKERRQ(ierr); 726 ierr = PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr); 727 if (!LSC_L) {ierr = PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr);} 728 if (LSC_L) {ierr = PetscObjectCompose((PetscObject)jac->schur,"LSC_L",(PetscObject)LSC_L);CHKERRQ(ierr);} 729 ierr = PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_Lp",ilink->splitname);CHKERRQ(ierr); 730 ierr = PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr); 731 if (!LSC_L) {ierr = PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr);} 732 if (LSC_L) {ierr = PetscObjectCompose((PetscObject)jac->schur,"LSC_Lp",(PetscObject)LSC_L);CHKERRQ(ierr);} 733 } else { 734 /* set up the individual splits' PCs */ 735 i = 0; 736 ilink = jac->head; 737 while (ilink) { 738 ierr = KSPSetOperators(ilink->ksp,jac->mat[i],jac->pmat[i],flag);CHKERRQ(ierr); 739 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 740 if (!jac->suboptionsset) {ierr = KSPSetFromOptions(ilink->ksp);CHKERRQ(ierr);} 741 i++; 742 ilink = ilink->next; 743 } 744 } 745 746 jac->suboptionsset = PETSC_TRUE; 747 PetscFunctionReturn(0); 748 } 749 750 #define FieldSplitSplitSolveAdd(ilink,xx,yy) \ 751 (VecScatterBegin(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 752 VecScatterEnd(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 753 KSPSolve(ilink->ksp,ilink->x,ilink->y) || \ 754 VecScatterBegin(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE) || \ 755 VecScatterEnd(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE)) 756 757 #undef __FUNCT__ 758 #define __FUNCT__ "PCApply_FieldSplit_Schur" 759 static PetscErrorCode PCApply_FieldSplit_Schur(PC pc,Vec x,Vec y) 760 { 761 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 762 PetscErrorCode ierr; 763 PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next; 764 KSP kspA = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper; 765 766 PetscFunctionBegin; 767 switch (jac->schurfactorization) { 768 case PC_FIELDSPLIT_SCHUR_FACT_DIAG: 769 /* [A00 0; 0 -S], positive definite, suitable for MINRES */ 770 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 771 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 772 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 773 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 774 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 775 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 776 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 777 ierr = VecScale(ilinkD->y,-1.);CHKERRQ(ierr); 778 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 779 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 780 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 781 break; 782 case PC_FIELDSPLIT_SCHUR_FACT_LOWER: 783 /* [A00 0; A10 S], suitable for left preconditioning */ 784 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 785 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 786 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 787 ierr = MatMult(jac->C,ilinkA->y,ilinkD->x);CHKERRQ(ierr); 788 ierr = VecScale(ilinkD->x,-1.);CHKERRQ(ierr); 789 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 790 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 791 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 792 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 793 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 794 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 795 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 796 break; 797 case PC_FIELDSPLIT_SCHUR_FACT_UPPER: 798 /* [A00 A01; 0 S], suitable for right preconditioning */ 799 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 800 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 801 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 802 ierr = MatMult(jac->B,ilinkD->y,ilinkA->x);CHKERRQ(ierr); 803 ierr = VecScale(ilinkA->x,-1.);CHKERRQ(ierr); 804 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 805 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 806 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 807 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 808 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 809 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 810 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 811 break; 812 case PC_FIELDSPLIT_SCHUR_FACT_FULL: 813 /* [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 */ 814 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 815 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 816 ierr = KSPSolve(kspLower,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 817 ierr = MatMult(jac->C,ilinkA->y,ilinkD->x);CHKERRQ(ierr); 818 ierr = VecScale(ilinkD->x,-1.0);CHKERRQ(ierr); 819 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 820 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 821 822 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 823 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 824 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 825 826 if (kspUpper == kspA) { 827 ierr = MatMult(jac->B,ilinkD->y,ilinkA->y);CHKERRQ(ierr); 828 ierr = VecAXPY(ilinkA->x,-1.0,ilinkA->y);CHKERRQ(ierr); 829 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 830 } else { 831 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 832 ierr = MatMult(jac->B,ilinkD->y,ilinkA->x);CHKERRQ(ierr); 833 ierr = KSPSolve(kspUpper,ilinkA->x,ilinkA->z);CHKERRQ(ierr); 834 ierr = VecAXPY(ilinkA->y,-1.0,ilinkA->z);CHKERRQ(ierr); 835 } 836 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 837 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 838 } 839 PetscFunctionReturn(0); 840 } 841 842 #undef __FUNCT__ 843 #define __FUNCT__ "PCApply_FieldSplit" 844 static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y) 845 { 846 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 847 PetscErrorCode ierr; 848 PC_FieldSplitLink ilink = jac->head; 849 PetscInt cnt,bs; 850 851 PetscFunctionBegin; 852 if (jac->type == PC_COMPOSITE_ADDITIVE) { 853 if (jac->defaultsplit) { 854 ierr = VecGetBlockSize(x,&bs);CHKERRQ(ierr); 855 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); 856 ierr = VecGetBlockSize(y,&bs);CHKERRQ(ierr); 857 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); 858 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 859 while (ilink) { 860 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 861 ilink = ilink->next; 862 } 863 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 864 } else { 865 ierr = VecSet(y,0.0);CHKERRQ(ierr); 866 while (ilink) { 867 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 868 ilink = ilink->next; 869 } 870 } 871 } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 872 if (!jac->w1) { 873 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 874 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 875 } 876 ierr = VecSet(y,0.0);CHKERRQ(ierr); 877 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 878 cnt = 1; 879 while (ilink->next) { 880 ilink = ilink->next; 881 /* compute the residual only over the part of the vector needed */ 882 ierr = MatMult(jac->Afield[cnt++],y,ilink->x);CHKERRQ(ierr); 883 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 884 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 885 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 886 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 887 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 888 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 889 } 890 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 891 cnt -= 2; 892 while (ilink->previous) { 893 ilink = ilink->previous; 894 /* compute the residual only over the part of the vector needed */ 895 ierr = MatMult(jac->Afield[cnt--],y,ilink->x);CHKERRQ(ierr); 896 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 897 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 898 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 899 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 900 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 901 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 902 } 903 } 904 } else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Unsupported or unknown composition",(int) jac->type); 905 PetscFunctionReturn(0); 906 } 907 908 #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \ 909 (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 910 VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 911 KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) || \ 912 VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \ 913 VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE)) 914 915 #undef __FUNCT__ 916 #define __FUNCT__ "PCApplyTranspose_FieldSplit" 917 static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y) 918 { 919 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 920 PetscErrorCode ierr; 921 PC_FieldSplitLink ilink = jac->head; 922 PetscInt bs; 923 924 PetscFunctionBegin; 925 if (jac->type == PC_COMPOSITE_ADDITIVE) { 926 if (jac->defaultsplit) { 927 ierr = VecGetBlockSize(x,&bs);CHKERRQ(ierr); 928 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); 929 ierr = VecGetBlockSize(y,&bs);CHKERRQ(ierr); 930 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); 931 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 932 while (ilink) { 933 ierr = KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 934 ilink = ilink->next; 935 } 936 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 937 } else { 938 ierr = VecSet(y,0.0);CHKERRQ(ierr); 939 while (ilink) { 940 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 941 ilink = ilink->next; 942 } 943 } 944 } else { 945 if (!jac->w1) { 946 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 947 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 948 } 949 ierr = VecSet(y,0.0);CHKERRQ(ierr); 950 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 951 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 952 while (ilink->next) { 953 ilink = ilink->next; 954 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 955 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 956 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 957 } 958 while (ilink->previous) { 959 ilink = ilink->previous; 960 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 961 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 962 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 963 } 964 } else { 965 while (ilink->next) { /* get to last entry in linked list */ 966 ilink = ilink->next; 967 } 968 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 969 while (ilink->previous) { 970 ilink = ilink->previous; 971 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 972 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 973 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 974 } 975 } 976 } 977 PetscFunctionReturn(0); 978 } 979 980 #undef __FUNCT__ 981 #define __FUNCT__ "PCReset_FieldSplit" 982 static PetscErrorCode PCReset_FieldSplit(PC pc) 983 { 984 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 985 PetscErrorCode ierr; 986 PC_FieldSplitLink ilink = jac->head,next; 987 988 PetscFunctionBegin; 989 while (ilink) { 990 ierr = KSPReset(ilink->ksp);CHKERRQ(ierr); 991 ierr = VecDestroy(&ilink->x);CHKERRQ(ierr); 992 ierr = VecDestroy(&ilink->y);CHKERRQ(ierr); 993 ierr = VecDestroy(&ilink->z);CHKERRQ(ierr); 994 ierr = VecScatterDestroy(&ilink->sctx);CHKERRQ(ierr); 995 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 996 ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr); 997 next = ilink->next; 998 ilink = next; 999 } 1000 ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr); 1001 if (jac->mat && jac->mat != jac->pmat) { 1002 ierr = MatDestroyMatrices(jac->nsplits,&jac->mat);CHKERRQ(ierr); 1003 } else if (jac->mat) { 1004 jac->mat = NULL; 1005 } 1006 if (jac->pmat) {ierr = MatDestroyMatrices(jac->nsplits,&jac->pmat);CHKERRQ(ierr);} 1007 if (jac->Afield) {ierr = MatDestroyMatrices(jac->nsplits,&jac->Afield);CHKERRQ(ierr);} 1008 ierr = VecDestroy(&jac->w1);CHKERRQ(ierr); 1009 ierr = VecDestroy(&jac->w2);CHKERRQ(ierr); 1010 ierr = MatDestroy(&jac->schur);CHKERRQ(ierr); 1011 ierr = MatDestroy(&jac->schur_user);CHKERRQ(ierr); 1012 ierr = KSPDestroy(&jac->kspschur);CHKERRQ(ierr); 1013 ierr = KSPDestroy(&jac->kspupper);CHKERRQ(ierr); 1014 ierr = MatDestroy(&jac->B);CHKERRQ(ierr); 1015 ierr = MatDestroy(&jac->C);CHKERRQ(ierr); 1016 jac->reset = PETSC_TRUE; 1017 PetscFunctionReturn(0); 1018 } 1019 1020 #undef __FUNCT__ 1021 #define __FUNCT__ "PCDestroy_FieldSplit" 1022 static PetscErrorCode PCDestroy_FieldSplit(PC pc) 1023 { 1024 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1025 PetscErrorCode ierr; 1026 PC_FieldSplitLink ilink = jac->head,next; 1027 1028 PetscFunctionBegin; 1029 ierr = PCReset_FieldSplit(pc);CHKERRQ(ierr); 1030 while (ilink) { 1031 ierr = KSPDestroy(&ilink->ksp);CHKERRQ(ierr); 1032 next = ilink->next; 1033 ierr = PetscFree(ilink->splitname);CHKERRQ(ierr); 1034 ierr = PetscFree(ilink->fields);CHKERRQ(ierr); 1035 ierr = PetscFree(ilink->fields_col);CHKERRQ(ierr); 1036 ierr = PetscFree(ilink);CHKERRQ(ierr); 1037 ilink = next; 1038 } 1039 ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr); 1040 ierr = PetscFree(pc->data);CHKERRQ(ierr); 1041 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",NULL);CHKERRQ(ierr); 1042 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",NULL);CHKERRQ(ierr); 1043 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",NULL);CHKERRQ(ierr); 1044 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",NULL);CHKERRQ(ierr); 1045 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",NULL);CHKERRQ(ierr); 1046 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSchurPrecondition_C",NULL);CHKERRQ(ierr); 1047 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",NULL);CHKERRQ(ierr); 1048 PetscFunctionReturn(0); 1049 } 1050 1051 #undef __FUNCT__ 1052 #define __FUNCT__ "PCSetFromOptions_FieldSplit" 1053 static PetscErrorCode PCSetFromOptions_FieldSplit(PC pc) 1054 { 1055 PetscErrorCode ierr; 1056 PetscInt bs; 1057 PetscBool flg,stokes = PETSC_FALSE; 1058 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1059 PCCompositeType ctype; 1060 1061 PetscFunctionBegin; 1062 ierr = PetscOptionsHead("FieldSplit options");CHKERRQ(ierr); 1063 ierr = PetscOptionsBool("-pc_fieldsplit_dm_splits","Whether to use DMCreateFieldDecomposition() for splits","PCFieldSplitSetDMSplits",jac->dm_splits,&jac->dm_splits,NULL);CHKERRQ(ierr); 1064 ierr = PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);CHKERRQ(ierr); 1065 if (flg) { 1066 ierr = PCFieldSplitSetBlockSize(pc,bs);CHKERRQ(ierr); 1067 } 1068 1069 ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,NULL);CHKERRQ(ierr); 1070 if (stokes) { 1071 ierr = PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);CHKERRQ(ierr); 1072 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_SELF; 1073 } 1074 1075 ierr = PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg);CHKERRQ(ierr); 1076 if (flg) { 1077 ierr = PCFieldSplitSetType(pc,ctype);CHKERRQ(ierr); 1078 } 1079 /* Only setup fields once */ 1080 if ((jac->bs > 0) && (jac->nsplits == 0)) { 1081 /* only allow user to set fields from command line if bs is already known. 1082 otherwise user can set them in PCFieldSplitSetDefaults() */ 1083 ierr = PCFieldSplitSetRuntimeSplits_Private(pc);CHKERRQ(ierr); 1084 if (jac->splitdefined) {ierr = PetscInfo(pc,"Splits defined using the options database\n");CHKERRQ(ierr);} 1085 } 1086 if (jac->type == PC_COMPOSITE_SCHUR) { 1087 ierr = PetscOptionsGetEnum(((PetscObject)pc)->prefix,"-pc_fieldsplit_schur_factorization_type",PCFieldSplitSchurFactTypes,(PetscEnum*)&jac->schurfactorization,&flg);CHKERRQ(ierr); 1088 if (flg) {ierr = PetscInfo(pc,"Deprecated use of -pc_fieldsplit_schur_factorization_type\n");CHKERRQ(ierr);} 1089 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); 1090 ierr = PetscOptionsEnum("-pc_fieldsplit_schur_precondition","How to build preconditioner for Schur complement","PCFieldSplitSchurPrecondition",PCFieldSplitSchurPreTypes,(PetscEnum)jac->schurpre,(PetscEnum*)&jac->schurpre,NULL);CHKERRQ(ierr); 1091 } 1092 ierr = PetscOptionsTail();CHKERRQ(ierr); 1093 PetscFunctionReturn(0); 1094 } 1095 1096 /*------------------------------------------------------------------------------------*/ 1097 1098 #undef __FUNCT__ 1099 #define __FUNCT__ "PCFieldSplitSetFields_FieldSplit" 1100 static PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 1101 { 1102 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1103 PetscErrorCode ierr; 1104 PC_FieldSplitLink ilink,next = jac->head; 1105 char prefix[128]; 1106 PetscInt i; 1107 1108 PetscFunctionBegin; 1109 if (jac->splitdefined) { 1110 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 1111 PetscFunctionReturn(0); 1112 } 1113 for (i=0; i<n; i++) { 1114 if (fields[i] >= jac->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs); 1115 if (fields[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]); 1116 } 1117 ierr = PetscNew(struct _PC_FieldSplitLink,&ilink);CHKERRQ(ierr); 1118 if (splitname) { 1119 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 1120 } else { 1121 ierr = PetscMalloc(3*sizeof(char),&ilink->splitname);CHKERRQ(ierr); 1122 ierr = PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);CHKERRQ(ierr); 1123 } 1124 ierr = PetscMalloc(n*sizeof(PetscInt),&ilink->fields);CHKERRQ(ierr); 1125 ierr = PetscMemcpy(ilink->fields,fields,n*sizeof(PetscInt));CHKERRQ(ierr); 1126 ierr = PetscMalloc(n*sizeof(PetscInt),&ilink->fields_col);CHKERRQ(ierr); 1127 ierr = PetscMemcpy(ilink->fields_col,fields_col,n*sizeof(PetscInt));CHKERRQ(ierr); 1128 1129 ilink->nfields = n; 1130 ilink->next = NULL; 1131 ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);CHKERRQ(ierr); 1132 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 1133 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 1134 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 1135 1136 ierr = PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);CHKERRQ(ierr); 1137 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 1138 1139 if (!next) { 1140 jac->head = ilink; 1141 ilink->previous = NULL; 1142 } else { 1143 while (next->next) { 1144 next = next->next; 1145 } 1146 next->next = ilink; 1147 ilink->previous = next; 1148 } 1149 jac->nsplits++; 1150 PetscFunctionReturn(0); 1151 } 1152 1153 #undef __FUNCT__ 1154 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit_Schur" 1155 static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp) 1156 { 1157 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1158 PetscErrorCode ierr; 1159 1160 PetscFunctionBegin; 1161 ierr = PetscMalloc(jac->nsplits*sizeof(KSP),subksp);CHKERRQ(ierr); 1162 ierr = MatSchurComplementGetKSP(jac->schur,*subksp);CHKERRQ(ierr); 1163 1164 (*subksp)[1] = jac->kspschur; 1165 if (n) *n = jac->nsplits; 1166 PetscFunctionReturn(0); 1167 } 1168 1169 #undef __FUNCT__ 1170 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit" 1171 static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp) 1172 { 1173 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1174 PetscErrorCode ierr; 1175 PetscInt cnt = 0; 1176 PC_FieldSplitLink ilink = jac->head; 1177 1178 PetscFunctionBegin; 1179 ierr = PetscMalloc(jac->nsplits*sizeof(KSP),subksp);CHKERRQ(ierr); 1180 while (ilink) { 1181 (*subksp)[cnt++] = ilink->ksp; 1182 ilink = ilink->next; 1183 } 1184 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); 1185 if (n) *n = jac->nsplits; 1186 PetscFunctionReturn(0); 1187 } 1188 1189 #undef __FUNCT__ 1190 #define __FUNCT__ "PCFieldSplitSetIS_FieldSplit" 1191 static PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is) 1192 { 1193 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1194 PetscErrorCode ierr; 1195 PC_FieldSplitLink ilink, next = jac->head; 1196 char prefix[128]; 1197 1198 PetscFunctionBegin; 1199 if (jac->splitdefined) { 1200 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 1201 PetscFunctionReturn(0); 1202 } 1203 ierr = PetscNew(struct _PC_FieldSplitLink,&ilink);CHKERRQ(ierr); 1204 if (splitname) { 1205 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 1206 } else { 1207 ierr = PetscMalloc(8*sizeof(char),&ilink->splitname);CHKERRQ(ierr); 1208 ierr = PetscSNPrintf(ilink->splitname,7,"%D",jac->nsplits);CHKERRQ(ierr); 1209 } 1210 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 1211 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 1212 ilink->is = is; 1213 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 1214 ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr); 1215 ilink->is_col = is; 1216 ilink->next = NULL; 1217 ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);CHKERRQ(ierr); 1218 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 1219 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 1220 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 1221 1222 ierr = PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);CHKERRQ(ierr); 1223 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 1224 1225 if (!next) { 1226 jac->head = ilink; 1227 ilink->previous = NULL; 1228 } else { 1229 while (next->next) { 1230 next = next->next; 1231 } 1232 next->next = ilink; 1233 ilink->previous = next; 1234 } 1235 jac->nsplits++; 1236 PetscFunctionReturn(0); 1237 } 1238 1239 #undef __FUNCT__ 1240 #define __FUNCT__ "PCFieldSplitSetFields" 1241 /*@ 1242 PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner 1243 1244 Logically Collective on PC 1245 1246 Input Parameters: 1247 + pc - the preconditioner context 1248 . splitname - name of this split, if NULL the number of the split is used 1249 . n - the number of fields in this split 1250 - fields - the fields in this split 1251 1252 Level: intermediate 1253 1254 Notes: Use PCFieldSplitSetIS() to set a completely general set of indices as a field. 1255 1256 The PCFieldSplitSetFields() is for defining fields as strided blocks. For example, if the block 1257 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 1258 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x.... 1259 where the numbered entries indicate what is in the field. 1260 1261 This function is called once per split (it creates a new split each time). Solve options 1262 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1263 1264 Developer Note: This routine does not actually create the IS representing the split, that is delayed 1265 until PCSetUp_FieldSplit(), because information about the vector/matrix layouts may not be 1266 available when this routine is called. 1267 1268 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize(), PCFieldSplitSetIS() 1269 1270 @*/ 1271 PetscErrorCode PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 1272 { 1273 PetscErrorCode ierr; 1274 1275 PetscFunctionBegin; 1276 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1277 PetscValidCharPointer(splitname,2); 1278 if (n < 1) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Provided number of fields %D in split \"%s\" not positive",n,splitname); 1279 PetscValidIntPointer(fields,3); 1280 ierr = PetscTryMethod(pc,"PCFieldSplitSetFields_C",(PC,const char[],PetscInt,const PetscInt*,const PetscInt*),(pc,splitname,n,fields,fields_col));CHKERRQ(ierr); 1281 PetscFunctionReturn(0); 1282 } 1283 1284 #undef __FUNCT__ 1285 #define __FUNCT__ "PCFieldSplitSetIS" 1286 /*@ 1287 PCFieldSplitSetIS - Sets the exact elements for field 1288 1289 Logically Collective on PC 1290 1291 Input Parameters: 1292 + pc - the preconditioner context 1293 . splitname - name of this split, if NULL the number of the split is used 1294 - is - the index set that defines the vector elements in this field 1295 1296 1297 Notes: 1298 Use PCFieldSplitSetFields(), for fields defined by strided types. 1299 1300 This function is called once per split (it creates a new split each time). Solve options 1301 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1302 1303 Level: intermediate 1304 1305 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize() 1306 1307 @*/ 1308 PetscErrorCode PCFieldSplitSetIS(PC pc,const char splitname[],IS is) 1309 { 1310 PetscErrorCode ierr; 1311 1312 PetscFunctionBegin; 1313 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1314 if (splitname) PetscValidCharPointer(splitname,2); 1315 PetscValidHeaderSpecific(is,IS_CLASSID,3); 1316 ierr = PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));CHKERRQ(ierr); 1317 PetscFunctionReturn(0); 1318 } 1319 1320 #undef __FUNCT__ 1321 #define __FUNCT__ "PCFieldSplitGetIS" 1322 /*@ 1323 PCFieldSplitGetIS - Retrieves the elements for a field as an IS 1324 1325 Logically Collective on PC 1326 1327 Input Parameters: 1328 + pc - the preconditioner context 1329 - splitname - name of this split 1330 1331 Output Parameter: 1332 - is - the index set that defines the vector elements in this field, or NULL if the field is not found 1333 1334 Level: intermediate 1335 1336 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetIS() 1337 1338 @*/ 1339 PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is) 1340 { 1341 PetscErrorCode ierr; 1342 1343 PetscFunctionBegin; 1344 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1345 PetscValidCharPointer(splitname,2); 1346 PetscValidPointer(is,3); 1347 { 1348 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 1349 PC_FieldSplitLink ilink = jac->head; 1350 PetscBool found; 1351 1352 *is = NULL; 1353 while (ilink) { 1354 ierr = PetscStrcmp(ilink->splitname, splitname, &found);CHKERRQ(ierr); 1355 if (found) { 1356 *is = ilink->is; 1357 break; 1358 } 1359 ilink = ilink->next; 1360 } 1361 } 1362 PetscFunctionReturn(0); 1363 } 1364 1365 #undef __FUNCT__ 1366 #define __FUNCT__ "PCFieldSplitSetBlockSize" 1367 /*@ 1368 PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the 1369 fieldsplit preconditioner. If not set the matrix block size is used. 1370 1371 Logically Collective on PC 1372 1373 Input Parameters: 1374 + pc - the preconditioner context 1375 - bs - the block size 1376 1377 Level: intermediate 1378 1379 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields() 1380 1381 @*/ 1382 PetscErrorCode PCFieldSplitSetBlockSize(PC pc,PetscInt bs) 1383 { 1384 PetscErrorCode ierr; 1385 1386 PetscFunctionBegin; 1387 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1388 PetscValidLogicalCollectiveInt(pc,bs,2); 1389 ierr = PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));CHKERRQ(ierr); 1390 PetscFunctionReturn(0); 1391 } 1392 1393 #undef __FUNCT__ 1394 #define __FUNCT__ "PCFieldSplitGetSubKSP" 1395 /*@C 1396 PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits 1397 1398 Collective on KSP 1399 1400 Input Parameter: 1401 . pc - the preconditioner context 1402 1403 Output Parameters: 1404 + n - the number of splits 1405 - pc - the array of KSP contexts 1406 1407 Note: 1408 After PCFieldSplitGetSubKSP() the array of KSPs IS to be freed by the user 1409 (not the KSP just the array that contains them). 1410 1411 You must call KSPSetUp() before calling PCFieldSplitGetSubKSP(). 1412 1413 Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs. 1414 You can call PCFieldSplitGetSubKSP(pc,n,NULL_OBJECT,ierr) to determine how large the 1415 KSP array must be. 1416 1417 1418 Level: advanced 1419 1420 .seealso: PCFIELDSPLIT 1421 @*/ 1422 PetscErrorCode PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[]) 1423 { 1424 PetscErrorCode ierr; 1425 1426 PetscFunctionBegin; 1427 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1428 if (n) PetscValidIntPointer(n,2); 1429 ierr = PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));CHKERRQ(ierr); 1430 PetscFunctionReturn(0); 1431 } 1432 1433 #undef __FUNCT__ 1434 #define __FUNCT__ "PCFieldSplitSchurPrecondition" 1435 /*@ 1436 PCFieldSplitSchurPrecondition - Indicates if the Schur complement is preconditioned by a preconditioner constructed by the 1437 A11 matrix. Otherwise no preconditioner is used. 1438 1439 Collective on PC 1440 1441 Input Parameters: 1442 + pc - the preconditioner context 1443 . ptype - which matrix to use for preconditioning the Schur complement, PC_FIELDSPLIT_SCHUR_PRE_A11 (diag) is default 1444 - userpre - matrix to use for preconditioning, or NULL 1445 1446 Options Database: 1447 . -pc_fieldsplit_schur_precondition <self,user,a11> default is a11 1448 1449 Notes: 1450 $ If ptype is 1451 $ user then the preconditioner for the Schur complement is generated by the provided matrix (pre argument 1452 $ to this function). 1453 $ a11 then the preconditioner for the Schur complement is generated by the block diagonal part of the original 1454 $ matrix associated with the Schur complement (i.e. A11) 1455 $ self the preconditioner for the Schur complement is generated from the Schur complement matrix itself: 1456 $ The only preconditioner that currently works directly with the Schur complement matrix object is the PCLSC 1457 $ preconditioner 1458 1459 When solving a saddle point problem, where the A11 block is identically zero, using a11 as the ptype only makes sense 1460 with the additional option -fieldsplit_1_pc_type none. Usually for saddle point problems one would use a ptype of self and 1461 -fieldsplit_1_pc_type lsc which uses the least squares commutator compute a preconditioner for the Schur complement. 1462 1463 Developer Notes: This is a terrible name, gives no good indication of what the function does and should also have Set in 1464 the name since it sets a proceedure to use. 1465 1466 Level: intermediate 1467 1468 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCLSC 1469 1470 @*/ 1471 PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1472 { 1473 PetscErrorCode ierr; 1474 1475 PetscFunctionBegin; 1476 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1477 ierr = PetscTryMethod(pc,"PCFieldSplitSchurPrecondition_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));CHKERRQ(ierr); 1478 PetscFunctionReturn(0); 1479 } 1480 1481 #undef __FUNCT__ 1482 #define __FUNCT__ "PCFieldSplitSchurPrecondition_FieldSplit" 1483 static PetscErrorCode PCFieldSplitSchurPrecondition_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1484 { 1485 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1486 PetscErrorCode ierr; 1487 1488 PetscFunctionBegin; 1489 jac->schurpre = ptype; 1490 if (pre) { 1491 ierr = MatDestroy(&jac->schur_user);CHKERRQ(ierr); 1492 jac->schur_user = pre; 1493 ierr = PetscObjectReference((PetscObject)jac->schur_user);CHKERRQ(ierr); 1494 } 1495 PetscFunctionReturn(0); 1496 } 1497 1498 #undef __FUNCT__ 1499 #define __FUNCT__ "PCFieldSplitSetSchurFactType" 1500 /*@ 1501 PCFieldSplitSetSchurFactType - sets which blocks of the approximate block factorization to retain 1502 1503 Collective on PC 1504 1505 Input Parameters: 1506 + pc - the preconditioner context 1507 - ftype - which blocks of factorization to retain, PC_FIELDSPLIT_SCHUR_FACT_FULL is default 1508 1509 Options Database: 1510 . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> default is full 1511 1512 1513 Level: intermediate 1514 1515 Notes: 1516 The FULL factorization is 1517 1518 $ (A B) = (1 0) (A 0) (1 Ainv*B) 1519 $ (C D) (C*Ainv 1) (0 S) (0 1 ) 1520 1521 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, 1522 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). 1523 1524 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 1525 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 1526 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, 1527 the preconditioned operator has three distinct nonzero eigenvalues and minimal polynomial of degree at most 4, so KSPGMRES converges in at most 4 iterations. 1528 1529 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 1530 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). 1531 1532 References: 1533 Murphy, Golub, and Wathen, A note on preconditioning indefinite linear systems, SIAM J. Sci. Comput., 21 (2000) pp. 1969-1972. 1534 1535 Ipsen, A note on preconditioning nonsymmetric matrices, SIAM J. Sci. Comput., 23 (2001), pp. 1050-1051. 1536 1537 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType 1538 @*/ 1539 PetscErrorCode PCFieldSplitSetSchurFactType(PC pc,PCFieldSplitSchurFactType ftype) 1540 { 1541 PetscErrorCode ierr; 1542 1543 PetscFunctionBegin; 1544 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1545 ierr = PetscTryMethod(pc,"PCFieldSplitSetSchurFactType_C",(PC,PCFieldSplitSchurFactType),(pc,ftype));CHKERRQ(ierr); 1546 PetscFunctionReturn(0); 1547 } 1548 1549 #undef __FUNCT__ 1550 #define __FUNCT__ "PCFieldSplitSetSchurFactType_FieldSplit" 1551 static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc,PCFieldSplitSchurFactType ftype) 1552 { 1553 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1554 1555 PetscFunctionBegin; 1556 jac->schurfactorization = ftype; 1557 PetscFunctionReturn(0); 1558 } 1559 1560 #undef __FUNCT__ 1561 #define __FUNCT__ "PCFieldSplitGetSchurBlocks" 1562 /*@C 1563 PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement 1564 1565 Collective on KSP 1566 1567 Input Parameter: 1568 . pc - the preconditioner context 1569 1570 Output Parameters: 1571 + A00 - the (0,0) block 1572 . A01 - the (0,1) block 1573 . A10 - the (1,0) block 1574 - A11 - the (1,1) block 1575 1576 Level: advanced 1577 1578 .seealso: PCFIELDSPLIT 1579 @*/ 1580 PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11) 1581 { 1582 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 1583 1584 PetscFunctionBegin; 1585 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1586 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach."); 1587 if (A00) *A00 = jac->pmat[0]; 1588 if (A01) *A01 = jac->B; 1589 if (A10) *A10 = jac->C; 1590 if (A11) *A11 = jac->pmat[1]; 1591 PetscFunctionReturn(0); 1592 } 1593 1594 #undef __FUNCT__ 1595 #define __FUNCT__ "PCFieldSplitSetType_FieldSplit" 1596 static PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type) 1597 { 1598 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1599 PetscErrorCode ierr; 1600 1601 PetscFunctionBegin; 1602 jac->type = type; 1603 if (type == PC_COMPOSITE_SCHUR) { 1604 pc->ops->apply = PCApply_FieldSplit_Schur; 1605 pc->ops->view = PCView_FieldSplit_Schur; 1606 1607 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit_Schur);CHKERRQ(ierr); 1608 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSchurPrecondition_C",PCFieldSplitSchurPrecondition_FieldSplit);CHKERRQ(ierr); 1609 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",PCFieldSplitSetSchurFactType_FieldSplit);CHKERRQ(ierr); 1610 1611 } else { 1612 pc->ops->apply = PCApply_FieldSplit; 1613 pc->ops->view = PCView_FieldSplit; 1614 1615 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 1616 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSchurPrecondition_C",0);CHKERRQ(ierr); 1617 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",0);CHKERRQ(ierr); 1618 } 1619 PetscFunctionReturn(0); 1620 } 1621 1622 #undef __FUNCT__ 1623 #define __FUNCT__ "PCFieldSplitSetBlockSize_FieldSplit" 1624 static PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs) 1625 { 1626 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1627 1628 PetscFunctionBegin; 1629 if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs); 1630 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); 1631 jac->bs = bs; 1632 PetscFunctionReturn(0); 1633 } 1634 1635 #undef __FUNCT__ 1636 #define __FUNCT__ "PCFieldSplitSetType" 1637 /*@ 1638 PCFieldSplitSetType - Sets the type of fieldsplit preconditioner. 1639 1640 Collective on PC 1641 1642 Input Parameter: 1643 . pc - the preconditioner context 1644 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 1645 1646 Options Database Key: 1647 . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type 1648 1649 Level: Intermediate 1650 1651 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 1652 1653 .seealso: PCCompositeSetType() 1654 1655 @*/ 1656 PetscErrorCode PCFieldSplitSetType(PC pc,PCCompositeType type) 1657 { 1658 PetscErrorCode ierr; 1659 1660 PetscFunctionBegin; 1661 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1662 ierr = PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));CHKERRQ(ierr); 1663 PetscFunctionReturn(0); 1664 } 1665 1666 #undef __FUNCT__ 1667 #define __FUNCT__ "PCFieldSplitGetType" 1668 /*@ 1669 PCFieldSplitGetType - Gets the type of fieldsplit preconditioner. 1670 1671 Not collective 1672 1673 Input Parameter: 1674 . pc - the preconditioner context 1675 1676 Output Parameter: 1677 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 1678 1679 Level: Intermediate 1680 1681 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 1682 .seealso: PCCompositeSetType() 1683 @*/ 1684 PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type) 1685 { 1686 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 1687 1688 PetscFunctionBegin; 1689 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1690 PetscValidIntPointer(type,2); 1691 *type = jac->type; 1692 PetscFunctionReturn(0); 1693 } 1694 1695 #undef __FUNCT__ 1696 #define __FUNCT__ "PCFieldSplitSetDMSplits" 1697 /*@ 1698 PCFieldSplitSetDMSplits - Flags whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible. 1699 1700 Logically Collective 1701 1702 Input Parameters: 1703 + pc - the preconditioner context 1704 - flg - boolean indicating whether to use field splits defined by the DM 1705 1706 Options Database Key: 1707 . -pc_fieldsplit_dm_splits 1708 1709 Level: Intermediate 1710 1711 .keywords: PC, DM, composite preconditioner, additive, multiplicative 1712 1713 .seealso: PCFieldSplitGetDMSplits() 1714 1715 @*/ 1716 PetscErrorCode PCFieldSplitSetDMSplits(PC pc,PetscBool flg) 1717 { 1718 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1719 PetscBool isfs; 1720 PetscErrorCode ierr; 1721 1722 PetscFunctionBegin; 1723 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1724 PetscValidLogicalCollectiveBool(pc,flg,2); 1725 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1726 if (isfs) { 1727 jac->dm_splits = flg; 1728 } 1729 PetscFunctionReturn(0); 1730 } 1731 1732 1733 #undef __FUNCT__ 1734 #define __FUNCT__ "PCFieldSplitGetDMSplits" 1735 /*@ 1736 PCFieldSplitGetDMSplits - Returns flag indicating whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible. 1737 1738 Logically Collective 1739 1740 Input Parameter: 1741 . pc - the preconditioner context 1742 1743 Output Parameter: 1744 . flg - boolean indicating whether to use field splits defined by the DM 1745 1746 Level: Intermediate 1747 1748 .keywords: PC, DM, composite preconditioner, additive, multiplicative 1749 1750 .seealso: PCFieldSplitSetDMSplits() 1751 1752 @*/ 1753 PetscErrorCode PCFieldSplitGetDMSplits(PC pc,PetscBool* flg) 1754 { 1755 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1756 PetscBool isfs; 1757 PetscErrorCode ierr; 1758 1759 PetscFunctionBegin; 1760 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1761 PetscValidPointer(flg,2); 1762 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1763 if (isfs) { 1764 if(flg) *flg = jac->dm_splits; 1765 } 1766 PetscFunctionReturn(0); 1767 } 1768 1769 /* -------------------------------------------------------------------------------------*/ 1770 /*MC 1771 PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual 1772 fields or groups of fields. See the users manual section "Solving Block Matrices" for more details. 1773 1774 To set options on the solvers for each block append -fieldsplit_ to all the PC 1775 options database keys. For example, -fieldsplit_pc_type ilu -fieldsplit_pc_factor_levels 1 1776 1777 To set the options on the solvers separate for each block call PCFieldSplitGetSubKSP() 1778 and set the options directly on the resulting KSP object 1779 1780 Level: intermediate 1781 1782 Options Database Keys: 1783 + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split 1784 . -pc_fieldsplit_default - automatically add any fields to additional splits that have not 1785 been supplied explicitly by -pc_fieldsplit_%d_fields 1786 . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields) 1787 . -pc_fieldsplit_type <additive,multiplicative,symmetric_multiplicative,schur> - type of relaxation or factorization splitting 1788 . -pc_fieldsplit_schur_precondition <self,user,a11> - default is a11 1789 . -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero or negative diagonal and uses Schur complement with no preconditioner as the solver 1790 1791 - Options prefix for inner solvers when using Schur complement preconditioner are -fieldsplit_0_ and -fieldsplit_1_ 1792 for all other solvers they are -fieldsplit_%d_ for the dth field, use -fieldsplit_ for all fields 1793 1794 Notes: 1795 Use PCFieldSplitSetFields() to set fields defined by "strided" entries and PCFieldSplitSetIS() 1796 to define a field by an arbitrary collection of entries. 1797 1798 If no fields are set the default is used. The fields are defined by entries strided by bs, 1799 beginning at 0 then 1, etc to bs-1. The block size can be set with PCFieldSplitSetBlockSize(), 1800 if this is not called the block size defaults to the blocksize of the second matrix passed 1801 to KSPSetOperators()/PCSetOperators(). 1802 1803 $ For the Schur complement preconditioner if J = ( A00 A01 ) 1804 $ ( A10 A11 ) 1805 $ the preconditioner using full factorization is 1806 $ ( I -A10 ksp(A00) ) ( inv(A00) 0 ) ( I 0 ) 1807 $ ( 0 I ) ( 0 ksp(S) ) ( -A10 ksp(A00) I ) 1808 where the action of inv(A00) is applied using the KSP solver with prefix -fieldsplit_0_. The action of 1809 ksp(S) is computed using the KSP solver with prefix -fieldsplit_splitname_ (where splitname was given 1810 in providing the SECOND split or 1 if not give). For PCFieldSplitGetKSP() when field number is 0, 1811 it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default 1812 A11 is used to construct a preconditioner for S, use PCFieldSplitSchurPrecondition() to turn on or off this 1813 option. You can use the preconditioner PCLSC to precondition the Schur complement with -fieldsplit_1_pc_type lsc. The 1814 factorization type is set using -pc_fieldsplit_schur_fact_type <diag, lower, upper, full>. The full is shown above, 1815 diag gives 1816 $ ( inv(A00) 0 ) 1817 $ ( 0 -ksp(S) ) 1818 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 1819 $ ( A00 0 ) 1820 $ ( A10 S ) 1821 where the inverses of A00 and S are applied using KSPs. The upper factorization is the inverse of 1822 $ ( A00 A01 ) 1823 $ ( 0 S ) 1824 where again the inverses of A00 and S are applied using KSPs. 1825 1826 If only one set of indices (one IS) is provided with PCFieldSplitSetIS() then the complement of that IS 1827 is used automatically for a second block. 1828 1829 The fieldsplit preconditioner cannot currently be used with the BAIJ or SBAIJ data formats if the blocksize is larger than 1. 1830 Generally it should be used with the AIJ format. 1831 1832 The forms of these preconditioners are closely related if not identical to forms derived as "Distributive Iterations", see, 1833 for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling. Note that one can also use PCFIELDSPLIT 1834 inside a smoother resulting in "Distributive Smoothers". 1835 1836 Concepts: physics based preconditioners, block preconditioners 1837 1838 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners, PCLSC, 1839 PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSchurPrecondition() 1840 M*/ 1841 1842 #undef __FUNCT__ 1843 #define __FUNCT__ "PCCreate_FieldSplit" 1844 PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc) 1845 { 1846 PetscErrorCode ierr; 1847 PC_FieldSplit *jac; 1848 1849 PetscFunctionBegin; 1850 ierr = PetscNewLog(pc,PC_FieldSplit,&jac);CHKERRQ(ierr); 1851 1852 jac->bs = -1; 1853 jac->nsplits = 0; 1854 jac->type = PC_COMPOSITE_MULTIPLICATIVE; 1855 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */ 1856 jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL; 1857 jac->dm_splits = PETSC_TRUE; 1858 1859 pc->data = (void*)jac; 1860 1861 pc->ops->apply = PCApply_FieldSplit; 1862 pc->ops->applytranspose = PCApplyTranspose_FieldSplit; 1863 pc->ops->setup = PCSetUp_FieldSplit; 1864 pc->ops->reset = PCReset_FieldSplit; 1865 pc->ops->destroy = PCDestroy_FieldSplit; 1866 pc->ops->setfromoptions = PCSetFromOptions_FieldSplit; 1867 pc->ops->view = PCView_FieldSplit; 1868 pc->ops->applyrichardson = 0; 1869 1870 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 1871 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",PCFieldSplitSetFields_FieldSplit);CHKERRQ(ierr); 1872 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",PCFieldSplitSetIS_FieldSplit);CHKERRQ(ierr); 1873 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",PCFieldSplitSetType_FieldSplit);CHKERRQ(ierr); 1874 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",PCFieldSplitSetBlockSize_FieldSplit);CHKERRQ(ierr); 1875 PetscFunctionReturn(0); 1876 } 1877 1878 1879 1880