1 /**************************************************************************
3 * Java Grande Forum Benchmark Suite - Thread Version 1.0 *
7 * Java Grande Benchmarking Project *
11 * Edinburgh Parallel Computing Centre *
13 * email: epcc-javagrande@epcc.ed.ac.uk *
15 * See below for the previous history of this code *
17 * This version copyright (c) The University of Edinburgh, 2001. *
18 * All rights reserved. *
20 **************************************************************************/
24 Modified 3/3/97 by David M. Doolin (dmd) doolin@cs.utk.edu
25 Fixed error in matgen() method. Added some comments.
27 Modified 1/22/97 by Paul McMahan mcmahan@cs.utk.edu
28 Added more MacOS options to form.
30 Optimized by Jonathan Hardwick (jch@cs.cmu.edu), 3/28/96
31 Compare to Linkpack.java.
32 Optimizations performed:
33 - added "final" modifier to performance-critical methods.
34 - changed lines of the form "a[i] = a[i] + x" to "a[i] += x".
35 - minimized array references using common subexpression elimination.
36 - eliminated unused variables.
37 - undid an unrolled loop.
38 - added temporary 1D arrays to hold frequently-used columns of 2D arrays.
39 - wrote my own abs() method
40 See http://www.cs.cmu.edu/~jch/java/linpack.html for more details.
43 Ported to Java by Reed Wade (wade@cs.utk.edu) 2/96
44 built using JDK 1.0 on solaris
45 using "javac -O Linpack.java"
48 Translated to C by Bonnie Toy 5/88
49 (modified on 2/25/94 to fix a problem with daxpy for
50 unequal increments or equal increments not equal to 1.
55 public class Linpack {
59 double ops,total,norma,normx;
62 int n,i,ntimes,info,lda,ldaa,kflops;
65 double abs (double d) {
68 //return (d >= 0) ? d : -d;
71 double matgen (double a[][], int lda, int n, double b[])
78 /* Next two for() statements switched. Solver wants
79 matrix in column order. --dmd 3/3/97
81 for (i = 0; i < n; i++) {
82 for (j = 0; j < n; j++) {
83 init = 3125*init % 65536;
84 a[j][i] = (init - 32768.0)/16384.0;
85 //norma = (a[j][i] > norma) ? a[j][i] : norma;
86 if (a[j][i] > norma) {
91 for (i = 0; i < n; i++) {
94 for (j = 0; j < n; j++) {
95 for (i = 0; i < n; i++) {
102 void dgesl( double a[][], int lda, int n, int ipvt[], double b[], int job)
109 // job = 0 , solve a * x = b. first solve l*y = b
111 for (k = 0; k < nm1; k++) {
119 daxpy(n-(kp1),t,a[k],kp1,1,b,kp1,1);
123 for (kb = 0; kb < n; kb++) {
127 daxpy(k,t,a[k],0,1,b,0,1);
130 // job = nonzero, solve trans(a) * x = b. first solve trans(u)*y = b
131 for (k = 0; k < n; k++) {
132 t = ddot(k,a[k],0,1,b,0,1);
133 b[k] = (b[k] - t)/a[k][k];
135 // now solve trans(l)*x = y
137 for (kb = 1; kb < nm1; kb++) {
140 b[k] += ddot(n-(kp1),a[k],kp1,1,b,kp1,1);
153 constant times a vector plus a vector.
154 jack dongarra, linpack, 3/11/78.
156 void daxpy( int n, double da, double dx[], int dx_off, int incx,
157 double dy[], int dy_off, int incy)
161 if ((n > 0) && (da != 0)) {
162 if (incx != 1 || incy != 1) {
163 // code for unequal increments or equal increments not equal to 1
166 if (incx < 0) ix = (-n+1)*incx;
167 if (incy < 0) iy = (-n+1)*incy;
168 for (i = 0;i < n; i++) {
169 dy[iy +dy_off] += da*dx[ix +dx_off];
175 // code for both increments equal to 1
176 for (i=0; i < n; i++)
177 dy[i +dy_off] += da*dx[i +dx_off];
183 forms the dot product of two vectors.
184 jack dongarra, linpack, 3/11/78.
186 double ddot( int n, double dx[], int dx_off, int incx, double dy[],
187 int dy_off, int incy)
193 if (incx != 1 || incy != 1) {
194 // code for unequal increments or equal increments not equal to 1
197 if (incx < 0) ix = (-n+1)*incx;
198 if (incy < 0) iy = (-n+1)*incy;
199 for (i = 0;i < n; i++) {
200 dtemp += dx[ix +dx_off]*dy[iy +dy_off];
205 // code for both increments equal to 1
207 dtemp += dx[i +dx_off]*dy[i +dy_off];
214 scales a vector by a constant.
215 jack dongarra, linpack, 3/11/78.
217 void dscal( int n, double da, double dx[], int dx_off, int incx)
222 // code for increment not equal to 1
224 for (i = 0; i < nincx; i += incx)
227 // code for increment equal to 1
228 for (i = 0; i < n; i++)
235 finds the index of element having max. absolute value.
236 jack dongarra, linpack, 3/11/78.
238 int idamax( int n, double dx[], int dx_off, int incx)
247 } else if (incx != 1) {
248 // code for increment not equal to 1
249 dmax = abs(dx[0 +dx_off]);
251 for (i = 1; i < n; i++) {
252 dtemp = abs(dx[ix + dx_off]);
260 // code for increment equal to 1
262 dmax = abs(dx[0 +dx_off]);
263 for (i = 1; i < n; i++) {
264 dtemp = abs(dx[i + dx_off]);
275 estimate unit roundoff in quantities of size x.
276 this program should function properly on all systems
277 satisfying the following two assumptions,
278 1. the base used in representing dfloating point
279 numbers is not a power of three.
280 2. the quantity a in statement 10 is represented to
281 the accuracy used in dfloating point variables
282 that are stored in memory.
283 the statement number 10 and the go to 10 are intended to
284 force optimizing compilers to generate code satisfying
286 under these assumptions, it should be true that,
287 a is not exactly equal to four-thirds,
288 b has a zero for its last bit or digit,
289 c is not exactly equal to one,
290 eps measures the separation of 1.0 from
291 the next larger dfloating point number.
292 the developers of eispack would appreciate being informed
293 about any systems where these assumptions do not hold.
295 *****************************************************************
296 this routine is one of the auxiliary routines used by eispack iii
297 to avoid machine dependencies.
298 *****************************************************************
300 this version dated 4/6/83.
302 double epslon (double x)
316 void dmxpy ( int n1, double y[], int n2, int ldm, double x[], double m[][])
319 // cleanup odd vector
320 for (j = 0; j < n2; j++) {
321 for (i = 0; i < n1; i++) {
322 y[i] += x[j]*m[j][i];