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