1 | #include <iostream> |
2 | #include <stdlib.h> |
3 | #include <cmath> |
4 | |
5 | #include <boost/numeric/ublas/vector.hpp> |
6 | #include <boost/numeric/ublas/matrix.hpp> |
7 | #include <boost/numeric/ublas/matrix_sparse.hpp> |
8 | #include <boost/numeric/ublas/triangular.hpp> |
9 | #include <boost/numeric/ublas/io.hpp> |
10 | |
11 | #include <boost/timer/timer.hpp> |
12 | |
13 | namespace ublas = boost::numeric::ublas; |
14 | |
15 | template<class mat, class vec> |
16 | double diff(const mat& A, const vec& x, const vec& b) { |
17 | vec temp(prod(A, x) - b); |
18 | double result = 0; |
19 | for (typename vec::size_type i=0; i<temp.size(); ++i) { |
20 | result += temp(i)*temp(i); |
21 | } |
22 | return std::sqrt(x: result); |
23 | } |
24 | |
25 | template<class mat, class vec> |
26 | double diff(const vec& x, const mat& A, const vec& b) { |
27 | return diff(trans(A), x, b); |
28 | } |
29 | |
30 | namespace ublas = boost::numeric::ublas; |
31 | |
32 | |
33 | int main() { |
34 | const int n=7000; |
35 | #if 1 |
36 | ublas::compressed_matrix<double, ublas::row_major> mat_row_upp(n, n); |
37 | ublas::compressed_matrix<double, ublas::column_major> mat_col_upp(n, n); |
38 | ublas::compressed_matrix<double, ublas::row_major> mat_row_low(n, n); |
39 | ublas::compressed_matrix<double, ublas::column_major> mat_col_low(n, n); |
40 | #else |
41 | ublas::matrix<double, ublas::row_major> mat_row_upp(n, n, 0); |
42 | ublas::matrix<double, ublas::column_major> mat_col_upp(n, n, 0); |
43 | ublas::matrix<double, ublas::row_major> mat_row_low(n, n, 0); |
44 | ublas::matrix<double, ublas::column_major> mat_col_low(n, n, 0); |
45 | #endif |
46 | ublas::vector<double> b(n, 1); |
47 | |
48 | std::cerr << "Constructing..." << std::endl; |
49 | for (int i=0; i<n; ++i) { |
50 | b(i) = std::rand() % 10; |
51 | double main = -10 + std::rand() % 20 ; |
52 | if (main == 0) main+=1; |
53 | double side = -10 + std::rand() % 20 ; |
54 | if (i-1>=0) { |
55 | mat_row_low(i, i-1) = side; |
56 | } |
57 | mat_row_low(i, i) = main; |
58 | |
59 | mat_col_low(i, i) = main; |
60 | if (i+1<n) { |
61 | mat_col_low(i+1, i) = side; |
62 | } |
63 | |
64 | mat_row_upp(i, i) = main; |
65 | if (i+1<n) { |
66 | mat_row_upp(i, i+1) = side; |
67 | } |
68 | |
69 | if (i-1>=0) { |
70 | mat_col_upp(i-1, i) = side; |
71 | } |
72 | mat_col_upp(i, i) = main; |
73 | } |
74 | |
75 | std::cerr << "Starting..." << std::endl; |
76 | { |
77 | boost::timer::auto_cpu_timer t(std::cerr, "col_low x: %t sec CPU, %w sec real\n" ); |
78 | ublas::vector<double> x(b); |
79 | ublas::inplace_solve(e1: mat_col_low, e2&: x, ublas::lower_tag()); |
80 | std::cerr << "delta: " << diff(A: mat_col_low, x, b) << "\n" ; |
81 | } |
82 | { |
83 | boost::timer::auto_cpu_timer t(std::cerr, "row_low x: %t sec CPU, %w sec real\n" ); |
84 | ublas::vector<double> x(b); |
85 | ublas::inplace_solve(e1: mat_row_low, e2&: x, ublas::lower_tag()); |
86 | std::cerr << "delta: " << diff(A: mat_row_low, x, b) << "\n" ; |
87 | } |
88 | |
89 | { |
90 | boost::timer::auto_cpu_timer t(std::cerr, "col_upp x: %t sec CPU, %w sec real\n" ); |
91 | ublas::vector<double> x(b); |
92 | ublas::inplace_solve(e1: mat_col_upp, e2&: x, ublas::upper_tag()); |
93 | std::cerr << "delta: " << diff(A: mat_col_upp, x, b) << "\n" ; |
94 | } |
95 | { |
96 | boost::timer::auto_cpu_timer t(std::cerr, "row_upp x: %t sec CPU, %w sec real\n" ); |
97 | ublas::vector<double> x(b); |
98 | ublas::inplace_solve(e1: mat_row_upp, e2&: x, ublas::upper_tag()); |
99 | std::cerr << "delta: " << diff(A: mat_row_upp, x, b) << "\n" ; |
100 | } |
101 | |
102 | { |
103 | boost::timer::auto_cpu_timer t(std::cerr, "x col_low: %t sec CPU, %w sec real\n" ); |
104 | ublas::vector<double> x(b); |
105 | ublas::inplace_solve(e1&: x, e2: mat_col_low, ublas::lower_tag()); |
106 | std::cerr << "delta: " << diff(x, A: mat_col_low, b) << "\n" ; |
107 | } |
108 | { |
109 | boost::timer::auto_cpu_timer t(std::cerr, "x row_low: %t sec CPU, %w sec real\n" ); |
110 | ublas::vector<double> x(b); |
111 | ublas::inplace_solve(e1&: x, e2: mat_row_low, ublas::lower_tag()); |
112 | std::cerr << "delta: " << diff(x, A: mat_row_low, b) << "\n" ; |
113 | } |
114 | |
115 | { |
116 | boost::timer::auto_cpu_timer t(std::cerr, "x col_upp: %t sec CPU, %w sec real\n" ); |
117 | ublas::vector<double> x(b); |
118 | ublas::inplace_solve(e1&: x, e2: mat_col_upp, ublas::upper_tag()); |
119 | std::cerr << "delta: " << diff(x, A: mat_col_upp, b) << "\n" ; |
120 | } |
121 | { |
122 | boost::timer::auto_cpu_timer t(std::cerr, "x row_upp: %t sec CPU, %w sec real\n" ); |
123 | ublas::vector<double> x(b); |
124 | ublas::inplace_solve(e1&: x, e2: mat_row_upp, ublas::upper_tag()); |
125 | std::cerr << "delta: " << diff(x, A: mat_row_upp, b) << "\n" ; |
126 | } |
127 | |
128 | |
129 | } |
130 | |