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