1 2 #include <private/pcimpl.h> /*I "petscpc.h" I*/ 3 #include <petscdmcomposite.h> /*I "petscdmcomposite.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 = A11 - A10 A00^{-1} A01 */ 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 A00 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 = A11 - A10 inv(A00) A01 \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_detect_saddle_point",&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 A01 and A10 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 /* set tabbing and options prefix of KSP inside the MatSchur */ 421 ierr = MatSchurComplementGetKSP(jac->schur,&ksp);CHKERRQ(ierr); 422 ierr = PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,2);CHKERRQ(ierr); 423 ierr = PetscSNPrintf(schurprefix,sizeof schurprefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",jac->head->splitname);CHKERRQ(ierr); 424 ierr = KSPSetOptionsPrefix(ksp,schurprefix);CHKERRQ(ierr); 425 ierr = MatSetFromOptions(jac->schur);CHKERRQ(ierr); 426 427 ierr = KSPCreate(((PetscObject)pc)->comm,&jac->kspschur);CHKERRQ(ierr); 428 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->kspschur);CHKERRQ(ierr); 429 ierr = PetscObjectIncrementTabLevel((PetscObject)jac->kspschur,(PetscObject)pc,1);CHKERRQ(ierr); 430 ierr = KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac),DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 431 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) { 432 PC pc; 433 ierr = KSPGetPC(jac->kspschur,&pc);CHKERRQ(ierr); 434 ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr); 435 /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */ 436 } 437 ierr = PetscSNPrintf(schurprefix,sizeof schurprefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",ilink->splitname);CHKERRQ(ierr); 438 ierr = KSPSetOptionsPrefix(jac->kspschur,schurprefix);CHKERRQ(ierr); 439 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 440 ierr = KSPSetFromOptions(jac->kspschur);CHKERRQ(ierr); 441 442 ierr = PetscMalloc2(2,Vec,&jac->x,2,Vec,&jac->y);CHKERRQ(ierr); 443 ierr = MatGetVecs(jac->pmat[0],&jac->x[0],&jac->y[0]);CHKERRQ(ierr); 444 ierr = MatGetVecs(jac->pmat[1],&jac->x[1],&jac->y[1]);CHKERRQ(ierr); 445 ilink = jac->head; 446 ilink->x = jac->x[0]; ilink->y = jac->y[0]; 447 ilink = ilink->next; 448 ilink->x = jac->x[1]; ilink->y = jac->y[1]; 449 } 450 } else { 451 /* set up the individual PCs */ 452 i = 0; 453 ilink = jac->head; 454 while (ilink) { 455 ierr = KSPSetOperators(ilink->ksp,jac->mat[i],jac->pmat[i],flag);CHKERRQ(ierr); 456 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 457 ierr = KSPSetFromOptions(ilink->ksp);CHKERRQ(ierr); 458 ierr = KSPSetUp(ilink->ksp);CHKERRQ(ierr); 459 i++; 460 ilink = ilink->next; 461 } 462 463 /* create work vectors for each split */ 464 if (!jac->x) { 465 ierr = PetscMalloc2(nsplit,Vec,&jac->x,nsplit,Vec,&jac->y);CHKERRQ(ierr); 466 ilink = jac->head; 467 for (i=0; i<nsplit; i++) { 468 Vec *vl,*vr; 469 470 ierr = KSPGetVecs(ilink->ksp,1,&vr,1,&vl);CHKERRQ(ierr); 471 ilink->x = *vr; 472 ilink->y = *vl; 473 ierr = PetscFree(vr);CHKERRQ(ierr); 474 ierr = PetscFree(vl);CHKERRQ(ierr); 475 jac->x[i] = ilink->x; 476 jac->y[i] = ilink->y; 477 ilink = ilink->next; 478 } 479 } 480 } 481 482 483 if (!jac->head->sctx) { 484 Vec xtmp; 485 486 /* compute scatter contexts needed by multiplicative versions and non-default splits */ 487 488 ilink = jac->head; 489 ierr = MatGetVecs(pc->pmat,&xtmp,PETSC_NULL);CHKERRQ(ierr); 490 for (i=0; i<nsplit; i++) { 491 ierr = VecScatterCreate(xtmp,ilink->is,jac->x[i],PETSC_NULL,&ilink->sctx);CHKERRQ(ierr); 492 ilink = ilink->next; 493 } 494 ierr = VecDestroy(xtmp);CHKERRQ(ierr); 495 } 496 PetscFunctionReturn(0); 497 } 498 499 #define FieldSplitSplitSolveAdd(ilink,xx,yy) \ 500 (VecScatterBegin(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 501 VecScatterEnd(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 502 KSPSolve(ilink->ksp,ilink->x,ilink->y) || \ 503 VecScatterBegin(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE) || \ 504 VecScatterEnd(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE)) 505 506 #undef __FUNCT__ 507 #define __FUNCT__ "PCApply_FieldSplit_Schur" 508 static PetscErrorCode PCApply_FieldSplit_Schur(PC pc,Vec x,Vec y) 509 { 510 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 511 PetscErrorCode ierr; 512 KSP ksp; 513 PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next; 514 515 PetscFunctionBegin; 516 ierr = MatSchurComplementGetKSP(jac->schur,&ksp);CHKERRQ(ierr); 517 518 switch (jac->schurfactorization) { 519 case PC_FIELDSPLIT_SCHUR_FACTORIZATION_DIAG: 520 /* [A00 0; 0 -S], positive definite, suitable for MINRES */ 521 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 522 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 523 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 524 ierr = KSPSolve(ksp,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 525 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 526 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 527 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 528 ierr = VecScale(ilinkD->y,-1.);CHKERRQ(ierr); 529 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 530 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 531 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 532 break; 533 case PC_FIELDSPLIT_SCHUR_FACTORIZATION_LOWER: 534 /* [A00 0; A10 S], suitable for left preconditioning */ 535 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 536 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 537 ierr = KSPSolve(ksp,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 538 ierr = MatMult(jac->C,ilinkA->y,ilinkD->x);CHKERRQ(ierr); 539 ierr = VecScale(ilinkD->x,-1.);CHKERRQ(ierr); 540 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 541 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 542 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 543 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 544 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 545 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 546 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 547 break; 548 case PC_FIELDSPLIT_SCHUR_FACTORIZATION_UPPER: 549 /* [A00 A01; 0 S], suitable for right preconditioning */ 550 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 551 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 552 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 553 ierr = MatMult(jac->B,ilinkD->y,ilinkA->x);CHKERRQ(ierr); 554 ierr = VecScale(ilinkA->x,-1.);CHKERRQ(ierr); 555 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 556 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 557 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 558 ierr = KSPSolve(ksp,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 559 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 560 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 561 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 562 break; 563 case PC_FIELDSPLIT_SCHUR_FACTORIZATION_FULL: 564 /* [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 */ 565 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 566 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 567 ierr = KSPSolve(ksp,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 568 ierr = MatMult(jac->C,ilinkA->y,ilinkD->x);CHKERRQ(ierr); 569 ierr = VecScale(ilinkD->x,-1.0);CHKERRQ(ierr); 570 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 571 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 572 573 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 574 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 575 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 576 577 ierr = MatMult(jac->B,ilinkD->y,ilinkA->y);CHKERRQ(ierr); 578 ierr = VecAXPY(ilinkA->x,-1.0,ilinkA->y);CHKERRQ(ierr); 579 ierr = KSPSolve(ksp,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 580 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 581 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 582 } 583 PetscFunctionReturn(0); 584 } 585 586 #undef __FUNCT__ 587 #define __FUNCT__ "PCApply_FieldSplit" 588 static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y) 589 { 590 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 591 PetscErrorCode ierr; 592 PC_FieldSplitLink ilink = jac->head; 593 PetscInt cnt; 594 595 PetscFunctionBegin; 596 CHKMEMQ; 597 ierr = VecSetBlockSize(x,jac->bs);CHKERRQ(ierr); 598 ierr = VecSetBlockSize(y,jac->bs);CHKERRQ(ierr); 599 600 if (jac->type == PC_COMPOSITE_ADDITIVE) { 601 if (jac->defaultsplit) { 602 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 603 while (ilink) { 604 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 605 ilink = ilink->next; 606 } 607 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 608 } else { 609 ierr = VecSet(y,0.0);CHKERRQ(ierr); 610 while (ilink) { 611 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 612 ilink = ilink->next; 613 } 614 } 615 } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 616 if (!jac->w1) { 617 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 618 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 619 } 620 ierr = VecSet(y,0.0);CHKERRQ(ierr); 621 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 622 cnt = 1; 623 while (ilink->next) { 624 ilink = ilink->next; 625 /* compute the residual only over the part of the vector needed */ 626 ierr = MatMult(jac->Afield[cnt++],y,ilink->x);CHKERRQ(ierr); 627 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 628 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 629 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 630 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 631 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 632 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 633 } 634 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 635 cnt -= 2; 636 while (ilink->previous) { 637 ilink = ilink->previous; 638 /* compute the residual only over the part of the vector needed */ 639 ierr = MatMult(jac->Afield[cnt--],y,ilink->x);CHKERRQ(ierr); 640 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 641 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 642 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 643 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 644 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 645 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 646 } 647 } 648 } else SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Unsupported or unknown composition",(int) jac->type); 649 CHKMEMQ; 650 PetscFunctionReturn(0); 651 } 652 653 #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \ 654 (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 655 VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 656 KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) || \ 657 VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \ 658 VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE)) 659 660 #undef __FUNCT__ 661 #define __FUNCT__ "PCApplyTranspose_FieldSplit" 662 static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y) 663 { 664 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 665 PetscErrorCode ierr; 666 PC_FieldSplitLink ilink = jac->head; 667 668 PetscFunctionBegin; 669 CHKMEMQ; 670 ierr = VecSetBlockSize(x,jac->bs);CHKERRQ(ierr); 671 ierr = VecSetBlockSize(y,jac->bs);CHKERRQ(ierr); 672 673 if (jac->type == PC_COMPOSITE_ADDITIVE) { 674 if (jac->defaultsplit) { 675 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 676 while (ilink) { 677 ierr = KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 678 ilink = ilink->next; 679 } 680 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 681 } else { 682 ierr = VecSet(y,0.0);CHKERRQ(ierr); 683 while (ilink) { 684 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 685 ilink = ilink->next; 686 } 687 } 688 } else { 689 if (!jac->w1) { 690 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 691 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 692 } 693 ierr = VecSet(y,0.0);CHKERRQ(ierr); 694 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 695 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 696 while (ilink->next) { 697 ilink = ilink->next; 698 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 699 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 700 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 701 } 702 while (ilink->previous) { 703 ilink = ilink->previous; 704 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 705 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 706 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 707 } 708 } else { 709 while (ilink->next) { /* get to last entry in linked list */ 710 ilink = ilink->next; 711 } 712 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 713 while (ilink->previous) { 714 ilink = ilink->previous; 715 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 716 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 717 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 718 } 719 } 720 } 721 CHKMEMQ; 722 PetscFunctionReturn(0); 723 } 724 725 #undef __FUNCT__ 726 #define __FUNCT__ "PCReset_FieldSplit" 727 static PetscErrorCode PCReset_FieldSplit(PC pc) 728 { 729 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 730 PetscErrorCode ierr; 731 PC_FieldSplitLink ilink = jac->head,next; 732 733 PetscFunctionBegin; 734 while (ilink) { 735 ierr = KSPReset(ilink->ksp);CHKERRQ(ierr); 736 if (ilink->x) {ierr = VecDestroy(ilink->x);CHKERRQ(ierr);} 737 if (ilink->y) {ierr = VecDestroy(ilink->y);CHKERRQ(ierr);} 738 if (ilink->sctx) {ierr = VecScatterDestroy(ilink->sctx);CHKERRQ(ierr);} 739 if (ilink->is) {ierr = ISDestroy(ilink->is);CHKERRQ(ierr);} 740 next = ilink->next; 741 ilink = next; 742 } 743 ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr); 744 if (jac->mat && jac->mat != jac->pmat) { 745 ierr = MatDestroyMatrices(jac->nsplits,&jac->mat);CHKERRQ(ierr); 746 } else if (jac->mat) { 747 jac->mat = PETSC_NULL; 748 } 749 if (jac->pmat) {ierr = MatDestroyMatrices(jac->nsplits,&jac->pmat);CHKERRQ(ierr);} 750 if (jac->Afield) {ierr = MatDestroyMatrices(jac->nsplits,&jac->Afield);CHKERRQ(ierr);} 751 if (jac->w1) {ierr = VecDestroy(jac->w1);CHKERRQ(ierr);} 752 if (jac->w2) {ierr = VecDestroy(jac->w2);CHKERRQ(ierr);} 753 if (jac->schur) {ierr = MatDestroy(jac->schur);CHKERRQ(ierr);} 754 if (jac->schur_user) {ierr = MatDestroy(jac->schur_user);CHKERRQ(ierr);} 755 if (jac->kspschur) {ierr = KSPDestroy(jac->kspschur);CHKERRQ(ierr);} 756 if (jac->B) {ierr = MatDestroy(jac->B);CHKERRQ(ierr);} 757 if (jac->C) {ierr = MatDestroy(jac->C);CHKERRQ(ierr);} 758 PetscFunctionReturn(0); 759 } 760 761 #undef __FUNCT__ 762 #define __FUNCT__ "PCDestroy_FieldSplit" 763 static PetscErrorCode PCDestroy_FieldSplit(PC pc) 764 { 765 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 766 PetscErrorCode ierr; 767 PC_FieldSplitLink ilink = jac->head,next; 768 769 PetscFunctionBegin; 770 ierr = PCReset_FieldSplit(pc);CHKERRQ(ierr); 771 while (ilink) { 772 ierr = KSPDestroy(ilink->ksp);CHKERRQ(ierr); 773 next = ilink->next; 774 ierr = PetscFree(ilink->splitname);CHKERRQ(ierr); 775 ierr = PetscFree(ilink->fields);CHKERRQ(ierr); 776 ierr = PetscFree(ilink);CHKERRQ(ierr); 777 ilink = next; 778 } 779 ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr); 780 ierr = PetscFree(pc->data);CHKERRQ(ierr); 781 PetscFunctionReturn(0); 782 } 783 784 #undef __FUNCT__ 785 #define __FUNCT__ "PCSetFromOptions_FieldSplit" 786 static PetscErrorCode PCSetFromOptions_FieldSplit(PC pc) 787 { 788 PetscErrorCode ierr; 789 PetscInt bs; 790 PetscBool flg,stokes = PETSC_FALSE; 791 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 792 PCCompositeType ctype; 793 794 PetscFunctionBegin; 795 ierr = PetscOptionsHead("FieldSplit options");CHKERRQ(ierr); 796 ierr = PetscOptionsBool("-pc_fieldsplit_real_diagonal","Use diagonal blocks of the operator","PCFieldSplitSetRealDiagonal",jac->realdiagonal,&jac->realdiagonal,PETSC_NULL);CHKERRQ(ierr); 797 ierr = PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);CHKERRQ(ierr); 798 if (flg) { 799 ierr = PCFieldSplitSetBlockSize(pc,bs);CHKERRQ(ierr); 800 } 801 802 ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,PETSC_NULL);CHKERRQ(ierr); 803 if (stokes) { 804 ierr = PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);CHKERRQ(ierr); 805 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_SELF; 806 } 807 808 ierr = PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg);CHKERRQ(ierr); 809 if (flg) { 810 ierr = PCFieldSplitSetType(pc,ctype);CHKERRQ(ierr); 811 } 812 813 /* Only setup fields once */ 814 if ((jac->bs > 0) && (jac->nsplits == 0)) { 815 /* only allow user to set fields from command line if bs is already known. 816 otherwise user can set them in PCFieldSplitSetDefaults() */ 817 ierr = PCFieldSplitSetRuntimeSplits_Private(pc);CHKERRQ(ierr); 818 if (jac->splitdefined) {ierr = PetscInfo(pc,"Splits defined using the options database\n");CHKERRQ(ierr);} 819 } 820 if (jac->type == PC_COMPOSITE_SCHUR) { 821 ierr = PetscOptionsEnum("-pc_fieldsplit_schur_factorization_type","Factorization to use","None",PCFieldSplitSchurFactorizationTypes,(PetscEnum)jac->schurfactorization,(PetscEnum*)&jac->schurfactorization,PETSC_NULL);CHKERRQ(ierr); 822 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); 823 } 824 ierr = PetscOptionsTail();CHKERRQ(ierr); 825 PetscFunctionReturn(0); 826 } 827 828 /*------------------------------------------------------------------------------------*/ 829 830 EXTERN_C_BEGIN 831 #undef __FUNCT__ 832 #define __FUNCT__ "PCFieldSplitSetFields_FieldSplit" 833 PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields) 834 { 835 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 836 PetscErrorCode ierr; 837 PC_FieldSplitLink ilink,next = jac->head; 838 char prefix[128]; 839 PetscInt i; 840 841 PetscFunctionBegin; 842 if (jac->splitdefined) { 843 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 844 PetscFunctionReturn(0); 845 } 846 for (i=0; i<n; i++) { 847 if (fields[i] >= jac->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs); 848 if (fields[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]); 849 } 850 ierr = PetscNew(struct _PC_FieldSplitLink,&ilink);CHKERRQ(ierr); 851 if (splitname) { 852 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 853 } else { 854 ierr = PetscMalloc(3*sizeof(char),&ilink->splitname);CHKERRQ(ierr); 855 ierr = PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);CHKERRQ(ierr); 856 } 857 ierr = PetscMalloc(n*sizeof(PetscInt),&ilink->fields);CHKERRQ(ierr); 858 ierr = PetscMemcpy(ilink->fields,fields,n*sizeof(PetscInt));CHKERRQ(ierr); 859 ilink->nfields = n; 860 ilink->next = PETSC_NULL; 861 ierr = KSPCreate(((PetscObject)pc)->comm,&ilink->ksp);CHKERRQ(ierr); 862 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 863 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 864 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 865 866 ierr = PetscSNPrintf(prefix,sizeof prefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",ilink->splitname);CHKERRQ(ierr); 867 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 868 869 if (!next) { 870 jac->head = ilink; 871 ilink->previous = PETSC_NULL; 872 } else { 873 while (next->next) { 874 next = next->next; 875 } 876 next->next = ilink; 877 ilink->previous = next; 878 } 879 jac->nsplits++; 880 PetscFunctionReturn(0); 881 } 882 EXTERN_C_END 883 884 EXTERN_C_BEGIN 885 #undef __FUNCT__ 886 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit_Schur" 887 PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp) 888 { 889 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 890 PetscErrorCode ierr; 891 892 PetscFunctionBegin; 893 ierr = PetscMalloc(jac->nsplits*sizeof(KSP),subksp);CHKERRQ(ierr); 894 ierr = MatSchurComplementGetKSP(jac->schur,*subksp);CHKERRQ(ierr); 895 (*subksp)[1] = jac->kspschur; 896 if (n) *n = jac->nsplits; 897 PetscFunctionReturn(0); 898 } 899 EXTERN_C_END 900 901 EXTERN_C_BEGIN 902 #undef __FUNCT__ 903 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit" 904 PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp) 905 { 906 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 907 PetscErrorCode ierr; 908 PetscInt cnt = 0; 909 PC_FieldSplitLink ilink = jac->head; 910 911 PetscFunctionBegin; 912 ierr = PetscMalloc(jac->nsplits*sizeof(KSP*),subksp);CHKERRQ(ierr); 913 while (ilink) { 914 (*subksp)[cnt++] = ilink->ksp; 915 ilink = ilink->next; 916 } 917 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); 918 if (n) *n = jac->nsplits; 919 PetscFunctionReturn(0); 920 } 921 EXTERN_C_END 922 923 EXTERN_C_BEGIN 924 #undef __FUNCT__ 925 #define __FUNCT__ "PCFieldSplitSetIS_FieldSplit" 926 PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is) 927 { 928 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 929 PetscErrorCode ierr; 930 PC_FieldSplitLink ilink, next = jac->head; 931 char prefix[128]; 932 933 PetscFunctionBegin; 934 if (jac->splitdefined) { 935 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 936 PetscFunctionReturn(0); 937 } 938 ierr = PetscNew(struct _PC_FieldSplitLink,&ilink);CHKERRQ(ierr); 939 if (splitname) { 940 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 941 } else { 942 ierr = PetscMalloc(3*sizeof(char),&ilink->splitname);CHKERRQ(ierr); 943 ierr = PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);CHKERRQ(ierr); 944 } 945 ilink->is = is; 946 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 947 ilink->next = PETSC_NULL; 948 ierr = KSPCreate(((PetscObject)pc)->comm,&ilink->ksp);CHKERRQ(ierr); 949 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 950 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 951 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 952 953 ierr = PetscSNPrintf(prefix,sizeof prefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",ilink->splitname);CHKERRQ(ierr); 954 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 955 956 if (!next) { 957 jac->head = ilink; 958 ilink->previous = PETSC_NULL; 959 } else { 960 while (next->next) { 961 next = next->next; 962 } 963 next->next = ilink; 964 ilink->previous = next; 965 } 966 jac->nsplits++; 967 968 PetscFunctionReturn(0); 969 } 970 EXTERN_C_END 971 972 #undef __FUNCT__ 973 #define __FUNCT__ "PCFieldSplitSetFields" 974 /*@ 975 PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner 976 977 Logically Collective on PC 978 979 Input Parameters: 980 + pc - the preconditioner context 981 . splitname - name of this split, if PETSC_NULL the number of the split is used 982 . n - the number of fields in this split 983 - fields - the fields in this split 984 985 Level: intermediate 986 987 Notes: Use PCFieldSplitSetIS() to set a completely general set of indices as a field. 988 989 The PCFieldSplitSetFields() is for defining fields as a strided blocks. For example, if the block 990 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 991 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x.... 992 where the numbered entries indicate what is in the field. 993 994 This function is called once per split (it creates a new split each time). Solve options 995 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 996 997 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize(), PCFieldSplitSetIS() 998 999 @*/ 1000 PetscErrorCode PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields) 1001 { 1002 PetscErrorCode ierr; 1003 1004 PetscFunctionBegin; 1005 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1006 PetscValidCharPointer(splitname,2); 1007 if (n < 1) SETERRQ2(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Provided number of fields %D in split \"%s\" not positive",n,splitname); 1008 PetscValidIntPointer(fields,3); 1009 ierr = PetscTryMethod(pc,"PCFieldSplitSetFields_C",(PC,const char[],PetscInt,const PetscInt *),(pc,splitname,n,fields));CHKERRQ(ierr); 1010 PetscFunctionReturn(0); 1011 } 1012 1013 #undef __FUNCT__ 1014 #define __FUNCT__ "PCFieldSplitSetIS" 1015 /*@ 1016 PCFieldSplitSetIS - Sets the exact elements for field 1017 1018 Logically Collective on PC 1019 1020 Input Parameters: 1021 + pc - the preconditioner context 1022 . splitname - name of this split, if PETSC_NULL the number of the split is used 1023 - is - the index set that defines the vector elements in this field 1024 1025 1026 Notes: 1027 Use PCFieldSplitSetFields(), for fields defined by strided types. 1028 1029 This function is called once per split (it creates a new split each time). Solve options 1030 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1031 1032 Level: intermediate 1033 1034 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize() 1035 1036 @*/ 1037 PetscErrorCode PCFieldSplitSetIS(PC pc,const char splitname[],IS is) 1038 { 1039 PetscErrorCode ierr; 1040 1041 PetscFunctionBegin; 1042 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1043 PetscValidCharPointer(splitname,2); 1044 PetscValidHeaderSpecific(is,IS_CLASSID,3); 1045 ierr = PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));CHKERRQ(ierr); 1046 PetscFunctionReturn(0); 1047 } 1048 1049 #undef __FUNCT__ 1050 #define __FUNCT__ "PCFieldSplitSetBlockSize" 1051 /*@ 1052 PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the 1053 fieldsplit preconditioner. If not set the matrix block size is used. 1054 1055 Logically Collective on PC 1056 1057 Input Parameters: 1058 + pc - the preconditioner context 1059 - bs - the block size 1060 1061 Level: intermediate 1062 1063 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields() 1064 1065 @*/ 1066 PetscErrorCode PCFieldSplitSetBlockSize(PC pc,PetscInt bs) 1067 { 1068 PetscErrorCode ierr; 1069 1070 PetscFunctionBegin; 1071 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1072 PetscValidLogicalCollectiveInt(pc,bs,2); 1073 ierr = PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));CHKERRQ(ierr); 1074 PetscFunctionReturn(0); 1075 } 1076 1077 #undef __FUNCT__ 1078 #define __FUNCT__ "PCFieldSplitGetSubKSP" 1079 /*@C 1080 PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits 1081 1082 Collective on KSP 1083 1084 Input Parameter: 1085 . pc - the preconditioner context 1086 1087 Output Parameters: 1088 + n - the number of splits 1089 - pc - the array of KSP contexts 1090 1091 Note: 1092 After PCFieldSplitGetSubKSP() the array of KSPs IS to be freed by the user 1093 (not the KSP just the array that contains them). 1094 1095 You must call KSPSetUp() before calling PCFieldSplitGetSubKSP(). 1096 1097 Level: advanced 1098 1099 .seealso: PCFIELDSPLIT 1100 @*/ 1101 PetscErrorCode PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[]) 1102 { 1103 PetscErrorCode ierr; 1104 1105 PetscFunctionBegin; 1106 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1107 if (n) PetscValidIntPointer(n,2); 1108 ierr = PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));CHKERRQ(ierr); 1109 PetscFunctionReturn(0); 1110 } 1111 1112 #undef __FUNCT__ 1113 #define __FUNCT__ "PCFieldSplitSchurPrecondition" 1114 /*@ 1115 PCFieldSplitSchurPrecondition - Indicates if the Schur complement is preconditioned by a preconditioner constructed by the 1116 A11 matrix. Otherwise no preconditioner is used. 1117 1118 Collective on PC 1119 1120 Input Parameters: 1121 + pc - the preconditioner context 1122 . ptype - which matrix to use for preconditioning the Schur complement 1123 - userpre - matrix to use for preconditioning, or PETSC_NULL 1124 1125 Notes: 1126 The default is to use the block on the diagonal of the preconditioning matrix. This is A11, in the (1,1) position. 1127 There are currently no preconditioners that work directly with the Schur complement so setting 1128 PC_FIELDSPLIT_SCHUR_PRE_SELF is observationally equivalent to -fieldsplit_1_pc_type none. 1129 1130 Options Database: 1131 . -pc_fieldsplit_schur_precondition <self,user,diag> default is diag 1132 1133 Level: intermediate 1134 1135 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType 1136 1137 @*/ 1138 PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1139 { 1140 PetscErrorCode ierr; 1141 1142 PetscFunctionBegin; 1143 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1144 ierr = PetscTryMethod(pc,"PCFieldSplitSchurPrecondition_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));CHKERRQ(ierr); 1145 PetscFunctionReturn(0); 1146 } 1147 1148 EXTERN_C_BEGIN 1149 #undef __FUNCT__ 1150 #define __FUNCT__ "PCFieldSplitSchurPrecondition_FieldSplit" 1151 PetscErrorCode PCFieldSplitSchurPrecondition_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1152 { 1153 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1154 PetscErrorCode ierr; 1155 1156 PetscFunctionBegin; 1157 jac->schurpre = ptype; 1158 if (pre) { 1159 if (jac->schur_user) {ierr = MatDestroy(jac->schur_user);CHKERRQ(ierr);} 1160 jac->schur_user = pre; 1161 ierr = PetscObjectReference((PetscObject)jac->schur_user);CHKERRQ(ierr); 1162 } 1163 PetscFunctionReturn(0); 1164 } 1165 EXTERN_C_END 1166 1167 #undef __FUNCT__ 1168 #define __FUNCT__ "PCFieldSplitGetSchurBlocks" 1169 /*@C 1170 PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement 1171 1172 Collective on KSP 1173 1174 Input Parameter: 1175 . pc - the preconditioner context 1176 1177 Output Parameters: 1178 + A00 - the (0,0) block 1179 . A01 - the (0,1) block 1180 . A10 - the (1,0) block 1181 - A11 - the (1,1) block 1182 1183 Level: advanced 1184 1185 .seealso: PCFIELDSPLIT 1186 @*/ 1187 PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11) 1188 { 1189 PC_FieldSplit *jac = (PC_FieldSplit *) pc->data; 1190 1191 PetscFunctionBegin; 1192 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1193 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach."); 1194 if (A00) *A00 = jac->pmat[0]; 1195 if (A01) *A01 = jac->B; 1196 if (A10) *A10 = jac->C; 1197 if (A11) *A11 = jac->pmat[1]; 1198 PetscFunctionReturn(0); 1199 } 1200 1201 EXTERN_C_BEGIN 1202 #undef __FUNCT__ 1203 #define __FUNCT__ "PCFieldSplitSetType_FieldSplit" 1204 PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type) 1205 { 1206 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1207 PetscErrorCode ierr; 1208 1209 PetscFunctionBegin; 1210 jac->type = type; 1211 if (type == PC_COMPOSITE_SCHUR) { 1212 pc->ops->apply = PCApply_FieldSplit_Schur; 1213 pc->ops->view = PCView_FieldSplit_Schur; 1214 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit_Schur",PCFieldSplitGetSubKSP_FieldSplit_Schur);CHKERRQ(ierr); 1215 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSchurPrecondition_C","PCFieldSplitSchurPrecondition_FieldSplit",PCFieldSplitSchurPrecondition_FieldSplit);CHKERRQ(ierr); 1216 1217 } else { 1218 pc->ops->apply = PCApply_FieldSplit; 1219 pc->ops->view = PCView_FieldSplit; 1220 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit",PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 1221 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSchurPrecondition_C","",0);CHKERRQ(ierr); 1222 } 1223 PetscFunctionReturn(0); 1224 } 1225 EXTERN_C_END 1226 1227 EXTERN_C_BEGIN 1228 #undef __FUNCT__ 1229 #define __FUNCT__ "PCFieldSplitSetBlockSize_FieldSplit" 1230 PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs) 1231 { 1232 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1233 1234 PetscFunctionBegin; 1235 if (bs < 1) SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs); 1236 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); 1237 jac->bs = bs; 1238 PetscFunctionReturn(0); 1239 } 1240 EXTERN_C_END 1241 1242 #undef __FUNCT__ 1243 #define __FUNCT__ "PCFieldSplitSetType" 1244 /*@ 1245 PCFieldSplitSetType - Sets the type of fieldsplit preconditioner. 1246 1247 Collective on PC 1248 1249 Input Parameter: 1250 . pc - the preconditioner context 1251 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 1252 1253 Options Database Key: 1254 . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type 1255 1256 Level: Developer 1257 1258 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 1259 1260 .seealso: PCCompositeSetType() 1261 1262 @*/ 1263 PetscErrorCode PCFieldSplitSetType(PC pc,PCCompositeType type) 1264 { 1265 PetscErrorCode ierr; 1266 1267 PetscFunctionBegin; 1268 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1269 ierr = PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));CHKERRQ(ierr); 1270 PetscFunctionReturn(0); 1271 } 1272 1273 /* -------------------------------------------------------------------------------------*/ 1274 /*MC 1275 PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual 1276 fields or groups of fields. See the users manual section "Solving Block Matrices" for more details. 1277 1278 To set options on the solvers for each block append -fieldsplit_ to all the PC 1279 options database keys. For example, -fieldsplit_pc_type ilu -fieldsplit_pc_factor_levels 1 1280 1281 To set the options on the solvers separate for each block call PCFieldSplitGetSubKSP() 1282 and set the options directly on the resulting KSP object 1283 1284 Level: intermediate 1285 1286 Options Database Keys: 1287 + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split 1288 . -pc_fieldsplit_default - automatically add any fields to additional splits that have not 1289 been supplied explicitly by -pc_fieldsplit_%d_fields 1290 . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields) 1291 . -pc_fieldsplit_type <additive,multiplicative,schur,symmetric_multiplicative> 1292 . -pc_fieldsplit_schur_precondition <true,false> default is true 1293 . -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero or negative diagonal and uses Schur complement with no preconditioner as the solver 1294 1295 - Options prefix for inner solvers when using Schur complement preconditioner are -fieldsplit_0_ and -fieldsplit_1_ 1296 for all other solvers they are -fieldsplit_%d_ for the dth field, use -fieldsplit_ for all fields 1297 1298 Notes: use PCFieldSplitSetFields() to set fields defined by "strided" entries and PCFieldSplitSetIS() 1299 to define a field by an arbitrary collection of entries. 1300 1301 If no fields are set the default is used. The fields are defined by entries strided by bs, 1302 beginning at 0 then 1, etc to bs-1. The block size can be set with PCFieldSplitSetBlockSize(), 1303 if this is not called the block size defaults to the blocksize of the second matrix passed 1304 to KSPSetOperators()/PCSetOperators(). 1305 1306 For the Schur complement preconditioner if J = ( A00 A01 ) 1307 ( A10 A11 ) 1308 the preconditioner is 1309 (I -B ksp(A00)) ( inv(A00) 0 (I 0 ) 1310 (0 I ) ( 0 ksp(S) ) (-A10 ksp(A00) I ) 1311 where the action of inv(A00) is applied using the KSP solver with prefix -fieldsplit_0_. The action of 1312 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). 1313 For PCFieldSplitGetKSP() when field number is 1314 0 it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default 1315 A11 is used to construct a preconditioner for S, use PCFieldSplitSchurPrecondition() to turn on or off this 1316 option. You can use the preconditioner PCLSC to precondition the Schur complement with -fieldsplit_1_pc_type lsc 1317 1318 If only one set of indices (one IS) is provided with PCFieldSplitSetIS() then the complement of that IS 1319 is used automatically for a second block. 1320 1321 Concepts: physics based preconditioners, block preconditioners 1322 1323 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners, PCLSC, 1324 PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSchurPrecondition() 1325 M*/ 1326 1327 EXTERN_C_BEGIN 1328 #undef __FUNCT__ 1329 #define __FUNCT__ "PCCreate_FieldSplit" 1330 PetscErrorCode PCCreate_FieldSplit(PC pc) 1331 { 1332 PetscErrorCode ierr; 1333 PC_FieldSplit *jac; 1334 1335 PetscFunctionBegin; 1336 ierr = PetscNewLog(pc,PC_FieldSplit,&jac);CHKERRQ(ierr); 1337 jac->bs = -1; 1338 jac->nsplits = 0; 1339 jac->type = PC_COMPOSITE_MULTIPLICATIVE; 1340 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */ 1341 jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACTORIZATION_FULL; 1342 1343 pc->data = (void*)jac; 1344 1345 pc->ops->apply = PCApply_FieldSplit; 1346 pc->ops->applytranspose = PCApplyTranspose_FieldSplit; 1347 pc->ops->setup = PCSetUp_FieldSplit; 1348 pc->ops->reset = PCReset_FieldSplit; 1349 pc->ops->destroy = PCDestroy_FieldSplit; 1350 pc->ops->setfromoptions = PCSetFromOptions_FieldSplit; 1351 pc->ops->view = PCView_FieldSplit; 1352 pc->ops->applyrichardson = 0; 1353 1354 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit", 1355 PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 1356 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetFields_C","PCFieldSplitSetFields_FieldSplit", 1357 PCFieldSplitSetFields_FieldSplit);CHKERRQ(ierr); 1358 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetIS_C","PCFieldSplitSetIS_FieldSplit", 1359 PCFieldSplitSetIS_FieldSplit);CHKERRQ(ierr); 1360 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetType_C","PCFieldSplitSetType_FieldSplit", 1361 PCFieldSplitSetType_FieldSplit);CHKERRQ(ierr); 1362 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetBlockSize_C","PCFieldSplitSetBlockSize_FieldSplit", 1363 PCFieldSplitSetBlockSize_FieldSplit);CHKERRQ(ierr); 1364 PetscFunctionReturn(0); 1365 } 1366 EXTERN_C_END 1367 1368 1369 1370