1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/os/ossrv/ossrv_pub/boost_apis/boost/numeric/ublas/lu.hpp Fri Jun 15 03:10:57 2012 +0200
1.3 @@ -0,0 +1,342 @@
1.4 +//
1.5 +// Copyright (c) 2000-2002
1.6 +// Joerg Walter, Mathias Koch
1.7 +//
1.8 +// Permission to use, copy, modify, distribute and sell this software
1.9 +// and its documentation for any purpose is hereby granted without fee,
1.10 +// provided that the above copyright notice appear in all copies and
1.11 +// that both that copyright notice and this permission notice appear
1.12 +// in supporting documentation. The authors make no representations
1.13 +// about the suitability of this software for any purpose.
1.14 +// It is provided "as is" without express or implied warranty.
1.15 +//
1.16 +// The authors gratefully acknowledge the support of
1.17 +// GeNeSys mbH & Co. KG in producing this work.
1.18 +//
1.19 +
1.20 +#ifndef _BOOST_UBLAS_LU_
1.21 +#define _BOOST_UBLAS_LU_
1.22 +
1.23 +#include <boost/numeric/ublas/operation.hpp>
1.24 +#include <boost/numeric/ublas/vector_proxy.hpp>
1.25 +#include <boost/numeric/ublas/matrix_proxy.hpp>
1.26 +#include <boost/numeric/ublas/vector.hpp>
1.27 +#include <boost/numeric/ublas/triangular.hpp>
1.28 +
1.29 +// LU factorizations in the spirit of LAPACK and Golub & van Loan
1.30 +
1.31 +namespace boost { namespace numeric { namespace ublas {
1.32 +
1.33 + template<class T = std::size_t, class A = unbounded_array<T> >
1.34 + class permutation_matrix:
1.35 + public vector<T, A> {
1.36 + public:
1.37 + typedef vector<T, A> vector_type;
1.38 + typedef typename vector_type::size_type size_type;
1.39 +
1.40 + // Construction and destruction
1.41 + BOOST_UBLAS_INLINE
1.42 + permutation_matrix (size_type size):
1.43 + vector<T, A> (size) {
1.44 + for (size_type i = 0; i < size; ++ i)
1.45 + (*this) (i) = i;
1.46 + }
1.47 + BOOST_UBLAS_INLINE
1.48 + ~permutation_matrix () {}
1.49 +
1.50 + // Assignment
1.51 + BOOST_UBLAS_INLINE
1.52 + permutation_matrix &operator = (const permutation_matrix &m) {
1.53 + vector_type::operator = (m);
1.54 + return *this;
1.55 + }
1.56 + };
1.57 +
1.58 + template<class PM, class MV>
1.59 + BOOST_UBLAS_INLINE
1.60 + void swap_rows (const PM &pm, MV &mv, vector_tag) {
1.61 + typedef typename PM::size_type size_type;
1.62 + typedef typename MV::value_type value_type;
1.63 +
1.64 + size_type size = pm.size ();
1.65 + for (size_type i = 0; i < size; ++ i) {
1.66 + if (i != pm (i))
1.67 + std::swap (mv (i), mv (pm (i)));
1.68 + }
1.69 + }
1.70 + template<class PM, class MV>
1.71 + BOOST_UBLAS_INLINE
1.72 + void swap_rows (const PM &pm, MV &mv, matrix_tag) {
1.73 + typedef typename PM::size_type size_type;
1.74 + typedef typename MV::value_type value_type;
1.75 +
1.76 + size_type size = pm.size ();
1.77 + for (size_type i = 0; i < size; ++ i) {
1.78 + if (i != pm (i))
1.79 + row (mv, i).swap (row (mv, pm (i)));
1.80 + }
1.81 + }
1.82 + // Dispatcher
1.83 + template<class PM, class MV>
1.84 + BOOST_UBLAS_INLINE
1.85 + void swap_rows (const PM &pm, MV &mv) {
1.86 + swap_rows (pm, mv, typename MV::type_category ());
1.87 + }
1.88 +
1.89 + // LU factorization without pivoting
1.90 + template<class M>
1.91 + typename M::size_type lu_factorize (M &m) {
1.92 + typedef M matrix_type;
1.93 + typedef typename M::size_type size_type;
1.94 + typedef typename M::value_type value_type;
1.95 +
1.96 +#if BOOST_UBLAS_TYPE_CHECK
1.97 + matrix_type cm (m);
1.98 +#endif
1.99 + int singular = 0;
1.100 + size_type size1 = m.size1 ();
1.101 + size_type size2 = m.size2 ();
1.102 + size_type size = (std::min) (size1, size2);
1.103 + for (size_type i = 0; i < size; ++ i) {
1.104 + matrix_column<M> mci (column (m, i));
1.105 + matrix_row<M> mri (row (m, i));
1.106 + if (m (i, i) != value_type/*zero*/()) {
1.107 + project (mci, range (i + 1, size1)) *= value_type (1) / m (i, i);
1.108 + } else if (singular == 0) {
1.109 + singular = i + 1;
1.110 + }
1.111 + project (m, range (i + 1, size1), range (i + 1, size2)).minus_assign (
1.112 + outer_prod (project (mci, range (i + 1, size1)),
1.113 + project (mri, range (i + 1, size2))));
1.114 + }
1.115 +#if BOOST_UBLAS_TYPE_CHECK
1.116 + BOOST_UBLAS_CHECK (singular != 0 ||
1.117 + detail::expression_type_check (prod (triangular_adaptor<matrix_type, unit_lower> (m),
1.118 + triangular_adaptor<matrix_type, upper> (m)),
1.119 + cm), internal_logic ());
1.120 +#endif
1.121 + return singular;
1.122 + }
1.123 +
1.124 + // LU factorization with partial pivoting
1.125 + template<class M, class PM>
1.126 + typename M::size_type lu_factorize (M &m, PM &pm) {
1.127 + typedef M matrix_type;
1.128 + typedef typename M::size_type size_type;
1.129 + typedef typename M::value_type value_type;
1.130 +
1.131 +#if BOOST_UBLAS_TYPE_CHECK
1.132 + matrix_type cm (m);
1.133 +#endif
1.134 + int singular = 0;
1.135 + size_type size1 = m.size1 ();
1.136 + size_type size2 = m.size2 ();
1.137 + size_type size = (std::min) (size1, size2);
1.138 + for (size_type i = 0; i < size; ++ i) {
1.139 + matrix_column<M> mci (column (m, i));
1.140 + matrix_row<M> mri (row (m, i));
1.141 + size_type i_norm_inf = i + index_norm_inf (project (mci, range (i, size1)));
1.142 + BOOST_UBLAS_CHECK (i_norm_inf < size1, external_logic ());
1.143 + if (m (i_norm_inf, i) != value_type/*zero*/()) {
1.144 + if (i_norm_inf != i) {
1.145 + pm (i) = i_norm_inf;
1.146 + row (m, i_norm_inf).swap (mri);
1.147 + } else {
1.148 + BOOST_UBLAS_CHECK (pm (i) == i_norm_inf, external_logic ());
1.149 + }
1.150 + project (mci, range (i + 1, size1)) *= value_type (1) / m (i, i);
1.151 + } else if (singular == 0) {
1.152 + singular = i + 1;
1.153 + }
1.154 + project (m, range (i + 1, size1), range (i + 1, size2)).minus_assign (
1.155 + outer_prod (project (mci, range (i + 1, size1)),
1.156 + project (mri, range (i + 1, size2))));
1.157 + }
1.158 +#if BOOST_UBLAS_TYPE_CHECK
1.159 + swap_rows (pm, cm);
1.160 + BOOST_UBLAS_CHECK (singular != 0 ||
1.161 + detail::expression_type_check (prod (triangular_adaptor<matrix_type, unit_lower> (m),
1.162 + triangular_adaptor<matrix_type, upper> (m)), cm), internal_logic ());
1.163 +#endif
1.164 + return singular;
1.165 + }
1.166 +
1.167 + template<class M, class PM>
1.168 + typename M::size_type axpy_lu_factorize (M &m, PM &pm) {
1.169 + typedef M matrix_type;
1.170 + typedef typename M::size_type size_type;
1.171 + typedef typename M::value_type value_type;
1.172 + typedef vector<value_type> vector_type;
1.173 +
1.174 +#if BOOST_UBLAS_TYPE_CHECK
1.175 + matrix_type cm (m);
1.176 +#endif
1.177 + int singular = 0;
1.178 + size_type size1 = m.size1 ();
1.179 + size_type size2 = m.size2 ();
1.180 + size_type size = (std::min) (size1, size2);
1.181 +#ifndef BOOST_UBLAS_LU_WITH_INPLACE_SOLVE
1.182 + matrix_type mr (m);
1.183 + mr.assign (zero_matrix<value_type> (size1, size2));
1.184 + vector_type v (size1);
1.185 + for (size_type i = 0; i < size; ++ i) {
1.186 + matrix_range<matrix_type> lrr (project (mr, range (0, i), range (0, i)));
1.187 + vector_range<matrix_column<matrix_type> > urr (project (column (mr, i), range (0, i)));
1.188 + urr.assign (solve (lrr, project (column (m, i), range (0, i)), unit_lower_tag ()));
1.189 + project (v, range (i, size1)).assign (
1.190 + project (column (m, i), range (i, size1)) -
1.191 + axpy_prod<vector_type> (project (mr, range (i, size1), range (0, i)), urr));
1.192 + size_type i_norm_inf = i + index_norm_inf (project (v, range (i, size1)));
1.193 + BOOST_UBLAS_CHECK (i_norm_inf < size1, external_logic ());
1.194 + if (v (i_norm_inf) != value_type/*zero*/()) {
1.195 + if (i_norm_inf != i) {
1.196 + pm (i) = i_norm_inf;
1.197 + std::swap (v (i_norm_inf), v (i));
1.198 + project (row (m, i_norm_inf), range (i + 1, size2)).swap (project (row (m, i), range (i + 1, size2)));
1.199 + } else {
1.200 + BOOST_UBLAS_CHECK (pm (i) == i_norm_inf, external_logic ());
1.201 + }
1.202 + project (column (mr, i), range (i + 1, size1)).assign (
1.203 + project (v, range (i + 1, size1)) / v (i));
1.204 + if (i_norm_inf != i) {
1.205 + project (row (mr, i_norm_inf), range (0, i)).swap (project (row (mr, i), range (0, i)));
1.206 + }
1.207 + } else if (singular == 0) {
1.208 + singular = i + 1;
1.209 + }
1.210 + mr (i, i) = v (i);
1.211 + }
1.212 + m.assign (mr);
1.213 +#else
1.214 + matrix_type lr (m);
1.215 + matrix_type ur (m);
1.216 + lr.assign (identity_matrix<value_type> (size1, size2));
1.217 + ur.assign (zero_matrix<value_type> (size1, size2));
1.218 + vector_type v (size1);
1.219 + for (size_type i = 0; i < size; ++ i) {
1.220 + matrix_range<matrix_type> lrr (project (lr, range (0, i), range (0, i)));
1.221 + vector_range<matrix_column<matrix_type> > urr (project (column (ur, i), range (0, i)));
1.222 + urr.assign (project (column (m, i), range (0, i)));
1.223 + inplace_solve (lrr, urr, unit_lower_tag ());
1.224 + project (v, range (i, size1)).assign (
1.225 + project (column (m, i), range (i, size1)) -
1.226 + axpy_prod<vector_type> (project (lr, range (i, size1), range (0, i)), urr));
1.227 + size_type i_norm_inf = i + index_norm_inf (project (v, range (i, size1)));
1.228 + BOOST_UBLAS_CHECK (i_norm_inf < size1, external_logic ());
1.229 + if (v (i_norm_inf) != value_type/*zero*/()) {
1.230 + if (i_norm_inf != i) {
1.231 + pm (i) = i_norm_inf;
1.232 + std::swap (v (i_norm_inf), v (i));
1.233 + project (row (m, i_norm_inf), range (i + 1, size2)).swap (project (row (m, i), range (i + 1, size2)));
1.234 + } else {
1.235 + BOOST_UBLAS_CHECK (pm (i) == i_norm_inf, external_logic ());
1.236 + }
1.237 + project (column (lr, i), range (i + 1, size1)).assign (
1.238 + project (v, range (i + 1, size1)) / v (i));
1.239 + if (i_norm_inf != i) {
1.240 + project (row (lr, i_norm_inf), range (0, i)).swap (project (row (lr, i), range (0, i)));
1.241 + }
1.242 + } else if (singular == 0) {
1.243 + singular = i + 1;
1.244 + }
1.245 + ur (i, i) = v (i);
1.246 + }
1.247 + m.assign (triangular_adaptor<matrix_type, strict_lower> (lr) +
1.248 + triangular_adaptor<matrix_type, upper> (ur));
1.249 +#endif
1.250 +#if BOOST_UBLAS_TYPE_CHECK
1.251 + swap_rows (pm, cm);
1.252 + BOOST_UBLAS_CHECK (singular != 0 ||
1.253 + detail::expression_type_check (prod (triangular_adaptor<matrix_type, unit_lower> (m),
1.254 + triangular_adaptor<matrix_type, upper> (m)), cm), internal_logic ());
1.255 +#endif
1.256 + return singular;
1.257 + }
1.258 +
1.259 + // LU substitution
1.260 + template<class M, class E>
1.261 + void lu_substitute (const M &m, vector_expression<E> &e) {
1.262 + typedef const M const_matrix_type;
1.263 + typedef vector<typename E::value_type> vector_type;
1.264 +
1.265 +#if BOOST_UBLAS_TYPE_CHECK
1.266 + vector_type cv1 (e);
1.267 +#endif
1.268 + inplace_solve (m, e, unit_lower_tag ());
1.269 +#if BOOST_UBLAS_TYPE_CHECK
1.270 + BOOST_UBLAS_CHECK (detail::expression_type_check (prod (triangular_adaptor<const_matrix_type, unit_lower> (m), e), cv1), internal_logic ());
1.271 + vector_type cv2 (e);
1.272 +#endif
1.273 + inplace_solve (m, e, upper_tag ());
1.274 +#if BOOST_UBLAS_TYPE_CHECK
1.275 + BOOST_UBLAS_CHECK (detail::expression_type_check (prod (triangular_adaptor<const_matrix_type, upper> (m), e), cv2), internal_logic ());
1.276 +#endif
1.277 + }
1.278 + template<class M, class E>
1.279 + void lu_substitute (const M &m, matrix_expression<E> &e) {
1.280 + typedef const M const_matrix_type;
1.281 + typedef matrix<typename E::value_type> matrix_type;
1.282 +
1.283 +#if BOOST_UBLAS_TYPE_CHECK
1.284 + matrix_type cm1 (e);
1.285 +#endif
1.286 + inplace_solve (m, e, unit_lower_tag ());
1.287 +#if BOOST_UBLAS_TYPE_CHECK
1.288 + BOOST_UBLAS_CHECK (detail::expression_type_check (prod (triangular_adaptor<const_matrix_type, unit_lower> (m), e), cm1), internal_logic ());
1.289 + matrix_type cm2 (e);
1.290 +#endif
1.291 + inplace_solve (m, e, upper_tag ());
1.292 +#if BOOST_UBLAS_TYPE_CHECK
1.293 + BOOST_UBLAS_CHECK (detail::expression_type_check (prod (triangular_adaptor<const_matrix_type, upper> (m), e), cm2), internal_logic ());
1.294 +#endif
1.295 + }
1.296 + template<class M, class PMT, class PMA, class MV>
1.297 + void lu_substitute (const M &m, const permutation_matrix<PMT, PMA> &pm, MV &mv) {
1.298 + swap_rows (pm, mv);
1.299 + lu_substitute (m, mv);
1.300 + }
1.301 + template<class E, class M>
1.302 + void lu_substitute (vector_expression<E> &e, const M &m) {
1.303 + typedef const M const_matrix_type;
1.304 + typedef vector<typename E::value_type> vector_type;
1.305 +
1.306 +#if BOOST_UBLAS_TYPE_CHECK
1.307 + vector_type cv1 (e);
1.308 +#endif
1.309 + inplace_solve (e, m, upper_tag ());
1.310 +#if BOOST_UBLAS_TYPE_CHECK
1.311 + BOOST_UBLAS_CHECK (detail::expression_type_check (prod (e, triangular_adaptor<const_matrix_type, upper> (m)), cv1), internal_logic ());
1.312 + vector_type cv2 (e);
1.313 +#endif
1.314 + inplace_solve (e, m, unit_lower_tag ());
1.315 +#if BOOST_UBLAS_TYPE_CHECK
1.316 + BOOST_UBLAS_CHECK (detail::expression_type_check (prod (e, triangular_adaptor<const_matrix_type, unit_lower> (m)), cv2), internal_logic ());
1.317 +#endif
1.318 + }
1.319 + template<class E, class M>
1.320 + void lu_substitute (matrix_expression<E> &e, const M &m) {
1.321 + typedef const M const_matrix_type;
1.322 + typedef matrix<typename E::value_type> matrix_type;
1.323 +
1.324 +#if BOOST_UBLAS_TYPE_CHECK
1.325 + matrix_type cm1 (e);
1.326 +#endif
1.327 + inplace_solve (e, m, upper_tag ());
1.328 +#if BOOST_UBLAS_TYPE_CHECK
1.329 + BOOST_UBLAS_CHECK (detail::expression_type_check (prod (e, triangular_adaptor<const_matrix_type, upper> (m)), cm1), internal_logic ());
1.330 + matrix_type cm2 (e);
1.331 +#endif
1.332 + inplace_solve (e, m, unit_lower_tag ());
1.333 +#if BOOST_UBLAS_TYPE_CHECK
1.334 + BOOST_UBLAS_CHECK (detail::expression_type_check (prod (e, triangular_adaptor<const_matrix_type, unit_lower> (m)), cm2), internal_logic ());
1.335 +#endif
1.336 + }
1.337 + template<class MV, class M, class PMT, class PMA>
1.338 + void lu_substitute (MV &mv, const M &m, const permutation_matrix<PMT, PMA> &pm) {
1.339 + swap_rows (pm, mv);
1.340 + lu_substitute (mv, m);
1.341 + }
1.342 +
1.343 +}}}
1.344 +
1.345 +#endif