29
29
30
30
#include " SquareMatrixMultiply.cpp"
31
31
32
+ template <typename T, typename F>
33
+ void parallel_for (T i1, T i2, F func) {
34
+ if (i1 == i2 - 1 )
35
+ func (i1);
36
+ else {
37
+ T mid = (i1 + i2) / 2 ;
38
+ std::thread t1 (parallel_for<T, F>, i1, mid, func);
39
+ parallel_for<T, F>(mid, i2, func);
40
+ t1.join ();
41
+ }
42
+ }
43
+
32
44
template <typename T>
33
45
void MatVecMainLoop (Matrix<T>* A, Matrix<T>* x, Matrix<T>* y, size_t n,
34
46
size_t i1, size_t i2) {
@@ -44,63 +56,60 @@ void MatVecMainLoop(Matrix<T>* A, Matrix<T>* x, Matrix<T>* y, size_t n,
44
56
}
45
57
46
58
template <typename T>
47
- void MatVec (Matrix<T>& A, Matrix<T>& x, Matrix<T>& y) {
59
+ void MatVecRecursive (Matrix<T>& A, Matrix<T>& x, Matrix<T>& y) {
48
60
size_t n = A.rows ;
49
61
MatVecMainLoop (&A, &x, &y, n, 0 , n);
50
62
}
51
63
52
64
template <typename T>
53
- void MatVecWrongLoop2 (Matrix<T>* A, Matrix<T>* x, Matrix<T>* y, size_t n,
54
- size_t i, size_t j1, size_t j2) {
55
- if (j1 == j2 - 1 )
56
- (*y)[i][0 ] += (*A)[i][j1] * (*x)[j1][0 ];
57
- else {
58
- size_t mid = (j1 + j2) / 2 ;
59
- std::thread t1 (MatVecWrongLoop2<T>, A, x, y, n, i, j1, mid);
60
- MatVecWrongLoop2 (A, x, y, n, i, mid, j2);
61
- t1.join ();
62
- }
63
- }
64
-
65
- template <typename T>
66
- void MatVecWrongLoop1 (Matrix<T>* A, Matrix<T>* x, Matrix<T>* y, size_t n,
67
- size_t i1, size_t i2) {
68
- if (i1 == i2 - 1 )
69
- MatVecWrongLoop2 (A, x, y, n, i1, 0 , n);
70
- else {
71
- size_t mid = (i1 + i2) / 2 ;
72
- std::thread t1 (MatVecWrongLoop1<T>, A, x, y, n, i1, mid);
73
- MatVecWrongLoop1 (A, x, y, n, mid, i2);
74
- t1.join ();
75
- }
65
+ void MatVec (Matrix<T>& A, Matrix<T>& x, Matrix<T>& y) {
66
+ size_t n = A.rows ;
67
+ parallel_for<size_t >(0 , n, [&](size_t i){
68
+ for (size_t j = 0 ; j < n; j++)
69
+ y[i][0 ] += A[i][j] * x[j][0 ];
70
+ });
76
71
}
77
72
78
73
template <typename T>
79
74
void MatVecWrong (Matrix<T>& A, Matrix<T>& x, Matrix<T>& y) {
80
75
size_t n = A.rows ;
81
- MatVecWrongLoop1 (&A, &x, &y, n, 0 , n);
76
+ parallel_for<size_t >(0 , n, [&](size_t i){
77
+ parallel_for<size_t >(0 , n, [&](size_t j){
78
+ y[i][0 ] += A[i][j] * x[j][0 ];
79
+ });
80
+ });
82
81
}
83
82
#endif
84
83
85
84
#ifdef MAIN_MatVec
86
85
int main (int argc, char *argv[]) {
87
86
size_t n = get_argv (argc, argv, 1 , 10 );
88
- size_t compute = get_argv (argc, argv, 2 , 1 );
87
+ size_t compute = get_argv (argc, argv, 2 , 7 );
89
88
std::vector<int > buf_A, buf_x;
90
89
random_integers (buf_A, 0 , n, n * n);
91
90
random_integers (buf_x, 0 , n, n);
92
91
Matrix<int > A (n, n, buf_A);
93
92
Matrix<int > x (n, 1 , buf_x);
94
- Matrix<int > y (n, 1 , 0 );
95
- MatVec (A, x, y);
96
93
std::cout << A << std::endl;
97
94
std::cout << x << std::endl;
98
- std::cout << y << std::endl;
99
- if (compute) {
100
- Matrix<int > yw (n, 1 , 0 );
101
- MatVecWrong (A, x, yw);
102
- std::cout << yw << std::endl;
103
- std::cout << std::boolalpha << (y == yw) << std::endl;
95
+ Matrix<int > y1 (n, 1 , 0 ), y2 (n, 1 , 0 ), y3 (n, 1 , 0 );
96
+ if (compute >> 0 & 1 ) {
97
+ MatVec (A, x, y1);
98
+ std::cout << y1 << std::endl;
99
+ }
100
+ if (compute >> 1 & 1 ) {
101
+ MatVecRecursive (A, x, y2);
102
+ std::cout << y2;
103
+ if (compute >> 0 & 1 )
104
+ std::cout << std::boolalpha << (y1 == y2) << std::endl;
105
+ std::cout << std::endl;
106
+ }
107
+ if (compute >> 2 & 1 ) {
108
+ MatVecWrong (A, x, y3);
109
+ std::cout << y3;
110
+ if (compute >> 0 & 1 )
111
+ std::cout << std::boolalpha << (y1 == y3) << std::endl;
112
+ std::cout << std::endl;
104
113
}
105
114
return 0 ;
106
115
}
0 commit comments