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