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     Barrier br;
  71     atomic {
  72       thobjects = global new LinpackRunner[numthreads];
  73       br = global new Barrier(numthreads);
  74     }
  75
  76     //JGFInstrumentor.startTimer("Section2:LUFact:Kernel", instr.timers);
  77     LinpackRunner tmp;
  78     int[] mid = new int[4];
  79     mid[0] = (128<<24)|(195<<16)|(175<<8)|73;
  80     mid[1] = (128<<24)|(195<<16)|(175<<8)|69;
  81     mid[2] = (128<<24)|(195<<16)|(175<<8)|78;
  82     mid[3] = (128<<24)|(195<<16)|(175<<8)|79;
  83     for(int i=1;i<numthreads;i++) {
  84       atomic {
  85         thobjects[i] = global new LinpackRunner(i,lub.a,lub.lda,lub.n,lub.ipvt,br,lub.nthreads);
  86         tmp = thobjects[i];
  87       }
  88       tmp.start(mid[i]);
  89     }
  90
  91     atomic {
  92       thobjects[0] = global new LinpackRunner(0,lub.a,lub.lda,lub.n,lub.ipvt,br,lub.nthreads);
  93       tmp = thobjects[0];
  94     }
  95     tmp.start(mid[0]);
  96     tmp.join();
  97
  98     for(int i=1;i<numthreads;i++) {
  99       atomic {
 100         tmp = thobjects[i];
 101       }
 102       tmp.join();
 103     }
 104
 105     atomic {
 106       lub.dgesl(lub.a,lub.lda,lub.n,lub.ipvt,lub.b,0);
 107     }
 108
 109     //JGFInstrumentor.stopTimer("Section2:LUFact:Kernel", instr.timers);
 110   }
 111
 112   public int JGFvalidate() {
 113     int i;
 114     double eps,residn;
 115     double[] ref;
 116
 117     ref = new double[3];
 118     ref[0] = 6.0;
 119     ref[1] = 12.0;
 120     ref[2] = 20.0;
 121
 122     for (i = 0; i < n; i++) {
 123       x[i] = b[i];
 124     }
 125     norma = matgen(a,lda,n,b);
 126     for (i = 0; i < n; i++) {
 127       b[i] = -(b[i]);
 128     }
 129
 130     dmxpy(n,b,n,lda,x,a);
 131     resid = 0.0;
 132     normx = 0.0;
 133     for (i = 0; i < n; i++) {
 134       //resid = (resid > abs(b[i])) ? resid : abs(b[i]);
 135       //normx = (normx > abs(x[i])) ? normx : abs(x[i]);
 136       if (resid <= abs(b[i])) resid = abs(b[i]);
 137       if (normx <= abs(x[i])) normx = abs(x[i]);
 138     }
 139     eps =  epslon((double)1.0);
 140     residn = resid/( n*norma*normx*eps );
 141
 142     /*******************Compare longs ***********/
 143     long lresidn, lref;
 144     lresidn = (long) residn * 1000000;
 145     lref = (long) ref[size] * 1000000;
 146
 147     if (lresidn > lref) {
 148       //System.printString("Validation failed");
 149       //System.printString("Computed Norm Res = " + (long) residn * 1000000);
 150       //System.printString("Reference Norm Res = " + (long) ref[size] * 1000000);
 151       return 1;
 152     } else {
 153       return 0;
 154     }
 155   }
 156
 157   double abs (double d) {
 158     if (d >= 0) return d;
 159     else return -d;
 160   }
 161
 162   double matgen (double a[][], int lda, int n, double b[])
 163   {
 164     double norma;
 165     int init, i, j;
 166
 167     init = 1325;
 168     norma = 0.0;
 169     /*  Next two for() statements switched.  Solver wants
 170         matrix in column order. --dmd 3/3/97
 171         */
 172     for (i = 0; i < n; i++) {
 173       for (j = 0; j < n; j++) {
 174         init = 3125*init % 65536;
 175         a[j][i] = (init - 32768.0)/16384.0;
 176         if (a[j][i] > norma) {
 177           norma = a[j][i];
 178         }
 179       }
 180     }
 181     for (i = 0; i < n; i++) {
 182       b[i] = 0.0;
 183     }
 184     for (j = 0; j < n; j++) {
 185       for (i = 0; i < n; i++) {
 186         b[i] += a[j][i];
 187       }
 188     }
 189     return norma;
 190   }
 191
 192   void dgesl( double a[][], int lda, int n, int ipvt[], double b[], int job)
 193   {
 194     double t;
 195     int k,kb,l,nm1,kp1;
 196
 197     nm1 = n - 1;
 198     if (job == 0) {
 199       // job = 0 , solve  a * x = b.  first solve  l*y = b
 200       if (nm1 >= 1) {
 201         for (k = 0; k < nm1; k++) {
 202           l = ipvt[k];
 203           t = b[l];
 204           if (l != k){
 205             b[l] = b[k];
 206             b[k] = t;
 207           }
 208           kp1 = k + 1;
 209           daxpy(n-(kp1),t,a[k],kp1,1,b,kp1,1);
 210         }
 211       }
 212       // now solve  u*x = y
 213       for (kb = 0; kb < n; kb++) {
 214         k = n - (kb + 1);
 215         b[k] /= a[k][k];
 216         t = -b[k];
 217         daxpy(k,t,a[k],0,1,b,0,1);
 218       }
 219     } else {
 220       // job = nonzero, solve  trans(a) * x = b.  first solve  trans(u)*y = b
 221       for (k = 0; k < n; k++) {
 222         t = ddot(k,a[k],0,1,b,0,1);
 223         b[k] = (b[k] - t)/a[k][k];
 224       }
 225       // now solve trans(l)*x = y
 226       if (nm1 >= 1) {
 227         for (kb = 1; kb < nm1; kb++) {
 228           k = n - (kb+1);
 229           kp1 = k + 1;
 230           b[k] += ddot(n-(kp1),a[k],kp1,1,b,kp1,1);
 231           l = ipvt[k];
 232           if (l != k) {
 233             t = b[l];
 234             b[l] = b[k];
 235             b[k] = t;
 236           }
 237         }
 238       }
 239     }
 240   }
 241
 242   /*
 243      constant times a vector plus a vector.
 244      jack dongarra, linpack, 3/11/78.
 245      */
 246   void daxpy( int n, double da, double dx[], int dx_off, int incx,
 247       double dy[], int dy_off, int incy)
 248   {
 249     int i,ix,iy;
 250
 251     if ((n > 0) && (da != 0)) {
 252       if (incx != 1 || incy != 1) {
 253         // code for unequal increments or equal increments not equal to 1
 254         ix = 0;
 255         iy = 0;
 256         if (incx < 0) ix = (-n+1)*incx;
 257         if (incy < 0) iy = (-n+1)*incy;
 258         for (i = 0;i < n; i++) {
 259           dy[iy +dy_off] += da*dx[ix +dx_off];
 260           ix += incx;
 261           iy += incy;
 262         }
 263         return;
 264       } else {
 265         // code for both increments equal to 1
 266         for (i=0; i < n; i++)
 267           dy[i +dy_off] += da*dx[i +dx_off];
 268       }
 269     }
 270   }
 271
 272   /*
 273      forms the dot product of two vectors.
 274      jack dongarra, linpack, 3/11/78.
 275      */
 276   double ddot( int n, double dx[], int dx_off, int incx, double dy[],
 277       int dy_off, int incy)
 278   {
 279     double dtemp;
 280     int i,ix,iy;
 281     dtemp = 0;
 282     if (n > 0) {
 283       if (incx != 1 || incy != 1) {
 284         // code for unequal increments or equal increments not equal to 1
 285         ix = 0;
 286         iy = 0;
 287         if (incx < 0) ix = (-n+1)*incx;
 288         if (incy < 0) iy = (-n+1)*incy;
 289         for (i = 0;i < n; i++) {
 290           dtemp += dx[ix +dx_off]*dy[iy +dy_off];
 291           ix += incx;
 292           iy += incy;
 293         }
 294       } else {
 295         // code for both increments equal to 1
 296         for (i=0;i < n; i++)
 297           dtemp += dx[i +dx_off]*dy[i +dy_off];
 298       }
 299     }
 300     return(dtemp);
 301   }
 302
 303   /*
 304      scales a vector by a constant.
 305      jack dongarra, linpack, 3/11/78.
 306      */
 307   void dscal( int n, double da, double dx[], int dx_off, int incx)
 308   {
 309     int i,nincx;
 310     if (n > 0) {
 311       if (incx != 1) {
 312         // code for increment not equal to 1
 313         nincx = n*incx;
 314         for (i = 0; i < nincx; i += incx)
 315           dx[i +dx_off] *= da;
 316       } else {
 317         // code for increment equal to 1
 318         for (i = 0; i < n; i++)
 319           dx[i +dx_off] *= da;
 320       }
 321     }
 322   }
 323
 324   /*
 325      finds the index of element having max. absolute value.
 326      jack dongarra, linpack, 3/11/78.
 327      */
 328   int idamax( int n, double dx[], int dx_off, int incx)
 329   {
 330     double dmax, dtemp;
 331     int i, ix, itemp=0;
 332
 333     if (n < 1) {
 334       itemp = -1;
 335     } else if (n ==1) {
 336       itemp = 0;
 337     } else if (incx != 1) {
 338       // code for increment not equal to 1
 339       dmax = abs(dx[0 +dx_off]);
 340       ix = 1 + incx;
 341       for (i = 1; i < n; i++) {
 342         dtemp = abs(dx[ix + dx_off]);
 343         if (dtemp > dmax)  {
 344           itemp = i;
 345           dmax = dtemp;
 346         }
 347         ix += incx;
 348       }
 349     } else {
 350       // code for increment equal to 1
 351       itemp = 0;
 352       dmax = abs(dx[0 +dx_off]);
 353       for (i = 1; i < n; i++) {
 354         dtemp = abs(dx[i + dx_off]);
 355         if (dtemp > dmax) {
 356           itemp = i;
 357           dmax = dtemp;
 358         }
 359       }
 360     }
 361     return (itemp);
 362   }
 363
 364   /*
 365      estimate unit roundoff in quantities of size x.
 366      this program should function properly on all systems
 367      satisfying the following two assumptions,
 368      1.  the base used in representing dfloating point
 369      numbers is not a power of three.
 370      2.  the quantity  a  in statement 10 is represented to
 371      the accuracy used in dfloating point variables
 372      that are stored in memory.
 373      the statement number 10 and the go to 10 are intended to
 374      force optimizing compilers to generate code satisfying
 375      assumption 2.
 376      under these assumptions, it should be true that,
 377      a  is not exactly equal to four-thirds,
 378      b  has a zero for its last bit or digit,
 379      c  is not exactly equal to one,
 380      eps  measures the separation of 1.0 from
 381      the next larger dfloating point number.
 382      the developers of eispack would appreciate being informed
 383      about any systems where these assumptions do not hold.
 384
 385    *****************************************************************
 386    this routine is one of the auxiliary routines used by eispack iii
 387    to avoid machine dependencies.
 388    *****************************************************************
 389
 390    this version dated 4/6/83.
 391    */
 392   double epslon (double x)
 393   {
 394     double a,b,c,eps;
 395
 396     a = 4.0e0/3.0e0;
 397     eps = 0;
 398     while (eps == 0) {
 399       b = a - 1.0;
 400       c = b + b + b;
 401       eps = abs(c-1.0);
 402     }
 403     return(eps*abs(x));
 404   }
 405
 406   void dmxpy ( int n1, double y[], int n2, int ldm, double x[], double m[][])
 407   {
 408     int j,i;
 409     // cleanup odd vector
 410     for (j = 0; j < n2; j++) {
 411       for (i = 0; i < n1; i++) {
 412         y[i] += x[j]*m[j][i];
 413       }
 414     }
 415   }
 416 }