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; /* additive or multiplicative */ 22 PetscTruth defaultsplit; 23 PetscInt bs; 24 PetscInt nsplits; 25 Vec *x,*y,w1,w2; 26 Mat *pmat; 27 PetscTruth issetup; 28 PC_FieldSplitLink head; 29 } PC_FieldSplit; 30 31 /* 32 Notes: there is no particular reason that pmat, x, and y are stored as arrays in PC_FieldSplit instead of 33 inside PC_FieldSplitLink, just historical. If you want to be able to add new fields after already using the 34 PC you could change this. 35 */ 36 #undef __FUNCT__ 37 #define __FUNCT__ "PCView_FieldSplit" 38 static PetscErrorCode PCView_FieldSplit(PC pc,PetscViewer viewer) 39 { 40 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 41 PetscErrorCode ierr; 42 PetscTruth iascii; 43 PetscInt i,j; 44 PC_FieldSplitLink ilink = jac->head; 45 46 PetscFunctionBegin; 47 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 48 if (iascii) { 49 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D, blocksize = %D\n",PCCompositeTypes[jac->type],jac->nsplits,jac->bs);CHKERRQ(ierr); 50 ierr = PetscViewerASCIIPrintf(viewer," Solver info for each split is in the following KSP objects:\n");CHKERRQ(ierr); 51 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 52 for (i=0; i<jac->nsplits; i++) { 53 if (ilink->fields) { 54 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);CHKERRQ(ierr); 55 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 56 for (j=0; j<ilink->nfields; j++) { 57 if (j > 0) { 58 ierr = PetscViewerASCIIPrintf(viewer,",");CHKERRQ(ierr); 59 } 60 ierr = PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);CHKERRQ(ierr); 61 } 62 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 63 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 64 } else { 65 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i);CHKERRQ(ierr); 66 } 67 ierr = KSPView(ilink->ksp,viewer);CHKERRQ(ierr); 68 ilink = ilink->next; 69 } 70 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 71 } else { 72 SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for PCFieldSplit",((PetscObject)viewer)->type_name); 73 } 74 PetscFunctionReturn(0); 75 } 76 77 #undef __FUNCT__ 78 #define __FUNCT__ "PCFieldSplitSetDefaults" 79 static PetscErrorCode PCFieldSplitSetDefaults(PC pc) 80 { 81 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 82 PetscErrorCode ierr; 83 PC_FieldSplitLink ilink = jac->head; 84 PetscInt i; 85 PetscTruth flg = PETSC_FALSE,*fields; 86 87 PetscFunctionBegin; 88 ierr = PetscOptionsGetTruth(((PetscObject)pc)->prefix,"-pc_fieldsplit_default",&flg,PETSC_NULL);CHKERRQ(ierr); 89 if (!ilink || flg) { 90 ierr = PetscInfo(pc,"Using default splitting of fields\n");CHKERRQ(ierr); 91 if (jac->bs <= 0) { 92 if (pc->pmat) { 93 ierr = MatGetBlockSize(pc->pmat,&jac->bs);CHKERRQ(ierr); 94 } else { 95 jac->bs = 1; 96 } 97 } 98 ierr = PetscMalloc(jac->bs*sizeof(PetscTruth),&fields);CHKERRQ(ierr); 99 ierr = PetscMemzero(fields,jac->bs*sizeof(PetscTruth));CHKERRQ(ierr); 100 while (ilink) { 101 for (i=0; i<ilink->nfields; i++) { 102 fields[ilink->fields[i]] = PETSC_TRUE; 103 } 104 ilink = ilink->next; 105 } 106 jac->defaultsplit = PETSC_TRUE; 107 for (i=0; i<jac->bs; i++) { 108 if (!fields[i]) { 109 ierr = PCFieldSplitSetFields(pc,1,&i);CHKERRQ(ierr); 110 } else { 111 jac->defaultsplit = PETSC_FALSE; 112 } 113 } 114 } 115 PetscFunctionReturn(0); 116 } 117 118 119 #undef __FUNCT__ 120 #define __FUNCT__ "PCSetUp_FieldSplit" 121 static PetscErrorCode PCSetUp_FieldSplit(PC pc) 122 { 123 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 124 PetscErrorCode ierr; 125 PC_FieldSplitLink ilink; 126 PetscInt i,nsplit,ccsize; 127 MatStructure flag = pc->flag; 128 PetscTruth sorted; 129 130 PetscFunctionBegin; 131 ierr = PCFieldSplitSetDefaults(pc);CHKERRQ(ierr); 132 nsplit = jac->nsplits; 133 ilink = jac->head; 134 135 /* get the matrices for each split */ 136 if (!jac->issetup) { 137 PetscInt rstart,rend,nslots,bs; 138 139 jac->issetup = PETSC_TRUE; 140 141 /* This is done here instead of in PCFieldSplitSetFields() because may not have matrix at that point */ 142 bs = jac->bs; 143 ierr = MatGetOwnershipRange(pc->pmat,&rstart,&rend);CHKERRQ(ierr); 144 ierr = MatGetLocalSize(pc->pmat,PETSC_NULL,&ccsize);CHKERRQ(ierr); 145 nslots = (rend - rstart)/bs; 146 for (i=0; i<nsplit; i++) { 147 if (jac->defaultsplit) { 148 ierr = ISCreateStride(((PetscObject)pc)->comm,nslots,rstart+i,nsplit,&ilink->is);CHKERRQ(ierr); 149 ilink->csize = ccsize/nsplit; 150 } else if (!ilink->is) { 151 if (ilink->nfields > 1) { 152 PetscInt *ii,j,k,nfields = ilink->nfields,*fields = ilink->fields; 153 ierr = PetscMalloc(ilink->nfields*nslots*sizeof(PetscInt),&ii);CHKERRQ(ierr); 154 for (j=0; j<nslots; j++) { 155 for (k=0; k<nfields; k++) { 156 ii[nfields*j + k] = rstart + bs*j + fields[k]; 157 } 158 } 159 ierr = ISCreateGeneral(((PetscObject)pc)->comm,nslots*nfields,ii,&ilink->is);CHKERRQ(ierr); 160 ierr = PetscFree(ii);CHKERRQ(ierr); 161 } else { 162 ierr = ISCreateStride(((PetscObject)pc)->comm,nslots,rstart+ilink->fields[0],bs,&ilink->is);CHKERRQ(ierr); 163 } 164 ilink->csize = (ccsize/bs)*ilink->nfields; 165 ierr = ISSorted(ilink->is,&sorted);CHKERRQ(ierr); 166 if (!sorted) SETERRQ(PETSC_ERR_USER,"Fields must be sorted when creating split"); 167 } 168 ierr = ISAllGather(ilink->is,&ilink->cis);CHKERRQ(ierr); 169 ilink = ilink->next; 170 } 171 } 172 173 ilink = jac->head; 174 if (!jac->pmat) { 175 ierr = PetscMalloc(nsplit*sizeof(Mat),&jac->pmat);CHKERRQ(ierr); 176 for (i=0; i<nsplit; i++) { 177 ierr = MatGetSubMatrix(pc->pmat,ilink->is,ilink->cis,ilink->csize,MAT_INITIAL_MATRIX,&jac->pmat[i]);CHKERRQ(ierr); 178 ilink = ilink->next; 179 } 180 } else { 181 for (i=0; i<nsplit; i++) { 182 ierr = MatGetSubMatrix(pc->pmat,ilink->is,ilink->cis,ilink->csize,MAT_REUSE_MATRIX,&jac->pmat[i]);CHKERRQ(ierr); 183 ilink = ilink->next; 184 } 185 } 186 187 /* set up the individual PCs */ 188 i = 0; 189 ilink = jac->head; 190 while (ilink) { 191 ierr = KSPSetOperators(ilink->ksp,jac->pmat[i],jac->pmat[i],flag);CHKERRQ(ierr); 192 ierr = KSPSetFromOptions(ilink->ksp);CHKERRQ(ierr); 193 ierr = KSPSetUp(ilink->ksp);CHKERRQ(ierr); 194 i++; 195 ilink = ilink->next; 196 } 197 198 /* create work vectors for each split */ 199 if (!jac->x) { 200 Vec xtmp; 201 ierr = PetscMalloc2(nsplit,Vec,&jac->x,nsplit,Vec,&jac->y);CHKERRQ(ierr); 202 ilink = jac->head; 203 for (i=0; i<nsplit; i++) { 204 Vec *vl,*vr; 205 206 ierr = KSPGetVecs(ilink->ksp,1,&vr,1,&vl);CHKERRQ(ierr); 207 ilink->x = *vr; 208 ilink->y = *vl; 209 ierr = PetscFree(vr);CHKERRQ(ierr); 210 ierr = PetscFree(vl);CHKERRQ(ierr); 211 jac->x[i] = ilink->x; 212 jac->y[i] = ilink->y; 213 ilink = ilink->next; 214 } 215 /* compute scatter contexts needed by multiplicative versions and non-default splits */ 216 217 ilink = jac->head; 218 ierr = MatGetVecs(pc->pmat,&xtmp,PETSC_NULL);CHKERRQ(ierr); 219 for (i=0; i<nsplit; i++) { 220 ierr = VecScatterCreate(xtmp,ilink->is,jac->x[i],PETSC_NULL,&ilink->sctx);CHKERRQ(ierr); 221 ilink = ilink->next; 222 } 223 ierr = VecDestroy(xtmp);CHKERRQ(ierr); 224 } 225 PetscFunctionReturn(0); 226 } 227 228 #define FieldSplitSplitSolveAdd(ilink,xx,yy) \ 229 (VecScatterBegin(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 230 VecScatterEnd(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 231 KSPSolve(ilink->ksp,ilink->x,ilink->y) || \ 232 VecScatterBegin(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE) || \ 233 VecScatterEnd(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE)) 234 235 #undef __FUNCT__ 236 #define __FUNCT__ "PCApply_FieldSplit" 237 static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y) 238 { 239 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 240 PetscErrorCode ierr; 241 PC_FieldSplitLink ilink = jac->head; 242 PetscInt bs; 243 244 PetscFunctionBegin; 245 CHKMEMQ; 246 ierr = VecGetBlockSize(x,&bs);CHKERRQ(ierr); 247 ierr = VecSetBlockSize(x,jac->bs);CHKERRQ(ierr); 248 ierr = VecSetBlockSize(y,jac->bs);CHKERRQ(ierr); 249 250 if (jac->type == PC_COMPOSITE_ADDITIVE) { 251 if (jac->defaultsplit) { 252 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 253 while (ilink) { 254 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 255 ilink = ilink->next; 256 } 257 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 258 } else { 259 ierr = VecSet(y,0.0);CHKERRQ(ierr); 260 while (ilink) { 261 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 262 ilink = ilink->next; 263 } 264 } 265 } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 266 if (!jac->w1) { 267 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 268 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 269 } 270 ierr = VecSet(y,0.0);CHKERRQ(ierr); 271 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 272 while (ilink->next) { 273 ilink = ilink->next; 274 ierr = MatMult(pc->mat,y,jac->w1);CHKERRQ(ierr); 275 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 276 ierr = FieldSplitSplitSolveAdd(ilink,jac->w2,y);CHKERRQ(ierr); 277 } 278 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 279 while (ilink->previous) { 280 ilink = ilink->previous; 281 ierr = MatMult(pc->mat,y,jac->w1);CHKERRQ(ierr); 282 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 283 ierr = FieldSplitSplitSolveAdd(ilink,jac->w2,y);CHKERRQ(ierr); 284 } 285 } 286 } else SETERRQ1(PETSC_ERR_SUP,"Unsupported or unknown composition",(int) jac->type); 287 CHKMEMQ; 288 PetscFunctionReturn(0); 289 } 290 291 #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \ 292 (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 293 VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 294 KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) || \ 295 VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \ 296 VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE)) 297 298 #undef __FUNCT__ 299 #define __FUNCT__ "PCApply_FieldSplit" 300 static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y) 301 { 302 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 303 PetscErrorCode ierr; 304 PC_FieldSplitLink ilink = jac->head; 305 PetscInt bs; 306 307 PetscFunctionBegin; 308 CHKMEMQ; 309 ierr = VecGetBlockSize(x,&bs);CHKERRQ(ierr); 310 ierr = VecSetBlockSize(x,jac->bs);CHKERRQ(ierr); 311 ierr = VecSetBlockSize(y,jac->bs);CHKERRQ(ierr); 312 313 if (jac->type == PC_COMPOSITE_ADDITIVE) { 314 if (jac->defaultsplit) { 315 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 316 while (ilink) { 317 ierr = KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 318 ilink = ilink->next; 319 } 320 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 321 } else { 322 ierr = VecSet(y,0.0);CHKERRQ(ierr); 323 while (ilink) { 324 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 325 ilink = ilink->next; 326 } 327 } 328 } else { 329 if (!jac->w1) { 330 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 331 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 332 } 333 ierr = VecSet(y,0.0);CHKERRQ(ierr); 334 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 335 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 336 while (ilink->next) { 337 ilink = ilink->next; 338 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 339 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 340 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 341 } 342 while (ilink->previous) { 343 ilink = ilink->previous; 344 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 345 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 346 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 347 } 348 } else { 349 while (ilink->next) { /* get to last entry in linked list */ 350 ilink = ilink->next; 351 } 352 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 353 while (ilink->previous) { 354 ilink = ilink->previous; 355 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 356 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 357 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 358 } 359 } 360 } 361 CHKMEMQ; 362 PetscFunctionReturn(0); 363 } 364 365 #undef __FUNCT__ 366 #define __FUNCT__ "PCDestroy_FieldSplit" 367 static PetscErrorCode PCDestroy_FieldSplit(PC pc) 368 { 369 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 370 PetscErrorCode ierr; 371 PC_FieldSplitLink ilink = jac->head,next; 372 373 PetscFunctionBegin; 374 while (ilink) { 375 ierr = KSPDestroy(ilink->ksp);CHKERRQ(ierr); 376 if (ilink->x) {ierr = VecDestroy(ilink->x);CHKERRQ(ierr);} 377 if (ilink->y) {ierr = VecDestroy(ilink->y);CHKERRQ(ierr);} 378 if (ilink->sctx) {ierr = VecScatterDestroy(ilink->sctx);CHKERRQ(ierr);} 379 if (ilink->is) {ierr = ISDestroy(ilink->is);CHKERRQ(ierr);} 380 if (ilink->cis) {ierr = ISDestroy(ilink->cis);CHKERRQ(ierr);} 381 next = ilink->next; 382 ierr = PetscFree(ilink->fields);CHKERRQ(ierr); 383 ierr = PetscFree(ilink);CHKERRQ(ierr); 384 ilink = next; 385 } 386 ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr); 387 if (jac->pmat) {ierr = MatDestroyMatrices(jac->nsplits,&jac->pmat);CHKERRQ(ierr);} 388 if (jac->w1) {ierr = VecDestroy(jac->w1);CHKERRQ(ierr);} 389 if (jac->w2) {ierr = VecDestroy(jac->w2);CHKERRQ(ierr);} 390 ierr = PetscFree(jac);CHKERRQ(ierr); 391 PetscFunctionReturn(0); 392 } 393 394 #undef __FUNCT__ 395 #define __FUNCT__ "PCSetFromOptions_FieldSplit" 396 static PetscErrorCode PCSetFromOptions_FieldSplit(PC pc) 397 { 398 PetscErrorCode ierr; 399 PetscInt i = 0,nfields,*fields,bs; 400 PetscTruth flg; 401 char optionname[128]; 402 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 403 404 PetscFunctionBegin; 405 ierr = PetscOptionsHead("FieldSplit options");CHKERRQ(ierr); 406 ierr = PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);CHKERRQ(ierr); 407 if (flg) { 408 ierr = PCFieldSplitSetBlockSize(pc,bs);CHKERRQ(ierr); 409 } 410 if (jac->bs <= 0) { 411 412 ierr = PCFieldSplitSetBlockSize(pc,1);CHKERRQ(ierr); 413 } 414 ierr = PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&jac->type,&flg);CHKERRQ(ierr); 415 ierr = PetscMalloc(jac->bs*sizeof(PetscInt),&fields);CHKERRQ(ierr); 416 while (PETSC_TRUE) { 417 sprintf(optionname,"-pc_fieldsplit_%d_fields",(int)i++); 418 nfields = jac->bs; 419 ierr = PetscOptionsIntArray(optionname,"Fields in this split","PCFieldSplitSetFields",fields,&nfields,&flg);CHKERRQ(ierr); 420 if (!flg) break; 421 if (!nfields) SETERRQ(PETSC_ERR_USER,"Cannot list zero fields"); 422 ierr = PCFieldSplitSetFields(pc,nfields,fields);CHKERRQ(ierr); 423 } 424 ierr = PetscFree(fields);CHKERRQ(ierr); 425 ierr = PetscOptionsTail();CHKERRQ(ierr); 426 PetscFunctionReturn(0); 427 } 428 429 /*------------------------------------------------------------------------------------*/ 430 431 EXTERN_C_BEGIN 432 #undef __FUNCT__ 433 #define __FUNCT__ "PCFieldSplitSetFields_FieldSplit" 434 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetFields_FieldSplit(PC pc,PetscInt n,PetscInt *fields) 435 { 436 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 437 PetscErrorCode ierr; 438 PC_FieldSplitLink ilink,next = jac->head; 439 char prefix[128]; 440 PetscInt i; 441 442 PetscFunctionBegin; 443 if (n <= 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Negative number of fields requested"); 444 for (i=0; i<n; i++) { 445 if (fields[i] >= jac->bs) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs); 446 if (fields[i] < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]); 447 } 448 ierr = PetscNew(struct _PC_FieldSplitLink,&ilink);CHKERRQ(ierr); 449 ierr = PetscMalloc(n*sizeof(PetscInt),&ilink->fields);CHKERRQ(ierr); 450 ierr = PetscMemcpy(ilink->fields,fields,n*sizeof(PetscInt));CHKERRQ(ierr); 451 ilink->nfields = n; 452 ilink->next = PETSC_NULL; 453 ierr = KSPCreate(((PetscObject)pc)->comm,&ilink->ksp);CHKERRQ(ierr); 454 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 455 456 if (((PetscObject)pc)->prefix) { 457 sprintf(prefix,"%sfieldsplit_%d_",((PetscObject)pc)->prefix,(int)jac->nsplits); 458 } else { 459 sprintf(prefix,"fieldsplit_%d_",(int)jac->nsplits); 460 } 461 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 462 463 if (!next) { 464 jac->head = ilink; 465 ilink->previous = PETSC_NULL; 466 } else { 467 while (next->next) { 468 next = next->next; 469 } 470 next->next = ilink; 471 ilink->previous = next; 472 } 473 jac->nsplits++; 474 PetscFunctionReturn(0); 475 } 476 EXTERN_C_END 477 478 479 EXTERN_C_BEGIN 480 #undef __FUNCT__ 481 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit" 482 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp) 483 { 484 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 485 PetscErrorCode ierr; 486 PetscInt cnt = 0; 487 PC_FieldSplitLink ilink = jac->head; 488 489 PetscFunctionBegin; 490 ierr = PetscMalloc(jac->nsplits*sizeof(KSP*),subksp);CHKERRQ(ierr); 491 while (ilink) { 492 (*subksp)[cnt++] = ilink->ksp; 493 ilink = ilink->next; 494 } 495 if (cnt != jac->nsplits) SETERRQ2(PETSC_ERR_PLIB,"Corrupt PCFIELDSPLIT object: number splits in linked list %D in object %D",cnt,jac->nsplits); 496 *n = jac->nsplits; 497 PetscFunctionReturn(0); 498 } 499 EXTERN_C_END 500 501 EXTERN_C_BEGIN 502 #undef __FUNCT__ 503 #define __FUNCT__ "PCFieldSplitSetIS_FieldSplit" 504 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetIS_FieldSplit(PC pc,IS is) 505 { 506 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 507 PetscErrorCode ierr; 508 PC_FieldSplitLink ilink, next = jac->head; 509 char prefix[128]; 510 511 PetscFunctionBegin; 512 ierr = PetscNew(struct _PC_FieldSplitLink,&ilink);CHKERRQ(ierr); 513 ilink->is = is; 514 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 515 ilink->next = PETSC_NULL; 516 ierr = KSPCreate(((PetscObject)pc)->comm,&ilink->ksp);CHKERRQ(ierr); 517 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 518 519 if (((PetscObject)pc)->prefix) { 520 sprintf(prefix,"%sfieldsplit_%d_",((PetscObject)pc)->prefix,(int)jac->nsplits); 521 } else { 522 sprintf(prefix,"fieldsplit_%d_",(int)jac->nsplits); 523 } 524 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 525 526 if (!next) { 527 jac->head = ilink; 528 ilink->previous = PETSC_NULL; 529 } else { 530 while (next->next) { 531 next = next->next; 532 } 533 next->next = ilink; 534 ilink->previous = next; 535 } 536 jac->nsplits++; 537 538 PetscFunctionReturn(0); 539 } 540 EXTERN_C_END 541 542 #undef __FUNCT__ 543 #define __FUNCT__ "PCFieldSplitSetFields" 544 /*@ 545 PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner 546 547 Collective on PC 548 549 Input Parameters: 550 + pc - the preconditioner context 551 . n - the number of fields in this split 552 . fields - the fields in this split 553 554 Level: intermediate 555 556 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize() 557 558 @*/ 559 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetFields(PC pc,PetscInt n, PetscInt *fields) 560 { 561 PetscErrorCode ierr,(*f)(PC,PetscInt,PetscInt *); 562 563 PetscFunctionBegin; 564 PetscValidHeaderSpecific(pc,PC_COOKIE,1); 565 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCFieldSplitSetFields_C",(void (**)(void))&f);CHKERRQ(ierr); 566 if (f) { 567 ierr = (*f)(pc,n,fields);CHKERRQ(ierr); 568 } 569 PetscFunctionReturn(0); 570 } 571 572 #undef __FUNCT__ 573 #define __FUNCT__ "PCFieldSplitSetIS" 574 /*@ 575 PCFieldSplitSetIS - Sets the exact elements for field 576 577 Collective on PC 578 579 Input Parameters: 580 + pc - the preconditioner context 581 . is - the index set that defines the vector elements in this field 582 583 Level: intermediate 584 585 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize() 586 587 @*/ 588 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetIS(PC pc,IS is) 589 { 590 PetscErrorCode ierr,(*f)(PC,IS); 591 592 PetscFunctionBegin; 593 PetscValidHeaderSpecific(pc,PC_COOKIE,1); 594 PetscValidHeaderSpecific(is,IS_COOKIE,1); 595 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCFieldSplitSetIS_C",(void (**)(void))&f);CHKERRQ(ierr); 596 if (f) { 597 ierr = (*f)(pc,is);CHKERRQ(ierr); 598 } 599 PetscFunctionReturn(0); 600 } 601 602 #undef __FUNCT__ 603 #define __FUNCT__ "PCFieldSplitSetBlockSize" 604 /*@ 605 PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the 606 fieldsplit preconditioner. If not set the matrix block size is used. 607 608 Collective on PC 609 610 Input Parameters: 611 + pc - the preconditioner context 612 - bs - the block size 613 614 Level: intermediate 615 616 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields() 617 618 @*/ 619 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetBlockSize(PC pc,PetscInt bs) 620 { 621 PetscErrorCode ierr,(*f)(PC,PetscInt); 622 623 PetscFunctionBegin; 624 PetscValidHeaderSpecific(pc,PC_COOKIE,1); 625 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",(void (**)(void))&f);CHKERRQ(ierr); 626 if (f) { 627 ierr = (*f)(pc,bs);CHKERRQ(ierr); 628 } 629 PetscFunctionReturn(0); 630 } 631 632 #undef __FUNCT__ 633 #define __FUNCT__ "PCFieldSplitGetSubKSP" 634 /*@C 635 PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits 636 637 Collective on KSP 638 639 Input Parameter: 640 . pc - the preconditioner context 641 642 Output Parameters: 643 + n - the number of split 644 - pc - the array of KSP contexts 645 646 Note: 647 After PCFieldSplitGetSubKSP() the array of KSPs IS to be freed 648 649 You must call KSPSetUp() before calling PCFieldSplitGetSubKSP(). 650 651 Level: advanced 652 653 .seealso: PCFIELDSPLIT 654 @*/ 655 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[]) 656 { 657 PetscErrorCode ierr,(*f)(PC,PetscInt*,KSP **); 658 659 PetscFunctionBegin; 660 PetscValidHeaderSpecific(pc,PC_COOKIE,1); 661 PetscValidIntPointer(n,2); 662 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",(void (**)(void))&f);CHKERRQ(ierr); 663 if (f) { 664 ierr = (*f)(pc,n,subksp);CHKERRQ(ierr); 665 } else { 666 SETERRQ(PETSC_ERR_ARG_WRONG,"Cannot get subksp for this type of PC"); 667 } 668 PetscFunctionReturn(0); 669 } 670 671 EXTERN_C_BEGIN 672 #undef __FUNCT__ 673 #define __FUNCT__ "PCFieldSplitSetType_FieldSplit" 674 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type) 675 { 676 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 677 678 PetscFunctionBegin; 679 jac->type = type; 680 PetscFunctionReturn(0); 681 } 682 EXTERN_C_END 683 684 EXTERN_C_BEGIN 685 #undef __FUNCT__ 686 #define __FUNCT__ "PCFieldSplitSetBlockSize_FieldSplit" 687 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs) 688 { 689 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 690 691 PetscFunctionBegin; 692 if (bs < 1) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs); 693 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); 694 jac->bs = bs; 695 PetscFunctionReturn(0); 696 } 697 EXTERN_C_END 698 699 #undef __FUNCT__ 700 #define __FUNCT__ "PCFieldSplitSetType" 701 /*@ 702 PCFieldSplitSetType - Sets the type of fieldsplit preconditioner. 703 704 Collective on PC 705 706 Input Parameter: 707 . pc - the preconditioner context 708 . type - PC_COMPOSITE_ADDITIVE (default), PC_COMPOSITE_MULTIPLICATIVE, PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE 709 710 Options Database Key: 711 . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special> - Sets fieldsplit preconditioner type 712 713 Level: Developer 714 715 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 716 717 .seealso: PCCompositeSetType() 718 719 @*/ 720 PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetType(PC pc,PCCompositeType type) 721 { 722 PetscErrorCode ierr,(*f)(PC,PCCompositeType); 723 724 PetscFunctionBegin; 725 PetscValidHeaderSpecific(pc,PC_COOKIE,1); 726 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCFieldSplitSetType_C",(void (**)(void))&f);CHKERRQ(ierr); 727 if (f) { 728 ierr = (*f)(pc,type);CHKERRQ(ierr); 729 } 730 PetscFunctionReturn(0); 731 } 732 733 /* -------------------------------------------------------------------------------------*/ 734 /*MC 735 PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual 736 fields or groups of fields 737 738 739 To set options on the solvers for each block append -sub_ to all the PC 740 options database keys. For example, -sub_pc_type ilu -sub_pc_factor_levels 1 741 742 To set the options on the solvers separate for each block call PCFieldSplitGetSubKSP() 743 and set the options directly on the resulting KSP object 744 745 Level: intermediate 746 747 Options Database Keys: 748 + -pc_splitfield_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split 749 . -pc_splitfield_default - automatically add any fields to additional splits that have not 750 been supplied explicitly by -pc_splitfield_%d_fields 751 . -pc_splitfield_block_size <bs> - size of block that defines fields (i.e. there are bs fields) 752 - -pc_splitfield_type <additive,multiplicative> 753 754 Concepts: physics based preconditioners 755 756 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, 757 PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(),PCFieldSplitSetType() 758 M*/ 759 760 EXTERN_C_BEGIN 761 #undef __FUNCT__ 762 #define __FUNCT__ "PCCreate_FieldSplit" 763 PetscErrorCode PETSCKSP_DLLEXPORT PCCreate_FieldSplit(PC pc) 764 { 765 PetscErrorCode ierr; 766 PC_FieldSplit *jac; 767 768 PetscFunctionBegin; 769 ierr = PetscNewLog(pc,PC_FieldSplit,&jac);CHKERRQ(ierr); 770 jac->bs = -1; 771 jac->nsplits = 0; 772 jac->type = PC_COMPOSITE_ADDITIVE; 773 774 pc->data = (void*)jac; 775 776 pc->ops->apply = PCApply_FieldSplit; 777 pc->ops->applytranspose = PCApplyTranspose_FieldSplit; 778 pc->ops->setup = PCSetUp_FieldSplit; 779 pc->ops->destroy = PCDestroy_FieldSplit; 780 pc->ops->setfromoptions = PCSetFromOptions_FieldSplit; 781 pc->ops->view = PCView_FieldSplit; 782 pc->ops->applyrichardson = 0; 783 784 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit", 785 PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 786 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetFields_C","PCFieldSplitSetFields_FieldSplit", 787 PCFieldSplitSetFields_FieldSplit);CHKERRQ(ierr); 788 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetIS_C","PCFieldSplitSetIS_FieldSplit", 789 PCFieldSplitSetIS_FieldSplit);CHKERRQ(ierr); 790 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetType_C","PCFieldSplitSetType_FieldSplit", 791 PCFieldSplitSetType_FieldSplit);CHKERRQ(ierr); 792 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetBlockSize_C","PCFieldSplitSetBlockSize_FieldSplit", 793 PCFieldSplitSetBlockSize_FieldSplit);CHKERRQ(ierr); 794 PetscFunctionReturn(0); 795 } 796 EXTERN_C_END 797 798 799