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