ublas: lu.hpp Source File

00001 //
00002 //  Copyright (c) 2000-2002
00003 //  Joerg Walter, Mathias Koch
00004 //
00005 //  Distributed under the Boost Software License, Version 1.0. (See
00006 //  accompanying file LICENSE_1_0.txt or copy at
00007 //  http://www.boost.org/LICENSE_1_0.txt)
00008 //
00009 //  The authors gratefully acknowledge the support of
00010 //  GeNeSys mbH & Co. KG in producing this work.
00011 //
00012 
00013 #ifndef _BOOST_UBLAS_LU_
00014 #define _BOOST_UBLAS_LU_
00015 
00016 #include <boost/numeric/ublas/operation.hpp>
00017 #include <boost/numeric/ublas/vector_proxy.hpp>
00018 #include <boost/numeric/ublas/matrix_proxy.hpp>
00019 #include <boost/numeric/ublas/vector.hpp>
00020 #include <boost/numeric/ublas/triangular.hpp>
00021 
00022 // LU factorizations in the spirit of LAPACK and Golub & van Loan
00023 
00024 namespace boost { namespace numeric { namespace ublas {
00025 
00031     template<class T = std::size_t, class A = unbounded_array<T> >
00032     class permutation_matrix:
00033         public vector<T, A> {
00034     public:
00035         typedef vector<T, A> vector_type;
00036         typedef typename vector_type::size_type size_type;
00037 
00038         // Construction and destruction
00039         BOOST_UBLAS_INLINE
00040         explicit
00041         permutation_matrix (size_type size):
00042             vector<T, A> (size) {
00043             for (size_type i = 0; i < size; ++ i)
00044                 (*this) (i) = i;
00045         }
00046         BOOST_UBLAS_INLINE
00047         explicit
00048         permutation_matrix (const vector_type & init) 
00049             : vector_type(init)
00050         { }
00051         BOOST_UBLAS_INLINE
00052         ~permutation_matrix () {}
00053 
00054         // Assignment
00055         BOOST_UBLAS_INLINE
00056         permutation_matrix &operator = (const permutation_matrix &m) {
00057             vector_type::operator = (m);
00058             return *this;
00059         }
00060     };
00061 
00062     template<class PM, class MV>
00063     BOOST_UBLAS_INLINE
00064     void swap_rows (const PM &pm, MV &mv, vector_tag) {
00065         typedef typename PM::size_type size_type;
00066         typedef typename MV::value_type value_type;
00067 
00068         size_type size = pm.size ();
00069         for (size_type i = 0; i < size; ++ i) {
00070             if (i != pm (i))
00071                 std::swap (mv (i), mv (pm (i)));
00072         }
00073     }
00074     template<class PM, class MV>
00075     BOOST_UBLAS_INLINE
00076     void swap_rows (const PM &pm, MV &mv, matrix_tag) {
00077         typedef typename PM::size_type size_type;
00078         typedef typename MV::value_type value_type;
00079 
00080         size_type size = pm.size ();
00081         for (size_type i = 0; i < size; ++ i) {
00082             if (i != pm (i))
00083                 row (mv, i).swap (row (mv, pm (i)));
00084         }
00085     }
00086     // Dispatcher
00087     template<class PM, class MV>
00088     BOOST_UBLAS_INLINE
00089     void swap_rows (const PM &pm, MV &mv) {
00090         swap_rows (pm, mv, typename MV::type_category ());
00091     }
00092 
00093     // LU factorization without pivoting
00094     template<class M>
00095     typename M::size_type lu_factorize (M &m) {
00096         typedef M matrix_type;
00097         typedef typename M::size_type size_type;
00098         typedef typename M::value_type value_type;
00099 
00100 #if BOOST_UBLAS_TYPE_CHECK
00101         matrix_type cm (m);
00102 #endif
00103         size_type singular = 0;
00104         size_type size1 = m.size1 ();
00105         size_type size2 = m.size2 ();
00106         size_type size = (std::min) (size1, size2);
00107         for (size_type i = 0; i < size; ++ i) {
00108             matrix_column<M> mci (column (m, i));
00109             matrix_row<M> mri (row (m, i));
00110             if (m (i, i) != value_type/*zero*/()) {
00111                 value_type m_inv = value_type (1) / m (i, i);
00112                 project (mci, range (i + 1, size1)) *= m_inv;
00113             } else if (singular == 0) {
00114                 singular = i + 1;
00115             }
00116             project (m, range (i + 1, size1), range (i + 1, size2)).minus_assign (
00117                 outer_prod (project (mci, range (i + 1, size1)),
00118                             project (mri, range (i + 1, size2))));
00119         }
00120 #if BOOST_UBLAS_TYPE_CHECK
00121         BOOST_UBLAS_CHECK (singular != 0 ||
00122                            detail::expression_type_check (prod (triangular_adaptor<matrix_type, unit_lower> (m),
00123                                                                 triangular_adaptor<matrix_type, upper> (m)), 
00124                                                           cm), internal_logic ());
00125 #endif
00126         return singular;
00127     }
00128 
00129     // LU factorization with partial pivoting
00130     template<class M, class PM>
00131     typename M::size_type lu_factorize (M &m, PM &pm) {
00132         typedef M matrix_type;
00133         typedef typename M::size_type size_type;
00134         typedef typename M::value_type value_type;
00135 
00136 #if BOOST_UBLAS_TYPE_CHECK
00137         matrix_type cm (m);
00138 #endif
00139         size_type singular = 0;
00140         size_type size1 = m.size1 ();
00141         size_type size2 = m.size2 ();
00142         size_type size = (std::min) (size1, size2);
00143         for (size_type i = 0; i < size; ++ i) {
00144             matrix_column<M> mci (column (m, i));
00145             matrix_row<M> mri (row (m, i));
00146             size_type i_norm_inf = i + index_norm_inf (project (mci, range (i, size1)));
00147             BOOST_UBLAS_CHECK (i_norm_inf < size1, external_logic ());
00148             if (m (i_norm_inf, i) != value_type/*zero*/()) {
00149                 if (i_norm_inf != i) {
00150                     pm (i) = i_norm_inf;
00151                     row (m, i_norm_inf).swap (mri);
00152                 } else {
00153                     BOOST_UBLAS_CHECK (pm (i) == i_norm_inf, external_logic ());
00154                 }
00155                 value_type m_inv = value_type (1) / m (i, i);
00156                 project (mci, range (i + 1, size1)) *= m_inv;
00157             } else if (singular == 0) {
00158                 singular = i + 1;
00159             }
00160             project (m, range (i + 1, size1), range (i + 1, size2)).minus_assign (
00161                 outer_prod (project (mci, range (i + 1, size1)),
00162                             project (mri, range (i + 1, size2))));
00163         }
00164 #if BOOST_UBLAS_TYPE_CHECK
00165         swap_rows (pm, cm);
00166         BOOST_UBLAS_CHECK (singular != 0 ||
00167                            detail::expression_type_check (prod (triangular_adaptor<matrix_type, unit_lower> (m),
00168                                                                 triangular_adaptor<matrix_type, upper> (m)), cm), internal_logic ());
00169 #endif
00170         return singular;
00171     }
00172 
00173     template<class M, class PM>
00174     typename M::size_type axpy_lu_factorize (M &m, PM &pm) {
00175         typedef M matrix_type;
00176         typedef typename M::size_type size_type;
00177         typedef typename M::value_type value_type;
00178         typedef vector<value_type> vector_type;
00179 
00180 #if BOOST_UBLAS_TYPE_CHECK
00181         matrix_type cm (m);
00182 #endif
00183         size_type singular = 0;
00184         size_type size1 = m.size1 ();
00185         size_type size2 = m.size2 ();
00186         size_type size = (std::min) (size1, size2);
00187 #ifndef BOOST_UBLAS_LU_WITH_INPLACE_SOLVE
00188         matrix_type mr (m);
00189         mr.assign (zero_matrix<value_type> (size1, size2));
00190         vector_type v (size1);
00191         for (size_type i = 0; i < size; ++ i) {
00192             matrix_range<matrix_type> lrr (project (mr, range (0, i), range (0, i)));
00193             vector_range<matrix_column<matrix_type> > urr (project (column (mr, i), range (0, i)));
00194             urr.assign (solve (lrr, project (column (m, i), range (0, i)), unit_lower_tag ()));
00195             project (v, range (i, size1)).assign (
00196                 project (column (m, i), range (i, size1)) -
00197                 axpy_prod<vector_type> (project (mr, range (i, size1), range (0, i)), urr));
00198             size_type i_norm_inf = i + index_norm_inf (project (v, range (i, size1)));
00199             BOOST_UBLAS_CHECK (i_norm_inf < size1, external_logic ());
00200             if (v (i_norm_inf) != value_type/*zero*/()) {
00201                 if (i_norm_inf != i) {
00202                     pm (i) = i_norm_inf;
00203                     std::swap (v (i_norm_inf), v (i));
00204                     project (row (m, i_norm_inf), range (i + 1, size2)).swap (project (row (m, i), range (i + 1, size2)));
00205                 } else {
00206                     BOOST_UBLAS_CHECK (pm (i) == i_norm_inf, external_logic ());
00207                 }
00208                 project (column (mr, i), range (i + 1, size1)).assign (
00209                     project (v, range (i + 1, size1)) / v (i));
00210                 if (i_norm_inf != i) {
00211                     project (row (mr, i_norm_inf), range (0, i)).swap (project (row (mr, i), range (0, i)));
00212                 }
00213             } else if (singular == 0) {
00214                 singular = i + 1;
00215             }
00216             mr (i, i) = v (i);
00217         }
00218         m.assign (mr);
00219 #else
00220         matrix_type lr (m);
00221         matrix_type ur (m);
00222         lr.assign (identity_matrix<value_type> (size1, size2));
00223         ur.assign (zero_matrix<value_type> (size1, size2));
00224         vector_type v (size1);
00225         for (size_type i = 0; i < size; ++ i) {
00226             matrix_range<matrix_type> lrr (project (lr, range (0, i), range (0, i)));
00227             vector_range<matrix_column<matrix_type> > urr (project (column (ur, i), range (0, i)));
00228             urr.assign (project (column (m, i), range (0, i)));
00229             inplace_solve (lrr, urr, unit_lower_tag ());
00230             project (v, range (i, size1)).assign (
00231                 project (column (m, i), range (i, size1)) -
00232                 axpy_prod<vector_type> (project (lr, range (i, size1), range (0, i)), urr));
00233             size_type i_norm_inf = i + index_norm_inf (project (v, range (i, size1)));
00234             BOOST_UBLAS_CHECK (i_norm_inf < size1, external_logic ());
00235             if (v (i_norm_inf) != value_type/*zero*/()) {
00236                 if (i_norm_inf != i) {
00237                     pm (i) = i_norm_inf;
00238                     std::swap (v (i_norm_inf), v (i));
00239                     project (row (m, i_norm_inf), range (i + 1, size2)).swap (project (row (m, i), range (i + 1, size2)));
00240                 } else {
00241                     BOOST_UBLAS_CHECK (pm (i) == i_norm_inf, external_logic ());
00242                 }
00243                 project (column (lr, i), range (i + 1, size1)).assign (
00244                     project (v, range (i + 1, size1)) / v (i));
00245                 if (i_norm_inf != i) {
00246                     project (row (lr, i_norm_inf), range (0, i)).swap (project (row (lr, i), range (0, i)));
00247                 }
00248             } else if (singular == 0) {
00249                 singular = i + 1;
00250             }
00251             ur (i, i) = v (i);
00252         }
00253         m.assign (triangular_adaptor<matrix_type, strict_lower> (lr) +
00254                   triangular_adaptor<matrix_type, upper> (ur));
00255 #endif
00256 #if BOOST_UBLAS_TYPE_CHECK
00257         swap_rows (pm, cm);
00258         BOOST_UBLAS_CHECK (singular != 0 ||
00259                            detail::expression_type_check (prod (triangular_adaptor<matrix_type, unit_lower> (m),
00260                                                                 triangular_adaptor<matrix_type, upper> (m)), cm), internal_logic ());
00261 #endif
00262         return singular;
00263     }
00264 
00265     // LU substitution
00266     template<class M, class E>
00267     void lu_substitute (const M &m, vector_expression<E> &e) {
00268         typedef const M const_matrix_type;
00269         typedef vector<typename E::value_type> vector_type;
00270 
00271 #if BOOST_UBLAS_TYPE_CHECK
00272         vector_type cv1 (e);
00273 #endif
00274         inplace_solve (m, e, unit_lower_tag ());
00275 #if BOOST_UBLAS_TYPE_CHECK
00276         BOOST_UBLAS_CHECK (detail::expression_type_check (prod (triangular_adaptor<const_matrix_type, unit_lower> (m), e), cv1), internal_logic ());
00277         vector_type cv2 (e);
00278 #endif
00279         inplace_solve (m, e, upper_tag ());
00280 #if BOOST_UBLAS_TYPE_CHECK
00281         BOOST_UBLAS_CHECK (detail::expression_type_check (prod (triangular_adaptor<const_matrix_type, upper> (m), e), cv2), internal_logic ());
00282 #endif
00283     }
00284     template<class M, class E>
00285     void lu_substitute (const M &m, matrix_expression<E> &e) {
00286         typedef const M const_matrix_type;
00287         typedef matrix<typename E::value_type> matrix_type;
00288 
00289 #if BOOST_UBLAS_TYPE_CHECK
00290         matrix_type cm1 (e);
00291 #endif
00292         inplace_solve (m, e, unit_lower_tag ());
00293 #if BOOST_UBLAS_TYPE_CHECK
00294         BOOST_UBLAS_CHECK (detail::expression_type_check (prod (triangular_adaptor<const_matrix_type, unit_lower> (m), e), cm1), internal_logic ());
00295         matrix_type cm2 (e);
00296 #endif
00297         inplace_solve (m, e, upper_tag ());
00298 #if BOOST_UBLAS_TYPE_CHECK
00299         BOOST_UBLAS_CHECK (detail::expression_type_check (prod (triangular_adaptor<const_matrix_type, upper> (m), e), cm2), internal_logic ());
00300 #endif
00301     }
00302     template<class M, class PMT, class PMA, class MV>
00303     void lu_substitute (const M &m, const permutation_matrix<PMT, PMA> &pm, MV &mv) {
00304         swap_rows (pm, mv);
00305         lu_substitute (m, mv);
00306     }
00307     template<class E, class M>
00308     void lu_substitute (vector_expression<E> &e, const M &m) {
00309         typedef const M const_matrix_type;
00310         typedef vector<typename E::value_type> vector_type;
00311 
00312 #if BOOST_UBLAS_TYPE_CHECK
00313         vector_type cv1 (e);
00314 #endif
00315         inplace_solve (e, m, upper_tag ());
00316 #if BOOST_UBLAS_TYPE_CHECK
00317         BOOST_UBLAS_CHECK (detail::expression_type_check (prod (e, triangular_adaptor<const_matrix_type, upper> (m)), cv1), internal_logic ());
00318         vector_type cv2 (e);
00319 #endif
00320         inplace_solve (e, m, unit_lower_tag ());
00321 #if BOOST_UBLAS_TYPE_CHECK
00322         BOOST_UBLAS_CHECK (detail::expression_type_check (prod (e, triangular_adaptor<const_matrix_type, unit_lower> (m)), cv2), internal_logic ());
00323 #endif
00324     }
00325     template<class E, class M>
00326     void lu_substitute (matrix_expression<E> &e, const M &m) {
00327         typedef const M const_matrix_type;
00328         typedef matrix<typename E::value_type> matrix_type;
00329 
00330 #if BOOST_UBLAS_TYPE_CHECK
00331         matrix_type cm1 (e);
00332 #endif
00333         inplace_solve (e, m, upper_tag ());
00334 #if BOOST_UBLAS_TYPE_CHECK
00335         BOOST_UBLAS_CHECK (detail::expression_type_check (prod (e, triangular_adaptor<const_matrix_type, upper> (m)), cm1), internal_logic ());
00336         matrix_type cm2 (e);
00337 #endif
00338         inplace_solve (e, m, unit_lower_tag ());
00339 #if BOOST_UBLAS_TYPE_CHECK
00340         BOOST_UBLAS_CHECK (detail::expression_type_check (prod (e, triangular_adaptor<const_matrix_type, unit_lower> (m)), cm2), internal_logic ());
00341 #endif
00342     }
00343     template<class MV, class M, class PMT, class PMA>
00344     void lu_substitute (MV &mv, const M &m, const permutation_matrix<PMT, PMA> &pm) {
00345         swap_rows (pm, mv);
00346         lu_substitute (mv, m);
00347     }
00348 
00349 }}}
00350 
00351 #endif