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