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 PetscFunctionReturn(0); 823 } 824 825 #undef __FUNCT__ 826 #define __FUNCT__ "PCSetFromOptions_FieldSplit" 827 static PetscErrorCode PCSetFromOptions_FieldSplit(PC pc) 828 { 829 PetscErrorCode ierr; 830 PetscInt bs; 831 PetscBool flg,stokes = PETSC_FALSE; 832 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 833 PCCompositeType ctype; 834 835 PetscFunctionBegin; 836 ierr = PetscOptionsHead("FieldSplit options");CHKERRQ(ierr); 837 ierr = PetscOptionsBool("-pc_fieldsplit_real_diagonal","Use diagonal blocks of the operator","PCFieldSplitSetRealDiagonal",jac->realdiagonal,&jac->realdiagonal,PETSC_NULL);CHKERRQ(ierr); 838 ierr = PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);CHKERRQ(ierr); 839 if (flg) { 840 ierr = PCFieldSplitSetBlockSize(pc,bs);CHKERRQ(ierr); 841 } 842 843 ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,PETSC_NULL);CHKERRQ(ierr); 844 if (stokes) { 845 ierr = PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);CHKERRQ(ierr); 846 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_SELF; 847 } 848 849 ierr = PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg);CHKERRQ(ierr); 850 if (flg) { 851 ierr = PCFieldSplitSetType(pc,ctype);CHKERRQ(ierr); 852 } 853 854 /* Only setup fields once */ 855 if ((jac->bs > 0) && (jac->nsplits == 0)) { 856 /* only allow user to set fields from command line if bs is already known. 857 otherwise user can set them in PCFieldSplitSetDefaults() */ 858 ierr = PCFieldSplitSetRuntimeSplits_Private(pc);CHKERRQ(ierr); 859 if (jac->splitdefined) {ierr = PetscInfo(pc,"Splits defined using the options database\n");CHKERRQ(ierr);} 860 } 861 if (jac->type == PC_COMPOSITE_SCHUR) { 862 ierr = PetscOptionsEnum("-pc_fieldsplit_schur_factorization_type","Factorization to use","None",PCFieldSplitSchurFactorizationTypes,(PetscEnum)jac->schurfactorization,(PetscEnum*)&jac->schurfactorization,PETSC_NULL);CHKERRQ(ierr); 863 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); 864 } 865 ierr = PetscOptionsTail();CHKERRQ(ierr); 866 PetscFunctionReturn(0); 867 } 868 869 /*------------------------------------------------------------------------------------*/ 870 871 EXTERN_C_BEGIN 872 #undef __FUNCT__ 873 #define __FUNCT__ "PCFieldSplitSetFields_FieldSplit" 874 PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields) 875 { 876 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 877 PetscErrorCode ierr; 878 PC_FieldSplitLink ilink,next = jac->head; 879 char prefix[128]; 880 PetscInt i; 881 882 PetscFunctionBegin; 883 if (jac->splitdefined) { 884 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 885 PetscFunctionReturn(0); 886 } 887 for (i=0; i<n; i++) { 888 if (fields[i] >= jac->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs); 889 if (fields[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]); 890 } 891 ierr = PetscNew(struct _PC_FieldSplitLink,&ilink);CHKERRQ(ierr); 892 if (splitname) { 893 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 894 } else { 895 ierr = PetscMalloc(3*sizeof(char),&ilink->splitname);CHKERRQ(ierr); 896 ierr = PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);CHKERRQ(ierr); 897 } 898 ierr = PetscMalloc(n*sizeof(PetscInt),&ilink->fields);CHKERRQ(ierr); 899 ierr = PetscMemcpy(ilink->fields,fields,n*sizeof(PetscInt));CHKERRQ(ierr); 900 ilink->nfields = n; 901 ilink->next = PETSC_NULL; 902 ierr = KSPCreate(((PetscObject)pc)->comm,&ilink->ksp);CHKERRQ(ierr); 903 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 904 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 905 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 906 907 ierr = PetscSNPrintf(prefix,sizeof prefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",ilink->splitname);CHKERRQ(ierr); 908 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 909 910 if (!next) { 911 jac->head = ilink; 912 ilink->previous = PETSC_NULL; 913 } else { 914 while (next->next) { 915 next = next->next; 916 } 917 next->next = ilink; 918 ilink->previous = next; 919 } 920 jac->nsplits++; 921 PetscFunctionReturn(0); 922 } 923 EXTERN_C_END 924 925 EXTERN_C_BEGIN 926 #undef __FUNCT__ 927 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit_Schur" 928 PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp) 929 { 930 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 931 PetscErrorCode ierr; 932 933 PetscFunctionBegin; 934 ierr = PetscMalloc(jac->nsplits*sizeof(KSP),subksp);CHKERRQ(ierr); 935 ierr = MatSchurComplementGetKSP(jac->schur,*subksp);CHKERRQ(ierr); 936 (*subksp)[1] = jac->kspschur; 937 if (n) *n = jac->nsplits; 938 PetscFunctionReturn(0); 939 } 940 EXTERN_C_END 941 942 EXTERN_C_BEGIN 943 #undef __FUNCT__ 944 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit" 945 PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp) 946 { 947 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 948 PetscErrorCode ierr; 949 PetscInt cnt = 0; 950 PC_FieldSplitLink ilink = jac->head; 951 952 PetscFunctionBegin; 953 ierr = PetscMalloc(jac->nsplits*sizeof(KSP),subksp);CHKERRQ(ierr); 954 while (ilink) { 955 (*subksp)[cnt++] = ilink->ksp; 956 ilink = ilink->next; 957 } 958 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); 959 if (n) *n = jac->nsplits; 960 PetscFunctionReturn(0); 961 } 962 EXTERN_C_END 963 964 EXTERN_C_BEGIN 965 #undef __FUNCT__ 966 #define __FUNCT__ "PCFieldSplitSetIS_FieldSplit" 967 PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is) 968 { 969 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 970 PetscErrorCode ierr; 971 PC_FieldSplitLink ilink, next = jac->head; 972 char prefix[128]; 973 974 PetscFunctionBegin; 975 if (jac->splitdefined) { 976 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 977 PetscFunctionReturn(0); 978 } 979 ierr = PetscNew(struct _PC_FieldSplitLink,&ilink);CHKERRQ(ierr); 980 if (splitname) { 981 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 982 } else { 983 ierr = PetscMalloc(3*sizeof(char),&ilink->splitname);CHKERRQ(ierr); 984 ierr = PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);CHKERRQ(ierr); 985 } 986 ilink->is = is; 987 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 988 ilink->next = PETSC_NULL; 989 ierr = KSPCreate(((PetscObject)pc)->comm,&ilink->ksp);CHKERRQ(ierr); 990 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 991 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 992 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 993 994 ierr = PetscSNPrintf(prefix,sizeof prefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",ilink->splitname);CHKERRQ(ierr); 995 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 996 997 if (!next) { 998 jac->head = ilink; 999 ilink->previous = PETSC_NULL; 1000 } else { 1001 while (next->next) { 1002 next = next->next; 1003 } 1004 next->next = ilink; 1005 ilink->previous = next; 1006 } 1007 jac->nsplits++; 1008 1009 PetscFunctionReturn(0); 1010 } 1011 EXTERN_C_END 1012 1013 #undef __FUNCT__ 1014 #define __FUNCT__ "PCFieldSplitSetFields" 1015 /*@ 1016 PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner 1017 1018 Logically Collective on PC 1019 1020 Input Parameters: 1021 + pc - the preconditioner context 1022 . splitname - name of this split, if PETSC_NULL the number of the split is used 1023 . n - the number of fields in this split 1024 - fields - the fields in this split 1025 1026 Level: intermediate 1027 1028 Notes: Use PCFieldSplitSetIS() to set a completely general set of indices as a field. 1029 1030 The PCFieldSplitSetFields() is for defining fields as a strided blocks. For example, if the block 1031 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 1032 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x.... 1033 where the numbered entries indicate what is in the field. 1034 1035 This function is called once per split (it creates a new split each time). Solve options 1036 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1037 1038 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize(), PCFieldSplitSetIS() 1039 1040 @*/ 1041 PetscErrorCode PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields) 1042 { 1043 PetscErrorCode ierr; 1044 1045 PetscFunctionBegin; 1046 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1047 PetscValidCharPointer(splitname,2); 1048 if (n < 1) SETERRQ2(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Provided number of fields %D in split \"%s\" not positive",n,splitname); 1049 PetscValidIntPointer(fields,3); 1050 ierr = PetscTryMethod(pc,"PCFieldSplitSetFields_C",(PC,const char[],PetscInt,const PetscInt *),(pc,splitname,n,fields));CHKERRQ(ierr); 1051 PetscFunctionReturn(0); 1052 } 1053 1054 #undef __FUNCT__ 1055 #define __FUNCT__ "PCFieldSplitSetIS" 1056 /*@ 1057 PCFieldSplitSetIS - Sets the exact elements for field 1058 1059 Logically Collective on PC 1060 1061 Input Parameters: 1062 + pc - the preconditioner context 1063 . splitname - name of this split, if PETSC_NULL the number of the split is used 1064 - is - the index set that defines the vector elements in this field 1065 1066 1067 Notes: 1068 Use PCFieldSplitSetFields(), for fields defined by strided types. 1069 1070 This function is called once per split (it creates a new split each time). Solve options 1071 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1072 1073 Level: intermediate 1074 1075 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize() 1076 1077 @*/ 1078 PetscErrorCode PCFieldSplitSetIS(PC pc,const char splitname[],IS is) 1079 { 1080 PetscErrorCode ierr; 1081 1082 PetscFunctionBegin; 1083 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1084 PetscValidCharPointer(splitname,2); 1085 PetscValidHeaderSpecific(is,IS_CLASSID,3); 1086 ierr = PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));CHKERRQ(ierr); 1087 PetscFunctionReturn(0); 1088 } 1089 1090 #undef __FUNCT__ 1091 #define __FUNCT__ "PCFieldSplitGetIS" 1092 /*@ 1093 PCFieldSplitGetIS - Retrieves the elements for a field as an IS 1094 1095 Logically Collective on PC 1096 1097 Input Parameters: 1098 + pc - the preconditioner context 1099 - splitname - name of this split 1100 1101 Output Parameter: 1102 - is - the index set that defines the vector elements in this field, or PETSC_NULL if the field is not found 1103 1104 Level: intermediate 1105 1106 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetIS() 1107 1108 @*/ 1109 PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is) 1110 { 1111 PetscErrorCode ierr; 1112 1113 PetscFunctionBegin; 1114 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1115 PetscValidCharPointer(splitname,2); 1116 PetscValidPointer(is,3); 1117 { 1118 PC_FieldSplit *jac = (PC_FieldSplit *) pc->data; 1119 PC_FieldSplitLink ilink = jac->head; 1120 PetscBool found; 1121 1122 *is = PETSC_NULL; 1123 while(ilink) { 1124 ierr = PetscStrcmp(ilink->splitname, splitname, &found);CHKERRQ(ierr); 1125 if (found) { 1126 *is = ilink->is; 1127 break; 1128 } 1129 ilink = ilink->next; 1130 } 1131 } 1132 PetscFunctionReturn(0); 1133 } 1134 1135 #undef __FUNCT__ 1136 #define __FUNCT__ "PCFieldSplitSetBlockSize" 1137 /*@ 1138 PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the 1139 fieldsplit preconditioner. If not set the matrix block size is used. 1140 1141 Logically Collective on PC 1142 1143 Input Parameters: 1144 + pc - the preconditioner context 1145 - bs - the block size 1146 1147 Level: intermediate 1148 1149 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields() 1150 1151 @*/ 1152 PetscErrorCode PCFieldSplitSetBlockSize(PC pc,PetscInt bs) 1153 { 1154 PetscErrorCode ierr; 1155 1156 PetscFunctionBegin; 1157 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1158 PetscValidLogicalCollectiveInt(pc,bs,2); 1159 ierr = PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));CHKERRQ(ierr); 1160 PetscFunctionReturn(0); 1161 } 1162 1163 #undef __FUNCT__ 1164 #define __FUNCT__ "PCFieldSplitGetSubKSP" 1165 /*@C 1166 PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits 1167 1168 Collective on KSP 1169 1170 Input Parameter: 1171 . pc - the preconditioner context 1172 1173 Output Parameters: 1174 + n - the number of splits 1175 - pc - the array of KSP contexts 1176 1177 Note: 1178 After PCFieldSplitGetSubKSP() the array of KSPs IS to be freed by the user 1179 (not the KSP just the array that contains them). 1180 1181 You must call KSPSetUp() before calling PCFieldSplitGetSubKSP(). 1182 1183 Level: advanced 1184 1185 .seealso: PCFIELDSPLIT 1186 @*/ 1187 PetscErrorCode PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[]) 1188 { 1189 PetscErrorCode ierr; 1190 1191 PetscFunctionBegin; 1192 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1193 if (n) PetscValidIntPointer(n,2); 1194 ierr = PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));CHKERRQ(ierr); 1195 PetscFunctionReturn(0); 1196 } 1197 1198 #undef __FUNCT__ 1199 #define __FUNCT__ "PCFieldSplitSchurPrecondition" 1200 /*@ 1201 PCFieldSplitSchurPrecondition - Indicates if the Schur complement is preconditioned by a preconditioner constructed by the 1202 A11 matrix. Otherwise no preconditioner is used. 1203 1204 Collective on PC 1205 1206 Input Parameters: 1207 + pc - the preconditioner context 1208 . ptype - which matrix to use for preconditioning the Schur complement, PC_FIELDSPLIT_SCHUR_PRE_DIAG (diag) is default 1209 - userpre - matrix to use for preconditioning, or PETSC_NULL 1210 1211 Options Database: 1212 . -pc_fieldsplit_schur_precondition <self,user,diag> default is diag 1213 1214 Notes: 1215 $ If ptype is 1216 $ user then the preconditioner for the Schur complement is generated by the provided matrix (pre argument 1217 $ to this function). 1218 $ diag then the preconditioner for the Schur complement is generated by the block diagonal part of the original 1219 $ matrix associated with the Schur complement (i.e. A11) 1220 $ self the preconditioner for the Schur complement is generated from the Schur complement matrix itself: 1221 $ The only preconditioner that currently works directly with the Schur complement matrix object is the PCLSC 1222 $ preconditioner 1223 1224 When solving a saddle point problem, where the A11 block is identically zero, using diag as the ptype only makes sense 1225 with the additional option -fieldsplit_1_pc_type none. Usually for saddle point problems one would use a ptype of self and 1226 -fieldsplit_1_pc_type lsc which uses the least squares commutator compute a preconditioner for the Schur complement. 1227 1228 Developer Notes: This is a terrible name, gives no good indication of what the function does and should also have Set in 1229 the name since it sets a proceedure to use. 1230 1231 Level: intermediate 1232 1233 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCLSC 1234 1235 @*/ 1236 PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1237 { 1238 PetscErrorCode ierr; 1239 1240 PetscFunctionBegin; 1241 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1242 ierr = PetscTryMethod(pc,"PCFieldSplitSchurPrecondition_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));CHKERRQ(ierr); 1243 PetscFunctionReturn(0); 1244 } 1245 1246 EXTERN_C_BEGIN 1247 #undef __FUNCT__ 1248 #define __FUNCT__ "PCFieldSplitSchurPrecondition_FieldSplit" 1249 PetscErrorCode PCFieldSplitSchurPrecondition_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1250 { 1251 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1252 PetscErrorCode ierr; 1253 1254 PetscFunctionBegin; 1255 jac->schurpre = ptype; 1256 if (pre) { 1257 ierr = MatDestroy(&jac->schur_user);CHKERRQ(ierr); 1258 jac->schur_user = pre; 1259 ierr = PetscObjectReference((PetscObject)jac->schur_user);CHKERRQ(ierr); 1260 } 1261 PetscFunctionReturn(0); 1262 } 1263 EXTERN_C_END 1264 1265 #undef __FUNCT__ 1266 #define __FUNCT__ "PCFieldSplitGetSchurBlocks" 1267 /*@C 1268 PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement 1269 1270 Collective on KSP 1271 1272 Input Parameter: 1273 . pc - the preconditioner context 1274 1275 Output Parameters: 1276 + A00 - the (0,0) block 1277 . A01 - the (0,1) block 1278 . A10 - the (1,0) block 1279 - A11 - the (1,1) block 1280 1281 Level: advanced 1282 1283 .seealso: PCFIELDSPLIT 1284 @*/ 1285 PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11) 1286 { 1287 PC_FieldSplit *jac = (PC_FieldSplit *) pc->data; 1288 1289 PetscFunctionBegin; 1290 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1291 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach."); 1292 if (A00) *A00 = jac->pmat[0]; 1293 if (A01) *A01 = jac->B; 1294 if (A10) *A10 = jac->C; 1295 if (A11) *A11 = jac->pmat[1]; 1296 PetscFunctionReturn(0); 1297 } 1298 1299 EXTERN_C_BEGIN 1300 #undef __FUNCT__ 1301 #define __FUNCT__ "PCFieldSplitSetType_FieldSplit" 1302 PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type) 1303 { 1304 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1305 PetscErrorCode ierr; 1306 1307 PetscFunctionBegin; 1308 jac->type = type; 1309 if (type == PC_COMPOSITE_SCHUR) { 1310 pc->ops->apply = PCApply_FieldSplit_Schur; 1311 pc->ops->view = PCView_FieldSplit_Schur; 1312 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit_Schur",PCFieldSplitGetSubKSP_FieldSplit_Schur);CHKERRQ(ierr); 1313 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSchurPrecondition_C","PCFieldSplitSchurPrecondition_FieldSplit",PCFieldSplitSchurPrecondition_FieldSplit);CHKERRQ(ierr); 1314 1315 } else { 1316 pc->ops->apply = PCApply_FieldSplit; 1317 pc->ops->view = PCView_FieldSplit; 1318 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit",PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 1319 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSchurPrecondition_C","",0);CHKERRQ(ierr); 1320 } 1321 PetscFunctionReturn(0); 1322 } 1323 EXTERN_C_END 1324 1325 EXTERN_C_BEGIN 1326 #undef __FUNCT__ 1327 #define __FUNCT__ "PCFieldSplitSetBlockSize_FieldSplit" 1328 PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs) 1329 { 1330 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1331 1332 PetscFunctionBegin; 1333 if (bs < 1) SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs); 1334 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); 1335 jac->bs = bs; 1336 PetscFunctionReturn(0); 1337 } 1338 EXTERN_C_END 1339 1340 #undef __FUNCT__ 1341 #define __FUNCT__ "PCFieldSplitSetType" 1342 /*@ 1343 PCFieldSplitSetType - Sets the type of fieldsplit preconditioner. 1344 1345 Collective on PC 1346 1347 Input Parameter: 1348 . pc - the preconditioner context 1349 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 1350 1351 Options Database Key: 1352 . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type 1353 1354 Level: Developer 1355 1356 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 1357 1358 .seealso: PCCompositeSetType() 1359 1360 @*/ 1361 PetscErrorCode PCFieldSplitSetType(PC pc,PCCompositeType type) 1362 { 1363 PetscErrorCode ierr; 1364 1365 PetscFunctionBegin; 1366 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1367 ierr = PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));CHKERRQ(ierr); 1368 PetscFunctionReturn(0); 1369 } 1370 1371 /* -------------------------------------------------------------------------------------*/ 1372 /*MC 1373 PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual 1374 fields or groups of fields. See the users manual section "Solving Block Matrices" for more details. 1375 1376 To set options on the solvers for each block append -fieldsplit_ to all the PC 1377 options database keys. For example, -fieldsplit_pc_type ilu -fieldsplit_pc_factor_levels 1 1378 1379 To set the options on the solvers separate for each block call PCFieldSplitGetSubKSP() 1380 and set the options directly on the resulting KSP object 1381 1382 Level: intermediate 1383 1384 Options Database Keys: 1385 + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split 1386 . -pc_fieldsplit_default - automatically add any fields to additional splits that have not 1387 been supplied explicitly by -pc_fieldsplit_%d_fields 1388 . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields) 1389 . -pc_fieldsplit_type <additive,multiplicative,schur,symmetric_multiplicative> 1390 . -pc_fieldsplit_schur_precondition <true,false> default is true 1391 . -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero or negative diagonal and uses Schur complement with no preconditioner as the solver 1392 1393 - Options prefix for inner solvers when using Schur complement preconditioner are -fieldsplit_0_ and -fieldsplit_1_ 1394 for all other solvers they are -fieldsplit_%d_ for the dth field, use -fieldsplit_ for all fields 1395 1396 Notes: use PCFieldSplitSetFields() to set fields defined by "strided" entries and PCFieldSplitSetIS() 1397 to define a field by an arbitrary collection of entries. 1398 1399 If no fields are set the default is used. The fields are defined by entries strided by bs, 1400 beginning at 0 then 1, etc to bs-1. The block size can be set with PCFieldSplitSetBlockSize(), 1401 if this is not called the block size defaults to the blocksize of the second matrix passed 1402 to KSPSetOperators()/PCSetOperators(). 1403 1404 For the Schur complement preconditioner if J = ( A00 A01 ) 1405 ( A10 A11 ) 1406 the preconditioner is 1407 (I -B ksp(A00)) ( inv(A00) 0 (I 0 ) 1408 (0 I ) ( 0 ksp(S) ) (-A10 ksp(A00) I ) 1409 where the action of inv(A00) is applied using the KSP solver with prefix -fieldsplit_0_. The action of 1410 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). 1411 For PCFieldSplitGetKSP() when field number is 1412 0 it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default 1413 A11 is used to construct a preconditioner for S, use PCFieldSplitSchurPrecondition() to turn on or off this 1414 option. You can use the preconditioner PCLSC to precondition the Schur complement with -fieldsplit_1_pc_type lsc 1415 1416 If only one set of indices (one IS) is provided with PCFieldSplitSetIS() then the complement of that IS 1417 is used automatically for a second block. 1418 1419 The fieldsplit preconditioner cannot currently be used with the BAIJ or SBAIJ data formats if the blocksize is larger than 1. 1420 Generally it should be used with the AIJ format. 1421 1422 The forms of these preconditioners are closely related if not identical to forms derived as "Distributive Iterations", see, 1423 for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling. Note that one can also use PCFIELDSPLIT 1424 inside a smoother resulting in "Distributive Smoothers". 1425 1426 Concepts: physics based preconditioners, block preconditioners 1427 1428 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners, PCLSC, 1429 PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSchurPrecondition() 1430 M*/ 1431 1432 EXTERN_C_BEGIN 1433 #undef __FUNCT__ 1434 #define __FUNCT__ "PCCreate_FieldSplit" 1435 PetscErrorCode PCCreate_FieldSplit(PC pc) 1436 { 1437 PetscErrorCode ierr; 1438 PC_FieldSplit *jac; 1439 1440 PetscFunctionBegin; 1441 ierr = PetscNewLog(pc,PC_FieldSplit,&jac);CHKERRQ(ierr); 1442 jac->bs = -1; 1443 jac->nsplits = 0; 1444 jac->type = PC_COMPOSITE_MULTIPLICATIVE; 1445 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */ 1446 jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACTORIZATION_FULL; 1447 1448 pc->data = (void*)jac; 1449 1450 pc->ops->apply = PCApply_FieldSplit; 1451 pc->ops->applytranspose = PCApplyTranspose_FieldSplit; 1452 pc->ops->setup = PCSetUp_FieldSplit; 1453 pc->ops->reset = PCReset_FieldSplit; 1454 pc->ops->destroy = PCDestroy_FieldSplit; 1455 pc->ops->setfromoptions = PCSetFromOptions_FieldSplit; 1456 pc->ops->view = PCView_FieldSplit; 1457 pc->ops->applyrichardson = 0; 1458 1459 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit", 1460 PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 1461 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetFields_C","PCFieldSplitSetFields_FieldSplit", 1462 PCFieldSplitSetFields_FieldSplit);CHKERRQ(ierr); 1463 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetIS_C","PCFieldSplitSetIS_FieldSplit", 1464 PCFieldSplitSetIS_FieldSplit);CHKERRQ(ierr); 1465 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetType_C","PCFieldSplitSetType_FieldSplit", 1466 PCFieldSplitSetType_FieldSplit);CHKERRQ(ierr); 1467 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetBlockSize_C","PCFieldSplitSetBlockSize_FieldSplit", 1468 PCFieldSplitSetBlockSize_FieldSplit);CHKERRQ(ierr); 1469 PetscFunctionReturn(0); 1470 } 1471 EXTERN_C_END 1472 1473 1474 1475