Chris@16: //
Chris@16: //  Copyright (c) 2000-2002
Chris@16: //  Joerg Walter, Mathias Koch
Chris@16: //
Chris@16: //  Distributed under the Boost Software License, Version 1.0. (See
Chris@16: //  accompanying file LICENSE_1_0.txt or copy at
Chris@16: //  http://www.boost.org/LICENSE_1_0.txt)
Chris@16: //
Chris@16: //  The authors gratefully acknowledge the support of
Chris@16: //  GeNeSys mbH & Co. KG in producing this work.
Chris@16: //
Chris@16: 
Chris@16: #ifndef _BOOST_UBLAS_LU_
Chris@16: #define _BOOST_UBLAS_LU_
Chris@16: 
Chris@16: #include <boost/numeric/ublas/operation.hpp>
Chris@16: #include <boost/numeric/ublas/vector_proxy.hpp>
Chris@16: #include <boost/numeric/ublas/matrix_proxy.hpp>
Chris@16: #include <boost/numeric/ublas/vector.hpp>
Chris@16: #include <boost/numeric/ublas/triangular.hpp>
Chris@16: 
Chris@16: // LU factorizations in the spirit of LAPACK and Golub & van Loan
Chris@16: 
Chris@16: namespace boost { namespace numeric { namespace ublas {
Chris@16: 
Chris@16:     /** \brief
Chris@16:      *
Chris@16:      * \tparam T
Chris@16:      * \tparam A
Chris@16:      */
Chris@16:     template<class T = std::size_t, class A = unbounded_array<T> >
Chris@16:     class permutation_matrix:
Chris@16:         public vector<T, A> {
Chris@16:     public:
Chris@16:         typedef vector<T, A> vector_type;
Chris@16:         typedef typename vector_type::size_type size_type;
Chris@16: 
Chris@16:         // Construction and destruction
Chris@16:         BOOST_UBLAS_INLINE
Chris@16:         explicit
Chris@16:         permutation_matrix (size_type size):
Chris@16:             vector<T, A> (size) {
Chris@16:             for (size_type i = 0; i < size; ++ i)
Chris@16:                 (*this) (i) = i;
Chris@16:         }
Chris@16:         BOOST_UBLAS_INLINE
Chris@16:         explicit
Chris@16:         permutation_matrix (const vector_type & init) 
Chris@16:             : vector_type(init)
Chris@16:         { }
Chris@16:         BOOST_UBLAS_INLINE
Chris@16:         ~permutation_matrix () {}
Chris@16: 
Chris@16:         // Assignment
Chris@16:         BOOST_UBLAS_INLINE
Chris@16:         permutation_matrix &operator = (const permutation_matrix &m) {
Chris@16:             vector_type::operator = (m);
Chris@16:             return *this;
Chris@16:         }
Chris@16:     };
Chris@16: 
Chris@16:     template<class PM, class MV>
Chris@16:     BOOST_UBLAS_INLINE
Chris@16:     void swap_rows (const PM &pm, MV &mv, vector_tag) {
Chris@16:         typedef typename PM::size_type size_type;
Chris@16: 
Chris@16:         size_type size = pm.size ();
Chris@16:         for (size_type i = 0; i < size; ++ i) {
Chris@16:             if (i != pm (i))
Chris@16:                 std::swap (mv (i), mv (pm (i)));
Chris@16:         }
Chris@16:     }
Chris@16:     template<class PM, class MV>
Chris@16:     BOOST_UBLAS_INLINE
Chris@16:     void swap_rows (const PM &pm, MV &mv, matrix_tag) {
Chris@16:         typedef typename PM::size_type size_type;
Chris@16: 
Chris@16:         size_type size = pm.size ();
Chris@16:         for (size_type i = 0; i < size; ++ i) {
Chris@16:             if (i != pm (i))
Chris@16:                 row (mv, i).swap (row (mv, pm (i)));
Chris@16:         }
Chris@16:     }
Chris@16:     // Dispatcher
Chris@16:     template<class PM, class MV>
Chris@16:     BOOST_UBLAS_INLINE
Chris@16:     void swap_rows (const PM &pm, MV &mv) {
Chris@16:         swap_rows (pm, mv, typename MV::type_category ());
Chris@16:     }
Chris@16: 
Chris@16:     // LU factorization without pivoting
Chris@16:     template<class M>
Chris@16:     typename M::size_type lu_factorize (M &m) {
Chris@101: 
Chris@16:         typedef typename M::size_type size_type;
Chris@16:         typedef typename M::value_type value_type;
Chris@16: 
Chris@16: #if BOOST_UBLAS_TYPE_CHECK
Chris@101:         typedef M matrix_type;
Chris@16:         matrix_type cm (m);
Chris@16: #endif
Chris@16:         size_type singular = 0;
Chris@16:         size_type size1 = m.size1 ();
Chris@16:         size_type size2 = m.size2 ();
Chris@16:         size_type size = (std::min) (size1, size2);
Chris@16:         for (size_type i = 0; i < size; ++ i) {
Chris@16:             matrix_column<M> mci (column (m, i));
Chris@16:             matrix_row<M> mri (row (m, i));
Chris@16:             if (m (i, i) != value_type/*zero*/()) {
Chris@16:                 value_type m_inv = value_type (1) / m (i, i);
Chris@16:                 project (mci, range (i + 1, size1)) *= m_inv;
Chris@16:             } else if (singular == 0) {
Chris@16:                 singular = i + 1;
Chris@16:             }
Chris@16:             project (m, range (i + 1, size1), range (i + 1, size2)).minus_assign (
Chris@16:                 outer_prod (project (mci, range (i + 1, size1)),
Chris@16:                             project (mri, range (i + 1, size2))));
Chris@16:         }
Chris@16: #if BOOST_UBLAS_TYPE_CHECK
Chris@16:         BOOST_UBLAS_CHECK (singular != 0 ||
Chris@16:                            detail::expression_type_check (prod (triangular_adaptor<matrix_type, unit_lower> (m),
Chris@16:                                                                 triangular_adaptor<matrix_type, upper> (m)), 
Chris@16:                                                           cm), internal_logic ());
Chris@16: #endif
Chris@16:         return singular;
Chris@16:     }
Chris@16: 
Chris@16:     // LU factorization with partial pivoting
Chris@16:     template<class M, class PM>
Chris@16:     typename M::size_type lu_factorize (M &m, PM &pm) {
Chris@16:         typedef typename M::size_type size_type;
Chris@16:         typedef typename M::value_type value_type;
Chris@16: 
Chris@16: #if BOOST_UBLAS_TYPE_CHECK
Chris@101:         typedef M matrix_type;
Chris@16:         matrix_type cm (m);
Chris@16: #endif
Chris@16:         size_type singular = 0;
Chris@16:         size_type size1 = m.size1 ();
Chris@16:         size_type size2 = m.size2 ();
Chris@16:         size_type size = (std::min) (size1, size2);
Chris@16:         for (size_type i = 0; i < size; ++ i) {
Chris@16:             matrix_column<M> mci (column (m, i));
Chris@16:             matrix_row<M> mri (row (m, i));
Chris@16:             size_type i_norm_inf = i + index_norm_inf (project (mci, range (i, size1)));
Chris@16:             BOOST_UBLAS_CHECK (i_norm_inf < size1, external_logic ());
Chris@16:             if (m (i_norm_inf, i) != value_type/*zero*/()) {
Chris@16:                 if (i_norm_inf != i) {
Chris@16:                     pm (i) = i_norm_inf;
Chris@16:                     row (m, i_norm_inf).swap (mri);
Chris@16:                 } else {
Chris@16:                     BOOST_UBLAS_CHECK (pm (i) == i_norm_inf, external_logic ());
Chris@16:                 }
Chris@16:                 value_type m_inv = value_type (1) / m (i, i);
Chris@16:                 project (mci, range (i + 1, size1)) *= m_inv;
Chris@16:             } else if (singular == 0) {
Chris@16:                 singular = i + 1;
Chris@16:             }
Chris@16:             project (m, range (i + 1, size1), range (i + 1, size2)).minus_assign (
Chris@16:                 outer_prod (project (mci, range (i + 1, size1)),
Chris@16:                             project (mri, range (i + 1, size2))));
Chris@16:         }
Chris@16: #if BOOST_UBLAS_TYPE_CHECK
Chris@16:         swap_rows (pm, cm);
Chris@16:         BOOST_UBLAS_CHECK (singular != 0 ||
Chris@16:                            detail::expression_type_check (prod (triangular_adaptor<matrix_type, unit_lower> (m),
Chris@16:                                                                 triangular_adaptor<matrix_type, upper> (m)), cm), internal_logic ());
Chris@16: #endif
Chris@16:         return singular;
Chris@16:     }
Chris@16: 
Chris@16:     template<class M, class PM>
Chris@16:     typename M::size_type axpy_lu_factorize (M &m, PM &pm) {
Chris@16:         typedef M matrix_type;
Chris@16:         typedef typename M::size_type size_type;
Chris@16:         typedef typename M::value_type value_type;
Chris@16:         typedef vector<value_type> vector_type;
Chris@16: 
Chris@16: #if BOOST_UBLAS_TYPE_CHECK
Chris@16:         matrix_type cm (m);
Chris@16: #endif
Chris@16:         size_type singular = 0;
Chris@16:         size_type size1 = m.size1 ();
Chris@16:         size_type size2 = m.size2 ();
Chris@16:         size_type size = (std::min) (size1, size2);
Chris@16: #ifndef BOOST_UBLAS_LU_WITH_INPLACE_SOLVE
Chris@16:         matrix_type mr (m);
Chris@16:         mr.assign (zero_matrix<value_type> (size1, size2));
Chris@16:         vector_type v (size1);
Chris@16:         for (size_type i = 0; i < size; ++ i) {
Chris@16:             matrix_range<matrix_type> lrr (project (mr, range (0, i), range (0, i)));
Chris@16:             vector_range<matrix_column<matrix_type> > urr (project (column (mr, i), range (0, i)));
Chris@16:             urr.assign (solve (lrr, project (column (m, i), range (0, i)), unit_lower_tag ()));
Chris@16:             project (v, range (i, size1)).assign (
Chris@16:                 project (column (m, i), range (i, size1)) -
Chris@16:                 axpy_prod<vector_type> (project (mr, range (i, size1), range (0, i)), urr));
Chris@16:             size_type i_norm_inf = i + index_norm_inf (project (v, range (i, size1)));
Chris@16:             BOOST_UBLAS_CHECK (i_norm_inf < size1, external_logic ());
Chris@16:             if (v (i_norm_inf) != value_type/*zero*/()) {
Chris@16:                 if (i_norm_inf != i) {
Chris@16:                     pm (i) = i_norm_inf;
Chris@16:                     std::swap (v (i_norm_inf), v (i));
Chris@16:                     project (row (m, i_norm_inf), range (i + 1, size2)).swap (project (row (m, i), range (i + 1, size2)));
Chris@16:                 } else {
Chris@16:                     BOOST_UBLAS_CHECK (pm (i) == i_norm_inf, external_logic ());
Chris@16:                 }
Chris@16:                 project (column (mr, i), range (i + 1, size1)).assign (
Chris@16:                     project (v, range (i + 1, size1)) / v (i));
Chris@16:                 if (i_norm_inf != i) {
Chris@16:                     project (row (mr, i_norm_inf), range (0, i)).swap (project (row (mr, i), range (0, i)));
Chris@16:                 }
Chris@16:             } else if (singular == 0) {
Chris@16:                 singular = i + 1;
Chris@16:             }
Chris@16:             mr (i, i) = v (i);
Chris@16:         }
Chris@16:         m.assign (mr);
Chris@16: #else
Chris@16:         matrix_type lr (m);
Chris@16:         matrix_type ur (m);
Chris@16:         lr.assign (identity_matrix<value_type> (size1, size2));
Chris@16:         ur.assign (zero_matrix<value_type> (size1, size2));
Chris@16:         vector_type v (size1);
Chris@16:         for (size_type i = 0; i < size; ++ i) {
Chris@16:             matrix_range<matrix_type> lrr (project (lr, range (0, i), range (0, i)));
Chris@16:             vector_range<matrix_column<matrix_type> > urr (project (column (ur, i), range (0, i)));
Chris@16:             urr.assign (project (column (m, i), range (0, i)));
Chris@16:             inplace_solve (lrr, urr, unit_lower_tag ());
Chris@16:             project (v, range (i, size1)).assign (
Chris@16:                 project (column (m, i), range (i, size1)) -
Chris@16:                 axpy_prod<vector_type> (project (lr, range (i, size1), range (0, i)), urr));
Chris@16:             size_type i_norm_inf = i + index_norm_inf (project (v, range (i, size1)));
Chris@16:             BOOST_UBLAS_CHECK (i_norm_inf < size1, external_logic ());
Chris@16:             if (v (i_norm_inf) != value_type/*zero*/()) {
Chris@16:                 if (i_norm_inf != i) {
Chris@16:                     pm (i) = i_norm_inf;
Chris@16:                     std::swap (v (i_norm_inf), v (i));
Chris@16:                     project (row (m, i_norm_inf), range (i + 1, size2)).swap (project (row (m, i), range (i + 1, size2)));
Chris@16:                 } else {
Chris@16:                     BOOST_UBLAS_CHECK (pm (i) == i_norm_inf, external_logic ());
Chris@16:                 }
Chris@16:                 project (column (lr, i), range (i + 1, size1)).assign (
Chris@16:                     project (v, range (i + 1, size1)) / v (i));
Chris@16:                 if (i_norm_inf != i) {
Chris@16:                     project (row (lr, i_norm_inf), range (0, i)).swap (project (row (lr, i), range (0, i)));
Chris@16:                 }
Chris@16:             } else if (singular == 0) {
Chris@16:                 singular = i + 1;
Chris@16:             }
Chris@16:             ur (i, i) = v (i);
Chris@16:         }
Chris@16:         m.assign (triangular_adaptor<matrix_type, strict_lower> (lr) +
Chris@16:                   triangular_adaptor<matrix_type, upper> (ur));
Chris@16: #endif
Chris@16: #if BOOST_UBLAS_TYPE_CHECK
Chris@16:         swap_rows (pm, cm);
Chris@16:         BOOST_UBLAS_CHECK (singular != 0 ||
Chris@16:                            detail::expression_type_check (prod (triangular_adaptor<matrix_type, unit_lower> (m),
Chris@16:                                                                 triangular_adaptor<matrix_type, upper> (m)), cm), internal_logic ());
Chris@16: #endif
Chris@16:         return singular;
Chris@16:     }
Chris@16: 
Chris@16:     // LU substitution
Chris@16:     template<class M, class E>
Chris@16:     void lu_substitute (const M &m, vector_expression<E> &e) {
Chris@101: #if BOOST_UBLAS_TYPE_CHECK
Chris@16:         typedef const M const_matrix_type;
Chris@16:         typedef vector<typename E::value_type> vector_type;
Chris@16: 
Chris@16:         vector_type cv1 (e);
Chris@16: #endif
Chris@16:         inplace_solve (m, e, unit_lower_tag ());
Chris@16: #if BOOST_UBLAS_TYPE_CHECK
Chris@16:         BOOST_UBLAS_CHECK (detail::expression_type_check (prod (triangular_adaptor<const_matrix_type, unit_lower> (m), e), cv1), internal_logic ());
Chris@16:         vector_type cv2 (e);
Chris@16: #endif
Chris@16:         inplace_solve (m, e, upper_tag ());
Chris@16: #if BOOST_UBLAS_TYPE_CHECK
Chris@16:         BOOST_UBLAS_CHECK (detail::expression_type_check (prod (triangular_adaptor<const_matrix_type, upper> (m), e), cv2), internal_logic ());
Chris@16: #endif
Chris@16:     }
Chris@16:     template<class M, class E>
Chris@16:     void lu_substitute (const M &m, matrix_expression<E> &e) {
Chris@101: #if BOOST_UBLAS_TYPE_CHECK
Chris@16:         typedef const M const_matrix_type;
Chris@16:         typedef matrix<typename E::value_type> matrix_type;
Chris@16: 
Chris@16:         matrix_type cm1 (e);
Chris@16: #endif
Chris@16:         inplace_solve (m, e, unit_lower_tag ());
Chris@16: #if BOOST_UBLAS_TYPE_CHECK
Chris@16:         BOOST_UBLAS_CHECK (detail::expression_type_check (prod (triangular_adaptor<const_matrix_type, unit_lower> (m), e), cm1), internal_logic ());
Chris@16:         matrix_type cm2 (e);
Chris@16: #endif
Chris@16:         inplace_solve (m, e, upper_tag ());
Chris@16: #if BOOST_UBLAS_TYPE_CHECK
Chris@16:         BOOST_UBLAS_CHECK (detail::expression_type_check (prod (triangular_adaptor<const_matrix_type, upper> (m), e), cm2), internal_logic ());
Chris@16: #endif
Chris@16:     }
Chris@16:     template<class M, class PMT, class PMA, class MV>
Chris@16:     void lu_substitute (const M &m, const permutation_matrix<PMT, PMA> &pm, MV &mv) {
Chris@16:         swap_rows (pm, mv);
Chris@16:         lu_substitute (m, mv);
Chris@16:     }
Chris@16:     template<class E, class M>
Chris@16:     void lu_substitute (vector_expression<E> &e, const M &m) {
Chris@101: #if BOOST_UBLAS_TYPE_CHECK
Chris@16:         typedef const M const_matrix_type;
Chris@16:         typedef vector<typename E::value_type> vector_type;
Chris@16: 
Chris@16:         vector_type cv1 (e);
Chris@16: #endif
Chris@16:         inplace_solve (e, m, upper_tag ());
Chris@16: #if BOOST_UBLAS_TYPE_CHECK
Chris@16:         BOOST_UBLAS_CHECK (detail::expression_type_check (prod (e, triangular_adaptor<const_matrix_type, upper> (m)), cv1), internal_logic ());
Chris@16:         vector_type cv2 (e);
Chris@16: #endif
Chris@16:         inplace_solve (e, m, unit_lower_tag ());
Chris@16: #if BOOST_UBLAS_TYPE_CHECK
Chris@16:         BOOST_UBLAS_CHECK (detail::expression_type_check (prod (e, triangular_adaptor<const_matrix_type, unit_lower> (m)), cv2), internal_logic ());
Chris@16: #endif
Chris@16:     }
Chris@16:     template<class E, class M>
Chris@16:     void lu_substitute (matrix_expression<E> &e, const M &m) {
Chris@101: #if BOOST_UBLAS_TYPE_CHECK
Chris@16:         typedef const M const_matrix_type;
Chris@16:         typedef matrix<typename E::value_type> matrix_type;
Chris@16: 
Chris@16:         matrix_type cm1 (e);
Chris@16: #endif
Chris@16:         inplace_solve (e, m, upper_tag ());
Chris@16: #if BOOST_UBLAS_TYPE_CHECK
Chris@16:         BOOST_UBLAS_CHECK (detail::expression_type_check (prod (e, triangular_adaptor<const_matrix_type, upper> (m)), cm1), internal_logic ());
Chris@16:         matrix_type cm2 (e);
Chris@16: #endif
Chris@16:         inplace_solve (e, m, unit_lower_tag ());
Chris@16: #if BOOST_UBLAS_TYPE_CHECK
Chris@16:         BOOST_UBLAS_CHECK (detail::expression_type_check (prod (e, triangular_adaptor<const_matrix_type, unit_lower> (m)), cm2), internal_logic ());
Chris@16: #endif
Chris@16:     }
Chris@16:     template<class MV, class M, class PMT, class PMA>
Chris@16:     void lu_substitute (MV &mv, const M &m, const permutation_matrix<PMT, PMA> &pm) {
Chris@16:         swap_rows (pm, mv);
Chris@16:         lu_substitute (mv, m);
Chris@16:     }
Chris@16: 
Chris@16: }}}
Chris@16: 
Chris@16: #endif