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