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