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 KSP ksp; 301 302 /* extract the B and C matrices */ 303 ierr = MatGetSize(pc->mat,PETSC_NULL,&N);CHKERRQ(ierr); 304 ilink = jac->head; 305 ierr = ISComplement(ilink->cis,N,&ccis);CHKERRQ(ierr); 306 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,PETSC_DECIDE,MAT_INITIAL_MATRIX,&jac->B);CHKERRQ(ierr); 307 ierr = ISDestroy(ccis);CHKERRQ(ierr); 308 ilink = ilink->next; 309 ierr = ISComplement(ilink->cis,N,&ccis);CHKERRQ(ierr); 310 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,PETSC_DECIDE,MAT_INITIAL_MATRIX,&jac->C);CHKERRQ(ierr); 311 ierr = ISDestroy(ccis);CHKERRQ(ierr); 312 ierr = MatCreateSchurComplement(jac->pmat[0],jac->B,jac->C,jac->pmat[1],&jac->schur);CHKERRQ(ierr); 313 ierr = MatSchurComplementGetKSP(jac->schur,&ksp);CHKERRQ(ierr); 314 ierr = PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,1);CHKERRQ(ierr); 315 ierr = MatSetFromOptions(jac->schur);CHKERRQ(ierr); 316 317 ierr = KSPCreate(((PetscObject)pc)->comm,&jac->kspschur);CHKERRQ(ierr); 318 ierr = PetscObjectIncrementTabLevel((PetscObject)jac->kspschur,(PetscObject)pc,1);CHKERRQ(ierr); 319 ierr = KSPSetOperators(jac->kspschur,jac->schur,jac->schur,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 320 ierr = KSPSetOptionsPrefix(jac->kspschur,((PetscObject)pc)->prefix);CHKERRQ(ierr); 321 ierr = KSPAppendOptionsPrefix(jac->kspschur,"schur_");CHKERRQ(ierr); 322 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 323 ierr = KSPSetFromOptions(jac->kspschur);CHKERRQ(ierr); 324 325 ierr = PetscMalloc2(2,Vec,&jac->x,2,Vec,&jac->y);CHKERRQ(ierr); 326 ierr = MatGetVecs(jac->pmat[0],&jac->x[0],&jac->y[0]);CHKERRQ(ierr); 327 ierr = MatGetVecs(jac->pmat[1],&jac->x[1],&jac->y[1]);CHKERRQ(ierr); 328 ilink = jac->head; 329 ilink->x = jac->x[0]; ilink->y = jac->y[0]; 330 ilink = ilink->next; 331 ilink->x = jac->x[1]; ilink->y = jac->y[1]; 332 } 333 } else { 334 /* set up the individual PCs */ 335 i = 0; 336 ilink = jac->head; 337 while (ilink) { 338 ierr = KSPSetOperators(ilink->ksp,jac->pmat[i],jac->pmat[i],flag);CHKERRQ(ierr); 339 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 340 ierr = KSPSetFromOptions(ilink->ksp);CHKERRQ(ierr); 341 ierr = KSPSetUp(ilink->ksp);CHKERRQ(ierr); 342 i++; 343 ilink = ilink->next; 344 } 345 346 /* create work vectors for each split */ 347 if (!jac->x) { 348 Vec xtmp; 349 ierr = PetscMalloc2(nsplit,Vec,&jac->x,nsplit,Vec,&jac->y);CHKERRQ(ierr); 350 ilink = jac->head; 351 for (i=0; i<nsplit; i++) { 352 Vec *vl,*vr; 353 354 ierr = KSPGetVecs(ilink->ksp,1,&vr,1,&vl);CHKERRQ(ierr); 355 ilink->x = *vr; 356 ilink->y = *vl; 357 ierr = PetscFree(vr);CHKERRQ(ierr); 358 ierr = PetscFree(vl);CHKERRQ(ierr); 359 jac->x[i] = ilink->x; 360 jac->y[i] = ilink->y; 361 ilink = ilink->next; 362 } 363 } 364 } 365 366 367 if (1) { 368 Vec xtmp; 369 370 /* compute scatter contexts needed by multiplicative versions and non-default splits */ 371 372 ilink = jac->head; 373 ierr = MatGetVecs(pc->pmat,&xtmp,PETSC_NULL);CHKERRQ(ierr); 374 for (i=0; i<nsplit; i++) { 375 ierr = VecScatterCreate(xtmp,ilink->is,jac->x[i],PETSC_NULL,&ilink->sctx);CHKERRQ(ierr); 376 ilink = ilink->next; 377 } 378 ierr = VecDestroy(xtmp);CHKERRQ(ierr); 379 } 380 PetscFunctionReturn(0); 381 } 382 383 #define FieldSplitSplitSolveAdd(ilink,xx,yy) \ 384 (VecScatterBegin(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 385 VecScatterEnd(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 386 KSPSolve(ilink->ksp,ilink->x,ilink->y) || \ 387 VecScatterBegin(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE) || \ 388 VecScatterEnd(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE)) 389 390 #undef __FUNCT__ 391 #define __FUNCT__ "PCApply_FieldSplit_Schur" 392 static PetscErrorCode PCApply_FieldSplit_Schur(PC pc,Vec x,Vec y) 393 { 394 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 395 PetscErrorCode ierr; 396 KSP ksp; 397 PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next; 398 399 PetscFunctionBegin; 400 ierr = MatSchurComplementGetKSP(jac->schur,&ksp);CHKERRQ(ierr); 401 402 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 403 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 404 ierr = KSPSolve(ksp,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 405 ierr = MatMult(jac->C,ilinkA->y,ilinkD->x);CHKERRQ(ierr); 406 ierr = VecScale(ilinkD->x,-1.0);CHKERRQ(ierr); 407 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 408 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 409 410 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 411 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 412 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 413 414 ierr = MatMult(jac->B,ilinkD->y,ilinkA->y);CHKERRQ(ierr); 415 ierr = VecAXPY(ilinkA->x,-1.0,ilinkA->y);CHKERRQ(ierr); 416 ierr = KSPSolve(ksp,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 417 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 418 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 419 420 PetscFunctionReturn(0); 421 } 422 423 #undef __FUNCT__ 424 #define __FUNCT__ "PCApply_FieldSplit" 425 static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y) 426 { 427 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 428 PetscErrorCode ierr; 429 PC_FieldSplitLink ilink = jac->head; 430 PetscInt bs,cnt; 431 432 PetscFunctionBegin; 433 CHKMEMQ; 434 ierr = VecGetBlockSize(x,&bs);CHKERRQ(ierr); 435 ierr = VecSetBlockSize(x,jac->bs);CHKERRQ(ierr); 436 ierr = VecSetBlockSize(y,jac->bs);CHKERRQ(ierr); 437 438 if (jac->type == PC_COMPOSITE_ADDITIVE) { 439 if (jac->defaultsplit) { 440 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 441 while (ilink) { 442 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 443 ilink = ilink->next; 444 } 445 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 446 } else { 447 ierr = VecSet(y,0.0);CHKERRQ(ierr); 448 while (ilink) { 449 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 450 ilink = ilink->next; 451 } 452 } 453 } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 454 if (!jac->w1) { 455 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 456 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 457 } 458 ierr = VecSet(y,0.0);CHKERRQ(ierr); 459 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 460 cnt = 1; 461 while (ilink->next) { 462 ilink = ilink->next; 463 if (jac->Afield) { 464 /* compute the residual only over the part of the vector needed */ 465 ierr = MatMult(jac->Afield[cnt++],y,ilink->x);CHKERRQ(ierr); 466 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 467 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 468 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 469 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 470 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 471 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 472 } else { 473 /* compute the residual over the entire vector */ 474 ierr = MatMult(pc->mat,y,jac->w1);CHKERRQ(ierr); 475 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 476 ierr = FieldSplitSplitSolveAdd(ilink,jac->w2,y);CHKERRQ(ierr); 477 } 478 } 479 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 480 cnt -= 2; 481 while (ilink->previous) { 482 ilink = ilink->previous; 483 if (jac->Afield) { 484 /* compute the residual only over the part of the vector needed */ 485 ierr = MatMult(jac->Afield[cnt--],y,ilink->x);CHKERRQ(ierr); 486 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 487 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 488 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 489 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 490 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 491 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 492 } else { 493 ierr = MatMult(pc->mat,y,jac->w1);CHKERRQ(ierr); 494 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 495 ierr = FieldSplitSplitSolveAdd(ilink,jac->w2,y);CHKERRQ(ierr); 496 } 497 } 498 } 499 } else SETERRQ1(PETSC_ERR_SUP,"Unsupported or unknown composition",(int) jac->type); 500 CHKMEMQ; 501 PetscFunctionReturn(0); 502 } 503 504 #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \ 505 (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 506 VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 507 KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) || \ 508 VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \ 509 VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE)) 510 511 #undef __FUNCT__ 512 #define __FUNCT__ "PCApply_FieldSplit" 513 static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y) 514 { 515 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 516 PetscErrorCode ierr; 517 PC_FieldSplitLink ilink = jac->head; 518 PetscInt bs; 519 520 PetscFunctionBegin; 521 CHKMEMQ; 522 ierr = VecGetBlockSize(x,&bs);CHKERRQ(ierr); 523 ierr = VecSetBlockSize(x,jac->bs);CHKERRQ(ierr); 524 ierr = VecSetBlockSize(y,jac->bs);CHKERRQ(ierr); 525 526 if (jac->type == PC_COMPOSITE_ADDITIVE) { 527 if (jac->defaultsplit) { 528 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 529 while (ilink) { 530 ierr = KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 531 ilink = ilink->next; 532 } 533 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 534 } else { 535 ierr = VecSet(y,0.0);CHKERRQ(ierr); 536 while (ilink) { 537 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 538 ilink = ilink->next; 539 } 540 } 541 } else { 542 if (!jac->w1) { 543 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 544 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 545 } 546 ierr = VecSet(y,0.0);CHKERRQ(ierr); 547 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 548 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 549 while (ilink->next) { 550 ilink = ilink->next; 551 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 552 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 553 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 554 } 555 while (ilink->previous) { 556 ilink = ilink->previous; 557 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 558 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 559 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 560 } 561 } else { 562 while (ilink->next) { /* get to last entry in linked list */ 563 ilink = ilink->next; 564 } 565 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 566 while (ilink->previous) { 567 ilink = ilink->previous; 568 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 569 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 570 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 571 } 572 } 573 } 574 CHKMEMQ; 575 PetscFunctionReturn(0); 576 } 577 578 #undef __FUNCT__ 579 #define __FUNCT__ "PCDestroy_FieldSplit" 580 static PetscErrorCode PCDestroy_FieldSplit(PC pc) 581 { 582 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 583 PetscErrorCode ierr; 584 PC_FieldSplitLink ilink = jac->head,next; 585 586 PetscFunctionBegin; 587 while (ilink) { 588 ierr = KSPDestroy(ilink->ksp);CHKERRQ(ierr); 589 if (ilink->x) {ierr = VecDestroy(ilink->x);CHKERRQ(ierr);} 590 if (ilink->y) {ierr = VecDestroy(ilink->y);CHKERRQ(ierr);} 591 if (ilink->sctx) {ierr = VecScatterDestroy(ilink->sctx);CHKERRQ(ierr);} 592 if (ilink->is) {ierr = ISDestroy(ilink->is);CHKERRQ(ierr);} 593 if (ilink->cis) {ierr = ISDestroy(ilink->cis);CHKERRQ(ierr);} 594 next = ilink->next; 595 ierr = PetscFree(ilink->fields);CHKERRQ(ierr); 596 ierr = PetscFree(ilink);CHKERRQ(ierr); 597 ilink = next; 598 } 599 ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr); 600 if (jac->pmat) {ierr = MatDestroyMatrices(jac->nsplits,&jac->pmat);CHKERRQ(ierr);} 601 if (jac->Afield) {ierr = MatDestroyMatrices(jac->nsplits,&jac->Afield);CHKERRQ(ierr);} 602 if (jac->w1) {ierr = VecDestroy(jac->w1);CHKERRQ(ierr);} 603 if (jac->w2) {ierr = VecDestroy(jac->w2);CHKERRQ(ierr);} 604 if (jac->schur) {ierr = MatDestroy(jac->schur);CHKERRQ(ierr);} 605 if (jac->B) {ierr = MatDestroy(jac->B);CHKERRQ(ierr);} 606 if (jac->C) {ierr = MatDestroy(jac->C);CHKERRQ(ierr);} 607 ierr = PetscFree(jac);CHKERRQ(ierr); 608 PetscFunctionReturn(0); 609 } 610 611 #undef __FUNCT__ 612 #define __FUNCT__ "PCSetFromOptions_FieldSplit" 613 static PetscErrorCode PCSetFromOptions_FieldSplit(PC pc) 614 { 615 PetscErrorCode ierr; 616 PetscInt i = 0,nfields,*fields,bs; 617 PetscTruth flg; 618 char optionname[128]; 619 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 620 PCCompositeType ctype; 621 622 PetscFunctionBegin; 623 ierr = PetscOptionsHead("FieldSplit options");CHKERRQ(ierr); 624 ierr = PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);CHKERRQ(ierr); 625 if (flg) { 626 ierr = PCFieldSplitSetBlockSize(pc,bs);CHKERRQ(ierr); 627 } 628 629 ierr = PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg);CHKERRQ(ierr); 630 if (flg) { 631 ierr = PCFieldSplitSetType(pc,ctype);CHKERRQ(ierr); 632 } 633 634 if (jac->bs > 0) { 635 /* only allow user to set fields from command line if bs is already known. 636 otherwise user can set them in PCFieldSplitSetDefaults() */ 637 ierr = PetscMalloc(jac->bs*sizeof(PetscInt),&fields);CHKERRQ(ierr); 638 while (PETSC_TRUE) { 639 sprintf(optionname,"-pc_fieldsplit_%d_fields",(int)i++); 640 nfields = jac->bs; 641 ierr = PetscOptionsIntArray(optionname,"Fields in this split","PCFieldSplitSetFields",fields,&nfields,&flg);CHKERRQ(ierr); 642 if (!flg) break; 643 if (!nfields) SETERRQ(PETSC_ERR_USER,"Cannot list zero fields"); 644 ierr = PCFieldSplitSetFields(pc,nfields,fields);CHKERRQ(ierr); 645 } 646 ierr = PetscFree(fields);CHKERRQ(ierr); 647 } 648 ierr = PetscOptionsTail();CHKERRQ(ierr); 649 PetscFunctionReturn(0); 650 } 651 652 /*------------------------------------------------------------------------------------*/ 653 654 EXTERN_C_BEGIN 655 #undef __FUNCT__ 656 #define __FUNCT__ "PCFieldSplitSetFields_FieldSplit" 657 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetFields_FieldSplit(PC pc,PetscInt n,PetscInt *fields) 658 { 659 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 660 PetscErrorCode ierr; 661 PC_FieldSplitLink ilink,next = jac->head; 662 char prefix[128]; 663 PetscInt i; 664 665 PetscFunctionBegin; 666 if (n <= 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Negative number of fields requested"); 667 for (i=0; i<n; i++) { 668 if (fields[i] >= jac->bs) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs); 669 if (fields[i] < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]); 670 } 671 ierr = PetscNew(struct _PC_FieldSplitLink,&ilink);CHKERRQ(ierr); 672 ierr = PetscMalloc(n*sizeof(PetscInt),&ilink->fields);CHKERRQ(ierr); 673 ierr = PetscMemcpy(ilink->fields,fields,n*sizeof(PetscInt));CHKERRQ(ierr); 674 ilink->nfields = n; 675 ilink->next = PETSC_NULL; 676 ierr = KSPCreate(((PetscObject)pc)->comm,&ilink->ksp);CHKERRQ(ierr); 677 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 678 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 679 680 if (((PetscObject)pc)->prefix) { 681 sprintf(prefix,"%sfieldsplit_%d_",((PetscObject)pc)->prefix,(int)jac->nsplits); 682 } else { 683 sprintf(prefix,"fieldsplit_%d_",(int)jac->nsplits); 684 } 685 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 686 687 if (!next) { 688 jac->head = ilink; 689 ilink->previous = PETSC_NULL; 690 } else { 691 while (next->next) { 692 next = next->next; 693 } 694 next->next = ilink; 695 ilink->previous = next; 696 } 697 jac->nsplits++; 698 PetscFunctionReturn(0); 699 } 700 EXTERN_C_END 701 702 703 EXTERN_C_BEGIN 704 #undef __FUNCT__ 705 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit" 706 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp) 707 { 708 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 709 PetscErrorCode ierr; 710 PetscInt cnt = 0; 711 PC_FieldSplitLink ilink = jac->head; 712 713 PetscFunctionBegin; 714 ierr = PetscMalloc(jac->nsplits*sizeof(KSP*),subksp);CHKERRQ(ierr); 715 while (ilink) { 716 (*subksp)[cnt++] = ilink->ksp; 717 ilink = ilink->next; 718 } 719 if (cnt != jac->nsplits) SETERRQ2(PETSC_ERR_PLIB,"Corrupt PCFIELDSPLIT object: number splits in linked list %D in object %D",cnt,jac->nsplits); 720 *n = jac->nsplits; 721 PetscFunctionReturn(0); 722 } 723 EXTERN_C_END 724 725 EXTERN_C_BEGIN 726 #undef __FUNCT__ 727 #define __FUNCT__ "PCFieldSplitSetIS_FieldSplit" 728 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetIS_FieldSplit(PC pc,IS is) 729 { 730 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 731 PetscErrorCode ierr; 732 PC_FieldSplitLink ilink, next = jac->head; 733 char prefix[128]; 734 735 PetscFunctionBegin; 736 ierr = PetscNew(struct _PC_FieldSplitLink,&ilink);CHKERRQ(ierr); 737 ilink->is = is; 738 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 739 ilink->next = PETSC_NULL; 740 ierr = KSPCreate(((PetscObject)pc)->comm,&ilink->ksp);CHKERRQ(ierr); 741 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 742 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 743 744 if (((PetscObject)pc)->prefix) { 745 sprintf(prefix,"%sfieldsplit_%d_",((PetscObject)pc)->prefix,(int)jac->nsplits); 746 } else { 747 sprintf(prefix,"fieldsplit_%d_",(int)jac->nsplits); 748 } 749 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 750 751 if (!next) { 752 jac->head = ilink; 753 ilink->previous = PETSC_NULL; 754 } else { 755 while (next->next) { 756 next = next->next; 757 } 758 next->next = ilink; 759 ilink->previous = next; 760 } 761 jac->nsplits++; 762 763 PetscFunctionReturn(0); 764 } 765 EXTERN_C_END 766 767 #undef __FUNCT__ 768 #define __FUNCT__ "PCFieldSplitSetFields" 769 /*@ 770 PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner 771 772 Collective on PC 773 774 Input Parameters: 775 + pc - the preconditioner context 776 . n - the number of fields in this split 777 . fields - the fields in this split 778 779 Level: intermediate 780 781 Notes: Use PCFieldSplitSetIS() to set a completely general set of indices as a field. 782 783 The PCFieldSplitSetFields() is for defining fields as a strided blocks. For example, if the block 784 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 785 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x.... 786 where the numbered entries indicate what is in the field. 787 788 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize(), PCFieldSplitSetIS() 789 790 @*/ 791 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetFields(PC pc,PetscInt n, PetscInt *fields) 792 { 793 PetscErrorCode ierr,(*f)(PC,PetscInt,PetscInt *); 794 795 PetscFunctionBegin; 796 PetscValidHeaderSpecific(pc,PC_COOKIE,1); 797 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCFieldSplitSetFields_C",(void (**)(void))&f);CHKERRQ(ierr); 798 if (f) { 799 ierr = (*f)(pc,n,fields);CHKERRQ(ierr); 800 } 801 PetscFunctionReturn(0); 802 } 803 804 #undef __FUNCT__ 805 #define __FUNCT__ "PCFieldSplitSetIS" 806 /*@ 807 PCFieldSplitSetIS - Sets the exact elements for field 808 809 Collective on PC 810 811 Input Parameters: 812 + pc - the preconditioner context 813 . is - the index set that defines the vector elements in this field 814 815 816 Notes: Use PCFieldSplitSetFields(), for fields defined by strided types. 817 818 Level: intermediate 819 820 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize() 821 822 @*/ 823 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetIS(PC pc,IS is) 824 { 825 PetscErrorCode ierr,(*f)(PC,IS); 826 827 PetscFunctionBegin; 828 PetscValidHeaderSpecific(pc,PC_COOKIE,1); 829 PetscValidHeaderSpecific(is,IS_COOKIE,1); 830 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCFieldSplitSetIS_C",(void (**)(void))&f);CHKERRQ(ierr); 831 if (f) { 832 ierr = (*f)(pc,is);CHKERRQ(ierr); 833 } 834 PetscFunctionReturn(0); 835 } 836 837 #undef __FUNCT__ 838 #define __FUNCT__ "PCFieldSplitSetBlockSize" 839 /*@ 840 PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the 841 fieldsplit preconditioner. If not set the matrix block size is used. 842 843 Collective on PC 844 845 Input Parameters: 846 + pc - the preconditioner context 847 - bs - the block size 848 849 Level: intermediate 850 851 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields() 852 853 @*/ 854 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetBlockSize(PC pc,PetscInt bs) 855 { 856 PetscErrorCode ierr,(*f)(PC,PetscInt); 857 858 PetscFunctionBegin; 859 PetscValidHeaderSpecific(pc,PC_COOKIE,1); 860 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",(void (**)(void))&f);CHKERRQ(ierr); 861 if (f) { 862 ierr = (*f)(pc,bs);CHKERRQ(ierr); 863 } 864 PetscFunctionReturn(0); 865 } 866 867 #undef __FUNCT__ 868 #define __FUNCT__ "PCFieldSplitGetSubKSP" 869 /*@C 870 PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits 871 872 Collective on KSP 873 874 Input Parameter: 875 . pc - the preconditioner context 876 877 Output Parameters: 878 + n - the number of split 879 - pc - the array of KSP contexts 880 881 Note: 882 After PCFieldSplitGetSubKSP() the array of KSPs IS to be freed by the user 883 (not the KSP just the array that contains them). 884 885 You must call KSPSetUp() before calling PCFieldSplitGetSubKSP(). 886 887 Level: advanced 888 889 .seealso: PCFIELDSPLIT 890 @*/ 891 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[]) 892 { 893 PetscErrorCode ierr,(*f)(PC,PetscInt*,KSP **); 894 895 PetscFunctionBegin; 896 PetscValidHeaderSpecific(pc,PC_COOKIE,1); 897 PetscValidIntPointer(n,2); 898 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",(void (**)(void))&f);CHKERRQ(ierr); 899 if (f) { 900 ierr = (*f)(pc,n,subksp);CHKERRQ(ierr); 901 } else { 902 SETERRQ(PETSC_ERR_ARG_WRONG,"Cannot get subksp for this type of PC"); 903 } 904 PetscFunctionReturn(0); 905 } 906 907 EXTERN_C_BEGIN 908 #undef __FUNCT__ 909 #define __FUNCT__ "PCFieldSplitSetType_FieldSplit" 910 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type) 911 { 912 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 913 914 PetscFunctionBegin; 915 jac->type = type; 916 if (type == PC_COMPOSITE_SCHUR) { 917 pc->ops->apply = PCApply_FieldSplit_Schur; 918 pc->ops->view = PCView_FieldSplit_Schur; 919 } else { 920 pc->ops->apply = PCApply_FieldSplit; 921 pc->ops->view = PCView_FieldSplit; 922 } 923 PetscFunctionReturn(0); 924 } 925 EXTERN_C_END 926 927 EXTERN_C_BEGIN 928 #undef __FUNCT__ 929 #define __FUNCT__ "PCFieldSplitSetBlockSize_FieldSplit" 930 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs) 931 { 932 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 933 934 PetscFunctionBegin; 935 if (bs < 1) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs); 936 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); 937 jac->bs = bs; 938 PetscFunctionReturn(0); 939 } 940 EXTERN_C_END 941 942 #undef __FUNCT__ 943 #define __FUNCT__ "PCFieldSplitSetType" 944 /*@ 945 PCFieldSplitSetType - Sets the type of fieldsplit preconditioner. 946 947 Collective on PC 948 949 Input Parameter: 950 . pc - the preconditioner context 951 . type - PC_COMPOSITE_ADDITIVE (default), PC_COMPOSITE_MULTIPLICATIVE, PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE 952 953 Options Database Key: 954 . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special> - Sets fieldsplit preconditioner type 955 956 Level: Developer 957 958 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 959 960 .seealso: PCCompositeSetType() 961 962 @*/ 963 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetType(PC pc,PCCompositeType type) 964 { 965 PetscErrorCode ierr,(*f)(PC,PCCompositeType); 966 967 PetscFunctionBegin; 968 PetscValidHeaderSpecific(pc,PC_COOKIE,1); 969 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCFieldSplitSetType_C",(void (**)(void))&f);CHKERRQ(ierr); 970 if (f) { 971 ierr = (*f)(pc,type);CHKERRQ(ierr); 972 } 973 PetscFunctionReturn(0); 974 } 975 976 /* -------------------------------------------------------------------------------------*/ 977 /*MC 978 PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual 979 fields or groups of fields 980 981 982 To set options on the solvers for each block append -sub_ to all the PC 983 options database keys. For example, -sub_pc_type ilu -sub_pc_factor_levels 1 984 985 To set the options on the solvers separate for each block call PCFieldSplitGetSubKSP() 986 and set the options directly on the resulting KSP object 987 988 Level: intermediate 989 990 Options Database Keys: 991 + -pc_splitfield_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split 992 . -pc_splitfield_default - automatically add any fields to additional splits that have not 993 been supplied explicitly by -pc_splitfield_%d_fields 994 . -pc_splitfield_block_size <bs> - size of block that defines fields (i.e. there are bs fields) 995 - -pc_splitfield_type <additive,multiplicative> 996 997 - Options prefix for inner solvers when using Schur complement preconditioner are -schurAblock_ and -schur, 998 for all other solvers they are -fieldsplit_%d_ for the dth field, use -fieldsplit_ for all fields 999 1000 1001 Notes: use PCFieldSplitSetFields() to set fields defined by "strided" entries and PCFieldSplitSetIS() 1002 to define a field by an arbitrary collection of entries. 1003 1004 If no fields are set the default is used. The fields are defined by entries strided by bs, 1005 beginning at 0 then 1, etc to bs-1. The block size can be set with PCFieldSplitSetBlockSize(), 1006 if this is not called the block size defaults to the blocksize of the second matrix passed 1007 to KSPSetOperators()/PCSetOperators(). 1008 1009 Currently for the multiplicative version, the updated residual needed for the next field 1010 solve is computed via a matrix vector product over the entire array. An optimization would be 1011 to update the residual only for the part of the right hand side associated with the next field 1012 solve. (This would involve more MatGetSubMatrix() calls or some other mechanism to compute the 1013 part of the matrix needed to just update part of the residual). 1014 1015 Concepts: physics based preconditioners 1016 1017 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, 1018 PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS() 1019 M*/ 1020 1021 EXTERN_C_BEGIN 1022 #undef __FUNCT__ 1023 #define __FUNCT__ "PCCreate_FieldSplit" 1024 PetscErrorCode PETSCKSP_DLLEXPORT PCCreate_FieldSplit(PC pc) 1025 { 1026 PetscErrorCode ierr; 1027 PC_FieldSplit *jac; 1028 1029 PetscFunctionBegin; 1030 ierr = PetscNewLog(pc,PC_FieldSplit,&jac);CHKERRQ(ierr); 1031 jac->bs = -1; 1032 jac->nsplits = 0; 1033 jac->type = PC_COMPOSITE_MULTIPLICATIVE; 1034 1035 pc->data = (void*)jac; 1036 1037 pc->ops->apply = PCApply_FieldSplit; 1038 pc->ops->applytranspose = PCApplyTranspose_FieldSplit; 1039 pc->ops->setup = PCSetUp_FieldSplit; 1040 pc->ops->destroy = PCDestroy_FieldSplit; 1041 pc->ops->setfromoptions = PCSetFromOptions_FieldSplit; 1042 pc->ops->view = PCView_FieldSplit; 1043 pc->ops->applyrichardson = 0; 1044 1045 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit", 1046 PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 1047 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetFields_C","PCFieldSplitSetFields_FieldSplit", 1048 PCFieldSplitSetFields_FieldSplit);CHKERRQ(ierr); 1049 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetIS_C","PCFieldSplitSetIS_FieldSplit", 1050 PCFieldSplitSetIS_FieldSplit);CHKERRQ(ierr); 1051 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetType_C","PCFieldSplitSetType_FieldSplit", 1052 PCFieldSplitSetType_FieldSplit);CHKERRQ(ierr); 1053 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetBlockSize_C","PCFieldSplitSetBlockSize_FieldSplit", 1054 PCFieldSplitSetBlockSize_FieldSplit);CHKERRQ(ierr); 1055 PetscFunctionReturn(0); 1056 } 1057 EXTERN_C_END 1058 1059 1060