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