From 6ea9939e0610a809f6f47d13ec68df00d1ca0afc Mon Sep 17 00:00:00 2001
From: Joshua Bakita <jbakita@cs.unc.edu>
Date: Fri, 16 Oct 2020 16:55:14 -0400
Subject: Move the DIS benchmarks up a directory and update hardcoded paths

Note that this repo does not attempt to keep a copy of the original
DIS benchmark distributions. UNC real-time has another repo for that.
---
 dis/Matrix/ver1/DISstressmarkRNG.h | 190 ++++++++++++
 dis/Matrix/ver1/matrix.c           | 600 +++++++++++++++++++++++++++++++++++++
 dis/Matrix/ver2/DISstressmarkRNG.h | 190 ++++++++++++
 dis/Matrix/ver2/matrix.c           | 594 ++++++++++++++++++++++++++++++++++++
 4 files changed, 1574 insertions(+)
 create mode 100755 dis/Matrix/ver1/DISstressmarkRNG.h
 create mode 100755 dis/Matrix/ver1/matrix.c
 create mode 100755 dis/Matrix/ver2/DISstressmarkRNG.h
 create mode 100755 dis/Matrix/ver2/matrix.c

(limited to 'dis/Matrix')

diff --git a/dis/Matrix/ver1/DISstressmarkRNG.h b/dis/Matrix/ver1/DISstressmarkRNG.h
new file mode 100755
index 0000000..4aa2620
--- /dev/null
+++ b/dis/Matrix/ver1/DISstressmarkRNG.h
@@ -0,0 +1,190 @@
+#include <math.h>
+
+#define IA 16807
+#define IM 2147483647
+#define AM (1.0/IM)
+#define IQ 127773
+#define IR 2836
+#define NTAB 32
+#define NDIV (1+(IM-1)/NTAB)
+#define EPS 1.2e-7
+#define RNMX (1.0-EPS)
+
+static long iy=0;
+static long iv[NTAB];
+static long iseed;
+
+int ABS(int x){
+  if (x>= 0) return x;
+  else
+    return (-x);
+}
+
+int sign(int x){
+  if (x >= 0) return 1;
+  else
+    return (-1);
+}
+
+int MAX(int x, int y){
+  if (x>= y) return x;
+  else
+    return y;
+}
+
+int MIN(int x, int y){
+  if (x<= y) return x;
+  else
+    return y;
+}
+
+void randInit(long idum)
+{
+  long j;
+  long k;
+
+  assert (idum <= 0);
+  assert (iy == 0);
+
+  iseed = idum;
+  if (-(iseed)<1){
+    iseed = 1;
+  }
+  else {
+    iseed = -(iseed);
+  }
+  for (j=NTAB+7; j>=0; j--){
+    k = (iseed)/IQ;
+    iseed = IA*(iseed-k*IQ)-IR*k;
+    if (iseed < 0){
+      iseed += IM;
+    }
+    if (j < NTAB){
+      iv[j] = iseed;
+    }
+  }
+  iy = iv[0];
+}
+
+float randNum()
+{
+  long j;
+  long k;
+  float temp;
+
+  assert (iy != 0);
+  
+  k = (iseed)/IQ;
+  iseed = IA*(iseed-k*IQ)-IR*k;
+
+  if (iseed < 0){
+    iseed += IM;
+  }
+  j = iy/NDIV; 
+  iy = iv[j];
+  iv[j] = iseed;
+
+  temp = AM * iy;
+
+  if (temp > RNMX){
+    return RNMX;
+  }
+  else {
+    return temp;
+  }
+}
+
+
+float randomFloat(float lowest_float, float highest_float)
+{
+ float value;
+ float range;
+
+assert (lowest_float < highest_float);
+
+range = highest_float - lowest_float;
+value = randNum()*(highest_float - lowest_float) + lowest_float;
+assert(value >= lowest_float);
+assert(value <= highest_float);
+
+return value;
+
+}
+
+float randomNonZeroFloat(float lowest_float, float highest_float, float epsilon)
+{
+
+  double range;
+  float value;
+
+
+  assert (lowest_float < 0);
+  assert (highest_float > 0);
+  assert (epsilon > 0);
+  assert ((epsilon < -lowest_float) && (epsilon < highest_float));
+  
+  range = highest_float - lowest_float;
+  value = (randNum() * range)+lowest_float;
+  
+  if (ABS(value) < epsilon) 
+    {
+      if (value > 0) value = value + epsilon;
+      else if (value < 0) value = value - epsilon;
+      
+    }
+
+  assert (value >= lowest_float);
+  assert (value <= highest_float);
+
+  return value;
+}
+
+unsigned int randomUInt(int lowest_uint, int highest_uint)
+{
+  float range;
+  unsigned int value;
+  float temp;
+
+  range =(float)(highest_uint - lowest_uint + 1);
+  temp = randNum();
+  value =(unsigned int)( floor(temp * range) + lowest_uint); 
+
+  assert (value >= lowest_uint);
+  assert (value <= highest_uint);
+
+  return value;
+}
+
+unsigned int randomNonZeroUInt(int lowest_uint, int highest_uint)
+{
+  float range;
+  unsigned int value;
+  float temp;
+
+  range =(float)(highest_uint - lowest_uint + 1);
+  value = 0;
+  while(value == 0){
+  temp = randNum();
+
+  value =(unsigned int)( floor(temp * range) + lowest_uint); 
+  }
+  
+  assert (value >= lowest_uint);
+  assert (value <= highest_uint);
+
+  return value;
+}
+
+int randInt(int lowest_uint,  int highest_uint)
+{
+  float range;
+  int value;
+  
+  range = highest_uint - lowest_uint + 1;
+  value = (int)(floor(randNum() * range) + lowest_uint);
+  
+  assert (value >= lowest_uint);
+  assert (value <= highest_uint);
+
+  return value;
+}			     
diff --git a/dis/Matrix/ver1/matrix.c b/dis/Matrix/ver1/matrix.c
new file mode 100755
index 0000000..518a638
--- /dev/null
+++ b/dis/Matrix/ver1/matrix.c
@@ -0,0 +1,600 @@
+/*
+ *  Sample code for the DIS Matrix Stressmark
+ *
+ * This source code is the completely correct source code based on 
+ * the example codes provided by Atlantic Aerospace Division, Titan 
+ * Systems Corporation, 2000.
+ * 
+ * If you just compile and generate the executables from this source 
+ * code, this code would be enough. However, if you wish to get a complete 
+ * understanding of this stressmark, it is strongly suggested that you
+ * read the Benchmark Analysis and Specifications Document Version 1.0
+ * before going on since the detailed comments are given in this documents.
+ * the comments are not repeated here.
+ */
+
+/*
+ *  The Sparse Matrix Storage is implemented by Compact Row Storage Scheme
+ *  In the code, the data is first generated by randomNonzeroFloat()
+ *  the data is first stored in a full-space matrix with size of dim*dim
+ *  then the data is transfered to the Compact Row Matrix, 
+ *  the data value is kept in *value,
+ *  the columns corresponding to the value are stored in *col_ind,
+ *  the start element of each row is stored in *row_start.
+ */
+ 
+/* 
+ * Please note: 
+ * the total number of data is numberNonzero +dim
+ * among which, NumberNonzero because this is symmetric matrix
+ * dim because the diagonal elements
+ */
+
+#include <stdio.h>
+#include <math.h>
+#include <stdlib.h>
+#include <time.h>
+#include <assert.h>
+#include "DISstressmarkRNG.h"
+
+#define MIN_SEED -2147483647
+#define MAX_SEED -1
+#define MIN_DIM  1
+#define MAX_DIM  32768
+#define MAX_ITERATIONS 65536
+#define MIN_TOLERANCE 0.000007
+#define MAX_TOLERANCE 0.5
+#define MIN_NUMBER   -3.4e10/dim
+#define MAX_NUMBER 3.4e10/dim
+#define EPSI   1.0e-10
+#define MIN_DIG_NUMBER 1.0e-10
+#define MAX_DIG_NUMBER 3.4e10
+
+/*
+ *  External variable, dimension
+ */
+
+static int dim;
+
+/*      
+ *  matrix * vector     
+ */
+
+double  *matrixMulvector(double *value, 
+			   int *col_ind, 
+			   int *row_start,
+			   double  *vector)
+{  
+  int l, ll;
+  double  *out;
+  double  sum;
+  int tmp_rs, tmp_re;
+  
+  out = (double *)malloc(dim*sizeof(double));
+ 
+  for (l=0; l<dim; l++){  
+     *(out + l) = 0;
+     tmp_rs = row_start[l];
+    
+     if (tmp_rs != -1){
+      tmp_re = row_start[l+1];   /*
+				  *get the start and ending elements of 
+				  *  each row
+				 */
+      for (ll=tmp_rs; ll<tmp_re; ll++){
+	*(out + l) += value[ll]*vector[col_ind[ll]];
+      }
+    }
+  }
+  return out;
+}
+
+
+/*
+ *    vector1 - vector2
+ */
+
+double  *vectorSub(double  *vector1, double  *vector2){
+
+  int l;
+  double  *vector;
+  vector = (double  *)malloc(dim*sizeof(double ));
+  for (l=0; l<dim; l++){
+    *(vector + l) = *(vector1 + l) - *(vector2 + l);
+  }
+  return vector;
+}
+
+
+/*
+ * vector1 + vector2
+ */
+
+double  *vectorAdd(double  *vector1, double  *vector2){
+
+  int l;
+  double  *vector;
+
+  vector = (double  *)malloc(dim*sizeof(double ));
+
+  for (l=0; l<dim; l++){
+    *(vector + l) = *(vector1 + l) + *(vector2 + l);
+  }
+  return vector;
+} 
+
+/* 
+ * vector1 * vector2
+ */
+
+double  vectorMul(double  *vector1, double  *vector2){
+
+  int l;
+  double  product;
+
+  product = 0;
+
+  for (l=0; l<dim; l++){
+    product += (*(vector1 + l))*(*(vector2 + l));
+
+  }
+  return product;
+} 
+
+/*
+ * /vector/
+ */
+
+double  vectorValue(double  *vector){
+
+  double  value;
+  int l;
+
+  value = 0;
+
+  for (l=0; l<dim; l++){
+    value += (*(vector + l)) * (*(vector + l));
+  }
+
+  return (sqrt(value));
+}
+
+/*
+ * transpose(vector)
+ * In fact, we return the original vector here
+ */
+
+double  *transpose(double  *vector){
+
+  double  *vect;
+  int l;
+
+  vect = (double  *)malloc(dim*sizeof(double ));
+
+  for (l=0; l<dim; l++){
+    *(vect+l) = *(vector+l);
+  }
+  return vect;
+}
+
+/*
+ * value * <vector>
+ */
+double  *valueMulvector(double  value, double  *vector){
+
+  int l;
+  double  *vect;
+    int lll;
+    double sum;
+
+  vect = (double  *) malloc(dim * sizeof(double ));
+
+  for (l=0; l<dim; l++){
+    *(vect + l) = (*(vector + l)) * value;
+  }
+
+  return vect;
+}
+  
+/*
+ * generate the data distributed sparsely in matrix
+ */
+
+void initMatrix(double  *matrix, int dim, int numberNonzero){
+  
+  int k, l, ll;
+  int i, j;
+
+  int lll;
+  double sum;
+
+  for (k=0; k< dim*dim; k++){
+    *(matrix + k) = 0;
+  }
+
+  for (l=0; l<numberNonzero/2; l++){
+
+    i = randomUInt(1, dim-1);
+    j = randomUInt(0, i-1);
+
+    while (*(matrix + i*dim + j) != 0){
+      
+     i++;
+       if (i == dim){
+       j++;
+       if (j == dim-1){
+	 j = 0;
+	 i = 1;
+       }
+       else{
+	 i = j+1;
+       }
+     }
+    }
+  
+    if (*(matrix + i*dim + j) == 0){
+      *(matrix + i*dim + j) = (double )randomNonZeroFloat(MIN_NUMBER, 
+							  MAX_NUMBER, 
+							  EPSI);
+      *(matrix + j*dim + i) = *(matrix + i*dim + j);
+    }
+  }
+ 
+  for (ll=0; ll<dim; ll++){
+
+
+    
+    *(matrix + ll*dim + ll) = (double )randomNonZeroFloat(-MAX_DIG_NUMBER,
+							  MAX_DIG_NUMBER, 
+							  MIN_DIG_NUMBER);
+    
+    sum = 0;
+
+    for (lll=0; lll<dim; lll++){
+      if (lll != ll){
+	sum += *(matrix + lll*dim + ll);
+      }
+    }
+    
+    if (*(matrix + ll*dim + ll) < sum ){
+      *(matrix + ll*dim + ll) += sum;
+   }
+  }
+
+  return;
+}
+
+/*
+ * generate the data value in the vectors
+ */
+ 
+void initVector(double *vector, int dim){
+
+  int l;
+  
+  for (l=0; l<dim; l++){
+    *(vector + l) = (double )randomFloat (MIN_NUMBER, MAX_NUMBER);
+  }
+
+  return;
+}
+
+/*
+ * make a vector contains value of zero
+ */
+
+void zeroVector(double *vector, int dim){
+  int l;
+  
+  for (l=0; l<dim; l++){
+    *(vector + l) = 0;
+  }
+  return;
+}
+
+/*
+ * return a vector which is the copy of the vect
+ */
+
+double *equalVector(double *vect){
+
+  int l;
+  double *vect1;
+ 
+  vect1 = (double *)malloc(dim*sizeof(double ));
+
+  for (l=0; l<dim; l++){
+    *(vect1+l) = *(vect+l);
+  }
+  return vect1;
+}
+
+
+
+void biConjugateGradient(double *value,
+			 int *col_ind,
+			 int *row_start,
+			 double *vectorB, 
+			 double *vectorX,
+			 double errorTolerance,
+			 int maxIterations,
+			 double *actualError,
+			 int *actualIteration,
+			 int dim)
+     /* 
+      * in the code, we use a lot of temparary vectors and variables
+      * this is just for simple and clear
+      * you can optimize these temporary variables and vectors 
+      * based on your need
+      *
+      */
+{
+  double *vectorR;
+  double *vectorP, *matrixAvectorP, *nextVectorR;
+  double  error;
+  int iteration;
+  double  alpha, beta;
+
+  double  *temp0, *temp1,*temp2, *temp3;
+  double  *temp4, temp5, temp6, *temp7 ;
+  double  *temp8, *temp9, temp10, *temp11;
+  double  temp12, *temp13, *temp14;
+  double  *temp15, *temp16, *temp17;
+
+  int l;
+  int ll;
+
+  alpha = 0;
+  beta = 0;
+
+  /*
+   * vectorR = vectorB - matrixA*vectorX
+   */
+  temp0 = matrixMulvector(value,col_ind, row_start, vectorX);
+  vectorR = vectorSub(vectorB, temp0);
+
+  /*
+   * vectorP = vectorR
+   */
+  vectorP = equalVector(vectorR);
+
+  /*
+   * error = |matrixA * vectorX - vectorB| / |vectorB|
+   */
+  temp2 = vectorSub(temp0, vectorB); 
+  error = vectorValue(temp2)/vectorValue(vectorB);
+
+  free(temp0);
+  free (temp2);
+  
+  iteration = 0;
+ 
+ 
+  while ((iteration < maxIterations) && (error > errorTolerance)){
+   
+    /* 
+     *   alpha = (transpose(vectorR) * vectorR) /
+     *           (transpose(vectorP) * (matrixA * vectorP)
+     */
+    temp1 = matrixMulvector(value, col_ind, row_start, vectorP);    
+    temp3 = transpose(vectorR);       
+    temp4 = transpose(vectorP);
+    temp5 = vectorMul(temp4, temp1);
+    temp6 = vectorMul(temp3, vectorR);
+    alpha = temp6/temp5;
+ 
+    /* 
+     * nextVectorR = vectorR - alpha*(matrixA * vectorP)
+     */
+
+    temp7 = valueMulvector(alpha, temp1); 
+    temp8 = vectorSub(vectorR, temp7);
+    nextVectorR = equalVector(temp8);
+
+    /* 
+     * beta = (transpose(nextVectorR) * nextVectorR) /
+     *           (transpose(vectorR) * vectorR)
+     */
+    temp9 = transpose(nextVectorR);
+    temp10 = vectorMul(temp9, nextVectorR);
+    temp11 = transpose(vectorR);
+    temp12 = vectorMul(temp11, vectorR);
+    beta = temp10/temp12;
+
+    /*
+     * vectorX = vectorX + alpha * vectorP
+     */
+    temp13 = valueMulvector(alpha, vectorP);       
+    vectorX = vectorAdd(vectorX,temp13);
+
+    /* 
+     *vectorP = nextVectorR + beta*vectorP
+     */       
+    temp14 = valueMulvector(beta, vectorP);   
+    temp17 = vectorAdd(nextVectorR, temp14);
+
+    for (ll=0; ll<dim; ll++){
+      *(vectorP + ll) = *(temp17 + ll);
+    }
+
+    /*
+     * vectorR = nextVectorR
+     */
+    
+    for (l=0; l<dim; l++){
+    *(vectorR+l) = *(nextVectorR+l);
+    }
+
+    /* 
+     * error = |matrixA * vectorX - vectorB| / |vectorB|
+     */
+    temp15 = matrixMulvector(value, col_ind,row_start, vectorX);
+    temp16 = vectorSub(temp15,vectorB);
+    error = vectorValue(temp16)/vectorValue(vectorB);
+
+    free(temp1);
+    free(temp3);
+    free(temp4);
+    free(temp7);
+    free(temp8);
+    free(temp9);
+    free(temp11);
+    free(nextVectorR);
+    free(temp13);
+    free(temp14);
+    free(temp15);
+    free(temp16);
+    free(temp17);
+
+    iteration++;
+  }
+
+  *actualError = error;
+  *actualIteration = iteration;
+
+  free(vectorR);
+  free(vectorP);
+
+  return;
+}
+  
+/*
+ * This is the function to transfer the data from the matrix of dense storage 
+ * to Compact Row Storage
+ */
+void create_CRS(double *matrixA,
+		double *value, 
+		int *col_ind, 
+		int *row_start,
+		int dim,
+		int numberNonzero)
+{
+
+  int i, j, k;
+  int cnt;
+  double tmp;
+
+  /* 
+   *initialize the row_start
+   */
+     
+  for(k=0; k<dim; k++){
+    row_start[k] = -1;
+  }
+  
+  /* 
+   * make the end of the last row to be numberNonzero + dim.
+   */
+
+  row_start[dim] = numberNonzero+dim;
+  
+  /*
+   * initialize the col_ind
+   */
+
+  for (k=0; k<numberNonzero+dim; k++){
+    col_ind[k] = -1;
+  }
+
+
+  cnt = 0;
+
+  for (i=0;  (cnt<numberNonzero+dim)&&(i<dim); i++){
+    for (j=0; (cnt<numberNonzero+dim)&&(j<dim); j++){
+      
+      tmp = *(matrixA + i*dim + j);
+
+         if (tmp!=0){
+
+	   value[cnt] = tmp;
+	   col_ind[cnt] = j;
+       
+	    if (row_start[i] == -1)
+	      row_start[i] = cnt;
+	    
+	    cnt += 1;
+	 }	
+    }
+  }
+  row_start[i] = cnt;
+
+  return;
+}
+
+
+
+int main()
+{
+  int seed;
+  int numberNonzero;
+  int maxIterations;
+  float errorTolerance;
+  double  actualError;
+  int actualIteration;
+  
+  time_t beginTime;
+  time_t endTime;
+
+  double  *matrixA;
+  double  *vectorB;
+  double  *vectorX;
+
+  double *value;
+  int    *col_ind;
+  int    *row_start;
+  int sum;
+  int k;
+
+  fscanf(stdin, "%d %d %d %d %f", 
+	 &seed, &dim, &numberNonzero,&maxIterations,&errorTolerance);
+  assert((seed > MIN_SEED) && (seed < MAX_SEED));
+  assert((dim > MIN_DIM) && (dim < MAX_DIM));
+  assert((numberNonzero > dim) && (numberNonzero < dim*dim));
+  assert((maxIterations > 0) && (maxIterations < MAX_ITERATIONS));
+  assert((errorTolerance > MIN_TOLERANCE) && (errorTolerance < MAX_TOLERANCE));
+  
+  matrixA = (double  *)malloc(dim*dim*sizeof(double ));
+  vectorB = (double *)malloc(dim*sizeof(double));
+  vectorX = (double *)malloc(dim*sizeof(double));
+
+  value = (double *)malloc((numberNonzero+dim)*sizeof(double));
+  col_ind = (int *)malloc((numberNonzero+dim)*sizeof(int));
+  row_start = (int *)malloc((dim+1)*sizeof(int));
+
+  randInit(seed);
+
+  initMatrix(matrixA, dim, numberNonzero);
+  
+  create_CRS(matrixA, value, col_ind, row_start, dim, numberNonzero);
+
+  initVector(vectorB, dim);
+  zeroVector(vectorX, dim);
+  printf(" after init\n");
+
+  beginTime = time(NULL);
+  
+  actualError = 0;
+  actualIteration = 0;
+   
+  biConjugateGradient(value, col_ind, row_start, vectorB, vectorX, errorTolerance,
+		      maxIterations,
+		      &actualError, &actualIteration, dim);
+  
+
+
+  endTime = time(NULL) - beginTime;
+  
+
+
+  sum = 0;
+  for (k=1; k<dim; k++){
+    sum += sum + *(vectorX + k);
+  }
+  
+  fprintf(stdout, "sum = %d, actualError = %e, actualIteration = %d\n", sum, actualError, actualIteration);
+  fprintf(stdout, "total time = %u sec. \n", (unsigned int)endTime);
+
+  return(0);
+    }
+
+    
diff --git a/dis/Matrix/ver2/DISstressmarkRNG.h b/dis/Matrix/ver2/DISstressmarkRNG.h
new file mode 100755
index 0000000..4aa2620
--- /dev/null
+++ b/dis/Matrix/ver2/DISstressmarkRNG.h
@@ -0,0 +1,190 @@
+#include <math.h>
+
+#define IA 16807
+#define IM 2147483647
+#define AM (1.0/IM)
+#define IQ 127773
+#define IR 2836
+#define NTAB 32
+#define NDIV (1+(IM-1)/NTAB)
+#define EPS 1.2e-7
+#define RNMX (1.0-EPS)
+
+static long iy=0;
+static long iv[NTAB];
+static long iseed;
+
+int ABS(int x){
+  if (x>= 0) return x;
+  else
+    return (-x);
+}
+
+int sign(int x){
+  if (x >= 0) return 1;
+  else
+    return (-1);
+}
+
+int MAX(int x, int y){
+  if (x>= y) return x;
+  else
+    return y;
+}
+
+int MIN(int x, int y){
+  if (x<= y) return x;
+  else
+    return y;
+}
+
+void randInit(long idum)
+{
+  long j;
+  long k;
+
+  assert (idum <= 0);
+  assert (iy == 0);
+
+  iseed = idum;
+  if (-(iseed)<1){
+    iseed = 1;
+  }
+  else {
+    iseed = -(iseed);
+  }
+  for (j=NTAB+7; j>=0; j--){
+    k = (iseed)/IQ;
+    iseed = IA*(iseed-k*IQ)-IR*k;
+    if (iseed < 0){
+      iseed += IM;
+    }
+    if (j < NTAB){
+      iv[j] = iseed;
+    }
+  }
+  iy = iv[0];
+}
+
+float randNum()
+{
+  long j;
+  long k;
+  float temp;
+
+  assert (iy != 0);
+  
+  k = (iseed)/IQ;
+  iseed = IA*(iseed-k*IQ)-IR*k;
+
+  if (iseed < 0){
+    iseed += IM;
+  }
+  j = iy/NDIV; 
+  iy = iv[j];
+  iv[j] = iseed;
+
+  temp = AM * iy;
+
+  if (temp > RNMX){
+    return RNMX;
+  }
+  else {
+    return temp;
+  }
+}
+
+
+float randomFloat(float lowest_float, float highest_float)
+{
+ float value;
+ float range;
+
+assert (lowest_float < highest_float);
+
+range = highest_float - lowest_float;
+value = randNum()*(highest_float - lowest_float) + lowest_float;
+assert(value >= lowest_float);
+assert(value <= highest_float);
+
+return value;
+
+}
+
+float randomNonZeroFloat(float lowest_float, float highest_float, float epsilon)
+{
+
+  double range;
+  float value;
+
+
+  assert (lowest_float < 0);
+  assert (highest_float > 0);
+  assert (epsilon > 0);
+  assert ((epsilon < -lowest_float) && (epsilon < highest_float));
+  
+  range = highest_float - lowest_float;
+  value = (randNum() * range)+lowest_float;
+  
+  if (ABS(value) < epsilon) 
+    {
+      if (value > 0) value = value + epsilon;
+      else if (value < 0) value = value - epsilon;
+      
+    }
+
+  assert (value >= lowest_float);
+  assert (value <= highest_float);
+
+  return value;
+}
+
+unsigned int randomUInt(int lowest_uint, int highest_uint)
+{
+  float range;
+  unsigned int value;
+  float temp;
+
+  range =(float)(highest_uint - lowest_uint + 1);
+  temp = randNum();
+  value =(unsigned int)( floor(temp * range) + lowest_uint); 
+
+  assert (value >= lowest_uint);
+  assert (value <= highest_uint);
+
+  return value;
+}
+
+unsigned int randomNonZeroUInt(int lowest_uint, int highest_uint)
+{
+  float range;
+  unsigned int value;
+  float temp;
+
+  range =(float)(highest_uint - lowest_uint + 1);
+  value = 0;
+  while(value == 0){
+  temp = randNum();
+
+  value =(unsigned int)( floor(temp * range) + lowest_uint); 
+  }
+  
+  assert (value >= lowest_uint);
+  assert (value <= highest_uint);
+
+  return value;
+}
+
+int randInt(int lowest_uint,  int highest_uint)
+{
+  float range;
+  int value;
+  
+  range = highest_uint - lowest_uint + 1;
+  value = (int)(floor(randNum() * range) + lowest_uint);
+  
+  assert (value >= lowest_uint);
+  assert (value <= highest_uint);
+
+  return value;
+}			     
diff --git a/dis/Matrix/ver2/matrix.c b/dis/Matrix/ver2/matrix.c
new file mode 100755
index 0000000..957d7c5
--- /dev/null
+++ b/dis/Matrix/ver2/matrix.c
@@ -0,0 +1,594 @@
+/* Please note:
+ * This code is the optimized version of the first version of Matrix 
+ * Stressmark. It uses less temporary vectors and vsariables, thus reduce 
+ * memory allocation/deallocation overhead. the simulation is faster 
+ */
+/*
+ *  Sample code for the DIS Matrix Stressmark
+ *
+ * This source code is the completely correct source code based on 
+ * the example codes provided by Atlantic Aerospace Division, Titan 
+ * Systems Corporation, 2000.
+ * 
+ * If you just compile and generate the executables from this source 
+ * code, this code would be enough. However, if you wish to get a complete 
+ * understanding of this stressmark, it is strongly suggested that you
+ * read the Benchmark Analysis and Specifications Document Version 1.0
+ * before going on since the detailed comments are given in this documents.
+ * the comments are not repeated here.
+ */
+
+/*
+ *  The Sparse Matrix Storage is implemented by Compact Row Storage Scheme
+ *  In the code, the data is first generated by randomNonzeroFloat()
+ *  the data is first stored in a full-space matrix with size of dim*dim
+ *  then the data is transfered to the Compact Row Matrix, 
+ *  the data value is kept in *value,
+ *  the columns corresponding to the value are stored in *col_ind,
+ *  the start element of each row is stored in *row_start.
+ */
+ 
+/* 
+ * Please note: 
+ * the total number of data is numberNonzero +dim
+ * among which, NumberNonzero because this is symmetric matrix
+ * dim because the diagonal elements
+ */
+
+#include <stdio.h>
+#include <math.h>
+#include <stdlib.h>
+#include <time.h>
+#include <assert.h>
+#include "DISstressmarkRNG.h"
+#include "extra.h"
+
+#define MIN_SEED -2147483647
+#define MAX_SEED -1
+#define MIN_DIM  1
+#define MAX_DIM  32768
+#define MAX_ITERATIONS 65536
+#define MIN_TOLERANCE 0.000007
+#define MAX_TOLERANCE 0.5
+#define MIN_NUMBER   -3.4e10/dim
+#define MAX_NUMBER 3.4e10/dim
+#define EPSI   1.0e-10
+#define MIN_DIG_NUMBER 1.0e-10
+#define MAX_DIG_NUMBER 3.4e10
+
+/*
+ *  External variable, dimension
+ */
+
+static int dim;
+int argc;
+char** argv;
+
+/*      
+ *  matrix * vector     
+ */
+
+void matrixMulvector(double *value, 
+		     int *col_ind, 
+		     int *row_start,
+		     double  *vector,
+		     double *out)
+{  
+  int l, ll;
+  double  sum;
+  int tmp_rs, tmp_re;
+ 
+  for (l=0; l<dim; l++){  
+     *(out + l) = 0;
+     tmp_rs = row_start[l];
+    
+     if (tmp_rs != -1){
+      tmp_re = row_start[l+1];   /*
+				  *get the start and ending elements of 
+				  *  each row
+				 */
+      for (ll=tmp_rs; ll<tmp_re; ll++){
+	*(out + l) += value[ll]*vector[col_ind[ll]];
+      }
+    }
+  }
+  return; 
+}
+
+
+/*
+ *    vector1 - vector2
+ */
+
+void  vectorSub(double  *vector1, double  *vector2, double *vector){
+
+  int l;
+
+  for (l=0; l<dim; l++){
+    *(vector + l) = *(vector1 + l) - *(vector2 + l);
+  }
+  return; 
+}
+
+
+/*
+ * vector1 + vector2
+ */
+
+void vectorAdd(double  *vector1, double  *vector2, double *vector){
+
+  int l;
+
+  for (l=0; l<dim; l++){
+    *(vector + l) = *(vector1 + l) + *(vector2 + l);
+  }
+  return; 
+} 
+
+/* 
+ * vector1 * vector2
+ */
+
+double  vectorMul(double  *vector1, double  *vector2){
+
+  int l;
+  double  product;
+
+  product = 0;
+
+  for (l=0; l<dim; l++){
+    product += (*(vector1 + l))*(*(vector2 + l));
+
+  }
+  return product;
+} 
+
+/*
+ * /vector/
+ */
+
+double  vectorValue(double  *vector){
+
+  double  value;
+  int l;
+
+  value = 0;
+
+  for (l=0; l<dim; l++){
+    value += (*(vector + l)) * (*(vector + l));
+  }
+
+  return (sqrt(value));
+}
+
+/*
+ * transpose(vector)
+ * In fact, we return the original vector here
+ */
+
+void  transpose(double  *vector, double *vect){
+
+  int l;
+
+  for (l=0; l<dim; l++){
+    *(vect+l) = *(vector+l);
+  }
+  return; 
+}
+
+/*
+ * value * <vector>
+ */
+void valueMulvector(double  value, double  *vector, double *vect){
+
+  int l;
+  int lll, i;
+  double tmp;
+
+  for (l=0; l<dim; l++){
+    *(vect + l) = *(vector + l) * value;
+  }
+  return;
+}
+  
+/*
+ * generate the data distributed sparsely in matrix
+ */
+
+void initMatrix(double  *matrix, int dim, int numberNonzero){
+  
+  int k, l, ll;
+  int i, j;
+
+  int lll;
+  double sum;
+
+  for (k=0; k< dim*dim; k++){
+    *(matrix + k) = 0;
+  }
+
+  for (l=0; l<numberNonzero/2; l++){
+
+    i = randomUInt(1, dim-1);
+    j = randomUInt(0, i-1);
+
+    while (*(matrix + i*dim + j) != 0){
+      
+     i++;
+       if (i == dim){
+       j++;
+       if (j == dim-1){
+	 j = 0;
+	 i = 1;
+       }
+       else{
+	 i = j+1;
+       }
+     }
+    }
+  
+    if (*(matrix + i*dim + j) == 0){
+      *(matrix + i*dim + j) = (double )randomNonZeroFloat(MIN_NUMBER, 
+							  MAX_NUMBER, 
+							  EPSI);
+      *(matrix + j*dim + i) = *(matrix + i*dim + j);
+    }
+  }
+ 
+  for (ll=0; ll<dim; ll++){
+
+
+    
+    *(matrix + ll*dim + ll) = (double )randomNonZeroFloat(-MAX_DIG_NUMBER,
+							  MAX_DIG_NUMBER, 
+							  MIN_DIG_NUMBER);
+    
+    sum = 0;
+
+    for (lll=0; lll<dim; lll++){
+      if (lll != ll){
+	sum += *(matrix + lll*dim + ll);
+      }
+    }
+    
+    if (*(matrix + ll*dim + ll) < sum ){
+      *(matrix + ll*dim + ll) += sum;
+   }
+  }
+
+  return;
+}
+
+/*
+ * generate the data value in the vectors
+ */
+ 
+void initVector(double *vector, int dim){
+
+  int l;
+  
+  for (l=0; l<dim; l++){
+    *(vector + l) = (double )randomFloat (MIN_NUMBER, MAX_NUMBER);
+  }
+
+  return;
+}
+
+/*
+ * make a vector contains value of zero
+ */
+
+void zeroVector(double *vector, int dim){
+  int l;
+  
+  for (l=0; l<dim; l++){
+    *(vector + l) = 0;
+  }
+  return;
+}
+
+/*
+ * return a vector which is the copy of the vect
+ */
+
+void equalVector(double *vect, double *vect1){
+
+  int l;
+
+  for (l=0; l<dim; l++){
+    *(vect1+l) = *(vect+l);
+  }
+  return; 
+}
+
+
+
+void biConjugateGradient(double *value,
+			 int *col_ind,
+			 int *row_start,
+			 double *vectorB, 
+			 double *vectorX,
+			 double errorTolerance,
+			 int maxIterations,
+			 double *actualError,
+			 int *actualIteration,
+			 int dim)
+     /* 
+      * in the code, we use a lot of temparary vectors and variables
+      * this is just for simple and clear
+      * you can optimize these temporary variables and vectors 
+      * based on your need
+      *
+      */
+{
+  double *vectorR;
+  double *vectorP, *matrixAvectorP, *nextVectorR;
+  double  error;
+  int iteration;
+  double  alpha, beta;
+
+  double  *tmpVector1, *tmpVector2, *tmpVector3;
+  double   tmpValue1, tmpValue2; 
+  int i;
+  int l;
+  int ll;
+  SET_UP
+
+  alpha = 0;
+  beta = 0;
+
+  vectorP = (double *)malloc(dim*sizeof(double));
+  vectorR = (double *)malloc(dim*sizeof(double));
+  nextVectorR = (double *)malloc(dim*sizeof(double));
+  vectorX = (double *)malloc(dim*sizeof(double));
+
+  tmpVector1 = (double *)malloc(dim*sizeof(double));
+  tmpVector2 = (double *)malloc(dim*sizeof(double));
+  tmpVector3 = (double *)malloc(dim*sizeof(double));
+
+  /*
+   * vectorR = vectorB - matrixA*vectorX
+   */
+  matrixMulvector(value,col_ind, row_start, vectorX, tmpVector1);
+
+  vectorSub(vectorB, tmpVector1, vectorR);
+
+  /*
+   * vectorP = vectorR
+   */
+
+  equalVector(vectorR, vectorP);
+
+  /*
+   * error = |matrixA * vectorX - vectorB| / |vectorB|
+   */
+  vectorSub(tmpVector1, vectorB, tmpVector1); 
+
+  error = vectorValue(tmpVector1)/vectorValue(vectorB);
+
+  iteration = 0;
+
+  while ((iteration < maxIterations) && (error > errorTolerance)){
+    START_LOOP
+   
+    /* 
+     *   alpha = (transpose(vectorR) * vectorR) /
+     *           (transpose(vectorP) * (matrixA * vectorP)
+     */
+
+    matrixMulvector(value, col_ind, row_start, vectorP, tmpVector1);  
+    transpose(vectorR, tmpVector2); 
+    transpose(vectorP, tmpVector3);
+    tmpValue1 = vectorMul(tmpVector3, tmpVector1);
+    tmpValue2 = vectorMul(tmpVector2, vectorR);
+    alpha = tmpValue2/tmpValue1;
+ 
+    /* 
+     * nextVectorR = vectorR - alpha*(matrixA * vectorP)
+     */
+
+    valueMulvector(alpha, tmpVector1, tmpVector2);
+    vectorSub(vectorR, tmpVector2, tmpVector1);
+    equalVector(tmpVector1, nextVectorR);
+ 
+    /* 
+     * beta = (transpose(nextVectorR) * nextVectorR) /
+     *           (transpose(vectorR) * vectorR)
+     */
+
+    transpose(nextVectorR, tmpVector3);
+    tmpValue1 = vectorMul(tmpVector3, nextVectorR);
+    transpose(vectorR, tmpVector2);
+    tmpValue2 = vectorMul(tmpVector2, vectorR);
+    beta = tmpValue1/tmpValue2;
+
+    /*
+     * vectorX = vectorX + alpha * vectorP
+     */
+    valueMulvector(alpha, vectorP, tmpVector1);       
+    vectorAdd(vectorX,tmpVector1, vectorX);
+
+    /* 
+     *vectorP = nextVectorR + beta*vectorP
+     */       
+    valueMulvector(beta, vectorP, tmpVector1);   
+    vectorAdd(nextVectorR, tmpVector1, tmpVector1);
+
+    for (ll=0; ll<dim; ll++){
+      *(vectorP + ll) = *(tmpVector1 + ll);
+    }
+
+    /*
+     * vectorR = nextVectorR
+     */
+    
+    for (l=0; l<dim; l++){
+    *(vectorR+l) = *(nextVectorR+l);
+    }
+
+    /* 
+     * error = |matrixA * vectorX - vectorB| / |vectorB|
+     */
+    matrixMulvector(value, col_ind,row_start, vectorX, tmpVector1);
+    vectorSub(tmpVector1,vectorB,tmpVector1);
+    error = vectorValue(tmpVector1)/vectorValue(vectorB);
+
+    iteration++;
+    STOP_LOOP
+  }
+
+  *actualError = error;
+  *actualIteration = iteration;
+
+  free(tmpVector1);
+  free(tmpVector2);
+  free(tmpVector3);
+
+  free(vectorR);
+  free(vectorP);
+  WRITE_TO_FILE
+
+  return;
+}
+  
+/*
+ * This is the function to transfer the data from the matrix of dense storage 
+ * to Compact Row Storage
+ */
+void create_CRS(double *matrixA,
+		double *value, 
+		int *col_ind, 
+		int *row_start,
+		int dim,
+		int numberNonzero)
+{
+
+  int i, j, k;
+  int cnt;
+  double tmp;
+
+  /* 
+   *initialize the row_start
+   */
+     
+  for(k=0; k<dim; k++){
+    row_start[k] = -1;
+  }
+  
+  /* 
+   * make the end of the last row to be numberNonzero + dim.
+   */
+
+  row_start[dim] = numberNonzero+dim;
+  
+  /*
+   * initialize the col_ind
+   */
+
+  for (k=0; k<numberNonzero+dim; k++){
+    col_ind[k] = -1;
+  }
+
+
+  cnt = 0;
+
+  for (i=0;  (cnt<numberNonzero+dim)&&(i<dim); i++){
+    for (j=0; (cnt<numberNonzero+dim)&&(j<dim); j++){
+      
+      tmp = *(matrixA + i*dim + j);
+
+         if (tmp!=0){
+
+	   value[cnt] = tmp;
+	   col_ind[cnt] = j;
+       
+	    if (row_start[i] == -1)
+	      row_start[i] = cnt;
+	    
+	    cnt += 1;
+	 }	
+    }
+  }
+  row_start[i] = cnt;
+
+  return;
+}
+
+
+int main(int _argc, char** _argv)
+{
+  argc = _argc;
+  argv = _argv;
+  int seed;
+  int numberNonzero;
+  int maxIterations;
+  float errorTolerance;
+  double  actualError;
+  int actualIteration;
+  
+  time_t beginTime;
+  time_t endTime;
+
+  double  *matrixA;
+  double  *vectorB;
+  double  *vectorX;
+
+  double *value;
+  int    *col_ind;
+  int    *row_start;
+  int sum;
+  int k;
+
+  fscanf(stdin, "%d %d %d %d %f", 
+	 &seed, &dim, &numberNonzero,&maxIterations,&errorTolerance);
+  assert((seed > MIN_SEED) && (seed < MAX_SEED));
+  assert((dim > MIN_DIM) && (dim < MAX_DIM));
+  assert((numberNonzero > dim) && (numberNonzero < dim*dim));
+  assert((maxIterations > 0) && (maxIterations < MAX_ITERATIONS));
+  assert((errorTolerance > MIN_TOLERANCE) && (errorTolerance < MAX_TOLERANCE));
+  
+  matrixA = (double  *)malloc(dim*dim*sizeof(double ));
+  vectorB = (double *)malloc(dim*sizeof(double));
+  vectorX = (double *)malloc(dim*sizeof(double));
+
+  value = (double *)malloc((numberNonzero+dim)*sizeof(double));
+  col_ind = (int *)malloc((numberNonzero+dim)*sizeof(int));
+  row_start = (int *)malloc((dim+1)*sizeof(int));
+
+  randInit(seed);
+
+  initMatrix(matrixA, dim, numberNonzero);
+  
+  create_CRS(matrixA, value, col_ind, row_start, dim, numberNonzero);
+
+  initVector(vectorB, dim);
+  zeroVector(vectorX, dim);
+  printf(" after init\n");
+
+  beginTime = time(NULL);
+  
+  actualError = 0;
+  actualIteration = 0;
+   
+  biConjugateGradient(value, col_ind, row_start, vectorB, vectorX, errorTolerance,
+		      maxIterations,
+		      &actualError, &actualIteration, dim);
+  
+
+
+  endTime = time(NULL) - beginTime;
+  
+
+
+  sum = 0;
+  for (k=1; k<dim; k++){
+    sum += sum + *(vectorX + k);
+  }
+  
+  fprintf(stdout, "sum = %d, actualError = %e, actualIteration = %d\n", sum, actualError, actualIteration);
+  fprintf(stdout, "total time = %u sec. \n", (unsigned int)endTime);
+
+  return(0);
+    }
+
+    
-- 
cgit v1.2.2