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