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