numeric/ublas/lu.hpp

0001 //
0002 //  Copyright (c) 2000-2002
0003 //  Joerg Walter, Mathias Koch
0004 //
0005 //  Distributed under the Boost Software License, Version 1.0. (See
0006 //  accompanying file LICENSE_1_0.txt or copy at
0007 //  http://www.boost.org/LICENSE_1_0.txt)
0008 //
0009 //  The authors gratefully acknowledge the support of
0010 //  GeNeSys mbH & Co. KG in producing this work.
0011 //
0012
0013 #ifndef _BOOST_UBLAS_LU_
0014 #define _BOOST_UBLAS_LU_
0015
0016 #include <boost/numeric/ublas/operation.hpp>
0017 #include <boost/numeric/ublas/vector_proxy.hpp>
0018 #include <boost/numeric/ublas/matrix_proxy.hpp>
0019 #include <boost/numeric/ublas/vector.hpp>
0020 #include <boost/numeric/ublas/triangular.hpp>
0021
0022 // LU factorizations in the spirit of LAPACK and Golub & van Loan
0023
0024 namespace boost { namespace numeric { namespace ublas {
0025
0026     /** \brief
0027      *
0028      * \tparam T
0029      * \tparam A
0030      */
0031     template<class T = std::size_t, class A = unbounded_array<T> >
0032     class permutation_matrix:
0033         public vector<T, A> {
0034     public:
0035         typedef vector<T, A> vector_type;
0036         typedef typename vector_type::size_type size_type;
0037
0038         // Construction and destruction
0039         BOOST_UBLAS_INLINE
0040         explicit
0041         permutation_matrix (size_type size):
0042             vector<T, A> (size) {
0043             for (size_type i = 0; i < size; ++ i)
0044                 (*this) (i) = i;
0045         }
0046         BOOST_UBLAS_INLINE
0047         explicit
0048         permutation_matrix (const vector_type & init)
0049             : vector_type(init)
0050         { }
0051         BOOST_UBLAS_INLINE
0052         ~permutation_matrix () {}
0053
0054         // Assignment
0055         BOOST_UBLAS_INLINE
0056         permutation_matrix &operator = (const permutation_matrix &m) {
0057             vector_type::operator = (m);
0058             return *this;
0059         }
0060     };
0061
0062     template<class PM, class MV>
0063     BOOST_UBLAS_INLINE
0064     void swap_rows (const PM &pm, MV &mv, vector_tag) {
0065         typedef typename PM::size_type size_type;
0066
0067         size_type size = pm.size ();
0068         for (size_type i = 0; i < size; ++ i) {
0069             if (i != pm (i))
0070                 std::swap (mv (i), mv (pm (i)));
0071         }
0072     }
0073     template<class PM, class MV>
0074     BOOST_UBLAS_INLINE
0075     void swap_rows (const PM &pm, MV &mv, matrix_tag) {
0076         typedef typename PM::size_type size_type;
0077
0078         size_type size = pm.size ();
0079         for (size_type i = 0; i < size; ++ i) {
0080             if (i != pm (i))
0081                 row (mv, i).swap (row (mv, pm (i)));
0082         }
0083     }
0084     // Dispatcher
0085     template<class PM, class MV>
0086     BOOST_UBLAS_INLINE
0087     void swap_rows (const PM &pm, MV &mv) {
0088         swap_rows (pm, mv, typename MV::type_category ());
0089     }
0090
0091     // LU factorization without pivoting
0092     template<class M>
0093     typename M::size_type lu_factorize (M &m) {
0094
0095         typedef typename M::size_type size_type;
0096         typedef typename M::value_type value_type;
0097
0098 #if BOOST_UBLAS_TYPE_CHECK
0099         typedef M matrix_type;
0100         matrix_type cm (m);
0101 #endif
0102         size_type singular = 0;
0103         size_type size1 = m.size1 ();
0104         size_type size2 = m.size2 ();
0105         size_type size = (std::min) (size1, size2);
0106         for (size_type i = 0; i < size; ++ i) {
0107             matrix_column<M> mci (column (m, i));
0108             matrix_row<M> mri (row (m, i));
0109             if (m (i, i) != value_type/*zero*/()) {
0110                 value_type m_inv = value_type (1) / m (i, i);
0111                 project (mci, range (i + 1, size1)) *= m_inv;
0112             } else if (singular == 0) {
0113                 singular = i + 1;
0114             }
0115             project (m, range (i + 1, size1), range (i + 1, size2)).minus_assign (
0116                 outer_prod (project (mci, range (i + 1, size1)),
0117                             project (mri, range (i + 1, size2))));
0118         }
0119 #if BOOST_UBLAS_TYPE_CHECK
0120         BOOST_UBLAS_CHECK (singular != 0 ||
0121                            detail::expression_type_check (prod (triangular_adaptor<matrix_type, unit_lower> (m),
0122                                                                 triangular_adaptor<matrix_type, upper> (m)),
0123                                                           cm), internal_logic ());
0124 #endif
0125         return singular;
0126     }
0127
0128     // LU factorization with partial pivoting
0129     template<class M, class PM>
0130     typename M::size_type lu_factorize (M &m, PM &pm) {
0131         typedef typename M::size_type size_type;
0132         typedef typename M::value_type value_type;
0133
0134 #if BOOST_UBLAS_TYPE_CHECK
0135         typedef M matrix_type;
0136         matrix_type cm (m);
0137 #endif
0138         size_type singular = 0;
0139         size_type size1 = m.size1 ();
0140         size_type size2 = m.size2 ();
0141         size_type size = (std::min) (size1, size2);
0142         for (size_type i = 0; i < size; ++ i) {
0143             matrix_column<M> mci (column (m, i));
0144             matrix_row<M> mri (row (m, i));
0145             size_type i_norm_inf = i + index_norm_inf (project (mci, range (i, size1)));
0146             BOOST_UBLAS_CHECK (i_norm_inf < size1, external_logic ());
0147             if (m (i_norm_inf, i) != value_type/*zero*/()) {
0148                 if (i_norm_inf != i) {
0149                     pm (i) = i_norm_inf;
0150                     row (m, i_norm_inf).swap (mri);
0151                 } else {
0152                     BOOST_UBLAS_CHECK (pm (i) == i_norm_inf, external_logic ());
0153                 }
0154                 value_type m_inv = value_type (1) / m (i, i);
0155                 project (mci, range (i + 1, size1)) *= m_inv;
0156             } else if (singular == 0) {
0157                 singular = i + 1;
0158             }
0159             project (m, range (i + 1, size1), range (i + 1, size2)).minus_assign (
0160                 outer_prod (project (mci, range (i + 1, size1)),
0161                             project (mri, range (i + 1, size2))));
0162         }
0163 #if BOOST_UBLAS_TYPE_CHECK
0164         swap_rows (pm, cm);
0165         BOOST_UBLAS_CHECK (singular != 0 ||
0166                            detail::expression_type_check (prod (triangular_adaptor<matrix_type, unit_lower> (m),
0167                                                                 triangular_adaptor<matrix_type, upper> (m)), cm), internal_logic ());
0168 #endif
0169         return singular;
0170     }
0171
0172     template<class M, class PM>
0173     typename M::size_type axpy_lu_factorize (M &m, PM &pm) {
0174         typedef M matrix_type;
0175         typedef typename M::size_type size_type;
0176         typedef typename M::value_type value_type;
0177         typedef vector<value_type> vector_type;
0178
0179 #if BOOST_UBLAS_TYPE_CHECK
0180         matrix_type cm (m);
0181 #endif
0182         size_type singular = 0;
0183         size_type size1 = m.size1 ();
0184         size_type size2 = m.size2 ();
0185         size_type size = (std::min) (size1, size2);
0186 #ifndef BOOST_UBLAS_LU_WITH_INPLACE_SOLVE
0187         matrix_type mr (m);
0188         mr.assign (zero_matrix<value_type> (size1, size2));
0189         vector_type v (size1);
0190         for (size_type i = 0; i < size; ++ i) {
0191             matrix_range<matrix_type> lrr (project (mr, range (0, i), range (0, i)));
0192             vector_range<matrix_column<matrix_type> > urr (project (column (mr, i), range (0, i)));
0193             urr.assign (solve (lrr, project (column (m, i), range (0, i)), unit_lower_tag ()));
0194             project (v, range (i, size1)).assign (
0195                 project (column (m, i), range (i, size1)) -
0196                 axpy_prod<vector_type> (project (mr, range (i, size1), range (0, i)), urr));
0197             size_type i_norm_inf = i + index_norm_inf (project (v, range (i, size1)));
0198             BOOST_UBLAS_CHECK (i_norm_inf < size1, external_logic ());
0199             if (v (i_norm_inf) != value_type/*zero*/()) {
0200                 if (i_norm_inf != i) {
0201                     pm (i) = i_norm_inf;
0202                     std::swap (v (i_norm_inf), v (i));
0203                     project (row (m, i_norm_inf), range (i + 1, size2)).swap (project (row (m, i), range (i + 1, size2)));
0204                 } else {
0205                     BOOST_UBLAS_CHECK (pm (i) == i_norm_inf, external_logic ());
0206                 }
0207                 project (column (mr, i), range (i + 1, size1)).assign (
0208                     project (v, range (i + 1, size1)) / v (i));
0209                 if (i_norm_inf != i) {
0210                     project (row (mr, i_norm_inf), range (0, i)).swap (project (row (mr, i), range (0, i)));
0211                 }
0212             } else if (singular == 0) {
0213                 singular = i + 1;
0214             }
0215             mr (i, i) = v (i);
0216         }
0217         m.assign (mr);
0218 #else
0219         matrix_type lr (m);
0220         matrix_type ur (m);
0221         lr.assign (identity_matrix<value_type> (size1, size2));
0222         ur.assign (zero_matrix<value_type> (size1, size2));
0223         vector_type v (size1);
0224         for (size_type i = 0; i < size; ++ i) {
0225             matrix_range<matrix_type> lrr (project (lr, range (0, i), range (0, i)));
0226             vector_range<matrix_column<matrix_type> > urr (project (column (ur, i), range (0, i)));
0227             urr.assign (project (column (m, i), range (0, i)));
0228             inplace_solve (lrr, urr, unit_lower_tag ());
0229             project (v, range (i, size1)).assign (
0230                 project (column (m, i), range (i, size1)) -
0231                 axpy_prod<vector_type> (project (lr, range (i, size1), range (0, i)), urr));
0232             size_type i_norm_inf = i + index_norm_inf (project (v, range (i, size1)));
0233             BOOST_UBLAS_CHECK (i_norm_inf < size1, external_logic ());
0234             if (v (i_norm_inf) != value_type/*zero*/()) {
0235                 if (i_norm_inf != i) {
0236                     pm (i) = i_norm_inf;
0237                     std::swap (v (i_norm_inf), v (i));
0238                     project (row (m, i_norm_inf), range (i + 1, size2)).swap (project (row (m, i), range (i + 1, size2)));
0239                 } else {
0240                     BOOST_UBLAS_CHECK (pm (i) == i_norm_inf, external_logic ());
0241                 }
0242                 project (column (lr, i), range (i + 1, size1)).assign (
0243                     project (v, range (i + 1, size1)) / v (i));
0244                 if (i_norm_inf != i) {
0245                     project (row (lr, i_norm_inf), range (0, i)).swap (project (row (lr, i), range (0, i)));
0246                 }
0247             } else if (singular == 0) {
0248                 singular = i + 1;
0249             }
0250             ur (i, i) = v (i);
0251         }
0252         m.assign (triangular_adaptor<matrix_type, strict_lower> (lr) +
0253                   triangular_adaptor<matrix_type, upper> (ur));
0254 #endif
0255 #if BOOST_UBLAS_TYPE_CHECK
0256         swap_rows (pm, cm);
0257         BOOST_UBLAS_CHECK (singular != 0 ||
0258                            detail::expression_type_check (prod (triangular_adaptor<matrix_type, unit_lower> (m),
0259                                                                 triangular_adaptor<matrix_type, upper> (m)), cm), internal_logic ());
0260 #endif
0261         return singular;
0262     }
0263
0264     // LU substitution
0265     template<class M, class E>
0266     void lu_substitute (const M &m, vector_expression<E> &e) {
0267 #if BOOST_UBLAS_TYPE_CHECK
0268         typedef const M const_matrix_type;
0269         typedef vector<typename E::value_type> vector_type;
0270
0271         vector_type cv1 (e);
0272 #endif
0273         inplace_solve (m, e, unit_lower_tag ());
0274 #if BOOST_UBLAS_TYPE_CHECK
0275         BOOST_UBLAS_CHECK (detail::expression_type_check (prod (triangular_adaptor<const_matrix_type, unit_lower> (m), e), cv1), internal_logic ());
0276         vector_type cv2 (e);
0277 #endif
0278         inplace_solve (m, e, upper_tag ());
0279 #if BOOST_UBLAS_TYPE_CHECK
0280         BOOST_UBLAS_CHECK (detail::expression_type_check (prod (triangular_adaptor<const_matrix_type, upper> (m), e), cv2), internal_logic ());
0281 #endif
0282     }
0283     template<class M, class E>
0284     void lu_substitute (const M &m, matrix_expression<E> &e) {
0285 #if BOOST_UBLAS_TYPE_CHECK
0286         typedef const M const_matrix_type;
0287         typedef matrix<typename E::value_type> matrix_type;
0288
0289         matrix_type cm1 (e);
0290 #endif
0291         inplace_solve (m, e, unit_lower_tag ());
0292 #if BOOST_UBLAS_TYPE_CHECK
0293         BOOST_UBLAS_CHECK (detail::expression_type_check (prod (triangular_adaptor<const_matrix_type, unit_lower> (m), e), cm1), internal_logic ());
0294         matrix_type cm2 (e);
0295 #endif
0296         inplace_solve (m, e, upper_tag ());
0297 #if BOOST_UBLAS_TYPE_CHECK
0298         BOOST_UBLAS_CHECK (detail::expression_type_check (prod (triangular_adaptor<const_matrix_type, upper> (m), e), cm2), internal_logic ());
0299 #endif
0300     }
0301     template<class M, class PMT, class PMA, class MV>
0302     void lu_substitute (const M &m, const permutation_matrix<PMT, PMA> &pm, MV &mv) {
0303         swap_rows (pm, mv);
0304         lu_substitute (m, mv);
0305     }
0306     template<class E, class M>
0307     void lu_substitute (vector_expression<E> &e, const M &m) {
0308 #if BOOST_UBLAS_TYPE_CHECK
0309         typedef const M const_matrix_type;
0310         typedef vector<typename E::value_type> vector_type;
0311
0312         vector_type cv1 (e);
0313 #endif
0314         inplace_solve (e, m, upper_tag ());
0315 #if BOOST_UBLAS_TYPE_CHECK
0316         BOOST_UBLAS_CHECK (detail::expression_type_check (prod (e, triangular_adaptor<const_matrix_type, upper> (m)), cv1), internal_logic ());
0317         vector_type cv2 (e);
0318 #endif
0319         inplace_solve (e, m, unit_lower_tag ());
0320 #if BOOST_UBLAS_TYPE_CHECK
0321         BOOST_UBLAS_CHECK (detail::expression_type_check (prod (e, triangular_adaptor<const_matrix_type, unit_lower> (m)), cv2), internal_logic ());
0322 #endif
0323     }
0324     template<class E, class M>
0325     void lu_substitute (matrix_expression<E> &e, const M &m) {
0326 #if BOOST_UBLAS_TYPE_CHECK
0327         typedef const M const_matrix_type;
0328         typedef matrix<typename E::value_type> matrix_type;
0329
0330         matrix_type cm1 (e);
0331 #endif
0332         inplace_solve (e, m, upper_tag ());
0333 #if BOOST_UBLAS_TYPE_CHECK
0334         BOOST_UBLAS_CHECK (detail::expression_type_check (prod (e, triangular_adaptor<const_matrix_type, upper> (m)), cm1), internal_logic ());
0335         matrix_type cm2 (e);
0336 #endif
0337         inplace_solve (e, m, unit_lower_tag ());
0338 #if BOOST_UBLAS_TYPE_CHECK
0339         BOOST_UBLAS_CHECK (detail::expression_type_check (prod (e, triangular_adaptor<const_matrix_type, unit_lower> (m)), cm2), internal_logic ());
0340 #endif
0341     }
0342     template<class MV, class M, class PMT, class PMA>
0343     void lu_substitute (MV &mv, const M &m, const permutation_matrix<PMT, PMA> &pm) {
0344         swap_rows (pm, mv);
0345         lu_substitute (mv, m);
0346     }
0347
0348 }}}
0349
0350 #endif