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