1 2 #include <private/pcimpl.h> /*I "petscpc.h" I*/ 3 #include <petscdmcomposite.h> /*I "petscdmcomposite.h" I*/ 4 5 const char *const PCFieldSplitSchurPreTypes[] = {"SELF","DIAG","USER","PCFieldSplitSchurPreType","PC_FIELDSPLIT_SCHUR_PRE_",0}; 6 const char *const PCFieldSplitSchurFactorizationTypes[] = {"DIAG","LOWER","UPPER","FULL","PCFieldSplitSchurFactorizationType","PC_FIELDSPLIT_SCHUR_FACTORIZATION_",0}; 7 8 typedef 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 PCFieldSplitSchurFactorizationType 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,PCFieldSplitSchurFactorizationTypes[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_FACTORIZATION_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_FACTORIZATION_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_FACTORIZATION_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_FACTORIZATION_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,"PCFieldSplitSetSchurFactorizationType_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 = PetscOptionsEnum("-pc_fieldsplit_schur_factorization_type","Which off-diagonal parts of the block factorization to use","PCFieldSplitSetSchurFactorizationType",PCFieldSplitSchurFactorizationTypes,(PetscEnum)jac->schurfactorization,(PetscEnum*)&jac->schurfactorization,PETSC_NULL);CHKERRQ(ierr); 915 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); 916 } 917 ierr = PetscOptionsTail();CHKERRQ(ierr); 918 PetscFunctionReturn(0); 919 } 920 921 /*------------------------------------------------------------------------------------*/ 922 923 EXTERN_C_BEGIN 924 #undef __FUNCT__ 925 #define __FUNCT__ "PCFieldSplitSetFields_FieldSplit" 926 PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 927 { 928 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 929 PetscErrorCode ierr; 930 PC_FieldSplitLink ilink,next = jac->head; 931 char prefix[128]; 932 PetscInt i; 933 934 PetscFunctionBegin; 935 if (jac->splitdefined) { 936 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 937 PetscFunctionReturn(0); 938 } 939 for (i=0; i<n; i++) { 940 if (fields[i] >= jac->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs); 941 if (fields[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]); 942 } 943 ierr = PetscNew(struct _PC_FieldSplitLink,&ilink);CHKERRQ(ierr); 944 if (splitname) { 945 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 946 } else { 947 ierr = PetscMalloc(3*sizeof(char),&ilink->splitname);CHKERRQ(ierr); 948 ierr = PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);CHKERRQ(ierr); 949 } 950 ierr = PetscMalloc(n*sizeof(PetscInt),&ilink->fields);CHKERRQ(ierr); 951 ierr = PetscMemcpy(ilink->fields,fields,n*sizeof(PetscInt));CHKERRQ(ierr); 952 ierr = PetscMalloc(n*sizeof(PetscInt),&ilink->fields_col);CHKERRQ(ierr); 953 ierr = PetscMemcpy(ilink->fields_col,fields_col,n*sizeof(PetscInt));CHKERRQ(ierr); 954 ilink->nfields = n; 955 ilink->next = PETSC_NULL; 956 ierr = KSPCreate(((PetscObject)pc)->comm,&ilink->ksp);CHKERRQ(ierr); 957 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 958 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 959 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 960 961 ierr = PetscSNPrintf(prefix,sizeof prefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",ilink->splitname);CHKERRQ(ierr); 962 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 963 964 if (!next) { 965 jac->head = ilink; 966 ilink->previous = PETSC_NULL; 967 } else { 968 while (next->next) { 969 next = next->next; 970 } 971 next->next = ilink; 972 ilink->previous = next; 973 } 974 jac->nsplits++; 975 PetscFunctionReturn(0); 976 } 977 EXTERN_C_END 978 979 EXTERN_C_BEGIN 980 #undef __FUNCT__ 981 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit_Schur" 982 PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp) 983 { 984 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 985 PetscErrorCode ierr; 986 987 PetscFunctionBegin; 988 ierr = PetscMalloc(jac->nsplits*sizeof(KSP),subksp);CHKERRQ(ierr); 989 ierr = MatSchurComplementGetKSP(jac->schur,*subksp);CHKERRQ(ierr); 990 (*subksp)[1] = jac->kspschur; 991 if (n) *n = jac->nsplits; 992 PetscFunctionReturn(0); 993 } 994 EXTERN_C_END 995 996 EXTERN_C_BEGIN 997 #undef __FUNCT__ 998 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit" 999 PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp) 1000 { 1001 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1002 PetscErrorCode ierr; 1003 PetscInt cnt = 0; 1004 PC_FieldSplitLink ilink = jac->head; 1005 1006 PetscFunctionBegin; 1007 ierr = PetscMalloc(jac->nsplits*sizeof(KSP),subksp);CHKERRQ(ierr); 1008 while (ilink) { 1009 (*subksp)[cnt++] = ilink->ksp; 1010 ilink = ilink->next; 1011 } 1012 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); 1013 if (n) *n = jac->nsplits; 1014 PetscFunctionReturn(0); 1015 } 1016 EXTERN_C_END 1017 1018 EXTERN_C_BEGIN 1019 #undef __FUNCT__ 1020 #define __FUNCT__ "PCFieldSplitSetIS_FieldSplit" 1021 PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is) 1022 { 1023 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1024 PetscErrorCode ierr; 1025 PC_FieldSplitLink ilink, next = jac->head; 1026 char prefix[128]; 1027 1028 PetscFunctionBegin; 1029 if (jac->splitdefined) { 1030 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 1031 PetscFunctionReturn(0); 1032 } 1033 ierr = PetscNew(struct _PC_FieldSplitLink,&ilink);CHKERRQ(ierr); 1034 if (splitname) { 1035 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 1036 } else { 1037 ierr = PetscMalloc(8*sizeof(char),&ilink->splitname);CHKERRQ(ierr); 1038 ierr = PetscSNPrintf(ilink->splitname,7,"%D",jac->nsplits);CHKERRQ(ierr); 1039 } 1040 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 1041 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 1042 ilink->is = is; 1043 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 1044 ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr); 1045 ilink->is_col = is; 1046 ilink->next = PETSC_NULL; 1047 ierr = KSPCreate(((PetscObject)pc)->comm,&ilink->ksp);CHKERRQ(ierr); 1048 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 1049 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 1050 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 1051 1052 ierr = PetscSNPrintf(prefix,sizeof prefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",ilink->splitname);CHKERRQ(ierr); 1053 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 1054 1055 if (!next) { 1056 jac->head = ilink; 1057 ilink->previous = PETSC_NULL; 1058 } else { 1059 while (next->next) { 1060 next = next->next; 1061 } 1062 next->next = ilink; 1063 ilink->previous = next; 1064 } 1065 jac->nsplits++; 1066 1067 PetscFunctionReturn(0); 1068 } 1069 EXTERN_C_END 1070 1071 #undef __FUNCT__ 1072 #define __FUNCT__ "PCFieldSplitSetFields" 1073 /*@ 1074 PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner 1075 1076 Logically Collective on PC 1077 1078 Input Parameters: 1079 + pc - the preconditioner context 1080 . splitname - name of this split, if PETSC_NULL the number of the split is used 1081 . n - the number of fields in this split 1082 - fields - the fields in this split 1083 1084 Level: intermediate 1085 1086 Notes: Use PCFieldSplitSetIS() to set a completely general set of indices as a field. 1087 1088 The PCFieldSplitSetFields() is for defining fields as a strided blocks. For example, if the block 1089 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 1090 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x.... 1091 where the numbered entries indicate what is in the field. 1092 1093 This function is called once per split (it creates a new split each time). Solve options 1094 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1095 1096 Developer Note: This routine does not actually create the IS representing the split, that is delayed 1097 until PCSetUp_FieldSplit(), because information about the vector/matrix layouts may not be 1098 available when this routine is called. 1099 1100 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize(), PCFieldSplitSetIS() 1101 1102 @*/ 1103 PetscErrorCode PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 1104 { 1105 PetscErrorCode ierr; 1106 1107 PetscFunctionBegin; 1108 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1109 PetscValidCharPointer(splitname,2); 1110 if (n < 1) SETERRQ2(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Provided number of fields %D in split \"%s\" not positive",n,splitname); 1111 PetscValidIntPointer(fields,3); 1112 ierr = PetscTryMethod(pc,"PCFieldSplitSetFields_C",(PC,const char[],PetscInt,const PetscInt *,const PetscInt *),(pc,splitname,n,fields,fields_col));CHKERRQ(ierr); 1113 PetscFunctionReturn(0); 1114 } 1115 1116 #undef __FUNCT__ 1117 #define __FUNCT__ "PCFieldSplitSetIS" 1118 /*@ 1119 PCFieldSplitSetIS - Sets the exact elements for field 1120 1121 Logically Collective on PC 1122 1123 Input Parameters: 1124 + pc - the preconditioner context 1125 . splitname - name of this split, if PETSC_NULL the number of the split is used 1126 - is - the index set that defines the vector elements in this field 1127 1128 1129 Notes: 1130 Use PCFieldSplitSetFields(), for fields defined by strided types. 1131 1132 This function is called once per split (it creates a new split each time). Solve options 1133 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1134 1135 Level: intermediate 1136 1137 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize() 1138 1139 @*/ 1140 PetscErrorCode PCFieldSplitSetIS(PC pc,const char splitname[],IS is) 1141 { 1142 PetscErrorCode ierr; 1143 1144 PetscFunctionBegin; 1145 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1146 if (splitname) PetscValidCharPointer(splitname,2); 1147 PetscValidHeaderSpecific(is,IS_CLASSID,3); 1148 ierr = PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));CHKERRQ(ierr); 1149 PetscFunctionReturn(0); 1150 } 1151 1152 #undef __FUNCT__ 1153 #define __FUNCT__ "PCFieldSplitGetIS" 1154 /*@ 1155 PCFieldSplitGetIS - Retrieves the elements for a field as an IS 1156 1157 Logically Collective on PC 1158 1159 Input Parameters: 1160 + pc - the preconditioner context 1161 - splitname - name of this split 1162 1163 Output Parameter: 1164 - is - the index set that defines the vector elements in this field, or PETSC_NULL if the field is not found 1165 1166 Level: intermediate 1167 1168 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetIS() 1169 1170 @*/ 1171 PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is) 1172 { 1173 PetscErrorCode ierr; 1174 1175 PetscFunctionBegin; 1176 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1177 PetscValidCharPointer(splitname,2); 1178 PetscValidPointer(is,3); 1179 { 1180 PC_FieldSplit *jac = (PC_FieldSplit *) pc->data; 1181 PC_FieldSplitLink ilink = jac->head; 1182 PetscBool found; 1183 1184 *is = PETSC_NULL; 1185 while(ilink) { 1186 ierr = PetscStrcmp(ilink->splitname, splitname, &found);CHKERRQ(ierr); 1187 if (found) { 1188 *is = ilink->is; 1189 break; 1190 } 1191 ilink = ilink->next; 1192 } 1193 } 1194 PetscFunctionReturn(0); 1195 } 1196 1197 #undef __FUNCT__ 1198 #define __FUNCT__ "PCFieldSplitSetBlockSize" 1199 /*@ 1200 PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the 1201 fieldsplit preconditioner. If not set the matrix block size is used. 1202 1203 Logically Collective on PC 1204 1205 Input Parameters: 1206 + pc - the preconditioner context 1207 - bs - the block size 1208 1209 Level: intermediate 1210 1211 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields() 1212 1213 @*/ 1214 PetscErrorCode PCFieldSplitSetBlockSize(PC pc,PetscInt bs) 1215 { 1216 PetscErrorCode ierr; 1217 1218 PetscFunctionBegin; 1219 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1220 PetscValidLogicalCollectiveInt(pc,bs,2); 1221 ierr = PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));CHKERRQ(ierr); 1222 PetscFunctionReturn(0); 1223 } 1224 1225 #undef __FUNCT__ 1226 #define __FUNCT__ "PCFieldSplitGetSubKSP" 1227 /*@C 1228 PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits 1229 1230 Collective on KSP 1231 1232 Input Parameter: 1233 . pc - the preconditioner context 1234 1235 Output Parameters: 1236 + n - the number of splits 1237 - pc - the array of KSP contexts 1238 1239 Note: 1240 After PCFieldSplitGetSubKSP() the array of KSPs IS to be freed by the user 1241 (not the KSP just the array that contains them). 1242 1243 You must call KSPSetUp() before calling PCFieldSplitGetSubKSP(). 1244 1245 Level: advanced 1246 1247 .seealso: PCFIELDSPLIT 1248 @*/ 1249 PetscErrorCode PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[]) 1250 { 1251 PetscErrorCode ierr; 1252 1253 PetscFunctionBegin; 1254 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1255 if (n) PetscValidIntPointer(n,2); 1256 ierr = PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));CHKERRQ(ierr); 1257 PetscFunctionReturn(0); 1258 } 1259 1260 #undef __FUNCT__ 1261 #define __FUNCT__ "PCFieldSplitSchurPrecondition" 1262 /*@ 1263 PCFieldSplitSchurPrecondition - Indicates if the Schur complement is preconditioned by a preconditioner constructed by the 1264 A11 matrix. Otherwise no preconditioner is used. 1265 1266 Collective on PC 1267 1268 Input Parameters: 1269 + pc - the preconditioner context 1270 . ptype - which matrix to use for preconditioning the Schur complement, PC_FIELDSPLIT_SCHUR_PRE_DIAG (diag) is default 1271 - userpre - matrix to use for preconditioning, or PETSC_NULL 1272 1273 Options Database: 1274 . -pc_fieldsplit_schur_precondition <self,user,diag> default is diag 1275 1276 Notes: 1277 $ If ptype is 1278 $ user then the preconditioner for the Schur complement is generated by the provided matrix (pre argument 1279 $ to this function). 1280 $ diag then the preconditioner for the Schur complement is generated by the block diagonal part of the original 1281 $ matrix associated with the Schur complement (i.e. A11) 1282 $ self the preconditioner for the Schur complement is generated from the Schur complement matrix itself: 1283 $ The only preconditioner that currently works directly with the Schur complement matrix object is the PCLSC 1284 $ preconditioner 1285 1286 When solving a saddle point problem, where the A11 block is identically zero, using diag as the ptype only makes sense 1287 with the additional option -fieldsplit_1_pc_type none. Usually for saddle point problems one would use a ptype of self and 1288 -fieldsplit_1_pc_type lsc which uses the least squares commutator compute a preconditioner for the Schur complement. 1289 1290 Developer Notes: This is a terrible name, gives no good indication of what the function does and should also have Set in 1291 the name since it sets a proceedure to use. 1292 1293 Level: intermediate 1294 1295 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCLSC 1296 1297 @*/ 1298 PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1299 { 1300 PetscErrorCode ierr; 1301 1302 PetscFunctionBegin; 1303 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1304 ierr = PetscTryMethod(pc,"PCFieldSplitSchurPrecondition_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));CHKERRQ(ierr); 1305 PetscFunctionReturn(0); 1306 } 1307 1308 EXTERN_C_BEGIN 1309 #undef __FUNCT__ 1310 #define __FUNCT__ "PCFieldSplitSchurPrecondition_FieldSplit" 1311 PetscErrorCode PCFieldSplitSchurPrecondition_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1312 { 1313 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1314 PetscErrorCode ierr; 1315 1316 PetscFunctionBegin; 1317 jac->schurpre = ptype; 1318 if (pre) { 1319 ierr = MatDestroy(&jac->schur_user);CHKERRQ(ierr); 1320 jac->schur_user = pre; 1321 ierr = PetscObjectReference((PetscObject)jac->schur_user);CHKERRQ(ierr); 1322 } 1323 PetscFunctionReturn(0); 1324 } 1325 EXTERN_C_END 1326 1327 #undef __FUNCT__ 1328 #define __FUNCT__ "PCFieldSplitSetSchurFactorizationType" 1329 /*@ 1330 PCFieldSplitSetSchurFactorizationType - sets which blocks of the approximate block factorization to retain 1331 1332 Collective on PC 1333 1334 Input Parameters: 1335 + pc - the preconditioner context 1336 - ftype - which blocks of factorization to retain, PC_FIELDSPLIT_SCHUR_FACTORIZATION_FULL is default 1337 1338 Options Database: 1339 . -pc_fieldsplit_schur_factorization_type <diag,lower,upper,full> default is full 1340 1341 1342 Level: intermediate 1343 1344 Notes: 1345 The FULL factorization is 1346 1347 $ (A B) = (1 0) (A 0) (1 Ainv*B) 1348 $ (C D) (C*Ainv 1) (0 S) (0 1 ) 1349 1350 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, 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). 1351 1352 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 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 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, the preconditioned operator has three distinct nonzero eigenvalues and minimal polynomial of degree at most 4, so KSPGMRES converges in at most 4 iterations. 1353 1354 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 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). 1355 1356 References: 1357 Murphy, Golub, and Wathen, A note on preconditioning indefinite linear systems, SIAM J. Sci. Comput., 21 (2000) pp. 1969-1972. 1358 1359 Ipsen, A note on preconditioning nonsymmetric matrices, SIAM J. Sci. Comput., 23 (2001), pp. 1050-1051. 1360 1361 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType 1362 @*/ 1363 PetscErrorCode PCFieldSplitSetSchurFactorizationType(PC pc,PCFieldSplitSchurFactorizationType ftype) 1364 { 1365 PetscErrorCode ierr; 1366 1367 PetscFunctionBegin; 1368 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1369 ierr = PetscTryMethod(pc,"PCFieldSplitSetSchurFactorizationType_C",(PC,PCFieldSplitSchurFactorizationType),(pc,ftype));CHKERRQ(ierr); 1370 PetscFunctionReturn(0); 1371 } 1372 1373 #undef __FUNCT__ 1374 #define __FUNCT__ "PCFieldSplitSetSchurFactorizationType_FieldSplit" 1375 PETSC_EXTERN_C PetscErrorCode PCFieldSplitSetSchurFactorizationType_FieldSplit(PC pc,PCFieldSplitSchurFactorizationType ftype) 1376 { 1377 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1378 1379 PetscFunctionBegin; 1380 jac->schurfactorization = ftype; 1381 PetscFunctionReturn(0); 1382 } 1383 1384 #undef __FUNCT__ 1385 #define __FUNCT__ "PCFieldSplitGetSchurBlocks" 1386 /*@C 1387 PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement 1388 1389 Collective on KSP 1390 1391 Input Parameter: 1392 . pc - the preconditioner context 1393 1394 Output Parameters: 1395 + A00 - the (0,0) block 1396 . A01 - the (0,1) block 1397 . A10 - the (1,0) block 1398 - A11 - the (1,1) block 1399 1400 Level: advanced 1401 1402 .seealso: PCFIELDSPLIT 1403 @*/ 1404 PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11) 1405 { 1406 PC_FieldSplit *jac = (PC_FieldSplit *) pc->data; 1407 1408 PetscFunctionBegin; 1409 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1410 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach."); 1411 if (A00) *A00 = jac->pmat[0]; 1412 if (A01) *A01 = jac->B; 1413 if (A10) *A10 = jac->C; 1414 if (A11) *A11 = jac->pmat[1]; 1415 PetscFunctionReturn(0); 1416 } 1417 1418 EXTERN_C_BEGIN 1419 #undef __FUNCT__ 1420 #define __FUNCT__ "PCFieldSplitSetType_FieldSplit" 1421 PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type) 1422 { 1423 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1424 PetscErrorCode ierr; 1425 1426 PetscFunctionBegin; 1427 jac->type = type; 1428 if (type == PC_COMPOSITE_SCHUR) { 1429 pc->ops->apply = PCApply_FieldSplit_Schur; 1430 pc->ops->view = PCView_FieldSplit_Schur; 1431 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit_Schur",PCFieldSplitGetSubKSP_FieldSplit_Schur);CHKERRQ(ierr); 1432 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSchurPrecondition_C","PCFieldSplitSchurPrecondition_FieldSplit",PCFieldSplitSchurPrecondition_FieldSplit);CHKERRQ(ierr); 1433 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetSchurFactorizationType_C","PCFieldSplitSetSchurFactorizationType_FieldSplit",PCFieldSplitSetSchurFactorizationType_FieldSplit);CHKERRQ(ierr); 1434 1435 } else { 1436 pc->ops->apply = PCApply_FieldSplit; 1437 pc->ops->view = PCView_FieldSplit; 1438 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit",PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 1439 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSchurPrecondition_C","",0);CHKERRQ(ierr); 1440 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetSchurFactorizationType_C","",0);CHKERRQ(ierr); 1441 } 1442 PetscFunctionReturn(0); 1443 } 1444 EXTERN_C_END 1445 1446 EXTERN_C_BEGIN 1447 #undef __FUNCT__ 1448 #define __FUNCT__ "PCFieldSplitSetBlockSize_FieldSplit" 1449 PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs) 1450 { 1451 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1452 1453 PetscFunctionBegin; 1454 if (bs < 1) SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs); 1455 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); 1456 jac->bs = bs; 1457 PetscFunctionReturn(0); 1458 } 1459 EXTERN_C_END 1460 1461 #undef __FUNCT__ 1462 #define __FUNCT__ "PCFieldSplitSetType" 1463 /*@ 1464 PCFieldSplitSetType - Sets the type of fieldsplit preconditioner. 1465 1466 Collective on PC 1467 1468 Input Parameter: 1469 . pc - the preconditioner context 1470 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 1471 1472 Options Database Key: 1473 . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type 1474 1475 Level: Developer 1476 1477 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 1478 1479 .seealso: PCCompositeSetType() 1480 1481 @*/ 1482 PetscErrorCode PCFieldSplitSetType(PC pc,PCCompositeType type) 1483 { 1484 PetscErrorCode ierr; 1485 1486 PetscFunctionBegin; 1487 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1488 ierr = PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));CHKERRQ(ierr); 1489 PetscFunctionReturn(0); 1490 } 1491 1492 /* -------------------------------------------------------------------------------------*/ 1493 /*MC 1494 PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual 1495 fields or groups of fields. See the users manual section "Solving Block Matrices" for more details. 1496 1497 To set options on the solvers for each block append -fieldsplit_ to all the PC 1498 options database keys. For example, -fieldsplit_pc_type ilu -fieldsplit_pc_factor_levels 1 1499 1500 To set the options on the solvers separate for each block call PCFieldSplitGetSubKSP() 1501 and set the options directly on the resulting KSP object 1502 1503 Level: intermediate 1504 1505 Options Database Keys: 1506 + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split 1507 . -pc_fieldsplit_default - automatically add any fields to additional splits that have not 1508 been supplied explicitly by -pc_fieldsplit_%d_fields 1509 . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields) 1510 . -pc_fieldsplit_type <additive,multiplicative,schur,symmetric_multiplicative> 1511 . -pc_fieldsplit_schur_precondition <true,false> default is true 1512 . -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero or negative diagonal and uses Schur complement with no preconditioner as the solver 1513 1514 - Options prefix for inner solvers when using Schur complement preconditioner are -fieldsplit_0_ and -fieldsplit_1_ 1515 for all other solvers they are -fieldsplit_%d_ for the dth field, use -fieldsplit_ for all fields 1516 1517 Notes: 1518 Use PCFieldSplitSetFields() to set fields defined by "strided" entries and PCFieldSplitSetIS() 1519 to define a field by an arbitrary collection of entries. 1520 1521 If no fields are set the default is used. The fields are defined by entries strided by bs, 1522 beginning at 0 then 1, etc to bs-1. The block size can be set with PCFieldSplitSetBlockSize(), 1523 if this is not called the block size defaults to the blocksize of the second matrix passed 1524 to KSPSetOperators()/PCSetOperators(). 1525 1526 $ For the Schur complement preconditioner if J = ( A00 A01 ) 1527 $ ( A10 A11 ) 1528 $ the preconditioner using full factorization is 1529 $ ( I -A10 ksp(A00) ) ( inv(A00) 0 ) ( I 0 ) 1530 $ ( 0 I ) ( 0 ksp(S) ) ( -A10 ksp(A00) I ) 1531 where the action of inv(A00) is applied using the KSP solver with prefix -fieldsplit_0_. The action of 1532 ksp(S) is computed using the KSP solver with prefix -fieldsplit_splitname_ (where splitname was given 1533 in providing the SECOND split or 1 if not give). For PCFieldSplitGetKSP() when field number is 0, 1534 it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default 1535 A11 is used to construct a preconditioner for S, use PCFieldSplitSchurPrecondition() to turn on or off this 1536 option. You can use the preconditioner PCLSC to precondition the Schur complement with -fieldsplit_1_pc_type lsc. The 1537 factorization type is set using -pc_fieldsplit_schur_factorization_type <diag, lower, upper, full>. The full is shown above, 1538 diag gives 1539 $ ( inv(A00) 0 ) 1540 $ ( 0 -ksp(S) ) 1541 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 1542 $ ( A00 0 ) 1543 $ ( A10 S ) 1544 where the inverses of A00 and S are applied using KSPs. The upper factorization is the inverse of 1545 $ ( A00 A01 ) 1546 $ ( 0 S ) 1547 where again the inverses of A00 and S are applied using KSPs. 1548 1549 If only one set of indices (one IS) is provided with PCFieldSplitSetIS() then the complement of that IS 1550 is used automatically for a second block. 1551 1552 The fieldsplit preconditioner cannot currently be used with the BAIJ or SBAIJ data formats if the blocksize is larger than 1. 1553 Generally it should be used with the AIJ format. 1554 1555 The forms of these preconditioners are closely related if not identical to forms derived as "Distributive Iterations", see, 1556 for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling. Note that one can also use PCFIELDSPLIT 1557 inside a smoother resulting in "Distributive Smoothers". 1558 1559 Concepts: physics based preconditioners, block preconditioners 1560 1561 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners, PCLSC, 1562 PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSchurPrecondition() 1563 M*/ 1564 1565 EXTERN_C_BEGIN 1566 #undef __FUNCT__ 1567 #define __FUNCT__ "PCCreate_FieldSplit" 1568 PetscErrorCode PCCreate_FieldSplit(PC pc) 1569 { 1570 PetscErrorCode ierr; 1571 PC_FieldSplit *jac; 1572 1573 PetscFunctionBegin; 1574 ierr = PetscNewLog(pc,PC_FieldSplit,&jac);CHKERRQ(ierr); 1575 jac->bs = -1; 1576 jac->nsplits = 0; 1577 jac->type = PC_COMPOSITE_MULTIPLICATIVE; 1578 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */ 1579 jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACTORIZATION_FULL; 1580 1581 pc->data = (void*)jac; 1582 1583 pc->ops->apply = PCApply_FieldSplit; 1584 pc->ops->applytranspose = PCApplyTranspose_FieldSplit; 1585 pc->ops->setup = PCSetUp_FieldSplit; 1586 pc->ops->reset = PCReset_FieldSplit; 1587 pc->ops->destroy = PCDestroy_FieldSplit; 1588 pc->ops->setfromoptions = PCSetFromOptions_FieldSplit; 1589 pc->ops->view = PCView_FieldSplit; 1590 pc->ops->applyrichardson = 0; 1591 1592 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit", 1593 PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 1594 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetFields_C","PCFieldSplitSetFields_FieldSplit", 1595 PCFieldSplitSetFields_FieldSplit);CHKERRQ(ierr); 1596 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetIS_C","PCFieldSplitSetIS_FieldSplit", 1597 PCFieldSplitSetIS_FieldSplit);CHKERRQ(ierr); 1598 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetType_C","PCFieldSplitSetType_FieldSplit", 1599 PCFieldSplitSetType_FieldSplit);CHKERRQ(ierr); 1600 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetBlockSize_C","PCFieldSplitSetBlockSize_FieldSplit", 1601 PCFieldSplitSetBlockSize_FieldSplit);CHKERRQ(ierr); 1602 PetscFunctionReturn(0); 1603 } 1604 EXTERN_C_END 1605 1606 1607 1608