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