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 496 ierr = KSPCreate(((PetscObject)pc)->comm,&jac->kspschur);CHKERRQ(ierr); 497 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->kspschur);CHKERRQ(ierr); 498 ierr = PetscObjectIncrementTabLevel((PetscObject)jac->kspschur,(PetscObject)pc,1);CHKERRQ(ierr); 499 ierr = KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac),DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 500 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) { 501 PC pc; 502 ierr = KSPGetPC(jac->kspschur,&pc);CHKERRQ(ierr); 503 ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr); 504 /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */ 505 } 506 ierr = PetscSNPrintf(schurprefix,sizeof schurprefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",ilink->splitname);CHKERRQ(ierr); 507 ierr = KSPSetOptionsPrefix(jac->kspschur,schurprefix);CHKERRQ(ierr); 508 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 509 /* need to call this every time, since the jac->kspschur is freshly created, otherwise its options never get set */ 510 ierr = KSPSetFromOptions(jac->kspschur);CHKERRQ(ierr); 511 512 ierr = PetscMalloc2(2,Vec,&jac->x,2,Vec,&jac->y);CHKERRQ(ierr); 513 ierr = MatGetVecs(jac->pmat[0],&jac->x[0],&jac->y[0]);CHKERRQ(ierr); 514 ierr = MatGetVecs(jac->pmat[1],&jac->x[1],&jac->y[1]);CHKERRQ(ierr); 515 ilink = jac->head; 516 ilink->x = jac->x[0]; ilink->y = jac->y[0]; 517 ilink = ilink->next; 518 ilink->x = jac->x[1]; ilink->y = jac->y[1]; 519 } 520 } else { 521 /* set up the individual PCs */ 522 i = 0; 523 ilink = jac->head; 524 while (ilink) { 525 ierr = KSPSetOperators(ilink->ksp,jac->mat[i],jac->pmat[i],flag);CHKERRQ(ierr); 526 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 527 if (!jac->suboptionsset) {ierr = KSPSetFromOptions(ilink->ksp);CHKERRQ(ierr);} 528 ierr = KSPSetUp(ilink->ksp);CHKERRQ(ierr); 529 i++; 530 ilink = ilink->next; 531 } 532 533 /* create work vectors for each split */ 534 if (!jac->x) { 535 ierr = PetscMalloc2(nsplit,Vec,&jac->x,nsplit,Vec,&jac->y);CHKERRQ(ierr); 536 ilink = jac->head; 537 for (i=0; i<nsplit; i++) { 538 Vec *vl,*vr; 539 540 ierr = KSPGetVecs(ilink->ksp,1,&vr,1,&vl);CHKERRQ(ierr); 541 ilink->x = *vr; 542 ilink->y = *vl; 543 ierr = PetscFree(vr);CHKERRQ(ierr); 544 ierr = PetscFree(vl);CHKERRQ(ierr); 545 jac->x[i] = ilink->x; 546 jac->y[i] = ilink->y; 547 ilink = ilink->next; 548 } 549 } 550 } 551 552 553 if (!jac->head->sctx) { 554 Vec xtmp; 555 556 /* compute scatter contexts needed by multiplicative versions and non-default splits */ 557 558 ilink = jac->head; 559 ierr = MatGetVecs(pc->pmat,&xtmp,PETSC_NULL);CHKERRQ(ierr); 560 for (i=0; i<nsplit; i++) { 561 ierr = VecScatterCreate(xtmp,ilink->is,jac->x[i],PETSC_NULL,&ilink->sctx);CHKERRQ(ierr); 562 ilink = ilink->next; 563 } 564 ierr = VecDestroy(&xtmp);CHKERRQ(ierr); 565 } 566 jac->suboptionsset = PETSC_TRUE; 567 PetscFunctionReturn(0); 568 } 569 570 #define FieldSplitSplitSolveAdd(ilink,xx,yy) \ 571 (VecScatterBegin(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 572 VecScatterEnd(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 573 KSPSolve(ilink->ksp,ilink->x,ilink->y) || \ 574 VecScatterBegin(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE) || \ 575 VecScatterEnd(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE)) 576 577 #undef __FUNCT__ 578 #define __FUNCT__ "PCApply_FieldSplit_Schur" 579 static PetscErrorCode PCApply_FieldSplit_Schur(PC pc,Vec x,Vec y) 580 { 581 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 582 PetscErrorCode ierr; 583 KSP ksp; 584 PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next; 585 586 PetscFunctionBegin; 587 ierr = MatSchurComplementGetKSP(jac->schur,&ksp);CHKERRQ(ierr); 588 589 switch (jac->schurfactorization) { 590 case PC_FIELDSPLIT_SCHUR_FACTORIZATION_DIAG: 591 /* [A00 0; 0 -S], positive definite, suitable for MINRES */ 592 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 593 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 594 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 595 ierr = KSPSolve(ksp,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 596 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 597 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 598 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 599 ierr = VecScale(ilinkD->y,-1.);CHKERRQ(ierr); 600 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 601 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 602 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 603 break; 604 case PC_FIELDSPLIT_SCHUR_FACTORIZATION_LOWER: 605 /* [A00 0; A10 S], suitable for left preconditioning */ 606 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 607 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 608 ierr = KSPSolve(ksp,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 609 ierr = MatMult(jac->C,ilinkA->y,ilinkD->x);CHKERRQ(ierr); 610 ierr = VecScale(ilinkD->x,-1.);CHKERRQ(ierr); 611 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 612 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 613 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 614 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 615 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 616 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 617 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 618 break; 619 case PC_FIELDSPLIT_SCHUR_FACTORIZATION_UPPER: 620 /* [A00 A01; 0 S], suitable for right preconditioning */ 621 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 622 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 623 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 624 ierr = MatMult(jac->B,ilinkD->y,ilinkA->x);CHKERRQ(ierr); 625 ierr = VecScale(ilinkA->x,-1.);CHKERRQ(ierr); 626 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 627 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 628 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 629 ierr = KSPSolve(ksp,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 630 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 631 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 632 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 633 break; 634 case PC_FIELDSPLIT_SCHUR_FACTORIZATION_FULL: 635 /* [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 */ 636 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 637 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 638 ierr = KSPSolve(ksp,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 639 ierr = MatMult(jac->C,ilinkA->y,ilinkD->x);CHKERRQ(ierr); 640 ierr = VecScale(ilinkD->x,-1.0);CHKERRQ(ierr); 641 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 642 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 643 644 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 645 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 646 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 647 648 ierr = MatMult(jac->B,ilinkD->y,ilinkA->y);CHKERRQ(ierr); 649 ierr = VecAXPY(ilinkA->x,-1.0,ilinkA->y);CHKERRQ(ierr); 650 ierr = KSPSolve(ksp,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 651 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 652 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 653 } 654 PetscFunctionReturn(0); 655 } 656 657 #undef __FUNCT__ 658 #define __FUNCT__ "PCApply_FieldSplit" 659 static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y) 660 { 661 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 662 PetscErrorCode ierr; 663 PC_FieldSplitLink ilink = jac->head; 664 PetscInt cnt; 665 666 PetscFunctionBegin; 667 CHKMEMQ; 668 ierr = VecSetBlockSize(x,jac->bs);CHKERRQ(ierr); 669 ierr = VecSetBlockSize(y,jac->bs);CHKERRQ(ierr); 670 671 if (jac->type == PC_COMPOSITE_ADDITIVE) { 672 if (jac->defaultsplit) { 673 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 674 while (ilink) { 675 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 676 ilink = ilink->next; 677 } 678 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 679 } else { 680 ierr = VecSet(y,0.0);CHKERRQ(ierr); 681 while (ilink) { 682 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 683 ilink = ilink->next; 684 } 685 } 686 } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 687 if (!jac->w1) { 688 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 689 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 690 } 691 ierr = VecSet(y,0.0);CHKERRQ(ierr); 692 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 693 cnt = 1; 694 while (ilink->next) { 695 ilink = ilink->next; 696 /* compute the residual only over the part of the vector needed */ 697 ierr = MatMult(jac->Afield[cnt++],y,ilink->x);CHKERRQ(ierr); 698 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 699 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 700 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 701 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 702 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 703 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 704 } 705 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 706 cnt -= 2; 707 while (ilink->previous) { 708 ilink = ilink->previous; 709 /* compute the residual only over the part of the vector needed */ 710 ierr = MatMult(jac->Afield[cnt--],y,ilink->x);CHKERRQ(ierr); 711 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 712 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 713 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 714 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 715 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 716 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 717 } 718 } 719 } else SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Unsupported or unknown composition",(int) jac->type); 720 CHKMEMQ; 721 PetscFunctionReturn(0); 722 } 723 724 #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \ 725 (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 726 VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 727 KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) || \ 728 VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \ 729 VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE)) 730 731 #undef __FUNCT__ 732 #define __FUNCT__ "PCApplyTranspose_FieldSplit" 733 static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y) 734 { 735 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 736 PetscErrorCode ierr; 737 PC_FieldSplitLink ilink = jac->head; 738 739 PetscFunctionBegin; 740 CHKMEMQ; 741 ierr = VecSetBlockSize(x,jac->bs);CHKERRQ(ierr); 742 ierr = VecSetBlockSize(y,jac->bs);CHKERRQ(ierr); 743 744 if (jac->type == PC_COMPOSITE_ADDITIVE) { 745 if (jac->defaultsplit) { 746 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 747 while (ilink) { 748 ierr = KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 749 ilink = ilink->next; 750 } 751 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 752 } else { 753 ierr = VecSet(y,0.0);CHKERRQ(ierr); 754 while (ilink) { 755 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 756 ilink = ilink->next; 757 } 758 } 759 } else { 760 if (!jac->w1) { 761 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 762 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 763 } 764 ierr = VecSet(y,0.0);CHKERRQ(ierr); 765 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 766 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 767 while (ilink->next) { 768 ilink = ilink->next; 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 while (ilink->previous) { 774 ilink = ilink->previous; 775 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 776 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 777 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 778 } 779 } else { 780 while (ilink->next) { /* get to last entry in linked list */ 781 ilink = ilink->next; 782 } 783 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 784 while (ilink->previous) { 785 ilink = ilink->previous; 786 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 787 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 788 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 789 } 790 } 791 } 792 CHKMEMQ; 793 PetscFunctionReturn(0); 794 } 795 796 #undef __FUNCT__ 797 #define __FUNCT__ "PCReset_FieldSplit" 798 static PetscErrorCode PCReset_FieldSplit(PC pc) 799 { 800 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 801 PetscErrorCode ierr; 802 PC_FieldSplitLink ilink = jac->head,next; 803 804 PetscFunctionBegin; 805 while (ilink) { 806 ierr = KSPReset(ilink->ksp);CHKERRQ(ierr); 807 ierr = VecDestroy(&ilink->x);CHKERRQ(ierr); 808 ierr = VecDestroy(&ilink->y);CHKERRQ(ierr); 809 ierr = VecScatterDestroy(&ilink->sctx);CHKERRQ(ierr); 810 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 811 next = ilink->next; 812 ilink = next; 813 } 814 ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr); 815 if (jac->mat && jac->mat != jac->pmat) { 816 ierr = MatDestroyMatrices(jac->nsplits,&jac->mat);CHKERRQ(ierr); 817 } else if (jac->mat) { 818 jac->mat = PETSC_NULL; 819 } 820 if (jac->pmat) {ierr = MatDestroyMatrices(jac->nsplits,&jac->pmat);CHKERRQ(ierr);} 821 if (jac->Afield) {ierr = MatDestroyMatrices(jac->nsplits,&jac->Afield);CHKERRQ(ierr);} 822 ierr = VecDestroy(&jac->w1);CHKERRQ(ierr); 823 ierr = VecDestroy(&jac->w2);CHKERRQ(ierr); 824 ierr = MatDestroy(&jac->schur);CHKERRQ(ierr); 825 ierr = MatDestroy(&jac->schur_user);CHKERRQ(ierr); 826 ierr = KSPDestroy(&jac->kspschur);CHKERRQ(ierr); 827 ierr = MatDestroy(&jac->B);CHKERRQ(ierr); 828 ierr = MatDestroy(&jac->C);CHKERRQ(ierr); 829 jac->reset = PETSC_TRUE; 830 PetscFunctionReturn(0); 831 } 832 833 #undef __FUNCT__ 834 #define __FUNCT__ "PCDestroy_FieldSplit" 835 static PetscErrorCode PCDestroy_FieldSplit(PC pc) 836 { 837 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 838 PetscErrorCode ierr; 839 PC_FieldSplitLink ilink = jac->head,next; 840 841 PetscFunctionBegin; 842 ierr = PCReset_FieldSplit(pc);CHKERRQ(ierr); 843 while (ilink) { 844 ierr = KSPDestroy(&ilink->ksp);CHKERRQ(ierr); 845 next = ilink->next; 846 ierr = PetscFree(ilink->splitname);CHKERRQ(ierr); 847 ierr = PetscFree(ilink->fields);CHKERRQ(ierr); 848 ierr = PetscFree(ilink);CHKERRQ(ierr); 849 ilink = next; 850 } 851 ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr); 852 ierr = PetscFree(pc->data);CHKERRQ(ierr); 853 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","",PETSC_NULL);CHKERRQ(ierr); 854 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetFields_C","",PETSC_NULL);CHKERRQ(ierr); 855 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetIS_C","",PETSC_NULL);CHKERRQ(ierr); 856 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetType_C","",PETSC_NULL);CHKERRQ(ierr); 857 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetBlockSize_C","",PETSC_NULL);CHKERRQ(ierr); 858 PetscFunctionReturn(0); 859 } 860 861 #undef __FUNCT__ 862 #define __FUNCT__ "PCSetFromOptions_FieldSplit" 863 static PetscErrorCode PCSetFromOptions_FieldSplit(PC pc) 864 { 865 PetscErrorCode ierr; 866 PetscInt bs; 867 PetscBool flg,stokes = PETSC_FALSE; 868 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 869 PCCompositeType ctype; 870 871 PetscFunctionBegin; 872 ierr = PetscOptionsHead("FieldSplit options");CHKERRQ(ierr); 873 ierr = PetscOptionsBool("-pc_fieldsplit_real_diagonal","Use diagonal blocks of the operator","PCFieldSplitSetRealDiagonal",jac->realdiagonal,&jac->realdiagonal,PETSC_NULL);CHKERRQ(ierr); 874 ierr = PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);CHKERRQ(ierr); 875 if (flg) { 876 ierr = PCFieldSplitSetBlockSize(pc,bs);CHKERRQ(ierr); 877 } 878 879 ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,PETSC_NULL);CHKERRQ(ierr); 880 if (stokes) { 881 ierr = PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);CHKERRQ(ierr); 882 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_SELF; 883 } 884 885 ierr = PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg);CHKERRQ(ierr); 886 if (flg) { 887 ierr = PCFieldSplitSetType(pc,ctype);CHKERRQ(ierr); 888 } 889 890 /* Only setup fields once */ 891 if ((jac->bs > 0) && (jac->nsplits == 0)) { 892 /* only allow user to set fields from command line if bs is already known. 893 otherwise user can set them in PCFieldSplitSetDefaults() */ 894 ierr = PCFieldSplitSetRuntimeSplits_Private(pc);CHKERRQ(ierr); 895 if (jac->splitdefined) {ierr = PetscInfo(pc,"Splits defined using the options database\n");CHKERRQ(ierr);} 896 } 897 if (jac->type == PC_COMPOSITE_SCHUR) { 898 ierr = PetscOptionsEnum("-pc_fieldsplit_schur_factorization_type","Factorization to use","None",PCFieldSplitSchurFactorizationTypes,(PetscEnum)jac->schurfactorization,(PetscEnum*)&jac->schurfactorization,PETSC_NULL);CHKERRQ(ierr); 899 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); 900 } 901 ierr = PetscOptionsTail();CHKERRQ(ierr); 902 PetscFunctionReturn(0); 903 } 904 905 /*------------------------------------------------------------------------------------*/ 906 907 EXTERN_C_BEGIN 908 #undef __FUNCT__ 909 #define __FUNCT__ "PCFieldSplitSetFields_FieldSplit" 910 PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields) 911 { 912 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 913 PetscErrorCode ierr; 914 PC_FieldSplitLink ilink,next = jac->head; 915 char prefix[128]; 916 PetscInt i; 917 918 PetscFunctionBegin; 919 if (jac->splitdefined) { 920 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 921 PetscFunctionReturn(0); 922 } 923 for (i=0; i<n; i++) { 924 if (fields[i] >= jac->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs); 925 if (fields[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]); 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 ierr = KSPCreate(((PetscObject)pc)->comm,&ilink->ksp);CHKERRQ(ierr); 939 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 940 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 941 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 942 943 ierr = PetscSNPrintf(prefix,sizeof prefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",ilink->splitname);CHKERRQ(ierr); 944 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 945 946 if (!next) { 947 jac->head = ilink; 948 ilink->previous = PETSC_NULL; 949 } else { 950 while (next->next) { 951 next = next->next; 952 } 953 next->next = ilink; 954 ilink->previous = next; 955 } 956 jac->nsplits++; 957 PetscFunctionReturn(0); 958 } 959 EXTERN_C_END 960 961 EXTERN_C_BEGIN 962 #undef __FUNCT__ 963 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit_Schur" 964 PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp) 965 { 966 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 967 PetscErrorCode ierr; 968 969 PetscFunctionBegin; 970 ierr = PetscMalloc(jac->nsplits*sizeof(KSP),subksp);CHKERRQ(ierr); 971 ierr = MatSchurComplementGetKSP(jac->schur,*subksp);CHKERRQ(ierr); 972 (*subksp)[1] = jac->kspschur; 973 if (n) *n = jac->nsplits; 974 PetscFunctionReturn(0); 975 } 976 EXTERN_C_END 977 978 EXTERN_C_BEGIN 979 #undef __FUNCT__ 980 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit" 981 PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp) 982 { 983 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 984 PetscErrorCode ierr; 985 PetscInt cnt = 0; 986 PC_FieldSplitLink ilink = jac->head; 987 988 PetscFunctionBegin; 989 ierr = PetscMalloc(jac->nsplits*sizeof(KSP),subksp);CHKERRQ(ierr); 990 while (ilink) { 991 (*subksp)[cnt++] = ilink->ksp; 992 ilink = ilink->next; 993 } 994 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); 995 if (n) *n = jac->nsplits; 996 PetscFunctionReturn(0); 997 } 998 EXTERN_C_END 999 1000 EXTERN_C_BEGIN 1001 #undef __FUNCT__ 1002 #define __FUNCT__ "PCFieldSplitSetIS_FieldSplit" 1003 PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is) 1004 { 1005 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1006 PetscErrorCode ierr; 1007 PC_FieldSplitLink ilink, next = jac->head; 1008 char prefix[128]; 1009 1010 PetscFunctionBegin; 1011 if (jac->splitdefined) { 1012 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 1013 PetscFunctionReturn(0); 1014 } 1015 ierr = PetscNew(struct _PC_FieldSplitLink,&ilink);CHKERRQ(ierr); 1016 if (splitname) { 1017 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 1018 } else { 1019 ierr = PetscMalloc(3*sizeof(char),&ilink->splitname);CHKERRQ(ierr); 1020 ierr = PetscSNPrintf(ilink->splitname,3,"%D",jac->nsplits);CHKERRQ(ierr); 1021 } 1022 ilink->is = is; 1023 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 1024 ilink->next = PETSC_NULL; 1025 ierr = KSPCreate(((PetscObject)pc)->comm,&ilink->ksp);CHKERRQ(ierr); 1026 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 1027 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 1028 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 1029 1030 ierr = PetscSNPrintf(prefix,sizeof prefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",ilink->splitname);CHKERRQ(ierr); 1031 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 1032 1033 if (!next) { 1034 jac->head = ilink; 1035 ilink->previous = PETSC_NULL; 1036 } else { 1037 while (next->next) { 1038 next = next->next; 1039 } 1040 next->next = ilink; 1041 ilink->previous = next; 1042 } 1043 jac->nsplits++; 1044 1045 PetscFunctionReturn(0); 1046 } 1047 EXTERN_C_END 1048 1049 #undef __FUNCT__ 1050 #define __FUNCT__ "PCFieldSplitSetFields" 1051 /*@ 1052 PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner 1053 1054 Logically Collective on PC 1055 1056 Input Parameters: 1057 + pc - the preconditioner context 1058 . splitname - name of this split, if PETSC_NULL the number of the split is used 1059 . n - the number of fields in this split 1060 - fields - the fields in this split 1061 1062 Level: intermediate 1063 1064 Notes: Use PCFieldSplitSetIS() to set a completely general set of indices as a field. 1065 1066 The PCFieldSplitSetFields() is for defining fields as a strided blocks. For example, if the block 1067 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 1068 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x.... 1069 where the numbered entries indicate what is in the field. 1070 1071 This function is called once per split (it creates a new split each time). Solve options 1072 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1073 1074 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize(), PCFieldSplitSetIS() 1075 1076 @*/ 1077 PetscErrorCode PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields) 1078 { 1079 PetscErrorCode ierr; 1080 1081 PetscFunctionBegin; 1082 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1083 PetscValidCharPointer(splitname,2); 1084 if (n < 1) SETERRQ2(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Provided number of fields %D in split \"%s\" not positive",n,splitname); 1085 PetscValidIntPointer(fields,3); 1086 ierr = PetscTryMethod(pc,"PCFieldSplitSetFields_C",(PC,const char[],PetscInt,const PetscInt *),(pc,splitname,n,fields));CHKERRQ(ierr); 1087 PetscFunctionReturn(0); 1088 } 1089 1090 #undef __FUNCT__ 1091 #define __FUNCT__ "PCFieldSplitSetIS" 1092 /*@ 1093 PCFieldSplitSetIS - Sets the exact elements for field 1094 1095 Logically Collective on PC 1096 1097 Input Parameters: 1098 + pc - the preconditioner context 1099 . splitname - name of this split, if PETSC_NULL the number of the split is used 1100 - is - the index set that defines the vector elements in this field 1101 1102 1103 Notes: 1104 Use PCFieldSplitSetFields(), for fields defined by strided types. 1105 1106 This function is called once per split (it creates a new split each time). Solve options 1107 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1108 1109 Level: intermediate 1110 1111 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize() 1112 1113 @*/ 1114 PetscErrorCode PCFieldSplitSetIS(PC pc,const char splitname[],IS is) 1115 { 1116 PetscErrorCode ierr; 1117 1118 PetscFunctionBegin; 1119 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1120 if (splitname) PetscValidCharPointer(splitname,2); 1121 PetscValidHeaderSpecific(is,IS_CLASSID,3); 1122 ierr = PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));CHKERRQ(ierr); 1123 PetscFunctionReturn(0); 1124 } 1125 1126 #undef __FUNCT__ 1127 #define __FUNCT__ "PCFieldSplitGetIS" 1128 /*@ 1129 PCFieldSplitGetIS - Retrieves the elements for a field as an IS 1130 1131 Logically Collective on PC 1132 1133 Input Parameters: 1134 + pc - the preconditioner context 1135 - splitname - name of this split 1136 1137 Output Parameter: 1138 - is - the index set that defines the vector elements in this field, or PETSC_NULL if the field is not found 1139 1140 Level: intermediate 1141 1142 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetIS() 1143 1144 @*/ 1145 PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is) 1146 { 1147 PetscErrorCode ierr; 1148 1149 PetscFunctionBegin; 1150 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1151 PetscValidCharPointer(splitname,2); 1152 PetscValidPointer(is,3); 1153 { 1154 PC_FieldSplit *jac = (PC_FieldSplit *) pc->data; 1155 PC_FieldSplitLink ilink = jac->head; 1156 PetscBool found; 1157 1158 *is = PETSC_NULL; 1159 while(ilink) { 1160 ierr = PetscStrcmp(ilink->splitname, splitname, &found);CHKERRQ(ierr); 1161 if (found) { 1162 *is = ilink->is; 1163 break; 1164 } 1165 ilink = ilink->next; 1166 } 1167 } 1168 PetscFunctionReturn(0); 1169 } 1170 1171 #undef __FUNCT__ 1172 #define __FUNCT__ "PCFieldSplitSetBlockSize" 1173 /*@ 1174 PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the 1175 fieldsplit preconditioner. If not set the matrix block size is used. 1176 1177 Logically Collective on PC 1178 1179 Input Parameters: 1180 + pc - the preconditioner context 1181 - bs - the block size 1182 1183 Level: intermediate 1184 1185 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields() 1186 1187 @*/ 1188 PetscErrorCode PCFieldSplitSetBlockSize(PC pc,PetscInt bs) 1189 { 1190 PetscErrorCode ierr; 1191 1192 PetscFunctionBegin; 1193 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1194 PetscValidLogicalCollectiveInt(pc,bs,2); 1195 ierr = PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));CHKERRQ(ierr); 1196 PetscFunctionReturn(0); 1197 } 1198 1199 #undef __FUNCT__ 1200 #define __FUNCT__ "PCFieldSplitGetSubKSP" 1201 /*@C 1202 PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits 1203 1204 Collective on KSP 1205 1206 Input Parameter: 1207 . pc - the preconditioner context 1208 1209 Output Parameters: 1210 + n - the number of splits 1211 - pc - the array of KSP contexts 1212 1213 Note: 1214 After PCFieldSplitGetSubKSP() the array of KSPs IS to be freed by the user 1215 (not the KSP just the array that contains them). 1216 1217 You must call KSPSetUp() before calling PCFieldSplitGetSubKSP(). 1218 1219 Level: advanced 1220 1221 .seealso: PCFIELDSPLIT 1222 @*/ 1223 PetscErrorCode PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[]) 1224 { 1225 PetscErrorCode ierr; 1226 1227 PetscFunctionBegin; 1228 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1229 if (n) PetscValidIntPointer(n,2); 1230 ierr = PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));CHKERRQ(ierr); 1231 PetscFunctionReturn(0); 1232 } 1233 1234 #undef __FUNCT__ 1235 #define __FUNCT__ "PCFieldSplitSchurPrecondition" 1236 /*@ 1237 PCFieldSplitSchurPrecondition - Indicates if the Schur complement is preconditioned by a preconditioner constructed by the 1238 A11 matrix. Otherwise no preconditioner is used. 1239 1240 Collective on PC 1241 1242 Input Parameters: 1243 + pc - the preconditioner context 1244 . ptype - which matrix to use for preconditioning the Schur complement, PC_FIELDSPLIT_SCHUR_PRE_DIAG (diag) is default 1245 - userpre - matrix to use for preconditioning, or PETSC_NULL 1246 1247 Options Database: 1248 . -pc_fieldsplit_schur_precondition <self,user,diag> default is diag 1249 1250 Notes: 1251 $ If ptype is 1252 $ user then the preconditioner for the Schur complement is generated by the provided matrix (pre argument 1253 $ to this function). 1254 $ diag then the preconditioner for the Schur complement is generated by the block diagonal part of the original 1255 $ matrix associated with the Schur complement (i.e. A11) 1256 $ self the preconditioner for the Schur complement is generated from the Schur complement matrix itself: 1257 $ The only preconditioner that currently works directly with the Schur complement matrix object is the PCLSC 1258 $ preconditioner 1259 1260 When solving a saddle point problem, where the A11 block is identically zero, using diag as the ptype only makes sense 1261 with the additional option -fieldsplit_1_pc_type none. Usually for saddle point problems one would use a ptype of self and 1262 -fieldsplit_1_pc_type lsc which uses the least squares commutator compute a preconditioner for the Schur complement. 1263 1264 Developer Notes: This is a terrible name, gives no good indication of what the function does and should also have Set in 1265 the name since it sets a proceedure to use. 1266 1267 Level: intermediate 1268 1269 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCLSC 1270 1271 @*/ 1272 PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1273 { 1274 PetscErrorCode ierr; 1275 1276 PetscFunctionBegin; 1277 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1278 ierr = PetscTryMethod(pc,"PCFieldSplitSchurPrecondition_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));CHKERRQ(ierr); 1279 PetscFunctionReturn(0); 1280 } 1281 1282 EXTERN_C_BEGIN 1283 #undef __FUNCT__ 1284 #define __FUNCT__ "PCFieldSplitSchurPrecondition_FieldSplit" 1285 PetscErrorCode PCFieldSplitSchurPrecondition_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1286 { 1287 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1288 PetscErrorCode ierr; 1289 1290 PetscFunctionBegin; 1291 jac->schurpre = ptype; 1292 if (pre) { 1293 ierr = MatDestroy(&jac->schur_user);CHKERRQ(ierr); 1294 jac->schur_user = pre; 1295 ierr = PetscObjectReference((PetscObject)jac->schur_user);CHKERRQ(ierr); 1296 } 1297 PetscFunctionReturn(0); 1298 } 1299 EXTERN_C_END 1300 1301 #undef __FUNCT__ 1302 #define __FUNCT__ "PCFieldSplitGetSchurBlocks" 1303 /*@C 1304 PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement 1305 1306 Collective on KSP 1307 1308 Input Parameter: 1309 . pc - the preconditioner context 1310 1311 Output Parameters: 1312 + A00 - the (0,0) block 1313 . A01 - the (0,1) block 1314 . A10 - the (1,0) block 1315 - A11 - the (1,1) block 1316 1317 Level: advanced 1318 1319 .seealso: PCFIELDSPLIT 1320 @*/ 1321 PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11) 1322 { 1323 PC_FieldSplit *jac = (PC_FieldSplit *) pc->data; 1324 1325 PetscFunctionBegin; 1326 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1327 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach."); 1328 if (A00) *A00 = jac->pmat[0]; 1329 if (A01) *A01 = jac->B; 1330 if (A10) *A10 = jac->C; 1331 if (A11) *A11 = jac->pmat[1]; 1332 PetscFunctionReturn(0); 1333 } 1334 1335 EXTERN_C_BEGIN 1336 #undef __FUNCT__ 1337 #define __FUNCT__ "PCFieldSplitSetType_FieldSplit" 1338 PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type) 1339 { 1340 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1341 PetscErrorCode ierr; 1342 1343 PetscFunctionBegin; 1344 jac->type = type; 1345 if (type == PC_COMPOSITE_SCHUR) { 1346 pc->ops->apply = PCApply_FieldSplit_Schur; 1347 pc->ops->view = PCView_FieldSplit_Schur; 1348 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit_Schur",PCFieldSplitGetSubKSP_FieldSplit_Schur);CHKERRQ(ierr); 1349 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSchurPrecondition_C","PCFieldSplitSchurPrecondition_FieldSplit",PCFieldSplitSchurPrecondition_FieldSplit);CHKERRQ(ierr); 1350 1351 } else { 1352 pc->ops->apply = PCApply_FieldSplit; 1353 pc->ops->view = PCView_FieldSplit; 1354 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit",PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 1355 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSchurPrecondition_C","",0);CHKERRQ(ierr); 1356 } 1357 PetscFunctionReturn(0); 1358 } 1359 EXTERN_C_END 1360 1361 EXTERN_C_BEGIN 1362 #undef __FUNCT__ 1363 #define __FUNCT__ "PCFieldSplitSetBlockSize_FieldSplit" 1364 PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs) 1365 { 1366 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1367 1368 PetscFunctionBegin; 1369 if (bs < 1) SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs); 1370 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); 1371 jac->bs = bs; 1372 PetscFunctionReturn(0); 1373 } 1374 EXTERN_C_END 1375 1376 #undef __FUNCT__ 1377 #define __FUNCT__ "PCFieldSplitSetType" 1378 /*@ 1379 PCFieldSplitSetType - Sets the type of fieldsplit preconditioner. 1380 1381 Collective on PC 1382 1383 Input Parameter: 1384 . pc - the preconditioner context 1385 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 1386 1387 Options Database Key: 1388 . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type 1389 1390 Level: Developer 1391 1392 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 1393 1394 .seealso: PCCompositeSetType() 1395 1396 @*/ 1397 PetscErrorCode PCFieldSplitSetType(PC pc,PCCompositeType type) 1398 { 1399 PetscErrorCode ierr; 1400 1401 PetscFunctionBegin; 1402 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1403 ierr = PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));CHKERRQ(ierr); 1404 PetscFunctionReturn(0); 1405 } 1406 1407 /* -------------------------------------------------------------------------------------*/ 1408 /*MC 1409 PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual 1410 fields or groups of fields. See the users manual section "Solving Block Matrices" for more details. 1411 1412 To set options on the solvers for each block append -fieldsplit_ to all the PC 1413 options database keys. For example, -fieldsplit_pc_type ilu -fieldsplit_pc_factor_levels 1 1414 1415 To set the options on the solvers separate for each block call PCFieldSplitGetSubKSP() 1416 and set the options directly on the resulting KSP object 1417 1418 Level: intermediate 1419 1420 Options Database Keys: 1421 + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split 1422 . -pc_fieldsplit_default - automatically add any fields to additional splits that have not 1423 been supplied explicitly by -pc_fieldsplit_%d_fields 1424 . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields) 1425 . -pc_fieldsplit_type <additive,multiplicative,schur,symmetric_multiplicative> 1426 . -pc_fieldsplit_schur_precondition <true,false> default is true 1427 . -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero or negative diagonal and uses Schur complement with no preconditioner as the solver 1428 1429 - Options prefix for inner solvers when using Schur complement preconditioner are -fieldsplit_0_ and -fieldsplit_1_ 1430 for all other solvers they are -fieldsplit_%d_ for the dth field, use -fieldsplit_ for all fields 1431 1432 Notes: use PCFieldSplitSetFields() to set fields defined by "strided" entries and PCFieldSplitSetIS() 1433 to define a field by an arbitrary collection of entries. 1434 1435 If no fields are set the default is used. The fields are defined by entries strided by bs, 1436 beginning at 0 then 1, etc to bs-1. The block size can be set with PCFieldSplitSetBlockSize(), 1437 if this is not called the block size defaults to the blocksize of the second matrix passed 1438 to KSPSetOperators()/PCSetOperators(). 1439 1440 For the Schur complement preconditioner if J = ( A00 A01 ) 1441 ( A10 A11 ) 1442 the preconditioner is 1443 (I -B ksp(A00)) ( inv(A00) 0 (I 0 ) 1444 (0 I ) ( 0 ksp(S) ) (-A10 ksp(A00) I ) 1445 where the action of inv(A00) is applied using the KSP solver with prefix -fieldsplit_0_. The action of 1446 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). 1447 For PCFieldSplitGetKSP() when field number is 1448 0 it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default 1449 A11 is used to construct a preconditioner for S, use PCFieldSplitSchurPrecondition() to turn on or off this 1450 option. You can use the preconditioner PCLSC to precondition the Schur complement with -fieldsplit_1_pc_type lsc 1451 1452 If only one set of indices (one IS) is provided with PCFieldSplitSetIS() then the complement of that IS 1453 is used automatically for a second block. 1454 1455 The fieldsplit preconditioner cannot currently be used with the BAIJ or SBAIJ data formats if the blocksize is larger than 1. 1456 Generally it should be used with the AIJ format. 1457 1458 The forms of these preconditioners are closely related if not identical to forms derived as "Distributive Iterations", see, 1459 for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling. Note that one can also use PCFIELDSPLIT 1460 inside a smoother resulting in "Distributive Smoothers". 1461 1462 Concepts: physics based preconditioners, block preconditioners 1463 1464 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners, PCLSC, 1465 PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSchurPrecondition() 1466 M*/ 1467 1468 EXTERN_C_BEGIN 1469 #undef __FUNCT__ 1470 #define __FUNCT__ "PCCreate_FieldSplit" 1471 PetscErrorCode PCCreate_FieldSplit(PC pc) 1472 { 1473 PetscErrorCode ierr; 1474 PC_FieldSplit *jac; 1475 1476 PetscFunctionBegin; 1477 ierr = PetscNewLog(pc,PC_FieldSplit,&jac);CHKERRQ(ierr); 1478 jac->bs = -1; 1479 jac->nsplits = 0; 1480 jac->type = PC_COMPOSITE_MULTIPLICATIVE; 1481 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */ 1482 jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACTORIZATION_FULL; 1483 1484 pc->data = (void*)jac; 1485 1486 pc->ops->apply = PCApply_FieldSplit; 1487 pc->ops->applytranspose = PCApplyTranspose_FieldSplit; 1488 pc->ops->setup = PCSetUp_FieldSplit; 1489 pc->ops->reset = PCReset_FieldSplit; 1490 pc->ops->destroy = PCDestroy_FieldSplit; 1491 pc->ops->setfromoptions = PCSetFromOptions_FieldSplit; 1492 pc->ops->view = PCView_FieldSplit; 1493 pc->ops->applyrichardson = 0; 1494 1495 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit", 1496 PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 1497 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetFields_C","PCFieldSplitSetFields_FieldSplit", 1498 PCFieldSplitSetFields_FieldSplit);CHKERRQ(ierr); 1499 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetIS_C","PCFieldSplitSetIS_FieldSplit", 1500 PCFieldSplitSetIS_FieldSplit);CHKERRQ(ierr); 1501 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetType_C","PCFieldSplitSetType_FieldSplit", 1502 PCFieldSplitSetType_FieldSplit);CHKERRQ(ierr); 1503 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetBlockSize_C","PCFieldSplitSetBlockSize_FieldSplit", 1504 PCFieldSplitSetBlockSize_FieldSplit);CHKERRQ(ierr); 1505 PetscFunctionReturn(0); 1506 } 1507 EXTERN_C_END 1508 1509 1510 1511