View Javadoc
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 }