Actual source code: pbjacobi.c
petsc-3.5.4 2015-05-23
2: /*
3: Include files needed for the PBJacobi preconditioner:
4: pcimpl.h - private include file intended for use by all preconditioners
5: */
7: #include <petsc-private/matimpl.h>
8: #include <petsc-private/pcimpl.h> /*I "petscpc.h" I*/
10: /*
11: Private context (data structure) for the PBJacobi preconditioner.
12: */
13: typedef struct {
14: const MatScalar *diag;
15: PetscInt bs,mbs;
16: } PC_PBJacobi;
21: static PetscErrorCode PCApply_PBJacobi_1(PC pc,Vec x,Vec y)
22: {
23: PC_PBJacobi *jac = (PC_PBJacobi*)pc->data;
24: PetscErrorCode ierr;
25: PetscInt i,m = jac->mbs;
26: const MatScalar *diag = jac->diag;
27: const PetscScalar *xx;
28: PetscScalar *yy;
31: VecGetArrayRead(x,&xx);
32: VecGetArray(y,&yy);
33: for (i=0; i<m; i++) yy[i] = diag[i]*xx[i];
34: VecRestoreArrayRead(x,&xx);
35: VecRestoreArray(y,&yy);
36: PetscLogFlops(2.0*m);
37: return(0);
38: }
42: static PetscErrorCode PCApply_PBJacobi_2(PC pc,Vec x,Vec y)
43: {
44: PC_PBJacobi *jac = (PC_PBJacobi*)pc->data;
45: PetscErrorCode ierr;
46: PetscInt i,m = jac->mbs;
47: const MatScalar *diag = jac->diag;
48: PetscScalar x0,x1,*xx,*yy;
51: VecGetArray(x,&xx);
52: VecGetArray(y,&yy);
53: for (i=0; i<m; i++) {
54: x0 = xx[2*i]; x1 = xx[2*i+1];
55: yy[2*i] = diag[0]*x0 + diag[2]*x1;
56: yy[2*i+1] = diag[1]*x0 + diag[3]*x1;
57: diag += 4;
58: }
59: VecRestoreArray(x,&xx);
60: VecRestoreArray(y,&yy);
61: PetscLogFlops(6.0*m);
62: return(0);
63: }
66: static PetscErrorCode PCApply_PBJacobi_3(PC pc,Vec x,Vec y)
67: {
68: PC_PBJacobi *jac = (PC_PBJacobi*)pc->data;
69: PetscErrorCode ierr;
70: PetscInt i,m = jac->mbs;
71: const MatScalar *diag = jac->diag;
72: PetscScalar x0,x1,x2,*xx,*yy;
75: VecGetArray(x,&xx);
76: VecGetArray(y,&yy);
77: for (i=0; i<m; i++) {
78: x0 = xx[3*i]; x1 = xx[3*i+1]; x2 = xx[3*i+2];
80: yy[3*i] = diag[0]*x0 + diag[3]*x1 + diag[6]*x2;
81: yy[3*i+1] = diag[1]*x0 + diag[4]*x1 + diag[7]*x2;
82: yy[3*i+2] = diag[2]*x0 + diag[5]*x1 + diag[8]*x2;
83: diag += 9;
84: }
85: VecRestoreArray(x,&xx);
86: VecRestoreArray(y,&yy);
87: PetscLogFlops(15.0*m);
88: return(0);
89: }
92: static PetscErrorCode PCApply_PBJacobi_4(PC pc,Vec x,Vec y)
93: {
94: PC_PBJacobi *jac = (PC_PBJacobi*)pc->data;
95: PetscErrorCode ierr;
96: PetscInt i,m = jac->mbs;
97: const MatScalar *diag = jac->diag;
98: PetscScalar x0,x1,x2,x3,*xx,*yy;
101: VecGetArray(x,&xx);
102: VecGetArray(y,&yy);
103: for (i=0; i<m; i++) {
104: x0 = xx[4*i]; x1 = xx[4*i+1]; x2 = xx[4*i+2]; x3 = xx[4*i+3];
106: yy[4*i] = diag[0]*x0 + diag[4]*x1 + diag[8]*x2 + diag[12]*x3;
107: yy[4*i+1] = diag[1]*x0 + diag[5]*x1 + diag[9]*x2 + diag[13]*x3;
108: yy[4*i+2] = diag[2]*x0 + diag[6]*x1 + diag[10]*x2 + diag[14]*x3;
109: yy[4*i+3] = diag[3]*x0 + diag[7]*x1 + diag[11]*x2 + diag[15]*x3;
110: diag += 16;
111: }
112: VecRestoreArray(x,&xx);
113: VecRestoreArray(y,&yy);
114: PetscLogFlops(28.0*m);
115: return(0);
116: }
119: static PetscErrorCode PCApply_PBJacobi_5(PC pc,Vec x,Vec y)
120: {
121: PC_PBJacobi *jac = (PC_PBJacobi*)pc->data;
122: PetscErrorCode ierr;
123: PetscInt i,m = jac->mbs;
124: const MatScalar *diag = jac->diag;
125: PetscScalar x0,x1,x2,x3,x4,*xx,*yy;
128: VecGetArray(x,&xx);
129: VecGetArray(y,&yy);
130: for (i=0; i<m; i++) {
131: x0 = xx[5*i]; x1 = xx[5*i+1]; x2 = xx[5*i+2]; x3 = xx[5*i+3]; x4 = xx[5*i+4];
133: yy[5*i] = diag[0]*x0 + diag[5]*x1 + diag[10]*x2 + diag[15]*x3 + diag[20]*x4;
134: yy[5*i+1] = diag[1]*x0 + diag[6]*x1 + diag[11]*x2 + diag[16]*x3 + diag[21]*x4;
135: yy[5*i+2] = diag[2]*x0 + diag[7]*x1 + diag[12]*x2 + diag[17]*x3 + diag[22]*x4;
136: yy[5*i+3] = diag[3]*x0 + diag[8]*x1 + diag[13]*x2 + diag[18]*x3 + diag[23]*x4;
137: yy[5*i+4] = diag[4]*x0 + diag[9]*x1 + diag[14]*x2 + diag[19]*x3 + diag[24]*x4;
138: diag += 25;
139: }
140: VecRestoreArray(x,&xx);
141: VecRestoreArray(y,&yy);
142: PetscLogFlops(45.0*m);
143: return(0);
144: }
147: static PetscErrorCode PCApply_PBJacobi_6(PC pc,Vec x,Vec y)
148: {
149: PC_PBJacobi *jac = (PC_PBJacobi*)pc->data;
150: PetscErrorCode ierr;
151: PetscInt i,m = jac->mbs;
152: const MatScalar *diag = jac->diag;
153: PetscScalar x0,x1,x2,x3,x4,x5,*xx,*yy;
156: VecGetArray(x,&xx);
157: VecGetArray(y,&yy);
158: for (i=0; i<m; i++) {
159: x0 = xx[6*i]; x1 = xx[6*i+1]; x2 = xx[6*i+2]; x3 = xx[6*i+3]; x4 = xx[6*i+4]; x5 = xx[6*i+5];
161: yy[6*i] = diag[0]*x0 + diag[6]*x1 + diag[12]*x2 + diag[18]*x3 + diag[24]*x4 + diag[30]*x5;
162: yy[6*i+1] = diag[1]*x0 + diag[7]*x1 + diag[13]*x2 + diag[19]*x3 + diag[25]*x4 + diag[31]*x5;
163: yy[6*i+2] = diag[2]*x0 + diag[8]*x1 + diag[14]*x2 + diag[20]*x3 + diag[26]*x4 + diag[32]*x5;
164: yy[6*i+3] = diag[3]*x0 + diag[9]*x1 + diag[15]*x2 + diag[21]*x3 + diag[27]*x4 + diag[33]*x5;
165: yy[6*i+4] = diag[4]*x0 + diag[10]*x1 + diag[16]*x2 + diag[22]*x3 + diag[28]*x4 + diag[34]*x5;
166: yy[6*i+5] = diag[5]*x0 + diag[11]*x1 + diag[17]*x2 + diag[23]*x3 + diag[29]*x4 + diag[35]*x5;
167: diag += 36;
168: }
169: VecRestoreArray(x,&xx);
170: VecRestoreArray(y,&yy);
171: PetscLogFlops(66.0*m);
172: return(0);
173: }
176: static PetscErrorCode PCApply_PBJacobi_7(PC pc,Vec x,Vec y)
177: {
178: PC_PBJacobi *jac = (PC_PBJacobi*)pc->data;
179: PetscErrorCode ierr;
180: PetscInt i,m = jac->mbs;
181: const MatScalar *diag = jac->diag;
182: PetscScalar x0,x1,x2,x3,x4,x5,x6,*xx,*yy;
185: VecGetArray(x,&xx);
186: VecGetArray(y,&yy);
187: for (i=0; i<m; i++) {
188: x0 = xx[7*i]; x1 = xx[7*i+1]; x2 = xx[7*i+2]; x3 = xx[7*i+3]; x4 = xx[7*i+4]; x5 = xx[7*i+5]; x6 = xx[7*i+6];
190: yy[7*i] = diag[0]*x0 + diag[7]*x1 + diag[14]*x2 + diag[21]*x3 + diag[28]*x4 + diag[35]*x5 + diag[42]*x6;
191: yy[7*i+1] = diag[1]*x0 + diag[8]*x1 + diag[15]*x2 + diag[22]*x3 + diag[29]*x4 + diag[36]*x5 + diag[43]*x6;
192: yy[7*i+2] = diag[2]*x0 + diag[9]*x1 + diag[16]*x2 + diag[23]*x3 + diag[30]*x4 + diag[37]*x5 + diag[44]*x6;
193: yy[7*i+3] = diag[3]*x0 + diag[10]*x1 + diag[17]*x2 + diag[24]*x3 + diag[31]*x4 + diag[38]*x5 + diag[45]*x6;
194: yy[7*i+4] = diag[4]*x0 + diag[11]*x1 + diag[18]*x2 + diag[25]*x3 + diag[32]*x4 + diag[39]*x5 + diag[46]*x6;
195: yy[7*i+5] = diag[5]*x0 + diag[12]*x1 + diag[19]*x2 + diag[26]*x3 + diag[33]*x4 + diag[40]*x5 + diag[47]*x6;
196: yy[7*i+6] = diag[6]*x0 + diag[13]*x1 + diag[20]*x2 + diag[27]*x3 + diag[34]*x4 + diag[41]*x5 + diag[48]*x6;
197: diag += 49;
198: }
199: VecRestoreArray(x,&xx);
200: VecRestoreArray(y,&yy);
201: PetscLogFlops(80.0*m);
202: return(0);
203: }
204: /* -------------------------------------------------------------------------- */
207: static PetscErrorCode PCSetUp_PBJacobi(PC pc)
208: {
209: PC_PBJacobi *jac = (PC_PBJacobi*)pc->data;
211: Mat A = pc->pmat;
214: if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Supported only for square matrices and square storage");
216: MatInvertBlockDiagonal(A,&jac->diag);
217: MatGetBlockSize(A,&jac->bs);
218: jac->mbs = A->rmap->n/jac->bs;
219: switch (jac->bs) {
220: case 1:
221: pc->ops->apply = PCApply_PBJacobi_1;
222: break;
223: case 2:
224: pc->ops->apply = PCApply_PBJacobi_2;
225: break;
226: case 3:
227: pc->ops->apply = PCApply_PBJacobi_3;
228: break;
229: case 4:
230: pc->ops->apply = PCApply_PBJacobi_4;
231: break;
232: case 5:
233: pc->ops->apply = PCApply_PBJacobi_5;
234: break;
235: case 6:
236: pc->ops->apply = PCApply_PBJacobi_6;
237: break;
238: case 7:
239: pc->ops->apply = PCApply_PBJacobi_7;
240: break;
241: default:
242: SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"not supported for block size %D",jac->bs);
243: }
244: return(0);
245: }
246: /* -------------------------------------------------------------------------- */
249: static PetscErrorCode PCDestroy_PBJacobi(PC pc)
250: {
254: /*
255: Free the private data structure that was hanging off the PC
256: */
257: PetscFree(pc->data);
258: return(0);
259: }
263: static PetscErrorCode PCView_PBJacobi(PC pc,PetscViewer viewer)
264: {
266: PC_PBJacobi *jac = (PC_PBJacobi*)pc->data;
267: PetscBool iascii;
270: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
271: if (iascii) {
272: PetscViewerASCIIPrintf(viewer," point-block Jacobi: block size %D\n",jac->bs);
273: }
274: return(0);
275: }
277: /* -------------------------------------------------------------------------- */
278: /*MC
279: PCPBJACOBI - Point block Jacobi
281: Level: beginner
283: Concepts: point block Jacobi
286: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC
288: M*/
292: PETSC_EXTERN PetscErrorCode PCCreate_PBJacobi(PC pc)
293: {
294: PC_PBJacobi *jac;
298: /*
299: Creates the private data structure for this preconditioner and
300: attach it to the PC object.
301: */
302: PetscNewLog(pc,&jac);
303: pc->data = (void*)jac;
305: /*
306: Initialize the pointers to vectors to ZERO; these will be used to store
307: diagonal entries of the matrix for fast preconditioner application.
308: */
309: jac->diag = 0;
311: /*
312: Set the pointers for the functions that are provided above.
313: Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
314: are called, they will automatically call these functions. Note we
315: choose not to provide a couple of these functions since they are
316: not needed.
317: */
318: pc->ops->apply = 0; /*set depending on the block size */
319: pc->ops->applytranspose = 0;
320: pc->ops->setup = PCSetUp_PBJacobi;
321: pc->ops->destroy = PCDestroy_PBJacobi;
322: pc->ops->setfromoptions = 0;
323: pc->ops->view = PCView_PBJacobi;
324: pc->ops->applyrichardson = 0;
325: pc->ops->applysymmetricleft = 0;
326: pc->ops->applysymmetricright = 0;
327: return(0);
328: }