1 | #include <iostream> |
2 | |
3 | #include <boost/numeric/ublas/vector.hpp> |
4 | #include <boost/numeric/ublas/matrix.hpp> |
5 | #include <boost/numeric/ublas/matrix_sparse.hpp> |
6 | #include <boost/numeric/ublas/triangular.hpp> |
7 | #include <boost/numeric/ublas/io.hpp> |
8 | |
9 | #include "utils.hpp" |
10 | |
11 | namespace ublas = boost::numeric::ublas; |
12 | |
13 | static const double TOL(1.0e-5); ///< Used for comparing two real numbers. |
14 | static const int n(10); ///< defines the test matrix size |
15 | |
16 | template<class mat, class vec> |
17 | double diff(const mat& A, const vec& x, const vec& b) { |
18 | return ublas::norm_2(prod(A, x) - b); |
19 | } |
20 | |
21 | // efficiently fill matrix depending on majority |
22 | template<class mat> |
23 | void fill_matrix(mat& A, ublas::column_major_tag) { |
24 | for (int i=0; i<n; ++i) { |
25 | if (i-1>=0) { |
26 | A(i-1, i) = -1; |
27 | } |
28 | A(i, i) = 1; |
29 | if (i+1<n) { |
30 | A(i+1, i) = -2; |
31 | } |
32 | } |
33 | } |
34 | template<class mat> |
35 | void fill_matrix(mat& A, ublas::row_major_tag) { |
36 | for (int i=0; i<n; ++i) { |
37 | if (i-1>=0) { |
38 | A(i, i-1) = -1; |
39 | } |
40 | A(i, i) = 1; |
41 | if (i+1<n) { |
42 | A(i, i+1) = -2; |
43 | } |
44 | } |
45 | } |
46 | |
47 | template<class mat> |
48 | BOOST_UBLAS_TEST_DEF ( test_inplace_solve ) |
49 | { |
50 | mat A(n, n); |
51 | A.clear(); |
52 | fill_matrix(A, typename mat::orientation_category()); |
53 | |
54 | ublas::vector<double> b(n, 1.0); |
55 | |
56 | // The test matrix is not triangular, but is interpreted that way by |
57 | // inplace_solve using the lower_tag/upper_tags. For checking, the |
58 | // triangular_adaptor makes A triangular for comparison. |
59 | { |
60 | ublas::vector<double> x(b); |
61 | ublas::inplace_solve(A, x, ublas::lower_tag()); |
62 | BOOST_UBLAS_TEST_CHECK(diff(ublas::triangular_adaptor<mat, ublas::lower>(A), x, b) < TOL); |
63 | } |
64 | { |
65 | ublas::vector<double> x(b); |
66 | ublas::inplace_solve(A, x, ublas::upper_tag()); |
67 | BOOST_UBLAS_TEST_CHECK(diff (ublas::triangular_adaptor<mat, ublas::upper>(A), x, b) < TOL); |
68 | } |
69 | { |
70 | ublas::vector<double> x(b); |
71 | ublas::inplace_solve(x, A, ublas::lower_tag()); |
72 | BOOST_UBLAS_TEST_CHECK(diff (trans(ublas::triangular_adaptor<mat, ublas::lower>(A)), x, b) < TOL); |
73 | } |
74 | { |
75 | ublas::vector<double> x(b); |
76 | ublas::inplace_solve(x, A, ublas::upper_tag()); |
77 | BOOST_UBLAS_TEST_CHECK(diff (trans(ublas::triangular_adaptor<mat, ublas::upper>(A)), x , b) < TOL); |
78 | } |
79 | } |
80 | |
81 | int main() { |
82 | |
83 | // typedefs are needed as macros do not work with "," in template arguments |
84 | |
85 | BOOST_UBLAS_TEST_BEGIN(); |
86 | |
87 | #ifdef USE_MATRIX |
88 | typedef ublas::matrix<double, ublas::row_major> mat_doub_rowmaj; |
89 | typedef ublas::matrix<double, ublas::column_major> mat_doub_colmaj; |
90 | BOOST_UBLAS_TEST_DO( test_inplace_solve<mat_doub_rowmaj> ); |
91 | BOOST_UBLAS_TEST_DO( test_inplace_solve<mat_doub_colmaj> ); |
92 | #endif |
93 | |
94 | #ifdef USE_COMPRESSED_MATRIX |
95 | typedef ublas::compressed_matrix<double, ublas::row_major> commat_doub_rowmaj; |
96 | typedef ublas::compressed_matrix<double, ublas::column_major> commat_doub_colmaj; |
97 | BOOST_UBLAS_TEST_DO( test_inplace_solve<commat_doub_rowmaj> ); |
98 | BOOST_UBLAS_TEST_DO( test_inplace_solve<commat_doub_colmaj> ); |
99 | #endif |
100 | |
101 | #ifdef USE_MAPPED_MATRIX |
102 | typedef ublas::mapped_matrix<double, ublas::row_major> mapmat_doub_rowmaj; |
103 | typedef ublas::mapped_matrix<double, ublas::column_major> mapmat_doub_colmaj; |
104 | BOOST_UBLAS_TEST_DO( test_inplace_solve<mapmat_doub_rowmaj> ); |
105 | BOOST_UBLAS_TEST_DO( test_inplace_solve<mapmat_doub_colmaj> ); |
106 | #endif |
107 | |
108 | #ifdef USE_COORDINATE_MATRIX |
109 | typedef ublas::coordinate_matrix<double, ublas::row_major> cormat_doub_rowmaj; |
110 | typedef ublas::coordinate_matrix<double, ublas::column_major> cormat_doub_colmaj; |
111 | BOOST_UBLAS_TEST_DO( test_inplace_solve<cormat_doub_rowmaj> ); |
112 | BOOST_UBLAS_TEST_DO( test_inplace_solve<cormat_doub_colmaj> ); |
113 | #endif |
114 | |
115 | #ifdef USE_MAPPED_VECTOR_OF_MAPPED_VECTOR |
116 | typedef ublas::mapped_vector_of_mapped_vector<double, ublas::row_major> mvmv_doub_rowmaj; |
117 | typedef ublas::mapped_vector_of_mapped_vector<double, ublas::column_major> mvmv_doub_colmaj; |
118 | BOOST_UBLAS_TEST_DO( test_inplace_solve<mvmv_doub_rowmaj> ); |
119 | BOOST_UBLAS_TEST_DO( test_inplace_solve<mvmv_doub_colmaj> ); |
120 | #endif |
121 | |
122 | BOOST_UBLAS_TEST_END(); |
123 | } |
124 | |