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