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