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; 669 670 PetscFunctionBegin; 671 CHKMEMQ; 672 ierr = VecSetBlockSize(x,jac->bs);CHKERRQ(ierr); 673 ierr = VecSetBlockSize(y,jac->bs);CHKERRQ(ierr); 674 675 if (jac->type == PC_COMPOSITE_ADDITIVE) { 676 if (jac->defaultsplit) { 677 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 678 while (ilink) { 679 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 680 ilink = ilink->next; 681 } 682 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 683 } else { 684 ierr = VecSet(y,0.0);CHKERRQ(ierr); 685 while (ilink) { 686 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 687 ilink = ilink->next; 688 } 689 } 690 } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 691 if (!jac->w1) { 692 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 693 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 694 } 695 ierr = VecSet(y,0.0);CHKERRQ(ierr); 696 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 697 cnt = 1; 698 while (ilink->next) { 699 ilink = ilink->next; 700 /* compute the residual only over the part of the vector needed */ 701 ierr = MatMult(jac->Afield[cnt++],y,ilink->x);CHKERRQ(ierr); 702 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 703 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 704 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 705 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 706 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 707 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 708 } 709 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 710 cnt -= 2; 711 while (ilink->previous) { 712 ilink = ilink->previous; 713 /* compute the residual only over the part of the vector needed */ 714 ierr = MatMult(jac->Afield[cnt--],y,ilink->x);CHKERRQ(ierr); 715 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 716 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 717 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 718 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 719 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 720 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 721 } 722 } 723 } else SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Unsupported or unknown composition",(int) jac->type); 724 CHKMEMQ; 725 PetscFunctionReturn(0); 726 } 727 728 #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \ 729 (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 730 VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 731 KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) || \ 732 VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \ 733 VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE)) 734 735 #undef __FUNCT__ 736 #define __FUNCT__ "PCApplyTranspose_FieldSplit" 737 static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y) 738 { 739 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 740 PetscErrorCode ierr; 741 PC_FieldSplitLink ilink = jac->head; 742 743 PetscFunctionBegin; 744 CHKMEMQ; 745 ierr = VecSetBlockSize(x,jac->bs);CHKERRQ(ierr); 746 ierr = VecSetBlockSize(y,jac->bs);CHKERRQ(ierr); 747 748 if (jac->type == PC_COMPOSITE_ADDITIVE) { 749 if (jac->defaultsplit) { 750 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 751 while (ilink) { 752 ierr = KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 753 ilink = ilink->next; 754 } 755 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 756 } else { 757 ierr = VecSet(y,0.0);CHKERRQ(ierr); 758 while (ilink) { 759 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 760 ilink = ilink->next; 761 } 762 } 763 } else { 764 if (!jac->w1) { 765 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 766 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 767 } 768 ierr = VecSet(y,0.0);CHKERRQ(ierr); 769 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 770 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 771 while (ilink->next) { 772 ilink = ilink->next; 773 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 774 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 775 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 776 } 777 while (ilink->previous) { 778 ilink = ilink->previous; 779 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 780 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 781 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 782 } 783 } else { 784 while (ilink->next) { /* get to last entry in linked list */ 785 ilink = ilink->next; 786 } 787 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 788 while (ilink->previous) { 789 ilink = ilink->previous; 790 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 791 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 792 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 793 } 794 } 795 } 796 CHKMEMQ; 797 PetscFunctionReturn(0); 798 } 799 800 #undef __FUNCT__ 801 #define __FUNCT__ "PCReset_FieldSplit" 802 static PetscErrorCode PCReset_FieldSplit(PC pc) 803 { 804 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 805 PetscErrorCode ierr; 806 PC_FieldSplitLink ilink = jac->head,next; 807 808 PetscFunctionBegin; 809 while (ilink) { 810 ierr = KSPReset(ilink->ksp);CHKERRQ(ierr); 811 ierr = VecDestroy(&ilink->x);CHKERRQ(ierr); 812 ierr = VecDestroy(&ilink->y);CHKERRQ(ierr); 813 ierr = VecScatterDestroy(&ilink->sctx);CHKERRQ(ierr); 814 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 815 ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr); 816 next = ilink->next; 817 ilink = next; 818 } 819 ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr); 820 if (jac->mat && jac->mat != jac->pmat) { 821 ierr = MatDestroyMatrices(jac->nsplits,&jac->mat);CHKERRQ(ierr); 822 } else if (jac->mat) { 823 jac->mat = PETSC_NULL; 824 } 825 if (jac->pmat) {ierr = MatDestroyMatrices(jac->nsplits,&jac->pmat);CHKERRQ(ierr);} 826 if (jac->Afield) {ierr = MatDestroyMatrices(jac->nsplits,&jac->Afield);CHKERRQ(ierr);} 827 ierr = VecDestroy(&jac->w1);CHKERRQ(ierr); 828 ierr = VecDestroy(&jac->w2);CHKERRQ(ierr); 829 ierr = MatDestroy(&jac->schur);CHKERRQ(ierr); 830 ierr = MatDestroy(&jac->schur_user);CHKERRQ(ierr); 831 ierr = KSPDestroy(&jac->kspschur);CHKERRQ(ierr); 832 ierr = MatDestroy(&jac->B);CHKERRQ(ierr); 833 ierr = MatDestroy(&jac->C);CHKERRQ(ierr); 834 jac->reset = PETSC_TRUE; 835 PetscFunctionReturn(0); 836 } 837 838 #undef __FUNCT__ 839 #define __FUNCT__ "PCDestroy_FieldSplit" 840 static PetscErrorCode PCDestroy_FieldSplit(PC pc) 841 { 842 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 843 PetscErrorCode ierr; 844 PC_FieldSplitLink ilink = jac->head,next; 845 846 PetscFunctionBegin; 847 ierr = PCReset_FieldSplit(pc);CHKERRQ(ierr); 848 while (ilink) { 849 ierr = KSPDestroy(&ilink->ksp);CHKERRQ(ierr); 850 next = ilink->next; 851 ierr = PetscFree(ilink->splitname);CHKERRQ(ierr); 852 ierr = PetscFree(ilink->fields);CHKERRQ(ierr); 853 ierr = PetscFree(ilink->fields_col);CHKERRQ(ierr); 854 ierr = PetscFree(ilink);CHKERRQ(ierr); 855 ilink = next; 856 } 857 ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr); 858 ierr = PetscFree(pc->data);CHKERRQ(ierr); 859 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","",PETSC_NULL);CHKERRQ(ierr); 860 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetFields_C","",PETSC_NULL);CHKERRQ(ierr); 861 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetIS_C","",PETSC_NULL);CHKERRQ(ierr); 862 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetType_C","",PETSC_NULL);CHKERRQ(ierr); 863 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetBlockSize_C","",PETSC_NULL);CHKERRQ(ierr); 864 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSchurPrecondition_C","",PETSC_NULL);CHKERRQ(ierr); 865 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetSchurFactType_C","",PETSC_NULL);CHKERRQ(ierr); 866 PetscFunctionReturn(0); 867 } 868 869 #undef __FUNCT__ 870 #define __FUNCT__ "PCSetFromOptions_FieldSplit" 871 static PetscErrorCode PCSetFromOptions_FieldSplit(PC pc) 872 { 873 PetscErrorCode ierr; 874 PetscInt bs; 875 PetscBool flg,stokes = PETSC_FALSE; 876 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 877 PCCompositeType ctype; 878 DM ddm; 879 char ddm_name[1025]; 880 881 PetscFunctionBegin; 882 ierr = PetscOptionsHead("FieldSplit options");CHKERRQ(ierr); 883 if(pc->dm) { 884 /* Allow the user to request a decomposition DM by name */ 885 ierr = PetscStrncpy(ddm_name, "", 1024); CHKERRQ(ierr); 886 ierr = PetscOptionsString("-pc_fieldsplit_decomposition_dm", "Name of the DM defining the composition", "PCSetDM", ddm_name, ddm_name,1024,&flg); CHKERRQ(ierr); 887 if(flg) { 888 ierr = DMCreateDecompositionDM(pc->dm, ddm_name, &ddm); CHKERRQ(ierr); 889 if(!ddm) { 890 SETERRQ1(((PetscObject)pc)->comm, PETSC_ERR_ARG_WRONGSTATE, "Uknown DM decomposition name %s", ddm_name); 891 } 892 ierr = PetscInfo(pc,"Using decomposition DM defined using options database\n");CHKERRQ(ierr); 893 ierr = PCSetDM(pc,ddm); CHKERRQ(ierr); 894 } 895 } 896 ierr = PetscOptionsBool("-pc_fieldsplit_real_diagonal","Use diagonal blocks of the operator","PCFieldSplitSetRealDiagonal",jac->realdiagonal,&jac->realdiagonal,PETSC_NULL);CHKERRQ(ierr); 897 ierr = PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);CHKERRQ(ierr); 898 if (flg) { 899 ierr = PCFieldSplitSetBlockSize(pc,bs);CHKERRQ(ierr); 900 } 901 902 ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,PETSC_NULL);CHKERRQ(ierr); 903 if (stokes) { 904 ierr = PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);CHKERRQ(ierr); 905 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_SELF; 906 } 907 908 ierr = PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg);CHKERRQ(ierr); 909 if (flg) { 910 ierr = PCFieldSplitSetType(pc,ctype);CHKERRQ(ierr); 911 } 912 /* Only setup fields once */ 913 if ((jac->bs > 0) && (jac->nsplits == 0)) { 914 /* only allow user to set fields from command line if bs is already known. 915 otherwise user can set them in PCFieldSplitSetDefaults() */ 916 ierr = PCFieldSplitSetRuntimeSplits_Private(pc);CHKERRQ(ierr); 917 if (jac->splitdefined) {ierr = PetscInfo(pc,"Splits defined using the options database\n");CHKERRQ(ierr);} 918 } 919 if (jac->type == PC_COMPOSITE_SCHUR) { 920 ierr = PetscOptionsGetEnum(((PetscObject)pc)->prefix,"-pc_fieldsplit_schur_factorization_type",PCFieldSplitSchurFactTypes,(PetscEnum*)&jac->schurfactorization,&flg);CHKERRQ(ierr); 921 if (flg) {ierr = PetscInfo(pc,"Deprecated use of -pc_fieldsplit_schur_factorization_type\n");CHKERRQ(ierr);} 922 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); 923 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); 924 } 925 ierr = PetscOptionsTail();CHKERRQ(ierr); 926 PetscFunctionReturn(0); 927 } 928 929 /*------------------------------------------------------------------------------------*/ 930 931 EXTERN_C_BEGIN 932 #undef __FUNCT__ 933 #define __FUNCT__ "PCFieldSplitSetFields_FieldSplit" 934 PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 935 { 936 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 937 PetscErrorCode ierr; 938 PC_FieldSplitLink ilink,next = jac->head; 939 char prefix[128]; 940 PetscInt i; 941 942 PetscFunctionBegin; 943 if (jac->splitdefined) { 944 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 945 PetscFunctionReturn(0); 946 } 947 for (i=0; i<n; i++) { 948 if (fields[i] >= jac->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs); 949 if (fields[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]); 950 } 951 ierr = PetscNew(struct _PC_FieldSplitLink,&ilink);CHKERRQ(ierr); 952 if (splitname) { 953 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 954 } else { 955 ierr = PetscMalloc(3*sizeof(char),&ilink->splitname);CHKERRQ(ierr); 956 ierr = PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);CHKERRQ(ierr); 957 } 958 ierr = PetscMalloc(n*sizeof(PetscInt),&ilink->fields);CHKERRQ(ierr); 959 ierr = PetscMemcpy(ilink->fields,fields,n*sizeof(PetscInt));CHKERRQ(ierr); 960 ierr = PetscMalloc(n*sizeof(PetscInt),&ilink->fields_col);CHKERRQ(ierr); 961 ierr = PetscMemcpy(ilink->fields_col,fields_col,n*sizeof(PetscInt));CHKERRQ(ierr); 962 ilink->nfields = n; 963 ilink->next = PETSC_NULL; 964 ierr = KSPCreate(((PetscObject)pc)->comm,&ilink->ksp);CHKERRQ(ierr); 965 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 966 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 967 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 968 969 ierr = PetscSNPrintf(prefix,sizeof prefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",ilink->splitname);CHKERRQ(ierr); 970 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 971 972 if (!next) { 973 jac->head = ilink; 974 ilink->previous = PETSC_NULL; 975 } else { 976 while (next->next) { 977 next = next->next; 978 } 979 next->next = ilink; 980 ilink->previous = next; 981 } 982 jac->nsplits++; 983 PetscFunctionReturn(0); 984 } 985 EXTERN_C_END 986 987 EXTERN_C_BEGIN 988 #undef __FUNCT__ 989 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit_Schur" 990 PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp) 991 { 992 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 993 PetscErrorCode ierr; 994 995 PetscFunctionBegin; 996 ierr = PetscMalloc(jac->nsplits*sizeof(KSP),subksp);CHKERRQ(ierr); 997 ierr = MatSchurComplementGetKSP(jac->schur,*subksp);CHKERRQ(ierr); 998 (*subksp)[1] = jac->kspschur; 999 if (n) *n = jac->nsplits; 1000 PetscFunctionReturn(0); 1001 } 1002 EXTERN_C_END 1003 1004 EXTERN_C_BEGIN 1005 #undef __FUNCT__ 1006 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit" 1007 PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp) 1008 { 1009 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1010 PetscErrorCode ierr; 1011 PetscInt cnt = 0; 1012 PC_FieldSplitLink ilink = jac->head; 1013 1014 PetscFunctionBegin; 1015 ierr = PetscMalloc(jac->nsplits*sizeof(KSP),subksp);CHKERRQ(ierr); 1016 while (ilink) { 1017 (*subksp)[cnt++] = ilink->ksp; 1018 ilink = ilink->next; 1019 } 1020 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); 1021 if (n) *n = jac->nsplits; 1022 PetscFunctionReturn(0); 1023 } 1024 EXTERN_C_END 1025 1026 EXTERN_C_BEGIN 1027 #undef __FUNCT__ 1028 #define __FUNCT__ "PCFieldSplitSetIS_FieldSplit" 1029 PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is) 1030 { 1031 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1032 PetscErrorCode ierr; 1033 PC_FieldSplitLink ilink, next = jac->head; 1034 char prefix[128]; 1035 1036 PetscFunctionBegin; 1037 if (jac->splitdefined) { 1038 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 1039 PetscFunctionReturn(0); 1040 } 1041 ierr = PetscNew(struct _PC_FieldSplitLink,&ilink);CHKERRQ(ierr); 1042 if (splitname) { 1043 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 1044 } else { 1045 ierr = PetscMalloc(8*sizeof(char),&ilink->splitname);CHKERRQ(ierr); 1046 ierr = PetscSNPrintf(ilink->splitname,7,"%D",jac->nsplits);CHKERRQ(ierr); 1047 } 1048 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 1049 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 1050 ilink->is = is; 1051 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 1052 ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr); 1053 ilink->is_col = is; 1054 ilink->next = PETSC_NULL; 1055 ierr = KSPCreate(((PetscObject)pc)->comm,&ilink->ksp);CHKERRQ(ierr); 1056 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 1057 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 1058 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 1059 1060 ierr = PetscSNPrintf(prefix,sizeof prefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",ilink->splitname);CHKERRQ(ierr); 1061 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 1062 1063 if (!next) { 1064 jac->head = ilink; 1065 ilink->previous = PETSC_NULL; 1066 } else { 1067 while (next->next) { 1068 next = next->next; 1069 } 1070 next->next = ilink; 1071 ilink->previous = next; 1072 } 1073 jac->nsplits++; 1074 1075 PetscFunctionReturn(0); 1076 } 1077 EXTERN_C_END 1078 1079 #undef __FUNCT__ 1080 #define __FUNCT__ "PCFieldSplitSetFields" 1081 /*@ 1082 PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner 1083 1084 Logically Collective on PC 1085 1086 Input Parameters: 1087 + pc - the preconditioner context 1088 . splitname - name of this split, if PETSC_NULL the number of the split is used 1089 . n - the number of fields in this split 1090 - fields - the fields in this split 1091 1092 Level: intermediate 1093 1094 Notes: Use PCFieldSplitSetIS() to set a completely general set of indices as a field. 1095 1096 The PCFieldSplitSetFields() is for defining fields as a strided blocks. For example, if the block 1097 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 1098 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x.... 1099 where the numbered entries indicate what is in the field. 1100 1101 This function is called once per split (it creates a new split each time). Solve options 1102 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1103 1104 Developer Note: This routine does not actually create the IS representing the split, that is delayed 1105 until PCSetUp_FieldSplit(), because information about the vector/matrix layouts may not be 1106 available when this routine is called. 1107 1108 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize(), PCFieldSplitSetIS() 1109 1110 @*/ 1111 PetscErrorCode PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 1112 { 1113 PetscErrorCode ierr; 1114 1115 PetscFunctionBegin; 1116 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1117 PetscValidCharPointer(splitname,2); 1118 if (n < 1) SETERRQ2(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Provided number of fields %D in split \"%s\" not positive",n,splitname); 1119 PetscValidIntPointer(fields,3); 1120 ierr = PetscTryMethod(pc,"PCFieldSplitSetFields_C",(PC,const char[],PetscInt,const PetscInt *,const PetscInt *),(pc,splitname,n,fields,fields_col));CHKERRQ(ierr); 1121 PetscFunctionReturn(0); 1122 } 1123 1124 #undef __FUNCT__ 1125 #define __FUNCT__ "PCFieldSplitSetIS" 1126 /*@ 1127 PCFieldSplitSetIS - Sets the exact elements for field 1128 1129 Logically Collective on PC 1130 1131 Input Parameters: 1132 + pc - the preconditioner context 1133 . splitname - name of this split, if PETSC_NULL the number of the split is used 1134 - is - the index set that defines the vector elements in this field 1135 1136 1137 Notes: 1138 Use PCFieldSplitSetFields(), for fields defined by strided types. 1139 1140 This function is called once per split (it creates a new split each time). Solve options 1141 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1142 1143 Level: intermediate 1144 1145 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize() 1146 1147 @*/ 1148 PetscErrorCode PCFieldSplitSetIS(PC pc,const char splitname[],IS is) 1149 { 1150 PetscErrorCode ierr; 1151 1152 PetscFunctionBegin; 1153 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1154 if (splitname) PetscValidCharPointer(splitname,2); 1155 PetscValidHeaderSpecific(is,IS_CLASSID,3); 1156 ierr = PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));CHKERRQ(ierr); 1157 PetscFunctionReturn(0); 1158 } 1159 1160 #undef __FUNCT__ 1161 #define __FUNCT__ "PCFieldSplitGetIS" 1162 /*@ 1163 PCFieldSplitGetIS - Retrieves the elements for a field as an IS 1164 1165 Logically Collective on PC 1166 1167 Input Parameters: 1168 + pc - the preconditioner context 1169 - splitname - name of this split 1170 1171 Output Parameter: 1172 - is - the index set that defines the vector elements in this field, or PETSC_NULL if the field is not found 1173 1174 Level: intermediate 1175 1176 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetIS() 1177 1178 @*/ 1179 PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is) 1180 { 1181 PetscErrorCode ierr; 1182 1183 PetscFunctionBegin; 1184 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1185 PetscValidCharPointer(splitname,2); 1186 PetscValidPointer(is,3); 1187 { 1188 PC_FieldSplit *jac = (PC_FieldSplit *) pc->data; 1189 PC_FieldSplitLink ilink = jac->head; 1190 PetscBool found; 1191 1192 *is = PETSC_NULL; 1193 while(ilink) { 1194 ierr = PetscStrcmp(ilink->splitname, splitname, &found);CHKERRQ(ierr); 1195 if (found) { 1196 *is = ilink->is; 1197 break; 1198 } 1199 ilink = ilink->next; 1200 } 1201 } 1202 PetscFunctionReturn(0); 1203 } 1204 1205 #undef __FUNCT__ 1206 #define __FUNCT__ "PCFieldSplitSetBlockSize" 1207 /*@ 1208 PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the 1209 fieldsplit preconditioner. If not set the matrix block size is used. 1210 1211 Logically Collective on PC 1212 1213 Input Parameters: 1214 + pc - the preconditioner context 1215 - bs - the block size 1216 1217 Level: intermediate 1218 1219 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields() 1220 1221 @*/ 1222 PetscErrorCode PCFieldSplitSetBlockSize(PC pc,PetscInt bs) 1223 { 1224 PetscErrorCode ierr; 1225 1226 PetscFunctionBegin; 1227 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1228 PetscValidLogicalCollectiveInt(pc,bs,2); 1229 ierr = PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));CHKERRQ(ierr); 1230 PetscFunctionReturn(0); 1231 } 1232 1233 #undef __FUNCT__ 1234 #define __FUNCT__ "PCFieldSplitGetSubKSP" 1235 /*@C 1236 PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits 1237 1238 Collective on KSP 1239 1240 Input Parameter: 1241 . pc - the preconditioner context 1242 1243 Output Parameters: 1244 + n - the number of splits 1245 - pc - the array of KSP contexts 1246 1247 Note: 1248 After PCFieldSplitGetSubKSP() the array of KSPs IS to be freed by the user 1249 (not the KSP just the array that contains them). 1250 1251 You must call KSPSetUp() before calling PCFieldSplitGetSubKSP(). 1252 1253 Level: advanced 1254 1255 .seealso: PCFIELDSPLIT 1256 @*/ 1257 PetscErrorCode PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[]) 1258 { 1259 PetscErrorCode ierr; 1260 1261 PetscFunctionBegin; 1262 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1263 if (n) PetscValidIntPointer(n,2); 1264 ierr = PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));CHKERRQ(ierr); 1265 PetscFunctionReturn(0); 1266 } 1267 1268 #undef __FUNCT__ 1269 #define __FUNCT__ "PCFieldSplitSchurPrecondition" 1270 /*@ 1271 PCFieldSplitSchurPrecondition - Indicates if the Schur complement is preconditioned by a preconditioner constructed by the 1272 A11 matrix. Otherwise no preconditioner is used. 1273 1274 Collective on PC 1275 1276 Input Parameters: 1277 + pc - the preconditioner context 1278 . ptype - which matrix to use for preconditioning the Schur complement, PC_FIELDSPLIT_SCHUR_PRE_DIAG (diag) is default 1279 - userpre - matrix to use for preconditioning, or PETSC_NULL 1280 1281 Options Database: 1282 . -pc_fieldsplit_schur_precondition <self,user,diag> default is diag 1283 1284 Notes: 1285 $ If ptype is 1286 $ user then the preconditioner for the Schur complement is generated by the provided matrix (pre argument 1287 $ to this function). 1288 $ diag then the preconditioner for the Schur complement is generated by the block diagonal part of the original 1289 $ matrix associated with the Schur complement (i.e. A11) 1290 $ self the preconditioner for the Schur complement is generated from the Schur complement matrix itself: 1291 $ The only preconditioner that currently works directly with the Schur complement matrix object is the PCLSC 1292 $ preconditioner 1293 1294 When solving a saddle point problem, where the A11 block is identically zero, using diag as the ptype only makes sense 1295 with the additional option -fieldsplit_1_pc_type none. Usually for saddle point problems one would use a ptype of self and 1296 -fieldsplit_1_pc_type lsc which uses the least squares commutator compute a preconditioner for the Schur complement. 1297 1298 Developer Notes: This is a terrible name, gives no good indication of what the function does and should also have Set in 1299 the name since it sets a proceedure to use. 1300 1301 Level: intermediate 1302 1303 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCLSC 1304 1305 @*/ 1306 PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1307 { 1308 PetscErrorCode ierr; 1309 1310 PetscFunctionBegin; 1311 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1312 ierr = PetscTryMethod(pc,"PCFieldSplitSchurPrecondition_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));CHKERRQ(ierr); 1313 PetscFunctionReturn(0); 1314 } 1315 1316 EXTERN_C_BEGIN 1317 #undef __FUNCT__ 1318 #define __FUNCT__ "PCFieldSplitSchurPrecondition_FieldSplit" 1319 PetscErrorCode PCFieldSplitSchurPrecondition_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1320 { 1321 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1322 PetscErrorCode ierr; 1323 1324 PetscFunctionBegin; 1325 jac->schurpre = ptype; 1326 if (pre) { 1327 ierr = MatDestroy(&jac->schur_user);CHKERRQ(ierr); 1328 jac->schur_user = pre; 1329 ierr = PetscObjectReference((PetscObject)jac->schur_user);CHKERRQ(ierr); 1330 } 1331 PetscFunctionReturn(0); 1332 } 1333 EXTERN_C_END 1334 1335 #undef __FUNCT__ 1336 #define __FUNCT__ "PCFieldSplitSetSchurFactType" 1337 /*@ 1338 PCFieldSplitSetSchurFactType - sets which blocks of the approximate block factorization to retain 1339 1340 Collective on PC 1341 1342 Input Parameters: 1343 + pc - the preconditioner context 1344 - ftype - which blocks of factorization to retain, PC_FIELDSPLIT_SCHUR_FACT_FULL is default 1345 1346 Options Database: 1347 . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> default is full 1348 1349 1350 Level: intermediate 1351 1352 Notes: 1353 The FULL factorization is 1354 1355 $ (A B) = (1 0) (A 0) (1 Ainv*B) 1356 $ (C D) (C*Ainv 1) (0 S) (0 1 ) 1357 1358 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, 1359 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). 1360 1361 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 1362 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 1363 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, 1364 the preconditioned operator has three distinct nonzero eigenvalues and minimal polynomial of degree at most 4, so KSPGMRES converges in at most 4 iterations. 1365 1366 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 1367 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). 1368 1369 References: 1370 Murphy, Golub, and Wathen, A note on preconditioning indefinite linear systems, SIAM J. Sci. Comput., 21 (2000) pp. 1969-1972. 1371 1372 Ipsen, A note on preconditioning nonsymmetric matrices, SIAM J. Sci. Comput., 23 (2001), pp. 1050-1051. 1373 1374 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType 1375 @*/ 1376 PetscErrorCode PCFieldSplitSetSchurFactType(PC pc,PCFieldSplitSchurFactType ftype) 1377 { 1378 PetscErrorCode ierr; 1379 1380 PetscFunctionBegin; 1381 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1382 ierr = PetscTryMethod(pc,"PCFieldSplitSetSchurFactType_C",(PC,PCFieldSplitSchurFactType),(pc,ftype));CHKERRQ(ierr); 1383 PetscFunctionReturn(0); 1384 } 1385 1386 #undef __FUNCT__ 1387 #define __FUNCT__ "PCFieldSplitSetSchurFactType_FieldSplit" 1388 PETSC_EXTERN_C PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc,PCFieldSplitSchurFactType ftype) 1389 { 1390 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1391 1392 PetscFunctionBegin; 1393 jac->schurfactorization = ftype; 1394 PetscFunctionReturn(0); 1395 } 1396 1397 #undef __FUNCT__ 1398 #define __FUNCT__ "PCFieldSplitGetSchurBlocks" 1399 /*@C 1400 PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement 1401 1402 Collective on KSP 1403 1404 Input Parameter: 1405 . pc - the preconditioner context 1406 1407 Output Parameters: 1408 + A00 - the (0,0) block 1409 . A01 - the (0,1) block 1410 . A10 - the (1,0) block 1411 - A11 - the (1,1) block 1412 1413 Level: advanced 1414 1415 .seealso: PCFIELDSPLIT 1416 @*/ 1417 PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11) 1418 { 1419 PC_FieldSplit *jac = (PC_FieldSplit *) pc->data; 1420 1421 PetscFunctionBegin; 1422 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1423 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach."); 1424 if (A00) *A00 = jac->pmat[0]; 1425 if (A01) *A01 = jac->B; 1426 if (A10) *A10 = jac->C; 1427 if (A11) *A11 = jac->pmat[1]; 1428 PetscFunctionReturn(0); 1429 } 1430 1431 EXTERN_C_BEGIN 1432 #undef __FUNCT__ 1433 #define __FUNCT__ "PCFieldSplitSetType_FieldSplit" 1434 PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type) 1435 { 1436 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1437 PetscErrorCode ierr; 1438 1439 PetscFunctionBegin; 1440 jac->type = type; 1441 if (type == PC_COMPOSITE_SCHUR) { 1442 pc->ops->apply = PCApply_FieldSplit_Schur; 1443 pc->ops->view = PCView_FieldSplit_Schur; 1444 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit_Schur",PCFieldSplitGetSubKSP_FieldSplit_Schur);CHKERRQ(ierr); 1445 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSchurPrecondition_C","PCFieldSplitSchurPrecondition_FieldSplit",PCFieldSplitSchurPrecondition_FieldSplit);CHKERRQ(ierr); 1446 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetSchurFactType_C","PCFieldSplitSetSchurFactType_FieldSplit",PCFieldSplitSetSchurFactType_FieldSplit);CHKERRQ(ierr); 1447 1448 } else { 1449 pc->ops->apply = PCApply_FieldSplit; 1450 pc->ops->view = PCView_FieldSplit; 1451 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit",PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 1452 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSchurPrecondition_C","",0);CHKERRQ(ierr); 1453 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetSchurFactType_C","",0);CHKERRQ(ierr); 1454 } 1455 PetscFunctionReturn(0); 1456 } 1457 EXTERN_C_END 1458 1459 EXTERN_C_BEGIN 1460 #undef __FUNCT__ 1461 #define __FUNCT__ "PCFieldSplitSetBlockSize_FieldSplit" 1462 PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs) 1463 { 1464 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1465 1466 PetscFunctionBegin; 1467 if (bs < 1) SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs); 1468 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); 1469 jac->bs = bs; 1470 PetscFunctionReturn(0); 1471 } 1472 EXTERN_C_END 1473 1474 #undef __FUNCT__ 1475 #define __FUNCT__ "PCFieldSplitSetType" 1476 /*@ 1477 PCFieldSplitSetType - Sets the type of fieldsplit preconditioner. 1478 1479 Collective on PC 1480 1481 Input Parameter: 1482 . pc - the preconditioner context 1483 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 1484 1485 Options Database Key: 1486 . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type 1487 1488 Level: Developer 1489 1490 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 1491 1492 .seealso: PCCompositeSetType() 1493 1494 @*/ 1495 PetscErrorCode PCFieldSplitSetType(PC pc,PCCompositeType type) 1496 { 1497 PetscErrorCode ierr; 1498 1499 PetscFunctionBegin; 1500 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1501 ierr = PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));CHKERRQ(ierr); 1502 PetscFunctionReturn(0); 1503 } 1504 1505 /* -------------------------------------------------------------------------------------*/ 1506 /*MC 1507 PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual 1508 fields or groups of fields. See the users manual section "Solving Block Matrices" for more details. 1509 1510 To set options on the solvers for each block append -fieldsplit_ to all the PC 1511 options database keys. For example, -fieldsplit_pc_type ilu -fieldsplit_pc_factor_levels 1 1512 1513 To set the options on the solvers separate for each block call PCFieldSplitGetSubKSP() 1514 and set the options directly on the resulting KSP object 1515 1516 Level: intermediate 1517 1518 Options Database Keys: 1519 + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split 1520 . -pc_fieldsplit_default - automatically add any fields to additional splits that have not 1521 been supplied explicitly by -pc_fieldsplit_%d_fields 1522 . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields) 1523 . -pc_fieldsplit_type <additive,multiplicative,schur,symmetric_multiplicative> 1524 . -pc_fieldsplit_schur_precondition <true,false> default is true 1525 . -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero or negative diagonal and uses Schur complement with no preconditioner as the solver 1526 1527 - Options prefix for inner solvers when using Schur complement preconditioner are -fieldsplit_0_ and -fieldsplit_1_ 1528 for all other solvers they are -fieldsplit_%d_ for the dth field, use -fieldsplit_ for all fields 1529 1530 Notes: 1531 Use PCFieldSplitSetFields() to set fields defined by "strided" entries and PCFieldSplitSetIS() 1532 to define a field by an arbitrary collection of entries. 1533 1534 If no fields are set the default is used. The fields are defined by entries strided by bs, 1535 beginning at 0 then 1, etc to bs-1. The block size can be set with PCFieldSplitSetBlockSize(), 1536 if this is not called the block size defaults to the blocksize of the second matrix passed 1537 to KSPSetOperators()/PCSetOperators(). 1538 1539 $ For the Schur complement preconditioner if J = ( A00 A01 ) 1540 $ ( A10 A11 ) 1541 $ the preconditioner using full factorization is 1542 $ ( I -A10 ksp(A00) ) ( inv(A00) 0 ) ( I 0 ) 1543 $ ( 0 I ) ( 0 ksp(S) ) ( -A10 ksp(A00) I ) 1544 where the action of inv(A00) is applied using the KSP solver with prefix -fieldsplit_0_. The action of 1545 ksp(S) is computed using the KSP solver with prefix -fieldsplit_splitname_ (where splitname was given 1546 in providing the SECOND split or 1 if not give). For PCFieldSplitGetKSP() when field number is 0, 1547 it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default 1548 A11 is used to construct a preconditioner for S, use PCFieldSplitSchurPrecondition() to turn on or off this 1549 option. You can use the preconditioner PCLSC to precondition the Schur complement with -fieldsplit_1_pc_type lsc. The 1550 factorization type is set using -pc_fieldsplit_schur_fact_type <diag, lower, upper, full>. The full is shown above, 1551 diag gives 1552 $ ( inv(A00) 0 ) 1553 $ ( 0 -ksp(S) ) 1554 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 1555 $ ( A00 0 ) 1556 $ ( A10 S ) 1557 where the inverses of A00 and S are applied using KSPs. The upper factorization is the inverse of 1558 $ ( A00 A01 ) 1559 $ ( 0 S ) 1560 where again the inverses of A00 and S are applied using KSPs. 1561 1562 If only one set of indices (one IS) is provided with PCFieldSplitSetIS() then the complement of that IS 1563 is used automatically for a second block. 1564 1565 The fieldsplit preconditioner cannot currently be used with the BAIJ or SBAIJ data formats if the blocksize is larger than 1. 1566 Generally it should be used with the AIJ format. 1567 1568 The forms of these preconditioners are closely related if not identical to forms derived as "Distributive Iterations", see, 1569 for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling. Note that one can also use PCFIELDSPLIT 1570 inside a smoother resulting in "Distributive Smoothers". 1571 1572 Concepts: physics based preconditioners, block preconditioners 1573 1574 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners, PCLSC, 1575 PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSchurPrecondition() 1576 M*/ 1577 1578 EXTERN_C_BEGIN 1579 #undef __FUNCT__ 1580 #define __FUNCT__ "PCCreate_FieldSplit" 1581 PetscErrorCode PCCreate_FieldSplit(PC pc) 1582 { 1583 PetscErrorCode ierr; 1584 PC_FieldSplit *jac; 1585 1586 PetscFunctionBegin; 1587 ierr = PetscNewLog(pc,PC_FieldSplit,&jac);CHKERRQ(ierr); 1588 jac->bs = -1; 1589 jac->nsplits = 0; 1590 jac->type = PC_COMPOSITE_MULTIPLICATIVE; 1591 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */ 1592 jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL; 1593 1594 pc->data = (void*)jac; 1595 1596 pc->ops->apply = PCApply_FieldSplit; 1597 pc->ops->applytranspose = PCApplyTranspose_FieldSplit; 1598 pc->ops->setup = PCSetUp_FieldSplit; 1599 pc->ops->reset = PCReset_FieldSplit; 1600 pc->ops->destroy = PCDestroy_FieldSplit; 1601 pc->ops->setfromoptions = PCSetFromOptions_FieldSplit; 1602 pc->ops->view = PCView_FieldSplit; 1603 pc->ops->applyrichardson = 0; 1604 1605 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit", 1606 PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 1607 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetFields_C","PCFieldSplitSetFields_FieldSplit", 1608 PCFieldSplitSetFields_FieldSplit);CHKERRQ(ierr); 1609 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetIS_C","PCFieldSplitSetIS_FieldSplit", 1610 PCFieldSplitSetIS_FieldSplit);CHKERRQ(ierr); 1611 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetType_C","PCFieldSplitSetType_FieldSplit", 1612 PCFieldSplitSetType_FieldSplit);CHKERRQ(ierr); 1613 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetBlockSize_C","PCFieldSplitSetBlockSize_FieldSplit", 1614 PCFieldSplitSetBlockSize_FieldSplit);CHKERRQ(ierr); 1615 PetscFunctionReturn(0); 1616 } 1617 EXTERN_C_END 1618 1619 1620 1621