1 #define PETSCKSP_DLL 2 3 #include "private/pcimpl.h" /*I "petscpc.h" I*/ 4 5 const char *const PCFieldSplitSchurPreTypes[] = {"SELF","DIAG","USER","PCFieldSplitSchurPreType","PC_FIELDSPLIT_SCHUR_PRE_",0}; 6 const char *const PCFieldSplitSchurFactorizationTypes[] = {"DIAG","LOWER","UPPER","FULL","PCFieldSplitSchurFactorizationType","PC_FIELDSPLIT_SCHUR_FACTORIZATION_",0}; 7 8 typedef enum { 9 PC_FIELDSPLIT_SCHUR_FACTORIZATION_DIAG, 10 PC_FIELDSPLIT_SCHUR_FACTORIZATION_LOWER, 11 PC_FIELDSPLIT_SCHUR_FACTORIZATION_UPPER, 12 PC_FIELDSPLIT_SCHUR_FACTORIZATION_FULL 13 } PCFieldSplitSchurFactorizationType; 14 15 typedef struct _PC_FieldSplitLink *PC_FieldSplitLink; 16 struct _PC_FieldSplitLink { 17 KSP ksp; 18 Vec x,y; 19 char *splitname; 20 PetscInt nfields; 21 PetscInt *fields; 22 VecScatter sctx; 23 IS is; 24 PC_FieldSplitLink next,previous; 25 }; 26 27 typedef struct { 28 PCCompositeType type; 29 PetscBool defaultsplit; /* Flag for a system with a set of 'k' scalar fields with the same layout (and bs = k) */ 30 PetscBool splitdefined; /* Flag is set after the splits have been defined, to prevent more splits from being added */ 31 PetscBool realdiagonal; /* Flag to use the diagonal blocks of mat preconditioned by pmat, instead of just pmat */ 32 PetscInt bs; /* Block size for IS and Mat structures */ 33 PetscInt nsplits; /* Number of field divisions defined */ 34 Vec *x,*y,w1,w2; 35 Mat *mat; /* The diagonal block for each split */ 36 Mat *pmat; /* The preconditioning diagonal block for each split */ 37 Mat *Afield; /* The rows of the matrix associated with each split */ 38 PetscBool issetup; 39 /* Only used when Schur complement preconditioning is used */ 40 Mat B; /* The (0,1) block */ 41 Mat C; /* The (1,0) block */ 42 Mat schur; /* The Schur complement S = D - C A^{-1} B */ 43 Mat schur_user; /* User-provided preconditioning matrix for the Schur complement */ 44 PCFieldSplitSchurPreType schurpre; /* Determines which preconditioning matrix is used for the Schur complement */ 45 PCFieldSplitSchurFactorizationType schurfactorization; 46 KSP kspschur; /* The solver for S */ 47 PC_FieldSplitLink head; 48 } PC_FieldSplit; 49 50 /* 51 Notes: there is no particular reason that pmat, x, and y are stored as arrays in PC_FieldSplit instead of 52 inside PC_FieldSplitLink, just historical. If you want to be able to add new fields after already using the 53 PC you could change this. 54 */ 55 56 /* This helper is so that setting a user-provided preconditioning matrix is orthogonal to choosing to use it. This way the 57 * application-provided FormJacobian can provide this matrix without interfering with the user's (command-line) choices. */ 58 static Mat FieldSplitSchurPre(PC_FieldSplit *jac) 59 { 60 switch (jac->schurpre) { 61 case PC_FIELDSPLIT_SCHUR_PRE_SELF: return jac->schur; 62 case PC_FIELDSPLIT_SCHUR_PRE_DIAG: return jac->pmat[1]; 63 case PC_FIELDSPLIT_SCHUR_PRE_USER: /* Use a user-provided matrix if it is given, otherwise diagonal block */ 64 default: 65 return jac->schur_user ? jac->schur_user : jac->pmat[1]; 66 } 67 } 68 69 70 #undef __FUNCT__ 71 #define __FUNCT__ "PCView_FieldSplit" 72 static PetscErrorCode PCView_FieldSplit(PC pc,PetscViewer viewer) 73 { 74 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 75 PetscErrorCode ierr; 76 PetscBool iascii; 77 PetscInt i,j; 78 PC_FieldSplitLink ilink = jac->head; 79 80 PetscFunctionBegin; 81 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 82 if (iascii) { 83 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D, blocksize = %D\n",PCCompositeTypes[jac->type],jac->nsplits,jac->bs);CHKERRQ(ierr); 84 ierr = PetscViewerASCIIPrintf(viewer," Solver info for each split is in the following KSP objects:\n");CHKERRQ(ierr); 85 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 86 for (i=0; i<jac->nsplits; i++) { 87 if (ilink->fields) { 88 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);CHKERRQ(ierr); 89 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 90 for (j=0; j<ilink->nfields; j++) { 91 if (j > 0) { 92 ierr = PetscViewerASCIIPrintf(viewer,",");CHKERRQ(ierr); 93 } 94 ierr = PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);CHKERRQ(ierr); 95 } 96 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 97 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 98 } else { 99 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i);CHKERRQ(ierr); 100 } 101 ierr = KSPView(ilink->ksp,viewer);CHKERRQ(ierr); 102 ilink = ilink->next; 103 } 104 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 105 } else { 106 SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for PCFieldSplit",((PetscObject)viewer)->type_name); 107 } 108 PetscFunctionReturn(0); 109 } 110 111 #undef __FUNCT__ 112 #define __FUNCT__ "PCView_FieldSplit_Schur" 113 static PetscErrorCode PCView_FieldSplit_Schur(PC pc,PetscViewer viewer) 114 { 115 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 116 PetscErrorCode ierr; 117 PetscBool iascii; 118 PetscInt i,j; 119 PC_FieldSplitLink ilink = jac->head; 120 KSP ksp; 121 122 PetscFunctionBegin; 123 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 124 if (iascii) { 125 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit with Schur preconditioner, blocksize = %D, factorization %s\n",jac->bs,PCFieldSplitSchurFactorizationTypes[jac->schurfactorization]);CHKERRQ(ierr); 126 ierr = PetscViewerASCIIPrintf(viewer," Split info:\n");CHKERRQ(ierr); 127 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 128 for (i=0; i<jac->nsplits; i++) { 129 if (ilink->fields) { 130 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);CHKERRQ(ierr); 131 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 132 for (j=0; j<ilink->nfields; j++) { 133 if (j > 0) { 134 ierr = PetscViewerASCIIPrintf(viewer,",");CHKERRQ(ierr); 135 } 136 ierr = PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);CHKERRQ(ierr); 137 } 138 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 139 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 140 } else { 141 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i);CHKERRQ(ierr); 142 } 143 ilink = ilink->next; 144 } 145 ierr = PetscViewerASCIIPrintf(viewer,"KSP solver for A block \n");CHKERRQ(ierr); 146 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 147 if (jac->schur) { 148 ierr = MatSchurComplementGetKSP(jac->schur,&ksp);CHKERRQ(ierr); 149 ierr = KSPView(ksp,viewer);CHKERRQ(ierr); 150 } else { 151 ierr = PetscViewerASCIIPrintf(viewer," not yet available\n");CHKERRQ(ierr); 152 } 153 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 154 ierr = PetscViewerASCIIPrintf(viewer,"KSP solver for S = D - C inv(A) B \n");CHKERRQ(ierr); 155 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 156 if (jac->kspschur) { 157 ierr = KSPView(jac->kspschur,viewer);CHKERRQ(ierr); 158 } else { 159 ierr = PetscViewerASCIIPrintf(viewer," not yet available\n");CHKERRQ(ierr); 160 } 161 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 162 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 163 } else { 164 SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for PCFieldSplit",((PetscObject)viewer)->type_name); 165 } 166 PetscFunctionReturn(0); 167 } 168 169 #undef __FUNCT__ 170 #define __FUNCT__ "PCFieldSplitSetRuntimeSplits_Private" 171 /* Precondition: jac->bs is set to a meaningful value */ 172 static PetscErrorCode PCFieldSplitSetRuntimeSplits_Private(PC pc) 173 { 174 PetscErrorCode ierr; 175 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 176 PetscInt i,nfields,*ifields; 177 PetscBool flg; 178 char optionname[128],splitname[8]; 179 180 PetscFunctionBegin; 181 ierr = PetscMalloc(jac->bs*sizeof(PetscInt),&ifields);CHKERRQ(ierr); 182 for (i=0,flg=PETSC_TRUE; ; i++) { 183 ierr = PetscSNPrintf(splitname,sizeof splitname,"%D",i);CHKERRQ(ierr); 184 ierr = PetscSNPrintf(optionname,sizeof optionname,"-pc_fieldsplit_%D_fields",i);CHKERRQ(ierr); 185 nfields = jac->bs; 186 ierr = PetscOptionsGetIntArray(((PetscObject)pc)->prefix,optionname,ifields,&nfields,&flg);CHKERRQ(ierr); 187 if (!flg) break; 188 if (!nfields) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot list zero fields"); 189 ierr = PCFieldSplitSetFields(pc,splitname,nfields,ifields);CHKERRQ(ierr); 190 } 191 if (i > 0) { 192 /* Makes command-line setting of splits take precedence over setting them in code. 193 Otherwise subsequent calls to PCFieldSplitSetIS() or PCFieldSplitSetFields() would 194 create new splits, which would probably not be what the user wanted. */ 195 jac->splitdefined = PETSC_TRUE; 196 } 197 ierr = PetscFree(ifields);CHKERRQ(ierr); 198 PetscFunctionReturn(0); 199 } 200 201 #undef __FUNCT__ 202 #define __FUNCT__ "PCFieldSplitSetDefaults" 203 static PetscErrorCode PCFieldSplitSetDefaults(PC pc) 204 { 205 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 206 PetscErrorCode ierr; 207 PC_FieldSplitLink ilink = jac->head; 208 PetscBool flg = PETSC_FALSE; 209 PetscInt i; 210 211 PetscFunctionBegin; 212 if (!ilink) { 213 if (pc->dm) { 214 PetscBool dmcomposite; 215 ierr = PetscTypeCompare((PetscObject)pc->dm,DMCOMPOSITE,&dmcomposite);CHKERRQ(ierr); 216 if (dmcomposite) { 217 PetscInt nDM; 218 IS *fields; 219 ierr = PetscInfo(pc,"Setting up physics based fieldsplit preconditioner using the embedded DM\n");CHKERRQ(ierr); 220 ierr = DMCompositeGetNumberDM(pc->dm,&nDM);CHKERRQ(ierr); 221 ierr = DMCompositeGetGlobalISs(pc->dm,&fields);CHKERRQ(ierr); 222 for (i=0; i<nDM; i++) { 223 char splitname[8]; 224 ierr = PetscSNPrintf(splitname,sizeof splitname,"%D",i);CHKERRQ(ierr); 225 ierr = PCFieldSplitSetIS(pc,splitname,fields[i]);CHKERRQ(ierr); 226 ierr = ISDestroy(fields[i]);CHKERRQ(ierr); 227 } 228 ierr = PetscFree(fields);CHKERRQ(ierr); 229 } 230 } else { 231 if (jac->bs <= 0) { 232 if (pc->pmat) { 233 ierr = MatGetBlockSize(pc->pmat,&jac->bs);CHKERRQ(ierr); 234 } else { 235 jac->bs = 1; 236 } 237 } 238 239 ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_default",&flg,PETSC_NULL);CHKERRQ(ierr); 240 if (!flg) { 241 /* Allow user to set fields from command line, if bs was known at the time of PCSetFromOptions_FieldSplit() 242 then it is set there. This is not ideal because we should only have options set in XXSetFromOptions(). */ 243 ierr = PCFieldSplitSetRuntimeSplits_Private(pc);CHKERRQ(ierr); 244 if (jac->splitdefined) {ierr = PetscInfo(pc,"Splits defined using the options database\n");CHKERRQ(ierr);} 245 } 246 if (flg || !jac->splitdefined) { 247 ierr = PetscInfo(pc,"Using default splitting of fields\n");CHKERRQ(ierr); 248 for (i=0; i<jac->bs; i++) { 249 char splitname[8]; 250 ierr = PetscSNPrintf(splitname,sizeof splitname,"%D",i);CHKERRQ(ierr); 251 ierr = PCFieldSplitSetFields(pc,splitname,1,&i);CHKERRQ(ierr); 252 } 253 jac->defaultsplit = PETSC_TRUE; 254 } 255 } 256 } else if (jac->nsplits == 1) { 257 if (ilink->is) { 258 IS is2; 259 PetscInt nmin,nmax; 260 261 ierr = MatGetOwnershipRange(pc->mat,&nmin,&nmax);CHKERRQ(ierr); 262 ierr = ISComplement(ilink->is,nmin,nmax,&is2);CHKERRQ(ierr); 263 ierr = PCFieldSplitSetIS(pc,"1",is2);CHKERRQ(ierr); 264 ierr = ISDestroy(is2);CHKERRQ(ierr); 265 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Must provide at least two sets of fields to PCFieldSplit()"); 266 } 267 if (jac->nsplits < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unhandled case, must have at least two fields"); 268 PetscFunctionReturn(0); 269 } 270 271 272 #undef __FUNCT__ 273 #define __FUNCT__ "PCSetUp_FieldSplit" 274 static PetscErrorCode PCSetUp_FieldSplit(PC pc) 275 { 276 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 277 PetscErrorCode ierr; 278 PC_FieldSplitLink ilink; 279 PetscInt i,nsplit,ccsize; 280 MatStructure flag = pc->flag; 281 PetscBool sorted; 282 283 PetscFunctionBegin; 284 ierr = PCFieldSplitSetDefaults(pc);CHKERRQ(ierr); 285 nsplit = jac->nsplits; 286 ilink = jac->head; 287 288 /* get the matrices for each split */ 289 if (!jac->issetup) { 290 PetscInt rstart,rend,nslots,bs; 291 292 jac->issetup = PETSC_TRUE; 293 294 /* This is done here instead of in PCFieldSplitSetFields() because may not have matrix at that point */ 295 bs = jac->bs; 296 ierr = MatGetOwnershipRange(pc->pmat,&rstart,&rend);CHKERRQ(ierr); 297 ierr = MatGetLocalSize(pc->pmat,PETSC_NULL,&ccsize);CHKERRQ(ierr); 298 nslots = (rend - rstart)/bs; 299 for (i=0; i<nsplit; i++) { 300 if (jac->defaultsplit) { 301 ierr = ISCreateStride(((PetscObject)pc)->comm,nslots,rstart+i,nsplit,&ilink->is);CHKERRQ(ierr); 302 } else if (!ilink->is) { 303 if (ilink->nfields > 1) { 304 PetscInt *ii,j,k,nfields = ilink->nfields,*fields = ilink->fields; 305 ierr = PetscMalloc(ilink->nfields*nslots*sizeof(PetscInt),&ii);CHKERRQ(ierr); 306 for (j=0; j<nslots; j++) { 307 for (k=0; k<nfields; k++) { 308 ii[nfields*j + k] = rstart + bs*j + fields[k]; 309 } 310 } 311 ierr = ISCreateGeneral(((PetscObject)pc)->comm,nslots*nfields,ii,PETSC_OWN_POINTER,&ilink->is);CHKERRQ(ierr); 312 } else { 313 ierr = ISCreateStride(((PetscObject)pc)->comm,nslots,rstart+ilink->fields[0],bs,&ilink->is);CHKERRQ(ierr); 314 } 315 } 316 ierr = ISSorted(ilink->is,&sorted);CHKERRQ(ierr); 317 if (!sorted) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Fields must be sorted when creating split"); 318 ilink = ilink->next; 319 } 320 } 321 322 ilink = jac->head; 323 if (!jac->pmat) { 324 ierr = PetscMalloc(nsplit*sizeof(Mat),&jac->pmat);CHKERRQ(ierr); 325 for (i=0; i<nsplit; i++) { 326 ierr = MatGetSubMatrix(pc->pmat,ilink->is,ilink->is,MAT_INITIAL_MATRIX,&jac->pmat[i]);CHKERRQ(ierr); 327 ilink = ilink->next; 328 } 329 } else { 330 for (i=0; i<nsplit; i++) { 331 ierr = MatGetSubMatrix(pc->pmat,ilink->is,ilink->is,MAT_REUSE_MATRIX,&jac->pmat[i]);CHKERRQ(ierr); 332 ilink = ilink->next; 333 } 334 } 335 if (jac->realdiagonal) { 336 ilink = jac->head; 337 if (!jac->mat) { 338 ierr = PetscMalloc(nsplit*sizeof(Mat),&jac->mat);CHKERRQ(ierr); 339 for (i=0; i<nsplit; i++) { 340 ierr = MatGetSubMatrix(pc->mat,ilink->is,ilink->is,MAT_INITIAL_MATRIX,&jac->mat[i]);CHKERRQ(ierr); 341 ilink = ilink->next; 342 } 343 } else { 344 for (i=0; i<nsplit; i++) { 345 ierr = MatGetSubMatrix(pc->mat,ilink->is,ilink->is,MAT_REUSE_MATRIX,&jac->mat[i]);CHKERRQ(ierr); 346 ilink = ilink->next; 347 } 348 } 349 } else { 350 jac->mat = jac->pmat; 351 } 352 353 if (jac->type != PC_COMPOSITE_ADDITIVE && jac->type != PC_COMPOSITE_SCHUR) { 354 /* extract the rows of the matrix associated with each field: used for efficient computation of residual inside algorithm */ 355 ilink = jac->head; 356 if (!jac->Afield) { 357 ierr = PetscMalloc(nsplit*sizeof(Mat),&jac->Afield);CHKERRQ(ierr); 358 for (i=0; i<nsplit; i++) { 359 ierr = MatGetSubMatrix(pc->mat,ilink->is,PETSC_NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);CHKERRQ(ierr); 360 ilink = ilink->next; 361 } 362 } else { 363 for (i=0; i<nsplit; i++) { 364 ierr = MatGetSubMatrix(pc->mat,ilink->is,PETSC_NULL,MAT_REUSE_MATRIX,&jac->Afield[i]);CHKERRQ(ierr); 365 ilink = ilink->next; 366 } 367 } 368 } 369 370 if (jac->type == PC_COMPOSITE_SCHUR) { 371 IS ccis; 372 PetscInt rstart,rend; 373 if (nsplit != 2) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_INCOMP,"To use Schur complement preconditioner you must have exactly 2 fields"); 374 375 /* When extracting off-diagonal submatrices, we take complements from this range */ 376 ierr = MatGetOwnershipRangeColumn(pc->mat,&rstart,&rend);CHKERRQ(ierr); 377 378 /* need to handle case when one is resetting up the preconditioner */ 379 if (jac->schur) { 380 ilink = jac->head; 381 ierr = ISComplement(ilink->is,rstart,rend,&ccis);CHKERRQ(ierr); 382 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->B);CHKERRQ(ierr); 383 ierr = ISDestroy(ccis);CHKERRQ(ierr); 384 ilink = ilink->next; 385 ierr = ISComplement(ilink->is,rstart,rend,&ccis);CHKERRQ(ierr); 386 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->C);CHKERRQ(ierr); 387 ierr = ISDestroy(ccis);CHKERRQ(ierr); 388 ierr = MatSchurComplementUpdate(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->pmat[1],pc->flag);CHKERRQ(ierr); 389 ierr = KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac),pc->flag);CHKERRQ(ierr); 390 391 } else { 392 KSP ksp; 393 char schurprefix[256]; 394 395 /* extract the B and C matrices */ 396 ilink = jac->head; 397 ierr = ISComplement(ilink->is,rstart,rend,&ccis);CHKERRQ(ierr); 398 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);CHKERRQ(ierr); 399 ierr = ISDestroy(ccis);CHKERRQ(ierr); 400 ilink = ilink->next; 401 ierr = ISComplement(ilink->is,rstart,rend,&ccis);CHKERRQ(ierr); 402 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);CHKERRQ(ierr); 403 ierr = ISDestroy(ccis);CHKERRQ(ierr); 404 /* Better would be to use 'mat[0]' (diagonal block of the real matrix) preconditioned by pmat[0] */ 405 ierr = MatCreateSchurComplement(jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1],&jac->schur);CHKERRQ(ierr); 406 ierr = MatSchurComplementGetKSP(jac->schur,&ksp);CHKERRQ(ierr); 407 ierr = PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,2);CHKERRQ(ierr); 408 ierr = MatSetFromOptions(jac->schur);CHKERRQ(ierr); 409 410 ierr = KSPCreate(((PetscObject)pc)->comm,&jac->kspschur);CHKERRQ(ierr); 411 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->kspschur);CHKERRQ(ierr); 412 ierr = PetscObjectIncrementTabLevel((PetscObject)jac->kspschur,(PetscObject)pc,1);CHKERRQ(ierr); 413 ierr = KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac),DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 414 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) { 415 PC pc; 416 ierr = KSPGetPC(jac->kspschur,&pc);CHKERRQ(ierr); 417 ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr); 418 /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */ 419 } 420 ierr = PetscSNPrintf(schurprefix,sizeof schurprefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",ilink->splitname);CHKERRQ(ierr); 421 ierr = KSPSetOptionsPrefix(jac->kspschur,schurprefix);CHKERRQ(ierr); 422 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 423 ierr = KSPSetFromOptions(jac->kspschur);CHKERRQ(ierr); 424 425 ierr = PetscMalloc2(2,Vec,&jac->x,2,Vec,&jac->y);CHKERRQ(ierr); 426 ierr = MatGetVecs(jac->pmat[0],&jac->x[0],&jac->y[0]);CHKERRQ(ierr); 427 ierr = MatGetVecs(jac->pmat[1],&jac->x[1],&jac->y[1]);CHKERRQ(ierr); 428 ilink = jac->head; 429 ilink->x = jac->x[0]; ilink->y = jac->y[0]; 430 ilink = ilink->next; 431 ilink->x = jac->x[1]; ilink->y = jac->y[1]; 432 } 433 } else { 434 /* set up the individual PCs */ 435 i = 0; 436 ilink = jac->head; 437 while (ilink) { 438 ierr = KSPSetOperators(ilink->ksp,jac->mat[i],jac->pmat[i],flag);CHKERRQ(ierr); 439 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 440 ierr = KSPSetFromOptions(ilink->ksp);CHKERRQ(ierr); 441 ierr = KSPSetUp(ilink->ksp);CHKERRQ(ierr); 442 i++; 443 ilink = ilink->next; 444 } 445 446 /* create work vectors for each split */ 447 if (!jac->x) { 448 ierr = PetscMalloc2(nsplit,Vec,&jac->x,nsplit,Vec,&jac->y);CHKERRQ(ierr); 449 ilink = jac->head; 450 for (i=0; i<nsplit; i++) { 451 Vec *vl,*vr; 452 453 ierr = KSPGetVecs(ilink->ksp,1,&vr,1,&vl);CHKERRQ(ierr); 454 ilink->x = *vr; 455 ilink->y = *vl; 456 ierr = PetscFree(vr);CHKERRQ(ierr); 457 ierr = PetscFree(vl);CHKERRQ(ierr); 458 jac->x[i] = ilink->x; 459 jac->y[i] = ilink->y; 460 ilink = ilink->next; 461 } 462 } 463 } 464 465 466 if (!jac->head->sctx) { 467 Vec xtmp; 468 469 /* compute scatter contexts needed by multiplicative versions and non-default splits */ 470 471 ilink = jac->head; 472 ierr = MatGetVecs(pc->pmat,&xtmp,PETSC_NULL);CHKERRQ(ierr); 473 for (i=0; i<nsplit; i++) { 474 ierr = VecScatterCreate(xtmp,ilink->is,jac->x[i],PETSC_NULL,&ilink->sctx);CHKERRQ(ierr); 475 ilink = ilink->next; 476 } 477 ierr = VecDestroy(xtmp);CHKERRQ(ierr); 478 } 479 PetscFunctionReturn(0); 480 } 481 482 #define FieldSplitSplitSolveAdd(ilink,xx,yy) \ 483 (VecScatterBegin(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 484 VecScatterEnd(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 485 KSPSolve(ilink->ksp,ilink->x,ilink->y) || \ 486 VecScatterBegin(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE) || \ 487 VecScatterEnd(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE)) 488 489 #undef __FUNCT__ 490 #define __FUNCT__ "PCApply_FieldSplit_Schur" 491 static PetscErrorCode PCApply_FieldSplit_Schur(PC pc,Vec x,Vec y) 492 { 493 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 494 PetscErrorCode ierr; 495 KSP ksp; 496 PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next; 497 498 PetscFunctionBegin; 499 ierr = MatSchurComplementGetKSP(jac->schur,&ksp);CHKERRQ(ierr); 500 501 switch (jac->schurfactorization) { 502 case PC_FIELDSPLIT_SCHUR_FACTORIZATION_DIAG: 503 /* [A 0; 0 -S], positive definite, suitable for MINRES */ 504 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 505 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 506 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 507 ierr = KSPSolve(ksp,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 508 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 509 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 510 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 511 ierr = VecScale(ilinkD->y,-1.);CHKERRQ(ierr); 512 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 513 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 514 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 515 break; 516 case PC_FIELDSPLIT_SCHUR_FACTORIZATION_LOWER: 517 /* [A 0; C S], suitable for left preconditioning */ 518 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 519 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 520 ierr = KSPSolve(ksp,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 521 ierr = MatMult(jac->C,ilinkA->y,ilinkD->x);CHKERRQ(ierr); 522 ierr = VecScale(ilinkD->x,-1.);CHKERRQ(ierr); 523 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 524 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 525 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 526 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 527 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 528 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 529 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 530 break; 531 case PC_FIELDSPLIT_SCHUR_FACTORIZATION_UPPER: 532 /* [A B; 0 S], suitable for right preconditioning */ 533 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 534 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 535 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 536 ierr = MatMult(jac->B,ilinkD->y,ilinkA->x);CHKERRQ(ierr); 537 ierr = VecScale(ilinkA->x,-1.);CHKERRQ(ierr); 538 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 539 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 540 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 541 ierr = KSPSolve(ksp,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 542 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 543 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 544 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 545 break; 546 case PC_FIELDSPLIT_SCHUR_FACTORIZATION_FULL: 547 /* [1 0; CA^{-1} 1] [A 0; 0 S] [1 A^{-1}B; 0 1], an exact solve if applied exactly, needs one extra solve with A */ 548 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 549 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 550 ierr = KSPSolve(ksp,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 551 ierr = MatMult(jac->C,ilinkA->y,ilinkD->x);CHKERRQ(ierr); 552 ierr = VecScale(ilinkD->x,-1.0);CHKERRQ(ierr); 553 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 554 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 555 556 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 557 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 558 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 559 560 ierr = MatMult(jac->B,ilinkD->y,ilinkA->y);CHKERRQ(ierr); 561 ierr = VecAXPY(ilinkA->x,-1.0,ilinkA->y);CHKERRQ(ierr); 562 ierr = KSPSolve(ksp,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 563 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 564 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 565 } 566 PetscFunctionReturn(0); 567 } 568 569 #undef __FUNCT__ 570 #define __FUNCT__ "PCApply_FieldSplit" 571 static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y) 572 { 573 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 574 PetscErrorCode ierr; 575 PC_FieldSplitLink ilink = jac->head; 576 PetscInt cnt; 577 578 PetscFunctionBegin; 579 CHKMEMQ; 580 ierr = VecSetBlockSize(x,jac->bs);CHKERRQ(ierr); 581 ierr = VecSetBlockSize(y,jac->bs);CHKERRQ(ierr); 582 583 if (jac->type == PC_COMPOSITE_ADDITIVE) { 584 if (jac->defaultsplit) { 585 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 586 while (ilink) { 587 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 588 ilink = ilink->next; 589 } 590 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 591 } else { 592 ierr = VecSet(y,0.0);CHKERRQ(ierr); 593 while (ilink) { 594 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 595 ilink = ilink->next; 596 } 597 } 598 } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 599 if (!jac->w1) { 600 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 601 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 602 } 603 ierr = VecSet(y,0.0);CHKERRQ(ierr); 604 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 605 cnt = 1; 606 while (ilink->next) { 607 ilink = ilink->next; 608 /* compute the residual only over the part of the vector needed */ 609 ierr = MatMult(jac->Afield[cnt++],y,ilink->x);CHKERRQ(ierr); 610 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 611 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 612 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 613 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 614 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 615 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 616 } 617 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 618 cnt -= 2; 619 while (ilink->previous) { 620 ilink = ilink->previous; 621 /* compute the residual only over the part of the vector needed */ 622 ierr = MatMult(jac->Afield[cnt--],y,ilink->x);CHKERRQ(ierr); 623 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 624 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 625 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 626 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 627 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 628 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 629 } 630 } 631 } else SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Unsupported or unknown composition",(int) jac->type); 632 CHKMEMQ; 633 PetscFunctionReturn(0); 634 } 635 636 #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \ 637 (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 638 VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 639 KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) || \ 640 VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \ 641 VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE)) 642 643 #undef __FUNCT__ 644 #define __FUNCT__ "PCApplyTranspose_FieldSplit" 645 static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y) 646 { 647 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 648 PetscErrorCode ierr; 649 PC_FieldSplitLink ilink = jac->head; 650 651 PetscFunctionBegin; 652 CHKMEMQ; 653 ierr = VecSetBlockSize(x,jac->bs);CHKERRQ(ierr); 654 ierr = VecSetBlockSize(y,jac->bs);CHKERRQ(ierr); 655 656 if (jac->type == PC_COMPOSITE_ADDITIVE) { 657 if (jac->defaultsplit) { 658 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 659 while (ilink) { 660 ierr = KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 661 ilink = ilink->next; 662 } 663 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 664 } else { 665 ierr = VecSet(y,0.0);CHKERRQ(ierr); 666 while (ilink) { 667 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 668 ilink = ilink->next; 669 } 670 } 671 } else { 672 if (!jac->w1) { 673 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 674 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 675 } 676 ierr = VecSet(y,0.0);CHKERRQ(ierr); 677 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 678 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 679 while (ilink->next) { 680 ilink = ilink->next; 681 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 682 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 683 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 684 } 685 while (ilink->previous) { 686 ilink = ilink->previous; 687 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 688 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 689 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 690 } 691 } else { 692 while (ilink->next) { /* get to last entry in linked list */ 693 ilink = ilink->next; 694 } 695 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 696 while (ilink->previous) { 697 ilink = ilink->previous; 698 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 699 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 700 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 701 } 702 } 703 } 704 CHKMEMQ; 705 PetscFunctionReturn(0); 706 } 707 708 #undef __FUNCT__ 709 #define __FUNCT__ "PCDestroy_FieldSplit" 710 static PetscErrorCode PCDestroy_FieldSplit(PC pc) 711 { 712 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 713 PetscErrorCode ierr; 714 PC_FieldSplitLink ilink = jac->head,next; 715 716 PetscFunctionBegin; 717 while (ilink) { 718 ierr = KSPDestroy(ilink->ksp);CHKERRQ(ierr); 719 if (ilink->x) {ierr = VecDestroy(ilink->x);CHKERRQ(ierr);} 720 if (ilink->y) {ierr = VecDestroy(ilink->y);CHKERRQ(ierr);} 721 if (ilink->sctx) {ierr = VecScatterDestroy(ilink->sctx);CHKERRQ(ierr);} 722 if (ilink->is) {ierr = ISDestroy(ilink->is);CHKERRQ(ierr);} 723 next = ilink->next; 724 ierr = PetscFree(ilink->splitname);CHKERRQ(ierr); 725 ierr = PetscFree(ilink->fields);CHKERRQ(ierr); 726 ierr = PetscFree(ilink);CHKERRQ(ierr); 727 ilink = next; 728 } 729 ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr); 730 if (jac->mat && jac->mat != jac->pmat) {ierr = MatDestroyMatrices(jac->nsplits,&jac->mat);CHKERRQ(ierr);} 731 if (jac->pmat) {ierr = MatDestroyMatrices(jac->nsplits,&jac->pmat);CHKERRQ(ierr);} 732 if (jac->Afield) {ierr = MatDestroyMatrices(jac->nsplits,&jac->Afield);CHKERRQ(ierr);} 733 if (jac->w1) {ierr = VecDestroy(jac->w1);CHKERRQ(ierr);} 734 if (jac->w2) {ierr = VecDestroy(jac->w2);CHKERRQ(ierr);} 735 if (jac->schur) {ierr = MatDestroy(jac->schur);CHKERRQ(ierr);} 736 if (jac->schur_user) {ierr = MatDestroy(jac->schur_user);CHKERRQ(ierr);} 737 if (jac->kspschur) {ierr = KSPDestroy(jac->kspschur);CHKERRQ(ierr);} 738 if (jac->B) {ierr = MatDestroy(jac->B);CHKERRQ(ierr);} 739 if (jac->C) {ierr = MatDestroy(jac->C);CHKERRQ(ierr);} 740 ierr = PetscFree(jac);CHKERRQ(ierr); 741 PetscFunctionReturn(0); 742 } 743 744 #undef __FUNCT__ 745 #define __FUNCT__ "PCSetFromOptions_FieldSplit" 746 static PetscErrorCode PCSetFromOptions_FieldSplit(PC pc) 747 { 748 PetscErrorCode ierr; 749 PetscInt bs; 750 PetscBool flg; 751 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 752 PCCompositeType ctype; 753 754 PetscFunctionBegin; 755 ierr = PetscOptionsHead("FieldSplit options");CHKERRQ(ierr); 756 ierr = PetscOptionsBool("-pc_fieldsplit_real_diagonal","Use diagonal blocks of the operator","PCFieldSplitSetRealDiagonal",jac->realdiagonal,&jac->realdiagonal,PETSC_NULL);CHKERRQ(ierr); 757 ierr = PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);CHKERRQ(ierr); 758 if (flg) { 759 ierr = PCFieldSplitSetBlockSize(pc,bs);CHKERRQ(ierr); 760 } 761 762 ierr = PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg);CHKERRQ(ierr); 763 if (flg) { 764 ierr = PCFieldSplitSetType(pc,ctype);CHKERRQ(ierr); 765 } 766 767 /* Only setup fields once */ 768 if ((jac->bs > 0) && (jac->nsplits == 0)) { 769 /* only allow user to set fields from command line if bs is already known. 770 otherwise user can set them in PCFieldSplitSetDefaults() */ 771 ierr = PCFieldSplitSetRuntimeSplits_Private(pc);CHKERRQ(ierr); 772 if (jac->splitdefined) {ierr = PetscInfo(pc,"Splits defined using the options database\n");CHKERRQ(ierr);} 773 } 774 if (jac->type == PC_COMPOSITE_SCHUR) { 775 ierr = PetscOptionsEnum("-pc_fieldsplit_schur_factorization_type","Factorization to use","None",PCFieldSplitSchurFactorizationTypes,(PetscEnum)jac->schurfactorization,(PetscEnum*)&jac->schurfactorization,PETSC_NULL);CHKERRQ(ierr); 776 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); 777 } 778 ierr = PetscOptionsTail();CHKERRQ(ierr); 779 PetscFunctionReturn(0); 780 } 781 782 /*------------------------------------------------------------------------------------*/ 783 784 EXTERN_C_BEGIN 785 #undef __FUNCT__ 786 #define __FUNCT__ "PCFieldSplitSetFields_FieldSplit" 787 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields) 788 { 789 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 790 PetscErrorCode ierr; 791 PC_FieldSplitLink ilink,next = jac->head; 792 char prefix[128]; 793 PetscInt i; 794 795 PetscFunctionBegin; 796 if (jac->splitdefined) { 797 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 798 PetscFunctionReturn(0); 799 } 800 for (i=0; i<n; i++) { 801 if (fields[i] >= jac->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs); 802 if (fields[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]); 803 } 804 ierr = PetscNew(struct _PC_FieldSplitLink,&ilink);CHKERRQ(ierr); 805 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 806 ierr = PetscMalloc(n*sizeof(PetscInt),&ilink->fields);CHKERRQ(ierr); 807 ierr = PetscMemcpy(ilink->fields,fields,n*sizeof(PetscInt));CHKERRQ(ierr); 808 ilink->nfields = n; 809 ilink->next = PETSC_NULL; 810 ierr = KSPCreate(((PetscObject)pc)->comm,&ilink->ksp);CHKERRQ(ierr); 811 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 812 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 813 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 814 815 ierr = PetscSNPrintf(prefix,sizeof prefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",splitname);CHKERRQ(ierr); 816 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 817 818 if (!next) { 819 jac->head = ilink; 820 ilink->previous = PETSC_NULL; 821 } else { 822 while (next->next) { 823 next = next->next; 824 } 825 next->next = ilink; 826 ilink->previous = next; 827 } 828 jac->nsplits++; 829 PetscFunctionReturn(0); 830 } 831 EXTERN_C_END 832 833 EXTERN_C_BEGIN 834 #undef __FUNCT__ 835 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit_Schur" 836 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp) 837 { 838 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 839 PetscErrorCode ierr; 840 841 PetscFunctionBegin; 842 ierr = PetscMalloc(jac->nsplits*sizeof(KSP),subksp);CHKERRQ(ierr); 843 ierr = MatSchurComplementGetKSP(jac->schur,*subksp);CHKERRQ(ierr); 844 (*subksp)[1] = jac->kspschur; 845 *n = jac->nsplits; 846 PetscFunctionReturn(0); 847 } 848 EXTERN_C_END 849 850 EXTERN_C_BEGIN 851 #undef __FUNCT__ 852 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit" 853 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp) 854 { 855 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 856 PetscErrorCode ierr; 857 PetscInt cnt = 0; 858 PC_FieldSplitLink ilink = jac->head; 859 860 PetscFunctionBegin; 861 ierr = PetscMalloc(jac->nsplits*sizeof(KSP*),subksp);CHKERRQ(ierr); 862 while (ilink) { 863 (*subksp)[cnt++] = ilink->ksp; 864 ilink = ilink->next; 865 } 866 if (cnt != jac->nsplits) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Corrupt PCFIELDSPLIT object: number splits in linked list %D in object %D",cnt,jac->nsplits); 867 *n = jac->nsplits; 868 PetscFunctionReturn(0); 869 } 870 EXTERN_C_END 871 872 EXTERN_C_BEGIN 873 #undef __FUNCT__ 874 #define __FUNCT__ "PCFieldSplitSetIS_FieldSplit" 875 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is) 876 { 877 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 878 PetscErrorCode ierr; 879 PC_FieldSplitLink ilink, next = jac->head; 880 char prefix[128]; 881 882 PetscFunctionBegin; 883 if (jac->splitdefined) { 884 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 885 PetscFunctionReturn(0); 886 } 887 ierr = PetscNew(struct _PC_FieldSplitLink,&ilink);CHKERRQ(ierr); 888 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 889 ilink->is = is; 890 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 891 ilink->next = PETSC_NULL; 892 ierr = KSPCreate(((PetscObject)pc)->comm,&ilink->ksp);CHKERRQ(ierr); 893 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 894 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 895 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 896 897 ierr = PetscSNPrintf(prefix,sizeof prefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",splitname);CHKERRQ(ierr); 898 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 899 900 if (!next) { 901 jac->head = ilink; 902 ilink->previous = PETSC_NULL; 903 } else { 904 while (next->next) { 905 next = next->next; 906 } 907 next->next = ilink; 908 ilink->previous = next; 909 } 910 jac->nsplits++; 911 912 PetscFunctionReturn(0); 913 } 914 EXTERN_C_END 915 916 #undef __FUNCT__ 917 #define __FUNCT__ "PCFieldSplitSetFields" 918 /*@ 919 PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner 920 921 Logically Collective on PC 922 923 Input Parameters: 924 + pc - the preconditioner context 925 . splitname - name of this split 926 . n - the number of fields in this split 927 - fields - the fields in this split 928 929 Level: intermediate 930 931 Notes: Use PCFieldSplitSetIS() to set a completely general set of indices as a field. 932 933 The PCFieldSplitSetFields() is for defining fields as a strided blocks. For example, if the block 934 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 935 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x.... 936 where the numbered entries indicate what is in the field. 937 938 This function is called once per split (it creates a new split each time). Solve options 939 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 940 941 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize(), PCFieldSplitSetIS() 942 943 @*/ 944 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields) 945 { 946 PetscErrorCode ierr; 947 948 PetscFunctionBegin; 949 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 950 PetscValidCharPointer(splitname,2); 951 if (n < 1) SETERRQ2(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Provided number of fields %D in split \"%s\" not positive",n,splitname); 952 PetscValidIntPointer(fields,3); 953 ierr = PetscTryMethod(pc,"PCFieldSplitSetFields_C",(PC,const char[],PetscInt,const PetscInt *),(pc,splitname,n,fields));CHKERRQ(ierr); 954 PetscFunctionReturn(0); 955 } 956 957 #undef __FUNCT__ 958 #define __FUNCT__ "PCFieldSplitSetIS" 959 /*@ 960 PCFieldSplitSetIS - Sets the exact elements for field 961 962 Logically Collective on PC 963 964 Input Parameters: 965 + pc - the preconditioner context 966 . splitname - name of this split 967 - is - the index set that defines the vector elements in this field 968 969 970 Notes: 971 Use PCFieldSplitSetFields(), for fields defined by strided types. 972 973 This function is called once per split (it creates a new split each time). Solve options 974 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 975 976 Level: intermediate 977 978 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize() 979 980 @*/ 981 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetIS(PC pc,const char splitname[],IS is) 982 { 983 PetscErrorCode ierr; 984 985 PetscFunctionBegin; 986 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 987 PetscValidCharPointer(splitname,2); 988 PetscValidHeaderSpecific(is,IS_CLASSID,3); 989 ierr = PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));CHKERRQ(ierr); 990 PetscFunctionReturn(0); 991 } 992 993 #undef __FUNCT__ 994 #define __FUNCT__ "PCFieldSplitSetBlockSize" 995 /*@ 996 PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the 997 fieldsplit preconditioner. If not set the matrix block size is used. 998 999 Logically Collective on PC 1000 1001 Input Parameters: 1002 + pc - the preconditioner context 1003 - bs - the block size 1004 1005 Level: intermediate 1006 1007 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields() 1008 1009 @*/ 1010 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetBlockSize(PC pc,PetscInt bs) 1011 { 1012 PetscErrorCode ierr; 1013 1014 PetscFunctionBegin; 1015 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1016 PetscValidLogicalCollectiveInt(pc,bs,2); 1017 ierr = PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));CHKERRQ(ierr); 1018 PetscFunctionReturn(0); 1019 } 1020 1021 #undef __FUNCT__ 1022 #define __FUNCT__ "PCFieldSplitGetSubKSP" 1023 /*@C 1024 PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits 1025 1026 Collective on KSP 1027 1028 Input Parameter: 1029 . pc - the preconditioner context 1030 1031 Output Parameters: 1032 + n - the number of split 1033 - pc - the array of KSP contexts 1034 1035 Note: 1036 After PCFieldSplitGetSubKSP() the array of KSPs IS to be freed by the user 1037 (not the KSP just the array that contains them). 1038 1039 You must call KSPSetUp() before calling PCFieldSplitGetSubKSP(). 1040 1041 Level: advanced 1042 1043 .seealso: PCFIELDSPLIT 1044 @*/ 1045 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[]) 1046 { 1047 PetscErrorCode ierr; 1048 1049 PetscFunctionBegin; 1050 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1051 PetscValidIntPointer(n,2); 1052 ierr = PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));CHKERRQ(ierr); 1053 PetscFunctionReturn(0); 1054 } 1055 1056 #undef __FUNCT__ 1057 #define __FUNCT__ "PCFieldSplitSchurPrecondition" 1058 /*@ 1059 PCFieldSplitSchurPrecondition - Indicates if the Schur complement is preconditioned by a preconditioner constructed by the 1060 D matrix. Otherwise no preconditioner is used. 1061 1062 Collective on PC 1063 1064 Input Parameters: 1065 + pc - the preconditioner context 1066 . ptype - which matrix to use for preconditioning the Schur complement 1067 - userpre - matrix to use for preconditioning, or PETSC_NULL 1068 1069 Notes: 1070 The default is to use the block on the diagonal of the preconditioning matrix. This is D, in the (1,1) position. 1071 There are currently no preconditioners that work directly with the Schur complement so setting 1072 PC_FIELDSPLIT_SCHUR_PRE_SELF is observationally equivalent to -fieldsplit_1_pc_type none. 1073 1074 Options Database: 1075 . -pc_fieldsplit_schur_precondition <self,user,diag> default is diag 1076 1077 Level: intermediate 1078 1079 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType 1080 1081 @*/ 1082 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1083 { 1084 PetscErrorCode ierr; 1085 1086 PetscFunctionBegin; 1087 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1088 ierr = PetscTryMethod(pc,"PCFieldSplitSchurPrecondition_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));CHKERRQ(ierr); 1089 PetscFunctionReturn(0); 1090 } 1091 1092 EXTERN_C_BEGIN 1093 #undef __FUNCT__ 1094 #define __FUNCT__ "PCFieldSplitSchurPrecondition_FieldSplit" 1095 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSchurPrecondition_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1096 { 1097 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1098 PetscErrorCode ierr; 1099 1100 PetscFunctionBegin; 1101 jac->schurpre = ptype; 1102 if (pre) { 1103 if (jac->schur_user) {ierr = MatDestroy(jac->schur_user);CHKERRQ(ierr);} 1104 jac->schur_user = pre; 1105 ierr = PetscObjectReference((PetscObject)jac->schur_user);CHKERRQ(ierr); 1106 } 1107 PetscFunctionReturn(0); 1108 } 1109 EXTERN_C_END 1110 1111 #undef __FUNCT__ 1112 #define __FUNCT__ "PCFieldSplitGetSchurBlocks" 1113 /*@C 1114 PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement 1115 1116 Collective on KSP 1117 1118 Input Parameter: 1119 . pc - the preconditioner context 1120 1121 Output Parameters: 1122 + A - the (0,0) block 1123 . B - the (0,1) block 1124 . C - the (1,0) block 1125 - D - the (1,1) block 1126 1127 Level: advanced 1128 1129 .seealso: PCFIELDSPLIT 1130 @*/ 1131 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitGetSchurBlocks(PC pc,Mat *A,Mat *B,Mat *C, Mat *D) 1132 { 1133 PC_FieldSplit *jac = (PC_FieldSplit *) pc->data; 1134 1135 PetscFunctionBegin; 1136 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1137 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach."); 1138 if (A) *A = jac->pmat[0]; 1139 if (B) *B = jac->B; 1140 if (C) *C = jac->C; 1141 if (D) *D = jac->pmat[1]; 1142 PetscFunctionReturn(0); 1143 } 1144 1145 EXTERN_C_BEGIN 1146 #undef __FUNCT__ 1147 #define __FUNCT__ "PCFieldSplitSetType_FieldSplit" 1148 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type) 1149 { 1150 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1151 PetscErrorCode ierr; 1152 1153 PetscFunctionBegin; 1154 jac->type = type; 1155 if (type == PC_COMPOSITE_SCHUR) { 1156 pc->ops->apply = PCApply_FieldSplit_Schur; 1157 pc->ops->view = PCView_FieldSplit_Schur; 1158 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit_Schur",PCFieldSplitGetSubKSP_FieldSplit_Schur);CHKERRQ(ierr); 1159 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSchurPrecondition_C","PCFieldSplitSchurPrecondition_FieldSplit",PCFieldSplitSchurPrecondition_FieldSplit);CHKERRQ(ierr); 1160 1161 } else { 1162 pc->ops->apply = PCApply_FieldSplit; 1163 pc->ops->view = PCView_FieldSplit; 1164 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit",PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 1165 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSchurPrecondition_C","",0);CHKERRQ(ierr); 1166 } 1167 PetscFunctionReturn(0); 1168 } 1169 EXTERN_C_END 1170 1171 EXTERN_C_BEGIN 1172 #undef __FUNCT__ 1173 #define __FUNCT__ "PCFieldSplitSetBlockSize_FieldSplit" 1174 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs) 1175 { 1176 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1177 1178 PetscFunctionBegin; 1179 if (bs < 1) SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs); 1180 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); 1181 jac->bs = bs; 1182 PetscFunctionReturn(0); 1183 } 1184 EXTERN_C_END 1185 1186 #undef __FUNCT__ 1187 #define __FUNCT__ "PCFieldSplitSetType" 1188 /*@ 1189 PCFieldSplitSetType - Sets the type of fieldsplit preconditioner. 1190 1191 Collective on PC 1192 1193 Input Parameter: 1194 . pc - the preconditioner context 1195 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 1196 1197 Options Database Key: 1198 . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type 1199 1200 Level: Developer 1201 1202 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 1203 1204 .seealso: PCCompositeSetType() 1205 1206 @*/ 1207 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetType(PC pc,PCCompositeType type) 1208 { 1209 PetscErrorCode ierr; 1210 1211 PetscFunctionBegin; 1212 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1213 ierr = PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));CHKERRQ(ierr); 1214 PetscFunctionReturn(0); 1215 } 1216 1217 /* -------------------------------------------------------------------------------------*/ 1218 /*MC 1219 PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual 1220 fields or groups of fields 1221 1222 1223 To set options on the solvers for each block append -fieldsplit_ to all the PC 1224 options database keys. For example, -fieldsplit_pc_type ilu -fieldsplit_pc_factor_levels 1 1225 1226 To set the options on the solvers separate for each block call PCFieldSplitGetSubKSP() 1227 and set the options directly on the resulting KSP object 1228 1229 Level: intermediate 1230 1231 Options Database Keys: 1232 + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split 1233 . -pc_fieldsplit_default - automatically add any fields to additional splits that have not 1234 been supplied explicitly by -pc_fieldsplit_%d_fields 1235 . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields) 1236 . -pc_fieldsplit_type <additive,multiplicative,schur,symmetric_multiplicative> 1237 . -pc_fieldsplit_schur_precondition <true,false> default is true 1238 1239 - Options prefix for inner solvers when using Schur complement preconditioner are -fieldsplit_0_ and -fieldsplit_1_ 1240 for all other solvers they are -fieldsplit_%d_ for the dth field, use -fieldsplit_ for all fields 1241 1242 1243 Notes: use PCFieldSplitSetFields() to set fields defined by "strided" entries and PCFieldSplitSetIS() 1244 to define a field by an arbitrary collection of entries. 1245 1246 If no fields are set the default is used. The fields are defined by entries strided by bs, 1247 beginning at 0 then 1, etc to bs-1. The block size can be set with PCFieldSplitSetBlockSize(), 1248 if this is not called the block size defaults to the blocksize of the second matrix passed 1249 to KSPSetOperators()/PCSetOperators(). 1250 1251 Currently for the multiplicative version, the updated residual needed for the next field 1252 solve is computed via a matrix vector product over the entire array. An optimization would be 1253 to update the residual only for the part of the right hand side associated with the next field 1254 solve. (This would involve more MatGetSubMatrix() calls or some other mechanism to compute the 1255 part of the matrix needed to just update part of the residual). 1256 1257 For the Schur complement preconditioner if J = ( A B ) 1258 ( C D ) 1259 the preconditioner is 1260 (I -B inv(A)) ( inv(A) 0 ) (I 0 ) 1261 (0 I ) ( 0 inv(S) ) (-C inv(A) I ) 1262 where the action of inv(A) is applied using the KSP solver with prefix -fieldsplit_0_. The action of 1263 inv(S) is computed using the KSP solver with prefix -schur_. For PCFieldSplitGetKSP() when field number is 1264 0 it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default 1265 D is used to construct a preconditioner for S, use PCFieldSplitSchurPrecondition() to turn on or off this 1266 option. 1267 1268 If only one set of indices (one IS) is provided with PCFieldSplitSetIS() then the complement of that IS 1269 is used automatically for a second block. 1270 1271 Concepts: physics based preconditioners, block preconditioners 1272 1273 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners 1274 PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSchurPrecondition() 1275 M*/ 1276 1277 EXTERN_C_BEGIN 1278 #undef __FUNCT__ 1279 #define __FUNCT__ "PCCreate_FieldSplit" 1280 PetscErrorCode PETSCKSP_DLLEXPORT PCCreate_FieldSplit(PC pc) 1281 { 1282 PetscErrorCode ierr; 1283 PC_FieldSplit *jac; 1284 1285 PetscFunctionBegin; 1286 ierr = PetscNewLog(pc,PC_FieldSplit,&jac);CHKERRQ(ierr); 1287 jac->bs = -1; 1288 jac->nsplits = 0; 1289 jac->type = PC_COMPOSITE_MULTIPLICATIVE; 1290 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */ 1291 jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACTORIZATION_FULL; 1292 1293 pc->data = (void*)jac; 1294 1295 pc->ops->apply = PCApply_FieldSplit; 1296 pc->ops->applytranspose = PCApplyTranspose_FieldSplit; 1297 pc->ops->setup = PCSetUp_FieldSplit; 1298 pc->ops->destroy = PCDestroy_FieldSplit; 1299 pc->ops->setfromoptions = PCSetFromOptions_FieldSplit; 1300 pc->ops->view = PCView_FieldSplit; 1301 pc->ops->applyrichardson = 0; 1302 1303 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit", 1304 PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 1305 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetFields_C","PCFieldSplitSetFields_FieldSplit", 1306 PCFieldSplitSetFields_FieldSplit);CHKERRQ(ierr); 1307 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetIS_C","PCFieldSplitSetIS_FieldSplit", 1308 PCFieldSplitSetIS_FieldSplit);CHKERRQ(ierr); 1309 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetType_C","PCFieldSplitSetType_FieldSplit", 1310 PCFieldSplitSetType_FieldSplit);CHKERRQ(ierr); 1311 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetBlockSize_C","PCFieldSplitSetBlockSize_FieldSplit", 1312 PCFieldSplitSetBlockSize_FieldSplit);CHKERRQ(ierr); 1313 PetscFunctionReturn(0); 1314 } 1315 EXTERN_C_END 1316 1317 1318 /*MC 1319 Block_Preconditioners - PETSc provides support for a variety of "block preconditioners", this provides an 1320 overview of these methods. 1321 1322 Consider the solution to ( A B ) (x_1) = (b_1) 1323 ( C D ) (x_2) (b_2) 1324 1325 Important special cases, the Stokes equation: C = B' and D = 0 (A B) (x_1) = (b_1) 1326 B' 0) (x_2) (b_2) 1327 1328 One of the goals of the PCFieldSplit preconditioner in PETSc is to provide a variety of preconditioners 1329 for this block system. 1330 1331 Consider an additional matrix (Ap Bp) 1332 (Cp Dp) where some or all of the entries may be the same as 1333 in the original matrix (for example Ap == A). 1334 1335 In the following, A^ denotes the approximate application of the inverse of A, possibly using Ap in the 1336 approximation. In PETSc this simply means one has called KSPSetOperators(ksp,A,Ap,...) or KSPSetOperators(ksp,Ap,Ap,...) 1337 1338 Block Jacobi: x_1 = A^ b_1 1339 x_2 = D^ b_2 1340 1341 Lower block Gauss-Seidel: x_1 = A^ b_1 1342 x_2 = D^ (b_2 - C x_1) variant x_2 = D^ (b_2 - Cp x_1) 1343 1344 Symmetric Gauss-Seidel: x_1 = x_1 + A^(b_1 - A x_1 - B x_2) variant x_1 = x_1 + A^(b_1 - Ap x_1 - Bp x_2) 1345 Interestingly this form is not actually a symmetric matrix, the symmetric version is 1346 x_1 = A^(b_1 - B x_2) variant x_1 = A^(b_1 - Bp x_2) 1347 1348 Level: intermediate 1349 1350 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCFIELDSPLIT 1351 PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSchurPrecondition() 1352 M*/ 1353