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