HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
matmult_native.h
Go to the documentation of this file.
1 
8 #ifndef _MATMULT_H_
9 #define _MATMULT_H_
10 
11 /* The following macro names conflict with PETSc stuff */
12 #ifndef _PETSC_INTERFACE_H_
13 
18 #define MatMult(N,A,X,Y) \
19  { \
20  int i,j,k; \
21  for (i = 0; i < (N); i++) { \
22  for (j = 0; j < (N); j++) { \
23  A[i*(N)+j] = 0.0; \
24  for (k = 0; k < (N); k++) A[i*(N)+j] += X[i*(N)+k]*Y[k*(N)+j]; \
25  } \
26  } \
27  }
28 
32 #define MatMult2(N,A,X,Y) \
33  { \
34  A[0] = X[0]*Y[0] + X[1]*Y[2]; \
35  A[1] = X[0]*Y[1] + X[1]*Y[3]; \
36  A[2] = X[2]*Y[0] + X[3]*Y[2]; \
37  A[3] = X[2]*Y[1] + X[3]*Y[3]; \
38  }
39 
43 #define MatMult3(N,A,X,Y) \
44  { \
45  A[0] = X[0]*Y[0] + X[1]*Y[3] + X[2]*Y[6]; \
46  A[1] = X[0]*Y[1] + X[1]*Y[4] + X[2]*Y[7]; \
47  A[2] = X[0]*Y[2] + X[1]*Y[5] + X[2]*Y[8]; \
48  A[3] = X[3]*Y[0] + X[4]*Y[3] + X[5]*Y[6]; \
49  A[4] = X[3]*Y[1] + X[4]*Y[4] + X[5]*Y[7]; \
50  A[5] = X[3]*Y[2] + X[4]*Y[5] + X[5]*Y[8]; \
51  A[6] = X[6]*Y[0] + X[7]*Y[3] + X[8]*Y[6]; \
52  A[7] = X[6]*Y[1] + X[7]*Y[4] + X[8]*Y[7]; \
53  A[8] = X[6]*Y[2] + X[7]*Y[5] + X[8]*Y[8]; \
54  }
55 
59 #define MatMult4(N,A,X,Y) \
60  { \
61  A[0] = X[0]*Y[0] + X[1]*Y[4] + X[2]*Y[8] + X[3]*Y[12]; \
62  A[1] = X[0]*Y[1] + X[1]*Y[5] + X[2]*Y[9] + X[3]*Y[13]; \
63  A[2] = X[0]*Y[2] + X[1]*Y[6] + X[2]*Y[10] + X[3]*Y[14]; \
64  A[3] = X[0]*Y[3] + X[1]*Y[7] + X[2]*Y[11] + X[3]*Y[15]; \
65  A[4] = X[4]*Y[0] + X[5]*Y[4] + X[6]*Y[8] + X[7]*Y[12]; \
66  A[5] = X[4]*Y[1] + X[5]*Y[5] + X[6]*Y[9] + X[7]*Y[13]; \
67  A[6] = X[4]*Y[2] + X[5]*Y[6] + X[6]*Y[10] + X[7]*Y[14]; \
68  A[7] = X[4]*Y[3] + X[5]*Y[7] + X[6]*Y[11] + X[7]*Y[15]; \
69  A[8] = X[8]*Y[0] + X[9]*Y[4] + X[10]*Y[8] + X[11]*Y[12];\
70  A[9] = X[8]*Y[1] + X[9]*Y[5] + X[10]*Y[9] + X[11]*Y[13];\
71  A[10] = X[8]*Y[2] + X[9]*Y[6] + X[10]*Y[10] + X[11]*Y[14];\
72  A[11] = X[8]*Y[3] + X[9]*Y[7] + X[10]*Y[11] + X[11]*Y[15];\
73  A[12] = X[12]*Y[0] + X[13]*Y[4] + X[14]*Y[8] + X[15]*Y[12];\
74  A[13] = X[12]*Y[1] + X[13]*Y[5] + X[14]*Y[9] + X[15]*Y[13];\
75  A[14] = X[12]*Y[2] + X[13]*Y[6] + X[14]*Y[10] + X[15]*Y[14];\
76  A[15] = X[12]*Y[3] + X[13]*Y[7] + X[14]*Y[11] + X[15]*Y[15];\
77  }
78 
82 #define MatMult5(N,A,X,Y) \
83  { \
84  A[0] = X[0]*Y[0] + X[1]*Y[5] + X[2]*Y[10] + X[3]*Y[15] + X[4]*Y[20]; \
85  A[1] = X[0]*Y[1] + X[1]*Y[6] + X[2]*Y[11] + X[3]*Y[16] + X[4]*Y[21]; \
86  A[2] = X[0]*Y[2] + X[1]*Y[7] + X[2]*Y[12] + X[3]*Y[17] + X[4]*Y[22]; \
87  A[3] = X[0]*Y[3] + X[1]*Y[8] + X[2]*Y[13] + X[3]*Y[18] + X[4]*Y[23]; \
88  A[4] = X[0]*Y[4] + X[1]*Y[9] + X[2]*Y[14] + X[3]*Y[19] + X[4]*Y[24]; \
89  A[5] = X[5]*Y[0] + X[6]*Y[5] + X[7]*Y[10] + X[8]*Y[15] + X[9]*Y[20]; \
90  A[6] = X[5]*Y[1] + X[6]*Y[6] + X[7]*Y[11] + X[8]*Y[16] + X[9]*Y[21]; \
91  A[7] = X[5]*Y[2] + X[6]*Y[7] + X[7]*Y[12] + X[8]*Y[17] + X[9]*Y[22]; \
92  A[8] = X[5]*Y[3] + X[6]*Y[8] + X[7]*Y[13] + X[8]*Y[18] + X[9]*Y[23]; \
93  A[9] = X[5]*Y[4] + X[6]*Y[9] + X[7]*Y[14] + X[8]*Y[19] + X[9]*Y[24]; \
94  A[10] = X[10]*Y[0] + X[11]*Y[5] + X[12]*Y[10] + X[13]*Y[15] + X[14]*Y[20];\
95  A[11] = X[10]*Y[1] + X[11]*Y[6] + X[12]*Y[11] + X[13]*Y[16] + X[14]*Y[21];\
96  A[12] = X[10]*Y[2] + X[11]*Y[7] + X[12]*Y[12] + X[13]*Y[17] + X[14]*Y[22];\
97  A[13] = X[10]*Y[3] + X[11]*Y[8] + X[12]*Y[13] + X[13]*Y[18] + X[14]*Y[23];\
98  A[14] = X[10]*Y[4] + X[11]*Y[9] + X[12]*Y[14] + X[13]*Y[19] + X[14]*Y[24];\
99  A[15] = X[15]*Y[0] + X[16]*Y[5] + X[17]*Y[10] + X[18]*Y[15] + X[19]*Y[20];\
100  A[16] = X[15]*Y[1] + X[16]*Y[6] + X[17]*Y[11] + X[18]*Y[16] + X[19]*Y[21];\
101  A[17] = X[15]*Y[2] + X[16]*Y[7] + X[17]*Y[12] + X[18]*Y[17] + X[19]*Y[22];\
102  A[18] = X[15]*Y[3] + X[16]*Y[8] + X[17]*Y[13] + X[18]*Y[18] + X[19]*Y[23];\
103  A[19] = X[15]*Y[4] + X[16]*Y[9] + X[17]*Y[14] + X[18]*Y[19] + X[19]*Y[24];\
104  A[20] = X[20]*Y[0] + X[21]*Y[5] + X[22]*Y[10] + X[23]*Y[15] + X[24]*Y[20];\
105  A[21] = X[20]*Y[1] + X[21]*Y[6] + X[22]*Y[11] + X[23]*Y[16] + X[24]*Y[21];\
106  A[22] = X[20]*Y[2] + X[21]*Y[7] + X[22]*Y[12] + X[23]*Y[17] + X[24]*Y[22];\
107  A[23] = X[20]*Y[3] + X[21]*Y[8] + X[22]*Y[13] + X[23]*Y[18] + X[24]*Y[23];\
108  A[24] = X[20]*Y[4] + X[21]*Y[9] + X[22]*Y[14] + X[23]*Y[19] + X[24]*Y[24];\
109  }
110 
116 #define MatVecMult(N,y,A,x) \
117  { \
118  int i,j; \
119  for (i = 0; i < (N); i++) { \
120  y[i] = 0; \
121  for (j = 0; j < (N); j++) y[i] += A[i*(N)+j]*x[j]; \
122  } \
123  }
124 
128 #define MatVecMult2(N,y,A,x) \
129  { \
130  y[0] = A[0]*x[0] + A[1]*x[1];\
131  y[1] = A[2]*x[0] + A[3]*x[1];\
132  }
133 
137 #define MatVecMult3(N,y,A,x) \
138  { \
139  y[0] = A[0]*x[0] + A[1]*x[1] + A[2]*x[2];\
140  y[1] = A[3]*x[0] + A[4]*x[1] + A[5]*x[2];\
141  y[2] = A[6]*x[0] + A[7]*x[1] + A[8]*x[2];\
142  }
143 
147 #define MatVecMult4(N,y,A,x) \
148  { \
149  y[0] = A[0]*x[0] + A[1]*x[1] + A[2]*x[2] + A[3]*x[3]; \
150  y[1] = A[4]*x[0] + A[5]*x[1] + A[6]*x[2] + A[7]*x[3]; \
151  y[2] = A[8]*x[0] + A[9]*x[1] + A[10]*x[2] + A[11]*x[3];\
152  y[3] = A[12]*x[0] + A[13]*x[1] + A[14]*x[2] + A[15]*x[3];\
153  }
154 
158 #define MatVecMult5(N,y,A,x) \
159  { \
160  y[0] = A[0]*x[0] + A[1]*x[1] + A[2]*x[2] + A[3]*x[3] + A[4]*x[4]; \
161  y[1] = A[5]*x[0] + A[6]*x[1] + A[7]*x[2] + A[8]*x[3] + A[9]*x[4]; \
162  y[2] = A[10]*x[0] + A[11]*x[1] + A[12]*x[2] + A[13]*x[3] + A[14]*x[4];\
163  y[3] = A[15]*x[0] + A[16]*x[1] + A[17]*x[2] + A[18]*x[3] + A[19]*x[4];\
164  y[4] = A[20]*x[0] + A[21]*x[1] + A[22]*x[2] + A[23]*x[3] + A[24]*x[4];\
165  }
166 
167 #endif
168 
169 #endif