Robust/src/Benchmarks/Prefetch/LUFact/dsm/LinpackRunner.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
  21 class LinpackRunner extends Thread {
  22   int id,lda,n,info,ipvt[];
  23   double a[][];
  24   Barrier br;
  25   int nthreads;
  26
  27   public LinpackRunner(int id, double a[][], int lda, int n, int ipvt[],Barrier br, int nthreads) {
  28     this.id = id;
  29     this.a=a;
  30     this.lda=lda;
  31     this.n=n;
  32     this.ipvt=ipvt;
  33     this.br=br;
  34     this.nthreads = nthreads;
  35   }
  36
  37   double abs (double d) {
  38     if (d >= 0) return d;
  39     else return -d;
  40   }
  41
  42   public void run() {
  43       double[] col_k, col_j;
  44       double t;
  45       int j,k,kp1,l,nm1;
  46       int info;
  47       int slice,ilow,iupper;
  48       // gaussian elimination with partial pivoting
  49       info = 0;
  50       int nlocal;
  51       Barrier tmpbr;
  52       int lid;
  53       atomic {
  54           nlocal=n;
  55           tmpbr=br;
  56           lid=id;
  57       }
  58
  59       nm1 = nlocal - 1;
  60       if (nm1 >=  0) {
  61           for (k = 0; k < nm1; k++) {
  62               atomic {
  63                   col_k = a[k];
  64                   kp1 = k + 1;
  65                   // find l = pivot index
  66                   l = idamax(nlocal-k,col_k,k,1) + k;
  67                   if(lid==0) {
  68                       ipvt[k] = l;
  69                   }
  70               }
  71               // synchronise threads
  72               Barrier.enterBarrier(tmpbr);
  73               System.clearPrefetchCache();
  74
  75               // zero pivot implies this column already triangularized
  76               boolean b;
  77               atomic {
  78                   b=col_k[l]!=0;
  79               }
  80               if (b) {
  81                   Barrier.enterBarrier(tmpbr);
  82                   System.clearPrefetchCache();
  83                   // interchange if necessary
  84                   atomic {
  85                       if(lid == 0 ) {
  86                           if (l != k) {
  87                               t = col_k[l];
  88                               col_k[l] = col_k[k];
  89                               col_k[k] = t;
  90                           }
  91                       }
  92                   }
  93                   // synchronise threads
  94                   Barrier.enterBarrier(tmpbr);
  95                   System.clearPrefetchCache();
  96                   // compute multipliers
  97                   atomic {
  98                       t = -1.0/col_k[k];
  99                       if(lid == 0) {
 100                           dscal(nlocal-(kp1),t,col_k,kp1,1);
 101                       }
 102                   }
 103                   // synchronise threads
 104                   Barrier.enterBarrier(tmpbr);
 105                   System.clearPrefetchCache();
 106                   // row elimination with column indexing
 107                   atomic {
 108                       slice = ((nlocal-kp1) + nthreads-1)/nthreads;
 109                       ilow = (lid*slice)+kp1;
 110                       iupper = ((lid+1)*slice)+kp1;
 111                       if (iupper > nlocal ) iupper=nlocal;
 112                       if (ilow > nlocal ) ilow=nlocal;
 113                       for (j = ilow; j < iupper; j++) {
 114                           col_j = a[j];
 115                           t = col_j[l];
 116                           if (l != k) {
 117                               col_j[l] = col_j[k];
 118                               col_j[k] = t;
 119                           }
 120                           daxpy(nlocal-(kp1),t,col_k,kp1,1,
 121                                 col_j,kp1,1);
 122                       }
 123                   }
 124                   // synchronise threads
 125                   Barrier.enterBarrier(tmpbr);
 126                   System.clearPrefetchCache();
 127               } else {
 128                   info = k;
 129               }
 130               Barrier.enterBarrier(tmpbr);
 131               System.clearPrefetchCache();
 132           }
 133       }
 134
 135       atomic {
 136           if(lid==0) {
 137               ipvt[nlocal-1] = nlocal-1;
 138           }
 139           if (a[(nlocal-1)][(nlocal-1)] == 0) info = nlocal-1;
 140       }
 141   }
 142
 143   /*
 144      finds the index of element having max. absolute value.
 145      jack dongarra, linpack, 3/11/78.
 146      */
 147   int idamax( int n, double dx[], int dx_off, int incx)
 148   {
 149     double dmax, dtemp;
 150     int i, ix, itemp=0;
 151
 152     if (n < 1) {
 153       itemp = -1;
 154     } else if (n ==1) {
 155       itemp = 0;
 156     } else if (incx != 1) {
 157       // code for increment not equal to 1
 158       dmax = abs(dx[0 +dx_off]);
 159       ix = 1 + incx;
 160       for (i = 1; i < n; i++) {
 161         dtemp = abs(dx[ix + dx_off]);
 162         if (dtemp > dmax)  {
 163           itemp = i;
 164           dmax = dtemp;
 165         }
 166         ix += incx;
 167       }
 168     } else {
 169       // code for increment equal to 1
 170       itemp = 0;
 171       dmax = abs(dx[0 +dx_off]);
 172       for (i = 1; i < n; i++) {
 173         dtemp = abs(dx[i + dx_off]);
 174         if (dtemp > dmax) {
 175           itemp = i;
 176           dmax = dtemp;
 177         }
 178       }
 179     }
 180     return (itemp);
 181   }
 182
 183   /*
 184      scales a vector by a constant.
 185      jack dongarra, linpack, 3/11/78.
 186      */
 187   void dscal( int n, double da, double dx[], int dx_off, int incx)
 188   {
 189     int i,nincx;
 190     if (n > 0) {
 191       if (incx != 1) {
 192         // code for increment not equal to 1
 193         nincx = n*incx;
 194         for (i = 0; i < nincx; i += incx)
 195           dx[i +dx_off] *= da;
 196       } else {
 197         // code for increment equal to 1
 198         for (i = 0; i < n; i++)
 199           dx[i +dx_off] *= da;
 200       }
 201     }
 202   }
 203
 204   /*
 205      constant times a vector plus a vector.
 206      jack dongarra, linpack, 3/11/78.
 207      */
 208   void daxpy( int n, double da, double dx[], int dx_off, int incx,
 209       double dy[], int dy_off, int incy)
 210   {
 211     int i,ix,iy;
 212     if ((n > 0) && (da != 0)) {
 213       if (incx != 1 || incy != 1) {
 214         // code for unequal increments or equal increments not equal to 1
 215         ix = 0;
 216         iy = 0;
 217         if (incx < 0) ix = (-n+1)*incx;
 218         if (incy < 0) iy = (-n+1)*incy;
 219         for (i = 0;i < n; i++) {
 220           dy[iy +dy_off] += da*dx[ix +dx_off];
 221           ix += incx;
 222           iy += incy;
 223         }
 224         return;
 225       } else {
 226         // code for both increments equal to 1
 227         for (i=0; i < n; i++)
 228           dy[i +dy_off] += da*dx[i +dx_off];
 229       }
 230     }
 231   }
 232 }
 233