diff options
Diffstat (limited to 'src/boost/libs/numeric/ublas/test/test_triangular.cpp')
-rw-r--r-- | src/boost/libs/numeric/ublas/test/test_triangular.cpp | 129 |
1 files changed, 129 insertions, 0 deletions
diff --git a/src/boost/libs/numeric/ublas/test/test_triangular.cpp b/src/boost/libs/numeric/ublas/test/test_triangular.cpp new file mode 100644 index 000000000..c6cba9134 --- /dev/null +++ b/src/boost/libs/numeric/ublas/test/test_triangular.cpp @@ -0,0 +1,129 @@ +#include <iostream> +#include <stdlib.h> +#include <cmath> + +#include <boost/numeric/ublas/vector.hpp> +#include <boost/numeric/ublas/matrix.hpp> +#include <boost/numeric/ublas/matrix_sparse.hpp> +#include <boost/numeric/ublas/triangular.hpp> +#include <boost/numeric/ublas/io.hpp> + +#include <boost/timer/timer.hpp> + +namespace ublas = boost::numeric::ublas; + +template<class mat, class vec> +double diff(const mat& A, const vec& x, const vec& b) { + vec temp(prod(A, x) - b); + double result = 0; + for (typename vec::size_type i=0; i<temp.size(); ++i) { + result += temp(i)*temp(i); + } + return std::sqrt(result); +} + +template<class mat, class vec> +double diff(const vec& x, const mat& A, const vec& b) { + return diff(trans(A), x, b); +} + +namespace ublas = boost::numeric::ublas; + + +int main() { + const int n=7000; +#if 1 + ublas::compressed_matrix<double, ublas::row_major> mat_row_upp(n, n); + ublas::compressed_matrix<double, ublas::column_major> mat_col_upp(n, n); + ublas::compressed_matrix<double, ublas::row_major> mat_row_low(n, n); + ublas::compressed_matrix<double, ublas::column_major> mat_col_low(n, n); +#else + ublas::matrix<double, ublas::row_major> mat_row_upp(n, n, 0); + ublas::matrix<double, ublas::column_major> mat_col_upp(n, n, 0); + ublas::matrix<double, ublas::row_major> mat_row_low(n, n, 0); + ublas::matrix<double, ublas::column_major> mat_col_low(n, n, 0); +#endif + ublas::vector<double> b(n, 1); + + std::cerr << "Constructing..." << std::endl; + for (int i=0; i<n; ++i) { + b(i) = std::rand() % 10; + double main = -10 + std::rand() % 20 ; + if (main == 0) main+=1; + double side = -10 + std::rand() % 20 ; + if (i-1>=0) { + mat_row_low(i, i-1) = side; + } + mat_row_low(i, i) = main; + + mat_col_low(i, i) = main; + if (i+1<n) { + mat_col_low(i+1, i) = side; + } + + mat_row_upp(i, i) = main; + if (i+1<n) { + mat_row_upp(i, i+1) = side; + } + + if (i-1>=0) { + mat_col_upp(i-1, i) = side; + } + mat_col_upp(i, i) = main; + } + + std::cerr << "Starting..." << std::endl; + { + boost::timer::auto_cpu_timer t(std::cerr, "col_low x: %t sec CPU, %w sec real\n"); + ublas::vector<double> x(b); + ublas::inplace_solve(mat_col_low, x, ublas::lower_tag()); + std::cerr << "delta: " << diff(mat_col_low, x, b) << "\n"; + } + { + boost::timer::auto_cpu_timer t(std::cerr, "row_low x: %t sec CPU, %w sec real\n"); + ublas::vector<double> x(b); + ublas::inplace_solve(mat_row_low, x, ublas::lower_tag()); + std::cerr << "delta: " << diff(mat_row_low, x, b) << "\n"; + } + + { + boost::timer::auto_cpu_timer t(std::cerr, "col_upp x: %t sec CPU, %w sec real\n"); + ublas::vector<double> x(b); + ublas::inplace_solve(mat_col_upp, x, ublas::upper_tag()); + std::cerr << "delta: " << diff(mat_col_upp, x, b) << "\n"; + } + { + boost::timer::auto_cpu_timer t(std::cerr, "row_upp x: %t sec CPU, %w sec real\n"); + ublas::vector<double> x(b); + ublas::inplace_solve(mat_row_upp, x, ublas::upper_tag()); + std::cerr << "delta: " << diff(mat_row_upp, x, b) << "\n"; + } + + { + boost::timer::auto_cpu_timer t(std::cerr, "x col_low: %t sec CPU, %w sec real\n"); + ublas::vector<double> x(b); + ublas::inplace_solve(x, mat_col_low, ublas::lower_tag()); + std::cerr << "delta: " << diff(x, mat_col_low, b) << "\n"; + } + { + boost::timer::auto_cpu_timer t(std::cerr, "x row_low: %t sec CPU, %w sec real\n"); + ublas::vector<double> x(b); + ublas::inplace_solve(x, mat_row_low, ublas::lower_tag()); + std::cerr << "delta: " << diff(x, mat_row_low, b) << "\n"; + } + + { + boost::timer::auto_cpu_timer t(std::cerr, "x col_upp: %t sec CPU, %w sec real\n"); + ublas::vector<double> x(b); + ublas::inplace_solve(x, mat_col_upp, ublas::upper_tag()); + std::cerr << "delta: " << diff(x, mat_col_upp, b) << "\n"; + } + { + boost::timer::auto_cpu_timer t(std::cerr, "x row_upp: %t sec CPU, %w sec real\n"); + ublas::vector<double> x(b); + ublas::inplace_solve(x, mat_row_upp, ublas::upper_tag()); + std::cerr << "delta: " << diff(x, mat_row_upp, b) << "\n"; + } + + +} |