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