source: NonGTP/Boost/boost/numeric/ublas/lu.hpp @ 857

Revision 857, 14.1 KB checked in by igarcia, 19 years ago (diff)
RevLine 
[857]1//
2//  Copyright (c) 2000-2002
3//  Joerg Walter, Mathias Koch
4//
5//  Permission to use, copy, modify, distribute and sell this software
6//  and its documentation for any purpose is hereby granted without fee,
7//  provided that the above copyright notice appear in all copies and
8//  that both that copyright notice and this permission notice appear
9//  in supporting documentation.  The authors make no representations
10//  about the suitability of this software for any purpose.
11//  It is provided "as is" without express or implied warranty.
12//
13//  The authors gratefully acknowledge the support of
14//  GeNeSys mbH & Co. KG in producing this work.
15//
16
17#ifndef _BOOST_UBLAS_LU_
18#define _BOOST_UBLAS_LU_
19
20#include <boost/numeric/ublas/operation.hpp>
21#include <boost/numeric/ublas/matrix_expression.hpp>
22#include <boost/numeric/ublas/vector.hpp>
23#include <boost/numeric/ublas/detail/matrix_assign.hpp>
24
25// LU factorizations in the spirit of LAPACK and Golub & van Loan
26
27namespace boost { namespace numeric { namespace ublas {
28
29    template<class T = std::size_t, class A = unbounded_array<T> >
30    class permutation_matrix:
31        public vector<T, A> {
32    public:
33        typedef vector<T, A> vector_type;
34        typedef typename vector_type::size_type size_type;
35
36        // Construction and destruction
37        BOOST_UBLAS_INLINE
38        permutation_matrix (size_type size):
39            vector<T, A> (size) {
40            for (size_type i = 0; i < size; ++ i)
41                (*this) (i) = i;
42        }
43        BOOST_UBLAS_INLINE
44        ~permutation_matrix () {}
45
46        // Assignment
47        BOOST_UBLAS_INLINE
48        permutation_matrix &operator = (const permutation_matrix &m) {
49            vector_type::operator = (m);
50            return *this;
51        }
52    };
53
54    template<class PM, class MV>
55    BOOST_UBLAS_INLINE
56    void swap_rows (const PM &pm, MV &mv, vector_tag) {
57        typedef typename PM::size_type size_type;
58        typedef typename MV::value_type value_type;
59
60        size_type size = pm.size ();
61        for (size_type i = 0; i < size; ++ i) {
62            if (i != pm (i))
63                std::swap (mv (i), mv (pm (i)));
64        }
65    }
66    template<class PM, class MV>
67    BOOST_UBLAS_INLINE
68    void swap_rows (const PM &pm, MV &mv, matrix_tag) {
69        typedef typename PM::size_type size_type;
70        typedef typename MV::value_type value_type;
71
72        size_type size = pm.size ();
73        for (size_type i = 0; i < size; ++ i) {
74            if (i != pm (i))
75                row (mv, i).swap (row (mv, pm (i)));
76        }
77    }
78    // Dispatcher
79    template<class PM, class MV>
80    BOOST_UBLAS_INLINE
81    void swap_rows (const PM &pm, MV &mv) {
82        swap_rows (pm, mv, typename MV::type_category ());
83    }
84
85    // LU factorization without pivoting
86    template<class M>
87    typename M::size_type lu_factorize (M &m) {
88        typedef M matrix_type;
89        typedef typename M::size_type size_type;
90        typedef typename M::value_type value_type;
91
92#if BOOST_UBLAS_TYPE_CHECK
93        matrix_type cm (m);
94#endif
95        int singular = 0;
96        size_type size1 = m.size1 ();
97        size_type size2 = m.size2 ();
98        size_type size = (std::min) (size1, size2);
99        for (size_type i = 0; i < size; ++ i) {
100            matrix_column<M> mci (column (m, i));
101            matrix_row<M> mri (row (m, i));
102            if (m (i, i) != value_type/*zero*/()) {
103                project (mci, range (i + 1, size1)) *= value_type (1) / m (i, i);
104            } else if (singular == 0) {
105                singular = i + 1;
106            }
107            project (m, range (i + 1, size1), range (i + 1, size2)).minus_assign (
108                outer_prod (project (mci, range (i + 1, size1)),
109                            project (mri, range (i + 1, size2))));
110        }
111#if BOOST_UBLAS_TYPE_CHECK
112        BOOST_UBLAS_CHECK (singular != 0 ||
113                           detail::expression_type_check (prod (triangular_adaptor<matrix_type, unit_lower> (m),
114                                                                triangular_adaptor<matrix_type, upper> (m)),
115                                                          cm), internal_logic ());
116#endif
117        return singular;
118    }
119
120    // LU factorization with partial pivoting
121    template<class M, class PM>
122    typename M::size_type lu_factorize (M &m, PM &pm) {
123        typedef M matrix_type;
124        typedef typename M::size_type size_type;
125        typedef typename M::value_type value_type;
126
127#if BOOST_UBLAS_TYPE_CHECK
128        matrix_type cm (m);
129#endif
130        int singular = 0;
131        size_type size1 = m.size1 ();
132        size_type size2 = m.size2 ();
133        size_type size = (std::min) (size1, size2);
134        for (size_type i = 0; i < size; ++ i) {
135            matrix_column<M> mci (column (m, i));
136            matrix_row<M> mri (row (m, i));
137            size_type i_norm_inf = i + index_norm_inf (project (mci, range (i, size1)));
138            BOOST_UBLAS_CHECK (i_norm_inf < size1, external_logic ());
139            if (m (i_norm_inf, i) != value_type/*zero*/()) {
140                if (i_norm_inf != i) {
141                    pm (i) = i_norm_inf;
142                    row (m, i_norm_inf).swap (mri);
143                } else {
144                    BOOST_UBLAS_CHECK (pm (i) == i_norm_inf, external_logic ());
145                }
146                project (mci, range (i + 1, size1)) *= value_type (1) / m (i, i);
147            } else if (singular == 0) {
148                singular = i + 1;
149            }
150            project (m, range (i + 1, size1), range (i + 1, size2)).minus_assign (
151                outer_prod (project (mci, range (i + 1, size1)),
152                            project (mri, range (i + 1, size2))));
153        }
154#if BOOST_UBLAS_TYPE_CHECK
155        swap_rows (pm, cm);
156        BOOST_UBLAS_CHECK (singular != 0 ||
157                           detail::expression_type_check (prod (triangular_adaptor<matrix_type, unit_lower> (m),
158                                                                triangular_adaptor<matrix_type, upper> (m)), cm), internal_logic ());
159#endif
160        return singular;
161    }
162
163    template<class M, class PM>
164    typename M::size_type axpy_lu_factorize (M &m, PM &pm) {
165        typedef M matrix_type;
166        typedef typename M::size_type size_type;
167        typedef typename M::value_type value_type;
168        typedef vector<value_type> vector_type;
169
170#if BOOST_UBLAS_TYPE_CHECK
171        matrix_type cm (m);
172#endif
173        int singular = 0;
174        size_type size1 = m.size1 ();
175        size_type size2 = m.size2 ();
176        size_type size = (std::min) (size1, size2);
177#ifndef BOOST_UBLAS_LU_WITH_INPLACE_SOLVE
178        matrix_type mr (m);
179        mr.assign (zero_matrix<value_type> (size1, size2));
180        vector_type v (size1);
181        for (size_type i = 0; i < size; ++ i) {
182            matrix_range<matrix_type> lrr (project (mr, range (0, i), range (0, i)));
183            vector_range<matrix_column<matrix_type> > urr (project (column (mr, i), range (0, i)));
184            urr.assign (solve (lrr, project (column (m, i), range (0, i)), unit_lower_tag ()));
185            project (v, range (i, size1)).assign (
186                project (column (m, i), range (i, size1)) -
187                axpy_prod<vector_type> (project (mr, range (i, size1), range (0, i)), urr));
188            size_type i_norm_inf = i + index_norm_inf (project (v, range (i, size1)));
189            BOOST_UBLAS_CHECK (i_norm_inf < size1, external_logic ());
190            if (v (i_norm_inf) != value_type/*zero*/()) {
191                if (i_norm_inf != i) {
192                    pm (i) = i_norm_inf;
193                    std::swap (v (i_norm_inf), v (i));
194                    project (row (m, i_norm_inf), range (i + 1, size2)).swap (project (row (m, i), range (i + 1, size2)));
195                } else {
196                    BOOST_UBLAS_CHECK (pm (i) == i_norm_inf, external_logic ());
197                }
198                project (column (mr, i), range (i + 1, size1)).assign (
199                    project (v, range (i + 1, size1)) / v (i));
200                if (i_norm_inf != i) {
201                    project (row (mr, i_norm_inf), range (0, i)).swap (project (row (mr, i), range (0, i)));
202                }
203            } else if (singular == 0) {
204                singular = i + 1;
205            }
206            mr (i, i) = v (i);
207        }
208        m.assign (mr);
209#else
210        matrix_type lr (m);
211        matrix_type ur (m);
212        lr.assign (identity_matrix<value_type> (size1, size2));
213        ur.assign (zero_matrix<value_type> (size1, size2));
214        vector_type v (size1);
215        for (size_type i = 0; i < size; ++ i) {
216            matrix_range<matrix_type> lrr (project (lr, range (0, i), range (0, i)));
217            vector_range<matrix_column<matrix_type> > urr (project (column (ur, i), range (0, i)));
218            urr.assign (project (column (m, i), range (0, i)));
219            inplace_solve (lrr, urr, unit_lower_tag ());
220            project (v, range (i, size1)).assign (
221                project (column (m, i), range (i, size1)) -
222                axpy_prod<vector_type> (project (lr, range (i, size1), range (0, i)), urr));
223            size_type i_norm_inf = i + index_norm_inf (project (v, range (i, size1)));
224            BOOST_UBLAS_CHECK (i_norm_inf < size1, external_logic ());
225            if (v (i_norm_inf) != value_type/*zero*/()) {
226                if (i_norm_inf != i) {
227                    pm (i) = i_norm_inf;
228                    std::swap (v (i_norm_inf), v (i));
229                    project (row (m, i_norm_inf), range (i + 1, size2)).swap (project (row (m, i), range (i + 1, size2)));
230                } else {
231                    BOOST_UBLAS_CHECK (pm (i) == i_norm_inf, external_logic ());
232                }
233                project (column (lr, i), range (i + 1, size1)).assign (
234                    project (v, range (i + 1, size1)) / v (i));
235                if (i_norm_inf != i) {
236                    project (row (lr, i_norm_inf), range (0, i)).swap (project (row (lr, i), range (0, i)));
237                }
238            } else if (singular == 0) {
239                singular = i + 1;
240            }
241            ur (i, i) = v (i);
242        }
243        m.assign (triangular_adaptor<matrix_type, strict_lower> (lr) +
244                  triangular_adaptor<matrix_type, upper> (ur));
245#endif
246#if BOOST_UBLAS_TYPE_CHECK
247        swap_rows (pm, cm);
248        BOOST_UBLAS_CHECK (singular != 0 ||
249                           detail::expression_type_check (prod (triangular_adaptor<matrix_type, unit_lower> (m),
250                                                                triangular_adaptor<matrix_type, upper> (m)), cm), internal_logic ());
251#endif
252        return singular;
253    }
254
255    // LU substitution
256    template<class M, class E>
257    void lu_substitute (const M &m, vector_expression<E> &e) {
258        typedef const M const_matrix_type;
259        typedef vector<typename E::value_type> vector_type;
260
261#if BOOST_UBLAS_TYPE_CHECK
262        vector_type cv1 (e);
263#endif
264        inplace_solve (m, e, unit_lower_tag ());
265#if BOOST_UBLAS_TYPE_CHECK
266        BOOST_UBLAS_CHECK (detail::expression_type_check (prod (triangular_adaptor<const_matrix_type, unit_lower> (m), e), cv1), internal_logic ());
267        vector_type cv2 (e);
268#endif
269        inplace_solve (m, e, upper_tag ());
270#if BOOST_UBLAS_TYPE_CHECK
271        BOOST_UBLAS_CHECK (detail::expression_type_check (prod (triangular_adaptor<const_matrix_type, upper> (m), e), cv2), internal_logic ());
272#endif
273    }
274    template<class M, class E>
275    void lu_substitute (const M &m, matrix_expression<E> &e) {
276        typedef const M const_matrix_type;
277        typedef matrix<typename E::value_type> matrix_type;
278
279#if BOOST_UBLAS_TYPE_CHECK
280        matrix_type cm1 (e);
281#endif
282        inplace_solve (m, e, unit_lower_tag ());
283#if BOOST_UBLAS_TYPE_CHECK
284        BOOST_UBLAS_CHECK (detail::expression_type_check (prod (triangular_adaptor<const_matrix_type, unit_lower> (m), e), cm1), internal_logic ());
285        matrix_type cm2 (e);
286#endif
287        inplace_solve (m, e, upper_tag ());
288#if BOOST_UBLAS_TYPE_CHECK
289        BOOST_UBLAS_CHECK (detail::expression_type_check (prod (triangular_adaptor<const_matrix_type, upper> (m), e), cm2), internal_logic ());
290#endif
291    }
292    template<class M, class PMT, class PMA, class MV>
293    void lu_substitute (const M &m, const permutation_matrix<PMT, PMA> &pm, MV &mv) {
294        swap_rows (pm, mv);
295        lu_substitute (m, mv);
296    }
297    template<class E, class M>
298    void lu_substitute (vector_expression<E> &e, const M &m) {
299        typedef const M const_matrix_type;
300        typedef vector<typename E::value_type> vector_type;
301
302#if BOOST_UBLAS_TYPE_CHECK
303        vector_type cv1 (e);
304#endif
305        inplace_solve (e, m, upper_tag ());
306#if BOOST_UBLAS_TYPE_CHECK
307        BOOST_UBLAS_CHECK (detail::expression_type_check (prod (e, triangular_adaptor<const_matrix_type, upper> (m)), cv1), internal_logic ());
308        vector_type cv2 (e);
309#endif
310        inplace_solve (e, m, unit_lower_tag ());
311#if BOOST_UBLAS_TYPE_CHECK
312        BOOST_UBLAS_CHECK (detail::expression_type_check (prod (e, triangular_adaptor<const_matrix_type, unit_lower> (m)), cv2), internal_logic ());
313#endif
314    }
315    template<class E, class M>
316    void lu_substitute (matrix_expression<E> &e, const M &m) {
317        typedef const M const_matrix_type;
318        typedef matrix<typename E::value_type> matrix_type;
319
320#if BOOST_UBLAS_TYPE_CHECK
321        matrix_type cm1 (e);
322#endif
323        inplace_solve (e, m, upper_tag ());
324#if BOOST_UBLAS_TYPE_CHECK
325        BOOST_UBLAS_CHECK (detail::expression_type_check (prod (e, triangular_adaptor<const_matrix_type, upper> (m)), cm1), internal_logic ());
326        matrix_type cm2 (e);
327#endif
328        inplace_solve (e, m, unit_lower_tag ());
329#if BOOST_UBLAS_TYPE_CHECK
330        BOOST_UBLAS_CHECK (detail::expression_type_check (prod (e, triangular_adaptor<const_matrix_type, unit_lower> (m)), cm2), internal_logic ());
331#endif
332    }
333    template<class MV, class M, class PMT, class PMA>
334    void lu_substitute (MV &mv, const M &m, const permutation_matrix<PMT, PMA> &pm) {
335        swap_rows (pm, mv);
336        lu_substitute (mv, m);
337    }
338
339}}}
340
341#endif
Note: See TracBrowser for help on using the repository browser.