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