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 PetscValidLogicalCollectiveInt(pc,bs,2); 1003 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",(void (**)(void))&f);CHKERRQ(ierr); 1004 if (f) { 1005 ierr = (*f)(pc,bs);CHKERRQ(ierr); 1006 } 1007 PetscFunctionReturn(0); 1008 } 1009 1010 #undef __FUNCT__ 1011 #define __FUNCT__ "PCFieldSplitGetSubKSP" 1012 /*@C 1013 PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits 1014 1015 Collective on KSP 1016 1017 Input Parameter: 1018 . pc - the preconditioner context 1019 1020 Output Parameters: 1021 + n - the number of split 1022 - pc - the array of KSP contexts 1023 1024 Note: 1025 After PCFieldSplitGetSubKSP() the array of KSPs IS to be freed by the user 1026 (not the KSP just the array that contains them). 1027 1028 You must call KSPSetUp() before calling PCFieldSplitGetSubKSP(). 1029 1030 Level: advanced 1031 1032 .seealso: PCFIELDSPLIT 1033 @*/ 1034 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[]) 1035 { 1036 PetscErrorCode ierr,(*f)(PC,PetscInt*,KSP **); 1037 1038 PetscFunctionBegin; 1039 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1040 PetscValidIntPointer(n,2); 1041 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",(void (**)(void))&f);CHKERRQ(ierr); 1042 if (f) { 1043 ierr = (*f)(pc,n,subksp);CHKERRQ(ierr); 1044 } else SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONG,"Cannot get subksp for this type of PC"); 1045 PetscFunctionReturn(0); 1046 } 1047 1048 #undef __FUNCT__ 1049 #define __FUNCT__ "PCFieldSplitSchurPrecondition" 1050 /*@ 1051 PCFieldSplitSchurPrecondition - Indicates if the Schur complement is preconditioned by a preconditioner constructed by the 1052 D matrix. Otherwise no preconditioner is used. 1053 1054 Collective on PC 1055 1056 Input Parameters: 1057 + pc - the preconditioner context 1058 . ptype - which matrix to use for preconditioning the Schur complement 1059 - userpre - matrix to use for preconditioning, or PETSC_NULL 1060 1061 Notes: 1062 The default is to use the block on the diagonal of the preconditioning matrix. This is D, in the (1,1) position. 1063 There are currently no preconditioners that work directly with the Schur complement so setting 1064 PC_FIELDSPLIT_SCHUR_PRE_SELF is observationally equivalent to -fieldsplit_1_pc_type none. 1065 1066 Options Database: 1067 . -pc_fieldsplit_schur_precondition <self,user,diag> default is diag 1068 1069 Level: intermediate 1070 1071 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType 1072 1073 @*/ 1074 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1075 { 1076 PetscErrorCode ierr,(*f)(PC,PCFieldSplitSchurPreType,Mat); 1077 1078 PetscFunctionBegin; 1079 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1080 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCFieldSplitSchurPrecondition_C",(void (**)(void))&f);CHKERRQ(ierr); 1081 if (f) { 1082 ierr = (*f)(pc,ptype,pre);CHKERRQ(ierr); 1083 } 1084 PetscFunctionReturn(0); 1085 } 1086 1087 EXTERN_C_BEGIN 1088 #undef __FUNCT__ 1089 #define __FUNCT__ "PCFieldSplitSchurPrecondition_FieldSplit" 1090 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSchurPrecondition_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1091 { 1092 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1093 PetscErrorCode ierr; 1094 1095 PetscFunctionBegin; 1096 jac->schurpre = ptype; 1097 if (pre) { 1098 if (jac->schur_user) {ierr = MatDestroy(jac->schur_user);CHKERRQ(ierr);} 1099 jac->schur_user = pre; 1100 ierr = PetscObjectReference((PetscObject)jac->schur_user);CHKERRQ(ierr); 1101 } 1102 PetscFunctionReturn(0); 1103 } 1104 EXTERN_C_END 1105 1106 #undef __FUNCT__ 1107 #define __FUNCT__ "PCFieldSplitGetSchurBlocks" 1108 /*@C 1109 PCFieldSplitGetSchurBlocks - Gets the all matrix blocks for the Schur complement 1110 1111 Collective on KSP 1112 1113 Input Parameter: 1114 . pc - the preconditioner context 1115 1116 Output Parameters: 1117 + A - the (0,0) block 1118 . B - the (0,1) block 1119 . C - the (1,0) block 1120 - D - the (1,1) block 1121 1122 Level: advanced 1123 1124 .seealso: PCFIELDSPLIT 1125 @*/ 1126 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitGetSchurBlocks(PC pc,Mat *A,Mat *B,Mat *C, Mat *D) 1127 { 1128 PC_FieldSplit *jac = (PC_FieldSplit *) pc->data; 1129 1130 PetscFunctionBegin; 1131 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1132 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach."); 1133 if (A) *A = jac->pmat[0]; 1134 if (B) *B = jac->B; 1135 if (C) *C = jac->C; 1136 if (D) *D = jac->pmat[1]; 1137 PetscFunctionReturn(0); 1138 } 1139 1140 EXTERN_C_BEGIN 1141 #undef __FUNCT__ 1142 #define __FUNCT__ "PCFieldSplitSetType_FieldSplit" 1143 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type) 1144 { 1145 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1146 PetscErrorCode ierr; 1147 1148 PetscFunctionBegin; 1149 jac->type = type; 1150 if (type == PC_COMPOSITE_SCHUR) { 1151 pc->ops->apply = PCApply_FieldSplit_Schur; 1152 pc->ops->view = PCView_FieldSplit_Schur; 1153 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit_Schur",PCFieldSplitGetSubKSP_FieldSplit_Schur);CHKERRQ(ierr); 1154 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSchurPrecondition_C","PCFieldSplitSchurPrecondition_FieldSplit",PCFieldSplitSchurPrecondition_FieldSplit);CHKERRQ(ierr); 1155 1156 } else { 1157 pc->ops->apply = PCApply_FieldSplit; 1158 pc->ops->view = PCView_FieldSplit; 1159 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit",PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 1160 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSchurPrecondition_C","",0);CHKERRQ(ierr); 1161 } 1162 PetscFunctionReturn(0); 1163 } 1164 EXTERN_C_END 1165 1166 EXTERN_C_BEGIN 1167 #undef __FUNCT__ 1168 #define __FUNCT__ "PCFieldSplitSetBlockSize_FieldSplit" 1169 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs) 1170 { 1171 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1172 1173 PetscFunctionBegin; 1174 if (bs < 1) SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs); 1175 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); 1176 jac->bs = bs; 1177 PetscFunctionReturn(0); 1178 } 1179 EXTERN_C_END 1180 1181 #undef __FUNCT__ 1182 #define __FUNCT__ "PCFieldSplitSetType" 1183 /*@ 1184 PCFieldSplitSetType - Sets the type of fieldsplit preconditioner. 1185 1186 Collective on PC 1187 1188 Input Parameter: 1189 . pc - the preconditioner context 1190 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 1191 1192 Options Database Key: 1193 . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type 1194 1195 Level: Developer 1196 1197 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 1198 1199 .seealso: PCCompositeSetType() 1200 1201 @*/ 1202 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetType(PC pc,PCCompositeType type) 1203 { 1204 PetscErrorCode ierr,(*f)(PC,PCCompositeType); 1205 1206 PetscFunctionBegin; 1207 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1208 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCFieldSplitSetType_C",(void (**)(void))&f);CHKERRQ(ierr); 1209 if (f) { 1210 ierr = (*f)(pc,type);CHKERRQ(ierr); 1211 } 1212 PetscFunctionReturn(0); 1213 } 1214 1215 /* -------------------------------------------------------------------------------------*/ 1216 /*MC 1217 PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual 1218 fields or groups of fields 1219 1220 1221 To set options on the solvers for each block append -fieldsplit_ to all the PC 1222 options database keys. For example, -fieldsplit_pc_type ilu -fieldsplit_pc_factor_levels 1 1223 1224 To set the options on the solvers separate for each block call PCFieldSplitGetSubKSP() 1225 and set the options directly on the resulting KSP object 1226 1227 Level: intermediate 1228 1229 Options Database Keys: 1230 + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split 1231 . -pc_fieldsplit_default - automatically add any fields to additional splits that have not 1232 been supplied explicitly by -pc_fieldsplit_%d_fields 1233 . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields) 1234 . -pc_fieldsplit_type <additive,multiplicative,schur,symmetric_multiplicative> 1235 . -pc_fieldsplit_schur_precondition <true,false> default is true 1236 1237 - Options prefix for inner solvers when using Schur complement preconditioner are -fieldsplit_0_ and -fieldsplit_1_ 1238 for all other solvers they are -fieldsplit_%d_ for the dth field, use -fieldsplit_ for all fields 1239 1240 1241 Notes: use PCFieldSplitSetFields() to set fields defined by "strided" entries and PCFieldSplitSetIS() 1242 to define a field by an arbitrary collection of entries. 1243 1244 If no fields are set the default is used. The fields are defined by entries strided by bs, 1245 beginning at 0 then 1, etc to bs-1. The block size can be set with PCFieldSplitSetBlockSize(), 1246 if this is not called the block size defaults to the blocksize of the second matrix passed 1247 to KSPSetOperators()/PCSetOperators(). 1248 1249 Currently for the multiplicative version, the updated residual needed for the next field 1250 solve is computed via a matrix vector product over the entire array. An optimization would be 1251 to update the residual only for the part of the right hand side associated with the next field 1252 solve. (This would involve more MatGetSubMatrix() calls or some other mechanism to compute the 1253 part of the matrix needed to just update part of the residual). 1254 1255 For the Schur complement preconditioner if J = ( A B ) 1256 ( C D ) 1257 the preconditioner is 1258 (I -B inv(A)) ( inv(A) 0 ) (I 0 ) 1259 (0 I ) ( 0 inv(S) ) (-C inv(A) I ) 1260 where the action of inv(A) is applied using the KSP solver with prefix -fieldsplit_0_. The action of 1261 inv(S) is computed using the KSP solver with prefix -schur_. For PCFieldSplitGetKSP() when field number is 1262 0 it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default 1263 D is used to construct a preconditioner for S, use PCFieldSplitSchurPrecondition() to turn on or off this 1264 option. 1265 1266 If only one set of indices (one IS) is provided with PCFieldSplitSetIS() then the complement of that IS 1267 is used automatically for a second block. 1268 1269 Concepts: physics based preconditioners, block preconditioners 1270 1271 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners 1272 PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSchurPrecondition() 1273 M*/ 1274 1275 EXTERN_C_BEGIN 1276 #undef __FUNCT__ 1277 #define __FUNCT__ "PCCreate_FieldSplit" 1278 PetscErrorCode PETSCKSP_DLLEXPORT PCCreate_FieldSplit(PC pc) 1279 { 1280 PetscErrorCode ierr; 1281 PC_FieldSplit *jac; 1282 1283 PetscFunctionBegin; 1284 ierr = PetscNewLog(pc,PC_FieldSplit,&jac);CHKERRQ(ierr); 1285 jac->bs = -1; 1286 jac->nsplits = 0; 1287 jac->type = PC_COMPOSITE_MULTIPLICATIVE; 1288 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */ 1289 jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACTORIZATION_FULL; 1290 1291 pc->data = (void*)jac; 1292 1293 pc->ops->apply = PCApply_FieldSplit; 1294 pc->ops->applytranspose = PCApplyTranspose_FieldSplit; 1295 pc->ops->setup = PCSetUp_FieldSplit; 1296 pc->ops->destroy = PCDestroy_FieldSplit; 1297 pc->ops->setfromoptions = PCSetFromOptions_FieldSplit; 1298 pc->ops->view = PCView_FieldSplit; 1299 pc->ops->applyrichardson = 0; 1300 1301 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit", 1302 PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 1303 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetFields_C","PCFieldSplitSetFields_FieldSplit", 1304 PCFieldSplitSetFields_FieldSplit);CHKERRQ(ierr); 1305 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetIS_C","PCFieldSplitSetIS_FieldSplit", 1306 PCFieldSplitSetIS_FieldSplit);CHKERRQ(ierr); 1307 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetType_C","PCFieldSplitSetType_FieldSplit", 1308 PCFieldSplitSetType_FieldSplit);CHKERRQ(ierr); 1309 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetBlockSize_C","PCFieldSplitSetBlockSize_FieldSplit", 1310 PCFieldSplitSetBlockSize_FieldSplit);CHKERRQ(ierr); 1311 PetscFunctionReturn(0); 1312 } 1313 EXTERN_C_END 1314 1315 1316 /*MC 1317 Block_Preconditioners - PETSc provides support for a variety of "block preconditioners", this provides an 1318 overview of these methods. 1319 1320 Consider the solution to ( A B ) (x_1) = (b_1) 1321 ( C D ) (x_2) (b_2) 1322 1323 Important special cases, the Stokes equation: C = B' and D = 0 (A B) (x_1) = (b_1) 1324 B' 0) (x_2) (b_2) 1325 1326 One of the goals of the PCFieldSplit preconditioner in PETSc is to provide a variety of preconditioners 1327 for this block system. 1328 1329 Consider an additional matrix (Ap Bp) 1330 (Cp Dp) where some or all of the entries may be the same as 1331 in the original matrix (for example Ap == A). 1332 1333 In the following, A^ denotes the approximate application of the inverse of A, possibly using Ap in the 1334 approximation. In PETSc this simply means one has called KSPSetOperators(ksp,A,Ap,...) or KSPSetOperators(ksp,Ap,Ap,...) 1335 1336 Block Jacobi: x_1 = A^ b_1 1337 x_2 = D^ b_2 1338 1339 Lower block Gauss-Seidel: x_1 = A^ b_1 1340 x_2 = D^ (b_2 - C x_1) variant x_2 = D^ (b_2 - Cp x_1) 1341 1342 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) 1343 Interestingly this form is not actually a symmetric matrix, the symmetric version is 1344 x_1 = A^(b_1 - B x_2) variant x_1 = A^(b_1 - Bp x_2) 1345 1346 Level: intermediate 1347 1348 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCFIELDSPLIT 1349 PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSchurPrecondition() 1350 M*/ 1351