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