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
13namespace ublas = boost::numeric::ublas;
14
15template<class mat, class vec>
16double 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
25template<class mat, class vec>
26double diff(const vec& x, const mat& A, const vec& b) {
27 return diff(trans(A), x, b);
28}
29
30namespace ublas = boost::numeric::ublas;
31
32
33int 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

source code of boost/libs/numeric/ublas/test/test_triangular.cpp