1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2007, 2008 Klaus Spanderen
5
6 This file is part of QuantLib, a free-software/open-source library
7 for financial quantitative analysts and developers - http://quantlib.org/
8
9 QuantLib is free software: you can redistribute it and/or modify it
10 under the terms of the QuantLib license. You should have received a
11 copy of the license along with this program; if not, please email
12 <quantlib-dev@lists.sf.net>. The license is also available online at
13 <http://quantlib.org/license.shtml>.
14
15 This program is distributed in the hope that it will be useful, but WITHOUT
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 FOR A PARTICULAR PURPOSE. See the license for more details.
18*/
19
20/*! \file matrix.hpp
21 \brief matrix used in linear algebra.
22*/
23
24#include <ql/math/matrix.hpp>
25#if defined(QL_PATCH_MSVC)
26#pragma warning(push)
27#pragma warning(disable:4180)
28#pragma warning(disable:4127)
29#endif
30
31#if BOOST_VERSION == 106400
32#include <boost/serialization/array_wrapper.hpp>
33#endif
34#include <boost/numeric/ublas/vector_proxy.hpp>
35#include <boost/numeric/ublas/triangular.hpp>
36#include <boost/numeric/ublas/lu.hpp>
37
38#if defined(QL_PATCH_MSVC)
39#pragma warning(pop)
40#endif
41
42namespace QuantLib {
43
44 Matrix inverse(const Matrix& m) {
45 QL_REQUIRE(m.rows() == m.columns(), "matrix is not square");
46
47 boost::numeric::ublas::matrix<Real> a(m.rows(), m.columns());
48
49 std::copy(first: m.begin(), last: m.end(), result: a.data().begin());
50
51 boost::numeric::ublas::permutation_matrix<Size> pert(m.rows());
52
53 // lu decomposition
54 Size singular = 1;
55 try {
56 singular = lu_factorize(m&: a, pm&: pert);
57 } catch (const boost::numeric::ublas::internal_logic& e) {
58 QL_FAIL("lu_factorize error: " << e.what());
59 } catch (const boost::numeric::ublas::external_logic& e) {
60 QL_FAIL("lu_factorize error: " << e.what());
61 }
62 QL_REQUIRE(singular == 0, "singular matrix given");
63
64 boost::numeric::ublas::matrix<Real>
65 inverse = boost::numeric::ublas::identity_matrix<Real>(m.rows());
66
67 // backsubstitution
68 try {
69 boost::numeric::ublas::lu_substitute(m: a, pm: pert, mv&: inverse);
70 } catch (const boost::numeric::ublas::internal_logic& e) {
71 QL_FAIL("lu_substitute error: " << e.what());
72 }
73
74 Matrix retVal(m.rows(), m.columns());
75 std::copy(first: inverse.data().begin(), last: inverse.data().end(),
76 result: retVal.begin());
77
78 return retVal;
79 }
80
81 Real determinant(const Matrix& m) {
82 QL_REQUIRE(m.rows() == m.columns(), "matrix is not square");
83
84 boost::numeric::ublas::matrix<Real> a(m.rows(), m.columns());
85 std::copy(first: m.begin(), last: m.end(), result: a.data().begin());
86
87
88 // lu decomposition
89 boost::numeric::ublas::permutation_matrix<Size> pert(m.rows());
90 /* const Size singular = */ lu_factorize(m&: a, pm&: pert);
91
92 Real retVal = 1.0;
93
94 for (Size i=0; i < m.rows(); ++i) {
95 if (pert[i] != i)
96 retVal *= -a(i,i);
97 else
98 retVal *= a(i,i);
99 }
100 return retVal;
101 }
102
103}
104

source code of quantlib/ql/math/matrix.cpp