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