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 permutation_matrix(size_type size)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 permutation_matrix(const vector_type & init)48 permutation_matrix (const vector_type & init) 49 : vector_type(init) 50 { } 51 BOOST_UBLAS_INLINE ~permutation_matrix()52 ~permutation_matrix () {} 53 54 // Assignment 55 BOOST_UBLAS_INLINE operator =(const permutation_matrix & m)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 swap_rows(const PM & pm,MV & mv,vector_tag)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 swap_rows(const PM & pm,MV & mv,matrix_tag)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 swap_rows(const PM & pm,MV & mv)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> lu_factorize(M & 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> lu_factorize(M & m,PM & 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> axpy_lu_factorize(M & m,PM & 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> lu_substitute(const M & m,vector_expression<E> & 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> lu_substitute(const M & m,matrix_expression<E> & 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> lu_substitute(const M & m,const permutation_matrix<PMT,PMA> & pm,MV & 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> lu_substitute(vector_expression<E> & e,const M & 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> lu_substitute(matrix_expression<E> & e,const M & 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> lu_substitute(MV & mv,const M & m,const permutation_matrix<PMT,PMA> & pm)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 351