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