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