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