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