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