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