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__ "PCReset_FieldSplit" 726 static PetscErrorCode PCReset_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 = KSPReset(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 ilink = next; 741 } 742 ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr); 743 if (jac->mat && jac->mat != jac->pmat) { 744 ierr = MatDestroyMatrices(jac->nsplits,&jac->mat);CHKERRQ(ierr); 745 } else if (jac->mat) { 746 jac->mat = PETSC_NULL; 747 } 748 if (jac->pmat) {ierr = MatDestroyMatrices(jac->nsplits,&jac->pmat);CHKERRQ(ierr);} 749 if (jac->Afield) {ierr = MatDestroyMatrices(jac->nsplits,&jac->Afield);CHKERRQ(ierr);} 750 if (jac->w1) {ierr = VecDestroy(jac->w1);CHKERRQ(ierr);} 751 if (jac->w2) {ierr = VecDestroy(jac->w2);CHKERRQ(ierr);} 752 if (jac->schur) {ierr = MatDestroy(jac->schur);CHKERRQ(ierr);} 753 if (jac->schur_user) {ierr = MatDestroy(jac->schur_user);CHKERRQ(ierr);} 754 if (jac->kspschur) {ierr = KSPDestroy(jac->kspschur);CHKERRQ(ierr);} 755 if (jac->B) {ierr = MatDestroy(jac->B);CHKERRQ(ierr);} 756 if (jac->C) {ierr = MatDestroy(jac->C);CHKERRQ(ierr);} 757 PetscFunctionReturn(0); 758 } 759 760 #undef __FUNCT__ 761 #define __FUNCT__ "PCDestroy_FieldSplit" 762 static PetscErrorCode PCDestroy_FieldSplit(PC pc) 763 { 764 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 765 PetscErrorCode ierr; 766 PC_FieldSplitLink ilink = jac->head,next; 767 768 PetscFunctionBegin; 769 ierr = PCReset_FieldSplit(pc);CHKERRQ(ierr); 770 while (ilink) { 771 ierr = KSPDestroy(ilink->ksp);CHKERRQ(ierr); 772 next = ilink->next; 773 ierr = PetscFree(ilink->splitname);CHKERRQ(ierr); 774 ierr = PetscFree(ilink->fields);CHKERRQ(ierr); 775 ierr = PetscFree(ilink);CHKERRQ(ierr); 776 ilink = next; 777 } 778 ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr); 779 ierr = PetscFree(pc->data);CHKERRQ(ierr); 780 PetscFunctionReturn(0); 781 } 782 783 #undef __FUNCT__ 784 #define __FUNCT__ "PCSetFromOptions_FieldSplit" 785 static PetscErrorCode PCSetFromOptions_FieldSplit(PC pc) 786 { 787 PetscErrorCode ierr; 788 PetscInt bs; 789 PetscBool flg,stokes = PETSC_FALSE; 790 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 791 PCCompositeType ctype; 792 793 PetscFunctionBegin; 794 ierr = PetscOptionsHead("FieldSplit options");CHKERRQ(ierr); 795 ierr = PetscOptionsBool("-pc_fieldsplit_real_diagonal","Use diagonal blocks of the operator","PCFieldSplitSetRealDiagonal",jac->realdiagonal,&jac->realdiagonal,PETSC_NULL);CHKERRQ(ierr); 796 ierr = PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);CHKERRQ(ierr); 797 if (flg) { 798 ierr = PCFieldSplitSetBlockSize(pc,bs);CHKERRQ(ierr); 799 } 800 801 ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,PETSC_NULL);CHKERRQ(ierr); 802 if (stokes) { 803 ierr = PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);CHKERRQ(ierr); 804 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_SELF; 805 } 806 807 ierr = PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg);CHKERRQ(ierr); 808 if (flg) { 809 ierr = PCFieldSplitSetType(pc,ctype);CHKERRQ(ierr); 810 } 811 812 /* Only setup fields once */ 813 if ((jac->bs > 0) && (jac->nsplits == 0)) { 814 /* only allow user to set fields from command line if bs is already known. 815 otherwise user can set them in PCFieldSplitSetDefaults() */ 816 ierr = PCFieldSplitSetRuntimeSplits_Private(pc);CHKERRQ(ierr); 817 if (jac->splitdefined) {ierr = PetscInfo(pc,"Splits defined using the options database\n");CHKERRQ(ierr);} 818 } 819 if (jac->type == PC_COMPOSITE_SCHUR) { 820 ierr = PetscOptionsEnum("-pc_fieldsplit_schur_factorization_type","Factorization to use","None",PCFieldSplitSchurFactorizationTypes,(PetscEnum)jac->schurfactorization,(PetscEnum*)&jac->schurfactorization,PETSC_NULL);CHKERRQ(ierr); 821 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); 822 } 823 ierr = PetscOptionsTail();CHKERRQ(ierr); 824 PetscFunctionReturn(0); 825 } 826 827 /*------------------------------------------------------------------------------------*/ 828 829 EXTERN_C_BEGIN 830 #undef __FUNCT__ 831 #define __FUNCT__ "PCFieldSplitSetFields_FieldSplit" 832 PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields) 833 { 834 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 835 PetscErrorCode ierr; 836 PC_FieldSplitLink ilink,next = jac->head; 837 char prefix[128]; 838 PetscInt i; 839 840 PetscFunctionBegin; 841 if (jac->splitdefined) { 842 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 843 PetscFunctionReturn(0); 844 } 845 for (i=0; i<n; i++) { 846 if (fields[i] >= jac->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs); 847 if (fields[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]); 848 } 849 ierr = PetscNew(struct _PC_FieldSplitLink,&ilink);CHKERRQ(ierr); 850 if (splitname) { 851 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 852 } else { 853 ierr = PetscMalloc(3*sizeof(char),&ilink->splitname);CHKERRQ(ierr); 854 ierr = PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);CHKERRQ(ierr); 855 } 856 ierr = PetscMalloc(n*sizeof(PetscInt),&ilink->fields);CHKERRQ(ierr); 857 ierr = PetscMemcpy(ilink->fields,fields,n*sizeof(PetscInt));CHKERRQ(ierr); 858 ilink->nfields = n; 859 ilink->next = PETSC_NULL; 860 ierr = KSPCreate(((PetscObject)pc)->comm,&ilink->ksp);CHKERRQ(ierr); 861 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 862 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 863 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 864 865 ierr = PetscSNPrintf(prefix,sizeof prefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",ilink->splitname);CHKERRQ(ierr); 866 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 867 868 if (!next) { 869 jac->head = ilink; 870 ilink->previous = PETSC_NULL; 871 } else { 872 while (next->next) { 873 next = next->next; 874 } 875 next->next = ilink; 876 ilink->previous = next; 877 } 878 jac->nsplits++; 879 PetscFunctionReturn(0); 880 } 881 EXTERN_C_END 882 883 EXTERN_C_BEGIN 884 #undef __FUNCT__ 885 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit_Schur" 886 PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp) 887 { 888 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 889 PetscErrorCode ierr; 890 891 PetscFunctionBegin; 892 ierr = PetscMalloc(jac->nsplits*sizeof(KSP),subksp);CHKERRQ(ierr); 893 ierr = MatSchurComplementGetKSP(jac->schur,*subksp);CHKERRQ(ierr); 894 (*subksp)[1] = jac->kspschur; 895 if (n) *n = jac->nsplits; 896 PetscFunctionReturn(0); 897 } 898 EXTERN_C_END 899 900 EXTERN_C_BEGIN 901 #undef __FUNCT__ 902 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit" 903 PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp) 904 { 905 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 906 PetscErrorCode ierr; 907 PetscInt cnt = 0; 908 PC_FieldSplitLink ilink = jac->head; 909 910 PetscFunctionBegin; 911 ierr = PetscMalloc(jac->nsplits*sizeof(KSP*),subksp);CHKERRQ(ierr); 912 while (ilink) { 913 (*subksp)[cnt++] = ilink->ksp; 914 ilink = ilink->next; 915 } 916 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); 917 if (n) *n = jac->nsplits; 918 PetscFunctionReturn(0); 919 } 920 EXTERN_C_END 921 922 EXTERN_C_BEGIN 923 #undef __FUNCT__ 924 #define __FUNCT__ "PCFieldSplitSetIS_FieldSplit" 925 PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is) 926 { 927 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 928 PetscErrorCode ierr; 929 PC_FieldSplitLink ilink, next = jac->head; 930 char prefix[128]; 931 932 PetscFunctionBegin; 933 if (jac->splitdefined) { 934 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 935 PetscFunctionReturn(0); 936 } 937 ierr = PetscNew(struct _PC_FieldSplitLink,&ilink);CHKERRQ(ierr); 938 if (splitname) { 939 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 940 } else { 941 ierr = PetscMalloc(3*sizeof(char),&ilink->splitname);CHKERRQ(ierr); 942 ierr = PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);CHKERRQ(ierr); 943 } 944 ilink->is = is; 945 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 946 ilink->next = PETSC_NULL; 947 ierr = KSPCreate(((PetscObject)pc)->comm,&ilink->ksp);CHKERRQ(ierr); 948 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 949 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 950 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 951 952 ierr = PetscSNPrintf(prefix,sizeof prefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",ilink->splitname);CHKERRQ(ierr); 953 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 954 955 if (!next) { 956 jac->head = ilink; 957 ilink->previous = PETSC_NULL; 958 } else { 959 while (next->next) { 960 next = next->next; 961 } 962 next->next = ilink; 963 ilink->previous = next; 964 } 965 jac->nsplits++; 966 967 PetscFunctionReturn(0); 968 } 969 EXTERN_C_END 970 971 #undef __FUNCT__ 972 #define __FUNCT__ "PCFieldSplitSetFields" 973 /*@ 974 PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner 975 976 Logically Collective on PC 977 978 Input Parameters: 979 + pc - the preconditioner context 980 . splitname - name of this split, if PETSC_NULL the number of the split is used 981 . n - the number of fields in this split 982 - fields - the fields in this split 983 984 Level: intermediate 985 986 Notes: Use PCFieldSplitSetIS() to set a completely general set of indices as a field. 987 988 The PCFieldSplitSetFields() is for defining fields as a strided blocks. For example, if the block 989 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 990 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x.... 991 where the numbered entries indicate what is in the field. 992 993 This function is called once per split (it creates a new split each time). Solve options 994 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 995 996 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize(), PCFieldSplitSetIS() 997 998 @*/ 999 PetscErrorCode PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields) 1000 { 1001 PetscErrorCode ierr; 1002 1003 PetscFunctionBegin; 1004 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1005 PetscValidCharPointer(splitname,2); 1006 if (n < 1) SETERRQ2(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Provided number of fields %D in split \"%s\" not positive",n,splitname); 1007 PetscValidIntPointer(fields,3); 1008 ierr = PetscTryMethod(pc,"PCFieldSplitSetFields_C",(PC,const char[],PetscInt,const PetscInt *),(pc,splitname,n,fields));CHKERRQ(ierr); 1009 PetscFunctionReturn(0); 1010 } 1011 1012 #undef __FUNCT__ 1013 #define __FUNCT__ "PCFieldSplitSetIS" 1014 /*@ 1015 PCFieldSplitSetIS - Sets the exact elements for field 1016 1017 Logically Collective on PC 1018 1019 Input Parameters: 1020 + pc - the preconditioner context 1021 . splitname - name of this split, if PETSC_NULL the number of the split is used 1022 - is - the index set that defines the vector elements in this field 1023 1024 1025 Notes: 1026 Use PCFieldSplitSetFields(), for fields defined by strided types. 1027 1028 This function is called once per split (it creates a new split each time). Solve options 1029 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1030 1031 Level: intermediate 1032 1033 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize() 1034 1035 @*/ 1036 PetscErrorCode PCFieldSplitSetIS(PC pc,const char splitname[],IS is) 1037 { 1038 PetscErrorCode ierr; 1039 1040 PetscFunctionBegin; 1041 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1042 PetscValidCharPointer(splitname,2); 1043 PetscValidHeaderSpecific(is,IS_CLASSID,3); 1044 ierr = PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));CHKERRQ(ierr); 1045 PetscFunctionReturn(0); 1046 } 1047 1048 #undef __FUNCT__ 1049 #define __FUNCT__ "PCFieldSplitSetBlockSize" 1050 /*@ 1051 PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the 1052 fieldsplit preconditioner. If not set the matrix block size is used. 1053 1054 Logically Collective on PC 1055 1056 Input Parameters: 1057 + pc - the preconditioner context 1058 - bs - the block size 1059 1060 Level: intermediate 1061 1062 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields() 1063 1064 @*/ 1065 PetscErrorCode PCFieldSplitSetBlockSize(PC pc,PetscInt bs) 1066 { 1067 PetscErrorCode ierr; 1068 1069 PetscFunctionBegin; 1070 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1071 PetscValidLogicalCollectiveInt(pc,bs,2); 1072 ierr = PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));CHKERRQ(ierr); 1073 PetscFunctionReturn(0); 1074 } 1075 1076 #undef __FUNCT__ 1077 #define __FUNCT__ "PCFieldSplitGetSubKSP" 1078 /*@C 1079 PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits 1080 1081 Collective on KSP 1082 1083 Input Parameter: 1084 . pc - the preconditioner context 1085 1086 Output Parameters: 1087 + n - the number of splits 1088 - pc - the array of KSP contexts 1089 1090 Note: 1091 After PCFieldSplitGetSubKSP() the array of KSPs IS to be freed by the user 1092 (not the KSP just the array that contains them). 1093 1094 You must call KSPSetUp() before calling PCFieldSplitGetSubKSP(). 1095 1096 Level: advanced 1097 1098 .seealso: PCFIELDSPLIT 1099 @*/ 1100 PetscErrorCode PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[]) 1101 { 1102 PetscErrorCode ierr; 1103 1104 PetscFunctionBegin; 1105 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1106 if (n) PetscValidIntPointer(n,2); 1107 ierr = PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));CHKERRQ(ierr); 1108 PetscFunctionReturn(0); 1109 } 1110 1111 #undef __FUNCT__ 1112 #define __FUNCT__ "PCFieldSplitSchurPrecondition" 1113 /*@ 1114 PCFieldSplitSchurPrecondition - Indicates if the Schur complement is preconditioned by a preconditioner constructed by the 1115 A11 matrix. Otherwise no preconditioner is used. 1116 1117 Collective on PC 1118 1119 Input Parameters: 1120 + pc - the preconditioner context 1121 . ptype - which matrix to use for preconditioning the Schur complement 1122 - userpre - matrix to use for preconditioning, or PETSC_NULL 1123 1124 Notes: 1125 The default is to use the block on the diagonal of the preconditioning matrix. This is A11, in the (1,1) position. 1126 There are currently no preconditioners that work directly with the Schur complement so setting 1127 PC_FIELDSPLIT_SCHUR_PRE_SELF is observationally equivalent to -fieldsplit_1_pc_type none. 1128 1129 Options Database: 1130 . -pc_fieldsplit_schur_precondition <self,user,diag> default is diag 1131 1132 Level: intermediate 1133 1134 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType 1135 1136 @*/ 1137 PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1138 { 1139 PetscErrorCode ierr; 1140 1141 PetscFunctionBegin; 1142 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1143 ierr = PetscTryMethod(pc,"PCFieldSplitSchurPrecondition_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));CHKERRQ(ierr); 1144 PetscFunctionReturn(0); 1145 } 1146 1147 EXTERN_C_BEGIN 1148 #undef __FUNCT__ 1149 #define __FUNCT__ "PCFieldSplitSchurPrecondition_FieldSplit" 1150 PetscErrorCode PCFieldSplitSchurPrecondition_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1151 { 1152 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1153 PetscErrorCode ierr; 1154 1155 PetscFunctionBegin; 1156 jac->schurpre = ptype; 1157 if (pre) { 1158 if (jac->schur_user) {ierr = MatDestroy(jac->schur_user);CHKERRQ(ierr);} 1159 jac->schur_user = pre; 1160 ierr = PetscObjectReference((PetscObject)jac->schur_user);CHKERRQ(ierr); 1161 } 1162 PetscFunctionReturn(0); 1163 } 1164 EXTERN_C_END 1165 1166 #undef __FUNCT__ 1167 #define __FUNCT__ "PCFieldSplitGetSchurBlocks" 1168 /*@C 1169 PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement 1170 1171 Collective on KSP 1172 1173 Input Parameter: 1174 . pc - the preconditioner context 1175 1176 Output Parameters: 1177 + A00 - the (0,0) block 1178 . A01 - the (0,1) block 1179 . A10 - the (1,0) block 1180 - A11 - the (1,1) block 1181 1182 Level: advanced 1183 1184 .seealso: PCFIELDSPLIT 1185 @*/ 1186 PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11) 1187 { 1188 PC_FieldSplit *jac = (PC_FieldSplit *) pc->data; 1189 1190 PetscFunctionBegin; 1191 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1192 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach."); 1193 if (A00) *A00 = jac->pmat[0]; 1194 if (A01) *A01 = jac->B; 1195 if (A10) *A10 = jac->C; 1196 if (A11) *A11 = jac->pmat[1]; 1197 PetscFunctionReturn(0); 1198 } 1199 1200 EXTERN_C_BEGIN 1201 #undef __FUNCT__ 1202 #define __FUNCT__ "PCFieldSplitSetType_FieldSplit" 1203 PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type) 1204 { 1205 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1206 PetscErrorCode ierr; 1207 1208 PetscFunctionBegin; 1209 jac->type = type; 1210 if (type == PC_COMPOSITE_SCHUR) { 1211 pc->ops->apply = PCApply_FieldSplit_Schur; 1212 pc->ops->view = PCView_FieldSplit_Schur; 1213 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit_Schur",PCFieldSplitGetSubKSP_FieldSplit_Schur);CHKERRQ(ierr); 1214 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSchurPrecondition_C","PCFieldSplitSchurPrecondition_FieldSplit",PCFieldSplitSchurPrecondition_FieldSplit);CHKERRQ(ierr); 1215 1216 } else { 1217 pc->ops->apply = PCApply_FieldSplit; 1218 pc->ops->view = PCView_FieldSplit; 1219 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit",PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 1220 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSchurPrecondition_C","",0);CHKERRQ(ierr); 1221 } 1222 PetscFunctionReturn(0); 1223 } 1224 EXTERN_C_END 1225 1226 EXTERN_C_BEGIN 1227 #undef __FUNCT__ 1228 #define __FUNCT__ "PCFieldSplitSetBlockSize_FieldSplit" 1229 PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs) 1230 { 1231 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1232 1233 PetscFunctionBegin; 1234 if (bs < 1) SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs); 1235 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); 1236 jac->bs = bs; 1237 PetscFunctionReturn(0); 1238 } 1239 EXTERN_C_END 1240 1241 #undef __FUNCT__ 1242 #define __FUNCT__ "PCFieldSplitSetType" 1243 /*@ 1244 PCFieldSplitSetType - Sets the type of fieldsplit preconditioner. 1245 1246 Collective on PC 1247 1248 Input Parameter: 1249 . pc - the preconditioner context 1250 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 1251 1252 Options Database Key: 1253 . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type 1254 1255 Level: Developer 1256 1257 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 1258 1259 .seealso: PCCompositeSetType() 1260 1261 @*/ 1262 PetscErrorCode PCFieldSplitSetType(PC pc,PCCompositeType type) 1263 { 1264 PetscErrorCode ierr; 1265 1266 PetscFunctionBegin; 1267 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1268 ierr = PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));CHKERRQ(ierr); 1269 PetscFunctionReturn(0); 1270 } 1271 1272 /* -------------------------------------------------------------------------------------*/ 1273 /*MC 1274 PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual 1275 fields or groups of fields. See the users manual section "Solving Block Matrices" for more details. 1276 1277 To set options on the solvers for each block append -fieldsplit_ to all the PC 1278 options database keys. For example, -fieldsplit_pc_type ilu -fieldsplit_pc_factor_levels 1 1279 1280 To set the options on the solvers separate for each block call PCFieldSplitGetSubKSP() 1281 and set the options directly on the resulting KSP object 1282 1283 Level: intermediate 1284 1285 Options Database Keys: 1286 + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split 1287 . -pc_fieldsplit_default - automatically add any fields to additional splits that have not 1288 been supplied explicitly by -pc_fieldsplit_%d_fields 1289 . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields) 1290 . -pc_fieldsplit_type <additive,multiplicative,schur,symmetric_multiplicative> 1291 . -pc_fieldsplit_schur_precondition <true,false> default is true 1292 . -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero or negative diagonal and uses Schur complement with no preconditioner as the solver 1293 1294 - Options prefix for inner solvers when using Schur complement preconditioner are -fieldsplit_0_ and -fieldsplit_1_ 1295 for all other solvers they are -fieldsplit_%d_ for the dth field, use -fieldsplit_ for all fields 1296 1297 Notes: use PCFieldSplitSetFields() to set fields defined by "strided" entries and PCFieldSplitSetIS() 1298 to define a field by an arbitrary collection of entries. 1299 1300 If no fields are set the default is used. The fields are defined by entries strided by bs, 1301 beginning at 0 then 1, etc to bs-1. The block size can be set with PCFieldSplitSetBlockSize(), 1302 if this is not called the block size defaults to the blocksize of the second matrix passed 1303 to KSPSetOperators()/PCSetOperators(). 1304 1305 For the Schur complement preconditioner if J = ( A00 A01 ) 1306 ( A10 A11 ) 1307 the preconditioner is 1308 (I -B ksp(A00)) ( inv(A00) 0 (I 0 ) 1309 (0 I ) ( 0 ksp(S) ) (-A10 ksp(A00) I ) 1310 where the action of inv(A00) is applied using the KSP solver with prefix -fieldsplit_0_. The action of 1311 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). 1312 For PCFieldSplitGetKSP() when field number is 1313 0 it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default 1314 A11 is used to construct a preconditioner for S, use PCFieldSplitSchurPrecondition() to turn on or off this 1315 option. You can use the preconditioner PCLSC to precondition the Schur complement with -fieldsplit_1_pc_type lsc 1316 1317 If only one set of indices (one IS) is provided with PCFieldSplitSetIS() then the complement of that IS 1318 is used automatically for a second block. 1319 1320 Concepts: physics based preconditioners, block preconditioners 1321 1322 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners, PCLSC, 1323 PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSchurPrecondition() 1324 M*/ 1325 1326 EXTERN_C_BEGIN 1327 #undef __FUNCT__ 1328 #define __FUNCT__ "PCCreate_FieldSplit" 1329 PetscErrorCode PCCreate_FieldSplit(PC pc) 1330 { 1331 PetscErrorCode ierr; 1332 PC_FieldSplit *jac; 1333 1334 PetscFunctionBegin; 1335 ierr = PetscNewLog(pc,PC_FieldSplit,&jac);CHKERRQ(ierr); 1336 jac->bs = -1; 1337 jac->nsplits = 0; 1338 jac->type = PC_COMPOSITE_MULTIPLICATIVE; 1339 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */ 1340 jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACTORIZATION_FULL; 1341 1342 pc->data = (void*)jac; 1343 1344 pc->ops->apply = PCApply_FieldSplit; 1345 pc->ops->applytranspose = PCApplyTranspose_FieldSplit; 1346 pc->ops->setup = PCSetUp_FieldSplit; 1347 pc->ops->reset = PCReset_FieldSplit; 1348 pc->ops->destroy = PCDestroy_FieldSplit; 1349 pc->ops->setfromoptions = PCSetFromOptions_FieldSplit; 1350 pc->ops->view = PCView_FieldSplit; 1351 pc->ops->applyrichardson = 0; 1352 1353 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit", 1354 PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 1355 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetFields_C","PCFieldSplitSetFields_FieldSplit", 1356 PCFieldSplitSetFields_FieldSplit);CHKERRQ(ierr); 1357 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetIS_C","PCFieldSplitSetIS_FieldSplit", 1358 PCFieldSplitSetIS_FieldSplit);CHKERRQ(ierr); 1359 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetType_C","PCFieldSplitSetType_FieldSplit", 1360 PCFieldSplitSetType_FieldSplit);CHKERRQ(ierr); 1361 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetBlockSize_C","PCFieldSplitSetBlockSize_FieldSplit", 1362 PCFieldSplitSetBlockSize_FieldSplit);CHKERRQ(ierr); 1363 PetscFunctionReturn(0); 1364 } 1365 EXTERN_C_END 1366 1367 1368 1369