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