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