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