1 files changed, 219 insertions, 249 deletions

diff --git a/dis/Matrix/ver2/matrix.c b/dis/Matrix/ver2/matrix.c
index ffa7cb7..2b075fb 100755
--- a/dis/Matrix/ver2/matrix.c
+++ b/dis/Matrix/ver2/matrix.c
@@ -1,17 +1,17 @@
 /* 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 
+ * 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 
+ * 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 
+ *
+ * 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.
@@ -22,37 +22,37 @@
  *  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, 
+ *  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: 
+
+/*
+ * 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 "DISstressmarkRNG.h"
+#include "extra.h"
+#include <assert.h>
 #include <math.h>
+#include <stdio.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 MIN_DIM 1
+#define MAX_DIM 32768
 #define MAX_ITERATIONS 65536
-#define MIN_TOLERANCE 1e-7//0.000007
+#define MIN_TOLERANCE 1e-7 // 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_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
 
@@ -62,97 +62,90 @@
 
 static int dim;
 
-/*      
- *  matrix * vector     
+/*
+ *  matrix * vector
  */
 
-void matrixMulvector(double *value, 
-                     int *col_ind, 
-                     int *row_start,
-                     double  *vector,
-                     double *out)
-{  
+void matrixMulvector(double *value, int *col_ind, int *row_start,
+                     double *vector, double *out) {
   int l, ll;
-  double  sum;
+  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]];
+
+  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; 
+  return;
 }
 
-
 /*
  *    vector1 - vector2
  */
 
-void  vectorSub(double  *vector1, double  *vector2, double *vector){
+void vectorSub(double *vector1, double *vector2, double *vector) {
 
   int l;
 
-  for (l=0; l<dim; l++){
+  for (l = 0; l < dim; l++) {
     *(vector + l) = *(vector1 + l) - *(vector2 + l);
   }
-  return; 
+  return;
 }
 
-
 /*
  * vector1 + vector2
  */
 
-void vectorAdd(double  *vector1, double  *vector2, double *vector){
+void vectorAdd(double *vector1, double *vector2, double *vector) {
 
   int l;
 
-  for (l=0; l<dim; l++){
+  for (l = 0; l < dim; l++) {
     *(vector + l) = *(vector1 + l) + *(vector2 + l);
   }
-  return; 
-} 
+  return;
+}
 
-/* 
+/*
  * vector1 * vector2
  */
 
-double  vectorMul(double  *vector1, double  *vector2){
+double vectorMul(double *vector1, double *vector2) {
 
   int l;
-  double  product;
+  double product;
 
   product = 0;
 
-  for (l=0; l<dim; l++){
-    product += (*(vector1 + l))*(*(vector2 + l));
-
+  for (l = 0; l < dim; l++) {
+    product += (*(vector1 + l)) * (*(vector2 + l));
   }
   return product;
-} 
+}
 
 /*
  * /vector/
  */
 
-double  vectorValue(double  *vector){
+double vectorValue(double *vector) {
 
-  double  value;
+  double value;
   int l;
 
   value = 0;
 
-  for (l=0; l<dim; l++){
+  for (l = 0; l < dim; l++) {
     value += (*(vector + l)) * (*(vector + l));
   }
 
@@ -164,94 +157,89 @@ double  vectorValue(double  *vector){
  * In fact, we return the original vector here
  */
 
-void  transpose(double  *vector, double *vect){
+void transpose(double *vector, double *vect) {
 
   int l;
 
-  for (l=0; l<dim; l++){
-    *(vect+l) = *(vector+l);
+  for (l = 0; l < dim; l++) {
+    *(vect + l) = *(vector + l);
   }
-  return; 
+  return;
 }
 
 /*
  * value * <vector>
  */
-void valueMulvector(double  value, double  *vector, double *vect){
+void valueMulvector(double value, double *vector, double *vect) {
 
   int l;
   int lll, i;
   double tmp;
 
-  for (l=0; l<dim; l++){
+  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){
-  
+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++){
+  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;
-       }
-     }
+  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);
+
+    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++){
 
+  for (ll = 0; ll < dim; ll++) {
+
+    *(matrix + ll * dim + ll) = (double)randomNonZeroFloat(
+        -MAX_DIG_NUMBER, MAX_DIG_NUMBER, MIN_DIG_NUMBER);
 
-    
-    *(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);
+    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;
-   }
+
+    if (*(matrix + ll * dim + ll) < sum) {
+      *(matrix + ll * dim + ll) += sum;
+    }
   }
 
   return;
@@ -260,13 +248,13 @@ void initMatrix(double  *matrix, int dim, int numberNonzero){
 /*
  * generate the data value in the vectors
  */
- 
-void initVector(double *vector, int dim){
+
+void initVector(double *vector, int dim) {
 
   int l;
-  
-  for (l=0; l<dim; l++){
-    *(vector + l) = (double )randomFloat (MIN_NUMBER, MAX_NUMBER);
+
+  for (l = 0; l < dim; l++) {
+    *(vector + l) = (double)randomFloat(MIN_NUMBER, MAX_NUMBER);
   }
 
   return;
@@ -276,10 +264,10 @@ void initVector(double *vector, int dim){
  * make a vector contains value of zero
  */
 
-void zeroVector(double *vector, int dim){
+void zeroVector(double *vector, int dim) {
   int l;
-  
-  for (l=0; l<dim; l++){
+
+  for (l = 0; l < dim; l++) {
     *(vector + l) = 0;
   }
   return;
@@ -289,48 +277,37 @@ void zeroVector(double *vector, int dim){
  * return a vector which is the copy of the vect
  */
 
-void equalVector(double *vect, double *vect1){
+void equalVector(double *vect, double *vect1) {
 
   int l;
 
-  for (l=0; l<dim; l++){
-    *(vect1+l) = *(vect+l);
+  for (l = 0; l < dim; l++) {
+    *(vect1 + l) = *(vect + l);
   }
-  return; 
+  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,
-                         // Start internal matrixes
-                         double *vectorP,
-                         double *vectorR,
-                         double *nextVectorR,
-                         double *tmpVector1,
-                         double *tmpVector2,
-                         double *tmpVector3)
-     /* 
-      * 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
-      *
-      */
+void biConjugateGradient(double *value, int *col_ind, int *row_start,
+                         double *vectorB, double *vectorX,
+                         double errorTolerance, int maxIterations,
+                         double *actualError, int *actualIteration, int dim,
+                         // Start internal matrixes
+                         double *vectorP, double *vectorR, double *nextVectorR,
+                         double *tmpVector1, double *tmpVector2,
+                         double *tmpVector3)
+/*
+ * 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  error;
+  double error;
   int iteration;
-  double  alpha, beta;
+  double alpha, beta;
 
-  double   tmpValue1, tmpValue2; 
+  double tmpValue1, tmpValue2;
   int i;
   int l;
   int ll;
@@ -341,7 +318,7 @@ void biConjugateGradient(double *value,
   /*
    * vectorR = vectorB - matrixA*vectorX
    */
-  matrixMulvector(value,col_ind, row_start, vectorX, tmpVector1);
+  matrixMulvector(value, col_ind, row_start, vectorX, tmpVector1);
 
   vectorSub(vectorB, tmpVector1, vectorR);
 
@@ -354,35 +331,35 @@ void biConjugateGradient(double *value,
   /*
    * error = |matrixA * vectorX - vectorB| / |vectorB|
    */
-  vectorSub(tmpVector1, vectorB, tmpVector1); 
+  vectorSub(tmpVector1, vectorB, tmpVector1);
 
-  error = vectorValue(tmpVector1)/vectorValue(vectorB);
+  error = vectorValue(tmpVector1) / vectorValue(vectorB);
 
   iteration = 0;
 
-  while ((iteration < maxIterations) && (error > errorTolerance)){
-   
-    /* 
+  while ((iteration < maxIterations) && (error > errorTolerance)) {
+
+    /*
      *   alpha = (transpose(vectorR) * vectorR) /
      *           (transpose(vectorP) * (matrixA * vectorP)
      */
 
-    matrixMulvector(value, col_ind, row_start, vectorP, tmpVector1);  
-    transpose(vectorR, tmpVector2); 
+    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;
- 
-    /* 
+    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)
      */
@@ -391,38 +368,38 @@ void biConjugateGradient(double *value,
     tmpValue1 = vectorMul(tmpVector3, nextVectorR);
     transpose(vectorR, tmpVector2);
     tmpValue2 = vectorMul(tmpVector2, vectorR);
-    beta = tmpValue1/tmpValue2;
+    beta = tmpValue1 / tmpValue2;
 
     /*
      * vectorX = vectorX + alpha * vectorP
      */
-    valueMulvector(alpha, vectorP, tmpVector1);       
-    vectorAdd(vectorX,tmpVector1, vectorX);
+    valueMulvector(alpha, vectorP, tmpVector1);
+    vectorAdd(vectorX, tmpVector1, vectorX);
 
-    /* 
+    /*
      *vectorP = nextVectorR + beta*vectorP
-     */       
-    valueMulvector(beta, vectorP, tmpVector1);   
+     */
+    valueMulvector(beta, vectorP, tmpVector1);
     vectorAdd(nextVectorR, tmpVector1, tmpVector1);
 
-    for (ll=0; ll<dim; ll++){
+    for (ll = 0; ll < dim; ll++) {
       *(vectorP + ll) = *(tmpVector1 + ll);
     }
 
     /*
      * vectorR = nextVectorR
      */
-    
-    for (l=0; l<dim; l++){
-    *(vectorR+l) = *(nextVectorR+l);
+
+    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);
+    matrixMulvector(value, col_ind, row_start, vectorX, tmpVector1);
+    vectorSub(tmpVector1, vectorB, tmpVector1);
+    error = vectorValue(tmpVector1) / vectorValue(vectorB);
 
     iteration++;
   }
@@ -439,63 +416,57 @@ void biConjugateGradient(double *value,
 
   return;
 }
-  
+
 /*
- * This is the function to transfer the data from the matrix of dense storage 
+ * 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)
-{
+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++){
+
+  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;
-  
+  row_start[dim] = numberNonzero + dim;
+
   /*
    * initialize the col_ind
    */
 
-  for (k=0; k<numberNonzero+dim; k++){
+  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;
-         }      
+  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;
@@ -503,56 +474,54 @@ void create_CRS(double *matrixA,
   return;
 }
 
-
-int main(int argc, char** argv)
-{
+int main(int argc, char **argv) {
   int seed;
   int numberNonzero;
   int maxIterations;
   float errorTolerance;
-  double  actualError;
+  double actualError;
   int actualIteration;
-  
+
   time_t beginTime;
   time_t endTime;
 
-  double  *matrixA;
-  double  *vectorB;
-  double  *vectorX;
+  double *matrixA;
+  double *vectorB;
+  double *vectorX;
 
   double *value;
-  int    *col_ind;
-  int    *row_start;
+  int *col_ind;
+  int *row_start;
   double sum;
   int k;
   // Internal vects for biConj
   double *vectorR, *vectorP, *nextVectorR;
-  double  *tmpVector1, *tmpVector2, *tmpVector3;
+  double *tmpVector1, *tmpVector2, *tmpVector3;
   SET_UP
 
-  assert(fscanf(stdin, "%d %d %d %d %f",
-         &seed, &dim, &numberNonzero,&maxIterations,&errorTolerance) == 5);
+  assert(fscanf(stdin, "%d %d %d %d %f", &seed, &dim, &numberNonzero,
+                &maxIterations, &errorTolerance) == 5);
   assert((seed > MIN_SEED) && (seed < MAX_SEED));
   assert((dim > MIN_DIM) && (dim < MAX_DIM));
-  assert((numberNonzero > dim) && (numberNonzero < dim*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));
+  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));
 
   // Internal matricies for biConj
-  vectorP = (double *)malloc(dim*sizeof(double));
-  vectorR = (double *)malloc(dim*sizeof(double));
-  nextVectorR = (double *)malloc(dim*sizeof(double));
-  tmpVector1 = (double *)malloc(dim*sizeof(double));
-  tmpVector2 = (double *)malloc(dim*sizeof(double));
-  tmpVector3 = (double *)malloc(dim*sizeof(double));
+  vectorP = (double *)malloc(dim * sizeof(double));
+  vectorR = (double *)malloc(dim * sizeof(double));
+  nextVectorR = (double *)malloc(dim * sizeof(double));
+  tmpVector1 = (double *)malloc(dim * sizeof(double));
+  tmpVector2 = (double *)malloc(dim * sizeof(double));
+  tmpVector3 = (double *)malloc(dim * sizeof(double));
 
   randInit(seed);
 
@@ -569,23 +538,24 @@ int main(int argc, char** argv)
   actualError = 0;
   actualIteration = 0;
 
-  biConjugateGradient(value, col_ind, row_start, vectorB, vectorX, errorTolerance,
-                      maxIterations,
-                      &actualError, &actualIteration, dim,
-              vectorP, vectorR, nextVectorR, tmpVector1, tmpVector2, tmpVector3);
+  biConjugateGradient(value, col_ind, row_start, vectorB, vectorX,
+                      errorTolerance, maxIterations, &actualError,
+                      &actualIteration, dim, vectorP, vectorR, nextVectorR,
+                      tmpVector1, tmpVector2, tmpVector3);
 
   STOP_LOOP
   endTime = time(NULL);
 
   sum = 0;
-  for (k=1; k<dim; k++){
+  for (k = 1; k < dim; k++) {
     sum += sum + *(vectorX + k);
   }
-  
-  fprintf(stdout, "sum = %f, actualError = %e, actualIteration = %d\n", sum, actualError, actualIteration);
-  fprintf(stderr, "time for matrix stressmark = %f seconds.\n", difftime(endTime, beginTime));
+
+  fprintf(stdout, "sum = %f, actualError = %e, actualIteration = %d\n", sum,
+          actualError, actualIteration);
+  fprintf(stderr, "time for matrix stressmark = %f seconds.\n",
+          difftime(endTime, beginTime));
 
   WRITE_TO_FILE
-  return(0);
+  return (0);
 }
-