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