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