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