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