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 PetscBool schurp_lump; /* Whether to lump A00 when forming Sp. */ 46 Mat schur_user; /* User-provided preconditioning matrix for the Schur complement */ 47 PCFieldSplitSchurPreType schurpre; /* Determines which preconditioning matrix is used for the Schur complement */ 48 PCFieldSplitSchurFactType schurfactorization; 49 KSP kspschur; /* The solver for S */ 50 KSP kspupper; /* The solver for A in the upper diagonal part of the factorization (H_2 in [El08]) */ 51 PC_FieldSplitLink head; 52 PetscBool reset; /* indicates PCReset() has been last called on this object, hack */ 53 PetscBool suboptionsset; /* Indicates that the KSPSetFromOptions() has been called on the sub-KSPs */ 54 PetscBool dm_splits; /* Whether to use DMCreateFieldDecomposition() whenever possible */ 55 } PC_FieldSplit; 56 57 /* 58 Notes: there is no particular reason that pmat, x, and y are stored as arrays in PC_FieldSplit instead of 59 inside PC_FieldSplitLink, just historical. If you want to be able to add new fields after already using the 60 PC you could change this. 61 */ 62 63 /* This helper is so that setting a user-provided preconditioning matrix is orthogonal to choosing to use it. This way the 64 * application-provided FormJacobian can provide this matrix without interfering with the user's (command-line) choices. */ 65 static Mat FieldSplitSchurPre(PC_FieldSplit *jac) 66 { 67 switch (jac->schurpre) { 68 case PC_FIELDSPLIT_SCHUR_PRE_SELF: return jac->schur; 69 case PC_FIELDSPLIT_SCHUR_PRE_SELFP: return jac->schurp; 70 case PC_FIELDSPLIT_SCHUR_PRE_A11: return jac->pmat[1]; 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_USER: 174 if (jac->schur_user) { 175 ierr = PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from user provided matrix\n");CHKERRQ(ierr); 176 } else { 177 ierr = PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from A11\n");CHKERRQ(ierr); 178 } 179 break; 180 default: 181 SETERRQ1(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Invalid Schur preconditioning type: %d", jac->schurpre); 182 } 183 ierr = PetscViewerASCIIPrintf(viewer," Split info:\n");CHKERRQ(ierr); 184 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 185 for (i=0; i<jac->nsplits; i++) { 186 if (ilink->fields) { 187 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);CHKERRQ(ierr); 188 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 189 for (j=0; j<ilink->nfields; j++) { 190 if (j > 0) { 191 ierr = PetscViewerASCIIPrintf(viewer,",");CHKERRQ(ierr); 192 } 193 ierr = PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);CHKERRQ(ierr); 194 } 195 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 196 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 197 } else { 198 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i);CHKERRQ(ierr); 199 } 200 ilink = ilink->next; 201 } 202 ierr = PetscViewerASCIIPrintf(viewer,"KSP solver for A00 block \n");CHKERRQ(ierr); 203 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 204 ierr = KSPView(jac->head->ksp,viewer);CHKERRQ(ierr); 205 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 206 if (jac->kspupper != jac->head->ksp) { 207 ierr = PetscViewerASCIIPrintf(viewer,"KSP solver for upper A00 in upper triangular factor \n");CHKERRQ(ierr); 208 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 209 if (jac->kspupper) {ierr = KSPView(jac->kspupper,viewer);CHKERRQ(ierr);} 210 else {ierr = PetscViewerASCIIPrintf(viewer," not yet available\n");CHKERRQ(ierr);} 211 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 212 } 213 ierr = PetscViewerASCIIPrintf(viewer,"KSP solver for S = A11 - A10 inv(A00) A01 \n");CHKERRQ(ierr); 214 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 215 if (jac->kspschur) { 216 ierr = KSPView(jac->kspschur,viewer);CHKERRQ(ierr); 217 } else { 218 ierr = PetscViewerASCIIPrintf(viewer," not yet available\n");CHKERRQ(ierr); 219 } 220 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 221 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 222 } else if (isdraw) { 223 PetscDraw draw; 224 PetscReal x,y,w,wd,h; 225 PetscInt cnt = 2; 226 char str[32]; 227 228 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 229 ierr = PetscDrawGetCurrentPoint(draw,&x,&y);CHKERRQ(ierr); 230 if (jac->kspupper != jac->head->ksp) cnt++; 231 w = 2*PetscMin(1.0 - x,x); 232 wd = w/(cnt + 1); 233 234 ierr = PetscSNPrintf(str,32,"Schur fact. %s",PCFieldSplitSchurFactTypes[jac->schurfactorization]);CHKERRQ(ierr); 235 ierr = PetscDrawBoxedString(draw,x,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);CHKERRQ(ierr); 236 y -= h; 237 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_USER && !jac->schur_user) { 238 ierr = PetscSNPrintf(str,32,"Prec. for Schur from %s",PCFieldSplitSchurPreTypes[PC_FIELDSPLIT_SCHUR_PRE_A11]);CHKERRQ(ierr); 239 } else { 240 ierr = PetscSNPrintf(str,32,"Prec. for Schur from %s",PCFieldSplitSchurPreTypes[jac->schurpre]);CHKERRQ(ierr); 241 } 242 ierr = PetscDrawBoxedString(draw,x+wd*(cnt-1)/2.0,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);CHKERRQ(ierr); 243 y -= h; 244 x = x - wd*(cnt-1)/2.0; 245 246 ierr = PetscDrawPushCurrentPoint(draw,x,y);CHKERRQ(ierr); 247 ierr = KSPView(jac->head->ksp,viewer);CHKERRQ(ierr); 248 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr); 249 if (jac->kspupper != jac->head->ksp) { 250 x += wd; 251 ierr = PetscDrawPushCurrentPoint(draw,x,y);CHKERRQ(ierr); 252 ierr = KSPView(jac->kspupper,viewer);CHKERRQ(ierr); 253 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr); 254 } 255 x += wd; 256 ierr = PetscDrawPushCurrentPoint(draw,x,y);CHKERRQ(ierr); 257 ierr = KSPView(jac->kspschur,viewer);CHKERRQ(ierr); 258 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr); 259 } 260 PetscFunctionReturn(0); 261 } 262 263 #undef __FUNCT__ 264 #define __FUNCT__ "PCFieldSplitSetRuntimeSplits_Private" 265 /* Precondition: jac->bs is set to a meaningful value */ 266 static PetscErrorCode PCFieldSplitSetRuntimeSplits_Private(PC pc) 267 { 268 PetscErrorCode ierr; 269 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 270 PetscInt i,nfields,*ifields,nfields_col,*ifields_col; 271 PetscBool flg,flg_col; 272 char optionname[128],splitname[8],optionname_col[128]; 273 274 PetscFunctionBegin; 275 ierr = PetscMalloc1(jac->bs,&ifields);CHKERRQ(ierr); 276 ierr = PetscMalloc1(jac->bs,&ifields_col);CHKERRQ(ierr); 277 for (i=0,flg=PETSC_TRUE;; i++) { 278 ierr = PetscSNPrintf(splitname,sizeof(splitname),"%D",i);CHKERRQ(ierr); 279 ierr = PetscSNPrintf(optionname,sizeof(optionname),"-pc_fieldsplit_%D_fields",i);CHKERRQ(ierr); 280 ierr = PetscSNPrintf(optionname_col,sizeof(optionname_col),"-pc_fieldsplit_%D_fields_col",i);CHKERRQ(ierr); 281 nfields = jac->bs; 282 nfields_col = jac->bs; 283 ierr = PetscOptionsGetIntArray(((PetscObject)pc)->prefix,optionname,ifields,&nfields,&flg);CHKERRQ(ierr); 284 ierr = PetscOptionsGetIntArray(((PetscObject)pc)->prefix,optionname_col,ifields_col,&nfields_col,&flg_col);CHKERRQ(ierr); 285 if (!flg) break; 286 else if (flg && !flg_col) { 287 if (!nfields) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot list zero fields"); 288 ierr = PCFieldSplitSetFields(pc,splitname,nfields,ifields,ifields);CHKERRQ(ierr); 289 } else { 290 if (!nfields || !nfields_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot list zero fields"); 291 if (nfields != nfields_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Number of row and column fields must match"); 292 ierr = PCFieldSplitSetFields(pc,splitname,nfields,ifields,ifields_col);CHKERRQ(ierr); 293 } 294 } 295 if (i > 0) { 296 /* Makes command-line setting of splits take precedence over setting them in code. 297 Otherwise subsequent calls to PCFieldSplitSetIS() or PCFieldSplitSetFields() would 298 create new splits, which would probably not be what the user wanted. */ 299 jac->splitdefined = PETSC_TRUE; 300 } 301 ierr = PetscFree(ifields);CHKERRQ(ierr); 302 ierr = PetscFree(ifields_col);CHKERRQ(ierr); 303 PetscFunctionReturn(0); 304 } 305 306 #undef __FUNCT__ 307 #define __FUNCT__ "PCFieldSplitSetDefaults" 308 static PetscErrorCode PCFieldSplitSetDefaults(PC pc) 309 { 310 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 311 PetscErrorCode ierr; 312 PC_FieldSplitLink ilink = jac->head; 313 PetscBool fieldsplit_default = PETSC_FALSE,stokes = PETSC_FALSE; 314 PetscInt i; 315 316 PetscFunctionBegin; 317 /* 318 Kinda messy, but at least this now uses DMCreateFieldDecomposition() even with jac->reset. 319 Should probably be rewritten. 320 */ 321 if (!ilink || jac->reset) { 322 ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,NULL);CHKERRQ(ierr); 323 if (pc->dm && jac->dm_splits && !stokes) { 324 PetscInt numFields, f, i, j; 325 char **fieldNames; 326 IS *fields; 327 DM *dms; 328 DM subdm[128]; 329 PetscBool flg; 330 331 ierr = DMCreateFieldDecomposition(pc->dm, &numFields, &fieldNames, &fields, &dms);CHKERRQ(ierr); 332 /* Allow the user to prescribe the splits */ 333 for (i = 0, flg = PETSC_TRUE;; i++) { 334 PetscInt ifields[128]; 335 IS compField; 336 char optionname[128], splitname[8]; 337 PetscInt nfields = numFields; 338 339 ierr = PetscSNPrintf(optionname, sizeof(optionname), "-pc_fieldsplit_%D_fields", i);CHKERRQ(ierr); 340 ierr = PetscOptionsGetIntArray(((PetscObject) pc)->prefix, optionname, ifields, &nfields, &flg);CHKERRQ(ierr); 341 if (!flg) break; 342 if (numFields > 128) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot currently support %d > 128 fields", numFields); 343 ierr = DMCreateSubDM(pc->dm, nfields, ifields, &compField, &subdm[i]);CHKERRQ(ierr); 344 if (nfields == 1) { 345 ierr = PCFieldSplitSetIS(pc, fieldNames[ifields[0]], compField);CHKERRQ(ierr); 346 /* ierr = PetscPrintf(PetscObjectComm((PetscObject)pc), "%s Field Indices:", fieldNames[ifields[0]]);CHKERRQ(ierr); 347 ierr = ISView(compField, NULL);CHKERRQ(ierr); */ 348 } else { 349 ierr = PetscSNPrintf(splitname, sizeof(splitname), "%D", i);CHKERRQ(ierr); 350 ierr = PCFieldSplitSetIS(pc, splitname, compField);CHKERRQ(ierr); 351 /* ierr = PetscPrintf(PetscObjectComm((PetscObject)pc), "%s Field Indices:", splitname);CHKERRQ(ierr); 352 ierr = ISView(compField, NULL);CHKERRQ(ierr); */ 353 } 354 ierr = ISDestroy(&compField);CHKERRQ(ierr); 355 for (j = 0; j < nfields; ++j) { 356 f = ifields[j]; 357 ierr = PetscFree(fieldNames[f]);CHKERRQ(ierr); 358 ierr = ISDestroy(&fields[f]);CHKERRQ(ierr); 359 } 360 } 361 if (i == 0) { 362 for (f = 0; f < numFields; ++f) { 363 ierr = PCFieldSplitSetIS(pc, fieldNames[f], fields[f]);CHKERRQ(ierr); 364 ierr = PetscFree(fieldNames[f]);CHKERRQ(ierr); 365 ierr = ISDestroy(&fields[f]);CHKERRQ(ierr); 366 } 367 } else { 368 for (j=0; j<numFields; j++) { 369 ierr = DMDestroy(dms+j);CHKERRQ(ierr); 370 } 371 ierr = PetscFree(dms);CHKERRQ(ierr); 372 ierr = PetscMalloc1(i, &dms);CHKERRQ(ierr); 373 for (j = 0; j < i; ++j) dms[j] = subdm[j]; 374 } 375 ierr = PetscFree(fieldNames);CHKERRQ(ierr); 376 ierr = PetscFree(fields);CHKERRQ(ierr); 377 if (dms) { 378 ierr = PetscInfo(pc, "Setting up physics based fieldsplit preconditioner using the embedded DM\n");CHKERRQ(ierr); 379 for (ilink = jac->head, i = 0; ilink; ilink = ilink->next, ++i) { 380 const char *prefix; 381 ierr = PetscObjectGetOptionsPrefix((PetscObject)(ilink->ksp),&prefix);CHKERRQ(ierr); 382 ierr = PetscObjectSetOptionsPrefix((PetscObject)(dms[i]), prefix);CHKERRQ(ierr); 383 ierr = KSPSetDM(ilink->ksp, dms[i]);CHKERRQ(ierr); 384 ierr = KSPSetDMActive(ilink->ksp, PETSC_FALSE);CHKERRQ(ierr); 385 ierr = PetscObjectIncrementTabLevel((PetscObject)dms[i],(PetscObject)ilink->ksp,0);CHKERRQ(ierr); 386 ierr = DMDestroy(&dms[i]);CHKERRQ(ierr); 387 } 388 ierr = PetscFree(dms);CHKERRQ(ierr); 389 } 390 } else { 391 if (jac->bs <= 0) { 392 if (pc->pmat) { 393 ierr = MatGetBlockSize(pc->pmat,&jac->bs);CHKERRQ(ierr); 394 } else jac->bs = 1; 395 } 396 397 if (stokes) { 398 IS zerodiags,rest; 399 PetscInt nmin,nmax; 400 401 ierr = MatGetOwnershipRange(pc->mat,&nmin,&nmax);CHKERRQ(ierr); 402 ierr = MatFindZeroDiagonals(pc->mat,&zerodiags);CHKERRQ(ierr); 403 ierr = ISComplement(zerodiags,nmin,nmax,&rest);CHKERRQ(ierr); 404 if (jac->reset) { 405 jac->head->is = rest; 406 jac->head->next->is = zerodiags; 407 } else { 408 ierr = PCFieldSplitSetIS(pc,"0",rest);CHKERRQ(ierr); 409 ierr = PCFieldSplitSetIS(pc,"1",zerodiags);CHKERRQ(ierr); 410 } 411 ierr = ISDestroy(&zerodiags);CHKERRQ(ierr); 412 ierr = ISDestroy(&rest);CHKERRQ(ierr); 413 } else { 414 if (jac->reset) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cases not yet handled when PCReset() was used"); 415 ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_default",&fieldsplit_default,NULL);CHKERRQ(ierr); 416 if (!fieldsplit_default) { 417 /* Allow user to set fields from command line, if bs was known at the time of PCSetFromOptions_FieldSplit() 418 then it is set there. This is not ideal because we should only have options set in XXSetFromOptions(). */ 419 ierr = PCFieldSplitSetRuntimeSplits_Private(pc);CHKERRQ(ierr); 420 if (jac->splitdefined) {ierr = PetscInfo(pc,"Splits defined using the options database\n");CHKERRQ(ierr);} 421 } 422 if (fieldsplit_default || !jac->splitdefined) { 423 ierr = PetscInfo(pc,"Using default splitting of fields\n");CHKERRQ(ierr); 424 for (i=0; i<jac->bs; i++) { 425 char splitname[8]; 426 ierr = PetscSNPrintf(splitname,sizeof(splitname),"%D",i);CHKERRQ(ierr); 427 ierr = PCFieldSplitSetFields(pc,splitname,1,&i,&i);CHKERRQ(ierr); 428 } 429 jac->defaultsplit = PETSC_TRUE; 430 } 431 } 432 } 433 } else if (jac->nsplits == 1) { 434 if (ilink->is) { 435 IS is2; 436 PetscInt nmin,nmax; 437 438 ierr = MatGetOwnershipRange(pc->mat,&nmin,&nmax);CHKERRQ(ierr); 439 ierr = ISComplement(ilink->is,nmin,nmax,&is2);CHKERRQ(ierr); 440 ierr = PCFieldSplitSetIS(pc,"1",is2);CHKERRQ(ierr); 441 ierr = ISDestroy(&is2);CHKERRQ(ierr); 442 } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Must provide at least two sets of fields to PCFieldSplit()"); 443 } 444 445 446 if (jac->nsplits < 2) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Unhandled case, must have at least two fields, not %d", jac->nsplits); 447 PetscFunctionReturn(0); 448 } 449 450 PETSC_EXTERN PetscErrorCode PetscOptionsFindPairPrefix_Private(const char pre[], const char name[], char *value[], PetscBool *flg); 451 452 #undef __FUNCT__ 453 #define __FUNCT__ "PCSetUp_FieldSplit" 454 static PetscErrorCode PCSetUp_FieldSplit(PC pc) 455 { 456 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 457 PetscErrorCode ierr; 458 PC_FieldSplitLink ilink; 459 PetscInt i,nsplit; 460 MatStructure flag = pc->flag; 461 PetscBool sorted, sorted_col; 462 463 PetscFunctionBegin; 464 ierr = PCFieldSplitSetDefaults(pc);CHKERRQ(ierr); 465 nsplit = jac->nsplits; 466 ilink = jac->head; 467 468 /* get the matrices for each split */ 469 if (!jac->issetup) { 470 PetscInt rstart,rend,nslots,bs; 471 472 jac->issetup = PETSC_TRUE; 473 474 /* This is done here instead of in PCFieldSplitSetFields() because may not have matrix at that point */ 475 if (jac->defaultsplit || !ilink->is) { 476 if (jac->bs <= 0) jac->bs = nsplit; 477 } 478 bs = jac->bs; 479 ierr = MatGetOwnershipRange(pc->pmat,&rstart,&rend);CHKERRQ(ierr); 480 nslots = (rend - rstart)/bs; 481 for (i=0; i<nsplit; i++) { 482 if (jac->defaultsplit) { 483 ierr = ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+i,nsplit,&ilink->is);CHKERRQ(ierr); 484 ierr = ISDuplicate(ilink->is,&ilink->is_col);CHKERRQ(ierr); 485 } else if (!ilink->is) { 486 if (ilink->nfields > 1) { 487 PetscInt *ii,*jj,j,k,nfields = ilink->nfields,*fields = ilink->fields,*fields_col = ilink->fields_col; 488 ierr = PetscMalloc1(ilink->nfields*nslots,&ii);CHKERRQ(ierr); 489 ierr = PetscMalloc1(ilink->nfields*nslots,&jj);CHKERRQ(ierr); 490 for (j=0; j<nslots; j++) { 491 for (k=0; k<nfields; k++) { 492 ii[nfields*j + k] = rstart + bs*j + fields[k]; 493 jj[nfields*j + k] = rstart + bs*j + fields_col[k]; 494 } 495 } 496 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,ii,PETSC_OWN_POINTER,&ilink->is);CHKERRQ(ierr); 497 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,jj,PETSC_OWN_POINTER,&ilink->is_col);CHKERRQ(ierr); 498 } else { 499 ierr = ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields[0],bs,&ilink->is);CHKERRQ(ierr); 500 ierr = ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields_col[0],bs,&ilink->is_col);CHKERRQ(ierr); 501 } 502 } 503 ierr = ISSorted(ilink->is,&sorted);CHKERRQ(ierr); 504 if (ilink->is_col) { ierr = ISSorted(ilink->is_col,&sorted_col);CHKERRQ(ierr); } 505 if (!sorted || !sorted_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Fields must be sorted when creating split"); 506 ilink = ilink->next; 507 } 508 } 509 510 ilink = jac->head; 511 if (!jac->pmat) { 512 Vec xtmp; 513 514 ierr = MatGetVecs(pc->pmat,&xtmp,NULL);CHKERRQ(ierr); 515 ierr = PetscMalloc1(nsplit,&jac->pmat);CHKERRQ(ierr); 516 ierr = PetscMalloc2(nsplit,&jac->x,nsplit,&jac->y);CHKERRQ(ierr); 517 for (i=0; i<nsplit; i++) { 518 MatNullSpace sp; 519 520 /* Check for preconditioning matrix attached to IS */ 521 ierr = PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &jac->pmat[i]);CHKERRQ(ierr); 522 if (jac->pmat[i]) { 523 ierr = PetscObjectReference((PetscObject) jac->pmat[i]);CHKERRQ(ierr); 524 if (jac->type == PC_COMPOSITE_SCHUR) { 525 jac->schur_user = jac->pmat[i]; 526 527 ierr = PetscObjectReference((PetscObject) jac->schur_user);CHKERRQ(ierr); 528 } 529 } else { 530 ierr = MatGetSubMatrix(pc->pmat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->pmat[i]);CHKERRQ(ierr); 531 } 532 /* create work vectors for each split */ 533 ierr = MatGetVecs(jac->pmat[i],&jac->x[i],&jac->y[i]);CHKERRQ(ierr); 534 ilink->x = jac->x[i]; ilink->y = jac->y[i]; ilink->z = NULL; 535 /* compute scatter contexts needed by multiplicative versions and non-default splits */ 536 ierr = VecScatterCreate(xtmp,ilink->is,jac->x[i],NULL,&ilink->sctx);CHKERRQ(ierr); 537 /* Check for null space attached to IS */ 538 ierr = PetscObjectQuery((PetscObject) ilink->is, "nullspace", (PetscObject*) &sp);CHKERRQ(ierr); 539 if (sp) { 540 ierr = MatSetNullSpace(jac->pmat[i], sp);CHKERRQ(ierr); 541 } 542 ierr = PetscObjectQuery((PetscObject) ilink->is, "nearnullspace", (PetscObject*) &sp);CHKERRQ(ierr); 543 if (sp) { 544 ierr = MatSetNearNullSpace(jac->pmat[i], sp);CHKERRQ(ierr); 545 } 546 ilink = ilink->next; 547 } 548 ierr = VecDestroy(&xtmp);CHKERRQ(ierr); 549 } else { 550 for (i=0; i<nsplit; i++) { 551 Mat pmat; 552 553 /* Check for preconditioning matrix attached to IS */ 554 ierr = PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &pmat);CHKERRQ(ierr); 555 if (!pmat) { 556 ierr = MatGetSubMatrix(pc->pmat,ilink->is,ilink->is_col,MAT_REUSE_MATRIX,&jac->pmat[i]);CHKERRQ(ierr); 557 } 558 ilink = ilink->next; 559 } 560 } 561 if (pc->useAmat) { 562 ilink = jac->head; 563 if (!jac->mat) { 564 ierr = PetscMalloc1(nsplit,&jac->mat);CHKERRQ(ierr); 565 for (i=0; i<nsplit; i++) { 566 ierr = MatGetSubMatrix(pc->mat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->mat[i]);CHKERRQ(ierr); 567 ilink = ilink->next; 568 } 569 } else { 570 for (i=0; i<nsplit; i++) { 571 if (jac->mat[i]) {ierr = MatGetSubMatrix(pc->mat,ilink->is,ilink->is_col,MAT_REUSE_MATRIX,&jac->mat[i]);CHKERRQ(ierr);} 572 ilink = ilink->next; 573 } 574 } 575 } else { 576 jac->mat = jac->pmat; 577 } 578 579 if (jac->type != PC_COMPOSITE_ADDITIVE && jac->type != PC_COMPOSITE_SCHUR) { 580 /* extract the rows of the matrix associated with each field: used for efficient computation of residual inside algorithm */ 581 ilink = jac->head; 582 if (!jac->Afield) { 583 ierr = PetscMalloc1(nsplit,&jac->Afield);CHKERRQ(ierr); 584 for (i=0; i<nsplit; i++) { 585 ierr = MatGetSubMatrix(pc->mat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);CHKERRQ(ierr); 586 ilink = ilink->next; 587 } 588 } else { 589 for (i=0; i<nsplit; i++) { 590 ierr = MatGetSubMatrix(pc->mat,ilink->is,NULL,MAT_REUSE_MATRIX,&jac->Afield[i]);CHKERRQ(ierr); 591 ilink = ilink->next; 592 } 593 } 594 } 595 596 if (jac->type == PC_COMPOSITE_SCHUR) { 597 IS ccis; 598 PetscInt rstart,rend; 599 char lscname[256]; 600 PetscObject LSC_L; 601 602 if (nsplit != 2) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"To use Schur complement preconditioner you must have exactly 2 fields"); 603 604 /* When extracting off-diagonal submatrices, we take complements from this range */ 605 ierr = MatGetOwnershipRangeColumn(pc->mat,&rstart,&rend);CHKERRQ(ierr); 606 607 /* need to handle case when one is resetting up the preconditioner */ 608 if (jac->schur) { 609 KSP kspA = jac->head->ksp, kspInner = NULL, kspUpper = jac->kspupper; 610 611 ierr = MatSchurComplementGetKSP(jac->schur, &kspInner);CHKERRQ(ierr); 612 ilink = jac->head; 613 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 614 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->B);CHKERRQ(ierr); 615 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 616 ilink = ilink->next; 617 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 618 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->C);CHKERRQ(ierr); 619 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 620 ierr = MatSchurComplementUpdate(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1],pc->flag);CHKERRQ(ierr); 621 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) { 622 ierr = MatDestroy(&jac->schurp);CHKERRQ(ierr); 623 ierr = MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);CHKERRQ(ierr); 624 } 625 if (kspA != kspInner) { 626 ierr = KSPSetOperators(kspA,jac->mat[0],jac->pmat[0],pc->flag);CHKERRQ(ierr); 627 } 628 if (kspUpper != kspA) { 629 ierr = KSPSetOperators(kspUpper,jac->mat[0],jac->pmat[0],pc->flag);CHKERRQ(ierr); 630 } 631 ierr = KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac),pc->flag);CHKERRQ(ierr); 632 } else { 633 const char *Dprefix; 634 char schurprefix[256]; 635 char schurtestoption[256]; 636 MatNullSpace sp; 637 PetscBool flg; 638 639 /* extract the A01 and A10 matrices */ 640 ilink = jac->head; 641 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 642 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);CHKERRQ(ierr); 643 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 644 ilink = ilink->next; 645 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 646 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);CHKERRQ(ierr); 647 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 648 649 /* Use mat[0] (diagonal block of Amat) preconditioned by pmat[0] to define Schur complement */ 650 ierr = MatCreate(((PetscObject)jac->mat[0])->comm,&jac->schur);CHKERRQ(ierr); 651 ierr = MatSetType(jac->schur,MATSCHURCOMPLEMENT);CHKERRQ(ierr); 652 ierr = MatSchurComplementSet(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);CHKERRQ(ierr); 653 654 ierr = MatGetNullSpace(jac->pmat[1], &sp);CHKERRQ(ierr); 655 if (sp) { 656 ierr = MatSetNullSpace(jac->schur, sp);CHKERRQ(ierr); 657 } 658 659 ierr = PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_inner_", ilink->splitname);CHKERRQ(ierr); 660 ierr = PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);CHKERRQ(ierr); 661 if (flg) { 662 DM dmInner; 663 KSP kspInner; 664 665 ierr = MatSchurComplementGetKSP(jac->schur, &kspInner);CHKERRQ(ierr); 666 ierr = PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_inner_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);CHKERRQ(ierr); 667 /* Indent this deeper to emphasize the "inner" nature of this solver. */ 668 ierr = PetscObjectIncrementTabLevel((PetscObject)kspInner, (PetscObject) pc, 2);CHKERRQ(ierr); 669 ierr = KSPSetOptionsPrefix(kspInner, schurprefix);CHKERRQ(ierr); 670 671 /* Set DM for new solver */ 672 ierr = KSPGetDM(jac->head->ksp, &dmInner);CHKERRQ(ierr); 673 ierr = KSPSetDM(kspInner, dmInner);CHKERRQ(ierr); 674 ierr = KSPSetDMActive(kspInner, PETSC_FALSE);CHKERRQ(ierr); 675 } else { 676 /* Use the outer solver for the inner solve, but revert the KSPPREONLY from PCFieldSplitSetFields_FieldSplit or 677 * PCFieldSplitSetIS_FieldSplit. We don't want KSPPREONLY because it makes the Schur complement inexact, 678 * preventing Schur complement reduction to be an accurate solve. Usually when an iterative solver is used for 679 * S = D - C A_inner^{-1} B, we expect S to be defined using an accurate definition of A_inner^{-1}, so we make 680 * GMRES the default. Note that it is also common to use PREONLY for S, in which case S may not be used 681 * directly, and the user is responsible for setting an inexact method for fieldsplit's A^{-1}. */ 682 ierr = KSPSetType(jac->head->ksp,KSPGMRES);CHKERRQ(ierr); 683 ierr = MatSchurComplementSetKSP(jac->schur,jac->head->ksp);CHKERRQ(ierr); 684 } 685 ierr = KSPSetOperators(jac->head->ksp,jac->mat[0],jac->pmat[0],flag);CHKERRQ(ierr); 686 ierr = KSPSetFromOptions(jac->head->ksp);CHKERRQ(ierr); 687 ierr = MatSetFromOptions(jac->schur);CHKERRQ(ierr); 688 689 ierr = PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_upper_", ilink->splitname);CHKERRQ(ierr); 690 ierr = PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);CHKERRQ(ierr); 691 if (flg) { 692 DM dmInner; 693 694 ierr = PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_upper_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);CHKERRQ(ierr); 695 ierr = KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspupper);CHKERRQ(ierr); 696 ierr = KSPSetOptionsPrefix(jac->kspupper, schurprefix);CHKERRQ(ierr); 697 ierr = KSPGetDM(jac->head->ksp, &dmInner);CHKERRQ(ierr); 698 ierr = KSPSetDM(jac->kspupper, dmInner);CHKERRQ(ierr); 699 ierr = KSPSetDMActive(jac->kspupper, PETSC_FALSE);CHKERRQ(ierr); 700 ierr = KSPSetFromOptions(jac->kspupper);CHKERRQ(ierr); 701 ierr = KSPSetOperators(jac->kspupper,jac->mat[0],jac->pmat[0],flag);CHKERRQ(ierr); 702 ierr = VecDuplicate(jac->head->x, &jac->head->z);CHKERRQ(ierr); 703 } else { 704 jac->kspupper = jac->head->ksp; 705 ierr = PetscObjectReference((PetscObject) jac->head->ksp);CHKERRQ(ierr); 706 } 707 708 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) { 709 ierr = MatSchurComplementSetPmatLumping(jac->schur,jac->schurp_lump);CHKERRQ(ierr); 710 ierr = MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);CHKERRQ(ierr); 711 } 712 ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&jac->kspschur);CHKERRQ(ierr); 713 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->kspschur);CHKERRQ(ierr); 714 ierr = PetscObjectIncrementTabLevel((PetscObject)jac->kspschur,(PetscObject)pc,1);CHKERRQ(ierr); 715 ierr = KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac),DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 716 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) { 717 PC pcschur; 718 ierr = KSPGetPC(jac->kspschur,&pcschur);CHKERRQ(ierr); 719 ierr = PCSetType(pcschur,PCNONE);CHKERRQ(ierr); 720 /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */ 721 } 722 ierr = KSPGetOptionsPrefix(jac->head->next->ksp, &Dprefix);CHKERRQ(ierr); 723 ierr = KSPSetOptionsPrefix(jac->kspschur, Dprefix);CHKERRQ(ierr); 724 /* propogate DM */ 725 { 726 DM sdm; 727 ierr = KSPGetDM(jac->head->next->ksp, &sdm);CHKERRQ(ierr); 728 if (sdm) { 729 ierr = KSPSetDM(jac->kspschur, sdm);CHKERRQ(ierr); 730 ierr = KSPSetDMActive(jac->kspschur, PETSC_FALSE);CHKERRQ(ierr); 731 } 732 } 733 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 734 /* need to call this every time, since the jac->kspschur is freshly created, otherwise its options never get set */ 735 ierr = KSPSetFromOptions(jac->kspschur);CHKERRQ(ierr); 736 } 737 738 /* HACK: special support to forward L and Lp matrices that might be used by PCLSC */ 739 ierr = PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_L",ilink->splitname);CHKERRQ(ierr); 740 ierr = PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr); 741 if (!LSC_L) {ierr = PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr);} 742 if (LSC_L) {ierr = PetscObjectCompose((PetscObject)jac->schur,"LSC_L",(PetscObject)LSC_L);CHKERRQ(ierr);} 743 ierr = PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_Lp",ilink->splitname);CHKERRQ(ierr); 744 ierr = PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr); 745 if (!LSC_L) {ierr = PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr);} 746 if (LSC_L) {ierr = PetscObjectCompose((PetscObject)jac->schur,"LSC_Lp",(PetscObject)LSC_L);CHKERRQ(ierr);} 747 } else { 748 /* set up the individual splits' PCs */ 749 i = 0; 750 ilink = jac->head; 751 while (ilink) { 752 ierr = KSPSetOperators(ilink->ksp,jac->mat[i],jac->pmat[i],flag);CHKERRQ(ierr); 753 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 754 if (!jac->suboptionsset) {ierr = KSPSetFromOptions(ilink->ksp);CHKERRQ(ierr);} 755 i++; 756 ilink = ilink->next; 757 } 758 } 759 760 jac->suboptionsset = PETSC_TRUE; 761 PetscFunctionReturn(0); 762 } 763 764 #define FieldSplitSplitSolveAdd(ilink,xx,yy) \ 765 (VecScatterBegin(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 766 VecScatterEnd(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 767 KSPSolve(ilink->ksp,ilink->x,ilink->y) || \ 768 VecScatterBegin(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE) || \ 769 VecScatterEnd(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE)) 770 771 #undef __FUNCT__ 772 #define __FUNCT__ "PCApply_FieldSplit_Schur" 773 static PetscErrorCode PCApply_FieldSplit_Schur(PC pc,Vec x,Vec y) 774 { 775 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 776 PetscErrorCode ierr; 777 PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next; 778 KSP kspA = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper; 779 780 PetscFunctionBegin; 781 switch (jac->schurfactorization) { 782 case PC_FIELDSPLIT_SCHUR_FACT_DIAG: 783 /* [A00 0; 0 -S], positive definite, suitable for MINRES */ 784 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 785 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 786 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 787 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 788 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 789 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 790 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 791 ierr = VecScale(ilinkD->y,-1.);CHKERRQ(ierr); 792 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 793 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 794 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 795 break; 796 case PC_FIELDSPLIT_SCHUR_FACT_LOWER: 797 /* [A00 0; A10 S], suitable for left preconditioning */ 798 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 799 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 800 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 801 ierr = MatMult(jac->C,ilinkA->y,ilinkD->x);CHKERRQ(ierr); 802 ierr = VecScale(ilinkD->x,-1.);CHKERRQ(ierr); 803 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 804 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 805 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 806 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 807 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 808 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 809 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 810 break; 811 case PC_FIELDSPLIT_SCHUR_FACT_UPPER: 812 /* [A00 A01; 0 S], suitable for right preconditioning */ 813 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 814 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 815 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 816 ierr = MatMult(jac->B,ilinkD->y,ilinkA->x);CHKERRQ(ierr); 817 ierr = VecScale(ilinkA->x,-1.);CHKERRQ(ierr); 818 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 819 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 820 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 821 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 822 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 823 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 824 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 825 break; 826 case PC_FIELDSPLIT_SCHUR_FACT_FULL: 827 /* [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 */ 828 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 829 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 830 ierr = KSPSolve(kspLower,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 831 ierr = MatMult(jac->C,ilinkA->y,ilinkD->x);CHKERRQ(ierr); 832 ierr = VecScale(ilinkD->x,-1.0);CHKERRQ(ierr); 833 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 834 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 835 836 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 837 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 838 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 839 840 if (kspUpper == kspA) { 841 ierr = MatMult(jac->B,ilinkD->y,ilinkA->y);CHKERRQ(ierr); 842 ierr = VecAXPY(ilinkA->x,-1.0,ilinkA->y);CHKERRQ(ierr); 843 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 844 } else { 845 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 846 ierr = MatMult(jac->B,ilinkD->y,ilinkA->x);CHKERRQ(ierr); 847 ierr = KSPSolve(kspUpper,ilinkA->x,ilinkA->z);CHKERRQ(ierr); 848 ierr = VecAXPY(ilinkA->y,-1.0,ilinkA->z);CHKERRQ(ierr); 849 } 850 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 851 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 852 } 853 PetscFunctionReturn(0); 854 } 855 856 #undef __FUNCT__ 857 #define __FUNCT__ "PCApply_FieldSplit" 858 static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y) 859 { 860 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 861 PetscErrorCode ierr; 862 PC_FieldSplitLink ilink = jac->head; 863 PetscInt cnt,bs; 864 865 PetscFunctionBegin; 866 if (jac->type == PC_COMPOSITE_ADDITIVE) { 867 if (jac->defaultsplit) { 868 ierr = VecGetBlockSize(x,&bs);CHKERRQ(ierr); 869 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); 870 ierr = VecGetBlockSize(y,&bs);CHKERRQ(ierr); 871 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); 872 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 873 while (ilink) { 874 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 875 ilink = ilink->next; 876 } 877 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 878 } else { 879 ierr = VecSet(y,0.0);CHKERRQ(ierr); 880 while (ilink) { 881 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 882 ilink = ilink->next; 883 } 884 } 885 } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 886 if (!jac->w1) { 887 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 888 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 889 } 890 ierr = VecSet(y,0.0);CHKERRQ(ierr); 891 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 892 cnt = 1; 893 while (ilink->next) { 894 ilink = ilink->next; 895 /* compute the residual only over the part of the vector needed */ 896 ierr = MatMult(jac->Afield[cnt++],y,ilink->x);CHKERRQ(ierr); 897 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 898 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 899 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 900 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 901 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 902 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 903 } 904 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 905 cnt -= 2; 906 while (ilink->previous) { 907 ilink = ilink->previous; 908 /* compute the residual only over the part of the vector needed */ 909 ierr = MatMult(jac->Afield[cnt--],y,ilink->x);CHKERRQ(ierr); 910 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 911 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 912 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 913 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 914 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 915 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 916 } 917 } 918 } else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Unsupported or unknown composition",(int) jac->type); 919 PetscFunctionReturn(0); 920 } 921 922 #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \ 923 (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 924 VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 925 KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) || \ 926 VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \ 927 VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE)) 928 929 #undef __FUNCT__ 930 #define __FUNCT__ "PCApplyTranspose_FieldSplit" 931 static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y) 932 { 933 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 934 PetscErrorCode ierr; 935 PC_FieldSplitLink ilink = jac->head; 936 PetscInt bs; 937 938 PetscFunctionBegin; 939 if (jac->type == PC_COMPOSITE_ADDITIVE) { 940 if (jac->defaultsplit) { 941 ierr = VecGetBlockSize(x,&bs);CHKERRQ(ierr); 942 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); 943 ierr = VecGetBlockSize(y,&bs);CHKERRQ(ierr); 944 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); 945 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 946 while (ilink) { 947 ierr = KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 948 ilink = ilink->next; 949 } 950 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 951 } else { 952 ierr = VecSet(y,0.0);CHKERRQ(ierr); 953 while (ilink) { 954 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 955 ilink = ilink->next; 956 } 957 } 958 } else { 959 if (!jac->w1) { 960 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 961 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 962 } 963 ierr = VecSet(y,0.0);CHKERRQ(ierr); 964 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 965 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 966 while (ilink->next) { 967 ilink = ilink->next; 968 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 969 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 970 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 971 } 972 while (ilink->previous) { 973 ilink = ilink->previous; 974 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 975 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 976 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 977 } 978 } else { 979 while (ilink->next) { /* get to last entry in linked list */ 980 ilink = ilink->next; 981 } 982 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 983 while (ilink->previous) { 984 ilink = ilink->previous; 985 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 986 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 987 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 988 } 989 } 990 } 991 PetscFunctionReturn(0); 992 } 993 994 #undef __FUNCT__ 995 #define __FUNCT__ "PCReset_FieldSplit" 996 static PetscErrorCode PCReset_FieldSplit(PC pc) 997 { 998 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 999 PetscErrorCode ierr; 1000 PC_FieldSplitLink ilink = jac->head,next; 1001 1002 PetscFunctionBegin; 1003 while (ilink) { 1004 ierr = KSPReset(ilink->ksp);CHKERRQ(ierr); 1005 ierr = VecDestroy(&ilink->x);CHKERRQ(ierr); 1006 ierr = VecDestroy(&ilink->y);CHKERRQ(ierr); 1007 ierr = VecDestroy(&ilink->z);CHKERRQ(ierr); 1008 ierr = VecScatterDestroy(&ilink->sctx);CHKERRQ(ierr); 1009 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 1010 ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr); 1011 next = ilink->next; 1012 ilink = next; 1013 } 1014 ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr); 1015 if (jac->mat && jac->mat != jac->pmat) { 1016 ierr = MatDestroyMatrices(jac->nsplits,&jac->mat);CHKERRQ(ierr); 1017 } else if (jac->mat) { 1018 jac->mat = NULL; 1019 } 1020 if (jac->pmat) {ierr = MatDestroyMatrices(jac->nsplits,&jac->pmat);CHKERRQ(ierr);} 1021 if (jac->Afield) {ierr = MatDestroyMatrices(jac->nsplits,&jac->Afield);CHKERRQ(ierr);} 1022 ierr = VecDestroy(&jac->w1);CHKERRQ(ierr); 1023 ierr = VecDestroy(&jac->w2);CHKERRQ(ierr); 1024 ierr = MatDestroy(&jac->schur);CHKERRQ(ierr); 1025 ierr = MatDestroy(&jac->schurp);CHKERRQ(ierr); 1026 ierr = MatDestroy(&jac->schur_user);CHKERRQ(ierr); 1027 ierr = KSPDestroy(&jac->kspschur);CHKERRQ(ierr); 1028 ierr = KSPDestroy(&jac->kspupper);CHKERRQ(ierr); 1029 ierr = MatDestroy(&jac->B);CHKERRQ(ierr); 1030 ierr = MatDestroy(&jac->C);CHKERRQ(ierr); 1031 jac->reset = PETSC_TRUE; 1032 PetscFunctionReturn(0); 1033 } 1034 1035 #undef __FUNCT__ 1036 #define __FUNCT__ "PCDestroy_FieldSplit" 1037 static PetscErrorCode PCDestroy_FieldSplit(PC pc) 1038 { 1039 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1040 PetscErrorCode ierr; 1041 PC_FieldSplitLink ilink = jac->head,next; 1042 1043 PetscFunctionBegin; 1044 ierr = PCReset_FieldSplit(pc);CHKERRQ(ierr); 1045 while (ilink) { 1046 ierr = KSPDestroy(&ilink->ksp);CHKERRQ(ierr); 1047 next = ilink->next; 1048 ierr = PetscFree(ilink->splitname);CHKERRQ(ierr); 1049 ierr = PetscFree(ilink->fields);CHKERRQ(ierr); 1050 ierr = PetscFree(ilink->fields_col);CHKERRQ(ierr); 1051 ierr = PetscFree(ilink);CHKERRQ(ierr); 1052 ilink = next; 1053 } 1054 ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr); 1055 ierr = PetscFree(pc->data);CHKERRQ(ierr); 1056 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",NULL);CHKERRQ(ierr); 1057 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",NULL);CHKERRQ(ierr); 1058 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",NULL);CHKERRQ(ierr); 1059 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",NULL);CHKERRQ(ierr); 1060 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",NULL);CHKERRQ(ierr); 1061 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSchurPrecondition_C",NULL);CHKERRQ(ierr); 1062 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",NULL);CHKERRQ(ierr); 1063 PetscFunctionReturn(0); 1064 } 1065 1066 #undef __FUNCT__ 1067 #define __FUNCT__ "PCSetFromOptions_FieldSplit" 1068 static PetscErrorCode PCSetFromOptions_FieldSplit(PC pc) 1069 { 1070 PetscErrorCode ierr; 1071 PetscInt bs; 1072 PetscBool flg,stokes = PETSC_FALSE; 1073 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1074 PCCompositeType ctype; 1075 1076 PetscFunctionBegin; 1077 ierr = PetscOptionsHead("FieldSplit options");CHKERRQ(ierr); 1078 ierr = PetscOptionsBool("-pc_fieldsplit_dm_splits","Whether to use DMCreateFieldDecomposition() for splits","PCFieldSplitSetDMSplits",jac->dm_splits,&jac->dm_splits,NULL);CHKERRQ(ierr); 1079 ierr = PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);CHKERRQ(ierr); 1080 if (flg) { 1081 ierr = PCFieldSplitSetBlockSize(pc,bs);CHKERRQ(ierr); 1082 } 1083 1084 ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,NULL);CHKERRQ(ierr); 1085 if (stokes) { 1086 ierr = PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);CHKERRQ(ierr); 1087 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_SELF; 1088 } 1089 1090 ierr = PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg);CHKERRQ(ierr); 1091 if (flg) { 1092 ierr = PCFieldSplitSetType(pc,ctype);CHKERRQ(ierr); 1093 } 1094 /* Only setup fields once */ 1095 if ((jac->bs > 0) && (jac->nsplits == 0)) { 1096 /* only allow user to set fields from command line if bs is already known. 1097 otherwise user can set them in PCFieldSplitSetDefaults() */ 1098 ierr = PCFieldSplitSetRuntimeSplits_Private(pc);CHKERRQ(ierr); 1099 if (jac->splitdefined) {ierr = PetscInfo(pc,"Splits defined using the options database\n");CHKERRQ(ierr);} 1100 } 1101 if (jac->type == PC_COMPOSITE_SCHUR) { 1102 ierr = PetscOptionsGetEnum(((PetscObject)pc)->prefix,"-pc_fieldsplit_schur_factorization_type",PCFieldSplitSchurFactTypes,(PetscEnum*)&jac->schurfactorization,&flg);CHKERRQ(ierr); 1103 if (flg) {ierr = PetscInfo(pc,"Deprecated use of -pc_fieldsplit_schur_factorization_type\n");CHKERRQ(ierr);} 1104 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); 1105 ierr = PetscOptionsEnum("-pc_fieldsplit_schur_precondition","How to build preconditioner for Schur complement","PCFieldSplitSchurPrecondition",PCFieldSplitSchurPreTypes,(PetscEnum)jac->schurpre,(PetscEnum*)&jac->schurpre,NULL);CHKERRQ(ierr); 1106 jac->schurp_lump = PETSC_FALSE; 1107 ierr = PetscOptionsBool("-pc_fieldsplit_schur_precondition_selfp_lump","Lump A00 when assembling selfp preconditioning matrix for Schur complement","PCFieldSplitSetSelfpLumping",jac->schurp_lump,&jac->schurp_lump,NULL);CHKERRQ(ierr); 1108 } 1109 ierr = PetscOptionsTail();CHKERRQ(ierr); 1110 PetscFunctionReturn(0); 1111 } 1112 1113 /*------------------------------------------------------------------------------------*/ 1114 1115 #undef __FUNCT__ 1116 #define __FUNCT__ "PCFieldSplitSetFields_FieldSplit" 1117 static PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 1118 { 1119 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1120 PetscErrorCode ierr; 1121 PC_FieldSplitLink ilink,next = jac->head; 1122 char prefix[128]; 1123 PetscInt i; 1124 1125 PetscFunctionBegin; 1126 if (jac->splitdefined) { 1127 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 1128 PetscFunctionReturn(0); 1129 } 1130 for (i=0; i<n; i++) { 1131 if (fields[i] >= jac->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs); 1132 if (fields[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]); 1133 } 1134 ierr = PetscNew(&ilink);CHKERRQ(ierr); 1135 if (splitname) { 1136 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 1137 } else { 1138 ierr = PetscMalloc1(3,&ilink->splitname);CHKERRQ(ierr); 1139 ierr = PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);CHKERRQ(ierr); 1140 } 1141 ierr = PetscMalloc1(n,&ilink->fields);CHKERRQ(ierr); 1142 ierr = PetscMemcpy(ilink->fields,fields,n*sizeof(PetscInt));CHKERRQ(ierr); 1143 ierr = PetscMalloc1(n,&ilink->fields_col);CHKERRQ(ierr); 1144 ierr = PetscMemcpy(ilink->fields_col,fields_col,n*sizeof(PetscInt));CHKERRQ(ierr); 1145 1146 ilink->nfields = n; 1147 ilink->next = NULL; 1148 ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);CHKERRQ(ierr); 1149 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 1150 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 1151 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 1152 1153 ierr = PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);CHKERRQ(ierr); 1154 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 1155 1156 if (!next) { 1157 jac->head = ilink; 1158 ilink->previous = NULL; 1159 } else { 1160 while (next->next) { 1161 next = next->next; 1162 } 1163 next->next = ilink; 1164 ilink->previous = next; 1165 } 1166 jac->nsplits++; 1167 PetscFunctionReturn(0); 1168 } 1169 1170 #undef __FUNCT__ 1171 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit_Schur" 1172 static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp) 1173 { 1174 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1175 PetscErrorCode ierr; 1176 1177 PetscFunctionBegin; 1178 ierr = PetscMalloc1(jac->nsplits,subksp);CHKERRQ(ierr); 1179 ierr = MatSchurComplementGetKSP(jac->schur,*subksp);CHKERRQ(ierr); 1180 1181 (*subksp)[1] = jac->kspschur; 1182 if (n) *n = jac->nsplits; 1183 PetscFunctionReturn(0); 1184 } 1185 1186 #undef __FUNCT__ 1187 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit" 1188 static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp) 1189 { 1190 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1191 PetscErrorCode ierr; 1192 PetscInt cnt = 0; 1193 PC_FieldSplitLink ilink = jac->head; 1194 1195 PetscFunctionBegin; 1196 ierr = PetscMalloc1(jac->nsplits,subksp);CHKERRQ(ierr); 1197 while (ilink) { 1198 (*subksp)[cnt++] = ilink->ksp; 1199 ilink = ilink->next; 1200 } 1201 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); 1202 if (n) *n = jac->nsplits; 1203 PetscFunctionReturn(0); 1204 } 1205 1206 #undef __FUNCT__ 1207 #define __FUNCT__ "PCFieldSplitSetIS_FieldSplit" 1208 static PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is) 1209 { 1210 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1211 PetscErrorCode ierr; 1212 PC_FieldSplitLink ilink, next = jac->head; 1213 char prefix[128]; 1214 1215 PetscFunctionBegin; 1216 if (jac->splitdefined) { 1217 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 1218 PetscFunctionReturn(0); 1219 } 1220 ierr = PetscNew(&ilink);CHKERRQ(ierr); 1221 if (splitname) { 1222 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 1223 } else { 1224 ierr = PetscMalloc1(8,&ilink->splitname);CHKERRQ(ierr); 1225 ierr = PetscSNPrintf(ilink->splitname,7,"%D",jac->nsplits);CHKERRQ(ierr); 1226 } 1227 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 1228 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 1229 ilink->is = is; 1230 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 1231 ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr); 1232 ilink->is_col = is; 1233 ilink->next = NULL; 1234 ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);CHKERRQ(ierr); 1235 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 1236 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 1237 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 1238 1239 ierr = PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);CHKERRQ(ierr); 1240 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 1241 1242 if (!next) { 1243 jac->head = ilink; 1244 ilink->previous = NULL; 1245 } else { 1246 while (next->next) { 1247 next = next->next; 1248 } 1249 next->next = ilink; 1250 ilink->previous = next; 1251 } 1252 jac->nsplits++; 1253 PetscFunctionReturn(0); 1254 } 1255 1256 #undef __FUNCT__ 1257 #define __FUNCT__ "PCFieldSplitSetFields" 1258 /*@ 1259 PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner 1260 1261 Logically Collective on PC 1262 1263 Input Parameters: 1264 + pc - the preconditioner context 1265 . splitname - name of this split, if NULL the number of the split is used 1266 . n - the number of fields in this split 1267 - fields - the fields in this split 1268 1269 Level: intermediate 1270 1271 Notes: Use PCFieldSplitSetIS() to set a completely general set of indices as a field. 1272 1273 The PCFieldSplitSetFields() is for defining fields as strided blocks. For example, if the block 1274 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 1275 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x.... 1276 where the numbered entries indicate what is in the field. 1277 1278 This function is called once per split (it creates a new split each time). Solve options 1279 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1280 1281 Developer Note: This routine does not actually create the IS representing the split, that is delayed 1282 until PCSetUp_FieldSplit(), because information about the vector/matrix layouts may not be 1283 available when this routine is called. 1284 1285 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize(), PCFieldSplitSetIS() 1286 1287 @*/ 1288 PetscErrorCode PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 1289 { 1290 PetscErrorCode ierr; 1291 1292 PetscFunctionBegin; 1293 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1294 PetscValidCharPointer(splitname,2); 1295 if (n < 1) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Provided number of fields %D in split \"%s\" not positive",n,splitname); 1296 PetscValidIntPointer(fields,3); 1297 ierr = PetscTryMethod(pc,"PCFieldSplitSetFields_C",(PC,const char[],PetscInt,const PetscInt*,const PetscInt*),(pc,splitname,n,fields,fields_col));CHKERRQ(ierr); 1298 PetscFunctionReturn(0); 1299 } 1300 1301 #undef __FUNCT__ 1302 #define __FUNCT__ "PCFieldSplitSetIS" 1303 /*@ 1304 PCFieldSplitSetIS - Sets the exact elements for field 1305 1306 Logically Collective on PC 1307 1308 Input Parameters: 1309 + pc - the preconditioner context 1310 . splitname - name of this split, if NULL the number of the split is used 1311 - is - the index set that defines the vector elements in this field 1312 1313 1314 Notes: 1315 Use PCFieldSplitSetFields(), for fields defined by strided types. 1316 1317 This function is called once per split (it creates a new split each time). Solve options 1318 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1319 1320 Level: intermediate 1321 1322 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize() 1323 1324 @*/ 1325 PetscErrorCode PCFieldSplitSetIS(PC pc,const char splitname[],IS is) 1326 { 1327 PetscErrorCode ierr; 1328 1329 PetscFunctionBegin; 1330 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1331 if (splitname) PetscValidCharPointer(splitname,2); 1332 PetscValidHeaderSpecific(is,IS_CLASSID,3); 1333 ierr = PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));CHKERRQ(ierr); 1334 PetscFunctionReturn(0); 1335 } 1336 1337 #undef __FUNCT__ 1338 #define __FUNCT__ "PCFieldSplitGetIS" 1339 /*@ 1340 PCFieldSplitGetIS - Retrieves the elements for a field as an IS 1341 1342 Logically Collective on PC 1343 1344 Input Parameters: 1345 + pc - the preconditioner context 1346 - splitname - name of this split 1347 1348 Output Parameter: 1349 - is - the index set that defines the vector elements in this field, or NULL if the field is not found 1350 1351 Level: intermediate 1352 1353 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetIS() 1354 1355 @*/ 1356 PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is) 1357 { 1358 PetscErrorCode ierr; 1359 1360 PetscFunctionBegin; 1361 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1362 PetscValidCharPointer(splitname,2); 1363 PetscValidPointer(is,3); 1364 { 1365 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 1366 PC_FieldSplitLink ilink = jac->head; 1367 PetscBool found; 1368 1369 *is = NULL; 1370 while (ilink) { 1371 ierr = PetscStrcmp(ilink->splitname, splitname, &found);CHKERRQ(ierr); 1372 if (found) { 1373 *is = ilink->is; 1374 break; 1375 } 1376 ilink = ilink->next; 1377 } 1378 } 1379 PetscFunctionReturn(0); 1380 } 1381 1382 #undef __FUNCT__ 1383 #define __FUNCT__ "PCFieldSplitSetBlockSize" 1384 /*@ 1385 PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the 1386 fieldsplit preconditioner. If not set the matrix block size is used. 1387 1388 Logically Collective on PC 1389 1390 Input Parameters: 1391 + pc - the preconditioner context 1392 - bs - the block size 1393 1394 Level: intermediate 1395 1396 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields() 1397 1398 @*/ 1399 PetscErrorCode PCFieldSplitSetBlockSize(PC pc,PetscInt bs) 1400 { 1401 PetscErrorCode ierr; 1402 1403 PetscFunctionBegin; 1404 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1405 PetscValidLogicalCollectiveInt(pc,bs,2); 1406 ierr = PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));CHKERRQ(ierr); 1407 PetscFunctionReturn(0); 1408 } 1409 1410 #undef __FUNCT__ 1411 #define __FUNCT__ "PCFieldSplitGetSubKSP" 1412 /*@C 1413 PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits 1414 1415 Collective on KSP 1416 1417 Input Parameter: 1418 . pc - the preconditioner context 1419 1420 Output Parameters: 1421 + n - the number of splits 1422 - pc - the array of KSP contexts 1423 1424 Note: 1425 After PCFieldSplitGetSubKSP() the array of KSPs IS to be freed by the user 1426 (not the KSP just the array that contains them). 1427 1428 You must call KSPSetUp() before calling PCFieldSplitGetSubKSP(). 1429 1430 Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs. 1431 You can call PCFieldSplitGetSubKSP(pc,n,NULL_OBJECT,ierr) to determine how large the 1432 KSP array must be. 1433 1434 1435 Level: advanced 1436 1437 .seealso: PCFIELDSPLIT 1438 @*/ 1439 PetscErrorCode PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[]) 1440 { 1441 PetscErrorCode ierr; 1442 1443 PetscFunctionBegin; 1444 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1445 if (n) PetscValidIntPointer(n,2); 1446 ierr = PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));CHKERRQ(ierr); 1447 PetscFunctionReturn(0); 1448 } 1449 1450 #undef __FUNCT__ 1451 #define __FUNCT__ "PCFieldSplitSchurPrecondition" 1452 /*@ 1453 PCFieldSplitSchurPrecondition - Indicates if the Schur complement is preconditioned by a preconditioner constructed by the 1454 A11 matrix. Otherwise no preconditioner is used. 1455 1456 Collective on PC 1457 1458 Input Parameters: 1459 + pc - the preconditioner context 1460 . 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 1461 - userpre - matrix to use for preconditioning, or NULL 1462 1463 Options Database: 1464 . -pc_fieldsplit_schur_precondition <self,selfp,user,a11> default is a11 1465 1466 Notes: 1467 $ If ptype is 1468 $ user then the preconditioner for the Schur complement is generated by the provided matrix (pre argument 1469 $ to this function). 1470 $ a11 then the preconditioner for the Schur complement is generated by the block diagonal part of the original 1471 $ matrix associated with the Schur complement (i.e. A11) 1472 $ self the preconditioner for the Schur complement is generated from the Schur complement matrix itself: 1473 $ The only preconditioner that currently works directly with the Schur complement matrix object is the PCLSC 1474 $ preconditioner 1475 $ selfp then the preconditioning matrix is an explicitly-assembled approximation Sp = A11 - A10 inv(diag(A00)) A01 1476 $ This is only a good preconditioner when diag(A00) is a good preconditioner for A00. 1477 1478 When solving a saddle point problem, where the A11 block is identically zero, using a11 as the ptype only makes sense 1479 with the additional option -fieldsplit_1_pc_type none. Usually for saddle point problems one would use a ptype of self and 1480 -fieldsplit_1_pc_type lsc which uses the least squares commutator to compute a preconditioner for the Schur complement. 1481 1482 Developer Notes: This is a terrible name, which gives no good indication of what the function does. 1483 Among other things, it should have 'Set' in the name, since it sets the type of matrix to use. 1484 1485 Level: intermediate 1486 1487 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCLSC 1488 1489 @*/ 1490 PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1491 { 1492 PetscErrorCode ierr; 1493 1494 PetscFunctionBegin; 1495 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1496 ierr = PetscTryMethod(pc,"PCFieldSplitSchurPrecondition_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));CHKERRQ(ierr); 1497 PetscFunctionReturn(0); 1498 } 1499 1500 #undef __FUNCT__ 1501 #define __FUNCT__ "PCFieldSplitSchurSetSelfpLumping" 1502 /*@ 1503 PCFieldSplitSchurSetSelfpLumping - Indicates whether the Schur complement preconditioner matrix 'selfp' is to be assembled 1504 using a row-lumped A00 matrix. 1505 1506 Not collective 1507 1508 Input Parameters: 1509 + pc - the preconditioner context 1510 - lump - whether to lump A00 when forming Sp = A11 - A10 inv(diag(A00)) A01 1511 1512 Options Database: 1513 + -pc_fieldsplit_schur_precondition_selfp_lump 1514 . Only used with -pc_fieldsplit_schur_precondition selfp or programmatic equivalent PCFieldSplitSchurSetSelfLumping() 1515 - Default is PETSC_FALSE. 1516 1517 1518 Level: advanced 1519 1520 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSchurPrecondition(), PCFieldSplitSchurGetSelfpLumping() 1521 1522 @*/ 1523 PetscErrorCode PCFieldSplitSchurSetSelfpLumping(PC pc,PetscBool lump) 1524 { 1525 PetscErrorCode ierr; 1526 const char* t; 1527 PetscBool isfs; 1528 PC_FieldSplit *jac; 1529 1530 PetscFunctionBegin; 1531 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1532 ierr = PetscObjectGetType((PetscObject)pc,&t);CHKERRQ(ierr); 1533 ierr = PetscStrcmp(t,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1534 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t); 1535 jac = (PC_FieldSplit*)pc->data; 1536 jac->schurp_lump = lump; 1537 PetscFunctionReturn(0); 1538 } 1539 1540 #undef __FUNCT__ 1541 #define __FUNCT__ "PCFieldSplitSchurGetSelfpLumping" 1542 /*@ 1543 PCFieldSplitSchurGetSelfpLumping - Indicates whether the Schur complement preconditioner matrix 'selfp' is to be assembled 1544 using a row-lumped A00 matrix. 1545 1546 Not collective 1547 1548 Input Parameter: 1549 . pc - the preconditioner context 1550 1551 Output Parameter: 1552 . lump - whether to lump A00 when forming Sp = A11 - A10 inv(diag(A00)) A01 1553 1554 Options Database: 1555 + Only used with -pc_fieldsplit_schur_precondition selfp or programmatic equivalent PCFieldSplitSchurSetSelfLumping() 1556 - Default is PETSC_FALSE. 1557 1558 1559 Level: advanced 1560 1561 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSchurPrecondition(), PCFieldSplitSchurSetSelfpLumping() 1562 1563 @*/ 1564 PetscErrorCode PCFieldSplitSchurGetSelfpLumping(PC pc,PetscBool *lump) 1565 { 1566 PetscErrorCode ierr; 1567 const char* t; 1568 PetscBool isfs; 1569 PC_FieldSplit *jac; 1570 1571 PetscFunctionBegin; 1572 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1573 ierr = PetscObjectGetType((PetscObject)pc,&t);CHKERRQ(ierr); 1574 ierr = PetscStrcmp(t,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1575 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t); 1576 jac = (PC_FieldSplit*)pc->data; 1577 if (lump) *lump = jac->schurp_lump; 1578 PetscFunctionReturn(0); 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> - 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 -A10 ksp(A00) ) ( 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_. The action of 1909 ksp(S) is computed using the KSP solver with prefix -fieldsplit_splitname_ (where splitname was given 1910 in providing the SECOND split or 1 if not given). For PCFieldSplitGetKSP() when field number is 0, 1911 it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default 1912 A11 is used to construct a preconditioner for S, use PCFieldSplitSchurPrecondition() to turn on or off this 1913 option. You can use the preconditioner PCLSC to precondition the Schur complement with -fieldsplit_1_pc_type lsc. 1914 When option -fieldsplit_schur_precondition selfp is given, an approximation to S is assembled -- 1915 Sp = A11 - A10 inv(diag(A00)) A01, which has type AIJ and can be used with a variety of preconditioners 1916 (e.g., -fieldsplit_1_pc_type asm). 1917 The factorization type is set using -pc_fieldsplit_schur_fact_type <diag, lower, upper, full>. The full is shown above, 1918 diag gives 1919 $ ( inv(A00) 0 ) 1920 $ ( 0 -ksp(S) ) 1921 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 1922 $ ( A00 0 ) 1923 $ ( A10 S ) 1924 where the inverses of A00 and S are applied using KSPs. The upper factorization is the inverse of 1925 $ ( A00 A01 ) 1926 $ ( 0 S ) 1927 where again the inverses of A00 and S are applied using KSPs. 1928 1929 If only one set of indices (one IS) is provided with PCFieldSplitSetIS() then the complement of that IS 1930 is used automatically for a second block. 1931 1932 The fieldsplit preconditioner cannot currently be used with the BAIJ or SBAIJ data formats if the blocksize is larger than 1. 1933 Generally it should be used with the AIJ format. 1934 1935 The forms of these preconditioners are closely related if not identical to forms derived as "Distributive Iterations", see, 1936 for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling. Note that one can also use PCFIELDSPLIT 1937 inside a smoother resulting in "Distributive Smoothers". 1938 1939 Concepts: physics based preconditioners, block preconditioners 1940 1941 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners, PCLSC, 1942 PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSchurPrecondition() 1943 M*/ 1944 1945 #undef __FUNCT__ 1946 #define __FUNCT__ "PCCreate_FieldSplit" 1947 PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc) 1948 { 1949 PetscErrorCode ierr; 1950 PC_FieldSplit *jac; 1951 1952 PetscFunctionBegin; 1953 ierr = PetscNewLog(pc,&jac);CHKERRQ(ierr); 1954 1955 jac->bs = -1; 1956 jac->nsplits = 0; 1957 jac->type = PC_COMPOSITE_MULTIPLICATIVE; 1958 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */ 1959 jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL; 1960 jac->dm_splits = PETSC_TRUE; 1961 1962 pc->data = (void*)jac; 1963 1964 pc->ops->apply = PCApply_FieldSplit; 1965 pc->ops->applytranspose = PCApplyTranspose_FieldSplit; 1966 pc->ops->setup = PCSetUp_FieldSplit; 1967 pc->ops->reset = PCReset_FieldSplit; 1968 pc->ops->destroy = PCDestroy_FieldSplit; 1969 pc->ops->setfromoptions = PCSetFromOptions_FieldSplit; 1970 pc->ops->view = PCView_FieldSplit; 1971 pc->ops->applyrichardson = 0; 1972 1973 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 1974 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",PCFieldSplitSetFields_FieldSplit);CHKERRQ(ierr); 1975 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",PCFieldSplitSetIS_FieldSplit);CHKERRQ(ierr); 1976 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",PCFieldSplitSetType_FieldSplit);CHKERRQ(ierr); 1977 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",PCFieldSplitSetBlockSize_FieldSplit);CHKERRQ(ierr); 1978 PetscFunctionReturn(0); 1979 } 1980 1981 1982 1983