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 *
16 * This version copyright (c) The University of Edinburgh, 2001. *
17 * All rights reserved. *
19 **************************************************************************/
20 public class JGFLUFactBench {
23 private int[] datasizes;
24 public int cachelinesize;
28 double ops,total,norma,normx;
31 int n,i,ntimes,info,lda,ldaa,kflops;
34 public JGFLUFactBench(int nthreads) {
35 this.nthreads=nthreads;
36 datasizes = global new int[3];
43 public void JGFsetsize(int size) {
47 public void JGFinitialise() {
52 a = global new double[ldaa][lda];
53 b = global new double [ldaa];
54 x = global new double [ldaa];
55 ipvt = global new int [ldaa];
57 long nl = (long) n; //avoid integer overflow
58 ops = (2.0*(nl*nl*nl))/3.0 + 2.0*(nl*nl);
59 norma = matgen(a,lda,n,b);
62 public static void JGFkernel(JGFLUFactBench lub) {
65 numthreads = lub.nthreads;
69 LinpackRunner[] thobjects;
72 thobjects = global new LinpackRunner[numthreads];
73 br = global new Barrier(numthreads);
76 //JGFInstrumentor.startTimer("Section2:LUFact:Kernel", instr.timers);
79 int mid = (128<<24)|(195<<16)|(175<<8)|73;
80 for(int i=1;i<numthreads;i++) {
82 thobjects[i] = global new LinpackRunner(i,lub.a,lub.lda,lub.n,lub.ipvt,br,lub.nthreads);
89 thobjects[0] = global new LinpackRunner(0,lub.a,lub.lda,lub.n,lub.ipvt,br,lub.nthreads);
95 for(int i=1;i<numthreads;i++) {
103 lub.dgesl(lub.a,lub.lda,lub.n,lub.ipvt,lub.b,0);
105 //JGFInstrumentor.stopTimer("Section2:LUFact:Kernel", instr.timers);
108 public int JGFvalidate() {
118 for (i = 0; i < n; i++) {
121 norma = matgen(a,lda,n,b);
122 for (i = 0; i < n; i++) {
126 dmxpy(n,b,n,lda,x,a);
129 for (i = 0; i < n; i++) {
130 //resid = (resid > abs(b[i])) ? resid : abs(b[i]);
131 //normx = (normx > abs(x[i])) ? normx : abs(x[i]);
132 if (resid <= abs(b[i])) resid = abs(b[i]);
133 if (normx <= abs(x[i])) normx = abs(x[i]);
135 eps = epslon((double)1.0);
136 residn = resid/( n*norma*normx*eps );
138 /*******************Compare longs ***********/
140 lresidn = (long) residn * 1000000;
141 lref = (long) ref[size] * 1000000;
143 if (lresidn > lref) {
144 //System.printString("Validation failed");
145 System.printString("Computed Norm Res = " + (long) residn * 1000000);
146 System.printString("Reference Norm Res = " + (long) ref[size] * 1000000);
153 double abs (double d) {
154 if (d >= 0) return d;
158 double matgen (double a[][], int lda, int n, double b[])
165 /* Next two for() statements switched. Solver wants
166 matrix in column order. --dmd 3/3/97
168 for (i = 0; i < n; i++) {
169 for (j = 0; j < n; j++) {
170 init = 3125*init % 65536;
171 a[j][i] = (init - 32768.0)/16384.0;
172 if (a[j][i] > norma) {
177 for (i = 0; i < n; i++) {
180 for (j = 0; j < n; j++) {
181 for (i = 0; i < n; i++) {
188 void dgesl( double a[][], int lda, int n, int ipvt[], double b[], int job)
195 // job = 0 , solve a * x = b. first solve l*y = b
197 for (k = 0; k < nm1; k++) {
205 daxpy(n-(kp1),t,a[k],kp1,1,b,kp1,1);
209 for (kb = 0; kb < n; kb++) {
213 daxpy(k,t,a[k],0,1,b,0,1);
216 // job = nonzero, solve trans(a) * x = b. first solve trans(u)*y = b
217 for (k = 0; k < n; k++) {
218 t = ddot(k,a[k],0,1,b,0,1);
219 b[k] = (b[k] - t)/a[k][k];
221 // now solve trans(l)*x = y
223 for (kb = 1; kb < nm1; kb++) {
226 b[k] += ddot(n-(kp1),a[k],kp1,1,b,kp1,1);
239 constant times a vector plus a vector.
240 jack dongarra, linpack, 3/11/78.
242 void daxpy( int n, double da, double dx[], int dx_off, int incx,
243 double dy[], int dy_off, int incy)
247 if ((n > 0) && (da != 0)) {
248 if (incx != 1 || incy != 1) {
249 // code for unequal increments or equal increments not equal to 1
252 if (incx < 0) ix = (-n+1)*incx;
253 if (incy < 0) iy = (-n+1)*incy;
254 for (i = 0;i < n; i++) {
255 dy[iy +dy_off] += da*dx[ix +dx_off];
261 // code for both increments equal to 1
262 for (i=0; i < n; i++)
263 dy[i +dy_off] += da*dx[i +dx_off];
269 forms the dot product of two vectors.
270 jack dongarra, linpack, 3/11/78.
272 double ddot( int n, double dx[], int dx_off, int incx, double dy[],
273 int dy_off, int incy)
279 if (incx != 1 || incy != 1) {
280 // code for unequal increments or equal increments not equal to 1
283 if (incx < 0) ix = (-n+1)*incx;
284 if (incy < 0) iy = (-n+1)*incy;
285 for (i = 0;i < n; i++) {
286 dtemp += dx[ix +dx_off]*dy[iy +dy_off];
291 // code for both increments equal to 1
293 dtemp += dx[i +dx_off]*dy[i +dy_off];
300 scales a vector by a constant.
301 jack dongarra, linpack, 3/11/78.
303 void dscal( int n, double da, double dx[], int dx_off, int incx)
308 // code for increment not equal to 1
310 for (i = 0; i < nincx; i += incx)
313 // code for increment equal to 1
314 for (i = 0; i < n; i++)
321 finds the index of element having max. absolute value.
322 jack dongarra, linpack, 3/11/78.
324 int idamax( int n, double dx[], int dx_off, int incx)
333 } else if (incx != 1) {
334 // code for increment not equal to 1
335 dmax = abs(dx[0 +dx_off]);
337 for (i = 1; i < n; i++) {
338 dtemp = abs(dx[ix + dx_off]);
346 // code for increment equal to 1
348 dmax = abs(dx[0 +dx_off]);
349 for (i = 1; i < n; i++) {
350 dtemp = abs(dx[i + dx_off]);
361 estimate unit roundoff in quantities of size x.
362 this program should function properly on all systems
363 satisfying the following two assumptions,
364 1. the base used in representing dfloating point
365 numbers is not a power of three.
366 2. the quantity a in statement 10 is represented to
367 the accuracy used in dfloating point variables
368 that are stored in memory.
369 the statement number 10 and the go to 10 are intended to
370 force optimizing compilers to generate code satisfying
372 under these assumptions, it should be true that,
373 a is not exactly equal to four-thirds,
374 b has a zero for its last bit or digit,
375 c is not exactly equal to one,
376 eps measures the separation of 1.0 from
377 the next larger dfloating point number.
378 the developers of eispack would appreciate being informed
379 about any systems where these assumptions do not hold.
381 *****************************************************************
382 this routine is one of the auxiliary routines used by eispack iii
383 to avoid machine dependencies.
384 *****************************************************************
386 this version dated 4/6/83.
388 double epslon (double x)
402 void dmxpy ( int n1, double y[], int n2, int ldm, double x[], double m[][])
405 // cleanup odd vector
406 for (j = 0; j < n2; j++) {
407 for (i = 0; i < n1; i++) {
408 y[i] += x[j]*m[j][i];
413 public static void JGFvalidate(JGFLUFactBench lub) {
424 for (i = 0; i < lub.n; i++) {
427 lub.norma = lub.matgen(lub.a,lub.lda,lub.n,lub.b);
428 for (i = 0; i < lub.n; i++) {
429 lub.b[i] = -(lub.b[i]);
432 lub.dmxpy(lub.n,lub.b,lub.n,lub.lda,lub.x,lub.a);
435 for (i = 0; i < lub.n; i++) {
436 //resid = (resid > abs(b[i])) ? resid : abs(b[i]);
437 //normx = (normx > abs(x[i])) ? normx : abs(x[i]);
438 if (lub.resid <= abs(lub.b[i])) lub.resid = lub.abs(lub.b[i]);
439 if (lub.normx <= abs(lub.x[i])) lub.normx = lub.abs(lub.x[i]);
441 eps = lub.epslon((double)1.0);
442 residn = lub.resid/( lub.n*lub.norma*lub.normx*eps );
445 if (residn > ref[size]) {
446 System.printString("Validation failed");
447 System.printString("Computed Norm Res = " + (long) residn);
448 System.printString("Reference Norm Res = " + (long) ref[size]);