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