1 2 #include <petsc-private/pcimpl.h> /*I "petscpc.h" I*/ 3 #include <petscdmcomposite.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 ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_default",&fieldsplit_default,NULL);CHKERRQ(ierr); 393 if (stokes) { 394 IS zerodiags,rest; 395 PetscInt nmin,nmax; 396 397 ierr = MatGetOwnershipRange(pc->mat,&nmin,&nmax);CHKERRQ(ierr); 398 ierr = MatFindZeroDiagonals(pc->mat,&zerodiags);CHKERRQ(ierr); 399 ierr = ISComplement(zerodiags,nmin,nmax,&rest);CHKERRQ(ierr); 400 if (jac->reset) { 401 jac->head->is = rest; 402 jac->head->next->is = zerodiags; 403 } else { 404 ierr = PCFieldSplitSetIS(pc,"0",rest);CHKERRQ(ierr); 405 ierr = PCFieldSplitSetIS(pc,"1",zerodiags);CHKERRQ(ierr); 406 } 407 ierr = ISDestroy(&zerodiags);CHKERRQ(ierr); 408 ierr = ISDestroy(&rest);CHKERRQ(ierr); 409 } else { 410 if (jac->reset) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cases not yet handled when PCReset() was used"); 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 if (jac->bs > 0) { 838 ierr = VecSetBlockSize(x,jac->bs);CHKERRQ(ierr); 839 ierr = VecSetBlockSize(y,jac->bs);CHKERRQ(ierr); 840 } 841 CHKMEMQ; 842 if (jac->type == PC_COMPOSITE_ADDITIVE) { 843 if (jac->defaultsplit) { 844 ierr = VecGetBlockSize(x,&bs);CHKERRQ(ierr); 845 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); 846 ierr = VecGetBlockSize(y,&bs);CHKERRQ(ierr); 847 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); 848 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 849 while (ilink) { 850 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 851 ilink = ilink->next; 852 } 853 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 854 } else { 855 ierr = VecSet(y,0.0);CHKERRQ(ierr); 856 while (ilink) { 857 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 858 ilink = ilink->next; 859 } 860 } 861 } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 862 if (!jac->w1) { 863 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 864 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 865 } 866 ierr = VecSet(y,0.0);CHKERRQ(ierr); 867 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 868 cnt = 1; 869 while (ilink->next) { 870 ilink = ilink->next; 871 /* compute the residual only over the part of the vector needed */ 872 ierr = MatMult(jac->Afield[cnt++],y,ilink->x);CHKERRQ(ierr); 873 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 874 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 875 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 876 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 877 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 878 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 879 } 880 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 881 cnt -= 2; 882 while (ilink->previous) { 883 ilink = ilink->previous; 884 /* compute the residual only over the part of the vector needed */ 885 ierr = MatMult(jac->Afield[cnt--],y,ilink->x);CHKERRQ(ierr); 886 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 887 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 888 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 889 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 890 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 891 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 892 } 893 } 894 } else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Unsupported or unknown composition",(int) jac->type); 895 CHKMEMQ; 896 PetscFunctionReturn(0); 897 } 898 899 #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \ 900 (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 901 VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 902 KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) || \ 903 VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \ 904 VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE)) 905 906 #undef __FUNCT__ 907 #define __FUNCT__ "PCApplyTranspose_FieldSplit" 908 static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y) 909 { 910 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 911 PetscErrorCode ierr; 912 PC_FieldSplitLink ilink = jac->head; 913 PetscInt bs; 914 915 PetscFunctionBegin; 916 CHKMEMQ; 917 if (jac->type == PC_COMPOSITE_ADDITIVE) { 918 if (jac->defaultsplit) { 919 ierr = VecGetBlockSize(x,&bs);CHKERRQ(ierr); 920 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); 921 ierr = VecGetBlockSize(y,&bs);CHKERRQ(ierr); 922 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); 923 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 924 while (ilink) { 925 ierr = KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 926 ilink = ilink->next; 927 } 928 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 929 } else { 930 ierr = VecSet(y,0.0);CHKERRQ(ierr); 931 while (ilink) { 932 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 933 ilink = ilink->next; 934 } 935 } 936 } else { 937 if (!jac->w1) { 938 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 939 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 940 } 941 ierr = VecSet(y,0.0);CHKERRQ(ierr); 942 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 943 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 944 while (ilink->next) { 945 ilink = ilink->next; 946 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 947 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 948 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 949 } 950 while (ilink->previous) { 951 ilink = ilink->previous; 952 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 953 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 954 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 955 } 956 } else { 957 while (ilink->next) { /* get to last entry in linked list */ 958 ilink = ilink->next; 959 } 960 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 961 while (ilink->previous) { 962 ilink = ilink->previous; 963 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 964 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 965 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 966 } 967 } 968 } 969 CHKMEMQ; 970 PetscFunctionReturn(0); 971 } 972 973 #undef __FUNCT__ 974 #define __FUNCT__ "PCReset_FieldSplit" 975 static PetscErrorCode PCReset_FieldSplit(PC pc) 976 { 977 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 978 PetscErrorCode ierr; 979 PC_FieldSplitLink ilink = jac->head,next; 980 981 PetscFunctionBegin; 982 while (ilink) { 983 ierr = KSPReset(ilink->ksp);CHKERRQ(ierr); 984 ierr = VecDestroy(&ilink->x);CHKERRQ(ierr); 985 ierr = VecDestroy(&ilink->y);CHKERRQ(ierr); 986 ierr = VecDestroy(&ilink->z);CHKERRQ(ierr); 987 ierr = VecScatterDestroy(&ilink->sctx);CHKERRQ(ierr); 988 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 989 ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr); 990 next = ilink->next; 991 ilink = next; 992 } 993 ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr); 994 if (jac->mat && jac->mat != jac->pmat) { 995 ierr = MatDestroyMatrices(jac->nsplits,&jac->mat);CHKERRQ(ierr); 996 } else if (jac->mat) { 997 jac->mat = NULL; 998 } 999 if (jac->pmat) {ierr = MatDestroyMatrices(jac->nsplits,&jac->pmat);CHKERRQ(ierr);} 1000 if (jac->Afield) {ierr = MatDestroyMatrices(jac->nsplits,&jac->Afield);CHKERRQ(ierr);} 1001 ierr = VecDestroy(&jac->w1);CHKERRQ(ierr); 1002 ierr = VecDestroy(&jac->w2);CHKERRQ(ierr); 1003 ierr = MatDestroy(&jac->schur);CHKERRQ(ierr); 1004 ierr = MatDestroy(&jac->schur_user);CHKERRQ(ierr); 1005 ierr = KSPDestroy(&jac->kspschur);CHKERRQ(ierr); 1006 ierr = KSPDestroy(&jac->kspupper);CHKERRQ(ierr); 1007 ierr = MatDestroy(&jac->B);CHKERRQ(ierr); 1008 ierr = MatDestroy(&jac->C);CHKERRQ(ierr); 1009 jac->reset = PETSC_TRUE; 1010 PetscFunctionReturn(0); 1011 } 1012 1013 #undef __FUNCT__ 1014 #define __FUNCT__ "PCDestroy_FieldSplit" 1015 static PetscErrorCode PCDestroy_FieldSplit(PC pc) 1016 { 1017 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1018 PetscErrorCode ierr; 1019 PC_FieldSplitLink ilink = jac->head,next; 1020 1021 PetscFunctionBegin; 1022 ierr = PCReset_FieldSplit(pc);CHKERRQ(ierr); 1023 while (ilink) { 1024 ierr = KSPDestroy(&ilink->ksp);CHKERRQ(ierr); 1025 next = ilink->next; 1026 ierr = PetscFree(ilink->splitname);CHKERRQ(ierr); 1027 ierr = PetscFree(ilink->fields);CHKERRQ(ierr); 1028 ierr = PetscFree(ilink->fields_col);CHKERRQ(ierr); 1029 ierr = PetscFree(ilink);CHKERRQ(ierr); 1030 ilink = next; 1031 } 1032 ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr); 1033 ierr = PetscFree(pc->data);CHKERRQ(ierr); 1034 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","",NULL);CHKERRQ(ierr); 1035 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetFields_C","",NULL);CHKERRQ(ierr); 1036 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetIS_C","",NULL);CHKERRQ(ierr); 1037 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetType_C","",NULL);CHKERRQ(ierr); 1038 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetBlockSize_C","",NULL);CHKERRQ(ierr); 1039 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSchurPrecondition_C","",NULL);CHKERRQ(ierr); 1040 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetSchurFactType_C","",NULL);CHKERRQ(ierr); 1041 PetscFunctionReturn(0); 1042 } 1043 1044 #undef __FUNCT__ 1045 #define __FUNCT__ "PCSetFromOptions_FieldSplit" 1046 static PetscErrorCode PCSetFromOptions_FieldSplit(PC pc) 1047 { 1048 PetscErrorCode ierr; 1049 PetscInt bs; 1050 PetscBool flg,stokes = PETSC_FALSE; 1051 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1052 PCCompositeType ctype; 1053 1054 PetscFunctionBegin; 1055 ierr = PetscOptionsHead("FieldSplit options");CHKERRQ(ierr); 1056 ierr = PetscOptionsBool("-pc_fieldsplit_real_diagonal","Use diagonal blocks of the operator","PCFieldSplitSetRealDiagonal",jac->realdiagonal,&jac->realdiagonal,NULL);CHKERRQ(ierr); 1057 ierr = PetscOptionsBool("-pc_fieldsplit_dm_splits","Whether to use DMCreateFieldDecomposition() for splits","PCFieldSplitSetDMSplits",jac->dm_splits,&jac->dm_splits,NULL);CHKERRQ(ierr); 1058 ierr = PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);CHKERRQ(ierr); 1059 if (flg) { 1060 ierr = PCFieldSplitSetBlockSize(pc,bs);CHKERRQ(ierr); 1061 } 1062 1063 ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,NULL);CHKERRQ(ierr); 1064 if (stokes) { 1065 ierr = PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);CHKERRQ(ierr); 1066 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_SELF; 1067 } 1068 1069 ierr = PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg);CHKERRQ(ierr); 1070 if (flg) { 1071 ierr = PCFieldSplitSetType(pc,ctype);CHKERRQ(ierr); 1072 } 1073 /* Only setup fields once */ 1074 if ((jac->bs > 0) && (jac->nsplits == 0)) { 1075 /* only allow user to set fields from command line if bs is already known. 1076 otherwise user can set them in PCFieldSplitSetDefaults() */ 1077 ierr = PCFieldSplitSetRuntimeSplits_Private(pc);CHKERRQ(ierr); 1078 if (jac->splitdefined) {ierr = PetscInfo(pc,"Splits defined using the options database\n");CHKERRQ(ierr);} 1079 } 1080 if (jac->type == PC_COMPOSITE_SCHUR) { 1081 ierr = PetscOptionsGetEnum(((PetscObject)pc)->prefix,"-pc_fieldsplit_schur_factorization_type",PCFieldSplitSchurFactTypes,(PetscEnum*)&jac->schurfactorization,&flg);CHKERRQ(ierr); 1082 if (flg) {ierr = PetscInfo(pc,"Deprecated use of -pc_fieldsplit_schur_factorization_type\n");CHKERRQ(ierr);} 1083 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); 1084 ierr = PetscOptionsEnum("-pc_fieldsplit_schur_precondition","How to build preconditioner for Schur complement","PCFieldSplitSchurPrecondition",PCFieldSplitSchurPreTypes,(PetscEnum)jac->schurpre,(PetscEnum*)&jac->schurpre,NULL);CHKERRQ(ierr); 1085 } 1086 ierr = PetscOptionsTail();CHKERRQ(ierr); 1087 PetscFunctionReturn(0); 1088 } 1089 1090 /*------------------------------------------------------------------------------------*/ 1091 1092 EXTERN_C_BEGIN 1093 #undef __FUNCT__ 1094 #define __FUNCT__ "PCFieldSplitSetFields_FieldSplit" 1095 PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 1096 { 1097 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1098 PetscErrorCode ierr; 1099 PC_FieldSplitLink ilink,next = jac->head; 1100 char prefix[128]; 1101 PetscInt i; 1102 1103 PetscFunctionBegin; 1104 if (jac->splitdefined) { 1105 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 1106 PetscFunctionReturn(0); 1107 } 1108 for (i=0; i<n; i++) { 1109 if (fields[i] >= jac->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs); 1110 if (fields[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]); 1111 } 1112 ierr = PetscNew(struct _PC_FieldSplitLink,&ilink);CHKERRQ(ierr); 1113 if (splitname) { 1114 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 1115 } else { 1116 ierr = PetscMalloc(3*sizeof(char),&ilink->splitname);CHKERRQ(ierr); 1117 ierr = PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);CHKERRQ(ierr); 1118 } 1119 ierr = PetscMalloc(n*sizeof(PetscInt),&ilink->fields);CHKERRQ(ierr); 1120 ierr = PetscMemcpy(ilink->fields,fields,n*sizeof(PetscInt));CHKERRQ(ierr); 1121 ierr = PetscMalloc(n*sizeof(PetscInt),&ilink->fields_col);CHKERRQ(ierr); 1122 ierr = PetscMemcpy(ilink->fields_col,fields_col,n*sizeof(PetscInt));CHKERRQ(ierr); 1123 1124 ilink->nfields = n; 1125 ilink->next = NULL; 1126 ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);CHKERRQ(ierr); 1127 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 1128 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 1129 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 1130 1131 ierr = PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);CHKERRQ(ierr); 1132 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 1133 1134 if (!next) { 1135 jac->head = ilink; 1136 ilink->previous = NULL; 1137 } else { 1138 while (next->next) { 1139 next = next->next; 1140 } 1141 next->next = ilink; 1142 ilink->previous = next; 1143 } 1144 jac->nsplits++; 1145 PetscFunctionReturn(0); 1146 } 1147 EXTERN_C_END 1148 1149 EXTERN_C_BEGIN 1150 #undef __FUNCT__ 1151 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit_Schur" 1152 PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp) 1153 { 1154 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1155 PetscErrorCode ierr; 1156 1157 PetscFunctionBegin; 1158 ierr = PetscMalloc(jac->nsplits*sizeof(KSP),subksp);CHKERRQ(ierr); 1159 ierr = MatSchurComplementGetKSP(jac->schur,*subksp);CHKERRQ(ierr); 1160 1161 (*subksp)[1] = jac->kspschur; 1162 if (n) *n = jac->nsplits; 1163 PetscFunctionReturn(0); 1164 } 1165 EXTERN_C_END 1166 1167 EXTERN_C_BEGIN 1168 #undef __FUNCT__ 1169 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit" 1170 PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp) 1171 { 1172 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1173 PetscErrorCode ierr; 1174 PetscInt cnt = 0; 1175 PC_FieldSplitLink ilink = jac->head; 1176 1177 PetscFunctionBegin; 1178 ierr = PetscMalloc(jac->nsplits*sizeof(KSP),subksp);CHKERRQ(ierr); 1179 while (ilink) { 1180 (*subksp)[cnt++] = ilink->ksp; 1181 ilink = ilink->next; 1182 } 1183 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); 1184 if (n) *n = jac->nsplits; 1185 PetscFunctionReturn(0); 1186 } 1187 EXTERN_C_END 1188 1189 EXTERN_C_BEGIN 1190 #undef __FUNCT__ 1191 #define __FUNCT__ "PCFieldSplitSetIS_FieldSplit" 1192 PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is) 1193 { 1194 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1195 PetscErrorCode ierr; 1196 PC_FieldSplitLink ilink, next = jac->head; 1197 char prefix[128]; 1198 1199 PetscFunctionBegin; 1200 if (jac->splitdefined) { 1201 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 1202 PetscFunctionReturn(0); 1203 } 1204 ierr = PetscNew(struct _PC_FieldSplitLink,&ilink);CHKERRQ(ierr); 1205 if (splitname) { 1206 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 1207 } else { 1208 ierr = PetscMalloc(8*sizeof(char),&ilink->splitname);CHKERRQ(ierr); 1209 ierr = PetscSNPrintf(ilink->splitname,7,"%D",jac->nsplits);CHKERRQ(ierr); 1210 } 1211 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 1212 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 1213 ilink->is = is; 1214 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 1215 ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr); 1216 ilink->is_col = is; 1217 ilink->next = NULL; 1218 ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);CHKERRQ(ierr); 1219 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 1220 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 1221 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 1222 1223 ierr = PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);CHKERRQ(ierr); 1224 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 1225 1226 if (!next) { 1227 jac->head = ilink; 1228 ilink->previous = NULL; 1229 } else { 1230 while (next->next) { 1231 next = next->next; 1232 } 1233 next->next = ilink; 1234 ilink->previous = next; 1235 } 1236 jac->nsplits++; 1237 PetscFunctionReturn(0); 1238 } 1239 EXTERN_C_END 1240 1241 #undef __FUNCT__ 1242 #define __FUNCT__ "PCFieldSplitSetFields" 1243 /*@ 1244 PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner 1245 1246 Logically Collective on PC 1247 1248 Input Parameters: 1249 + pc - the preconditioner context 1250 . splitname - name of this split, if NULL the number of the split is used 1251 . n - the number of fields in this split 1252 - fields - the fields in this split 1253 1254 Level: intermediate 1255 1256 Notes: Use PCFieldSplitSetIS() to set a completely general set of indices as a field. 1257 1258 The PCFieldSplitSetFields() is for defining fields as strided blocks. For example, if the block 1259 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 1260 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x.... 1261 where the numbered entries indicate what is in the field. 1262 1263 This function is called once per split (it creates a new split each time). Solve options 1264 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1265 1266 Developer Note: This routine does not actually create the IS representing the split, that is delayed 1267 until PCSetUp_FieldSplit(), because information about the vector/matrix layouts may not be 1268 available when this routine is called. 1269 1270 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize(), PCFieldSplitSetIS() 1271 1272 @*/ 1273 PetscErrorCode PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 1274 { 1275 PetscErrorCode ierr; 1276 1277 PetscFunctionBegin; 1278 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1279 PetscValidCharPointer(splitname,2); 1280 if (n < 1) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Provided number of fields %D in split \"%s\" not positive",n,splitname); 1281 PetscValidIntPointer(fields,3); 1282 ierr = PetscTryMethod(pc,"PCFieldSplitSetFields_C",(PC,const char[],PetscInt,const PetscInt*,const PetscInt*),(pc,splitname,n,fields,fields_col));CHKERRQ(ierr); 1283 PetscFunctionReturn(0); 1284 } 1285 1286 #undef __FUNCT__ 1287 #define __FUNCT__ "PCFieldSplitSetIS" 1288 /*@ 1289 PCFieldSplitSetIS - Sets the exact elements for field 1290 1291 Logically Collective on PC 1292 1293 Input Parameters: 1294 + pc - the preconditioner context 1295 . splitname - name of this split, if NULL the number of the split is used 1296 - is - the index set that defines the vector elements in this field 1297 1298 1299 Notes: 1300 Use PCFieldSplitSetFields(), for fields defined by strided types. 1301 1302 This function is called once per split (it creates a new split each time). Solve options 1303 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1304 1305 Level: intermediate 1306 1307 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize() 1308 1309 @*/ 1310 PetscErrorCode PCFieldSplitSetIS(PC pc,const char splitname[],IS is) 1311 { 1312 PetscErrorCode ierr; 1313 1314 PetscFunctionBegin; 1315 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1316 if (splitname) PetscValidCharPointer(splitname,2); 1317 PetscValidHeaderSpecific(is,IS_CLASSID,3); 1318 ierr = PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));CHKERRQ(ierr); 1319 PetscFunctionReturn(0); 1320 } 1321 1322 #undef __FUNCT__ 1323 #define __FUNCT__ "PCFieldSplitGetIS" 1324 /*@ 1325 PCFieldSplitGetIS - Retrieves the elements for a field as an IS 1326 1327 Logically Collective on PC 1328 1329 Input Parameters: 1330 + pc - the preconditioner context 1331 - splitname - name of this split 1332 1333 Output Parameter: 1334 - is - the index set that defines the vector elements in this field, or NULL if the field is not found 1335 1336 Level: intermediate 1337 1338 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetIS() 1339 1340 @*/ 1341 PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is) 1342 { 1343 PetscErrorCode ierr; 1344 1345 PetscFunctionBegin; 1346 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1347 PetscValidCharPointer(splitname,2); 1348 PetscValidPointer(is,3); 1349 { 1350 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 1351 PC_FieldSplitLink ilink = jac->head; 1352 PetscBool found; 1353 1354 *is = NULL; 1355 while (ilink) { 1356 ierr = PetscStrcmp(ilink->splitname, splitname, &found);CHKERRQ(ierr); 1357 if (found) { 1358 *is = ilink->is; 1359 break; 1360 } 1361 ilink = ilink->next; 1362 } 1363 } 1364 PetscFunctionReturn(0); 1365 } 1366 1367 #undef __FUNCT__ 1368 #define __FUNCT__ "PCFieldSplitSetBlockSize" 1369 /*@ 1370 PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the 1371 fieldsplit preconditioner. If not set the matrix block size is used. 1372 1373 Logically Collective on PC 1374 1375 Input Parameters: 1376 + pc - the preconditioner context 1377 - bs - the block size 1378 1379 Level: intermediate 1380 1381 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields() 1382 1383 @*/ 1384 PetscErrorCode PCFieldSplitSetBlockSize(PC pc,PetscInt bs) 1385 { 1386 PetscErrorCode ierr; 1387 1388 PetscFunctionBegin; 1389 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1390 PetscValidLogicalCollectiveInt(pc,bs,2); 1391 ierr = PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));CHKERRQ(ierr); 1392 PetscFunctionReturn(0); 1393 } 1394 1395 #undef __FUNCT__ 1396 #define __FUNCT__ "PCFieldSplitGetSubKSP" 1397 /*@C 1398 PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits 1399 1400 Collective on KSP 1401 1402 Input Parameter: 1403 . pc - the preconditioner context 1404 1405 Output Parameters: 1406 + n - the number of splits 1407 - pc - the array of KSP contexts 1408 1409 Note: 1410 After PCFieldSplitGetSubKSP() the array of KSPs IS to be freed by the user 1411 (not the KSP just the array that contains them). 1412 1413 You must call KSPSetUp() before calling PCFieldSplitGetSubKSP(). 1414 1415 Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs. 1416 You can call PCFieldSplitGetSubKSP(pc,n,NULL_OBJECT,ierr) to determine how large the 1417 KSP array must be. 1418 1419 1420 Level: advanced 1421 1422 .seealso: PCFIELDSPLIT 1423 @*/ 1424 PetscErrorCode PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[]) 1425 { 1426 PetscErrorCode ierr; 1427 1428 PetscFunctionBegin; 1429 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1430 if (n) PetscValidIntPointer(n,2); 1431 ierr = PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));CHKERRQ(ierr); 1432 PetscFunctionReturn(0); 1433 } 1434 1435 #undef __FUNCT__ 1436 #define __FUNCT__ "PCFieldSplitSchurPrecondition" 1437 /*@ 1438 PCFieldSplitSchurPrecondition - Indicates if the Schur complement is preconditioned by a preconditioner constructed by the 1439 A11 matrix. Otherwise no preconditioner is used. 1440 1441 Collective on PC 1442 1443 Input Parameters: 1444 + pc - the preconditioner context 1445 . ptype - which matrix to use for preconditioning the Schur complement, PC_FIELDSPLIT_SCHUR_PRE_A11 (diag) is default 1446 - userpre - matrix to use for preconditioning, or NULL 1447 1448 Options Database: 1449 . -pc_fieldsplit_schur_precondition <self,user,a11> default is a11 1450 1451 Notes: 1452 $ If ptype is 1453 $ user then the preconditioner for the Schur complement is generated by the provided matrix (pre argument 1454 $ to this function). 1455 $ a11 then the preconditioner for the Schur complement is generated by the block diagonal part of the original 1456 $ matrix associated with the Schur complement (i.e. A11) 1457 $ self the preconditioner for the Schur complement is generated from the Schur complement matrix itself: 1458 $ The only preconditioner that currently works directly with the Schur complement matrix object is the PCLSC 1459 $ preconditioner 1460 1461 When solving a saddle point problem, where the A11 block is identically zero, using a11 as the ptype only makes sense 1462 with the additional option -fieldsplit_1_pc_type none. Usually for saddle point problems one would use a ptype of self and 1463 -fieldsplit_1_pc_type lsc which uses the least squares commutator compute a preconditioner for the Schur complement. 1464 1465 Developer Notes: This is a terrible name, gives no good indication of what the function does and should also have Set in 1466 the name since it sets a proceedure to use. 1467 1468 Level: intermediate 1469 1470 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCLSC 1471 1472 @*/ 1473 PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1474 { 1475 PetscErrorCode ierr; 1476 1477 PetscFunctionBegin; 1478 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1479 ierr = PetscTryMethod(pc,"PCFieldSplitSchurPrecondition_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));CHKERRQ(ierr); 1480 PetscFunctionReturn(0); 1481 } 1482 1483 EXTERN_C_BEGIN 1484 #undef __FUNCT__ 1485 #define __FUNCT__ "PCFieldSplitSchurPrecondition_FieldSplit" 1486 PetscErrorCode PCFieldSplitSchurPrecondition_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1487 { 1488 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1489 PetscErrorCode ierr; 1490 1491 PetscFunctionBegin; 1492 jac->schurpre = ptype; 1493 if (pre) { 1494 ierr = MatDestroy(&jac->schur_user);CHKERRQ(ierr); 1495 jac->schur_user = pre; 1496 ierr = PetscObjectReference((PetscObject)jac->schur_user);CHKERRQ(ierr); 1497 } 1498 PetscFunctionReturn(0); 1499 } 1500 EXTERN_C_END 1501 1502 #undef __FUNCT__ 1503 #define __FUNCT__ "PCFieldSplitSetSchurFactType" 1504 /*@ 1505 PCFieldSplitSetSchurFactType - sets which blocks of the approximate block factorization to retain 1506 1507 Collective on PC 1508 1509 Input Parameters: 1510 + pc - the preconditioner context 1511 - ftype - which blocks of factorization to retain, PC_FIELDSPLIT_SCHUR_FACT_FULL is default 1512 1513 Options Database: 1514 . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> default is full 1515 1516 1517 Level: intermediate 1518 1519 Notes: 1520 The FULL factorization is 1521 1522 $ (A B) = (1 0) (A 0) (1 Ainv*B) 1523 $ (C D) (C*Ainv 1) (0 S) (0 1 ) 1524 1525 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, 1526 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). 1527 1528 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 1529 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 1530 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, 1531 the preconditioned operator has three distinct nonzero eigenvalues and minimal polynomial of degree at most 4, so KSPGMRES converges in at most 4 iterations. 1532 1533 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 1534 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). 1535 1536 References: 1537 Murphy, Golub, and Wathen, A note on preconditioning indefinite linear systems, SIAM J. Sci. Comput., 21 (2000) pp. 1969-1972. 1538 1539 Ipsen, A note on preconditioning nonsymmetric matrices, SIAM J. Sci. Comput., 23 (2001), pp. 1050-1051. 1540 1541 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType 1542 @*/ 1543 PetscErrorCode PCFieldSplitSetSchurFactType(PC pc,PCFieldSplitSchurFactType ftype) 1544 { 1545 PetscErrorCode ierr; 1546 1547 PetscFunctionBegin; 1548 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1549 ierr = PetscTryMethod(pc,"PCFieldSplitSetSchurFactType_C",(PC,PCFieldSplitSchurFactType),(pc,ftype));CHKERRQ(ierr); 1550 PetscFunctionReturn(0); 1551 } 1552 1553 #undef __FUNCT__ 1554 #define __FUNCT__ "PCFieldSplitSetSchurFactType_FieldSplit" 1555 PETSC_EXTERN_C PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc,PCFieldSplitSchurFactType ftype) 1556 { 1557 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1558 1559 PetscFunctionBegin; 1560 jac->schurfactorization = ftype; 1561 PetscFunctionReturn(0); 1562 } 1563 1564 #undef __FUNCT__ 1565 #define __FUNCT__ "PCFieldSplitGetSchurBlocks" 1566 /*@C 1567 PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement 1568 1569 Collective on KSP 1570 1571 Input Parameter: 1572 . pc - the preconditioner context 1573 1574 Output Parameters: 1575 + A00 - the (0,0) block 1576 . A01 - the (0,1) block 1577 . A10 - the (1,0) block 1578 - A11 - the (1,1) block 1579 1580 Level: advanced 1581 1582 .seealso: PCFIELDSPLIT 1583 @*/ 1584 PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11) 1585 { 1586 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 1587 1588 PetscFunctionBegin; 1589 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1590 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach."); 1591 if (A00) *A00 = jac->pmat[0]; 1592 if (A01) *A01 = jac->B; 1593 if (A10) *A10 = jac->C; 1594 if (A11) *A11 = jac->pmat[1]; 1595 PetscFunctionReturn(0); 1596 } 1597 1598 EXTERN_C_BEGIN 1599 #undef __FUNCT__ 1600 #define __FUNCT__ "PCFieldSplitSetType_FieldSplit" 1601 PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type) 1602 { 1603 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1604 PetscErrorCode ierr; 1605 1606 PetscFunctionBegin; 1607 jac->type = type; 1608 if (type == PC_COMPOSITE_SCHUR) { 1609 pc->ops->apply = PCApply_FieldSplit_Schur; 1610 pc->ops->view = PCView_FieldSplit_Schur; 1611 1612 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit_Schur",PCFieldSplitGetSubKSP_FieldSplit_Schur);CHKERRQ(ierr); 1613 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSchurPrecondition_C","PCFieldSplitSchurPrecondition_FieldSplit",PCFieldSplitSchurPrecondition_FieldSplit);CHKERRQ(ierr); 1614 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetSchurFactType_C","PCFieldSplitSetSchurFactType_FieldSplit",PCFieldSplitSetSchurFactType_FieldSplit);CHKERRQ(ierr); 1615 1616 } else { 1617 pc->ops->apply = PCApply_FieldSplit; 1618 pc->ops->view = PCView_FieldSplit; 1619 1620 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit",PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 1621 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSchurPrecondition_C","",0);CHKERRQ(ierr); 1622 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetSchurFactType_C","",0);CHKERRQ(ierr); 1623 } 1624 PetscFunctionReturn(0); 1625 } 1626 EXTERN_C_END 1627 1628 EXTERN_C_BEGIN 1629 #undef __FUNCT__ 1630 #define __FUNCT__ "PCFieldSplitSetBlockSize_FieldSplit" 1631 PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs) 1632 { 1633 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1634 1635 PetscFunctionBegin; 1636 if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs); 1637 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); 1638 jac->bs = bs; 1639 PetscFunctionReturn(0); 1640 } 1641 EXTERN_C_END 1642 1643 #undef __FUNCT__ 1644 #define __FUNCT__ "PCFieldSplitSetType" 1645 /*@ 1646 PCFieldSplitSetType - Sets the type of fieldsplit preconditioner. 1647 1648 Collective on PC 1649 1650 Input Parameter: 1651 . pc - the preconditioner context 1652 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 1653 1654 Options Database Key: 1655 . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type 1656 1657 Level: Intermediate 1658 1659 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 1660 1661 .seealso: PCCompositeSetType() 1662 1663 @*/ 1664 PetscErrorCode PCFieldSplitSetType(PC pc,PCCompositeType type) 1665 { 1666 PetscErrorCode ierr; 1667 1668 PetscFunctionBegin; 1669 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1670 ierr = PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));CHKERRQ(ierr); 1671 PetscFunctionReturn(0); 1672 } 1673 1674 #undef __FUNCT__ 1675 #define __FUNCT__ "PCFieldSplitGetType" 1676 /*@ 1677 PCFieldSplitGetType - Gets the type of fieldsplit preconditioner. 1678 1679 Not collective 1680 1681 Input Parameter: 1682 . pc - the preconditioner context 1683 1684 Output Parameter: 1685 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 1686 1687 Level: Intermediate 1688 1689 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 1690 .seealso: PCCompositeSetType() 1691 @*/ 1692 PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type) 1693 { 1694 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 1695 1696 PetscFunctionBegin; 1697 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1698 PetscValidIntPointer(type,2); 1699 *type = jac->type; 1700 PetscFunctionReturn(0); 1701 } 1702 1703 #undef __FUNCT__ 1704 #define __FUNCT__ "PCFieldSplitSetDMSplits" 1705 /*@ 1706 PCFieldSplitSetDMSplits - Flags whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible. 1707 1708 Logically Collective 1709 1710 Input Parameters: 1711 + pc - the preconditioner context 1712 - flg - boolean indicating whether to use field splits defined by the DM 1713 1714 Options Database Key: 1715 . -pc_fieldsplit_dm_splits 1716 1717 Level: Intermediate 1718 1719 .keywords: PC, DM, composite preconditioner, additive, multiplicative 1720 1721 .seealso: PCFieldSplitGetDMSplits() 1722 1723 @*/ 1724 PetscErrorCode PCFieldSplitSetDMSplits(PC pc,PetscBool flg) 1725 { 1726 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1727 PetscBool isfs; 1728 PetscErrorCode ierr; 1729 1730 PetscFunctionBegin; 1731 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1732 PetscValidLogicalCollectiveBool(pc,flg,2); 1733 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1734 if (isfs) { 1735 jac->dm_splits = flg; 1736 } 1737 PetscFunctionReturn(0); 1738 } 1739 1740 1741 #undef __FUNCT__ 1742 #define __FUNCT__ "PCFieldSplitGetDMSplits" 1743 /*@ 1744 PCFieldSplitGetDMSplits - Returns flag indicating whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible. 1745 1746 Logically Collective 1747 1748 Input Parameter: 1749 . pc - the preconditioner context 1750 1751 Output Parameter: 1752 . flg - boolean indicating whether to use field splits defined by the DM 1753 1754 Level: Intermediate 1755 1756 .keywords: PC, DM, composite preconditioner, additive, multiplicative 1757 1758 .seealso: PCFieldSplitSetDMSplits() 1759 1760 @*/ 1761 PetscErrorCode PCFieldSplitGetDMSplits(PC pc,PetscBool* flg) 1762 { 1763 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1764 PetscBool isfs; 1765 PetscErrorCode ierr; 1766 1767 PetscFunctionBegin; 1768 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1769 PetscValidPointer(flg,2); 1770 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1771 if (isfs) { 1772 if(flg) *flg = jac->dm_splits; 1773 } 1774 PetscFunctionReturn(0); 1775 } 1776 1777 /* -------------------------------------------------------------------------------------*/ 1778 /*MC 1779 PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual 1780 fields or groups of fields. See the users manual section "Solving Block Matrices" for more details. 1781 1782 To set options on the solvers for each block append -fieldsplit_ to all the PC 1783 options database keys. For example, -fieldsplit_pc_type ilu -fieldsplit_pc_factor_levels 1 1784 1785 To set the options on the solvers separate for each block call PCFieldSplitGetSubKSP() 1786 and set the options directly on the resulting KSP object 1787 1788 Level: intermediate 1789 1790 Options Database Keys: 1791 + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split 1792 . -pc_fieldsplit_default - automatically add any fields to additional splits that have not 1793 been supplied explicitly by -pc_fieldsplit_%d_fields 1794 . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields) 1795 . -pc_fieldsplit_type <additive,multiplicative,symmetric_multiplicative,schur> - type of relaxation or factorization splitting 1796 . -pc_fieldsplit_schur_precondition <self,user,a11> - default is a11 1797 . -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero or negative diagonal and uses Schur complement with no preconditioner as the solver 1798 1799 - Options prefix for inner solvers when using Schur complement preconditioner are -fieldsplit_0_ and -fieldsplit_1_ 1800 for all other solvers they are -fieldsplit_%d_ for the dth field, use -fieldsplit_ for all fields 1801 1802 Notes: 1803 Use PCFieldSplitSetFields() to set fields defined by "strided" entries and PCFieldSplitSetIS() 1804 to define a field by an arbitrary collection of entries. 1805 1806 If no fields are set the default is used. The fields are defined by entries strided by bs, 1807 beginning at 0 then 1, etc to bs-1. The block size can be set with PCFieldSplitSetBlockSize(), 1808 if this is not called the block size defaults to the blocksize of the second matrix passed 1809 to KSPSetOperators()/PCSetOperators(). 1810 1811 $ For the Schur complement preconditioner if J = ( A00 A01 ) 1812 $ ( A10 A11 ) 1813 $ the preconditioner using full factorization is 1814 $ ( I -A10 ksp(A00) ) ( inv(A00) 0 ) ( I 0 ) 1815 $ ( 0 I ) ( 0 ksp(S) ) ( -A10 ksp(A00) I ) 1816 where the action of inv(A00) is applied using the KSP solver with prefix -fieldsplit_0_. The action of 1817 ksp(S) is computed using the KSP solver with prefix -fieldsplit_splitname_ (where splitname was given 1818 in providing the SECOND split or 1 if not give). For PCFieldSplitGetKSP() when field number is 0, 1819 it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default 1820 A11 is used to construct a preconditioner for S, use PCFieldSplitSchurPrecondition() to turn on or off this 1821 option. You can use the preconditioner PCLSC to precondition the Schur complement with -fieldsplit_1_pc_type lsc. The 1822 factorization type is set using -pc_fieldsplit_schur_fact_type <diag, lower, upper, full>. The full is shown above, 1823 diag gives 1824 $ ( inv(A00) 0 ) 1825 $ ( 0 -ksp(S) ) 1826 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 1827 $ ( A00 0 ) 1828 $ ( A10 S ) 1829 where the inverses of A00 and S are applied using KSPs. The upper factorization is the inverse of 1830 $ ( A00 A01 ) 1831 $ ( 0 S ) 1832 where again the inverses of A00 and S are applied using KSPs. 1833 1834 If only one set of indices (one IS) is provided with PCFieldSplitSetIS() then the complement of that IS 1835 is used automatically for a second block. 1836 1837 The fieldsplit preconditioner cannot currently be used with the BAIJ or SBAIJ data formats if the blocksize is larger than 1. 1838 Generally it should be used with the AIJ format. 1839 1840 The forms of these preconditioners are closely related if not identical to forms derived as "Distributive Iterations", see, 1841 for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling. Note that one can also use PCFIELDSPLIT 1842 inside a smoother resulting in "Distributive Smoothers". 1843 1844 Concepts: physics based preconditioners, block preconditioners 1845 1846 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners, PCLSC, 1847 PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSchurPrecondition() 1848 M*/ 1849 1850 EXTERN_C_BEGIN 1851 #undef __FUNCT__ 1852 #define __FUNCT__ "PCCreate_FieldSplit" 1853 PetscErrorCode PCCreate_FieldSplit(PC pc) 1854 { 1855 PetscErrorCode ierr; 1856 PC_FieldSplit *jac; 1857 1858 PetscFunctionBegin; 1859 ierr = PetscNewLog(pc,PC_FieldSplit,&jac);CHKERRQ(ierr); 1860 1861 jac->bs = -1; 1862 jac->nsplits = 0; 1863 jac->type = PC_COMPOSITE_MULTIPLICATIVE; 1864 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */ 1865 jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL; 1866 jac->dm_splits = PETSC_FALSE; 1867 1868 pc->data = (void*)jac; 1869 1870 pc->ops->apply = PCApply_FieldSplit; 1871 pc->ops->applytranspose = PCApplyTranspose_FieldSplit; 1872 pc->ops->setup = PCSetUp_FieldSplit; 1873 pc->ops->reset = PCReset_FieldSplit; 1874 pc->ops->destroy = PCDestroy_FieldSplit; 1875 pc->ops->setfromoptions = PCSetFromOptions_FieldSplit; 1876 pc->ops->view = PCView_FieldSplit; 1877 pc->ops->applyrichardson = 0; 1878 1879 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit", 1880 PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 1881 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetFields_C","PCFieldSplitSetFields_FieldSplit", 1882 PCFieldSplitSetFields_FieldSplit);CHKERRQ(ierr); 1883 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetIS_C","PCFieldSplitSetIS_FieldSplit", 1884 PCFieldSplitSetIS_FieldSplit);CHKERRQ(ierr); 1885 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetType_C","PCFieldSplitSetType_FieldSplit", 1886 PCFieldSplitSetType_FieldSplit);CHKERRQ(ierr); 1887 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetBlockSize_C","PCFieldSplitSetBlockSize_FieldSplit", 1888 PCFieldSplitSetBlockSize_FieldSplit);CHKERRQ(ierr); 1889 PetscFunctionReturn(0); 1890 } 1891 EXTERN_C_END 1892 1893 1894 1895