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