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