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 PetscFunctionReturn(0); 854 } 855 856 #undef __FUNCT__ 857 #define __FUNCT__ "PCSetFromOptions_FieldSplit" 858 static PetscErrorCode PCSetFromOptions_FieldSplit(PC pc) 859 { 860 PetscErrorCode ierr; 861 PetscInt bs; 862 PetscBool flg,stokes = PETSC_FALSE; 863 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 864 PCCompositeType ctype; 865 866 PetscFunctionBegin; 867 ierr = PetscOptionsHead("FieldSplit options");CHKERRQ(ierr); 868 ierr = PetscOptionsBool("-pc_fieldsplit_real_diagonal","Use diagonal blocks of the operator","PCFieldSplitSetRealDiagonal",jac->realdiagonal,&jac->realdiagonal,PETSC_NULL);CHKERRQ(ierr); 869 ierr = PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);CHKERRQ(ierr); 870 if (flg) { 871 ierr = PCFieldSplitSetBlockSize(pc,bs);CHKERRQ(ierr); 872 } 873 874 ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,PETSC_NULL);CHKERRQ(ierr); 875 if (stokes) { 876 ierr = PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);CHKERRQ(ierr); 877 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_SELF; 878 } 879 880 ierr = PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg);CHKERRQ(ierr); 881 if (flg) { 882 ierr = PCFieldSplitSetType(pc,ctype);CHKERRQ(ierr); 883 } 884 885 /* Only setup fields once */ 886 if ((jac->bs > 0) && (jac->nsplits == 0)) { 887 /* only allow user to set fields from command line if bs is already known. 888 otherwise user can set them in PCFieldSplitSetDefaults() */ 889 ierr = PCFieldSplitSetRuntimeSplits_Private(pc);CHKERRQ(ierr); 890 if (jac->splitdefined) {ierr = PetscInfo(pc,"Splits defined using the options database\n");CHKERRQ(ierr);} 891 } 892 if (jac->type == PC_COMPOSITE_SCHUR) { 893 ierr = PetscOptionsEnum("-pc_fieldsplit_schur_factorization_type","Factorization to use","None",PCFieldSplitSchurFactorizationTypes,(PetscEnum)jac->schurfactorization,(PetscEnum*)&jac->schurfactorization,PETSC_NULL);CHKERRQ(ierr); 894 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); 895 } 896 ierr = PetscOptionsTail();CHKERRQ(ierr); 897 PetscFunctionReturn(0); 898 } 899 900 /*------------------------------------------------------------------------------------*/ 901 902 EXTERN_C_BEGIN 903 #undef __FUNCT__ 904 #define __FUNCT__ "PCFieldSplitSetFields_FieldSplit" 905 PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields) 906 { 907 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 908 PetscErrorCode ierr; 909 PC_FieldSplitLink ilink,next = jac->head; 910 char prefix[128]; 911 PetscInt i; 912 913 PetscFunctionBegin; 914 if (jac->splitdefined) { 915 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 916 PetscFunctionReturn(0); 917 } 918 for (i=0; i<n; i++) { 919 if (fields[i] >= jac->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs); 920 if (fields[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]); 921 } 922 ierr = PetscNew(struct _PC_FieldSplitLink,&ilink);CHKERRQ(ierr); 923 if (splitname) { 924 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 925 } else { 926 ierr = PetscMalloc(3*sizeof(char),&ilink->splitname);CHKERRQ(ierr); 927 ierr = PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);CHKERRQ(ierr); 928 } 929 ierr = PetscMalloc(n*sizeof(PetscInt),&ilink->fields);CHKERRQ(ierr); 930 ierr = PetscMemcpy(ilink->fields,fields,n*sizeof(PetscInt));CHKERRQ(ierr); 931 ilink->nfields = n; 932 ilink->next = PETSC_NULL; 933 ierr = KSPCreate(((PetscObject)pc)->comm,&ilink->ksp);CHKERRQ(ierr); 934 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 935 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 936 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 937 938 ierr = PetscSNPrintf(prefix,sizeof prefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",ilink->splitname);CHKERRQ(ierr); 939 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 940 941 if (!next) { 942 jac->head = ilink; 943 ilink->previous = PETSC_NULL; 944 } else { 945 while (next->next) { 946 next = next->next; 947 } 948 next->next = ilink; 949 ilink->previous = next; 950 } 951 jac->nsplits++; 952 PetscFunctionReturn(0); 953 } 954 EXTERN_C_END 955 956 EXTERN_C_BEGIN 957 #undef __FUNCT__ 958 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit_Schur" 959 PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp) 960 { 961 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 962 PetscErrorCode ierr; 963 964 PetscFunctionBegin; 965 ierr = PetscMalloc(jac->nsplits*sizeof(KSP),subksp);CHKERRQ(ierr); 966 ierr = MatSchurComplementGetKSP(jac->schur,*subksp);CHKERRQ(ierr); 967 (*subksp)[1] = jac->kspschur; 968 if (n) *n = jac->nsplits; 969 PetscFunctionReturn(0); 970 } 971 EXTERN_C_END 972 973 EXTERN_C_BEGIN 974 #undef __FUNCT__ 975 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit" 976 PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp) 977 { 978 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 979 PetscErrorCode ierr; 980 PetscInt cnt = 0; 981 PC_FieldSplitLink ilink = jac->head; 982 983 PetscFunctionBegin; 984 ierr = PetscMalloc(jac->nsplits*sizeof(KSP),subksp);CHKERRQ(ierr); 985 while (ilink) { 986 (*subksp)[cnt++] = ilink->ksp; 987 ilink = ilink->next; 988 } 989 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); 990 if (n) *n = jac->nsplits; 991 PetscFunctionReturn(0); 992 } 993 EXTERN_C_END 994 995 EXTERN_C_BEGIN 996 #undef __FUNCT__ 997 #define __FUNCT__ "PCFieldSplitSetIS_FieldSplit" 998 PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is) 999 { 1000 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1001 PetscErrorCode ierr; 1002 PC_FieldSplitLink ilink, next = jac->head; 1003 char prefix[128]; 1004 1005 PetscFunctionBegin; 1006 if (jac->splitdefined) { 1007 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 1008 PetscFunctionReturn(0); 1009 } 1010 ierr = PetscNew(struct _PC_FieldSplitLink,&ilink);CHKERRQ(ierr); 1011 if (splitname) { 1012 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 1013 } else { 1014 ierr = PetscMalloc(3*sizeof(char),&ilink->splitname);CHKERRQ(ierr); 1015 ierr = PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);CHKERRQ(ierr); 1016 } 1017 ilink->is = is; 1018 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 1019 ilink->next = PETSC_NULL; 1020 ierr = KSPCreate(((PetscObject)pc)->comm,&ilink->ksp);CHKERRQ(ierr); 1021 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 1022 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 1023 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 1024 1025 ierr = PetscSNPrintf(prefix,sizeof prefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",ilink->splitname);CHKERRQ(ierr); 1026 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 1027 1028 if (!next) { 1029 jac->head = ilink; 1030 ilink->previous = PETSC_NULL; 1031 } else { 1032 while (next->next) { 1033 next = next->next; 1034 } 1035 next->next = ilink; 1036 ilink->previous = next; 1037 } 1038 jac->nsplits++; 1039 1040 PetscFunctionReturn(0); 1041 } 1042 EXTERN_C_END 1043 1044 #undef __FUNCT__ 1045 #define __FUNCT__ "PCFieldSplitSetFields" 1046 /*@ 1047 PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner 1048 1049 Logically Collective on PC 1050 1051 Input Parameters: 1052 + pc - the preconditioner context 1053 . splitname - name of this split, if PETSC_NULL the number of the split is used 1054 . n - the number of fields in this split 1055 - fields - the fields in this split 1056 1057 Level: intermediate 1058 1059 Notes: Use PCFieldSplitSetIS() to set a completely general set of indices as a field. 1060 1061 The PCFieldSplitSetFields() is for defining fields as a strided blocks. For example, if the block 1062 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 1063 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x.... 1064 where the numbered entries indicate what is in the field. 1065 1066 This function is called once per split (it creates a new split each time). Solve options 1067 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1068 1069 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize(), PCFieldSplitSetIS() 1070 1071 @*/ 1072 PetscErrorCode PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields) 1073 { 1074 PetscErrorCode ierr; 1075 1076 PetscFunctionBegin; 1077 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1078 PetscValidCharPointer(splitname,2); 1079 if (n < 1) SETERRQ2(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Provided number of fields %D in split \"%s\" not positive",n,splitname); 1080 PetscValidIntPointer(fields,3); 1081 ierr = PetscTryMethod(pc,"PCFieldSplitSetFields_C",(PC,const char[],PetscInt,const PetscInt *),(pc,splitname,n,fields));CHKERRQ(ierr); 1082 PetscFunctionReturn(0); 1083 } 1084 1085 #undef __FUNCT__ 1086 #define __FUNCT__ "PCFieldSplitSetIS" 1087 /*@ 1088 PCFieldSplitSetIS - Sets the exact elements for field 1089 1090 Logically Collective on PC 1091 1092 Input Parameters: 1093 + pc - the preconditioner context 1094 . splitname - name of this split, if PETSC_NULL the number of the split is used 1095 - is - the index set that defines the vector elements in this field 1096 1097 1098 Notes: 1099 Use PCFieldSplitSetFields(), for fields defined by strided types. 1100 1101 This function is called once per split (it creates a new split each time). Solve options 1102 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1103 1104 Level: intermediate 1105 1106 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize() 1107 1108 @*/ 1109 PetscErrorCode PCFieldSplitSetIS(PC pc,const char splitname[],IS is) 1110 { 1111 PetscErrorCode ierr; 1112 1113 PetscFunctionBegin; 1114 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1115 PetscValidCharPointer(splitname,2); 1116 PetscValidHeaderSpecific(is,IS_CLASSID,3); 1117 ierr = PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));CHKERRQ(ierr); 1118 PetscFunctionReturn(0); 1119 } 1120 1121 #undef __FUNCT__ 1122 #define __FUNCT__ "PCFieldSplitGetIS" 1123 /*@ 1124 PCFieldSplitGetIS - Retrieves the elements for a field as an IS 1125 1126 Logically Collective on PC 1127 1128 Input Parameters: 1129 + pc - the preconditioner context 1130 - splitname - name of this split 1131 1132 Output Parameter: 1133 - is - the index set that defines the vector elements in this field, or PETSC_NULL if the field is not found 1134 1135 Level: intermediate 1136 1137 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetIS() 1138 1139 @*/ 1140 PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is) 1141 { 1142 PetscErrorCode ierr; 1143 1144 PetscFunctionBegin; 1145 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1146 PetscValidCharPointer(splitname,2); 1147 PetscValidPointer(is,3); 1148 { 1149 PC_FieldSplit *jac = (PC_FieldSplit *) pc->data; 1150 PC_FieldSplitLink ilink = jac->head; 1151 PetscBool found; 1152 1153 *is = PETSC_NULL; 1154 while(ilink) { 1155 ierr = PetscStrcmp(ilink->splitname, splitname, &found);CHKERRQ(ierr); 1156 if (found) { 1157 *is = ilink->is; 1158 break; 1159 } 1160 ilink = ilink->next; 1161 } 1162 } 1163 PetscFunctionReturn(0); 1164 } 1165 1166 #undef __FUNCT__ 1167 #define __FUNCT__ "PCFieldSplitSetBlockSize" 1168 /*@ 1169 PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the 1170 fieldsplit preconditioner. If not set the matrix block size is used. 1171 1172 Logically Collective on PC 1173 1174 Input Parameters: 1175 + pc - the preconditioner context 1176 - bs - the block size 1177 1178 Level: intermediate 1179 1180 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields() 1181 1182 @*/ 1183 PetscErrorCode PCFieldSplitSetBlockSize(PC pc,PetscInt bs) 1184 { 1185 PetscErrorCode ierr; 1186 1187 PetscFunctionBegin; 1188 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1189 PetscValidLogicalCollectiveInt(pc,bs,2); 1190 ierr = PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));CHKERRQ(ierr); 1191 PetscFunctionReturn(0); 1192 } 1193 1194 #undef __FUNCT__ 1195 #define __FUNCT__ "PCFieldSplitGetSubKSP" 1196 /*@C 1197 PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits 1198 1199 Collective on KSP 1200 1201 Input Parameter: 1202 . pc - the preconditioner context 1203 1204 Output Parameters: 1205 + n - the number of splits 1206 - pc - the array of KSP contexts 1207 1208 Note: 1209 After PCFieldSplitGetSubKSP() the array of KSPs IS to be freed by the user 1210 (not the KSP just the array that contains them). 1211 1212 You must call KSPSetUp() before calling PCFieldSplitGetSubKSP(). 1213 1214 Level: advanced 1215 1216 .seealso: PCFIELDSPLIT 1217 @*/ 1218 PetscErrorCode PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[]) 1219 { 1220 PetscErrorCode ierr; 1221 1222 PetscFunctionBegin; 1223 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1224 if (n) PetscValidIntPointer(n,2); 1225 ierr = PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));CHKERRQ(ierr); 1226 PetscFunctionReturn(0); 1227 } 1228 1229 #undef __FUNCT__ 1230 #define __FUNCT__ "PCFieldSplitSchurPrecondition" 1231 /*@ 1232 PCFieldSplitSchurPrecondition - Indicates if the Schur complement is preconditioned by a preconditioner constructed by the 1233 A11 matrix. Otherwise no preconditioner is used. 1234 1235 Collective on PC 1236 1237 Input Parameters: 1238 + pc - the preconditioner context 1239 . ptype - which matrix to use for preconditioning the Schur complement, PC_FIELDSPLIT_SCHUR_PRE_DIAG (diag) is default 1240 - userpre - matrix to use for preconditioning, or PETSC_NULL 1241 1242 Options Database: 1243 . -pc_fieldsplit_schur_precondition <self,user,diag> default is diag 1244 1245 Notes: 1246 $ If ptype is 1247 $ user then the preconditioner for the Schur complement is generated by the provided matrix (pre argument 1248 $ to this function). 1249 $ diag then the preconditioner for the Schur complement is generated by the block diagonal part of the original 1250 $ matrix associated with the Schur complement (i.e. A11) 1251 $ self the preconditioner for the Schur complement is generated from the Schur complement matrix itself: 1252 $ The only preconditioner that currently works directly with the Schur complement matrix object is the PCLSC 1253 $ preconditioner 1254 1255 When solving a saddle point problem, where the A11 block is identically zero, using diag as the ptype only makes sense 1256 with the additional option -fieldsplit_1_pc_type none. Usually for saddle point problems one would use a ptype of self and 1257 -fieldsplit_1_pc_type lsc which uses the least squares commutator compute a preconditioner for the Schur complement. 1258 1259 Developer Notes: This is a terrible name, gives no good indication of what the function does and should also have Set in 1260 the name since it sets a proceedure to use. 1261 1262 Level: intermediate 1263 1264 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCLSC 1265 1266 @*/ 1267 PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1268 { 1269 PetscErrorCode ierr; 1270 1271 PetscFunctionBegin; 1272 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1273 ierr = PetscTryMethod(pc,"PCFieldSplitSchurPrecondition_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));CHKERRQ(ierr); 1274 PetscFunctionReturn(0); 1275 } 1276 1277 EXTERN_C_BEGIN 1278 #undef __FUNCT__ 1279 #define __FUNCT__ "PCFieldSplitSchurPrecondition_FieldSplit" 1280 PetscErrorCode PCFieldSplitSchurPrecondition_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1281 { 1282 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1283 PetscErrorCode ierr; 1284 1285 PetscFunctionBegin; 1286 jac->schurpre = ptype; 1287 if (pre) { 1288 ierr = MatDestroy(&jac->schur_user);CHKERRQ(ierr); 1289 jac->schur_user = pre; 1290 ierr = PetscObjectReference((PetscObject)jac->schur_user);CHKERRQ(ierr); 1291 } 1292 PetscFunctionReturn(0); 1293 } 1294 EXTERN_C_END 1295 1296 #undef __FUNCT__ 1297 #define __FUNCT__ "PCFieldSplitGetSchurBlocks" 1298 /*@C 1299 PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement 1300 1301 Collective on KSP 1302 1303 Input Parameter: 1304 . pc - the preconditioner context 1305 1306 Output Parameters: 1307 + A00 - the (0,0) block 1308 . A01 - the (0,1) block 1309 . A10 - the (1,0) block 1310 - A11 - the (1,1) block 1311 1312 Level: advanced 1313 1314 .seealso: PCFIELDSPLIT 1315 @*/ 1316 PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11) 1317 { 1318 PC_FieldSplit *jac = (PC_FieldSplit *) pc->data; 1319 1320 PetscFunctionBegin; 1321 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1322 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach."); 1323 if (A00) *A00 = jac->pmat[0]; 1324 if (A01) *A01 = jac->B; 1325 if (A10) *A10 = jac->C; 1326 if (A11) *A11 = jac->pmat[1]; 1327 PetscFunctionReturn(0); 1328 } 1329 1330 EXTERN_C_BEGIN 1331 #undef __FUNCT__ 1332 #define __FUNCT__ "PCFieldSplitSetType_FieldSplit" 1333 PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type) 1334 { 1335 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1336 PetscErrorCode ierr; 1337 1338 PetscFunctionBegin; 1339 jac->type = type; 1340 if (type == PC_COMPOSITE_SCHUR) { 1341 pc->ops->apply = PCApply_FieldSplit_Schur; 1342 pc->ops->view = PCView_FieldSplit_Schur; 1343 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit_Schur",PCFieldSplitGetSubKSP_FieldSplit_Schur);CHKERRQ(ierr); 1344 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSchurPrecondition_C","PCFieldSplitSchurPrecondition_FieldSplit",PCFieldSplitSchurPrecondition_FieldSplit);CHKERRQ(ierr); 1345 1346 } else { 1347 pc->ops->apply = PCApply_FieldSplit; 1348 pc->ops->view = PCView_FieldSplit; 1349 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit",PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 1350 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSchurPrecondition_C","",0);CHKERRQ(ierr); 1351 } 1352 PetscFunctionReturn(0); 1353 } 1354 EXTERN_C_END 1355 1356 EXTERN_C_BEGIN 1357 #undef __FUNCT__ 1358 #define __FUNCT__ "PCFieldSplitSetBlockSize_FieldSplit" 1359 PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs) 1360 { 1361 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1362 1363 PetscFunctionBegin; 1364 if (bs < 1) SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs); 1365 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); 1366 jac->bs = bs; 1367 PetscFunctionReturn(0); 1368 } 1369 EXTERN_C_END 1370 1371 #undef __FUNCT__ 1372 #define __FUNCT__ "PCFieldSplitSetType" 1373 /*@ 1374 PCFieldSplitSetType - Sets the type of fieldsplit preconditioner. 1375 1376 Collective on PC 1377 1378 Input Parameter: 1379 . pc - the preconditioner context 1380 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 1381 1382 Options Database Key: 1383 . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type 1384 1385 Level: Developer 1386 1387 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 1388 1389 .seealso: PCCompositeSetType() 1390 1391 @*/ 1392 PetscErrorCode PCFieldSplitSetType(PC pc,PCCompositeType type) 1393 { 1394 PetscErrorCode ierr; 1395 1396 PetscFunctionBegin; 1397 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1398 ierr = PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));CHKERRQ(ierr); 1399 PetscFunctionReturn(0); 1400 } 1401 1402 /* -------------------------------------------------------------------------------------*/ 1403 /*MC 1404 PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual 1405 fields or groups of fields. See the users manual section "Solving Block Matrices" for more details. 1406 1407 To set options on the solvers for each block append -fieldsplit_ to all the PC 1408 options database keys. For example, -fieldsplit_pc_type ilu -fieldsplit_pc_factor_levels 1 1409 1410 To set the options on the solvers separate for each block call PCFieldSplitGetSubKSP() 1411 and set the options directly on the resulting KSP object 1412 1413 Level: intermediate 1414 1415 Options Database Keys: 1416 + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split 1417 . -pc_fieldsplit_default - automatically add any fields to additional splits that have not 1418 been supplied explicitly by -pc_fieldsplit_%d_fields 1419 . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields) 1420 . -pc_fieldsplit_type <additive,multiplicative,schur,symmetric_multiplicative> 1421 . -pc_fieldsplit_schur_precondition <true,false> default is true 1422 . -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero or negative diagonal and uses Schur complement with no preconditioner as the solver 1423 1424 - Options prefix for inner solvers when using Schur complement preconditioner are -fieldsplit_0_ and -fieldsplit_1_ 1425 for all other solvers they are -fieldsplit_%d_ for the dth field, use -fieldsplit_ for all fields 1426 1427 Notes: use PCFieldSplitSetFields() to set fields defined by "strided" entries and PCFieldSplitSetIS() 1428 to define a field by an arbitrary collection of entries. 1429 1430 If no fields are set the default is used. The fields are defined by entries strided by bs, 1431 beginning at 0 then 1, etc to bs-1. The block size can be set with PCFieldSplitSetBlockSize(), 1432 if this is not called the block size defaults to the blocksize of the second matrix passed 1433 to KSPSetOperators()/PCSetOperators(). 1434 1435 For the Schur complement preconditioner if J = ( A00 A01 ) 1436 ( A10 A11 ) 1437 the preconditioner is 1438 (I -B ksp(A00)) ( inv(A00) 0 (I 0 ) 1439 (0 I ) ( 0 ksp(S) ) (-A10 ksp(A00) I ) 1440 where the action of inv(A00) is applied using the KSP solver with prefix -fieldsplit_0_. The action of 1441 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). 1442 For PCFieldSplitGetKSP() when field number is 1443 0 it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default 1444 A11 is used to construct a preconditioner for S, use PCFieldSplitSchurPrecondition() to turn on or off this 1445 option. You can use the preconditioner PCLSC to precondition the Schur complement with -fieldsplit_1_pc_type lsc 1446 1447 If only one set of indices (one IS) is provided with PCFieldSplitSetIS() then the complement of that IS 1448 is used automatically for a second block. 1449 1450 The fieldsplit preconditioner cannot currently be used with the BAIJ or SBAIJ data formats if the blocksize is larger than 1. 1451 Generally it should be used with the AIJ format. 1452 1453 The forms of these preconditioners are closely related if not identical to forms derived as "Distributive Iterations", see, 1454 for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling. Note that one can also use PCFIELDSPLIT 1455 inside a smoother resulting in "Distributive Smoothers". 1456 1457 Concepts: physics based preconditioners, block preconditioners 1458 1459 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners, PCLSC, 1460 PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSchurPrecondition() 1461 M*/ 1462 1463 EXTERN_C_BEGIN 1464 #undef __FUNCT__ 1465 #define __FUNCT__ "PCCreate_FieldSplit" 1466 PetscErrorCode PCCreate_FieldSplit(PC pc) 1467 { 1468 PetscErrorCode ierr; 1469 PC_FieldSplit *jac; 1470 1471 PetscFunctionBegin; 1472 ierr = PetscNewLog(pc,PC_FieldSplit,&jac);CHKERRQ(ierr); 1473 jac->bs = -1; 1474 jac->nsplits = 0; 1475 jac->type = PC_COMPOSITE_MULTIPLICATIVE; 1476 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */ 1477 jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACTORIZATION_FULL; 1478 1479 pc->data = (void*)jac; 1480 1481 pc->ops->apply = PCApply_FieldSplit; 1482 pc->ops->applytranspose = PCApplyTranspose_FieldSplit; 1483 pc->ops->setup = PCSetUp_FieldSplit; 1484 pc->ops->reset = PCReset_FieldSplit; 1485 pc->ops->destroy = PCDestroy_FieldSplit; 1486 pc->ops->setfromoptions = PCSetFromOptions_FieldSplit; 1487 pc->ops->view = PCView_FieldSplit; 1488 pc->ops->applyrichardson = 0; 1489 1490 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit", 1491 PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 1492 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetFields_C","PCFieldSplitSetFields_FieldSplit", 1493 PCFieldSplitSetFields_FieldSplit);CHKERRQ(ierr); 1494 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetIS_C","PCFieldSplitSetIS_FieldSplit", 1495 PCFieldSplitSetIS_FieldSplit);CHKERRQ(ierr); 1496 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetType_C","PCFieldSplitSetType_FieldSplit", 1497 PCFieldSplitSetType_FieldSplit);CHKERRQ(ierr); 1498 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetBlockSize_C","PCFieldSplitSetBlockSize_FieldSplit", 1499 PCFieldSplitSetBlockSize_FieldSplit);CHKERRQ(ierr); 1500 PetscFunctionReturn(0); 1501 } 1502 EXTERN_C_END 1503 1504 1505 1506