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