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