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