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