diff options
Diffstat (limited to 'dis/Matrix/ver2')
| -rwxr-xr-x | dis/Matrix/ver2/DISstressmarkRNG.h | 1 | ||||
| -rwxr-xr-x | dis/Matrix/ver2/matrix.c | 468 |
2 files changed, 220 insertions, 249 deletions
diff --git a/dis/Matrix/ver2/DISstressmarkRNG.h b/dis/Matrix/ver2/DISstressmarkRNG.h index 4aa2620..6243606 100755 --- a/dis/Matrix/ver2/DISstressmarkRNG.h +++ b/dis/Matrix/ver2/DISstressmarkRNG.h | |||
| @@ -1,3 +1,4 @@ | |||
| 1 | #include <assert.h> | ||
| 1 | #include <math.h> | 2 | #include <math.h> |
| 2 | 3 | ||
| 3 | #define IA 16807 | 4 | #define IA 16807 |
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 @@ | |||
| 1 | /* Please note: | 1 | /* Please note: |
| 2 | * This code is the optimized version of the first version of Matrix | 2 | * This code is the optimized version of the first version of Matrix |
| 3 | * Stressmark. It uses less temporary vectors and vsariables, thus reduce | 3 | * Stressmark. It uses less temporary vectors and vsariables, thus reduce |
| 4 | * memory allocation/deallocation overhead. the simulation is faster | 4 | * memory allocation/deallocation overhead. the simulation is faster |
| 5 | */ | 5 | */ |
| 6 | /* | 6 | /* |
| 7 | * Sample code for the DIS Matrix Stressmark | 7 | * Sample code for the DIS Matrix Stressmark |
| 8 | * | 8 | * |
| 9 | * This source code is the completely correct source code based on | 9 | * This source code is the completely correct source code based on |
| 10 | * the example codes provided by Atlantic Aerospace Division, Titan | 10 | * the example codes provided by Atlantic Aerospace Division, Titan |
| 11 | * Systems Corporation, 2000. | 11 | * Systems Corporation, 2000. |
| 12 | * | 12 | * |
| 13 | * If you just compile and generate the executables from this source | 13 | * If you just compile and generate the executables from this source |
| 14 | * code, this code would be enough. However, if you wish to get a complete | 14 | * code, this code would be enough. However, if you wish to get a complete |
| 15 | * understanding of this stressmark, it is strongly suggested that you | 15 | * understanding of this stressmark, it is strongly suggested that you |
| 16 | * read the Benchmark Analysis and Specifications Document Version 1.0 | 16 | * read the Benchmark Analysis and Specifications Document Version 1.0 |
| 17 | * before going on since the detailed comments are given in this documents. | 17 | * before going on since the detailed comments are given in this documents. |
| @@ -22,37 +22,37 @@ | |||
| 22 | * The Sparse Matrix Storage is implemented by Compact Row Storage Scheme | 22 | * The Sparse Matrix Storage is implemented by Compact Row Storage Scheme |
| 23 | * In the code, the data is first generated by randomNonzeroFloat() | 23 | * In the code, the data is first generated by randomNonzeroFloat() |
| 24 | * the data is first stored in a full-space matrix with size of dim*dim | 24 | * the data is first stored in a full-space matrix with size of dim*dim |
| 25 | * then the data is transfered to the Compact Row Matrix, | 25 | * then the data is transfered to the Compact Row Matrix, |
| 26 | * the data value is kept in *value, | 26 | * the data value is kept in *value, |
| 27 | * the columns corresponding to the value are stored in *col_ind, | 27 | * the columns corresponding to the value are stored in *col_ind, |
| 28 | * the start element of each row is stored in *row_start. | 28 | * the start element of each row is stored in *row_start. |
| 29 | */ | 29 | */ |
| 30 | 30 | ||
| 31 | /* | 31 | /* |
| 32 | * Please note: | 32 | * Please note: |
| 33 | * the total number of data is numberNonzero +dim | 33 | * the total number of data is numberNonzero +dim |
| 34 | * among which, NumberNonzero because this is symmetric matrix | 34 | * among which, NumberNonzero because this is symmetric matrix |
| 35 | * dim because the diagonal elements | 35 | * dim because the diagonal elements |
| 36 | */ | 36 | */ |
| 37 | 37 | ||
| 38 | #include <stdio.h> | 38 | #include "DISstressmarkRNG.h" |
| 39 | #include "extra.h" | ||
| 40 | #include <assert.h> | ||
| 39 | #include <math.h> | 41 | #include <math.h> |
| 42 | #include <stdio.h> | ||
| 40 | #include <stdlib.h> | 43 | #include <stdlib.h> |
| 41 | #include <time.h> | 44 | #include <time.h> |
| 42 | #include <assert.h> | ||
| 43 | #include "DISstressmarkRNG.h" | ||
| 44 | #include "extra.h" | ||
| 45 | 45 | ||
| 46 | #define MIN_SEED -2147483647 | 46 | #define MIN_SEED -2147483647 |
| 47 | #define MAX_SEED -1 | 47 | #define MAX_SEED -1 |
| 48 | #define MIN_DIM 1 | 48 | #define MIN_DIM 1 |
| 49 | #define MAX_DIM 32768 | 49 | #define MAX_DIM 32768 |
| 50 | #define MAX_ITERATIONS 65536 | 50 | #define MAX_ITERATIONS 65536 |
| 51 | #define MIN_TOLERANCE 1e-7//0.000007 | 51 | #define MIN_TOLERANCE 1e-7 // 0.000007 |
| 52 | #define MAX_TOLERANCE 0.5 | 52 | #define MAX_TOLERANCE 0.5 |
| 53 | #define MIN_NUMBER -3.4e10/dim | 53 | #define MIN_NUMBER -3.4e10 / dim |
| 54 | #define MAX_NUMBER 3.4e10/dim | 54 | #define MAX_NUMBER 3.4e10 / dim |
| 55 | #define EPSI 1.0e-10 | 55 | #define EPSI 1.0e-10 |
| 56 | #define MIN_DIG_NUMBER 1.0e-10 | 56 | #define MIN_DIG_NUMBER 1.0e-10 |
| 57 | #define MAX_DIG_NUMBER 3.4e10 | 57 | #define MAX_DIG_NUMBER 3.4e10 |
| 58 | 58 | ||
| @@ -62,97 +62,90 @@ | |||
| 62 | 62 | ||
| 63 | static int dim; | 63 | static int dim; |
| 64 | 64 | ||
| 65 | /* | 65 | /* |
| 66 | * matrix * vector | 66 | * matrix * vector |
| 67 | */ | 67 | */ |
| 68 | 68 | ||
| 69 | void matrixMulvector(double *value, | 69 | void matrixMulvector(double *value, int *col_ind, int *row_start, |
| 70 | int *col_ind, | 70 | double *vector, double *out) { |
| 71 | int *row_start, | ||
| 72 | double *vector, | ||
| 73 | double *out) | ||
| 74 | { | ||
| 75 | int l, ll; | 71 | int l, ll; |
| 76 | double sum; | 72 | double sum; |
| 77 | int tmp_rs, tmp_re; | 73 | int tmp_rs, tmp_re; |
| 78 | 74 | ||
| 79 | for (l=0; l<dim; l++){ | 75 | for (l = 0; l < dim; l++) { |
| 80 | *(out + l) = 0; | 76 | *(out + l) = 0; |
| 81 | tmp_rs = row_start[l]; | 77 | tmp_rs = row_start[l]; |
| 82 | 78 | ||
| 83 | if (tmp_rs != -1){ | 79 | if (tmp_rs != -1) { |
| 84 | tmp_re = row_start[l+1]; /* | 80 | tmp_re = row_start[l + 1]; /* |
| 85 | *get the start and ending elements of | 81 | *get the start and ending elements of |
| 86 | * each row | 82 | * each row |
| 87 | */ | 83 | */ |
| 88 | for (ll=tmp_rs; ll<tmp_re; ll++){ | 84 | for (ll = tmp_rs; ll < tmp_re; ll++) { |
| 89 | *(out + l) += value[ll]*vector[col_ind[ll]]; | 85 | *(out + l) += value[ll] * vector[col_ind[ll]]; |
| 90 | } | 86 | } |
| 91 | } | 87 | } |
| 92 | } | 88 | } |
| 93 | return; | 89 | return; |
| 94 | } | 90 | } |
| 95 | 91 | ||
| 96 | |||
| 97 | /* | 92 | /* |
| 98 | * vector1 - vector2 | 93 | * vector1 - vector2 |
| 99 | */ | 94 | */ |
| 100 | 95 | ||
| 101 | void vectorSub(double *vector1, double *vector2, double *vector){ | 96 | void vectorSub(double *vector1, double *vector2, double *vector) { |
| 102 | 97 | ||
| 103 | int l; | 98 | int l; |
| 104 | 99 | ||
| 105 | for (l=0; l<dim; l++){ | 100 | for (l = 0; l < dim; l++) { |
| 106 | *(vector + l) = *(vector1 + l) - *(vector2 + l); | 101 | *(vector + l) = *(vector1 + l) - *(vector2 + l); |
| 107 | } | 102 | } |
| 108 | return; | 103 | return; |
| 109 | } | 104 | } |
| 110 | 105 | ||
| 111 | |||
| 112 | /* | 106 | /* |
| 113 | * vector1 + vector2 | 107 | * vector1 + vector2 |
| 114 | */ | 108 | */ |
| 115 | 109 | ||
| 116 | void vectorAdd(double *vector1, double *vector2, double *vector){ | 110 | void vectorAdd(double *vector1, double *vector2, double *vector) { |
| 117 | 111 | ||
| 118 | int l; | 112 | int l; |
| 119 | 113 | ||
| 120 | for (l=0; l<dim; l++){ | 114 | for (l = 0; l < dim; l++) { |
| 121 | *(vector + l) = *(vector1 + l) + *(vector2 + l); | 115 | *(vector + l) = *(vector1 + l) + *(vector2 + l); |
| 122 | } | 116 | } |
| 123 | return; | 117 | return; |
| 124 | } | 118 | } |
| 125 | 119 | ||
| 126 | /* | 120 | /* |
| 127 | * vector1 * vector2 | 121 | * vector1 * vector2 |
| 128 | */ | 122 | */ |
| 129 | 123 | ||
| 130 | double vectorMul(double *vector1, double *vector2){ | 124 | double vectorMul(double *vector1, double *vector2) { |
| 131 | 125 | ||
| 132 | int l; | 126 | int l; |
| 133 | double product; | 127 | double product; |
| 134 | 128 | ||
| 135 | product = 0; | 129 | product = 0; |
| 136 | 130 | ||
| 137 | for (l=0; l<dim; l++){ | 131 | for (l = 0; l < dim; l++) { |
| 138 | product += (*(vector1 + l))*(*(vector2 + l)); | 132 | product += (*(vector1 + l)) * (*(vector2 + l)); |
| 139 | |||
| 140 | } | 133 | } |
| 141 | return product; | 134 | return product; |
| 142 | } | 135 | } |
| 143 | 136 | ||
| 144 | /* | 137 | /* |
| 145 | * /vector/ | 138 | * /vector/ |
| 146 | */ | 139 | */ |
| 147 | 140 | ||
| 148 | double vectorValue(double *vector){ | 141 | double vectorValue(double *vector) { |
| 149 | 142 | ||
| 150 | double value; | 143 | double value; |
| 151 | int l; | 144 | int l; |
| 152 | 145 | ||
| 153 | value = 0; | 146 | value = 0; |
| 154 | 147 | ||
| 155 | for (l=0; l<dim; l++){ | 148 | for (l = 0; l < dim; l++) { |
| 156 | value += (*(vector + l)) * (*(vector + l)); | 149 | value += (*(vector + l)) * (*(vector + l)); |
| 157 | } | 150 | } |
| 158 | 151 | ||
| @@ -164,94 +157,89 @@ double vectorValue(double *vector){ | |||
| 164 | * In fact, we return the original vector here | 157 | * In fact, we return the original vector here |
| 165 | */ | 158 | */ |
| 166 | 159 | ||
| 167 | void transpose(double *vector, double *vect){ | 160 | void transpose(double *vector, double *vect) { |
| 168 | 161 | ||
| 169 | int l; | 162 | int l; |
| 170 | 163 | ||
| 171 | for (l=0; l<dim; l++){ | 164 | for (l = 0; l < dim; l++) { |
| 172 | *(vect+l) = *(vector+l); | 165 | *(vect + l) = *(vector + l); |
| 173 | } | 166 | } |
| 174 | return; | 167 | return; |
| 175 | } | 168 | } |
| 176 | 169 | ||
| 177 | /* | 170 | /* |
| 178 | * value * <vector> | 171 | * value * <vector> |
| 179 | */ | 172 | */ |
| 180 | void valueMulvector(double value, double *vector, double *vect){ | 173 | void valueMulvector(double value, double *vector, double *vect) { |
| 181 | 174 | ||
| 182 | int l; | 175 | int l; |
| 183 | int lll, i; | 176 | int lll, i; |
| 184 | double tmp; | 177 | double tmp; |
| 185 | 178 | ||
| 186 | for (l=0; l<dim; l++){ | 179 | for (l = 0; l < dim; l++) { |
| 187 | *(vect + l) = *(vector + l) * value; | 180 | *(vect + l) = *(vector + l) * value; |
| 188 | } | 181 | } |
| 189 | return; | 182 | return; |
| 190 | } | 183 | } |
| 191 | 184 | ||
| 192 | /* | 185 | /* |
| 193 | * generate the data distributed sparsely in matrix | 186 | * generate the data distributed sparsely in matrix |
| 194 | */ | 187 | */ |
| 195 | 188 | ||
| 196 | void initMatrix(double *matrix, int dim, int numberNonzero){ | 189 | void initMatrix(double *matrix, int dim, int numberNonzero) { |
| 197 | 190 | ||
| 198 | int k, l, ll; | 191 | int k, l, ll; |
| 199 | int i, j; | 192 | int i, j; |
| 200 | 193 | ||
| 201 | int lll; | 194 | int lll; |
| 202 | double sum; | 195 | double sum; |
| 203 | 196 | ||
| 204 | for (k=0; k< dim*dim; k++){ | 197 | for (k = 0; k < dim * dim; k++) { |
| 205 | *(matrix + k) = 0; | 198 | *(matrix + k) = 0; |
| 206 | } | 199 | } |
| 207 | 200 | ||
| 208 | for (l=0; l<numberNonzero/2; l++){ | 201 | for (l = 0; l < numberNonzero / 2; l++) { |
| 209 | 202 | ||
| 210 | i = randomUInt(1, dim-1); | 203 | i = randomUInt(1, dim - 1); |
| 211 | j = randomUInt(0, i-1); | 204 | j = randomUInt(0, i - 1); |
| 212 | 205 | ||
| 213 | while (*(matrix + i*dim + j) != 0){ | 206 | while (*(matrix + i * dim + j) != 0) { |
| 214 | 207 | ||
| 215 | i++; | 208 | i++; |
| 216 | if (i == dim){ | 209 | if (i == dim) { |
| 217 | j++; | 210 | j++; |
| 218 | if (j == dim-1){ | 211 | if (j == dim - 1) { |
| 219 | j = 0; | 212 | j = 0; |
| 220 | i = 1; | 213 | i = 1; |
| 221 | } | 214 | } else { |
| 222 | else{ | 215 | i = j + 1; |
| 223 | i = j+1; | 216 | } |
| 224 | } | 217 | } |
| 225 | } | ||
| 226 | } | 218 | } |
| 227 | 219 | ||
| 228 | if (*(matrix + i*dim + j) == 0){ | 220 | if (*(matrix + i * dim + j) == 0) { |
| 229 | *(matrix + i*dim + j) = (double )randomNonZeroFloat(MIN_NUMBER, | 221 | *(matrix + i * dim + j) = |
| 230 | MAX_NUMBER, | 222 | (double)randomNonZeroFloat(MIN_NUMBER, MAX_NUMBER, EPSI); |
| 231 | EPSI); | 223 | *(matrix + j * dim + i) = *(matrix + i * dim + j); |
| 232 | *(matrix + j*dim + i) = *(matrix + i*dim + j); | ||
| 233 | } | 224 | } |
| 234 | } | 225 | } |
| 235 | |||
| 236 | for (ll=0; ll<dim; ll++){ | ||
| 237 | 226 | ||
| 227 | for (ll = 0; ll < dim; ll++) { | ||
| 228 | |||
| 229 | *(matrix + ll * dim + ll) = (double)randomNonZeroFloat( | ||
| 230 | -MAX_DIG_NUMBER, MAX_DIG_NUMBER, MIN_DIG_NUMBER); | ||
| 238 | 231 | ||
| 239 | |||
| 240 | *(matrix + ll*dim + ll) = (double )randomNonZeroFloat(-MAX_DIG_NUMBER, | ||
| 241 | MAX_DIG_NUMBER, | ||
| 242 | MIN_DIG_NUMBER); | ||
| 243 | |||
| 244 | sum = 0; | 232 | sum = 0; |
| 245 | 233 | ||
| 246 | for (lll=0; lll<dim; lll++){ | 234 | for (lll = 0; lll < dim; lll++) { |
| 247 | if (lll != ll){ | 235 | if (lll != ll) { |
| 248 | sum += *(matrix + lll*dim + ll); | 236 | sum += *(matrix + lll * dim + ll); |
| 249 | } | 237 | } |
| 250 | } | 238 | } |
| 251 | 239 | ||
| 252 | if (*(matrix + ll*dim + ll) < sum ){ | 240 | if (*(matrix + ll * dim + ll) < sum) { |
| 253 | *(matrix + ll*dim + ll) += sum; | 241 | *(matrix + ll * dim + ll) += sum; |
| 254 | } | 242 | } |
| 255 | } | 243 | } |
| 256 | 244 | ||
| 257 | return; | 245 | return; |
| @@ -260,13 +248,13 @@ void initMatrix(double *matrix, int dim, int numberNonzero){ | |||
| 260 | /* | 248 | /* |
| 261 | * generate the data value in the vectors | 249 | * generate the data value in the vectors |
| 262 | */ | 250 | */ |
| 263 | 251 | ||
| 264 | void initVector(double *vector, int dim){ | 252 | void initVector(double *vector, int dim) { |
| 265 | 253 | ||
| 266 | int l; | 254 | int l; |
| 267 | 255 | ||
| 268 | for (l=0; l<dim; l++){ | 256 | for (l = 0; l < dim; l++) { |
| 269 | *(vector + l) = (double )randomFloat (MIN_NUMBER, MAX_NUMBER); | 257 | *(vector + l) = (double)randomFloat(MIN_NUMBER, MAX_NUMBER); |
| 270 | } | 258 | } |
| 271 | 259 | ||
| 272 | return; | 260 | return; |
| @@ -276,10 +264,10 @@ void initVector(double *vector, int dim){ | |||
| 276 | * make a vector contains value of zero | 264 | * make a vector contains value of zero |
| 277 | */ | 265 | */ |
| 278 | 266 | ||
| 279 | void zeroVector(double *vector, int dim){ | 267 | void zeroVector(double *vector, int dim) { |
| 280 | int l; | 268 | int l; |
| 281 | 269 | ||
| 282 | for (l=0; l<dim; l++){ | 270 | for (l = 0; l < dim; l++) { |
| 283 | *(vector + l) = 0; | 271 | *(vector + l) = 0; |
| 284 | } | 272 | } |
| 285 | return; | 273 | return; |
| @@ -289,48 +277,37 @@ void zeroVector(double *vector, int dim){ | |||
| 289 | * return a vector which is the copy of the vect | 277 | * return a vector which is the copy of the vect |
| 290 | */ | 278 | */ |
| 291 | 279 | ||
| 292 | void equalVector(double *vect, double *vect1){ | 280 | void equalVector(double *vect, double *vect1) { |
| 293 | 281 | ||
| 294 | int l; | 282 | int l; |
| 295 | 283 | ||
| 296 | for (l=0; l<dim; l++){ | 284 | for (l = 0; l < dim; l++) { |
| 297 | *(vect1+l) = *(vect+l); | 285 | *(vect1 + l) = *(vect + l); |
| 298 | } | 286 | } |
| 299 | return; | 287 | return; |
| 300 | } | 288 | } |
| 301 | 289 | ||
| 302 | 290 | void biConjugateGradient(double *value, int *col_ind, int *row_start, | |
| 303 | 291 | double *vectorB, double *vectorX, | |
| 304 | void biConjugateGradient(double *value, | 292 | double errorTolerance, int maxIterations, |
| 305 | int *col_ind, | 293 | double *actualError, int *actualIteration, int dim, |
| 306 | int *row_start, | 294 | // Start internal matrixes |
| 307 | double *vectorB, | 295 | double *vectorP, double *vectorR, double *nextVectorR, |
| 308 | double *vectorX, | 296 | double *tmpVector1, double *tmpVector2, |
| 309 | double errorTolerance, | 297 | double *tmpVector3) |
| 310 | int maxIterations, | 298 | /* |
| 311 | double *actualError, | 299 | * in the code, we use a lot of temparary vectors and variables |
| 312 | int *actualIteration, | 300 | * this is just for simple and clear |
| 313 | int dim, | 301 | * you can optimize these temporary variables and vectors |
| 314 | // Start internal matrixes | 302 | * based on your need |
| 315 | double *vectorP, | 303 | * |
| 316 | double *vectorR, | 304 | */ |
| 317 | double *nextVectorR, | ||
| 318 | double *tmpVector1, | ||
| 319 | double *tmpVector2, | ||
| 320 | double *tmpVector3) | ||
| 321 | /* | ||
| 322 | * in the code, we use a lot of temparary vectors and variables | ||
| 323 | * this is just for simple and clear | ||
| 324 | * you can optimize these temporary variables and vectors | ||
| 325 | * based on your need | ||
| 326 | * | ||
| 327 | */ | ||
| 328 | { | 305 | { |
| 329 | double error; | 306 | double error; |
| 330 | int iteration; | 307 | int iteration; |
| 331 | double alpha, beta; | 308 | double alpha, beta; |
| 332 | 309 | ||
| 333 | double tmpValue1, tmpValue2; | 310 | double tmpValue1, tmpValue2; |
| 334 | int i; | 311 | int i; |
| 335 | int l; | 312 | int l; |
| 336 | int ll; | 313 | int ll; |
| @@ -341,7 +318,7 @@ void biConjugateGradient(double *value, | |||
| 341 | /* | 318 | /* |
| 342 | * vectorR = vectorB - matrixA*vectorX | 319 | * vectorR = vectorB - matrixA*vectorX |
| 343 | */ | 320 | */ |
| 344 | matrixMulvector(value,col_ind, row_start, vectorX, tmpVector1); | 321 | matrixMulvector(value, col_ind, row_start, vectorX, tmpVector1); |
| 345 | 322 | ||
| 346 | vectorSub(vectorB, tmpVector1, vectorR); | 323 | vectorSub(vectorB, tmpVector1, vectorR); |
| 347 | 324 | ||
| @@ -354,35 +331,35 @@ void biConjugateGradient(double *value, | |||
| 354 | /* | 331 | /* |
| 355 | * error = |matrixA * vectorX - vectorB| / |vectorB| | 332 | * error = |matrixA * vectorX - vectorB| / |vectorB| |
| 356 | */ | 333 | */ |
| 357 | vectorSub(tmpVector1, vectorB, tmpVector1); | 334 | vectorSub(tmpVector1, vectorB, tmpVector1); |
| 358 | 335 | ||
| 359 | error = vectorValue(tmpVector1)/vectorValue(vectorB); | 336 | error = vectorValue(tmpVector1) / vectorValue(vectorB); |
| 360 | 337 | ||
| 361 | iteration = 0; | 338 | iteration = 0; |
| 362 | 339 | ||
| 363 | while ((iteration < maxIterations) && (error > errorTolerance)){ | 340 | while ((iteration < maxIterations) && (error > errorTolerance)) { |
| 364 | 341 | ||
| 365 | /* | 342 | /* |
| 366 | * alpha = (transpose(vectorR) * vectorR) / | 343 | * alpha = (transpose(vectorR) * vectorR) / |
| 367 | * (transpose(vectorP) * (matrixA * vectorP) | 344 | * (transpose(vectorP) * (matrixA * vectorP) |
| 368 | */ | 345 | */ |
| 369 | 346 | ||
| 370 | matrixMulvector(value, col_ind, row_start, vectorP, tmpVector1); | 347 | matrixMulvector(value, col_ind, row_start, vectorP, tmpVector1); |
| 371 | transpose(vectorR, tmpVector2); | 348 | transpose(vectorR, tmpVector2); |
| 372 | transpose(vectorP, tmpVector3); | 349 | transpose(vectorP, tmpVector3); |
| 373 | tmpValue1 = vectorMul(tmpVector3, tmpVector1); | 350 | tmpValue1 = vectorMul(tmpVector3, tmpVector1); |
| 374 | tmpValue2 = vectorMul(tmpVector2, vectorR); | 351 | tmpValue2 = vectorMul(tmpVector2, vectorR); |
| 375 | alpha = tmpValue2/tmpValue1; | 352 | alpha = tmpValue2 / tmpValue1; |
| 376 | 353 | ||
| 377 | /* | 354 | /* |
| 378 | * nextVectorR = vectorR - alpha*(matrixA * vectorP) | 355 | * nextVectorR = vectorR - alpha*(matrixA * vectorP) |
| 379 | */ | 356 | */ |
| 380 | 357 | ||
| 381 | valueMulvector(alpha, tmpVector1, tmpVector2); | 358 | valueMulvector(alpha, tmpVector1, tmpVector2); |
| 382 | vectorSub(vectorR, tmpVector2, tmpVector1); | 359 | vectorSub(vectorR, tmpVector2, tmpVector1); |
| 383 | equalVector(tmpVector1, nextVectorR); | 360 | equalVector(tmpVector1, nextVectorR); |
| 384 | 361 | ||
| 385 | /* | 362 | /* |
| 386 | * beta = (transpose(nextVectorR) * nextVectorR) / | 363 | * beta = (transpose(nextVectorR) * nextVectorR) / |
| 387 | * (transpose(vectorR) * vectorR) | 364 | * (transpose(vectorR) * vectorR) |
| 388 | */ | 365 | */ |
| @@ -391,38 +368,38 @@ void biConjugateGradient(double *value, | |||
| 391 | tmpValue1 = vectorMul(tmpVector3, nextVectorR); | 368 | tmpValue1 = vectorMul(tmpVector3, nextVectorR); |
| 392 | transpose(vectorR, tmpVector2); | 369 | transpose(vectorR, tmpVector2); |
| 393 | tmpValue2 = vectorMul(tmpVector2, vectorR); | 370 | tmpValue2 = vectorMul(tmpVector2, vectorR); |
| 394 | beta = tmpValue1/tmpValue2; | 371 | beta = tmpValue1 / tmpValue2; |
| 395 | 372 | ||
| 396 | /* | 373 | /* |
| 397 | * vectorX = vectorX + alpha * vectorP | 374 | * vectorX = vectorX + alpha * vectorP |
| 398 | */ | 375 | */ |
| 399 | valueMulvector(alpha, vectorP, tmpVector1); | 376 | valueMulvector(alpha, vectorP, tmpVector1); |
| 400 | vectorAdd(vectorX,tmpVector1, vectorX); | 377 | vectorAdd(vectorX, tmpVector1, vectorX); |
| 401 | 378 | ||
| 402 | /* | 379 | /* |
| 403 | *vectorP = nextVectorR + beta*vectorP | 380 | *vectorP = nextVectorR + beta*vectorP |
| 404 | */ | 381 | */ |
| 405 | valueMulvector(beta, vectorP, tmpVector1); | 382 | valueMulvector(beta, vectorP, tmpVector1); |
| 406 | vectorAdd(nextVectorR, tmpVector1, tmpVector1); | 383 | vectorAdd(nextVectorR, tmpVector1, tmpVector1); |
| 407 | 384 | ||
| 408 | for (ll=0; ll<dim; ll++){ | 385 | for (ll = 0; ll < dim; ll++) { |
| 409 | *(vectorP + ll) = *(tmpVector1 + ll); | 386 | *(vectorP + ll) = *(tmpVector1 + ll); |
| 410 | } | 387 | } |
| 411 | 388 | ||
| 412 | /* | 389 | /* |
| 413 | * vectorR = nextVectorR | 390 | * vectorR = nextVectorR |
| 414 | */ | 391 | */ |
| 415 | 392 | ||
| 416 | for (l=0; l<dim; l++){ | 393 | for (l = 0; l < dim; l++) { |
| 417 | *(vectorR+l) = *(nextVectorR+l); | 394 | *(vectorR + l) = *(nextVectorR + l); |
| 418 | } | 395 | } |
| 419 | 396 | ||
| 420 | /* | 397 | /* |
| 421 | * error = |matrixA * vectorX - vectorB| / |vectorB| | 398 | * error = |matrixA * vectorX - vectorB| / |vectorB| |
| 422 | */ | 399 | */ |
| 423 | matrixMulvector(value, col_ind,row_start, vectorX, tmpVector1); | 400 | matrixMulvector(value, col_ind, row_start, vectorX, tmpVector1); |
| 424 | vectorSub(tmpVector1,vectorB,tmpVector1); | 401 | vectorSub(tmpVector1, vectorB, tmpVector1); |
| 425 | error = vectorValue(tmpVector1)/vectorValue(vectorB); | 402 | error = vectorValue(tmpVector1) / vectorValue(vectorB); |
| 426 | 403 | ||
| 427 | iteration++; | 404 | iteration++; |
| 428 | } | 405 | } |
| @@ -439,63 +416,57 @@ void biConjugateGradient(double *value, | |||
| 439 | 416 | ||
| 440 | return; | 417 | return; |
| 441 | } | 418 | } |
| 442 | 419 | ||
| 443 | /* | 420 | /* |
| 444 | * This is the function to transfer the data from the matrix of dense storage | 421 | * This is the function to transfer the data from the matrix of dense storage |
| 445 | * to Compact Row Storage | 422 | * to Compact Row Storage |
| 446 | */ | 423 | */ |
| 447 | void create_CRS(double *matrixA, | 424 | void create_CRS(double *matrixA, double *value, int *col_ind, int *row_start, |
| 448 | double *value, | 425 | int dim, int numberNonzero) { |
| 449 | int *col_ind, | ||
| 450 | int *row_start, | ||
| 451 | int dim, | ||
| 452 | int numberNonzero) | ||
| 453 | { | ||
| 454 | 426 | ||
| 455 | int i, j, k; | 427 | int i, j, k; |
| 456 | int cnt; | 428 | int cnt; |
| 457 | double tmp; | 429 | double tmp; |
| 458 | 430 | ||
| 459 | /* | 431 | /* |
| 460 | *initialize the row_start | 432 | *initialize the row_start |
| 461 | */ | 433 | */ |
| 462 | 434 | ||
| 463 | for(k=0; k<dim; k++){ | 435 | for (k = 0; k < dim; k++) { |
| 464 | row_start[k] = -1; | 436 | row_start[k] = -1; |
| 465 | } | 437 | } |
| 466 | 438 | ||
| 467 | /* | 439 | /* |
| 468 | * make the end of the last row to be numberNonzero + dim. | 440 | * make the end of the last row to be numberNonzero + dim. |
| 469 | */ | 441 | */ |
| 470 | 442 | ||
| 471 | row_start[dim] = numberNonzero+dim; | 443 | row_start[dim] = numberNonzero + dim; |
| 472 | 444 | ||
| 473 | /* | 445 | /* |
| 474 | * initialize the col_ind | 446 | * initialize the col_ind |
| 475 | */ | 447 | */ |
| 476 | 448 | ||
| 477 | for (k=0; k<numberNonzero+dim; k++){ | 449 | for (k = 0; k < numberNonzero + dim; k++) { |
| 478 | col_ind[k] = -1; | 450 | col_ind[k] = -1; |
| 479 | } | 451 | } |
| 480 | 452 | ||
| 481 | |||
| 482 | cnt = 0; | 453 | cnt = 0; |
| 483 | 454 | ||
| 484 | for (i=0; (cnt<numberNonzero+dim)&&(i<dim); i++){ | 455 | for (i = 0; (cnt < numberNonzero + dim) && (i < dim); i++) { |
| 485 | for (j=0; (cnt<numberNonzero+dim)&&(j<dim); j++){ | 456 | for (j = 0; (cnt < numberNonzero + dim) && (j < dim); j++) { |
| 486 | 457 | ||
| 487 | tmp = *(matrixA + i*dim + j); | 458 | tmp = *(matrixA + i * dim + j); |
| 488 | 459 | ||
| 489 | if (tmp!=0){ | 460 | if (tmp != 0) { |
| 490 | 461 | ||
| 491 | value[cnt] = tmp; | 462 | value[cnt] = tmp; |
| 492 | col_ind[cnt] = j; | 463 | col_ind[cnt] = j; |
| 493 | 464 | ||
| 494 | if (row_start[i] == -1) | 465 | if (row_start[i] == -1) |
| 495 | row_start[i] = cnt; | 466 | row_start[i] = cnt; |
| 496 | 467 | ||
| 497 | cnt += 1; | 468 | cnt += 1; |
| 498 | } | 469 | } |
| 499 | } | 470 | } |
| 500 | } | 471 | } |
| 501 | row_start[i] = cnt; | 472 | row_start[i] = cnt; |
| @@ -503,56 +474,54 @@ void create_CRS(double *matrixA, | |||
| 503 | return; | 474 | return; |
| 504 | } | 475 | } |
| 505 | 476 | ||
| 506 | 477 | int main(int argc, char **argv) { | |
| 507 | int main(int argc, char** argv) | ||
| 508 | { | ||
| 509 | int seed; | 478 | int seed; |
| 510 | int numberNonzero; | 479 | int numberNonzero; |
| 511 | int maxIterations; | 480 | int maxIterations; |
| 512 | float errorTolerance; | 481 | float errorTolerance; |
| 513 | double actualError; | 482 | double actualError; |
| 514 | int actualIteration; | 483 | int actualIteration; |
| 515 | 484 | ||
| 516 | time_t beginTime; | 485 | time_t beginTime; |
| 517 | time_t endTime; | 486 | time_t endTime; |
| 518 | 487 | ||
| 519 | double *matrixA; | 488 | double *matrixA; |
| 520 | double *vectorB; | 489 | double *vectorB; |
| 521 | double *vectorX; | 490 | double *vectorX; |
| 522 | 491 | ||
| 523 | double *value; | 492 | double *value; |
| 524 | int *col_ind; | 493 | int *col_ind; |
| 525 | int *row_start; | 494 | int *row_start; |
| 526 | double sum; | 495 | double sum; |
| 527 | int k; | 496 | int k; |
| 528 | // Internal vects for biConj | 497 | // Internal vects for biConj |
| 529 | double *vectorR, *vectorP, *nextVectorR; | 498 | double *vectorR, *vectorP, *nextVectorR; |
| 530 | double *tmpVector1, *tmpVector2, *tmpVector3; | 499 | double *tmpVector1, *tmpVector2, *tmpVector3; |
| 531 | SET_UP | 500 | SET_UP |
| 532 | 501 | ||
| 533 | assert(fscanf(stdin, "%d %d %d %d %f", | 502 | assert(fscanf(stdin, "%d %d %d %d %f", &seed, &dim, &numberNonzero, |
| 534 | &seed, &dim, &numberNonzero,&maxIterations,&errorTolerance) == 5); | 503 | &maxIterations, &errorTolerance) == 5); |
| 535 | assert((seed > MIN_SEED) && (seed < MAX_SEED)); | 504 | assert((seed > MIN_SEED) && (seed < MAX_SEED)); |
| 536 | assert((dim > MIN_DIM) && (dim < MAX_DIM)); | 505 | assert((dim > MIN_DIM) && (dim < MAX_DIM)); |
| 537 | assert((numberNonzero > dim) && (numberNonzero < dim*dim)); | 506 | assert((numberNonzero > dim) && (numberNonzero < dim * dim)); |
| 538 | assert((maxIterations > 0) && (maxIterations < MAX_ITERATIONS)); | 507 | assert((maxIterations > 0) && (maxIterations < MAX_ITERATIONS)); |
| 539 | assert((errorTolerance > MIN_TOLERANCE) && (errorTolerance < MAX_TOLERANCE)); | 508 | assert((errorTolerance > MIN_TOLERANCE) && (errorTolerance < MAX_TOLERANCE)); |
| 540 | |||
| 541 | matrixA = (double *)malloc(dim*dim*sizeof(double)); | ||
| 542 | vectorB = (double *)malloc(dim*sizeof(double)); | ||
| 543 | vectorX = (double *)malloc(dim*sizeof(double)); | ||
| 544 | 509 | ||
| 545 | value = (double *)malloc((numberNonzero+dim)*sizeof(double)); | 510 | matrixA = (double *)malloc(dim * dim * sizeof(double)); |
| 546 | col_ind = (int *)malloc((numberNonzero+dim)*sizeof(int)); | 511 | vectorB = (double *)malloc(dim * sizeof(double)); |
| 547 | row_start = (int *)malloc((dim+1)*sizeof(int)); | 512 | vectorX = (double *)malloc(dim * sizeof(double)); |
| 513 | |||
| 514 | value = (double *)malloc((numberNonzero + dim) * sizeof(double)); | ||
| 515 | col_ind = (int *)malloc((numberNonzero + dim) * sizeof(int)); | ||
| 516 | row_start = (int *)malloc((dim + 1) * sizeof(int)); | ||
| 548 | 517 | ||
| 549 | // Internal matricies for biConj | 518 | // Internal matricies for biConj |
| 550 | vectorP = (double *)malloc(dim*sizeof(double)); | 519 | vectorP = (double *)malloc(dim * sizeof(double)); |
| 551 | vectorR = (double *)malloc(dim*sizeof(double)); | 520 | vectorR = (double *)malloc(dim * sizeof(double)); |
| 552 | nextVectorR = (double *)malloc(dim*sizeof(double)); | 521 | nextVectorR = (double *)malloc(dim * sizeof(double)); |
| 553 | tmpVector1 = (double *)malloc(dim*sizeof(double)); | 522 | tmpVector1 = (double *)malloc(dim * sizeof(double)); |
| 554 | tmpVector2 = (double *)malloc(dim*sizeof(double)); | 523 | tmpVector2 = (double *)malloc(dim * sizeof(double)); |
| 555 | tmpVector3 = (double *)malloc(dim*sizeof(double)); | 524 | tmpVector3 = (double *)malloc(dim * sizeof(double)); |
| 556 | 525 | ||
| 557 | randInit(seed); | 526 | randInit(seed); |
| 558 | 527 | ||
| @@ -569,23 +538,24 @@ int main(int argc, char** argv) | |||
| 569 | actualError = 0; | 538 | actualError = 0; |
| 570 | actualIteration = 0; | 539 | actualIteration = 0; |
| 571 | 540 | ||
| 572 | biConjugateGradient(value, col_ind, row_start, vectorB, vectorX, errorTolerance, | 541 | biConjugateGradient(value, col_ind, row_start, vectorB, vectorX, |
| 573 | maxIterations, | 542 | errorTolerance, maxIterations, &actualError, |
| 574 | &actualError, &actualIteration, dim, | 543 | &actualIteration, dim, vectorP, vectorR, nextVectorR, |
| 575 | vectorP, vectorR, nextVectorR, tmpVector1, tmpVector2, tmpVector3); | 544 | tmpVector1, tmpVector2, tmpVector3); |
| 576 | 545 | ||
| 577 | STOP_LOOP | 546 | STOP_LOOP |
| 578 | endTime = time(NULL); | 547 | endTime = time(NULL); |
| 579 | 548 | ||
| 580 | sum = 0; | 549 | sum = 0; |
| 581 | for (k=1; k<dim; k++){ | 550 | for (k = 1; k < dim; k++) { |
| 582 | sum += sum + *(vectorX + k); | 551 | sum += sum + *(vectorX + k); |
| 583 | } | 552 | } |
| 584 | 553 | ||
| 585 | fprintf(stdout, "sum = %f, actualError = %e, actualIteration = %d\n", sum, actualError, actualIteration); | 554 | fprintf(stdout, "sum = %f, actualError = %e, actualIteration = %d\n", sum, |
| 586 | fprintf(stderr, "time for matrix stressmark = %f seconds.\n", difftime(endTime, beginTime)); | 555 | actualError, actualIteration); |
| 556 | fprintf(stderr, "time for matrix stressmark = %f seconds.\n", | ||
| 557 | difftime(endTime, beginTime)); | ||
| 587 | 558 | ||
| 588 | WRITE_TO_FILE | 559 | WRITE_TO_FILE |
| 589 | return(0); | 560 | return (0); |
| 590 | } | 561 | } |
| 591 | |||
