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