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