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