1 // ******************************************************************************
2 //
3 // Title: Force Field X.
4 // Description: Force Field X - Software for Molecular Biophysics.
5 // Copyright: Copyright (c) Michael J. Schnieders 2001-2025.
6 //
7 // This file is part of Force Field X.
8 //
9 // Force Field X is free software; you can redistribute it and/or modify it
10 // under the terms of the GNU General Public License version 3 as published by
11 // the Free Software Foundation.
12 //
13 // Force Field X is distributed in the hope that it will be useful, but WITHOUT
14 // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
15 // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
16 // details.
17 //
18 // You should have received a copy of the GNU General Public License along with
19 // Force Field X; if not, write to the Free Software Foundation, Inc., 59 Temple
20 // Place, Suite 330, Boston, MA 02111-1307 USA
21 //
22 // Linking this library statically or dynamically with other modules is making a
23 // combined work based on this library. Thus, the terms and conditions of the
24 // GNU General Public License cover the whole combination.
25 //
26 // As a special exception, the copyright holders of this library give you
27 // permission to link this library with independent modules to produce an
28 // executable, regardless of the license terms of these independent modules, and
29 // to copy and distribute the resulting executable under terms of your choice,
30 // provided that you also meet, for each linked independent module, the terms
31 // and conditions of the license of that module. An independent module is a
32 // module which is not derived from or based on this library. If you modify this
33 // library, you may extend this exception to your version of the library, but
34 // you are not obligated to do so. If you do not wish to do so, delete this
35 // exception statement from your version.
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 * java -cp target/numerics-1.0.0.jar -XX:+UnlockDiagnosticVMOptions -XX:+PrintAssembly
46 * -Djava.library.path=hsdis-amd64.dylib ffx.numerics.math.SSETest
47 *
48 * @author M. J. Schnieders
49 */
50 public class SSETest {
51
52 /**
53 * A matrix of double values.
54 */
55 public final double[][] A;
56 /**
57 * A vector of double values.
58 */
59 public final double[] x;
60 /**
61 * A flattened matrix of double values.
62 */
63 private final double[] flatA;
64
65 /**
66 * Constructor for SSETest.
67 *
68 * @param m the number of rows.
69 * @param n the number of columns.
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 * main.
79 *
80 * @param args an array of {@link java.lang.String} objects.
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 // Initialize A and flatA.
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 // Initialize X.
140 for (int j = 0; j < n; j++) {
141 x[j] = r.nextDouble();
142 }
143 }
144
145 /**
146 * matVec.
147 *
148 * @param A an array of double values.
149 * @param x an array of double values.
150 * @param m the number of rows.
151 * @param n the number of columns.
152 * @return an array of double values.
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 // y[i] = fma(A[i][j], x[j], y[i]);
159 y[i] += A[i][j] * x[j];
160 }
161 }
162 return y;
163 }
164
165 /**
166 * matVec.
167 *
168 * @param A an array of double values.
169 * @param x an array of double values.
170 * @param m the number of rows.
171 * @param n the number of columns.
172 * @return an array of double values.
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 }