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