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