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