summaryrefslogtreecommitdiffstats
path: root/dis/Matrix/ver2/matrix.c
diff options
context:
space:
mode:
authorJoshua Bakita <jbakita@cs.unc.edu>2020-10-16 16:55:14 -0400
committerJoshua Bakita <jbakita@cs.unc.edu>2020-10-16 16:55:14 -0400
commit6ea9939e0610a809f6f47d13ec68df00d1ca0afc (patch)
treefe4a2eee3ddcf77e2367309dcd75a232b76dcd62 /dis/Matrix/ver2/matrix.c
parente9285d0cdea756a2830f0ace378e4197b36869aa (diff)
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.
Diffstat (limited to 'dis/Matrix/ver2/matrix.c')
-rwxr-xr-xdis/Matrix/ver2/matrix.c594
1 files changed, 594 insertions, 0 deletions
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