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
24namespace 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

source code of include/boost/numeric/ublas/lu.hpp