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