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