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