- Main Page
- Related Pages
- Classes
- Files

- File List

# lu.hpp

    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