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,PETSC_OWN_POINTER,&ilink->is);CHKERRQ(ierr); 294 } else { 295 ierr = ISCreateStride(((PetscObject)pc)->comm,nslots,rstart+ilink->fields[0],bs,&ilink->is);CHKERRQ(ierr); 296 } 297 } 298 ierr = ISSorted(ilink->is,&sorted);CHKERRQ(ierr); 299 if (!sorted) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Fields must be sorted when creating split"); 300 ilink = ilink->next; 301 } 302 } 303 304 ilink = jac->head; 305 if (!jac->pmat) { 306 ierr = PetscMalloc(nsplit*sizeof(Mat),&jac->pmat);CHKERRQ(ierr); 307 for (i=0; i<nsplit; i++) { 308 ierr = MatGetSubMatrix(pc->pmat,ilink->is,ilink->is,MAT_INITIAL_MATRIX,&jac->pmat[i]);CHKERRQ(ierr); 309 ilink = ilink->next; 310 } 311 } else { 312 for (i=0; i<nsplit; i++) { 313 ierr = MatGetSubMatrix(pc->pmat,ilink->is,ilink->is,MAT_REUSE_MATRIX,&jac->pmat[i]);CHKERRQ(ierr); 314 ilink = ilink->next; 315 } 316 } 317 if (jac->realdiagonal) { 318 ilink = jac->head; 319 if (!jac->mat) { 320 ierr = PetscMalloc(nsplit*sizeof(Mat),&jac->mat);CHKERRQ(ierr); 321 for (i=0; i<nsplit; i++) { 322 ierr = MatGetSubMatrix(pc->mat,ilink->is,ilink->is,MAT_INITIAL_MATRIX,&jac->mat[i]);CHKERRQ(ierr); 323 ilink = ilink->next; 324 } 325 } else { 326 for (i=0; i<nsplit; i++) { 327 ierr = MatGetSubMatrix(pc->mat,ilink->is,ilink->is,MAT_REUSE_MATRIX,&jac->mat[i]);CHKERRQ(ierr); 328 ilink = ilink->next; 329 } 330 } 331 } else { 332 jac->mat = jac->pmat; 333 } 334 335 if (jac->type != PC_COMPOSITE_ADDITIVE && jac->type != PC_COMPOSITE_SCHUR) { 336 /* extract the rows of the matrix associated with each field: used for efficient computation of residual inside algorithm */ 337 ilink = jac->head; 338 if (!jac->Afield) { 339 ierr = PetscMalloc(nsplit*sizeof(Mat),&jac->Afield);CHKERRQ(ierr); 340 for (i=0; i<nsplit; i++) { 341 ierr = MatGetSubMatrix(pc->mat,ilink->is,PETSC_NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);CHKERRQ(ierr); 342 ilink = ilink->next; 343 } 344 } else { 345 for (i=0; i<nsplit; i++) { 346 ierr = MatGetSubMatrix(pc->mat,ilink->is,PETSC_NULL,MAT_REUSE_MATRIX,&jac->Afield[i]);CHKERRQ(ierr); 347 ilink = ilink->next; 348 } 349 } 350 } 351 352 if (jac->type == PC_COMPOSITE_SCHUR) { 353 IS ccis; 354 PetscInt rstart,rend; 355 if (nsplit != 2) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_INCOMP,"To use Schur complement preconditioner you must have exactly 2 fields"); 356 357 /* When extracting off-diagonal submatrices, we take complements from this range */ 358 ierr = MatGetOwnershipRangeColumn(pc->mat,&rstart,&rend);CHKERRQ(ierr); 359 360 /* need to handle case when one is resetting up the preconditioner */ 361 if (jac->schur) { 362 ilink = jac->head; 363 ierr = ISComplement(ilink->is,rstart,rend,&ccis);CHKERRQ(ierr); 364 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->B);CHKERRQ(ierr); 365 ierr = ISDestroy(ccis);CHKERRQ(ierr); 366 ilink = ilink->next; 367 ierr = ISComplement(ilink->is,rstart,rend,&ccis);CHKERRQ(ierr); 368 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->C);CHKERRQ(ierr); 369 ierr = ISDestroy(ccis);CHKERRQ(ierr); 370 ierr = MatSchurComplementUpdate(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->pmat[1],pc->flag);CHKERRQ(ierr); 371 ierr = KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac),pc->flag);CHKERRQ(ierr); 372 373 } else { 374 KSP ksp; 375 char schurprefix[256]; 376 377 /* extract the B and C matrices */ 378 ilink = jac->head; 379 ierr = ISComplement(ilink->is,rstart,rend,&ccis);CHKERRQ(ierr); 380 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);CHKERRQ(ierr); 381 ierr = ISDestroy(ccis);CHKERRQ(ierr); 382 ilink = ilink->next; 383 ierr = ISComplement(ilink->is,rstart,rend,&ccis);CHKERRQ(ierr); 384 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);CHKERRQ(ierr); 385 ierr = ISDestroy(ccis);CHKERRQ(ierr); 386 /* Better would be to use 'mat[0]' (diagonal block of the real matrix) preconditioned by pmat[0] */ 387 ierr = MatCreateSchurComplement(jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1],&jac->schur);CHKERRQ(ierr); 388 ierr = MatSchurComplementGetKSP(jac->schur,&ksp);CHKERRQ(ierr); 389 ierr = PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,2);CHKERRQ(ierr); 390 ierr = MatSetFromOptions(jac->schur);CHKERRQ(ierr); 391 392 ierr = KSPCreate(((PetscObject)pc)->comm,&jac->kspschur);CHKERRQ(ierr); 393 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)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 PetscBool 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 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 796 797 ierr = PetscSNPrintf(prefix,sizeof prefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",splitname);CHKERRQ(ierr); 798 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 799 800 if (!next) { 801 jac->head = ilink; 802 ilink->previous = PETSC_NULL; 803 } else { 804 while (next->next) { 805 next = next->next; 806 } 807 next->next = ilink; 808 ilink->previous = next; 809 } 810 jac->nsplits++; 811 PetscFunctionReturn(0); 812 } 813 EXTERN_C_END 814 815 EXTERN_C_BEGIN 816 #undef __FUNCT__ 817 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit_Schur" 818 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp) 819 { 820 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 821 PetscErrorCode ierr; 822 823 PetscFunctionBegin; 824 ierr = PetscMalloc(jac->nsplits*sizeof(KSP),subksp);CHKERRQ(ierr); 825 ierr = MatSchurComplementGetKSP(jac->schur,*subksp);CHKERRQ(ierr); 826 (*subksp)[1] = jac->kspschur; 827 *n = jac->nsplits; 828 PetscFunctionReturn(0); 829 } 830 EXTERN_C_END 831 832 EXTERN_C_BEGIN 833 #undef __FUNCT__ 834 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit" 835 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp) 836 { 837 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 838 PetscErrorCode ierr; 839 PetscInt cnt = 0; 840 PC_FieldSplitLink ilink = jac->head; 841 842 PetscFunctionBegin; 843 ierr = PetscMalloc(jac->nsplits*sizeof(KSP*),subksp);CHKERRQ(ierr); 844 while (ilink) { 845 (*subksp)[cnt++] = ilink->ksp; 846 ilink = ilink->next; 847 } 848 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); 849 *n = jac->nsplits; 850 PetscFunctionReturn(0); 851 } 852 EXTERN_C_END 853 854 EXTERN_C_BEGIN 855 #undef __FUNCT__ 856 #define __FUNCT__ "PCFieldSplitSetIS_FieldSplit" 857 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is) 858 { 859 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 860 PetscErrorCode ierr; 861 PC_FieldSplitLink ilink, next = jac->head; 862 char prefix[128]; 863 864 PetscFunctionBegin; 865 if (jac->splitdefined) { 866 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 867 PetscFunctionReturn(0); 868 } 869 ierr = PetscNew(struct _PC_FieldSplitLink,&ilink);CHKERRQ(ierr); 870 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 871 ilink->is = is; 872 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 873 ilink->next = PETSC_NULL; 874 ierr = KSPCreate(((PetscObject)pc)->comm,&ilink->ksp);CHKERRQ(ierr); 875 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 876 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 877 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 878 879 ierr = PetscSNPrintf(prefix,sizeof prefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",splitname);CHKERRQ(ierr); 880 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 881 882 if (!next) { 883 jac->head = ilink; 884 ilink->previous = PETSC_NULL; 885 } else { 886 while (next->next) { 887 next = next->next; 888 } 889 next->next = ilink; 890 ilink->previous = next; 891 } 892 jac->nsplits++; 893 894 PetscFunctionReturn(0); 895 } 896 EXTERN_C_END 897 898 #undef __FUNCT__ 899 #define __FUNCT__ "PCFieldSplitSetFields" 900 /*@ 901 PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner 902 903 Logically Collective on PC 904 905 Input Parameters: 906 + pc - the preconditioner context 907 . splitname - name of this split 908 . n - the number of fields in this split 909 - fields - the fields in this split 910 911 Level: intermediate 912 913 Notes: Use PCFieldSplitSetIS() to set a completely general set of indices as a field. 914 915 The PCFieldSplitSetFields() is for defining fields as a strided blocks. For example, if the block 916 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 917 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x.... 918 where the numbered entries indicate what is in the field. 919 920 This function is called once per split (it creates a new split each time). Solve options 921 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 922 923 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize(), PCFieldSplitSetIS() 924 925 @*/ 926 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields) 927 { 928 PetscErrorCode ierr,(*f)(PC,const char[],PetscInt,const PetscInt *); 929 930 PetscFunctionBegin; 931 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 932 PetscValidCharPointer(splitname,2); 933 if (n < 1) SETERRQ2(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Provided number of fields %D in split \"%s\" not positive",n,splitname); 934 PetscValidIntPointer(fields,3); 935 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCFieldSplitSetFields_C",(void (**)(void))&f);CHKERRQ(ierr); 936 if (f) { 937 ierr = (*f)(pc,splitname,n,fields);CHKERRQ(ierr); 938 } 939 PetscFunctionReturn(0); 940 } 941 942 #undef __FUNCT__ 943 #define __FUNCT__ "PCFieldSplitSetIS" 944 /*@ 945 PCFieldSplitSetIS - Sets the exact elements for field 946 947 Logically Collective on PC 948 949 Input Parameters: 950 + pc - the preconditioner context 951 . splitname - name of this split 952 - is - the index set that defines the vector elements in this field 953 954 955 Notes: 956 Use PCFieldSplitSetFields(), for fields defined by strided types. 957 958 This function is called once per split (it creates a new split each time). Solve options 959 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 960 961 Level: intermediate 962 963 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize() 964 965 @*/ 966 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetIS(PC pc,const char splitname[],IS is) 967 { 968 PetscErrorCode ierr,(*f)(PC,const char[],IS); 969 970 PetscFunctionBegin; 971 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 972 PetscValidCharPointer(splitname,2); 973 PetscValidHeaderSpecific(is,IS_CLASSID,3); 974 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCFieldSplitSetIS_C",(void (**)(void))&f);CHKERRQ(ierr); 975 if (f) { 976 ierr = (*f)(pc,splitname,is);CHKERRQ(ierr); 977 } 978 PetscFunctionReturn(0); 979 } 980 981 #undef __FUNCT__ 982 #define __FUNCT__ "PCFieldSplitSetBlockSize" 983 /*@ 984 PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the 985 fieldsplit preconditioner. If not set the matrix block size is used. 986 987 Logically Collective on PC 988 989 Input Parameters: 990 + pc - the preconditioner context 991 - bs - the block size 992 993 Level: intermediate 994 995 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields() 996 997 @*/ 998 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetBlockSize(PC pc,PetscInt bs) 999 { 1000 PetscErrorCode ierr,(*f)(PC,PetscInt); 1001 1002 PetscFunctionBegin; 1003 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1004 PetscValidLogicalCollectiveInt(pc,bs,2); 1005 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",(void (**)(void))&f);CHKERRQ(ierr); 1006 if (f) { 1007 ierr = (*f)(pc,bs);CHKERRQ(ierr); 1008 } 1009 PetscFunctionReturn(0); 1010 } 1011 1012 #undef __FUNCT__ 1013 #define __FUNCT__ "PCFieldSplitGetSubKSP" 1014 /*@C 1015 PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits 1016 1017 Collective on KSP 1018 1019 Input Parameter: 1020 . pc - the preconditioner context 1021 1022 Output Parameters: 1023 + n - the number of split 1024 - pc - the array of KSP contexts 1025 1026 Note: 1027 After PCFieldSplitGetSubKSP() the array of KSPs IS to be freed by the user 1028 (not the KSP just the array that contains them). 1029 1030 You must call KSPSetUp() before calling PCFieldSplitGetSubKSP(). 1031 1032 Level: advanced 1033 1034 .seealso: PCFIELDSPLIT 1035 @*/ 1036 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[]) 1037 { 1038 PetscErrorCode ierr,(*f)(PC,PetscInt*,KSP **); 1039 1040 PetscFunctionBegin; 1041 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1042 PetscValidIntPointer(n,2); 1043 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",(void (**)(void))&f);CHKERRQ(ierr); 1044 if (f) { 1045 ierr = (*f)(pc,n,subksp);CHKERRQ(ierr); 1046 } else SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONG,"Cannot get subksp for this type of PC"); 1047 PetscFunctionReturn(0); 1048 } 1049 1050 #undef __FUNCT__ 1051 #define __FUNCT__ "PCFieldSplitSchurPrecondition" 1052 /*@ 1053 PCFieldSplitSchurPrecondition - Indicates if the Schur complement is preconditioned by a preconditioner constructed by the 1054 D matrix. Otherwise no preconditioner is used. 1055 1056 Collective on PC 1057 1058 Input Parameters: 1059 + pc - the preconditioner context 1060 . ptype - which matrix to use for preconditioning the Schur complement 1061 - userpre - matrix to use for preconditioning, or PETSC_NULL 1062 1063 Notes: 1064 The default is to use the block on the diagonal of the preconditioning matrix. This is D, in the (1,1) position. 1065 There are currently no preconditioners that work directly with the Schur complement so setting 1066 PC_FIELDSPLIT_SCHUR_PRE_SELF is observationally equivalent to -fieldsplit_1_pc_type none. 1067 1068 Options Database: 1069 . -pc_fieldsplit_schur_precondition <self,user,diag> default is diag 1070 1071 Level: intermediate 1072 1073 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType 1074 1075 @*/ 1076 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1077 { 1078 PetscErrorCode ierr,(*f)(PC,PCFieldSplitSchurPreType,Mat); 1079 1080 PetscFunctionBegin; 1081 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1082 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCFieldSplitSchurPrecondition_C",(void (**)(void))&f);CHKERRQ(ierr); 1083 if (f) { 1084 ierr = (*f)(pc,ptype,pre);CHKERRQ(ierr); 1085 } 1086 PetscFunctionReturn(0); 1087 } 1088 1089 EXTERN_C_BEGIN 1090 #undef __FUNCT__ 1091 #define __FUNCT__ "PCFieldSplitSchurPrecondition_FieldSplit" 1092 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSchurPrecondition_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1093 { 1094 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1095 PetscErrorCode ierr; 1096 1097 PetscFunctionBegin; 1098 jac->schurpre = ptype; 1099 if (pre) { 1100 if (jac->schur_user) {ierr = MatDestroy(jac->schur_user);CHKERRQ(ierr);} 1101 jac->schur_user = pre; 1102 ierr = PetscObjectReference((PetscObject)jac->schur_user);CHKERRQ(ierr); 1103 } 1104 PetscFunctionReturn(0); 1105 } 1106 EXTERN_C_END 1107 1108 #undef __FUNCT__ 1109 #define __FUNCT__ "PCFieldSplitGetSchurBlocks" 1110 /*@C 1111 PCFieldSplitGetSchurBlocks - Gets the all matrix blocks for the Schur complement 1112 1113 Collective on KSP 1114 1115 Input Parameter: 1116 . pc - the preconditioner context 1117 1118 Output Parameters: 1119 + A - the (0,0) block 1120 . B - the (0,1) block 1121 . C - the (1,0) block 1122 - D - the (1,1) block 1123 1124 Level: advanced 1125 1126 .seealso: PCFIELDSPLIT 1127 @*/ 1128 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitGetSchurBlocks(PC pc,Mat *A,Mat *B,Mat *C, Mat *D) 1129 { 1130 PC_FieldSplit *jac = (PC_FieldSplit *) pc->data; 1131 1132 PetscFunctionBegin; 1133 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1134 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach."); 1135 if (A) *A = jac->pmat[0]; 1136 if (B) *B = jac->B; 1137 if (C) *C = jac->C; 1138 if (D) *D = jac->pmat[1]; 1139 PetscFunctionReturn(0); 1140 } 1141 1142 EXTERN_C_BEGIN 1143 #undef __FUNCT__ 1144 #define __FUNCT__ "PCFieldSplitSetType_FieldSplit" 1145 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type) 1146 { 1147 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1148 PetscErrorCode ierr; 1149 1150 PetscFunctionBegin; 1151 jac->type = type; 1152 if (type == PC_COMPOSITE_SCHUR) { 1153 pc->ops->apply = PCApply_FieldSplit_Schur; 1154 pc->ops->view = PCView_FieldSplit_Schur; 1155 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit_Schur",PCFieldSplitGetSubKSP_FieldSplit_Schur);CHKERRQ(ierr); 1156 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSchurPrecondition_C","PCFieldSplitSchurPrecondition_FieldSplit",PCFieldSplitSchurPrecondition_FieldSplit);CHKERRQ(ierr); 1157 1158 } else { 1159 pc->ops->apply = PCApply_FieldSplit; 1160 pc->ops->view = PCView_FieldSplit; 1161 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit",PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 1162 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSchurPrecondition_C","",0);CHKERRQ(ierr); 1163 } 1164 PetscFunctionReturn(0); 1165 } 1166 EXTERN_C_END 1167 1168 EXTERN_C_BEGIN 1169 #undef __FUNCT__ 1170 #define __FUNCT__ "PCFieldSplitSetBlockSize_FieldSplit" 1171 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs) 1172 { 1173 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1174 1175 PetscFunctionBegin; 1176 if (bs < 1) SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs); 1177 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); 1178 jac->bs = bs; 1179 PetscFunctionReturn(0); 1180 } 1181 EXTERN_C_END 1182 1183 #undef __FUNCT__ 1184 #define __FUNCT__ "PCFieldSplitSetType" 1185 /*@ 1186 PCFieldSplitSetType - Sets the type of fieldsplit preconditioner. 1187 1188 Collective on PC 1189 1190 Input Parameter: 1191 . pc - the preconditioner context 1192 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 1193 1194 Options Database Key: 1195 . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type 1196 1197 Level: Developer 1198 1199 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 1200 1201 .seealso: PCCompositeSetType() 1202 1203 @*/ 1204 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetType(PC pc,PCCompositeType type) 1205 { 1206 PetscErrorCode ierr,(*f)(PC,PCCompositeType); 1207 1208 PetscFunctionBegin; 1209 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1210 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCFieldSplitSetType_C",(void (**)(void))&f);CHKERRQ(ierr); 1211 if (f) { 1212 ierr = (*f)(pc,type);CHKERRQ(ierr); 1213 } 1214 PetscFunctionReturn(0); 1215 } 1216 1217 /* -------------------------------------------------------------------------------------*/ 1218 /*MC 1219 PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual 1220 fields or groups of fields 1221 1222 1223 To set options on the solvers for each block append -fieldsplit_ to all the PC 1224 options database keys. For example, -fieldsplit_pc_type ilu -fieldsplit_pc_factor_levels 1 1225 1226 To set the options on the solvers separate for each block call PCFieldSplitGetSubKSP() 1227 and set the options directly on the resulting KSP object 1228 1229 Level: intermediate 1230 1231 Options Database Keys: 1232 + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split 1233 . -pc_fieldsplit_default - automatically add any fields to additional splits that have not 1234 been supplied explicitly by -pc_fieldsplit_%d_fields 1235 . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields) 1236 . -pc_fieldsplit_type <additive,multiplicative,schur,symmetric_multiplicative> 1237 . -pc_fieldsplit_schur_precondition <true,false> default is true 1238 1239 - Options prefix for inner solvers when using Schur complement preconditioner are -fieldsplit_0_ and -fieldsplit_1_ 1240 for all other solvers they are -fieldsplit_%d_ for the dth field, use -fieldsplit_ for all fields 1241 1242 1243 Notes: use PCFieldSplitSetFields() to set fields defined by "strided" entries and PCFieldSplitSetIS() 1244 to define a field by an arbitrary collection of entries. 1245 1246 If no fields are set the default is used. The fields are defined by entries strided by bs, 1247 beginning at 0 then 1, etc to bs-1. The block size can be set with PCFieldSplitSetBlockSize(), 1248 if this is not called the block size defaults to the blocksize of the second matrix passed 1249 to KSPSetOperators()/PCSetOperators(). 1250 1251 Currently for the multiplicative version, the updated residual needed for the next field 1252 solve is computed via a matrix vector product over the entire array. An optimization would be 1253 to update the residual only for the part of the right hand side associated with the next field 1254 solve. (This would involve more MatGetSubMatrix() calls or some other mechanism to compute the 1255 part of the matrix needed to just update part of the residual). 1256 1257 For the Schur complement preconditioner if J = ( A B ) 1258 ( C D ) 1259 the preconditioner is 1260 (I -B inv(A)) ( inv(A) 0 ) (I 0 ) 1261 (0 I ) ( 0 inv(S) ) (-C inv(A) I ) 1262 where the action of inv(A) is applied using the KSP solver with prefix -fieldsplit_0_. The action of 1263 inv(S) is computed using the KSP solver with prefix -schur_. For PCFieldSplitGetKSP() when field number is 1264 0 it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default 1265 D is used to construct a preconditioner for S, use PCFieldSplitSchurPrecondition() to turn on or off this 1266 option. 1267 1268 If only one set of indices (one IS) is provided with PCFieldSplitSetIS() then the complement of that IS 1269 is used automatically for a second block. 1270 1271 Concepts: physics based preconditioners, block preconditioners 1272 1273 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners 1274 PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSchurPrecondition() 1275 M*/ 1276 1277 EXTERN_C_BEGIN 1278 #undef __FUNCT__ 1279 #define __FUNCT__ "PCCreate_FieldSplit" 1280 PetscErrorCode PETSCKSP_DLLEXPORT PCCreate_FieldSplit(PC pc) 1281 { 1282 PetscErrorCode ierr; 1283 PC_FieldSplit *jac; 1284 1285 PetscFunctionBegin; 1286 ierr = PetscNewLog(pc,PC_FieldSplit,&jac);CHKERRQ(ierr); 1287 jac->bs = -1; 1288 jac->nsplits = 0; 1289 jac->type = PC_COMPOSITE_MULTIPLICATIVE; 1290 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */ 1291 jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACTORIZATION_FULL; 1292 1293 pc->data = (void*)jac; 1294 1295 pc->ops->apply = PCApply_FieldSplit; 1296 pc->ops->applytranspose = PCApplyTranspose_FieldSplit; 1297 pc->ops->setup = PCSetUp_FieldSplit; 1298 pc->ops->destroy = PCDestroy_FieldSplit; 1299 pc->ops->setfromoptions = PCSetFromOptions_FieldSplit; 1300 pc->ops->view = PCView_FieldSplit; 1301 pc->ops->applyrichardson = 0; 1302 1303 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit", 1304 PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 1305 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetFields_C","PCFieldSplitSetFields_FieldSplit", 1306 PCFieldSplitSetFields_FieldSplit);CHKERRQ(ierr); 1307 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetIS_C","PCFieldSplitSetIS_FieldSplit", 1308 PCFieldSplitSetIS_FieldSplit);CHKERRQ(ierr); 1309 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetType_C","PCFieldSplitSetType_FieldSplit", 1310 PCFieldSplitSetType_FieldSplit);CHKERRQ(ierr); 1311 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetBlockSize_C","PCFieldSplitSetBlockSize_FieldSplit", 1312 PCFieldSplitSetBlockSize_FieldSplit);CHKERRQ(ierr); 1313 PetscFunctionReturn(0); 1314 } 1315 EXTERN_C_END 1316 1317 1318 /*MC 1319 Block_Preconditioners - PETSc provides support for a variety of "block preconditioners", this provides an 1320 overview of these methods. 1321 1322 Consider the solution to ( A B ) (x_1) = (b_1) 1323 ( C D ) (x_2) (b_2) 1324 1325 Important special cases, the Stokes equation: C = B' and D = 0 (A B) (x_1) = (b_1) 1326 B' 0) (x_2) (b_2) 1327 1328 One of the goals of the PCFieldSplit preconditioner in PETSc is to provide a variety of preconditioners 1329 for this block system. 1330 1331 Consider an additional matrix (Ap Bp) 1332 (Cp Dp) where some or all of the entries may be the same as 1333 in the original matrix (for example Ap == A). 1334 1335 In the following, A^ denotes the approximate application of the inverse of A, possibly using Ap in the 1336 approximation. In PETSc this simply means one has called KSPSetOperators(ksp,A,Ap,...) or KSPSetOperators(ksp,Ap,Ap,...) 1337 1338 Block Jacobi: x_1 = A^ b_1 1339 x_2 = D^ b_2 1340 1341 Lower block Gauss-Seidel: x_1 = A^ b_1 1342 x_2 = D^ (b_2 - C x_1) variant x_2 = D^ (b_2 - Cp x_1) 1343 1344 Symmetric Gauss-Seidel: x_1 = x_1 + A^(b_1 - A x_1 - B x_2) variant x_1 = x_1 + A^(b_1 - Ap x_1 - Bp x_2) 1345 Interestingly this form is not actually a symmetric matrix, the symmetric version is 1346 x_1 = A^(b_1 - B x_2) variant x_1 = A^(b_1 - Bp x_2) 1347 1348 Level: intermediate 1349 1350 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCFIELDSPLIT 1351 PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSchurPrecondition() 1352 M*/ 1353