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 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(8*sizeof(char),&ilink->splitname);CHKERRQ(ierr); 1018 ierr = PetscSNPrintf(ilink->splitname,7,"%D",jac->nsplits);CHKERRQ(ierr); 1019 } 1020 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 1021 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 1022 ilink->is = is; 1023 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 1024 ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr); 1025 ilink->is_col = is; 1026 ilink->next = PETSC_NULL; 1027 ierr = KSPCreate(((PetscObject)pc)->comm,&ilink->ksp);CHKERRQ(ierr); 1028 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 1029 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 1030 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 1031 1032 ierr = PetscSNPrintf(prefix,sizeof prefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",ilink->splitname);CHKERRQ(ierr); 1033 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 1034 1035 if (!next) { 1036 jac->head = ilink; 1037 ilink->previous = PETSC_NULL; 1038 } else { 1039 while (next->next) { 1040 next = next->next; 1041 } 1042 next->next = ilink; 1043 ilink->previous = next; 1044 } 1045 jac->nsplits++; 1046 1047 PetscFunctionReturn(0); 1048 } 1049 EXTERN_C_END 1050 1051 #undef __FUNCT__ 1052 #define __FUNCT__ "PCFieldSplitSetFields" 1053 /*@ 1054 PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner 1055 1056 Logically Collective on PC 1057 1058 Input Parameters: 1059 + pc - the preconditioner context 1060 . splitname - name of this split, if PETSC_NULL the number of the split is used 1061 . n - the number of fields in this split 1062 - fields - the fields in this split 1063 1064 Level: intermediate 1065 1066 Notes: Use PCFieldSplitSetIS() to set a completely general set of indices as a field. 1067 1068 The PCFieldSplitSetFields() is for defining fields as a strided blocks. For example, if the block 1069 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 1070 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x.... 1071 where the numbered entries indicate what is in the field. 1072 1073 This function is called once per split (it creates a new split each time). Solve options 1074 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1075 1076 Developer Note: This routine does not actually create the IS representing the split, that is delayed 1077 until PCSetUp_FieldSplit(), because information about the vector/matrix layouts may not be 1078 available when this routine is called. 1079 1080 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize(), PCFieldSplitSetIS() 1081 1082 @*/ 1083 PetscErrorCode PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 1084 { 1085 PetscErrorCode ierr; 1086 1087 PetscFunctionBegin; 1088 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1089 PetscValidCharPointer(splitname,2); 1090 if (n < 1) SETERRQ2(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Provided number of fields %D in split \"%s\" not positive",n,splitname); 1091 PetscValidIntPointer(fields,3); 1092 ierr = PetscTryMethod(pc,"PCFieldSplitSetFields_C",(PC,const char[],PetscInt,const PetscInt *,const PetscInt *),(pc,splitname,n,fields,fields_col));CHKERRQ(ierr); 1093 PetscFunctionReturn(0); 1094 } 1095 1096 #undef __FUNCT__ 1097 #define __FUNCT__ "PCFieldSplitSetIS" 1098 /*@ 1099 PCFieldSplitSetIS - Sets the exact elements for field 1100 1101 Logically Collective on PC 1102 1103 Input Parameters: 1104 + pc - the preconditioner context 1105 . splitname - name of this split, if PETSC_NULL the number of the split is used 1106 - is - the index set that defines the vector elements in this field 1107 1108 1109 Notes: 1110 Use PCFieldSplitSetFields(), for fields defined by strided types. 1111 1112 This function is called once per split (it creates a new split each time). Solve options 1113 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1114 1115 Level: intermediate 1116 1117 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize() 1118 1119 @*/ 1120 PetscErrorCode PCFieldSplitSetIS(PC pc,const char splitname[],IS is) 1121 { 1122 PetscErrorCode ierr; 1123 1124 PetscFunctionBegin; 1125 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1126 PetscValidCharPointer(splitname,2); 1127 PetscValidHeaderSpecific(is,IS_CLASSID,3); 1128 ierr = PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));CHKERRQ(ierr); 1129 PetscFunctionReturn(0); 1130 } 1131 1132 #undef __FUNCT__ 1133 #define __FUNCT__ "PCFieldSplitGetIS" 1134 /*@ 1135 PCFieldSplitGetIS - Retrieves the elements for a field as an IS 1136 1137 Logically Collective on PC 1138 1139 Input Parameters: 1140 + pc - the preconditioner context 1141 - splitname - name of this split 1142 1143 Output Parameter: 1144 - is - the index set that defines the vector elements in this field, or PETSC_NULL if the field is not found 1145 1146 Level: intermediate 1147 1148 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetIS() 1149 1150 @*/ 1151 PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is) 1152 { 1153 PetscErrorCode ierr; 1154 1155 PetscFunctionBegin; 1156 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1157 PetscValidCharPointer(splitname,2); 1158 PetscValidPointer(is,3); 1159 { 1160 PC_FieldSplit *jac = (PC_FieldSplit *) pc->data; 1161 PC_FieldSplitLink ilink = jac->head; 1162 PetscBool found; 1163 1164 *is = PETSC_NULL; 1165 while(ilink) { 1166 ierr = PetscStrcmp(ilink->splitname, splitname, &found);CHKERRQ(ierr); 1167 if (found) { 1168 *is = ilink->is; 1169 break; 1170 } 1171 ilink = ilink->next; 1172 } 1173 } 1174 PetscFunctionReturn(0); 1175 } 1176 1177 #undef __FUNCT__ 1178 #define __FUNCT__ "PCFieldSplitSetBlockSize" 1179 /*@ 1180 PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the 1181 fieldsplit preconditioner. If not set the matrix block size is used. 1182 1183 Logically Collective on PC 1184 1185 Input Parameters: 1186 + pc - the preconditioner context 1187 - bs - the block size 1188 1189 Level: intermediate 1190 1191 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields() 1192 1193 @*/ 1194 PetscErrorCode PCFieldSplitSetBlockSize(PC pc,PetscInt bs) 1195 { 1196 PetscErrorCode ierr; 1197 1198 PetscFunctionBegin; 1199 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1200 PetscValidLogicalCollectiveInt(pc,bs,2); 1201 ierr = PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));CHKERRQ(ierr); 1202 PetscFunctionReturn(0); 1203 } 1204 1205 #undef __FUNCT__ 1206 #define __FUNCT__ "PCFieldSplitGetSubKSP" 1207 /*@C 1208 PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits 1209 1210 Collective on KSP 1211 1212 Input Parameter: 1213 . pc - the preconditioner context 1214 1215 Output Parameters: 1216 + n - the number of splits 1217 - pc - the array of KSP contexts 1218 1219 Note: 1220 After PCFieldSplitGetSubKSP() the array of KSPs IS to be freed by the user 1221 (not the KSP just the array that contains them). 1222 1223 You must call KSPSetUp() before calling PCFieldSplitGetSubKSP(). 1224 1225 Level: advanced 1226 1227 .seealso: PCFIELDSPLIT 1228 @*/ 1229 PetscErrorCode PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[]) 1230 { 1231 PetscErrorCode ierr; 1232 1233 PetscFunctionBegin; 1234 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1235 if (n) PetscValidIntPointer(n,2); 1236 ierr = PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));CHKERRQ(ierr); 1237 PetscFunctionReturn(0); 1238 } 1239 1240 #undef __FUNCT__ 1241 #define __FUNCT__ "PCFieldSplitSchurPrecondition" 1242 /*@ 1243 PCFieldSplitSchurPrecondition - Indicates if the Schur complement is preconditioned by a preconditioner constructed by the 1244 A11 matrix. Otherwise no preconditioner is used. 1245 1246 Collective on PC 1247 1248 Input Parameters: 1249 + pc - the preconditioner context 1250 . ptype - which matrix to use for preconditioning the Schur complement, PC_FIELDSPLIT_SCHUR_PRE_DIAG (diag) is default 1251 - userpre - matrix to use for preconditioning, or PETSC_NULL 1252 1253 Options Database: 1254 . -pc_fieldsplit_schur_precondition <self,user,diag> default is diag 1255 1256 Notes: 1257 $ If ptype is 1258 $ user then the preconditioner for the Schur complement is generated by the provided matrix (pre argument 1259 $ to this function). 1260 $ diag then the preconditioner for the Schur complement is generated by the block diagonal part of the original 1261 $ matrix associated with the Schur complement (i.e. A11) 1262 $ self the preconditioner for the Schur complement is generated from the Schur complement matrix itself: 1263 $ The only preconditioner that currently works directly with the Schur complement matrix object is the PCLSC 1264 $ preconditioner 1265 1266 When solving a saddle point problem, where the A11 block is identically zero, using diag as the ptype only makes sense 1267 with the additional option -fieldsplit_1_pc_type none. Usually for saddle point problems one would use a ptype of self and 1268 -fieldsplit_1_pc_type lsc which uses the least squares commutator compute a preconditioner for the Schur complement. 1269 1270 Developer Notes: This is a terrible name, gives no good indication of what the function does and should also have Set in 1271 the name since it sets a proceedure to use. 1272 1273 Level: intermediate 1274 1275 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCLSC 1276 1277 @*/ 1278 PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1279 { 1280 PetscErrorCode ierr; 1281 1282 PetscFunctionBegin; 1283 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1284 ierr = PetscTryMethod(pc,"PCFieldSplitSchurPrecondition_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));CHKERRQ(ierr); 1285 PetscFunctionReturn(0); 1286 } 1287 1288 EXTERN_C_BEGIN 1289 #undef __FUNCT__ 1290 #define __FUNCT__ "PCFieldSplitSchurPrecondition_FieldSplit" 1291 PetscErrorCode PCFieldSplitSchurPrecondition_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1292 { 1293 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1294 PetscErrorCode ierr; 1295 1296 PetscFunctionBegin; 1297 jac->schurpre = ptype; 1298 if (pre) { 1299 ierr = MatDestroy(&jac->schur_user);CHKERRQ(ierr); 1300 jac->schur_user = pre; 1301 ierr = PetscObjectReference((PetscObject)jac->schur_user);CHKERRQ(ierr); 1302 } 1303 PetscFunctionReturn(0); 1304 } 1305 EXTERN_C_END 1306 1307 #undef __FUNCT__ 1308 #define __FUNCT__ "PCFieldSplitGetSchurBlocks" 1309 /*@C 1310 PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement 1311 1312 Collective on KSP 1313 1314 Input Parameter: 1315 . pc - the preconditioner context 1316 1317 Output Parameters: 1318 + A00 - the (0,0) block 1319 . A01 - the (0,1) block 1320 . A10 - the (1,0) block 1321 - A11 - the (1,1) block 1322 1323 Level: advanced 1324 1325 .seealso: PCFIELDSPLIT 1326 @*/ 1327 PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11) 1328 { 1329 PC_FieldSplit *jac = (PC_FieldSplit *) pc->data; 1330 1331 PetscFunctionBegin; 1332 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1333 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach."); 1334 if (A00) *A00 = jac->pmat[0]; 1335 if (A01) *A01 = jac->B; 1336 if (A10) *A10 = jac->C; 1337 if (A11) *A11 = jac->pmat[1]; 1338 PetscFunctionReturn(0); 1339 } 1340 1341 EXTERN_C_BEGIN 1342 #undef __FUNCT__ 1343 #define __FUNCT__ "PCFieldSplitSetType_FieldSplit" 1344 PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type) 1345 { 1346 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1347 PetscErrorCode ierr; 1348 1349 PetscFunctionBegin; 1350 jac->type = type; 1351 if (type == PC_COMPOSITE_SCHUR) { 1352 pc->ops->apply = PCApply_FieldSplit_Schur; 1353 pc->ops->view = PCView_FieldSplit_Schur; 1354 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit_Schur",PCFieldSplitGetSubKSP_FieldSplit_Schur);CHKERRQ(ierr); 1355 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSchurPrecondition_C","PCFieldSplitSchurPrecondition_FieldSplit",PCFieldSplitSchurPrecondition_FieldSplit);CHKERRQ(ierr); 1356 1357 } else { 1358 pc->ops->apply = PCApply_FieldSplit; 1359 pc->ops->view = PCView_FieldSplit; 1360 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit",PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 1361 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSchurPrecondition_C","",0);CHKERRQ(ierr); 1362 } 1363 PetscFunctionReturn(0); 1364 } 1365 EXTERN_C_END 1366 1367 EXTERN_C_BEGIN 1368 #undef __FUNCT__ 1369 #define __FUNCT__ "PCFieldSplitSetBlockSize_FieldSplit" 1370 PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs) 1371 { 1372 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1373 1374 PetscFunctionBegin; 1375 if (bs < 1) SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs); 1376 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); 1377 jac->bs = bs; 1378 PetscFunctionReturn(0); 1379 } 1380 EXTERN_C_END 1381 1382 #undef __FUNCT__ 1383 #define __FUNCT__ "PCFieldSplitSetType" 1384 /*@ 1385 PCFieldSplitSetType - Sets the type of fieldsplit preconditioner. 1386 1387 Collective on PC 1388 1389 Input Parameter: 1390 . pc - the preconditioner context 1391 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 1392 1393 Options Database Key: 1394 . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type 1395 1396 Level: Developer 1397 1398 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 1399 1400 .seealso: PCCompositeSetType() 1401 1402 @*/ 1403 PetscErrorCode PCFieldSplitSetType(PC pc,PCCompositeType type) 1404 { 1405 PetscErrorCode ierr; 1406 1407 PetscFunctionBegin; 1408 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1409 ierr = PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));CHKERRQ(ierr); 1410 PetscFunctionReturn(0); 1411 } 1412 1413 /* -------------------------------------------------------------------------------------*/ 1414 /*MC 1415 PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual 1416 fields or groups of fields. See the users manual section "Solving Block Matrices" for more details. 1417 1418 To set options on the solvers for each block append -fieldsplit_ to all the PC 1419 options database keys. For example, -fieldsplit_pc_type ilu -fieldsplit_pc_factor_levels 1 1420 1421 To set the options on the solvers separate for each block call PCFieldSplitGetSubKSP() 1422 and set the options directly on the resulting KSP object 1423 1424 Level: intermediate 1425 1426 Options Database Keys: 1427 + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split 1428 . -pc_fieldsplit_default - automatically add any fields to additional splits that have not 1429 been supplied explicitly by -pc_fieldsplit_%d_fields 1430 . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields) 1431 . -pc_fieldsplit_type <additive,multiplicative,schur,symmetric_multiplicative> 1432 . -pc_fieldsplit_schur_precondition <true,false> default is true 1433 . -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero or negative diagonal and uses Schur complement with no preconditioner as the solver 1434 1435 - Options prefix for inner solvers when using Schur complement preconditioner are -fieldsplit_0_ and -fieldsplit_1_ 1436 for all other solvers they are -fieldsplit_%d_ for the dth field, use -fieldsplit_ for all fields 1437 1438 Notes: 1439 Use PCFieldSplitSetFields() to set fields defined by "strided" entries and PCFieldSplitSetIS() 1440 to define a field by an arbitrary collection of entries. 1441 1442 If no fields are set the default is used. The fields are defined by entries strided by bs, 1443 beginning at 0 then 1, etc to bs-1. The block size can be set with PCFieldSplitSetBlockSize(), 1444 if this is not called the block size defaults to the blocksize of the second matrix passed 1445 to KSPSetOperators()/PCSetOperators(). 1446 1447 $ For the Schur complement preconditioner if J = ( A00 A01 ) 1448 $ ( A10 A11 ) 1449 $ the preconditioner using full factorization is 1450 $ ( I -A10 ksp(A00) ) ( inv(A00) 0 ) ( I 0 ) 1451 $ ( 0 I ) ( 0 ksp(S) ) ( -A10 ksp(A00) I ) 1452 where the action of inv(A00) is applied using the KSP solver with prefix -fieldsplit_0_. The action of 1453 ksp(S) is computed using the KSP solver with prefix -fieldsplit_splitname_ (where splitname was given 1454 in providing the SECOND split or 1 if not give). For PCFieldSplitGetKSP() when field number is 0, 1455 it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default 1456 A11 is used to construct a preconditioner for S, use PCFieldSplitSchurPrecondition() to turn on or off this 1457 option. You can use the preconditioner PCLSC to precondition the Schur complement with -fieldsplit_1_pc_type lsc. The 1458 factorization type is set using -pc_fieldsplit_schur_factorization_type <diag, lower, upper, full>. The full is shown above, 1459 but diag gives 1460 $ ( inv(A00) 0 ) 1461 $ ( 0 -ksp(S) ) 1462 so that the preconditioner is positive definite. The lower factorization is the inverse of 1463 $ ( A00 0 ) 1464 $ ( A10 S ) 1465 where the inverses of A00 and S are applied using KSPs. The upper factorization is the inverse of 1466 $ ( A00 A01 ) 1467 $ ( 0 S ) 1468 where again the inverses of A00 and S are applied using KSPs. 1469 1470 If only one set of indices (one IS) is provided with PCFieldSplitSetIS() then the complement of that IS 1471 is used automatically for a second block. 1472 1473 The fieldsplit preconditioner cannot currently be used with the BAIJ or SBAIJ data formats if the blocksize is larger than 1. 1474 Generally it should be used with the AIJ format. 1475 1476 The forms of these preconditioners are closely related if not identical to forms derived as "Distributive Iterations", see, 1477 for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling. Note that one can also use PCFIELDSPLIT 1478 inside a smoother resulting in "Distributive Smoothers". 1479 1480 Concepts: physics based preconditioners, block preconditioners 1481 1482 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners, PCLSC, 1483 PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSchurPrecondition() 1484 M*/ 1485 1486 EXTERN_C_BEGIN 1487 #undef __FUNCT__ 1488 #define __FUNCT__ "PCCreate_FieldSplit" 1489 PetscErrorCode PCCreate_FieldSplit(PC pc) 1490 { 1491 PetscErrorCode ierr; 1492 PC_FieldSplit *jac; 1493 1494 PetscFunctionBegin; 1495 ierr = PetscNewLog(pc,PC_FieldSplit,&jac);CHKERRQ(ierr); 1496 jac->bs = -1; 1497 jac->nsplits = 0; 1498 jac->type = PC_COMPOSITE_MULTIPLICATIVE; 1499 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */ 1500 jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACTORIZATION_FULL; 1501 1502 pc->data = (void*)jac; 1503 1504 pc->ops->apply = PCApply_FieldSplit; 1505 pc->ops->applytranspose = PCApplyTranspose_FieldSplit; 1506 pc->ops->setup = PCSetUp_FieldSplit; 1507 pc->ops->reset = PCReset_FieldSplit; 1508 pc->ops->destroy = PCDestroy_FieldSplit; 1509 pc->ops->setfromoptions = PCSetFromOptions_FieldSplit; 1510 pc->ops->view = PCView_FieldSplit; 1511 pc->ops->applyrichardson = 0; 1512 1513 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit", 1514 PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 1515 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetFields_C","PCFieldSplitSetFields_FieldSplit", 1516 PCFieldSplitSetFields_FieldSplit);CHKERRQ(ierr); 1517 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetIS_C","PCFieldSplitSetIS_FieldSplit", 1518 PCFieldSplitSetIS_FieldSplit);CHKERRQ(ierr); 1519 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetType_C","PCFieldSplitSetType_FieldSplit", 1520 PCFieldSplitSetType_FieldSplit);CHKERRQ(ierr); 1521 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetBlockSize_C","PCFieldSplitSetBlockSize_FieldSplit", 1522 PCFieldSplitSetBlockSize_FieldSplit);CHKERRQ(ierr); 1523 PetscFunctionReturn(0); 1524 } 1525 EXTERN_C_END 1526 1527 1528 1529