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