1
- //
2
- // Created by pkarakal on 12/6/20.
3
- //
1
+ #include < iostream>
2
+ #include < cstdlib>
3
+ #include < chrono>
4
+ #include < vector>
5
+ #include < pthread.h>
6
+ #include < algorithm>
7
+ #include " Graph.hpp"
4
8
9
+ extern " C" {
10
+ #include " mmio.h"
11
+ #include " coo2csc.h"
12
+ }
13
+
14
+
15
+ int main (int argc, char ** argv){
16
+ int ret_code;
17
+ MM_typecode matcode;
18
+ FILE *f;
19
+ int M, N, nnz;
20
+ uint32_t i;
21
+ std::vector<uint32_t > I, J;
22
+ std::vector<double > val;
23
+
24
+ if (argc < 2 )
25
+ {
26
+ fprintf (stderr, " Usage: %s [martix-market-filename] [0 for binary or 1 for non binary]\n " , argv[0 ]);
27
+ exit (1 );
28
+ }
29
+ else
30
+ {
31
+ if ((f = fopen (argv[1 ], " r" )) == nullptr )
32
+ exit (1 );
33
+ }
34
+
35
+ if (mm_read_banner (f, &matcode) != 0 )
36
+ {
37
+ printf (" Could not process Matrix Market banner.\n " );
38
+ exit (1 );
39
+ }
40
+
41
+ int threads{};
42
+ if (argc == 3 ){
43
+ threads= atoi (argv[2 ]);
44
+ } else {
45
+ threads = 4 ;
46
+ }
47
+
48
+
49
+ /* This is how one can screen matrix types if their application */
50
+ /* only supports a subset of the Matrix Market data types. */
51
+
52
+ if (mm_is_complex (matcode) && mm_is_matrix (matcode) &&
53
+ mm_is_sparse (matcode) )
54
+ {
55
+ printf (" Sorry, this application does not support " );
56
+ printf (" Market Market type: [%s]\n " , mm_typecode_to_str (matcode));
57
+ exit (1 );
58
+ }
59
+
60
+ /* find out size of sparse matrix .... */
61
+
62
+ if ((ret_code = mm_read_mtx_crd_size (f, &M, &N, &nnz)) != 0 )
63
+ exit (1 );
64
+
65
+
66
+ I = std::vector<uint32_t >(nnz);
67
+ J = std::vector<uint32_t >(nnz);
68
+ val = std::vector<double >(nnz);
69
+
70
+ std::vector<uint32_t > cscRow = std::vector<uint32_t >(2 *nnz);
71
+ std::vector<uint32_t > cscColumn = std::vector<uint32_t >(N+1 );
72
+ std::vector<uint32_t > c_values = std::vector<uint32_t >(0 );
73
+
74
+ /* NOTE: when reading in doubles, ANSI C requires the use of the "l" */
75
+ /* specifier as in "%lg", "%lf", "%le", otherwise errors will occur */
76
+ /* (ANSI C X3.159-1989, Sec. 4.9.6.2, p. 136 lines 13-15) */
77
+
78
+ /* Replace missing val column with 1s and change the fscanf to match pattern matrices*/
79
+
80
+ if (!mm_is_pattern (matcode)) {
81
+ for (i = 0 ; i < nnz; i++) {
82
+ fscanf (f, " %d %d %lg\n " , &I[i], &J[i], &val[i]);
83
+ I[i]--; /* adjust from 1-based to 0-based */
84
+ J[i]--;
85
+ }
86
+ } else {
87
+ for (i = 0 ; i < nnz; i++) {
88
+ fscanf (f, " %d %d\n " , &I[i], &J[i]);
89
+ val[i] = 1 ;
90
+ I[i]--; /* adjust from 1-based to 0-based */
91
+ J[i]--;
92
+ }
93
+ }
94
+
95
+ if (f !=stdin) fclose (f);
96
+
97
+ if (M != N) {
98
+ printf (" COO matrix' columns and rows are not the same" );
99
+ }
100
+
101
+ // create symmetric values
102
+ std::vector<uint32_t > temp1 = std::vector<uint32_t >(I.begin (), I.end ());
103
+ I.insert (std::end (I), std::begin (J), std::end (J));
104
+ J.insert (std::end (J), std::begin (temp1), std::end (temp1));
105
+ temp1.clear ();
106
+
107
+ if (I[0 ] < J[0 ]) {
108
+ coo2csc (&cscRow[0 ], &cscColumn[0 ], &I[0 ], &J[0 ], 2 * nnz, M, 0 );
109
+ } else {
110
+ coo2csc (&cscRow[0 ], &cscColumn[0 ], &J[0 ], &I[0 ], 2 * nnz, N, 0 );
111
+ }
112
+
113
+ std::sort (cscColumn.begin (), cscColumn.end ());
114
+
115
+ std::vector<int >c3 = std::vector<int >(0 );
116
+ std::vector<int >ones = std::vector<int >(N, 1 );
117
+ std::vector<int >result_vector = std::vector<int >(N, 0 );
118
+ c_values = std::vector<uint32_t >(2 *nnz);
119
+
120
+ auto start = std::chrono::high_resolution_clock::now ();
121
+
122
+ Graph graphs[threads];
123
+
124
+ pthread_t *pthreads;
125
+ pthreads = (pthread_t *)malloc (threads*sizeof (pthread_t ));
126
+
127
+ int chunk = 1 ;
128
+ if (threads > 0 ) {
129
+ chunk = N / (threads);
130
+ }
131
+
132
+ for (i = 0 ; i < threads-1 ; i++) {
133
+ graphs[i] = Graph (cscRow, cscColumn, c_values, i*chunk, (i+1 )*chunk, nnz, i);
134
+ pthread_create (&pthreads[i], NULL , Graph::statAdjMatMul, &graphs[i]);
135
+ }
136
+ // The last thread is left out so as to calculate the mod of the chunk division!
137
+
138
+ graphs[threads -1 ] = Graph (cscRow,cscColumn, c_values, (threads-1 )*chunk, (threads)*chunk + (N%threads), nnz, threads -1 );
139
+
140
+ pthread_create (&pthreads[threads - 1 ], NULL , Graph::statAdjMatMul, &graphs[threads - 1 ]);
141
+
142
+ for (i = 0 ; i < threads; i++) {
143
+ pthread_join (pthreads[i], NULL );
144
+ }
145
+
146
+ for (i = 0 ; i < N; i++) {
147
+ for (int j = 0 ; j < cscColumn.at (i+1 ) - cscColumn.at (i); j++) {
148
+ int row = cscRow.at (cscColumn.at (i) + j);
149
+ int col = i;
150
+ int value = c_values.at (cscColumn.at (i) + j);
151
+ result_vector.at (row) += value * ones.at (col);
152
+ }
153
+ }
154
+
155
+ for (int res: result_vector){
156
+ c3.push_back (res/2 );
157
+ }
158
+
159
+ auto stop = std::chrono::high_resolution_clock::now ();
160
+
161
+ std::chrono::duration<double > elapsed = stop - start;
162
+ std::cout<<" Took " << elapsed.count () <<std::endl;
163
+
164
+ // for(int item: c3){
165
+ // std::cout<< item << " ";
166
+ // }
167
+ std::cout<<std::endl;
168
+
169
+ return 0 ;
170
+ }
0 commit comments