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 |
351 | |