ceph/src/boost/libs/numeric/ublas/include/boost/numeric/ublas/lu.hpp

   1 //
   2 //  Copyright (c) 2000-2002
   3 //  Joerg Walter, Mathias Koch
   4 //
   5 //  Distributed under the Boost Software License, Version 1.0. (See
   6 //  accompanying file LICENSE_1_0.txt or copy at
   7 //  http://www.boost.org/LICENSE_1_0.txt)
   8 //
   9 //  The authors gratefully acknowledge the support of
  10 //  GeNeSys mbH & Co. KG in producing this work.
  11 //
  12
  13 #ifndef _BOOST_UBLAS_LU_
  14 #define _BOOST_UBLAS_LU_
  15
  16 #include <boost/numeric/ublas/operation.hpp>
  17 #include <boost/numeric/ublas/vector_proxy.hpp>
  18 #include <boost/numeric/ublas/matrix_proxy.hpp>
  19 #include <boost/numeric/ublas/vector.hpp>
  20 #include <boost/numeric/ublas/triangular.hpp>
  21
  22 // LU factorizations in the spirit of LAPACK and Golub & van Loan
  23
  24 namespace boost { namespace numeric { namespace ublas {
  25
  26     /** \brief
  27      *
  28      * \tparam T
  29      * \tparam A
  30      */
  31     template<class T = std::size_t, class A = unbounded_array<T> >
  32     class permutation_matrix:
  33         public vector<T, A> {
  34     public:
  35         typedef vector<T, A> vector_type;
  36         typedef typename vector_type::size_type size_type;
  37
  38         // Construction and destruction
  39         BOOST_UBLAS_INLINE
  40         explicit
  41         permutation_matrix (size_type size):
  42             vector<T, A> (size) {
  43             for (size_type i = 0; i < size; ++ i)
  44                 (*this) (i) = i;
  45         }
  46         BOOST_UBLAS_INLINE
  47         explicit
  48         permutation_matrix (const vector_type & init)
  49             : vector_type(init)
  50         { }
  51         BOOST_UBLAS_INLINE
  52         ~permutation_matrix () {}
  53
  54         // Assignment
  55         BOOST_UBLAS_INLINE
  56         permutation_matrix &operator = (const permutation_matrix &m) {
  57             vector_type::operator = (m);
  58             return *this;
  59         }
  60     };
  61
  62     template<class PM, class MV>
  63     BOOST_UBLAS_INLINE
  64     void swap_rows (const PM &pm, MV &mv, vector_tag) {
  65         typedef typename PM::size_type size_type;
  66
  67         size_type size = pm.size ();
  68         for (size_type i = 0; i < size; ++ i) {
  69             if (i != pm (i))
  70                 std::swap (mv (i), mv (pm (i)));
  71         }
  72     }
  73     template<class PM, class MV>
  74     BOOST_UBLAS_INLINE
  75     void swap_rows (const PM &pm, MV &mv, matrix_tag) {
  76         typedef typename PM::size_type size_type;
  77
  78         size_type size = pm.size ();
  79         for (size_type i = 0; i < size; ++ i) {
  80             if (i != pm (i))
  81                 row (mv, i).swap (row (mv, pm (i)));
  82         }
  83     }
  84     // Dispatcher
  85     template<class PM, class MV>
  86     BOOST_UBLAS_INLINE
  87     void swap_rows (const PM &pm, MV &mv) {
  88         swap_rows (pm, mv, typename MV::type_category ());
  89     }
  90
  91     // LU factorization without pivoting
  92     template<class M>
  93     typename M::size_type lu_factorize (M &m) {
  94
  95         typedef typename M::size_type size_type;
  96         typedef typename M::value_type value_type;
  97
  98 #if BOOST_UBLAS_TYPE_CHECK
  99         typedef M matrix_type;
 100         matrix_type cm (m);
 101 #endif
 102         size_type singular = 0;
 103         size_type size1 = m.size1 ();
 104         size_type size2 = m.size2 ();
 105         size_type size = (std::min) (size1, size2);
 106         for (size_type i = 0; i < size; ++ i) {
 107             matrix_column<M> mci (column (m, i));
 108             matrix_row<M> mri (row (m, i));
 109             if (m (i, i) != value_type/*zero*/()) {
 110                 value_type m_inv = value_type (1) / m (i, i);
 111                 project (mci, range (i + 1, size1)) *= m_inv;
 112             } else if (singular == 0) {
 113                 singular = i + 1;
 114             }
 115             project (m, range (i + 1, size1), range (i + 1, size2)).minus_assign (
 116                 outer_prod (project (mci, range (i + 1, size1)),
 117                             project (mri, range (i + 1, size2))));
 118         }
 119 #if BOOST_UBLAS_TYPE_CHECK
 120         BOOST_UBLAS_CHECK (singular != 0 ||
 121                            detail::expression_type_check (prod (triangular_adaptor<matrix_type, unit_lower> (m),
 122                                                                 triangular_adaptor<matrix_type, upper> (m)),
 123                                                           cm), internal_logic ());
 124 #endif
 125         return singular;
 126     }
 127
 128     // LU factorization with partial pivoting
 129     template<class M, class PM>
 130     typename M::size_type lu_factorize (M &m, PM &pm) {
 131         typedef typename M::size_type size_type;
 132         typedef typename M::value_type value_type;
 133
 134 #if BOOST_UBLAS_TYPE_CHECK
 135         typedef M matrix_type;
 136         matrix_type cm (m);
 137 #endif
 138         size_type singular = 0;
 139         size_type size1 = m.size1 ();
 140         size_type size2 = m.size2 ();
 141         size_type size = (std::min) (size1, size2);
 142         for (size_type i = 0; i < size; ++ i) {
 143             matrix_column<M> mci (column (m, i));
 144             matrix_row<M> mri (row (m, i));
 145             size_type i_norm_inf = i + index_norm_inf (project (mci, range (i, size1)));
 146             BOOST_UBLAS_CHECK (i_norm_inf < size1, external_logic ());
 147             if (m (i_norm_inf, i) != value_type/*zero*/()) {
 148                 if (i_norm_inf != i) {
 149                     pm (i) = i_norm_inf;
 150                     row (m, i_norm_inf).swap (mri);
 151                 } else {
 152                     BOOST_UBLAS_CHECK (pm (i) == i_norm_inf, external_logic ());
 153                 }
 154                 value_type m_inv = value_type (1) / m (i, i);
 155                 project (mci, range (i + 1, size1)) *= m_inv;
 156             } else if (singular == 0) {
 157                 singular = i + 1;
 158             }
 159             project (m, range (i + 1, size1), range (i + 1, size2)).minus_assign (
 160                 outer_prod (project (mci, range (i + 1, size1)),
 161                             project (mri, range (i + 1, size2))));
 162         }
 163 #if BOOST_UBLAS_TYPE_CHECK
 164         swap_rows (pm, cm);
 165         BOOST_UBLAS_CHECK (singular != 0 ||
 166                            detail::expression_type_check (prod (triangular_adaptor<matrix_type, unit_lower> (m),
 167                                                                 triangular_adaptor<matrix_type, upper> (m)), cm), internal_logic ());
 168 #endif
 169         return singular;
 170     }
 171
 172     template<class M, class PM>
 173     typename M::size_type axpy_lu_factorize (M &m, PM &pm) {
 174         typedef M matrix_type;
 175         typedef typename M::size_type size_type;
 176         typedef typename M::value_type value_type;
 177         typedef vector<value_type> vector_type;
 178
 179 #if BOOST_UBLAS_TYPE_CHECK
 180         matrix_type cm (m);
 181 #endif
 182         size_type singular = 0;
 183         size_type size1 = m.size1 ();
 184         size_type size2 = m.size2 ();
 185         size_type size = (std::min) (size1, size2);
 186 #ifndef BOOST_UBLAS_LU_WITH_INPLACE_SOLVE
 187         matrix_type mr (m);
 188         mr.assign (zero_matrix<value_type> (size1, size2));
 189         vector_type v (size1);
 190         for (size_type i = 0; i < size; ++ i) {
 191             matrix_range<matrix_type> lrr (project (mr, range (0, i), range (0, i)));
 192             vector_range<matrix_column<matrix_type> > urr (project (column (mr, i), range (0, i)));
 193             urr.assign (solve (lrr, project (column (m, i), range (0, i)), unit_lower_tag ()));
 194             project (v, range (i, size1)).assign (
 195                 project (column (m, i), range (i, size1)) -
 196                 axpy_prod<vector_type> (project (mr, range (i, size1), range (0, i)), urr));
 197             size_type i_norm_inf = i + index_norm_inf (project (v, range (i, size1)));
 198             BOOST_UBLAS_CHECK (i_norm_inf < size1, external_logic ());
 199             if (v (i_norm_inf) != value_type/*zero*/()) {
 200                 if (i_norm_inf != i) {
 201                     pm (i) = i_norm_inf;
 202                     std::swap (v (i_norm_inf), v (i));
 203                     project (row (m, i_norm_inf), range (i + 1, size2)).swap (project (row (m, i), range (i + 1, size2)));
 204                 } else {
 205                     BOOST_UBLAS_CHECK (pm (i) == i_norm_inf, external_logic ());
 206                 }
 207                 project (column (mr, i), range (i + 1, size1)).assign (
 208                     project (v, range (i + 1, size1)) / v (i));
 209                 if (i_norm_inf != i) {
 210                     project (row (mr, i_norm_inf), range (0, i)).swap (project (row (mr, i), range (0, i)));
 211                 }
 212             } else if (singular == 0) {
 213                 singular = i + 1;
 214             }
 215             mr (i, i) = v (i);
 216         }
 217         m.assign (mr);
 218 #else
 219         matrix_type lr (m);
 220         matrix_type ur (m);
 221         lr.assign (identity_matrix<value_type> (size1, size2));
 222         ur.assign (zero_matrix<value_type> (size1, size2));
 223         vector_type v (size1);
 224         for (size_type i = 0; i < size; ++ i) {
 225             matrix_range<matrix_type> lrr (project (lr, range (0, i), range (0, i)));
 226             vector_range<matrix_column<matrix_type> > urr (project (column (ur, i), range (0, i)));
 227             urr.assign (project (column (m, i), range (0, i)));
 228             inplace_solve (lrr, urr, unit_lower_tag ());
 229             project (v, range (i, size1)).assign (
 230                 project (column (m, i), range (i, size1)) -
 231                 axpy_prod<vector_type> (project (lr, range (i, size1), range (0, i)), urr));
 232             size_type i_norm_inf = i + index_norm_inf (project (v, range (i, size1)));
 233             BOOST_UBLAS_CHECK (i_norm_inf < size1, external_logic ());
 234             if (v (i_norm_inf) != value_type/*zero*/()) {
 235                 if (i_norm_inf != i) {
 236                     pm (i) = i_norm_inf;
 237                     std::swap (v (i_norm_inf), v (i));
 238                     project (row (m, i_norm_inf), range (i + 1, size2)).swap (project (row (m, i), range (i + 1, size2)));
 239                 } else {
 240                     BOOST_UBLAS_CHECK (pm (i) == i_norm_inf, external_logic ());
 241                 }
 242                 project (column (lr, i), range (i + 1, size1)).assign (
 243                     project (v, range (i + 1, size1)) / v (i));
 244                 if (i_norm_inf != i) {
 245                     project (row (lr, i_norm_inf), range (0, i)).swap (project (row (lr, i), range (0, i)));
 246                 }
 247             } else if (singular == 0) {
 248                 singular = i + 1;
 249             }
 250             ur (i, i) = v (i);
 251         }
 252         m.assign (triangular_adaptor<matrix_type, strict_lower> (lr) +
 253                   triangular_adaptor<matrix_type, upper> (ur));
 254 #endif
 255 #if BOOST_UBLAS_TYPE_CHECK
 256         swap_rows (pm, cm);
 257         BOOST_UBLAS_CHECK (singular != 0 ||
 258                            detail::expression_type_check (prod (triangular_adaptor<matrix_type, unit_lower> (m),
 259                                                                 triangular_adaptor<matrix_type, upper> (m)), cm), internal_logic ());
 260 #endif
 261         return singular;
 262     }
 263
 264     // LU substitution
 265     template<class M, class E>
 266     void lu_substitute (const M &m, vector_expression<E> &e) {
 267 #if BOOST_UBLAS_TYPE_CHECK
 268         typedef const M const_matrix_type;
 269         typedef vector<typename E::value_type> vector_type;
 270
 271         vector_type cv1 (e);
 272 #endif
 273         inplace_solve (m, e, unit_lower_tag ());
 274 #if BOOST_UBLAS_TYPE_CHECK
 275         BOOST_UBLAS_CHECK (detail::expression_type_check (prod (triangular_adaptor<const_matrix_type, unit_lower> (m), e), cv1), internal_logic ());
 276         vector_type cv2 (e);
 277 #endif
 278         inplace_solve (m, e, upper_tag ());
 279 #if BOOST_UBLAS_TYPE_CHECK
 280         BOOST_UBLAS_CHECK (detail::expression_type_check (prod (triangular_adaptor<const_matrix_type, upper> (m), e), cv2), internal_logic ());
 281 #endif
 282     }
 283     template<class M, class E>
 284     void lu_substitute (const M &m, matrix_expression<E> &e) {
 285 #if BOOST_UBLAS_TYPE_CHECK
 286         typedef const M const_matrix_type;
 287         typedef matrix<typename E::value_type> matrix_type;
 288
 289         matrix_type cm1 (e);
 290 #endif
 291         inplace_solve (m, e, unit_lower_tag ());
 292 #if BOOST_UBLAS_TYPE_CHECK
 293         BOOST_UBLAS_CHECK (detail::expression_type_check (prod (triangular_adaptor<const_matrix_type, unit_lower> (m), e), cm1), internal_logic ());
 294         matrix_type cm2 (e);
 295 #endif
 296         inplace_solve (m, e, upper_tag ());
 297 #if BOOST_UBLAS_TYPE_CHECK
 298         BOOST_UBLAS_CHECK (detail::expression_type_check (prod (triangular_adaptor<const_matrix_type, upper> (m), e), cm2), internal_logic ());
 299 #endif
 300     }
 301     template<class M, class PMT, class PMA, class MV>
 302     void lu_substitute (const M &m, const permutation_matrix<PMT, PMA> &pm, MV &mv) {
 303         swap_rows (pm, mv);
 304         lu_substitute (m, mv);
 305     }
 306     template<class E, class M>
 307     void lu_substitute (vector_expression<E> &e, const M &m) {
 308 #if BOOST_UBLAS_TYPE_CHECK
 309         typedef const M const_matrix_type;
 310         typedef vector<typename E::value_type> vector_type;
 311
 312         vector_type cv1 (e);
 313 #endif
 314         inplace_solve (e, m, upper_tag ());
 315 #if BOOST_UBLAS_TYPE_CHECK
 316         BOOST_UBLAS_CHECK (detail::expression_type_check (prod (e, triangular_adaptor<const_matrix_type, upper> (m)), cv1), internal_logic ());
 317         vector_type cv2 (e);
 318 #endif
 319         inplace_solve (e, m, unit_lower_tag ());
 320 #if BOOST_UBLAS_TYPE_CHECK
 321         BOOST_UBLAS_CHECK (detail::expression_type_check (prod (e, triangular_adaptor<const_matrix_type, unit_lower> (m)), cv2), internal_logic ());
 322 #endif
 323     }
 324     template<class E, class M>
 325     void lu_substitute (matrix_expression<E> &e, const M &m) {
 326 #if BOOST_UBLAS_TYPE_CHECK
 327         typedef const M const_matrix_type;
 328         typedef matrix<typename E::value_type> matrix_type;
 329
 330         matrix_type cm1 (e);
 331 #endif
 332         inplace_solve (e, m, upper_tag ());
 333 #if BOOST_UBLAS_TYPE_CHECK
 334         BOOST_UBLAS_CHECK (detail::expression_type_check (prod (e, triangular_adaptor<const_matrix_type, upper> (m)), cm1), internal_logic ());
 335         matrix_type cm2 (e);
 336 #endif
 337         inplace_solve (e, m, unit_lower_tag ());
 338 #if BOOST_UBLAS_TYPE_CHECK
 339         BOOST_UBLAS_CHECK (detail::expression_type_check (prod (e, triangular_adaptor<const_matrix_type, unit_lower> (m)), cm2), internal_logic ());
 340 #endif
 341     }
 342     template<class MV, class M, class PMT, class PMA>
 343     void lu_substitute (MV &mv, const M &m, const permutation_matrix<PMT, PMA> &pm) {
 344         swap_rows (pm, mv);
 345         lu_substitute (mv, m);
 346     }
 347
 348 }}}
 349
 350 #endif