...one of the most highly
regarded and expertly designed C++ library projects in the
world.
— Herb Sutter and Andrei
Alexandrescu, C++
Coding Standards
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