Added in everything!

Nikhil Sarda [2012-08-12 13:35:53]
Added in everything!
Filename
testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/iobench/FileRead.java
testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/Constants.java
testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/FFT.java
testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/Jacobi.java
testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/LU.java
testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/MonteCarlo.java
testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/README
testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/Random.java
testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/SOR.java
testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/SparseCompRow.java
testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/Stopwatch.java
testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/commandline.java
testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/kernel.java
diff --git a/testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/iobench/FileRead.java b/testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/iobench/FileRead.java
new file mode 100644
index 0000000..6718a84
--- /dev/null
+++ b/testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/iobench/FileRead.java
@@ -0,0 +1,34 @@
+package edu.columbia.cs.psl.invivo.wallace.iobench;
+
+import java.io.BufferedReader;
+import java.io.FileNotFoundException;
+import java.io.FileReader;
+import java.io.IOException;
+import java.util.ArrayList;
+import java.util.List;
+
+public class FileRead {
+	public static void main(String[] args) {
+		long start = System.currentTimeMillis();
+		FileReader fileReader = null;
+		try {
+			fileReader = new FileReader("/home/nikhil/Downloads/alldownloads.tar.gz");
+		} catch (FileNotFoundException e) {
+			// TODO Auto-generated catch block
+			e.printStackTrace();
+		}
+		BufferedReader bufferedReader = new BufferedReader(fileReader);
+		List<String> lines = new ArrayList<String>();
+		String line = null;
+		try {
+			while ((line = bufferedReader.readLine()) != null) {
+			}
+			bufferedReader.close();
+		} catch (IOException e) {
+			// TODO Auto-generated catch block
+			e.printStackTrace();
+		}
+		long end = System.currentTimeMillis();
+		System.out.println((end-start));
+	}
+}
diff --git a/testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/Constants.java b/testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/Constants.java
new file mode 100644
index 0000000..9ec1951
--- /dev/null
+++ b/testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/Constants.java
@@ -0,0 +1,37 @@
+package edu.columbia.cs.psl.invivo.wallace.numericalbench;
+
+public class Constants
+{
+
+	public static final double RESOLUTION_DEFAULT = 2.0;  /*secs*/
+	public static final int RANDOM_SEED = 101010;
+
+	// default: small (cache-contained) problem sizes
+	//
+	public static final int FFT_SIZE = 1024;  // must be a power of two
+	public static final int SOR_SIZE =100; // NxN grid
+	public static final int SPARSE_SIZE_M = 1000;
+	public static final int SPARSE_SIZE_nz = 5000;
+	public static final int LU_SIZE = 100;
+
+	// large (out-of-cache) problem sizes
+	//
+	public static final int LG_FFT_SIZE = 1048576;  // must be a power of two
+	public static final int LG_SOR_SIZE =1000; // NxN grid
+	public static final int LG_SPARSE_SIZE_M = 100000;
+	public static final int LG_SPARSE_SIZE_nz =1000000;
+	public static final int LG_LU_SIZE = 1000;
+
+	// tiny problem sizes (used to mainly to preload network classes
+	//                     for applet, so that network download times
+	//                     are factored out of benchmark.)
+	//
+	public static final int TINY_FFT_SIZE = 16;  // must be a power of two
+	public static final int TINY_SOR_SIZE =10; // NxN grid
+	public static final int TINY_SPARSE_SIZE_M = 10;
+	public static final int TINY_SPARSE_SIZE_N = 10;
+	public static final int TINY_SPARSE_SIZE_nz = 50;
+	public static final int TINY_LU_SIZE = 10;
+
+}
+
diff --git a/testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/FFT.java b/testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/FFT.java
new file mode 100644
index 0000000..e07bb9f
--- /dev/null
+++ b/testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/FFT.java
@@ -0,0 +1,192 @@
+package edu.columbia.cs.psl.invivo.wallace.numericalbench;
+
+/** Computes FFT's of complex, double precision data where n is an integer power of 2.
+  * This appears to be slower than the Radix2 method,
+  * but the code is smaller and simpler, and it requires no extra storage.
+  * <P>
+  *
+  * @author Bruce R. Miller bruce.miller@nist.gov,
+  * @author Derived from GSL (Gnu Scientific Library),
+  * @author GSL's FFT Code by Brian Gough bjg@vvv.lanl.gov
+  */
+
+  /* See {@link ComplexDoubleFFT ComplexDoubleFFT} for details of data layout.
+   */
+
+public class FFT {
+
+  public static final double num_flops(int N)
+  {
+	 double Nd = (double) N;
+	 double logN = (double) log2(N);
+
+	 return (5.0*Nd-2)*logN + 2*(Nd+1);
+   }
+
+
+  /** Compute Fast Fourier Transform of (complex) data, in place.*/
+  public static void transform (double data[]) {
+    transform_internal(data, -1); }
+
+  /** Compute Inverse Fast Fourier Transform of (complex) data, in place.*/
+  public static void inverse (double data[]) {
+    transform_internal(data, +1);
+    // Normalize
+    int nd=data.length;
+    int n =nd/2;
+    double norm=1/((double) n);
+    for(int i=0; i<nd; i++)
+      data[i] *= norm;
+  }
+
+  /** Accuracy check on FFT of data. Make a copy of data, Compute the FFT, then
+    * the inverse and compare to the original.  Returns the rms difference.*/
+  public static double test(double data[]){
+    int nd = data.length;
+    // Make duplicate for comparison
+    double copy[] = new double[nd];
+    System.arraycopy(data,0,copy,0,nd);
+    // Transform & invert
+    transform(data);
+    inverse(data);
+    // Compute RMS difference.
+    double diff = 0.0;
+    for(int i=0; i<nd; i++) {
+      double d = data[i]-copy[i];
+      diff += d*d; }
+    return Math.sqrt(diff/nd); }
+
+  /** Make a random array of n (complex) elements. */
+  public static double[] makeRandom(int n){
+    int nd = 2*n;
+    double data[] = new double[nd];
+    for(int i=0; i<nd; i++)
+      data[i]= Math.random();
+    return data; }
+
+  /** Simple Test routine. */
+  public static void main(String args[]){
+    if (args.length == 0) {
+      int n = 1024;
+      System.out.println("n="+n+" => RMS Error="+test(makeRandom(n))); }
+    for(int i=0; i<args.length; i++) {
+      int n = Integer.parseInt(args[i]);
+      System.out.println("n="+n+" => RMS Error="+test(makeRandom(n))); }
+  }
+  /* ______________________________________________________________________ */
+
+  protected static int log2 (int n){
+    int log = 0;
+    for(int k=1; k < n; k *= 2, log++);
+    if (n != (1 << log))
+      throw new Error("FFT: Data length is not a power of 2!: "+n);
+    return log; }
+
+  protected static void transform_internal (double data[], int direction) {
+	if (data.length == 0) return;
+	int n = data.length/2;
+    if (n == 1) return;         // Identity operation!
+    int logn = log2(n);
+
+    /* bit reverse the input data for decimation in time algorithm */
+    bitreverse(data) ;
+
+    /* apply fft recursion */
+	/* this loop executed log2(N) times */
+    for (int bit = 0, dual = 1; bit < logn; bit++, dual *= 2) {
+      double w_real = 1.0;
+      double w_imag = 0.0;
+
+      double theta = 2.0 * direction * Math.PI / (2.0 * (double) dual);
+      double s = Math.sin(theta);
+      double t = Math.sin(theta / 2.0);
+      double s2 = 2.0 * t * t;
+
+      /* a = 0 */
+      for (int b = 0; b < n; b += 2 * dual) {
+        int i = 2*b ;
+        int j = 2*(b + dual);
+
+        double wd_real = data[j] ;
+        double wd_imag = data[j+1] ;
+
+        data[j]   = data[i]   - wd_real;
+        data[j+1] = data[i+1] - wd_imag;
+        data[i]  += wd_real;
+        data[i+1]+= wd_imag;
+      }
+
+      /* a = 1 .. (dual-1) */
+      for (int a = 1; a < dual; a++) {
+        /* trignometric recurrence for w-> exp(i theta) w */
+        {
+          double tmp_real = w_real - s * w_imag - s2 * w_real;
+          double tmp_imag = w_imag + s * w_real - s2 * w_imag;
+          w_real = tmp_real;
+          w_imag = tmp_imag;
+        }
+        for (int b = 0; b < n; b += 2 * dual) {
+          int i = 2*(b + a);
+          int j = 2*(b + a + dual);
+
+          double z1_real = data[j];
+          double z1_imag = data[j+1];
+
+          double wd_real = w_real * z1_real - w_imag * z1_imag;
+          double wd_imag = w_real * z1_imag + w_imag * z1_real;
+
+          data[j]   = data[i]   - wd_real;
+          data[j+1] = data[i+1] - wd_imag;
+          data[i]  += wd_real;
+          data[i+1]+= wd_imag;
+        }
+      }
+    }
+  }
+
+
+  protected static void bitreverse(double data[]) {
+    /* This is the Goldrader bit-reversal algorithm */
+    int n=data.length/2;
+	int nm1 = n-1;
+	int i=0;
+	int j=0;
+    for (; i < nm1; i++) {
+
+      //int ii = 2*i;
+      int ii = i << 1;
+
+      //int jj = 2*j;
+      int jj = j << 1;
+
+      //int k = n / 2 ;
+      int k = n >> 1;
+
+      if (i < j) {
+        double tmp_real    = data[ii];
+        double tmp_imag    = data[ii+1];
+        data[ii]   = data[jj];
+        data[ii+1] = data[jj+1];
+        data[jj]   = tmp_real;
+        data[jj+1] = tmp_imag; }
+
+      while (k <= j)
+	  {
+        //j = j - k ;
+		j -= k;
+
+        //k = k / 2 ;
+        k >>= 1 ;
+	  }
+      j += k ;
+    }
+  }
+}
+
+
+
+
+
+
+
+
diff --git a/testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/Jacobi.java b/testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/Jacobi.java
new file mode 100644
index 0000000..5f71b26
--- /dev/null
+++ b/testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/Jacobi.java
@@ -0,0 +1,40 @@
+package edu.columbia.cs.psl.invivo.wallace.numericalbench;
+
+public class Jacobi
+{
+	public static final double num_flops(int M, int N, int num_iterations)
+	{
+		double Md = (double) M;
+		double Nd = (double) N;
+		double num_iterD = (double) num_iterations;
+
+		return (Md-1)*(Nd-1)*num_iterD*6.0;
+	}
+
+	public static final void SOR(double omega, double G[][], int num_iterations)
+	{
+		int M = G.length;
+		int N = G[0].length;
+
+		double omega_over_four = omega * 0.25;
+		double one_minus_omega = 1.0 - omega;
+
+		// update interior points
+		//
+		int Mm1 = M-1;
+		int Nm1 = N-1;
+		for (int p=0; p<num_iterations; p++)
+		{
+			for (int i=1; i<Mm1; i++)
+			{
+				double[] Gi = G[i];
+				double[] Gim1 = G[i-1];
+				double[] Gip1 = G[i+1];
+				for (int j=1; j<Nm1; j++)
+					Gi[j] = omega_over_four * (Gim1[j] + Gip1[j] + Gi[j-1]
+								+ Gi[j+1]) + one_minus_omega * Gi[j];
+			}
+		}
+	}
+}
+
diff --git a/testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/LU.java b/testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/LU.java
new file mode 100644
index 0000000..e41ae7d
--- /dev/null
+++ b/testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/LU.java
@@ -0,0 +1,290 @@
+package edu.columbia.cs.psl.invivo.wallace.numericalbench;
+
+/**
+    LU matrix factorization. (Based on TNT implementation.)
+    Decomposes a matrix A  into a triangular lower triangular
+    factor (L) and an upper triangular factor (U) such that
+    A = L*U.  By convnetion, the main diagonal of L consists
+    of 1's so that L and U can be stored compactly in
+    a NxN matrix.
+
+
+*/
+public class LU
+{
+    /**
+        Returns a <em>copy</em> of the compact LU factorization.
+        (useful mainly for debugging.)
+
+        @return the compact LU factorization.  The U factor
+        is stored in the upper triangular portion, and the L
+        factor is stored in the lower triangular portion.
+        The main diagonal of L consists (by convention) of
+        ones, and is not explicitly stored.
+    */
+
+
+	public static final double num_flops(int N)
+	{
+		// rougly 2/3*N^3
+
+		double Nd = (double) N;
+
+		return (2.0 * Nd *Nd *Nd/ 3.0);
+	}
+
+    protected static double[] new_copy(double x[])
+    {
+        int N = x.length;
+        double T[] = new double[N];
+        for (int i=0; i<N; i++)
+            T[i] = x[i];
+        return T;
+    }
+
+
+    protected static double[][] new_copy(double A[][])
+    {
+        int M = A.length;
+        int N = A[0].length;
+
+        double T[][] = new double[M][N];
+
+        for (int i=0; i<M; i++)
+        {
+            double Ti[] = T[i];
+            double Ai[] = A[i];
+            for (int j=0; j<N; j++)
+                Ti[j] = Ai[j];
+        }
+
+        return T;
+    }
+
+
+
+    public static int[] new_copy(int x[])
+    {
+        int N = x.length;
+        int T[] = new int[N];
+        for (int i=0; i<N; i++)
+            T[i] = x[i];
+        return T;
+    }
+
+    protected static final void insert_copy(double B[][], double A[][])
+    {
+        int M = A.length;
+        int N = A[0].length;
+
+		int remainder = N & 3;		 // N mod 4;
+
+        for (int i=0; i<M; i++)
+        {
+            double Bi[] = B[i];
+            double Ai[] = A[i];
+			for (int j=0; j<remainder; j++)
+                Bi[j] = Ai[j];
+            for (int j=remainder; j<N; j+=4)
+			{
+				Bi[j] = Ai[j];
+				Bi[j+1] = Ai[j+1];
+				Bi[j+2] = Ai[j+2];
+				Bi[j+3] = Ai[j+3];
+			}
+		}
+
+    }
+    public double[][] getLU()
+    {
+        return new_copy(LU_);
+    }
+
+    /**
+        Returns a <em>copy</em> of the pivot vector.
+
+        @return the pivot vector used in obtaining the
+        LU factorzation.  Subsequent solutions must
+        permute the right-hand side by this vector.
+
+    */
+    public int[] getPivot()
+    {
+        return new_copy(pivot_);
+    }
+
+    /**
+        Initalize LU factorization from matrix.
+
+        @param A (in) the matrix to associate with this
+                factorization.
+    */
+    public LU( double A[][] )
+    {
+        int M = A.length;
+        int N = A[0].length;
+
+        //if ( LU_ == null || LU_.length != M || LU_[0].length != N)
+            LU_ = new double[M][N];
+
+        insert_copy(LU_, A);
+
+        //if (pivot_.length != M)
+            pivot_ = new int[M];
+
+        factor(LU_, pivot_);
+    }
+
+    /**
+        Solve a linear system, with pre-computed factorization.
+
+        @param b (in) the right-hand side.
+        @return solution vector.
+    */
+    public double[] solve(double b[])
+    {
+        double x[] = new_copy(b);
+
+        solve(LU_, pivot_, x);
+        return x;
+    }
+
+
+/**
+    LU factorization (in place).
+
+    @param A (in/out) On input, the matrix to be factored.
+        On output, the compact LU factorization.
+
+    @param pivit (out) The pivot vector records the
+        reordering of the rows of A during factorization.
+
+    @return 0, if OK, nozero value, othewise.
+*/
+public static int factor(double A[][],  int pivot[])
+{
+
+
+
+    int N = A.length;
+    int M = A[0].length;
+
+    int minMN = Math.min(M,N);
+
+    for (int j=0; j<minMN; j++)
+    {
+        // find pivot in column j and  test for singularity.
+
+        int jp=j;
+
+        double t = Math.abs(A[j][j]);
+        for (int i=j+1; i<M; i++)
+        {
+            double ab = Math.abs(A[i][j]);
+            if ( ab > t)
+            {
+                jp = i;
+                t = ab;
+            }
+        }
+
+        pivot[j] = jp;
+
+        // jp now has the index of maximum element
+        // of column j, below the diagonal
+
+        if ( A[jp][j] == 0 )
+            return 1;       // factorization failed because of zero pivot
+
+
+        if (jp != j)
+        {
+            // swap rows j and jp
+            double tA[] = A[j];
+            A[j] = A[jp];
+            A[jp] = tA;
+        }
+
+        if (j<M-1)                // compute elements j+1:M of jth column
+        {
+            // note A(j,j), was A(jp,p) previously which was
+            // guarranteed not to be zero (Label #1)
+            //
+            double recp =  1.0 / A[j][j];
+
+            for (int k=j+1; k<M; k++)
+                A[k][j] *= recp;
+        }
+
+
+        if (j < minMN-1)
+        {
+            // rank-1 update to trailing submatrix:   E = E - x*y;
+            //
+            // E is the region A(j+1:M, j+1:N)
+            // x is the column vector A(j+1:M,j)
+            // y is row vector A(j,j+1:N)
+
+
+            for (int ii=j+1; ii<M; ii++)
+            {
+                double Aii[] = A[ii];
+                double Aj[] = A[j];
+                double AiiJ = Aii[j];
+                for (int jj=j+1; jj<N; jj++)
+                  Aii[jj] -= AiiJ * Aj[jj];
+
+            }
+        }
+    }
+
+    return 0;
+}
+
+
+    /**
+        Solve a linear system, using a prefactored matrix
+            in LU form.
+
+
+        @param LU (in) the factored matrix in LU form.
+        @param pivot (in) the pivot vector which lists
+            the reordering used during the factorization
+            stage.
+        @param b    (in/out) On input, the right-hand side.
+                    On output, the solution vector.
+    */
+    public static void solve(double LU[][], int pvt[], double b[])
+    {
+        int M = LU.length;
+        int N = LU[0].length;
+        int ii=0;
+
+        for (int i=0; i<M; i++)
+        {
+            int ip = pvt[i];
+            double sum = b[ip];
+
+            b[ip] = b[i];
+            if (ii==0)
+                for (int j=ii; j<i; j++)
+                    sum -= LU[i][j] * b[j];
+            else
+                if (sum == 0.0)
+                    ii = i;
+            b[i] = sum;
+        }
+
+        for (int i=N-1; i>=0; i--)
+        {
+            double sum = b[i];
+            for (int j=i+1; j<N; j++)
+                sum -= LU[i][j] * b[j];
+            b[i] = sum / LU[i][i];
+        }
+    }
+
+
+    private double LU_[][];
+    private int pivot_[];
+}
diff --git a/testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/MonteCarlo.java b/testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/MonteCarlo.java
new file mode 100644
index 0000000..68790a7
--- /dev/null
+++ b/testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/MonteCarlo.java
@@ -0,0 +1,67 @@
+package edu.columbia.cs.psl.invivo.wallace.numericalbench;
+
+/**
+ Estimate Pi by approximating the area of a circle.
+
+ How: generate N random numbers in the unit square, (0,0) to (1,1)
+ and see how are within a radius of 1 or less, i.e.
+ <pre>
+
+ sqrt(x^2 + y^2) < r
+
+ </pre>
+  since the radius is 1.0, we can square both sides
+  and avoid a sqrt() computation:
+  <pre>
+
+    x^2 + y^2 <= 1.0
+
+  </pre>
+  this area under the curve is (Pi * r^2)/ 4.0,
+  and the area of the unit of square is 1.0,
+  so Pi can be approximated by
+  <pre>
+		        # points with x^2+y^2 < 1
+     Pi =~ 		--------------------------  * 4.0
+		             total # points
+
+  </pre>
+
+*/
+
+public class MonteCarlo
+{
+	final static int SEED = 113;
+
+	public static final double num_flops(int Num_samples)
+	{
+		// 3 flops in x^2+y^2 and 1 flop in random routine
+
+		return ((double) Num_samples)* 4.0;
+
+	}
+
+
+
+	public static final double integrate(int Num_samples)
+	{
+
+		Random R = new Random(SEED);
+
+
+		int under_curve = 0;
+		for (int count=0; count<Num_samples; count++)
+		{
+			double x= R.nextDouble();
+			double y= R.nextDouble();
+
+			if ( x*x + y*y <= 1.0)
+				 under_curve ++;
+
+		}
+
+		return ((double) under_curve / Num_samples) * 4.0;
+	}
+
+
+}
diff --git a/testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/README b/testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/README
new file mode 100644
index 0000000..63b4ba8
--- /dev/null
+++ b/testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/README
@@ -0,0 +1,67 @@
+
+        SciMark 2.0  Java Numerical Benchmark
+
+        Roldan Pozo, Bruce Miller
+
+        NIST
+
+SciMark 2.0 is a composite Java benchmark measuring the  performance of
+numerical kernels occurring in scientific and engineering applications.
+It consists of five kernels which typify computational routines
+commonly found in numeric codes: Fast Fourier Transforms (FFTs),
+Jacobi Successive Over-relaxation (SOR), Sparse matrix-multiply,
+Monte Carlo integration, and dense LU matrix factorization.
+
+(See http://www.math.nist.gov/scimark for further information
+and latest updates.)
+
+
+1) INSTALLATION
+
+Unpack the contents of archive into a subdirectory on your
+CLASSPATH.  Be sure to keep the directory structure of the
+file contents.
+
+2) COMPILING THE BENCHMARKS (optional)
+
+From the directory above this one, issue the command:
+
+    >javac -O commandline.java
+
+This should compile main benchmark driver and dependent files.
+
+3) RUNNING THE BENCHMARKS
+
+From the directory above this one, issue the command:
+
+    >java jnt.scimark2.commandline
+
+or
+    >java jnt.scimark2.commandline -large
+
+to run the large problem size version. (Note that this one
+takes considerably longer to run.)
+
+After a few minutes, the program should respond with
+the benchmark results, e.g.
+
+    >javac jnt.scimark2.commandline
+
+	SciMark 2.0a
+
+    Composite Score: 20.791595999749727
+    FFT (4096): 30.260047144878346
+    Jacobi SOR (100x100):   33.074935359763934
+    Monte Carlo (25000): 11.510791361970528
+    Sparse matmult (nz=25000), 10 iterations: 8.007507030681996
+    LU (100x100): 21.104699101453836
+
+    java.vendor: Sun Microsystems Inc.
+    java.version: 1.2
+    os.arch: x86
+    os.name: Windows NT
+    os.version: 4.0
+
+One can send these results to "pozo@nist.gov".
+
+/* ----------------------  END OF README -----------------------------*/
diff --git a/testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/Random.java b/testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/Random.java
new file mode 100644
index 0000000..26173cf
--- /dev/null
+++ b/testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/Random.java
@@ -0,0 +1,268 @@
+package edu.columbia.cs.psl.invivo.wallace.numericalbench;
+
+/* Random.java based on Java Numerical Toolkit (JNT) Random.UniformSequence
+	class.  We do not use Java's own java.util.Random so that we can compare
+	results with equivalent C and Fortran coces.
+*/
+
+public class Random {
+
+
+/* ------------------------------------------------------------------------------
+                               CLASS VARIABLES
+   ------------------------------------------------------------------------------ */
+
+  int seed = 0;
+
+  private int m[];
+  private int i = 4;
+  private int j = 16;
+
+  private final int mdig = 32;
+  private final int one = 1;
+  private final int m1 = (one << mdig-2) + ((one << mdig-2)-one);
+  private final int m2 = one << mdig/2;
+
+  /* For mdig = 32 : m1 =          2147483647, m2 =      65536
+     For mdig = 64 : m1 = 9223372036854775807, m2 = 4294967296
+  */
+
+  private double dm1 = 1.0 / (double) m1;
+
+  private boolean haveRange = false;
+  private double left  = 0.0;
+  private double right = 1.0;
+  private double width = 1.0;
+
+
+/* ------------------------------------------------------------------------------
+                                CONSTRUCTORS
+   ------------------------------------------------------------------------------ */
+
+/**
+   Initializes a sequence of uniformly distributed quasi random numbers with a
+   seed based on the system clock.
+*/
+  public Random () {
+    initialize( (int) System.currentTimeMillis());
+  }
+
+/**
+   Initializes a sequence of uniformly distributed quasi random numbers on a
+   given half-open interval [left,right) with a seed based on the system
+   clock.
+
+@param <B>left</B> (double)<BR>
+
+       The left endpoint of the half-open interval [left,right).
+
+@param <B>right</B> (double)<BR>
+
+       The right endpoint of the half-open interval [left,right).
+*/
+  public Random ( double left, double right) {
+    initialize( (int) System.currentTimeMillis() );
+    this.left = left;
+    this.right = right;
+    width = right - left;
+    haveRange = true;
+  }
+
+/**
+   Initializes a sequence of uniformly distributed quasi random numbers with a
+   given seed.
+
+@param <B>seed</B> (int)<BR>
+
+       The seed of the random number generator.  Two sequences with the same
+       seed will be identical.
+*/
+  public Random (int seed) {
+    initialize( seed);
+  }
+
+/**
+   Initializes a sequence of uniformly distributed quasi random numbers
+   with a given seed on a given half-open interval [left,right).
+
+@param <B>seed</B> (int)<BR>
+
+       The seed of the random number generator.  Two sequences with the same
+       seed will be identical.
+
+@param <B>left</B> (double)<BR>
+
+       The left endpoint of the half-open interval [left,right).
+
+@param <B>right</B> (double)<BR>
+
+       The right endpoint of the half-open interval [left,right).
+*/
+  public Random (int seed, double left, double right) {
+    initialize( seed);
+    this.left = left;
+    this.right = right;
+    width = right - left;
+    haveRange = true;
+  }
+
+/* ------------------------------------------------------------------------------
+                             PUBLIC METHODS
+   ------------------------------------------------------------------------------ */
+
+/**
+   Returns the next random number in the sequence.
+*/
+  public final synchronized double nextDouble () {
+
+    int k;
+    double nextValue;
+
+    k = m[i] - m[j];
+    if (k < 0) k += m1;
+    m[j] = k;
+
+    if (i == 0)
+		i = 16;
+	else i--;
+
+    if (j == 0)
+		j = 16 ;
+	else j--;
+
+    if (haveRange)
+		return  left +  dm1 * (double) k * width;
+	else
+		return dm1 * (double) k;
+
+  }
+
+/**
+   Returns the next N random numbers in the sequence, as
+   a vector.
+*/
+  public final synchronized void nextDoubles (double x[])
+  {
+
+	int N = x.length;
+	int remainder = N & 3;		// N mod 4
+
+	if (haveRange)
+	{
+		for (int count=0; count<N; count++)
+		{
+     		int k = m[i] - m[j];
+
+     		if (i == 0) i = 16;
+	 			else i--;
+
+     		if (k < 0) k += m1;
+     		m[j] = k;
+
+     		if (j == 0) j = 16;
+	 			else j--;
+
+     		x[count] = left + dm1 * (double) k * width;
+		}
+
+	}
+	else
+	{
+
+		for (int count=0; count<remainder; count++)
+		{
+     		int k = m[i] - m[j];
+
+     		if (i == 0) i = 16;
+	 			else i--;
+
+     		if (k < 0) k += m1;
+     		m[j] = k;
+
+     		if (j == 0) j = 16;
+	 			else j--;
+
+
+     		x[count] = dm1 * (double) k;
+		}
+
+		for (int count=remainder; count<N; count+=4)
+		{
+     		int k = m[i] - m[j];
+     		if (i == 0) i = 16;
+	 			else i--;
+     		if (k < 0) k += m1;
+     		m[j] = k;
+     		if (j == 0) j = 16;
+	 			else j--;
+     		x[count] = dm1 * (double) k;
+
+
+     		k = m[i] - m[j];
+     		if (i == 0) i = 16;
+	 			else i--;
+     		if (k < 0) k += m1;
+     		m[j] = k;
+     		if (j == 0) j = 16;
+	 			else j--;
+     		x[count+1] = dm1 * (double) k;
+
+
+     		k = m[i] - m[j];
+     		if (i == 0) i = 16;
+	 			else i--;
+     		if (k < 0) k += m1;
+     		m[j] = k;
+     		if (j == 0) j = 16;
+	 			else j--;
+     		x[count+2] = dm1 * (double) k;
+
+
+     		k = m[i] - m[j];
+     		if (i == 0) i = 16;
+	 			else i--;
+     		if (k < 0) k += m1;
+     		m[j] = k;
+     		if (j == 0) j = 16;
+	 			else j--;
+     		x[count+3] = dm1 * (double) k;
+		}
+	}
+  }
+
+
+
+
+
+
+/*----------------------------------------------------------------------------
+                           PRIVATE METHODS
+  ------------------------------------------------------------------------ */
+
+  private void initialize (int seed) {
+
+    int jseed, k0, k1, j0, j1, iloop;
+
+    this.seed = seed;
+
+    m = new int[17];
+
+    jseed = Math.min(Math.abs(seed),m1);
+    if (jseed % 2 == 0) --jseed;
+    k0 = 9069 % m2;
+    k1 = 9069 / m2;
+    j0 = jseed % m2;
+    j1 = jseed / m2;
+    for (iloop = 0; iloop < 17; ++iloop)
+	{
+		jseed = j0 * k0;
+		j1 = (jseed / m2 + j0 * k1 + j1 * k0) % (m2 / 2);
+		j0 = jseed % m2;
+		m[iloop] = j0 + m2 * j1;
+    }
+    i = 4;
+    j = 16;
+
+  }
+
+}
diff --git a/testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/SOR.java b/testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/SOR.java
new file mode 100644
index 0000000..a5cd61f
--- /dev/null
+++ b/testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/SOR.java
@@ -0,0 +1,41 @@
+package edu.columbia.cs.psl.invivo.wallace.numericalbench;
+
+public class SOR
+{
+	public static final double num_flops(int M, int N, int num_iterations)
+	{
+		double Md = (double) M;
+		double Nd = (double) N;
+		double num_iterD = (double) num_iterations;
+
+		return (Md-1)*(Nd-1)*num_iterD*6.0;
+	}
+
+	public static final void execute(double omega, double G[][], int
+			num_iterations)
+	{
+		int M = G.length;
+		int N = G[0].length;
+
+		double omega_over_four = omega * 0.25;
+		double one_minus_omega = 1.0 - omega;
+
+		// update interior points
+		//
+		int Mm1 = M-1;
+		int Nm1 = N-1;
+		for (int p=0; p<num_iterations; p++)
+		{
+			for (int i=1; i<Mm1; i++)
+			{
+				double[] Gi = G[i];
+				double[] Gim1 = G[i-1];
+				double[] Gip1 = G[i+1];
+				for (int j=1; j<Nm1; j++)
+					Gi[j] = omega_over_four * (Gim1[j] + Gip1[j] + Gi[j-1]
+								+ Gi[j+1]) + one_minus_omega * Gi[j];
+			}
+		}
+	}
+}
+
diff --git a/testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/SparseCompRow.java b/testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/SparseCompRow.java
new file mode 100644
index 0000000..f4db878
--- /dev/null
+++ b/testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/SparseCompRow.java
@@ -0,0 +1,46 @@
+package edu.columbia.cs.psl.invivo.wallace.numericalbench;
+
+public class SparseCompRow
+{
+	/* multiple iterations used to make kernel have roughly
+		same granulairty as other Scimark kernels. */
+
+	public static double num_flops(int N, int nz, int num_iterations)
+	{
+		/* Note that if nz does not divide N evenly, then the
+		   actual number of nonzeros used is adjusted slightly.
+		*/
+		int actual_nz = (nz/N) * N;
+		return ((double)actual_nz) * 2.0 * ((double) num_iterations);
+	}
+
+
+	/* computes  a matrix-vector multiply with a sparse matrix
+		held in compress-row format.  If the size of the matrix
+		in MxN with nz nonzeros, then the val[] is the nz nonzeros,
+		with its ith entry in column col[i].  The integer vector row[]
+		is of size M+1 and row[i] points to the begining of the
+		ith row in col[].
+	*/
+
+	public static void matmult( double y[], double val[], int row[],
+		int col[], double x[], int NUM_ITERATIONS)
+	{
+		int M = row.length - 1;
+
+		for (int reps=0; reps<NUM_ITERATIONS; reps++)
+		{
+
+			for (int r=0; r<M; r++)
+			{
+				double sum = 0.0;
+				int rowR = row[r];
+				int rowRp1 = row[r+1];
+				for (int i=rowR; i<rowRp1; i++)
+					sum += x[ col[i] ] * val[i];
+				y[r] = sum;
+			}
+		}
+	}
+
+}
diff --git a/testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/Stopwatch.java b/testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/Stopwatch.java
new file mode 100644
index 0000000..54ff9d3
--- /dev/null
+++ b/testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/Stopwatch.java
@@ -0,0 +1,121 @@
+package edu.columbia.cs.psl.invivo.wallace.numericalbench;
+
+/**
+
+	Provides a stopwatch to measure elapsed time.
+
+<P>
+<DL>
+<DT><B>Example of use:</B></DT>
+<DD>
+<p>
+<pre>
+	Stopwatch Q = new Stopwatch;
+<p>
+	Q.start();
+	//
+	// code to be timed here ...
+	//
+	Q.stop();
+	System.out.println("elapsed time was: " + Q.read() + " seconds.");
+</pre>
+
+@author Roldan Pozo
+@version 14 October 1997, revised 1999-04-24
+*/
+public class Stopwatch
+{
+    private boolean running;
+    private double last_time;
+	private double total;
+
+
+/**
+	Return system time (in seconds)
+
+*/
+	public final static double seconds()
+	{
+		return (System.currentTimeMillis() * 0.001);
+	}
+
+/**
+	Return system time (in seconds)
+
+*/
+	public void reset()
+	{
+		running = false;
+		last_time = 0.0;
+		total=0.0;
+	}
+
+    public Stopwatch()
+	{
+		reset();
+	}
+
+
+/**
+	Start (and reset) timer
+
+*/
+  	public  void start()
+	{
+		if (!running)
+		{
+			running = true;
+			total = 0.0;
+			last_time = seconds();
+		}
+	}
+
+/**
+	Resume timing, after stopping.  (Does not wipe out
+		accumulated times.)
+
+*/
+  	public  void resume()
+	{
+		if (!running)
+		{
+			last_time = seconds();
+			running = true;
+		}
+	}
+
+
+/**
+	Stop timer
+
+*/
+   public  double stop()
+	{
+		if (running)
+        {
+			total += seconds() - last_time;
+            running = false;
+        }
+        return total;
+    }
+
+
+/**
+	Display the elapsed time (in seconds)
+
+*/
+   public  double read()
+	{
+		if (running)
+       	{
+			total += seconds() - last_time;
+			last_time = seconds();
+		}
+        return total;
+	}
+
+}
+
+
+
+
diff --git a/testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/commandline.java b/testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/commandline.java
new file mode 100644
index 0000000..03b7749
--- /dev/null
+++ b/testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/commandline.java
@@ -0,0 +1,120 @@
+package edu.columbia.cs.psl.invivo.wallace.numericalbench;
+
+import java.util.Properties;
+
+/**
+	SciMark2: A Java numerical benchmark measuring performance
+	of computational kernels for FFTs, Monte Carlo simulation,
+	sparse matrix computations, Jacobi SOR, and dense LU matrix
+	factorizations.
+
+
+*/
+
+
+public class commandline
+{
+
+  /* Benchmark 5 kernels with individual Mflops.
+	 "results[0]" has the average Mflop rate.
+
+  */
+
+
+	public static void main(String args[])
+	{
+		// default to the (small) cache-contained version
+
+		double min_time = Constants.RESOLUTION_DEFAULT;
+
+		int FFT_size = Constants.FFT_SIZE;
+		int SOR_size =  Constants.SOR_SIZE;
+		int Sparse_size_M = Constants.SPARSE_SIZE_M;
+		int Sparse_size_nz = Constants.SPARSE_SIZE_nz;
+		int LU_size = Constants.LU_SIZE;
+
+		// look for runtime options
+
+        if (args.length > 0)
+        {
+
+			if (args[0].equalsIgnoreCase("-h") ||
+						args[0].equalsIgnoreCase("-help"))
+			{
+				System.out.println("Usage: [-large] [minimum_time]");
+				return;
+			}
+
+			int current_arg = 0;
+			if (args[current_arg].equalsIgnoreCase("-large"))
+			{
+				FFT_size = Constants.LG_FFT_SIZE;
+				SOR_size =  Constants.LG_SOR_SIZE;
+				Sparse_size_M = Constants.LG_SPARSE_SIZE_M;
+				Sparse_size_nz = Constants.LG_SPARSE_SIZE_nz;
+				LU_size = Constants.LG_LU_SIZE;
+
+				current_arg++;
+			}
+
+			if (args.length > current_arg)
+        		min_time = Double.valueOf(args[current_arg]).doubleValue();
+        }
+
+
+		// run the benchmark
+
+		double res[] = new double[6];
+		Random R = new Random(Constants.RANDOM_SEED);
+
+		res[1] = kernel.measureFFT( FFT_size, min_time, R);
+		res[2] = kernel.measureSOR( SOR_size, min_time, R);
+		res[3] = kernel.measureMonteCarlo(min_time, R);
+		res[4] = kernel.measureSparseMatmult( Sparse_size_M,
+					Sparse_size_nz, min_time, R);
+		res[5] = kernel.measureLU( LU_size, min_time, R);
+
+
+		res[0] = (res[1] + res[2] + res[3] + res[4] + res[5]) / 5;
+
+
+	    // print out results
+
+		System.out.println();
+		System.out.println("SciMark 2.0a");
+		System.out.println();
+		System.out.println("Composite Score: " + res[0]);
+		System.out.print("FFT ("+FFT_size+"): ");
+		if (res[1]==0.0)
+			System.out.println(" ERROR, INVALID NUMERICAL RESULT!");
+		else
+			System.out.println(res[1]);
+
+		System.out.println("SOR ("+SOR_size+"x"+ SOR_size+"): "
+				+ "  " + res[2]);
+		System.out.println("Monte Carlo : " + res[3]);
+		System.out.println("Sparse matmult (N="+ Sparse_size_M+
+				", nz=" + Sparse_size_nz + "): " + res[4]);
+		System.out.print("LU (" + LU_size + "x" + LU_size + "): ");
+		if (res[5]==0.0)
+			System.out.println(" ERROR, INVALID NUMERICAL RESULT!");
+		else
+			System.out.println(res[5]);
+
+		// print out System info
+		System.out.println();
+		System.out.println("java.vendor: " +
+				System.getProperty("java.vendor"));
+		System.out.println("java.version: " +
+				System.getProperty("java.version"));
+		System.out.println("os.arch: " +
+				System.getProperty("os.arch"));
+		System.out.println("os.name: " +
+				System.getProperty("os.name"));
+		System.out.println("os.version: " +
+				System.getProperty("os.version"));
+
+
+	}
+
+}
diff --git a/testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/kernel.java b/testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/kernel.java
new file mode 100644
index 0000000..50d3dd2
--- /dev/null
+++ b/testcase-generation-tester/src/edu/columbia/cs/psl/invivo/wallace/numericalbench/kernel.java
@@ -0,0 +1,292 @@
+package edu.columbia.cs.psl.invivo.wallace.numericalbench;
+
+public class kernel
+{
+	// each measurement returns approx Mflops
+
+
+	public static double measureFFT(int N, double mintime, Random R)
+	{
+		// initialize FFT data as complex (N real/img pairs)
+
+		double x[] = RandomVector(2*N, R);
+		double oldx[] = NewVectorCopy(x);
+		long cycles = 1;
+		Stopwatch Q = new Stopwatch();
+
+		while(true)
+		{
+			Q.start();
+			for (int i=0; i<cycles; i++)
+			{
+				FFT.transform(x);	// forward transform
+				FFT.inverse(x);		// backward transform
+			}
+			Q.stop();
+			if (Q.read() >= mintime)
+				break;
+
+			cycles *= 2;
+		}
+		// approx Mflops
+
+		final double EPS = 1.0e-10;
+		if ( FFT.test(x) / N > EPS )
+			return 0.0;
+
+		return FFT.num_flops(N)*cycles/ Q.read() * 1.0e-6;
+	}
+
+
+	public static double measureSOR(int N, double min_time, Random R)
+	{
+		double G[][] = RandomMatrix(N, N, R);
+
+		Stopwatch Q = new Stopwatch();
+		int cycles=1;
+		while(true)
+		{
+			Q.start();
+			SOR.execute(1.25, G, cycles);
+			Q.stop();
+			if (Q.read() >= min_time) break;
+
+			cycles *= 2;
+		}
+		// approx Mflops
+		return SOR.num_flops(N, N, cycles) / Q.read() * 1.0e-6;
+	}
+
+	public static double measureMonteCarlo(double min_time, Random R)
+	{
+		Stopwatch Q = new Stopwatch();
+
+		int cycles=1;
+		while(true)
+		{
+			Q.start();
+			MonteCarlo.integrate(cycles);
+			Q.stop();
+			if (Q.read() >= min_time) break;
+
+			cycles *= 2;
+		}
+		// approx Mflops
+		return MonteCarlo.num_flops(cycles) / Q.read() * 1.0e-6;
+	}
+
+
+	public static double measureSparseMatmult(int N, int nz,
+			double min_time, Random R)
+	{
+		// initialize vector multipliers and storage for result
+		// y = A*y;
+
+		double x[] = RandomVector(N, R);
+		double y[] = new double[N];
+
+		// initialize square sparse matrix
+		//
+		// for this test, we create a sparse matrix wit M/nz nonzeros
+		// per row, with spaced-out evenly between the begining of the
+		// row to the main diagonal.  Thus, the resulting pattern looks
+		// like
+		//             +-----------------+
+		//             +*                +
+		//             +***              +
+		//             +* * *            +
+		//             +** *  *          +
+		//             +**  *   *        +
+		//             +* *   *   *      +
+		//             +*  *   *    *    +
+		//             +*   *    *    *  +
+		//             +-----------------+
+		//
+		// (as best reproducible with integer artihmetic)
+		// Note that the first nr rows will have elements past
+		// the diagonal.
+
+		int nr = nz/N; 		// average number of nonzeros per row
+		int anz = nr *N;   // _actual_ number of nonzeros
+
+
+		double val[] = RandomVector(anz, R);
+		int col[] = new int[anz];
+		int row[] = new int[N+1];
+
+		row[0] = 0;
+		for (int r=0; r<N; r++)
+		{
+			// initialize elements for row r
+
+			int rowr = row[r];
+			row[r+1] = rowr + nr;
+			int step = r/ nr;
+			if (step < 1) step = 1;   // take at least unit steps
+
+
+			for (int i=0; i<nr; i++)
+				col[rowr+i] = i*step;
+
+		}
+
+		Stopwatch Q = new Stopwatch();
+
+		int cycles=1;
+		while(true)
+		{
+			Q.start();
+			SparseCompRow.matmult(y, val, row, col, x, cycles);
+			Q.stop();
+			if (Q.read() >= min_time) break;
+
+			cycles *= 2;
+		}
+		// approx Mflops
+		return SparseCompRow.num_flops(N, nz, cycles) / Q.read() * 1.0e-6;
+	}
+
+
+	public static double measureLU(int N, double min_time, Random R)
+	{
+		// compute approx Mlfops, or O if LU yields large errors
+
+		double A[][] = RandomMatrix(N, N,  R);
+		double lu[][] = new double[N][N];
+		int pivot[] = new int[N];
+
+		Stopwatch Q = new Stopwatch();
+
+		int cycles=1;
+		while(true)
+		{
+			Q.start();
+			for (int i=0; i<cycles; i++)
+			{
+				CopyMatrix(lu, A);
+				LU.factor(lu, pivot);
+			}
+			Q.stop();
+			if (Q.read() >= min_time) break;
+
+			cycles *= 2;
+		}
+
+
+		// verify that LU is correct
+		double b[] = RandomVector(N, R);
+		double x[] = NewVectorCopy(b);
+
+		LU.solve(lu, pivot, x);
+
+		final double EPS = 1.0e-12;
+		if ( normabs(b, matvec(A,x)) / N > EPS )
+			return 0.0;
+
+
+		// else return approx Mflops
+		//
+		return LU.num_flops(N) * cycles / Q.read() * 1.0e-6;
+	}
+
+
+  private static double[] NewVectorCopy(double x[])
+  {
+		int N = x.length;
+
+		double y[] = new double[N];
+		for (int i=0; i<N; i++)
+			y[i] = x[i];
+
+		return y;
+  }
+
+  private static void CopyVector(double B[], double A[])
+  {
+		int N = A.length;
+
+		for (int i=0; i<N; i++)
+			B[i] = A[i];
+  }
+
+
+  private static double normabs(double x[], double y[])
+  {
+		int N = x.length;
+		double sum = 0.0;
+
+		for (int i=0; i<N; i++)
+			sum += Math.abs(x[i]-y[i]);
+
+		return sum;
+  }
+
+  private static void CopyMatrix(double B[][], double A[][])
+  {
+        int M = A.length;
+        int N = A[0].length;
+
+		int remainder = N & 3;		 // N mod 4;
+
+        for (int i=0; i<M; i++)
+        {
+            double Bi[] = B[i];
+            double Ai[] = A[i];
+			for (int j=0; j<remainder; j++)
+                Bi[j] = Ai[j];
+            for (int j=remainder; j<N; j+=4)
+			{
+				Bi[j] = Ai[j];
+				Bi[j+1] = Ai[j+1];
+				Bi[j+2] = Ai[j+2];
+				Bi[j+3] = Ai[j+3];
+			}
+		}
+  }
+
+  private static double[][] RandomMatrix(int M, int N, Random R)
+  {
+  		double A[][] = new double[M][N];
+
+        for (int i=0; i<N; i++)
+			for (int j=0; j<N; j++)
+            	A[i][j] = R.nextDouble();
+		return A;
+	}
+
+	private static double[] RandomVector(int N, Random R)
+	{
+		double A[] = new double[N];
+
+		for (int i=0; i<N; i++)
+			A[i] = R.nextDouble();
+		return A;
+	}
+
+	private static double[] matvec(double A[][], double x[])
+	{
+		int N = x.length;
+		double y[] = new double[N];
+
+		matvec(A, x, y);
+
+		return y;
+	}
+
+	private static void matvec(double A[][], double x[], double y[])
+	{
+		int M = A.length;
+		int N = A[0].length;
+
+		for (int i=0; i<M; i++)
+		{
+			double sum = 0.0;
+			double Ai[] = A[i];
+			for (int j=0; j<N; j++)
+				sum += Ai[j] * x[j];
+
+			y[i] = sum;
+		}
+	}
+
+}