]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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 |