1 #define PETSCKSP_DLL 2 3 /* 4 5 */ 6 #include "private/pcimpl.h" /*I "petscpc.h" I*/ 7 8 typedef struct _PC_FieldSplitLink *PC_FieldSplitLink; 9 struct _PC_FieldSplitLink { 10 KSP ksp; 11 Vec x,y; 12 PetscInt nfields; 13 PetscInt *fields; 14 VecScatter sctx; 15 IS is,cis; 16 PetscInt csize; 17 PC_FieldSplitLink next,previous; 18 }; 19 20 typedef struct { 21 PCCompositeType type; 22 PetscTruth defaultsplit; 23 PetscInt bs; 24 PetscInt nsplits; 25 Vec *x,*y,w1,w2; 26 Mat *pmat; 27 Mat *Afield; /* the rows of the matrix associated with each field */ 28 PetscTruth issetup; 29 Mat B,C,schur; /* only used when Schur complement preconditioning is used */ 30 KSP kspschur; 31 PC_FieldSplitLink head; 32 } PC_FieldSplit; 33 34 /* 35 Notes: there is no particular reason that pmat, x, and y are stored as arrays in PC_FieldSplit instead of 36 inside PC_FieldSplitLink, just historical. If you want to be able to add new fields after already using the 37 PC you could change this. 38 */ 39 #undef __FUNCT__ 40 #define __FUNCT__ "PCView_FieldSplit" 41 static PetscErrorCode PCView_FieldSplit(PC pc,PetscViewer viewer) 42 { 43 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 44 PetscErrorCode ierr; 45 PetscTruth iascii; 46 PetscInt i,j; 47 PC_FieldSplitLink ilink = jac->head; 48 49 PetscFunctionBegin; 50 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 51 if (iascii) { 52 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D, blocksize = %D\n",PCCompositeTypes[jac->type],jac->nsplits,jac->bs);CHKERRQ(ierr); 53 ierr = PetscViewerASCIIPrintf(viewer," Solver info for each split is in the following KSP objects:\n");CHKERRQ(ierr); 54 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 55 for (i=0; i<jac->nsplits; i++) { 56 if (ilink->fields) { 57 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);CHKERRQ(ierr); 58 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 59 for (j=0; j<ilink->nfields; j++) { 60 if (j > 0) { 61 ierr = PetscViewerASCIIPrintf(viewer,",");CHKERRQ(ierr); 62 } 63 ierr = PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);CHKERRQ(ierr); 64 } 65 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 66 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 67 } else { 68 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i);CHKERRQ(ierr); 69 } 70 ierr = KSPView(ilink->ksp,viewer);CHKERRQ(ierr); 71 ilink = ilink->next; 72 } 73 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 74 } else { 75 SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for PCFieldSplit",((PetscObject)viewer)->type_name); 76 } 77 PetscFunctionReturn(0); 78 } 79 80 #undef __FUNCT__ 81 #define __FUNCT__ "PCView_FieldSplit_Schur" 82 static PetscErrorCode PCView_FieldSplit_Schur(PC pc,PetscViewer viewer) 83 { 84 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 85 PetscErrorCode ierr; 86 PetscTruth iascii; 87 PetscInt i,j; 88 PC_FieldSplitLink ilink = jac->head; 89 KSP ksp; 90 91 PetscFunctionBegin; 92 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 93 if (iascii) { 94 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit with Schur preconditioner, blocksize = %D\n",jac->bs);CHKERRQ(ierr); 95 ierr = PetscViewerASCIIPrintf(viewer," Split info:\n");CHKERRQ(ierr); 96 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 97 for (i=0; i<jac->nsplits; i++) { 98 if (ilink->fields) { 99 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);CHKERRQ(ierr); 100 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 101 for (j=0; j<ilink->nfields; j++) { 102 if (j > 0) { 103 ierr = PetscViewerASCIIPrintf(viewer,",");CHKERRQ(ierr); 104 } 105 ierr = PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);CHKERRQ(ierr); 106 } 107 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 108 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 109 } else { 110 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i);CHKERRQ(ierr); 111 } 112 ilink = ilink->next; 113 } 114 ierr = PetscViewerASCIIPrintf(viewer,"KSP solver for A block \n");CHKERRQ(ierr); 115 ierr = MatSchurComplementGetKSP(jac->schur,&ksp);CHKERRQ(ierr); 116 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 117 ierr = KSPView(ksp,viewer);CHKERRQ(ierr); 118 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 119 ierr = PetscViewerASCIIPrintf(viewer,"KSP solver for S = D - C inv(A) B \n");CHKERRQ(ierr); 120 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 121 ierr = KSPView(jac->kspschur,viewer);CHKERRQ(ierr); 122 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 123 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 124 } else { 125 SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for PCFieldSplit",((PetscObject)viewer)->type_name); 126 } 127 PetscFunctionReturn(0); 128 } 129 130 #undef __FUNCT__ 131 #define __FUNCT__ "PCFieldSplitSetDefaults" 132 static PetscErrorCode PCFieldSplitSetDefaults(PC pc) 133 { 134 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 135 PetscErrorCode ierr; 136 PC_FieldSplitLink ilink = jac->head; 137 PetscInt i = 0,*ifields,nfields; 138 PetscTruth flg = PETSC_FALSE,*fields,flg2; 139 char optionname[128]; 140 141 PetscFunctionBegin; 142 if (!ilink) { 143 144 if (jac->bs <= 0) { 145 if (pc->pmat) { 146 ierr = MatGetBlockSize(pc->pmat,&jac->bs);CHKERRQ(ierr); 147 } else { 148 jac->bs = 1; 149 } 150 } 151 152 ierr = PetscOptionsGetTruth(((PetscObject)pc)->prefix,"-pc_fieldsplit_default",&flg,PETSC_NULL);CHKERRQ(ierr); 153 if (!flg) { 154 /* Allow user to set fields from command line, if bs was known at the time of PCSetFromOptions_FieldSplit() 155 then it is set there. This is not ideal because we should only have options set in XXSetFromOptions(). */ 156 flg = PETSC_TRUE; /* switched off automatically if user sets fields manually here */ 157 ierr = PetscMalloc(jac->bs*sizeof(PetscInt),&ifields);CHKERRQ(ierr); 158 while (PETSC_TRUE) { 159 sprintf(optionname,"-pc_fieldsplit_%d_fields",(int)i++); 160 nfields = jac->bs; 161 ierr = PetscOptionsGetIntArray(((PetscObject)pc)->prefix,optionname,ifields,&nfields,&flg2);CHKERRQ(ierr); 162 if (!flg2) break; 163 if (!nfields) SETERRQ(PETSC_ERR_USER,"Cannot list zero fields"); 164 flg = PETSC_FALSE; 165 ierr = PCFieldSplitSetFields(pc,nfields,ifields);CHKERRQ(ierr); 166 } 167 ierr = PetscFree(ifields);CHKERRQ(ierr); 168 } 169 170 if (flg) { 171 ierr = PetscInfo(pc,"Using default splitting of fields\n");CHKERRQ(ierr); 172 ierr = PetscMalloc(jac->bs*sizeof(PetscTruth),&fields);CHKERRQ(ierr); 173 ierr = PetscMemzero(fields,jac->bs*sizeof(PetscTruth));CHKERRQ(ierr); 174 while (ilink) { 175 for (i=0; i<ilink->nfields; i++) { 176 fields[ilink->fields[i]] = PETSC_TRUE; 177 } 178 ilink = ilink->next; 179 } 180 jac->defaultsplit = PETSC_TRUE; 181 for (i=0; i<jac->bs; i++) { 182 if (!fields[i]) { 183 ierr = PCFieldSplitSetFields(pc,1,&i);CHKERRQ(ierr); 184 } else { 185 jac->defaultsplit = PETSC_FALSE; 186 } 187 } 188 ierr = PetscFree(fields);CHKERRQ(ierr); 189 } 190 } 191 PetscFunctionReturn(0); 192 } 193 194 195 #undef __FUNCT__ 196 #define __FUNCT__ "PCSetUp_FieldSplit" 197 static PetscErrorCode PCSetUp_FieldSplit(PC pc) 198 { 199 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 200 PetscErrorCode ierr; 201 PC_FieldSplitLink ilink; 202 PetscInt i,nsplit,ccsize; 203 MatStructure flag = pc->flag; 204 PetscTruth sorted,getsub; 205 206 PetscFunctionBegin; 207 ierr = PCFieldSplitSetDefaults(pc);CHKERRQ(ierr); 208 nsplit = jac->nsplits; 209 ilink = jac->head; 210 211 /* get the matrices for each split */ 212 if (!jac->issetup) { 213 PetscInt rstart,rend,nslots,bs; 214 215 jac->issetup = PETSC_TRUE; 216 217 /* This is done here instead of in PCFieldSplitSetFields() because may not have matrix at that point */ 218 bs = jac->bs; 219 ierr = MatGetOwnershipRange(pc->pmat,&rstart,&rend);CHKERRQ(ierr); 220 ierr = MatGetLocalSize(pc->pmat,PETSC_NULL,&ccsize);CHKERRQ(ierr); 221 nslots = (rend - rstart)/bs; 222 for (i=0; i<nsplit; i++) { 223 if (jac->defaultsplit) { 224 ierr = ISCreateStride(((PetscObject)pc)->comm,nslots,rstart+i,nsplit,&ilink->is);CHKERRQ(ierr); 225 } else if (!ilink->is) { 226 if (ilink->nfields > 1) { 227 PetscInt *ii,j,k,nfields = ilink->nfields,*fields = ilink->fields; 228 ierr = PetscMalloc(ilink->nfields*nslots*sizeof(PetscInt),&ii);CHKERRQ(ierr); 229 for (j=0; j<nslots; j++) { 230 for (k=0; k<nfields; k++) { 231 ii[nfields*j + k] = rstart + bs*j + fields[k]; 232 } 233 } 234 ierr = ISCreateGeneral(((PetscObject)pc)->comm,nslots*nfields,ii,&ilink->is);CHKERRQ(ierr); 235 ierr = PetscFree(ii);CHKERRQ(ierr); 236 } else { 237 ierr = ISCreateStride(((PetscObject)pc)->comm,nslots,rstart+ilink->fields[0],bs,&ilink->is);CHKERRQ(ierr); 238 } 239 } 240 ierr = ISGetLocalSize(ilink->is,&ilink->csize);CHKERRQ(ierr); 241 ierr = ISSorted(ilink->is,&sorted);CHKERRQ(ierr); 242 if (!sorted) SETERRQ(PETSC_ERR_USER,"Fields must be sorted when creating split"); 243 ierr = ISAllGather(ilink->is,&ilink->cis);CHKERRQ(ierr); 244 ilink = ilink->next; 245 } 246 } 247 248 ilink = jac->head; 249 if (!jac->pmat) { 250 ierr = PetscMalloc(nsplit*sizeof(Mat),&jac->pmat);CHKERRQ(ierr); 251 for (i=0; i<nsplit; i++) { 252 ierr = MatGetSubMatrix(pc->pmat,ilink->is,ilink->cis,ilink->csize,MAT_INITIAL_MATRIX,&jac->pmat[i]);CHKERRQ(ierr); 253 ilink = ilink->next; 254 } 255 } else { 256 for (i=0; i<nsplit; i++) { 257 ierr = MatGetSubMatrix(pc->pmat,ilink->is,ilink->cis,ilink->csize,MAT_REUSE_MATRIX,&jac->pmat[i]);CHKERRQ(ierr); 258 ilink = ilink->next; 259 } 260 } 261 262 /* extract the rows of the matrix associated with each field: used for efficient computation of residual inside algorithm */ 263 ierr = MatHasOperation(pc->mat,MATOP_GET_SUBMATRIX,&getsub);CHKERRQ(ierr); 264 if (getsub && jac->type != PC_COMPOSITE_ADDITIVE && jac->type != PC_COMPOSITE_SCHUR) { 265 ilink = jac->head; 266 if (!jac->Afield) { 267 ierr = PetscMalloc(nsplit*sizeof(Mat),&jac->Afield);CHKERRQ(ierr); 268 for (i=0; i<nsplit; i++) { 269 ierr = MatGetSubMatrix(pc->mat,ilink->is,PETSC_NULL,PETSC_DECIDE,MAT_INITIAL_MATRIX,&jac->Afield[i]);CHKERRQ(ierr); 270 ilink = ilink->next; 271 } 272 } else { 273 for (i=0; i<nsplit; i++) { 274 ierr = MatGetSubMatrix(pc->mat,ilink->is,PETSC_NULL,PETSC_DECIDE,MAT_REUSE_MATRIX,&jac->Afield[i]);CHKERRQ(ierr); 275 ilink = ilink->next; 276 } 277 } 278 } 279 280 if (jac->type == PC_COMPOSITE_SCHUR) { 281 IS ccis; 282 PetscInt N; 283 if (nsplit != 2) SETERRQ(PETSC_ERR_ARG_INCOMP,"To use Schur complement preconditioner you must have exactly 2 fields"); 284 285 /* need to handle case when one is resetting up the preconditioner */ 286 if (jac->schur) { 287 ierr = MatGetSize(pc->mat,PETSC_NULL,&N);CHKERRQ(ierr); 288 ilink = jac->head; 289 ierr = ISComplement(ilink->cis,N,&ccis);CHKERRQ(ierr); 290 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,PETSC_DECIDE,MAT_REUSE_MATRIX,&jac->B);CHKERRQ(ierr); 291 ierr = ISDestroy(ccis);CHKERRQ(ierr); 292 ilink = ilink->next; 293 ierr = ISComplement(ilink->cis,N,&ccis);CHKERRQ(ierr); 294 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,PETSC_DECIDE,MAT_REUSE_MATRIX,&jac->C);CHKERRQ(ierr); 295 ierr = ISDestroy(ccis);CHKERRQ(ierr); 296 ierr = MatSchurComplementUpdate(jac->schur,jac->pmat[0],jac->B,jac->C,jac->pmat[1],pc->flag);CHKERRQ(ierr); 297 ierr = KSPSetOperators(jac->kspschur,jac->schur,jac->schur,pc->flag);CHKERRQ(ierr); 298 299 } else { 300 301 /* extract the B and C matrices */ 302 ierr = MatGetSize(pc->mat,PETSC_NULL,&N);CHKERRQ(ierr); 303 ilink = jac->head; 304 ierr = ISComplement(ilink->cis,N,&ccis);CHKERRQ(ierr); 305 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,PETSC_DECIDE,MAT_INITIAL_MATRIX,&jac->B);CHKERRQ(ierr); 306 ierr = ISDestroy(ccis);CHKERRQ(ierr); 307 ilink = ilink->next; 308 ierr = ISComplement(ilink->cis,N,&ccis);CHKERRQ(ierr); 309 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,PETSC_DECIDE,MAT_INITIAL_MATRIX,&jac->C);CHKERRQ(ierr); 310 ierr = ISDestroy(ccis);CHKERRQ(ierr); 311 ierr = MatCreateSchurComplement(jac->pmat[0],jac->B,jac->C,jac->pmat[1],&jac->schur);CHKERRQ(ierr); 312 ierr = MatSetFromOptions(jac->schur);CHKERRQ(ierr); 313 314 ierr = KSPCreate(((PetscObject)pc)->comm,&jac->kspschur);CHKERRQ(ierr); 315 ierr = KSPSetOperators(jac->kspschur,jac->schur,jac->schur,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 316 ierr = KSPSetOptionsPrefix(jac->kspschur,((PetscObject)pc)->prefix);CHKERRQ(ierr); 317 ierr = KSPAppendOptionsPrefix(jac->kspschur,"schur_");CHKERRQ(ierr); 318 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 319 ierr = KSPSetFromOptions(jac->kspschur);CHKERRQ(ierr); 320 321 ierr = PetscMalloc2(2,Vec,&jac->x,2,Vec,&jac->y);CHKERRQ(ierr); 322 ierr = MatGetVecs(jac->pmat[0],&jac->x[0],&jac->y[0]);CHKERRQ(ierr); 323 ierr = MatGetVecs(jac->pmat[1],&jac->x[1],&jac->y[1]);CHKERRQ(ierr); 324 ilink = jac->head; 325 ilink->x = jac->x[0]; ilink->y = jac->y[0]; 326 ilink = ilink->next; 327 ilink->x = jac->x[1]; ilink->y = jac->y[1]; 328 } 329 } else { 330 /* set up the individual PCs */ 331 i = 0; 332 ilink = jac->head; 333 while (ilink) { 334 ierr = KSPSetOperators(ilink->ksp,jac->pmat[i],jac->pmat[i],flag);CHKERRQ(ierr); 335 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 336 ierr = KSPSetFromOptions(ilink->ksp);CHKERRQ(ierr); 337 ierr = KSPSetUp(ilink->ksp);CHKERRQ(ierr); 338 i++; 339 ilink = ilink->next; 340 } 341 342 /* create work vectors for each split */ 343 if (!jac->x) { 344 Vec xtmp; 345 ierr = PetscMalloc2(nsplit,Vec,&jac->x,nsplit,Vec,&jac->y);CHKERRQ(ierr); 346 ilink = jac->head; 347 for (i=0; i<nsplit; i++) { 348 Vec *vl,*vr; 349 350 ierr = KSPGetVecs(ilink->ksp,1,&vr,1,&vl);CHKERRQ(ierr); 351 ilink->x = *vr; 352 ilink->y = *vl; 353 ierr = PetscFree(vr);CHKERRQ(ierr); 354 ierr = PetscFree(vl);CHKERRQ(ierr); 355 jac->x[i] = ilink->x; 356 jac->y[i] = ilink->y; 357 ilink = ilink->next; 358 } 359 } 360 } 361 362 363 if (1) { 364 Vec xtmp; 365 366 /* compute scatter contexts needed by multiplicative versions and non-default splits */ 367 368 ilink = jac->head; 369 ierr = MatGetVecs(pc->pmat,&xtmp,PETSC_NULL);CHKERRQ(ierr); 370 for (i=0; i<nsplit; i++) { 371 ierr = VecScatterCreate(xtmp,ilink->is,jac->x[i],PETSC_NULL,&ilink->sctx);CHKERRQ(ierr); 372 ilink = ilink->next; 373 } 374 ierr = VecDestroy(xtmp);CHKERRQ(ierr); 375 } 376 PetscFunctionReturn(0); 377 } 378 379 #define FieldSplitSplitSolveAdd(ilink,xx,yy) \ 380 (VecScatterBegin(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 381 VecScatterEnd(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 382 KSPSolve(ilink->ksp,ilink->x,ilink->y) || \ 383 VecScatterBegin(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE) || \ 384 VecScatterEnd(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE)) 385 386 #undef __FUNCT__ 387 #define __FUNCT__ "PCApply_FieldSplit_Schur" 388 static PetscErrorCode PCApply_FieldSplit_Schur(PC pc,Vec x,Vec y) 389 { 390 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 391 PetscErrorCode ierr; 392 KSP ksp; 393 PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next; 394 395 PetscFunctionBegin; 396 ierr = MatSchurComplementGetKSP(jac->schur,&ksp);CHKERRQ(ierr); 397 398 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 399 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 400 ierr = KSPSolve(ksp,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 401 ierr = MatMult(jac->C,ilinkA->y,ilinkD->x);CHKERRQ(ierr); 402 ierr = VecScale(ilinkD->x,-1.0);CHKERRQ(ierr); 403 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 404 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 405 406 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 407 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 408 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 409 410 ierr = MatMult(jac->B,ilinkD->y,ilinkA->y);CHKERRQ(ierr); 411 ierr = VecAXPY(ilinkA->x,-1.0,ilinkA->y);CHKERRQ(ierr); 412 ierr = KSPSolve(ksp,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 413 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 414 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 415 416 PetscFunctionReturn(0); 417 } 418 419 #undef __FUNCT__ 420 #define __FUNCT__ "PCApply_FieldSplit" 421 static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y) 422 { 423 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 424 PetscErrorCode ierr; 425 PC_FieldSplitLink ilink = jac->head; 426 PetscInt bs,cnt; 427 428 PetscFunctionBegin; 429 CHKMEMQ; 430 ierr = VecGetBlockSize(x,&bs);CHKERRQ(ierr); 431 ierr = VecSetBlockSize(x,jac->bs);CHKERRQ(ierr); 432 ierr = VecSetBlockSize(y,jac->bs);CHKERRQ(ierr); 433 434 if (jac->type == PC_COMPOSITE_ADDITIVE) { 435 if (jac->defaultsplit) { 436 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 437 while (ilink) { 438 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 439 ilink = ilink->next; 440 } 441 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 442 } else { 443 ierr = VecSet(y,0.0);CHKERRQ(ierr); 444 while (ilink) { 445 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 446 ilink = ilink->next; 447 } 448 } 449 } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 450 if (!jac->w1) { 451 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 452 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 453 } 454 ierr = VecSet(y,0.0);CHKERRQ(ierr); 455 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 456 cnt = 1; 457 while (ilink->next) { 458 ilink = ilink->next; 459 if (jac->Afield) { 460 /* compute the residual only over the part of the vector needed */ 461 ierr = MatMult(jac->Afield[cnt++],y,ilink->x);CHKERRQ(ierr); 462 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 463 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 464 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 465 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 466 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 467 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 468 } else { 469 /* compute the residual over the entire vector */ 470 ierr = MatMult(pc->mat,y,jac->w1);CHKERRQ(ierr); 471 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 472 ierr = FieldSplitSplitSolveAdd(ilink,jac->w2,y);CHKERRQ(ierr); 473 } 474 } 475 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 476 cnt -= 2; 477 while (ilink->previous) { 478 ilink = ilink->previous; 479 if (jac->Afield) { 480 /* compute the residual only over the part of the vector needed */ 481 ierr = MatMult(jac->Afield[cnt--],y,ilink->x);CHKERRQ(ierr); 482 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 483 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 484 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 485 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 486 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 487 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 488 } else { 489 ierr = MatMult(pc->mat,y,jac->w1);CHKERRQ(ierr); 490 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 491 ierr = FieldSplitSplitSolveAdd(ilink,jac->w2,y);CHKERRQ(ierr); 492 } 493 } 494 } 495 } else SETERRQ1(PETSC_ERR_SUP,"Unsupported or unknown composition",(int) jac->type); 496 CHKMEMQ; 497 PetscFunctionReturn(0); 498 } 499 500 #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \ 501 (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 502 VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 503 KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) || \ 504 VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \ 505 VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE)) 506 507 #undef __FUNCT__ 508 #define __FUNCT__ "PCApply_FieldSplit" 509 static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y) 510 { 511 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 512 PetscErrorCode ierr; 513 PC_FieldSplitLink ilink = jac->head; 514 PetscInt bs; 515 516 PetscFunctionBegin; 517 CHKMEMQ; 518 ierr = VecGetBlockSize(x,&bs);CHKERRQ(ierr); 519 ierr = VecSetBlockSize(x,jac->bs);CHKERRQ(ierr); 520 ierr = VecSetBlockSize(y,jac->bs);CHKERRQ(ierr); 521 522 if (jac->type == PC_COMPOSITE_ADDITIVE) { 523 if (jac->defaultsplit) { 524 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 525 while (ilink) { 526 ierr = KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 527 ilink = ilink->next; 528 } 529 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 530 } else { 531 ierr = VecSet(y,0.0);CHKERRQ(ierr); 532 while (ilink) { 533 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 534 ilink = ilink->next; 535 } 536 } 537 } else { 538 if (!jac->w1) { 539 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 540 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 541 } 542 ierr = VecSet(y,0.0);CHKERRQ(ierr); 543 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 544 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 545 while (ilink->next) { 546 ilink = ilink->next; 547 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 548 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 549 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 550 } 551 while (ilink->previous) { 552 ilink = ilink->previous; 553 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 554 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 555 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 556 } 557 } else { 558 while (ilink->next) { /* get to last entry in linked list */ 559 ilink = ilink->next; 560 } 561 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 562 while (ilink->previous) { 563 ilink = ilink->previous; 564 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 565 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 566 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 567 } 568 } 569 } 570 CHKMEMQ; 571 PetscFunctionReturn(0); 572 } 573 574 #undef __FUNCT__ 575 #define __FUNCT__ "PCDestroy_FieldSplit" 576 static PetscErrorCode PCDestroy_FieldSplit(PC pc) 577 { 578 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 579 PetscErrorCode ierr; 580 PC_FieldSplitLink ilink = jac->head,next; 581 582 PetscFunctionBegin; 583 while (ilink) { 584 ierr = KSPDestroy(ilink->ksp);CHKERRQ(ierr); 585 if (ilink->x) {ierr = VecDestroy(ilink->x);CHKERRQ(ierr);} 586 if (ilink->y) {ierr = VecDestroy(ilink->y);CHKERRQ(ierr);} 587 if (ilink->sctx) {ierr = VecScatterDestroy(ilink->sctx);CHKERRQ(ierr);} 588 if (ilink->is) {ierr = ISDestroy(ilink->is);CHKERRQ(ierr);} 589 if (ilink->cis) {ierr = ISDestroy(ilink->cis);CHKERRQ(ierr);} 590 next = ilink->next; 591 ierr = PetscFree(ilink->fields);CHKERRQ(ierr); 592 ierr = PetscFree(ilink);CHKERRQ(ierr); 593 ilink = next; 594 } 595 ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr); 596 if (jac->pmat) {ierr = MatDestroyMatrices(jac->nsplits,&jac->pmat);CHKERRQ(ierr);} 597 if (jac->Afield) {ierr = MatDestroyMatrices(jac->nsplits,&jac->Afield);CHKERRQ(ierr);} 598 if (jac->w1) {ierr = VecDestroy(jac->w1);CHKERRQ(ierr);} 599 if (jac->w2) {ierr = VecDestroy(jac->w2);CHKERRQ(ierr);} 600 if (jac->schur) {ierr = MatDestroy(jac->schur);CHKERRQ(ierr);} 601 if (jac->B) {ierr = MatDestroy(jac->B);CHKERRQ(ierr);} 602 if (jac->C) {ierr = MatDestroy(jac->C);CHKERRQ(ierr);} 603 ierr = PetscFree(jac);CHKERRQ(ierr); 604 PetscFunctionReturn(0); 605 } 606 607 #undef __FUNCT__ 608 #define __FUNCT__ "PCSetFromOptions_FieldSplit" 609 static PetscErrorCode PCSetFromOptions_FieldSplit(PC pc) 610 { 611 PetscErrorCode ierr; 612 PetscInt i = 0,nfields,*fields,bs; 613 PetscTruth flg; 614 char optionname[128]; 615 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 616 PCCompositeType ctype; 617 618 PetscFunctionBegin; 619 ierr = PetscOptionsHead("FieldSplit options");CHKERRQ(ierr); 620 ierr = PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);CHKERRQ(ierr); 621 if (flg) { 622 ierr = PCFieldSplitSetBlockSize(pc,bs);CHKERRQ(ierr); 623 } 624 625 ierr = PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg);CHKERRQ(ierr); 626 if (flg) { 627 ierr = PCFieldSplitSetType(pc,ctype);CHKERRQ(ierr); 628 } 629 630 if (jac->bs > 0) { 631 /* only allow user to set fields from command line if bs is already known. 632 otherwise user can set them in PCFieldSplitSetDefaults() */ 633 ierr = PetscMalloc(jac->bs*sizeof(PetscInt),&fields);CHKERRQ(ierr); 634 while (PETSC_TRUE) { 635 sprintf(optionname,"-pc_fieldsplit_%d_fields",(int)i++); 636 nfields = jac->bs; 637 ierr = PetscOptionsIntArray(optionname,"Fields in this split","PCFieldSplitSetFields",fields,&nfields,&flg);CHKERRQ(ierr); 638 if (!flg) break; 639 if (!nfields) SETERRQ(PETSC_ERR_USER,"Cannot list zero fields"); 640 ierr = PCFieldSplitSetFields(pc,nfields,fields);CHKERRQ(ierr); 641 } 642 ierr = PetscFree(fields);CHKERRQ(ierr); 643 } 644 ierr = PetscOptionsTail();CHKERRQ(ierr); 645 PetscFunctionReturn(0); 646 } 647 648 /*------------------------------------------------------------------------------------*/ 649 650 EXTERN_C_BEGIN 651 #undef __FUNCT__ 652 #define __FUNCT__ "PCFieldSplitSetFields_FieldSplit" 653 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetFields_FieldSplit(PC pc,PetscInt n,PetscInt *fields) 654 { 655 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 656 PetscErrorCode ierr; 657 PC_FieldSplitLink ilink,next = jac->head; 658 char prefix[128]; 659 PetscInt i; 660 661 PetscFunctionBegin; 662 if (n <= 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Negative number of fields requested"); 663 for (i=0; i<n; i++) { 664 if (fields[i] >= jac->bs) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs); 665 if (fields[i] < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]); 666 } 667 ierr = PetscNew(struct _PC_FieldSplitLink,&ilink);CHKERRQ(ierr); 668 ierr = PetscMalloc(n*sizeof(PetscInt),&ilink->fields);CHKERRQ(ierr); 669 ierr = PetscMemcpy(ilink->fields,fields,n*sizeof(PetscInt));CHKERRQ(ierr); 670 ilink->nfields = n; 671 ilink->next = PETSC_NULL; 672 ierr = KSPCreate(((PetscObject)pc)->comm,&ilink->ksp);CHKERRQ(ierr); 673 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 674 675 if (((PetscObject)pc)->prefix) { 676 sprintf(prefix,"%sfieldsplit_%d_",((PetscObject)pc)->prefix,(int)jac->nsplits); 677 } else { 678 sprintf(prefix,"fieldsplit_%d_",(int)jac->nsplits); 679 } 680 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 681 682 if (!next) { 683 jac->head = ilink; 684 ilink->previous = PETSC_NULL; 685 } else { 686 while (next->next) { 687 next = next->next; 688 } 689 next->next = ilink; 690 ilink->previous = next; 691 } 692 jac->nsplits++; 693 PetscFunctionReturn(0); 694 } 695 EXTERN_C_END 696 697 698 EXTERN_C_BEGIN 699 #undef __FUNCT__ 700 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit" 701 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp) 702 { 703 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 704 PetscErrorCode ierr; 705 PetscInt cnt = 0; 706 PC_FieldSplitLink ilink = jac->head; 707 708 PetscFunctionBegin; 709 ierr = PetscMalloc(jac->nsplits*sizeof(KSP*),subksp);CHKERRQ(ierr); 710 while (ilink) { 711 (*subksp)[cnt++] = ilink->ksp; 712 ilink = ilink->next; 713 } 714 if (cnt != jac->nsplits) SETERRQ2(PETSC_ERR_PLIB,"Corrupt PCFIELDSPLIT object: number splits in linked list %D in object %D",cnt,jac->nsplits); 715 *n = jac->nsplits; 716 PetscFunctionReturn(0); 717 } 718 EXTERN_C_END 719 720 EXTERN_C_BEGIN 721 #undef __FUNCT__ 722 #define __FUNCT__ "PCFieldSplitSetIS_FieldSplit" 723 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetIS_FieldSplit(PC pc,IS is) 724 { 725 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 726 PetscErrorCode ierr; 727 PC_FieldSplitLink ilink, next = jac->head; 728 char prefix[128]; 729 730 PetscFunctionBegin; 731 ierr = PetscNew(struct _PC_FieldSplitLink,&ilink);CHKERRQ(ierr); 732 ilink->is = is; 733 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 734 ilink->next = PETSC_NULL; 735 ierr = KSPCreate(((PetscObject)pc)->comm,&ilink->ksp);CHKERRQ(ierr); 736 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 737 738 if (((PetscObject)pc)->prefix) { 739 sprintf(prefix,"%sfieldsplit_%d_",((PetscObject)pc)->prefix,(int)jac->nsplits); 740 } else { 741 sprintf(prefix,"fieldsplit_%d_",(int)jac->nsplits); 742 } 743 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 744 745 if (!next) { 746 jac->head = ilink; 747 ilink->previous = PETSC_NULL; 748 } else { 749 while (next->next) { 750 next = next->next; 751 } 752 next->next = ilink; 753 ilink->previous = next; 754 } 755 jac->nsplits++; 756 757 PetscFunctionReturn(0); 758 } 759 EXTERN_C_END 760 761 #undef __FUNCT__ 762 #define __FUNCT__ "PCFieldSplitSetFields" 763 /*@ 764 PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner 765 766 Collective on PC 767 768 Input Parameters: 769 + pc - the preconditioner context 770 . n - the number of fields in this split 771 . fields - the fields in this split 772 773 Level: intermediate 774 775 Notes: Use PCFieldSplitSetIS() to set a completely general set of indices as a field. 776 777 The PCFieldSplitSetFields() is for defining fields as a strided blocks. For example, if the block 778 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 779 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x.... 780 where the numbered entries indicate what is in the field. 781 782 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize(), PCFieldSplitSetIS() 783 784 @*/ 785 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetFields(PC pc,PetscInt n, PetscInt *fields) 786 { 787 PetscErrorCode ierr,(*f)(PC,PetscInt,PetscInt *); 788 789 PetscFunctionBegin; 790 PetscValidHeaderSpecific(pc,PC_COOKIE,1); 791 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCFieldSplitSetFields_C",(void (**)(void))&f);CHKERRQ(ierr); 792 if (f) { 793 ierr = (*f)(pc,n,fields);CHKERRQ(ierr); 794 } 795 PetscFunctionReturn(0); 796 } 797 798 #undef __FUNCT__ 799 #define __FUNCT__ "PCFieldSplitSetIS" 800 /*@ 801 PCFieldSplitSetIS - Sets the exact elements for field 802 803 Collective on PC 804 805 Input Parameters: 806 + pc - the preconditioner context 807 . is - the index set that defines the vector elements in this field 808 809 810 Notes: Use PCFieldSplitSetFields(), for fields defined by strided types. 811 812 Level: intermediate 813 814 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize() 815 816 @*/ 817 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetIS(PC pc,IS is) 818 { 819 PetscErrorCode ierr,(*f)(PC,IS); 820 821 PetscFunctionBegin; 822 PetscValidHeaderSpecific(pc,PC_COOKIE,1); 823 PetscValidHeaderSpecific(is,IS_COOKIE,1); 824 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCFieldSplitSetIS_C",(void (**)(void))&f);CHKERRQ(ierr); 825 if (f) { 826 ierr = (*f)(pc,is);CHKERRQ(ierr); 827 } 828 PetscFunctionReturn(0); 829 } 830 831 #undef __FUNCT__ 832 #define __FUNCT__ "PCFieldSplitSetBlockSize" 833 /*@ 834 PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the 835 fieldsplit preconditioner. If not set the matrix block size is used. 836 837 Collective on PC 838 839 Input Parameters: 840 + pc - the preconditioner context 841 - bs - the block size 842 843 Level: intermediate 844 845 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields() 846 847 @*/ 848 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetBlockSize(PC pc,PetscInt bs) 849 { 850 PetscErrorCode ierr,(*f)(PC,PetscInt); 851 852 PetscFunctionBegin; 853 PetscValidHeaderSpecific(pc,PC_COOKIE,1); 854 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",(void (**)(void))&f);CHKERRQ(ierr); 855 if (f) { 856 ierr = (*f)(pc,bs);CHKERRQ(ierr); 857 } 858 PetscFunctionReturn(0); 859 } 860 861 #undef __FUNCT__ 862 #define __FUNCT__ "PCFieldSplitGetSubKSP" 863 /*@C 864 PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits 865 866 Collective on KSP 867 868 Input Parameter: 869 . pc - the preconditioner context 870 871 Output Parameters: 872 + n - the number of split 873 - pc - the array of KSP contexts 874 875 Note: 876 After PCFieldSplitGetSubKSP() the array of KSPs IS to be freed by the user 877 (not the KSP just the array that contains them). 878 879 You must call KSPSetUp() before calling PCFieldSplitGetSubKSP(). 880 881 Level: advanced 882 883 .seealso: PCFIELDSPLIT 884 @*/ 885 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[]) 886 { 887 PetscErrorCode ierr,(*f)(PC,PetscInt*,KSP **); 888 889 PetscFunctionBegin; 890 PetscValidHeaderSpecific(pc,PC_COOKIE,1); 891 PetscValidIntPointer(n,2); 892 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",(void (**)(void))&f);CHKERRQ(ierr); 893 if (f) { 894 ierr = (*f)(pc,n,subksp);CHKERRQ(ierr); 895 } else { 896 SETERRQ(PETSC_ERR_ARG_WRONG,"Cannot get subksp for this type of PC"); 897 } 898 PetscFunctionReturn(0); 899 } 900 901 EXTERN_C_BEGIN 902 #undef __FUNCT__ 903 #define __FUNCT__ "PCFieldSplitSetType_FieldSplit" 904 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type) 905 { 906 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 907 908 PetscFunctionBegin; 909 jac->type = type; 910 if (type == PC_COMPOSITE_SCHUR) { 911 pc->ops->apply = PCApply_FieldSplit_Schur; 912 pc->ops->view = PCView_FieldSplit_Schur; 913 } else { 914 pc->ops->apply = PCApply_FieldSplit; 915 pc->ops->view = PCView_FieldSplit; 916 } 917 PetscFunctionReturn(0); 918 } 919 EXTERN_C_END 920 921 EXTERN_C_BEGIN 922 #undef __FUNCT__ 923 #define __FUNCT__ "PCFieldSplitSetBlockSize_FieldSplit" 924 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs) 925 { 926 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 927 928 PetscFunctionBegin; 929 if (bs < 1) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs); 930 if (jac->bs > 0 && jac->bs != bs) SETERRQ2(PETSC_ERR_ARG_WRONGSTATE,"Cannot change fieldsplit blocksize from %D to %D after it has been set",jac->bs,bs); 931 jac->bs = bs; 932 PetscFunctionReturn(0); 933 } 934 EXTERN_C_END 935 936 #undef __FUNCT__ 937 #define __FUNCT__ "PCFieldSplitSetType" 938 /*@ 939 PCFieldSplitSetType - Sets the type of fieldsplit preconditioner. 940 941 Collective on PC 942 943 Input Parameter: 944 . pc - the preconditioner context 945 . type - PC_COMPOSITE_ADDITIVE (default), PC_COMPOSITE_MULTIPLICATIVE, PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE 946 947 Options Database Key: 948 . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special> - Sets fieldsplit preconditioner type 949 950 Level: Developer 951 952 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 953 954 .seealso: PCCompositeSetType() 955 956 @*/ 957 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetType(PC pc,PCCompositeType type) 958 { 959 PetscErrorCode ierr,(*f)(PC,PCCompositeType); 960 961 PetscFunctionBegin; 962 PetscValidHeaderSpecific(pc,PC_COOKIE,1); 963 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCFieldSplitSetType_C",(void (**)(void))&f);CHKERRQ(ierr); 964 if (f) { 965 ierr = (*f)(pc,type);CHKERRQ(ierr); 966 } 967 PetscFunctionReturn(0); 968 } 969 970 /* -------------------------------------------------------------------------------------*/ 971 /*MC 972 PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual 973 fields or groups of fields 974 975 976 To set options on the solvers for each block append -sub_ to all the PC 977 options database keys. For example, -sub_pc_type ilu -sub_pc_factor_levels 1 978 979 To set the options on the solvers separate for each block call PCFieldSplitGetSubKSP() 980 and set the options directly on the resulting KSP object 981 982 Level: intermediate 983 984 Options Database Keys: 985 + -pc_splitfield_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split 986 . -pc_splitfield_default - automatically add any fields to additional splits that have not 987 been supplied explicitly by -pc_splitfield_%d_fields 988 . -pc_splitfield_block_size <bs> - size of block that defines fields (i.e. there are bs fields) 989 - -pc_splitfield_type <additive,multiplicative> 990 991 - Options prefix for inner solvers when using Schur complement preconditioner are -schurAblock_ and -schur, 992 for all other solvers they are -fieldsplit_%d_ for the dth field, use -fieldsplit_ for all fields 993 994 995 Notes: use PCFieldSplitSetFields() to set fields defined by "strided" entries and PCFieldSplitSetIS() 996 to define a field by an arbitrary collection of entries. 997 998 If no fields are set the default is used. The fields are defined by entries strided by bs, 999 beginning at 0 then 1, etc to bs-1. The block size can be set with PCFieldSplitSetBlockSize(), 1000 if this is not called the block size defaults to the blocksize of the second matrix passed 1001 to KSPSetOperators()/PCSetOperators(). 1002 1003 Currently for the multiplicative version, the updated residual needed for the next field 1004 solve is computed via a matrix vector product over the entire array. An optimization would be 1005 to update the residual only for the part of the right hand side associated with the next field 1006 solve. (This would involve more MatGetSubMatrix() calls or some other mechanism to compute the 1007 part of the matrix needed to just update part of the residual). 1008 1009 Concepts: physics based preconditioners 1010 1011 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, 1012 PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS() 1013 M*/ 1014 1015 EXTERN_C_BEGIN 1016 #undef __FUNCT__ 1017 #define __FUNCT__ "PCCreate_FieldSplit" 1018 PetscErrorCode PETSCKSP_DLLEXPORT PCCreate_FieldSplit(PC pc) 1019 { 1020 PetscErrorCode ierr; 1021 PC_FieldSplit *jac; 1022 1023 PetscFunctionBegin; 1024 ierr = PetscNewLog(pc,PC_FieldSplit,&jac);CHKERRQ(ierr); 1025 jac->bs = -1; 1026 jac->nsplits = 0; 1027 jac->type = PC_COMPOSITE_MULTIPLICATIVE; 1028 1029 pc->data = (void*)jac; 1030 1031 pc->ops->apply = PCApply_FieldSplit; 1032 pc->ops->applytranspose = PCApplyTranspose_FieldSplit; 1033 pc->ops->setup = PCSetUp_FieldSplit; 1034 pc->ops->destroy = PCDestroy_FieldSplit; 1035 pc->ops->setfromoptions = PCSetFromOptions_FieldSplit; 1036 pc->ops->view = PCView_FieldSplit; 1037 pc->ops->applyrichardson = 0; 1038 1039 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit", 1040 PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 1041 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetFields_C","PCFieldSplitSetFields_FieldSplit", 1042 PCFieldSplitSetFields_FieldSplit);CHKERRQ(ierr); 1043 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetIS_C","PCFieldSplitSetIS_FieldSplit", 1044 PCFieldSplitSetIS_FieldSplit);CHKERRQ(ierr); 1045 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetType_C","PCFieldSplitSetType_FieldSplit", 1046 PCFieldSplitSetType_FieldSplit);CHKERRQ(ierr); 1047 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetBlockSize_C","PCFieldSplitSetBlockSize_FieldSplit", 1048 PCFieldSplitSetBlockSize_FieldSplit);CHKERRQ(ierr); 1049 PetscFunctionReturn(0); 1050 } 1051 EXTERN_C_END 1052 1053 1054