1 | // Copyright 2008 Gunter Winkler <guwi17@gmx.de> |
2 | // Distributed under the Boost Software License, Version 1.0. (See |
3 | // accompanying file LICENSE_1_0.txt or copy at |
4 | // http://www.boost.org/LICENSE_1_0.txt) |
5 | |
6 | // switch automatic singular check off |
7 | #define BOOST_UBLAS_TYPE_CHECK 0 |
8 | |
9 | #include <boost/numeric/ublas/io.hpp> |
10 | #include <boost/numeric/ublas/lu.hpp> |
11 | #include <boost/cstdlib.hpp> |
12 | |
13 | #include "common/testhelper.hpp" |
14 | |
15 | #include <iostream> |
16 | #include <sstream> |
17 | |
18 | using namespace boost::numeric::ublas; |
19 | using std::string; |
20 | |
21 | static const string matrix_IN = "[3,3]((1,2,2),(2,3,3),(3,4,6))\0" ; |
22 | static const string matrix_LU = "[3,3]((3,4,6),(3.33333343e-01,6.66666627e-01,0),(6.66666687e-01,4.99999911e-01,-1))\0" ; |
23 | static const string matrix_INV= "[3,3]((-3,2,-7.94728621e-08),(1.50000012,0,-5.00000060e-01),(4.99999911e-01,-1,5.00000060e-01))\0" ; |
24 | static const string matrix_PM = "[3](2,2,2)" ; |
25 | |
26 | int main () { |
27 | |
28 | typedef float TYPE; |
29 | |
30 | typedef matrix<TYPE> MATRIX; |
31 | |
32 | MATRIX A; |
33 | MATRIX LU; |
34 | MATRIX INV; |
35 | |
36 | { |
37 | std::istringstream is(matrix_IN); |
38 | is >> A; |
39 | } |
40 | { |
41 | std::istringstream is(matrix_LU); |
42 | is >> LU; |
43 | } |
44 | { |
45 | std::istringstream is(matrix_INV); |
46 | is >> INV; |
47 | } |
48 | permutation_matrix<>::vector_type temp; |
49 | { |
50 | std::istringstream is(matrix_PM); |
51 | is >> temp; |
52 | } |
53 | permutation_matrix<> PM(temp); |
54 | |
55 | permutation_matrix<> pm(3); |
56 | |
57 | std::size_t result = lu_factorize<MATRIX, permutation_matrix<> >(m&: A, pm); |
58 | |
59 | assertTrue(message: "factorization completed: " , condition: 0 == result); |
60 | assertTrue(message: "LU factors are correct: " , condition: compare(m1: A, m2: LU)); |
61 | assertTrue(message: "permutation is correct: " , condition: compare(m1: pm, m2: PM)); |
62 | |
63 | MATRIX B = identity_matrix<TYPE>(A.size2()); |
64 | |
65 | lu_substitute(m: A, pm, mv&: B); |
66 | |
67 | assertTrue(message: "inverse is correct: " , condition: compare(m1: B, m2: INV)); |
68 | |
69 | return (getResults().second > 0) ? boost::exit_failure : boost::exit_success; |
70 | } |
71 | |