Actual source code: pbjacobi.c
1: #define PETSCKSP_DLL
3: /*
4: Include files needed for the PBJacobi preconditioner:
5: pcimpl.h - private include file intended for use by all preconditioners
6: */
8: #include private/pcimpl.h
10: /*
11: Private context (data structure) for the PBJacobi preconditioner.
12: */
13: typedef struct {
14: PetscScalar *diag;
15: PetscInt bs,mbs;
16: } PC_PBJacobi;
18: /*
19: Currently only implemented for baij matrices and directly access baij
20: data structures.
21: */
22: #include src/mat/impls/baij/mpi/mpibaij.h
23: #include src/inline/ilu.h
27: static PetscErrorCode PCApply_PBJacobi_2(PC pc,Vec x,Vec y)
28: {
29: PC_PBJacobi *jac = (PC_PBJacobi*)pc->data;
31: PetscInt i,m = jac->mbs;
32: PetscScalar *diag = jac->diag,x0,x1,*xx,*yy;
33:
35: VecGetArray(x,&xx);
36: VecGetArray(y,&yy);
37: for (i=0; i<m; i++) {
38: x0 = xx[2*i]; x1 = xx[2*i+1];
39: yy[2*i] = diag[0]*x0 + diag[2]*x1;
40: yy[2*i+1] = diag[1]*x0 + diag[3]*x1;
41: diag += 4;
42: }
43: VecRestoreArray(x,&xx);
44: VecRestoreArray(y,&yy);
45: PetscLogFlops(6*m);
46: return(0);
47: }
50: static PetscErrorCode PCApply_PBJacobi_3(PC pc,Vec x,Vec y)
51: {
52: PC_PBJacobi *jac = (PC_PBJacobi*)pc->data;
54: PetscInt i,m = jac->mbs;
55: PetscScalar *diag = jac->diag,x0,x1,x2,*xx,*yy;
56:
58: VecGetArray(x,&xx);
59: VecGetArray(y,&yy);
60: for (i=0; i<m; i++) {
61: x0 = xx[3*i]; x1 = xx[3*i+1]; x2 = xx[3*i+2];
62: yy[3*i] = diag[0]*x0 + diag[3]*x1 + diag[6]*x2;
63: yy[3*i+1] = diag[1]*x0 + diag[4]*x1 + diag[7]*x2;
64: yy[3*i+2] = diag[2]*x0 + diag[5]*x1 + diag[8]*x2;
65: diag += 9;
66: }
67: VecRestoreArray(x,&xx);
68: VecRestoreArray(y,&yy);
69: PetscLogFlops(15*m);
70: return(0);
71: }
74: static PetscErrorCode PCApply_PBJacobi_4(PC pc,Vec x,Vec y)
75: {
76: PC_PBJacobi *jac = (PC_PBJacobi*)pc->data;
78: PetscInt i,m = jac->mbs;
79: PetscScalar *diag = jac->diag,x0,x1,x2,x3,*xx,*yy;
80:
82: VecGetArray(x,&xx);
83: VecGetArray(y,&yy);
84: for (i=0; i<m; i++) {
85: x0 = xx[4*i]; x1 = xx[4*i+1]; x2 = xx[4*i+2]; x3 = xx[4*i+3];
86: yy[4*i] = diag[0]*x0 + diag[4]*x1 + diag[8]*x2 + diag[12]*x3;
87: yy[4*i+1] = diag[1]*x0 + diag[5]*x1 + diag[9]*x2 + diag[13]*x3;
88: yy[4*i+2] = diag[2]*x0 + diag[6]*x1 + diag[10]*x2 + diag[14]*x3;
89: yy[4*i+3] = diag[3]*x0 + diag[7]*x1 + diag[11]*x2 + diag[15]*x3;
90: diag += 16;
91: }
92: VecRestoreArray(x,&xx);
93: VecRestoreArray(y,&yy);
94: PetscLogFlops(28*m);
95: return(0);
96: }
99: static PetscErrorCode PCApply_PBJacobi_5(PC pc,Vec x,Vec y)
100: {
101: PC_PBJacobi *jac = (PC_PBJacobi*)pc->data;
103: PetscInt i,m = jac->mbs;
104: PetscScalar *diag = jac->diag,x0,x1,x2,x3,x4,*xx,*yy;
105:
107: VecGetArray(x,&xx);
108: VecGetArray(y,&yy);
109: for (i=0; i<m; i++) {
110: x0 = xx[5*i]; x1 = xx[5*i+1]; x2 = xx[5*i+2]; x3 = xx[5*i+3]; x4 = xx[5*i+4];
111: yy[5*i] = diag[0]*x0 + diag[5]*x1 + diag[10]*x2 + diag[15]*x3 + diag[20]*x4;
112: yy[5*i+1] = diag[1]*x0 + diag[6]*x1 + diag[11]*x2 + diag[16]*x3 + diag[21]*x4;
113: yy[5*i+2] = diag[2]*x0 + diag[7]*x1 + diag[12]*x2 + diag[17]*x3 + diag[22]*x4;
114: yy[5*i+3] = diag[3]*x0 + diag[8]*x1 + diag[13]*x2 + diag[18]*x3 + diag[23]*x4;
115: yy[5*i+4] = diag[4]*x0 + diag[9]*x1 + diag[14]*x2 + diag[19]*x3 + diag[24]*x4;
116: diag += 25;
117: }
118: VecRestoreArray(x,&xx);
119: VecRestoreArray(y,&yy);
120: PetscLogFlops(45*m);
121: return(0);
122: }
123: /* -------------------------------------------------------------------------- */
126: static PetscErrorCode PCSetUp_PBJacobi(PC pc)
127: {
128: PC_PBJacobi *jac = (PC_PBJacobi*)pc->data;
130: PetscMPIInt size;
131: PetscTruth seqbaij,mpibaij,baij;
132: Mat A = pc->pmat;
133: Mat_SeqBAIJ *a;
136: PetscTypeCompare((PetscObject)pc->pmat,MATSEQBAIJ,&seqbaij);
137: PetscTypeCompare((PetscObject)pc->pmat,MATMPIBAIJ,&mpibaij);
138: PetscTypeCompare((PetscObject)pc->pmat,MATBAIJ,&baij);
139: if (!seqbaij && !mpibaij && !baij) {
140: SETERRQ(PETSC_ERR_SUP,"Currently only supports BAIJ matrices");
141: }
142: MPI_Comm_size(pc->comm,&size);
143: if (mpibaij || (baij && (size > 1))) A = ((Mat_MPIBAIJ*)A->data)->A;
144: if (A->rmap.n != A->cmap.n) SETERRQ(PETSC_ERR_SUP,"Supported only for square matrices and square storage");
146: MatSeqBAIJInvertBlockDiagonal(A);
147: a = (Mat_SeqBAIJ*)A->data;
148: jac->diag = a->idiag;
149: jac->bs = A->rmap.bs;
150: jac->mbs = a->mbs;
151: switch (jac->bs){
152: case 2:
153: pc->ops->apply = PCApply_PBJacobi_2;
154: break;
155: case 3:
156: pc->ops->apply = PCApply_PBJacobi_3;
157: break;
158: case 4:
159: pc->ops->apply = PCApply_PBJacobi_4;
160: break;
161: case 5:
162: pc->ops->apply = PCApply_PBJacobi_5;
163: break;
164: default:
165: SETERRQ1(PETSC_ERR_SUP,"not supported for block size %D",jac->bs);
166: }
168: return(0);
169: }
170: /* -------------------------------------------------------------------------- */
173: static PetscErrorCode PCDestroy_PBJacobi(PC pc)
174: {
175: PC_PBJacobi *jac = (PC_PBJacobi*)pc->data;
179: /*
180: Free the private data structure that was hanging off the PC
181: */
182: PetscFree(jac);
183: return(0);
184: }
185: /* -------------------------------------------------------------------------- */
186: /*MC
187: PCPBJACOBI - Point block Jacobi
189: Level: beginner
191: Concepts: point block Jacobi
193: Notes: Only implemented for the BAIJ matrix formats.
195: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC
197: M*/
202: PetscErrorCode PCCreate_PBJacobi(PC pc)
203: {
204: PC_PBJacobi *jac;
209: /*
210: Creates the private data structure for this preconditioner and
211: attach it to the PC object.
212: */
213: PetscNew(PC_PBJacobi,&jac);
214: pc->data = (void*)jac;
216: /*
217: Logs the memory usage; this is not needed but allows PETSc to
218: monitor how much memory is being used for various purposes.
219: */
220: PetscLogObjectMemory(pc,sizeof(PC_PBJacobi));
222: /*
223: Initialize the pointers to vectors to ZERO; these will be used to store
224: diagonal entries of the matrix for fast preconditioner application.
225: */
226: jac->diag = 0;
228: /*
229: Set the pointers for the functions that are provided above.
230: Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
231: are called, they will automatically call these functions. Note we
232: choose not to provide a couple of these functions since they are
233: not needed.
234: */
235: pc->ops->apply = 0; /*set depending on the block size */
236: pc->ops->applytranspose = 0;
237: pc->ops->setup = PCSetUp_PBJacobi;
238: pc->ops->destroy = PCDestroy_PBJacobi;
239: pc->ops->setfromoptions = 0;
240: pc->ops->view = 0;
241: pc->ops->applyrichardson = 0;
242: pc->ops->applysymmetricleft = 0;
243: pc->ops->applysymmetricright = 0;
244: return(0);
245: }