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