Robust/src/Benchmarks/Prefetch/LUFact/dsm/Linpack.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  *            See below for the previous history of this code              *
  16  *                                                                         *
  17  *      This version copyright (c) The University of Edinburgh, 2001.      *
  18  *                         All rights reserved.                            *
  19  *                                                                         *
  20  **************************************************************************/
  21
  22 /*
  23
  24    Modified 3/3/97 by David M. Doolin (dmd) doolin@cs.utk.edu
  25    Fixed error in matgen() method. Added some comments.
  26
  27    Modified 1/22/97 by Paul McMahan mcmahan@cs.utk.edu
  28    Added more MacOS options to form.
  29
  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.
  41
  42
  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"
  46
  47
  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.
  51    Jack Dongarra)
  52
  53 */
  54
  55 public class Linpack {
  56   double a[][];
  57   double b[];
  58   double x[];
  59   double ops,total,norma,normx;
  60   double resid,time;
  61   double kf;
  62   int n,i,ntimes,info,lda,ldaa,kflops;
  63   int ipvt[];
  64
  65   double abs (double d) {
  66     if (d >= 0) return d;
  67     else return -d;
  68     //return (d >= 0) ? d : -d;
  69   }
  70
  71   double matgen (double a[][], int lda, int n, double b[])
  72   {
  73     double norma;
  74     int init, i, j;
  75
  76     init = 1325;
  77     norma = 0.0;
  78     /*  Next two for() statements switched.  Solver wants
  79         matrix in column order. --dmd 3/3/97
  80         */
  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) {
  87           norma = a[j][i];
  88         }
  89       }
  90     }
  91     for (i = 0; i < n; i++) {
  92       b[i] = 0.0;
  93     }
  94     for (j = 0; j < n; j++) {
  95       for (i = 0; i < n; i++) {
  96         b[i] += a[j][i];
  97       }
  98     }
  99     return norma;
 100   }
 101
 102   void dgesl( double a[][], int lda, int n, int ipvt[], double b[], int job)
 103   {
 104     double t;
 105     int k,kb,l,nm1,kp1;
 106
 107     nm1 = n - 1;
 108     if (job == 0) {
 109       // job = 0 , solve  a * x = b.  first solve  l*y = b
 110       if (nm1 >= 1) {
 111         for (k = 0; k < nm1; k++) {
 112           l = ipvt[k];
 113           t = b[l];
 114           if (l != k){
 115             b[l] = b[k];
 116             b[k] = t;
 117           }
 118           kp1 = k + 1;
 119           daxpy(n-(kp1),t,a[k],kp1,1,b,kp1,1);
 120         }
 121       }
 122       // now solve  u*x = y
 123       for (kb = 0; kb < n; kb++) {
 124         k = n - (kb + 1);
 125         b[k] /= a[k][k];
 126         t = -b[k];
 127         daxpy(k,t,a[k],0,1,b,0,1);
 128       }
 129     } else {
 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];
 134       }
 135       // now solve trans(l)*x = y
 136       if (nm1 >= 1) {
 137         for (kb = 1; kb < nm1; kb++) {
 138           k = n - (kb+1);
 139           kp1 = k + 1;
 140           b[k] += ddot(n-(kp1),a[k],kp1,1,b,kp1,1);
 141           l = ipvt[k];
 142           if (l != k) {
 143             t = b[l];
 144             b[l] = b[k];
 145             b[k] = t;
 146           }
 147         }
 148       }
 149     }
 150   }
 151
 152   /*
 153      constant times a vector plus a vector.
 154      jack dongarra, linpack, 3/11/78.
 155      */
 156   void daxpy( int n, double da, double dx[], int dx_off, int incx,
 157       double dy[], int dy_off, int incy)
 158   {
 159     int i,ix,iy;
 160
 161     if ((n > 0) && (da != 0)) {
 162       if (incx != 1 || incy != 1) {
 163         // code for unequal increments or equal increments not equal to 1
 164         ix = 0;
 165         iy = 0;
 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];
 170           ix += incx;
 171           iy += incy;
 172         }
 173         return;
 174       } else {
 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];
 178       }
 179     }
 180   }
 181
 182   /*
 183      forms the dot product of two vectors.
 184      jack dongarra, linpack, 3/11/78.
 185      */
 186   double ddot( int n, double dx[], int dx_off, int incx, double dy[],
 187       int dy_off, int incy)
 188   {
 189     double dtemp;
 190     int i,ix,iy;
 191     dtemp = 0;
 192     if (n > 0) {
 193       if (incx != 1 || incy != 1) {
 194         // code for unequal increments or equal increments not equal to 1
 195         ix = 0;
 196         iy = 0;
 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];
 201           ix += incx;
 202           iy += incy;
 203         }
 204       } else {
 205         // code for both increments equal to 1
 206         for (i=0;i < n; i++)
 207           dtemp += dx[i +dx_off]*dy[i +dy_off];
 208       }
 209     }
 210     return(dtemp);
 211   }
 212
 213   /*
 214      scales a vector by a constant.
 215      jack dongarra, linpack, 3/11/78.
 216      */
 217   void dscal( int n, double da, double dx[], int dx_off, int incx)
 218   {
 219     int i,nincx;
 220     if (n > 0) {
 221       if (incx != 1) {
 222         // code for increment not equal to 1
 223         nincx = n*incx;
 224         for (i = 0; i < nincx; i += incx)
 225           dx[i +dx_off] *= da;
 226       } else {
 227         // code for increment equal to 1
 228         for (i = 0; i < n; i++)
 229           dx[i +dx_off] *= da;
 230       }
 231     }
 232   }
 233
 234   /*
 235      finds the index of element having max. absolute value.
 236      jack dongarra, linpack, 3/11/78.
 237      */
 238   int idamax( int n, double dx[], int dx_off, int incx)
 239   {
 240     double dmax, dtemp;
 241     int i, ix, itemp=0;
 242
 243     if (n < 1) {
 244       itemp = -1;
 245     } else if (n ==1) {
 246       itemp = 0;
 247     } else if (incx != 1) {
 248       // code for increment not equal to 1
 249       dmax = abs(dx[0 +dx_off]);
 250       ix = 1 + incx;
 251       for (i = 1; i < n; i++) {
 252         dtemp = abs(dx[ix + dx_off]);
 253         if (dtemp > dmax)  {
 254           itemp = i;
 255           dmax = dtemp;
 256         }
 257         ix += incx;
 258       }
 259     } else {
 260       // code for increment equal to 1
 261       itemp = 0;
 262       dmax = abs(dx[0 +dx_off]);
 263       for (i = 1; i < n; i++) {
 264         dtemp = abs(dx[i + dx_off]);
 265         if (dtemp > dmax) {
 266           itemp = i;
 267           dmax = dtemp;
 268         }
 269       }
 270     }
 271     return (itemp);
 272   }
 273
 274   /*
 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
 285      assumption 2.
 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.
 294
 295    *****************************************************************
 296    this routine is one of the auxiliary routines used by eispack iii
 297    to avoid machine dependencies.
 298    *****************************************************************
 299
 300    this version dated 4/6/83.
 301    */
 302   double epslon (double x)
 303   {
 304     double a,b,c,eps;
 305
 306     a = 4.0e0/3.0e0;
 307     eps = 0;
 308     while (eps == 0) {
 309       b = a - 1.0;
 310       c = b + b + b;
 311       eps = abs(c-1.0);
 312     }
 313     return(eps*abs(x));
 314   }
 315
 316   void dmxpy ( int n1, double y[], int n2, int ldm, double x[], double m[][])
 317   {
 318     int j,i;
 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];
 323       }
 324     }
 325   }
 326 }
 327