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