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