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