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