Robust/src/Benchmarks/Prefetch/LUFact/dsm/JGFLUFactBench.java

   1 /**************************************************************************
   2  *                                                                         *
   3  *         Java Grande Forum Benchmark Suite - Thread Version 1.0          *
   4  *                                                                         *
   5  *                            produced by                                  *
   6  *                                                                         *
   7  *                  Java Grande Benchmarking Project                       *
   8  *                                                                         *
   9  *                                at                                       *
  10  *                                                                         *
  11  *                Edinburgh Parallel Computing Centre                      *
  12  *                                                                         *
  13  *                email: epcc-javagrande@epcc.ed.ac.uk                     *
  14  *                                                                         *
  15  *                                                                         *
  16  *      This version copyright (c) The University of Edinburgh, 2001.      *
  17  *                         All rights reserved.                            *
  18  *                                                                         *
  19  **************************************************************************/
  20 public class JGFLUFactBench {
  21   public int nthreads;
  22   private int size;
  23   private int[] datasizes;
  24   public int cachelinesize;
  25   double a[][];
  26   double b[];
  27   double x[];
  28   double ops,total,norma,normx;
  29   double resid,time;
  30   double kf;
  31   int n,i,ntimes,info,lda,ldaa,kflops;
  32   int ipvt[];
  33
  34   public JGFLUFactBench(int nthreads) {
  35     this.nthreads=nthreads;
  36     datasizes = global new int[3];
  37     datasizes[0] = 500;
  38     datasizes[1] = 1000;
  39     datasizes[2] = 2000;
  40     cachelinesize = 128;
  41   }
  42
  43   public void JGFsetsize(int size) {
  44     this.size = size;
  45   }
  46
  47   public void JGFinitialise() {
  48     n = datasizes[size];
  49     ldaa = n;
  50     lda = ldaa + 1;
  51
  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];
  56
  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);
  60   }
  61
  62   public static void JGFkernel(JGFLUFactBench lub) {
  63     int numthreads;
  64     atomic {
  65       numthreads = lub.nthreads;
  66     }
  67
  68     /* spawn threads */
  69     LinpackRunner[] thobjects;
  70     TournamentBarrier br;
  71     atomic {
  72       thobjects = global new LinpackRunner[numthreads];
  73       br = global new TournamentBarrier(numthreads);
  74     }
  75
  76     //JGFInstrumentor.startTimer("Section2:LUFact:Kernel", instr.timers);
  77
  78     LinpackRunner tmp;
  79     int mid = (128<<24)|(195<<16)|(175<<8)|73;
  80     for(int i=1;i<numthreads;i++) {
  81       atomic {
  82         thobjects[i] = global new LinpackRunner(i,lub.a,lub.lda,lub.n,lub.ipvt,br,lub.nthreads);
  83         tmp = thobjects[i];
  84       }
  85       tmp.start(mid);
  86     }
  87
  88     atomic {
  89       thobjects[0] = global new LinpackRunner(0,lub.a,lub.lda,lub.n,lub.ipvt,br,lub.nthreads);
  90       tmp = thobjects[0];
  91     }
  92     tmp.start(mid);
  93     tmp.join();
  94
  95     for(int i=1;i<numthreads;i++) {
  96       atomic {
  97         tmp = thobjects[i];
  98       }
  99       tmp.join();
 100     }
 101
 102     atomic {
 103       lub.dgesl(lub.a,lub.lda,lub.n,lub.ipvt,lub.b,0);
 104     }
 105     //JGFInstrumentor.stopTimer("Section2:LUFact:Kernel", instr.timers);
 106   }
 107
 108   public int JGFvalidate() {
 109     int i;
 110     double eps,residn;
 111     double[] ref;
 112
 113     ref = new double[3];
 114     ref[0] = 6.0;
 115     ref[1] = 12.0;
 116     ref[2] = 20.0;
 117
 118     for (i = 0; i < n; i++) {
 119       x[i] = b[i];
 120     }
 121     norma = matgen(a,lda,n,b);
 122     for (i = 0; i < n; i++) {
 123       b[i] = -(b[i]);
 124     }
 125
 126     dmxpy(n,b,n,lda,x,a);
 127     resid = 0.0;
 128     normx = 0.0;
 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]);
 134     }
 135     eps =  epslon((double)1.0);
 136     residn = resid/( n*norma*normx*eps );
 137
 138     /*******************Compare longs ***********/
 139     long lresidn, lref;
 140     lresidn = (long) residn * 1000000;
 141     lref = (long) ref[size] * 1000000;
 142
 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);
 147       return 1;
 148     } else {
 149       return 0;
 150     }
 151   }
 152
 153   double abs (double d) {
 154     if (d >= 0) return d;
 155     else return -d;
 156   }
 157
 158   double matgen (double a[][], int lda, int n, double b[])
 159   {
 160     double norma;
 161     int init, i, j;
 162
 163     init = 1325;
 164     norma = 0.0;
 165     /*  Next two for() statements switched.  Solver wants
 166         matrix in column order. --dmd 3/3/97
 167         */
 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) {
 173           norma = a[j][i];
 174         }
 175       }
 176     }
 177     for (i = 0; i < n; i++) {
 178       b[i] = 0.0;
 179     }
 180     for (j = 0; j < n; j++) {
 181       for (i = 0; i < n; i++) {
 182         b[i] += a[j][i];
 183       }
 184     }
 185     return norma;
 186   }
 187
 188   void dgesl( double a[][], int lda, int n, int ipvt[], double b[], int job)
 189   {
 190     double t;
 191     int k,kb,l,nm1,kp1;
 192
 193     nm1 = n - 1;
 194     if (job == 0) {
 195       // job = 0 , solve  a * x = b.  first solve  l*y = b
 196       if (nm1 >= 1) {
 197         for (k = 0; k < nm1; k++) {
 198           l = ipvt[k];
 199           t = b[l];
 200           if (l != k){
 201             b[l] = b[k];
 202             b[k] = t;
 203           }
 204           kp1 = k + 1;
 205           daxpy(n-(kp1),t,a[k],kp1,1,b,kp1,1);
 206         }
 207       }
 208       // now solve  u*x = y
 209       for (kb = 0; kb < n; kb++) {
 210         k = n - (kb + 1);
 211         b[k] /= a[k][k];
 212         t = -b[k];
 213         daxpy(k,t,a[k],0,1,b,0,1);
 214       }
 215     } else {
 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];
 220       }
 221       // now solve trans(l)*x = y
 222       if (nm1 >= 1) {
 223         for (kb = 1; kb < nm1; kb++) {
 224           k = n - (kb+1);
 225           kp1 = k + 1;
 226           b[k] += ddot(n-(kp1),a[k],kp1,1,b,kp1,1);
 227           l = ipvt[k];
 228           if (l != k) {
 229             t = b[l];
 230             b[l] = b[k];
 231             b[k] = t;
 232           }
 233         }
 234       }
 235     }
 236   }
 237
 238   /*
 239      constant times a vector plus a vector.
 240      jack dongarra, linpack, 3/11/78.
 241      */
 242   void daxpy( int n, double da, double dx[], int dx_off, int incx,
 243       double dy[], int dy_off, int incy)
 244   {
 245     int i,ix,iy;
 246
 247     if ((n > 0) && (da != 0)) {
 248       if (incx != 1 || incy != 1) {
 249         // code for unequal increments or equal increments not equal to 1
 250         ix = 0;
 251         iy = 0;
 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];
 256           ix += incx;
 257           iy += incy;
 258         }
 259         return;
 260       } else {
 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];
 264       }
 265     }
 266   }
 267
 268   /*
 269      forms the dot product of two vectors.
 270      jack dongarra, linpack, 3/11/78.
 271      */
 272   double ddot( int n, double dx[], int dx_off, int incx, double dy[],
 273       int dy_off, int incy)
 274   {
 275     double dtemp;
 276     int i,ix,iy;
 277     dtemp = 0;
 278     if (n > 0) {
 279       if (incx != 1 || incy != 1) {
 280         // code for unequal increments or equal increments not equal to 1
 281         ix = 0;
 282         iy = 0;
 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];
 287           ix += incx;
 288           iy += incy;
 289         }
 290       } else {
 291         // code for both increments equal to 1
 292         for (i=0;i < n; i++)
 293           dtemp += dx[i +dx_off]*dy[i +dy_off];
 294       }
 295     }
 296     return(dtemp);
 297   }
 298
 299   /*
 300      scales a vector by a constant.
 301      jack dongarra, linpack, 3/11/78.
 302      */
 303   void dscal( int n, double da, double dx[], int dx_off, int incx)
 304   {
 305     int i,nincx;
 306     if (n > 0) {
 307       if (incx != 1) {
 308         // code for increment not equal to 1
 309         nincx = n*incx;
 310         for (i = 0; i < nincx; i += incx)
 311           dx[i +dx_off] *= da;
 312       } else {
 313         // code for increment equal to 1
 314         for (i = 0; i < n; i++)
 315           dx[i +dx_off] *= da;
 316       }
 317     }
 318   }
 319
 320   /*
 321      finds the index of element having max. absolute value.
 322      jack dongarra, linpack, 3/11/78.
 323      */
 324   int idamax( int n, double dx[], int dx_off, int incx)
 325   {
 326     double dmax, dtemp;
 327     int i, ix, itemp=0;
 328
 329     if (n < 1) {
 330       itemp = -1;
 331     } else if (n ==1) {
 332       itemp = 0;
 333     } else if (incx != 1) {
 334       // code for increment not equal to 1
 335       dmax = abs(dx[0 +dx_off]);
 336       ix = 1 + incx;
 337       for (i = 1; i < n; i++) {
 338         dtemp = abs(dx[ix + dx_off]);
 339         if (dtemp > dmax)  {
 340           itemp = i;
 341           dmax = dtemp;
 342         }
 343         ix += incx;
 344       }
 345     } else {
 346       // code for increment equal to 1
 347       itemp = 0;
 348       dmax = abs(dx[0 +dx_off]);
 349       for (i = 1; i < n; i++) {
 350         dtemp = abs(dx[i + dx_off]);
 351         if (dtemp > dmax) {
 352           itemp = i;
 353           dmax = dtemp;
 354         }
 355       }
 356     }
 357     return (itemp);
 358   }
 359
 360   /*
 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
 371      assumption 2.
 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.
 380
 381    *****************************************************************
 382    this routine is one of the auxiliary routines used by eispack iii
 383    to avoid machine dependencies.
 384    *****************************************************************
 385
 386    this version dated 4/6/83.
 387    */
 388   double epslon (double x)
 389   {
 390     double a,b,c,eps;
 391
 392     a = 4.0e0/3.0e0;
 393     eps = 0;
 394     while (eps == 0) {
 395       b = a - 1.0;
 396       c = b + b + b;
 397       eps = abs(c-1.0);
 398     }
 399     return(eps*abs(x));
 400   }
 401
 402   void dmxpy ( int n1, double y[], int n2, int ldm, double x[], double m[][])
 403   {
 404     int j,i;
 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];
 409       }
 410     }
 411   }
 412   /*
 413      public static void JGFvalidate(JGFLUFactBench lub) {
 414      int i;
 415      double eps,residn;
 416      double[] ref;
 417
 418      ref = new double[3];
 419      ref[0] = 6.0;
 420      ref[1] = 12.0;
 421      ref[2] = 20.0;
 422
 423      atomic {
 424      for (i = 0; i < lub.n; i++) {
 425      lub.x[i] = lub.b[i];
 426      }
 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]);
 430      }
 431
 432      lub.dmxpy(lub.n,lub.b,lub.n,lub.lda,lub.x,lub.a);
 433      lub.resid = 0.0;
 434      lub.normx = 0.0;
 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]);
 440 }
 441 eps =  lub.epslon((double)1.0);
 442 residn = lub.resid/( lub.n*lub.norma*lub.normx*eps );
 443 }
 444
 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]);
 449 }
 450 }
 451 */
 452
 453 }