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