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
main()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<> >(A, pm);
58
59 assertTrue("factorization completed: ", 0 == result);
60 assertTrue("LU factors are correct: ", compare(A, LU));
61 assertTrue("permutation is correct: ", compare(pm, PM));
62
63 MATRIX B = identity_matrix<TYPE>(A.size2());
64
65 lu_substitute(A, pm, B);
66
67 assertTrue("inverse is correct: ", compare(B, INV));
68
69 return (getResults().second > 0) ? boost::exit_failure : boost::exit_success;
70 }
71