1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38 package ffx.numerics.math;
39
40 import java.util.Random;
41
42 import static org.apache.commons.math3.util.FastMath.floor;
43
44
45
46
47
48
49
50 public class SSETest {
51
52
53
54
55 public final double[][] A;
56
57
58
59 public final double[] x;
60
61
62
63 private final double[] flatA;
64
65
66
67
68
69
70
71 private SSETest(int m, int n) {
72 A = new double[m][n];
73 flatA = new double[m * n];
74 x = new double[n];
75 }
76
77
78
79
80
81
82 public static void main(String[] args) {
83
84 int m = 500;
85 int n = 500;
86 if (args != null && args.length > 0) {
87 m = Integer.parseInt(args[0]);
88 n = Integer.parseInt(args[1]);
89 }
90
91 SSETest test = new SSETest(m, n);
92
93 int nLoops = 100000;
94
95 test.init(m, n);
96 Random r = new Random(0);
97 double temp = 0;
98 long time = 0;
99 for (int i = 1; i <= nLoops; i++) {
100 time -= System.nanoTime();
101 double[] y = test.matVec(test.A, test.x, m, n);
102 time += System.nanoTime();
103 int j = (int) floor(m * r.nextDouble());
104 temp += y[j];
105 if (i % 10000 == 0) {
106 System.out.println(" Nested: " + temp + " " + time * 1.0e-9);
107 time = 0;
108 }
109 }
110
111 test.init(m, n);
112 r = new Random(0);
113 temp = 0;
114 time = 0;
115 for (int i = 1; i <= nLoops; i++) {
116 time -= System.nanoTime();
117 double[] y = test.matVec(test.flatA, test.x, m, n);
118 time += System.nanoTime();
119 int j = (int) floor(m * r.nextDouble());
120 temp += y[j];
121 if (i % 10000 == 0) {
122 System.out.println(" Flat: " + temp + " " + time * 1.0e-9);
123 time = 0;
124 }
125 }
126 }
127
128 private void init(int n, int m) {
129 Random r = new Random(0);
130
131
132 for (int i = 0; i < m; i++) {
133 int idx = i * n;
134 for (int j = 0; j < n; j++) {
135 A[i][j] = r.nextDouble();
136 flatA[idx + j] = A[i][j];
137 }
138 }
139
140 for (int j = 0; j < n; j++) {
141 x[j] = r.nextDouble();
142 }
143 }
144
145
146
147
148
149
150
151
152
153
154 private double[] matVec(final double[][] A, final double[] x, final int m, final int n) {
155 double[] y = new double[m];
156 for (int i = 0; i < m; i++) {
157 for (int j = 0; j < n; j++) {
158
159 y[i] += A[i][j] * x[j];
160 }
161 }
162 return y;
163 }
164
165
166
167
168
169
170
171
172
173
174 private double[] matVec(final double[] A, final double[] x, final int m, final int n) {
175 double[] y = new double[m];
176 final int extra = n - n % 8;
177 final int ub = ((n / 8) * 8) - 1;
178 int idx = 0;
179 for (int i = 0; i < m; i++) {
180 double acc = 0;
181 for (int j = 0; j < ub; j += 8) {
182 int ptr = idx + j;
183 y[i] +=
184 A[ptr] * x[j]
185 + A[ptr + 1] * x[j + 1]
186 + A[ptr + 2] * x[j + 2]
187 + A[ptr + 3] * x[j + 3];
188 acc += A[ptr + 4] * x[j + 4]
189 + A[ptr + 5] * x[j + 5]
190 + A[ptr + 6] * x[j + 6]
191 + A[ptr + 7] * x[j + 7];
192 }
193 y[i] += acc;
194 for (int j = extra; j < n; j++) {
195 y[i] += A[idx + j] * x[j];
196 }
197 idx += n;
198 }
199 return y;
200 }
201 }