source: vendor/nemo/current/NEMOGCM/EXTERNAL/XIOS/extern/boost/include/boost/numeric/ublas/lu.hpp @ 44

Last change on this file since 44 was 44, checked in by cholod, 12 years ago

Load NEMO_TMP into vendor/nemo/current.

File size: 13.8 KB
Line 
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        typedef typename MV::value_type value_type;
67
68        size_type size = pm.size ();
69        for (size_type i = 0; i < size; ++ i) {
70            if (i != pm (i))
71                std::swap (mv (i), mv (pm (i)));
72        }
73    }
74    template<class PM, class MV>
75    BOOST_UBLAS_INLINE
76    void swap_rows (const PM &pm, MV &mv, matrix_tag) {
77        typedef typename PM::size_type size_type;
78        typedef typename MV::value_type value_type;
79
80        size_type size = pm.size ();
81        for (size_type i = 0; i < size; ++ i) {
82            if (i != pm (i))
83                row (mv, i).swap (row (mv, pm (i)));
84        }
85    }
86    // Dispatcher
87    template<class PM, class MV>
88    BOOST_UBLAS_INLINE
89    void swap_rows (const PM &pm, MV &mv) {
90        swap_rows (pm, mv, typename MV::type_category ());
91    }
92
93    // LU factorization without pivoting
94    template<class M>
95    typename M::size_type lu_factorize (M &m) {
96        typedef M matrix_type;
97        typedef typename M::size_type size_type;
98        typedef typename M::value_type value_type;
99
100#if BOOST_UBLAS_TYPE_CHECK
101        matrix_type cm (m);
102#endif
103        size_type singular = 0;
104        size_type size1 = m.size1 ();
105        size_type size2 = m.size2 ();
106        size_type size = (std::min) (size1, size2);
107        for (size_type i = 0; i < size; ++ i) {
108            matrix_column<M> mci (column (m, i));
109            matrix_row<M> mri (row (m, i));
110            if (m (i, i) != value_type/*zero*/()) {
111                value_type m_inv = value_type (1) / m (i, i);
112                project (mci, range (i + 1, size1)) *= m_inv;
113            } else if (singular == 0) {
114                singular = i + 1;
115            }
116            project (m, range (i + 1, size1), range (i + 1, size2)).minus_assign (
117                outer_prod (project (mci, range (i + 1, size1)),
118                            project (mri, range (i + 1, size2))));
119        }
120#if BOOST_UBLAS_TYPE_CHECK
121        BOOST_UBLAS_CHECK (singular != 0 ||
122                           detail::expression_type_check (prod (triangular_adaptor<matrix_type, unit_lower> (m),
123                                                                triangular_adaptor<matrix_type, upper> (m)), 
124                                                          cm), internal_logic ());
125#endif
126        return singular;
127    }
128
129    // LU factorization with partial pivoting
130    template<class M, class PM>
131    typename M::size_type lu_factorize (M &m, PM &pm) {
132        typedef M matrix_type;
133        typedef typename M::size_type size_type;
134        typedef typename M::value_type value_type;
135
136#if BOOST_UBLAS_TYPE_CHECK
137        matrix_type cm (m);
138#endif
139        size_type singular = 0;
140        size_type size1 = m.size1 ();
141        size_type size2 = m.size2 ();
142        size_type size = (std::min) (size1, size2);
143        for (size_type i = 0; i < size; ++ i) {
144            matrix_column<M> mci (column (m, i));
145            matrix_row<M> mri (row (m, i));
146            size_type i_norm_inf = i + index_norm_inf (project (mci, range (i, size1)));
147            BOOST_UBLAS_CHECK (i_norm_inf < size1, external_logic ());
148            if (m (i_norm_inf, i) != value_type/*zero*/()) {
149                if (i_norm_inf != i) {
150                    pm (i) = i_norm_inf;
151                    row (m, i_norm_inf).swap (mri);
152                } else {
153                    BOOST_UBLAS_CHECK (pm (i) == i_norm_inf, external_logic ());
154                }
155                value_type m_inv = value_type (1) / m (i, i);
156                project (mci, range (i + 1, size1)) *= m_inv;
157            } else if (singular == 0) {
158                singular = i + 1;
159            }
160            project (m, range (i + 1, size1), range (i + 1, size2)).minus_assign (
161                outer_prod (project (mci, range (i + 1, size1)),
162                            project (mri, range (i + 1, size2))));
163        }
164#if BOOST_UBLAS_TYPE_CHECK
165        swap_rows (pm, cm);
166        BOOST_UBLAS_CHECK (singular != 0 ||
167                           detail::expression_type_check (prod (triangular_adaptor<matrix_type, unit_lower> (m),
168                                                                triangular_adaptor<matrix_type, upper> (m)), cm), internal_logic ());
169#endif
170        return singular;
171    }
172
173    template<class M, class PM>
174    typename M::size_type axpy_lu_factorize (M &m, PM &pm) {
175        typedef M matrix_type;
176        typedef typename M::size_type size_type;
177        typedef typename M::value_type value_type;
178        typedef vector<value_type> vector_type;
179
180#if BOOST_UBLAS_TYPE_CHECK
181        matrix_type cm (m);
182#endif
183        size_type singular = 0;
184        size_type size1 = m.size1 ();
185        size_type size2 = m.size2 ();
186        size_type size = (std::min) (size1, size2);
187#ifndef BOOST_UBLAS_LU_WITH_INPLACE_SOLVE
188        matrix_type mr (m);
189        mr.assign (zero_matrix<value_type> (size1, size2));
190        vector_type v (size1);
191        for (size_type i = 0; i < size; ++ i) {
192            matrix_range<matrix_type> lrr (project (mr, range (0, i), range (0, i)));
193            vector_range<matrix_column<matrix_type> > urr (project (column (mr, i), range (0, i)));
194            urr.assign (solve (lrr, project (column (m, i), range (0, i)), unit_lower_tag ()));
195            project (v, range (i, size1)).assign (
196                project (column (m, i), range (i, size1)) -
197                axpy_prod<vector_type> (project (mr, range (i, size1), range (0, i)), urr));
198            size_type i_norm_inf = i + index_norm_inf (project (v, range (i, size1)));
199            BOOST_UBLAS_CHECK (i_norm_inf < size1, external_logic ());
200            if (v (i_norm_inf) != value_type/*zero*/()) {
201                if (i_norm_inf != i) {
202                    pm (i) = i_norm_inf;
203                    std::swap (v (i_norm_inf), v (i));
204                    project (row (m, i_norm_inf), range (i + 1, size2)).swap (project (row (m, i), range (i + 1, size2)));
205                } else {
206                    BOOST_UBLAS_CHECK (pm (i) == i_norm_inf, external_logic ());
207                }
208                project (column (mr, i), range (i + 1, size1)).assign (
209                    project (v, range (i + 1, size1)) / v (i));
210                if (i_norm_inf != i) {
211                    project (row (mr, i_norm_inf), range (0, i)).swap (project (row (mr, i), range (0, i)));
212                }
213            } else if (singular == 0) {
214                singular = i + 1;
215            }
216            mr (i, i) = v (i);
217        }
218        m.assign (mr);
219#else
220        matrix_type lr (m);
221        matrix_type ur (m);
222        lr.assign (identity_matrix<value_type> (size1, size2));
223        ur.assign (zero_matrix<value_type> (size1, size2));
224        vector_type v (size1);
225        for (size_type i = 0; i < size; ++ i) {
226            matrix_range<matrix_type> lrr (project (lr, range (0, i), range (0, i)));
227            vector_range<matrix_column<matrix_type> > urr (project (column (ur, i), range (0, i)));
228            urr.assign (project (column (m, i), range (0, i)));
229            inplace_solve (lrr, urr, unit_lower_tag ());
230            project (v, range (i, size1)).assign (
231                project (column (m, i), range (i, size1)) -
232                axpy_prod<vector_type> (project (lr, range (i, size1), range (0, i)), urr));
233            size_type i_norm_inf = i + index_norm_inf (project (v, range (i, size1)));
234            BOOST_UBLAS_CHECK (i_norm_inf < size1, external_logic ());
235            if (v (i_norm_inf) != value_type/*zero*/()) {
236                if (i_norm_inf != i) {
237                    pm (i) = i_norm_inf;
238                    std::swap (v (i_norm_inf), v (i));
239                    project (row (m, i_norm_inf), range (i + 1, size2)).swap (project (row (m, i), range (i + 1, size2)));
240                } else {
241                    BOOST_UBLAS_CHECK (pm (i) == i_norm_inf, external_logic ());
242                }
243                project (column (lr, i), range (i + 1, size1)).assign (
244                    project (v, range (i + 1, size1)) / v (i));
245                if (i_norm_inf != i) {
246                    project (row (lr, i_norm_inf), range (0, i)).swap (project (row (lr, i), range (0, i)));
247                }
248            } else if (singular == 0) {
249                singular = i + 1;
250            }
251            ur (i, i) = v (i);
252        }
253        m.assign (triangular_adaptor<matrix_type, strict_lower> (lr) +
254                  triangular_adaptor<matrix_type, upper> (ur));
255#endif
256#if BOOST_UBLAS_TYPE_CHECK
257        swap_rows (pm, cm);
258        BOOST_UBLAS_CHECK (singular != 0 ||
259                           detail::expression_type_check (prod (triangular_adaptor<matrix_type, unit_lower> (m),
260                                                                triangular_adaptor<matrix_type, upper> (m)), cm), internal_logic ());
261#endif
262        return singular;
263    }
264
265    // LU substitution
266    template<class M, class E>
267    void lu_substitute (const M &m, vector_expression<E> &e) {
268        typedef const M const_matrix_type;
269        typedef vector<typename E::value_type> vector_type;
270
271#if BOOST_UBLAS_TYPE_CHECK
272        vector_type cv1 (e);
273#endif
274        inplace_solve (m, e, unit_lower_tag ());
275#if BOOST_UBLAS_TYPE_CHECK
276        BOOST_UBLAS_CHECK (detail::expression_type_check (prod (triangular_adaptor<const_matrix_type, unit_lower> (m), e), cv1), internal_logic ());
277        vector_type cv2 (e);
278#endif
279        inplace_solve (m, e, upper_tag ());
280#if BOOST_UBLAS_TYPE_CHECK
281        BOOST_UBLAS_CHECK (detail::expression_type_check (prod (triangular_adaptor<const_matrix_type, upper> (m), e), cv2), internal_logic ());
282#endif
283    }
284    template<class M, class E>
285    void lu_substitute (const M &m, matrix_expression<E> &e) {
286        typedef const M const_matrix_type;
287        typedef matrix<typename E::value_type> matrix_type;
288
289#if BOOST_UBLAS_TYPE_CHECK
290        matrix_type cm1 (e);
291#endif
292        inplace_solve (m, e, unit_lower_tag ());
293#if BOOST_UBLAS_TYPE_CHECK
294        BOOST_UBLAS_CHECK (detail::expression_type_check (prod (triangular_adaptor<const_matrix_type, unit_lower> (m), e), cm1), internal_logic ());
295        matrix_type cm2 (e);
296#endif
297        inplace_solve (m, e, upper_tag ());
298#if BOOST_UBLAS_TYPE_CHECK
299        BOOST_UBLAS_CHECK (detail::expression_type_check (prod (triangular_adaptor<const_matrix_type, upper> (m), e), cm2), internal_logic ());
300#endif
301    }
302    template<class M, class PMT, class PMA, class MV>
303    void lu_substitute (const M &m, const permutation_matrix<PMT, PMA> &pm, MV &mv) {
304        swap_rows (pm, mv);
305        lu_substitute (m, mv);
306    }
307    template<class E, class M>
308    void lu_substitute (vector_expression<E> &e, const M &m) {
309        typedef const M const_matrix_type;
310        typedef vector<typename E::value_type> vector_type;
311
312#if BOOST_UBLAS_TYPE_CHECK
313        vector_type cv1 (e);
314#endif
315        inplace_solve (e, m, upper_tag ());
316#if BOOST_UBLAS_TYPE_CHECK
317        BOOST_UBLAS_CHECK (detail::expression_type_check (prod (e, triangular_adaptor<const_matrix_type, upper> (m)), cv1), internal_logic ());
318        vector_type cv2 (e);
319#endif
320        inplace_solve (e, m, unit_lower_tag ());
321#if BOOST_UBLAS_TYPE_CHECK
322        BOOST_UBLAS_CHECK (detail::expression_type_check (prod (e, triangular_adaptor<const_matrix_type, unit_lower> (m)), cv2), internal_logic ());
323#endif
324    }
325    template<class E, class M>
326    void lu_substitute (matrix_expression<E> &e, const M &m) {
327        typedef const M const_matrix_type;
328        typedef matrix<typename E::value_type> matrix_type;
329
330#if BOOST_UBLAS_TYPE_CHECK
331        matrix_type cm1 (e);
332#endif
333        inplace_solve (e, m, upper_tag ());
334#if BOOST_UBLAS_TYPE_CHECK
335        BOOST_UBLAS_CHECK (detail::expression_type_check (prod (e, triangular_adaptor<const_matrix_type, upper> (m)), cm1), internal_logic ());
336        matrix_type cm2 (e);
337#endif
338        inplace_solve (e, m, unit_lower_tag ());
339#if BOOST_UBLAS_TYPE_CHECK
340        BOOST_UBLAS_CHECK (detail::expression_type_check (prod (e, triangular_adaptor<const_matrix_type, unit_lower> (m)), cm2), internal_logic ());
341#endif
342    }
343    template<class MV, class M, class PMT, class PMA>
344    void lu_substitute (MV &mv, const M &m, const permutation_matrix<PMT, PMA> &pm) {
345        swap_rows (pm, mv);
346        lu_substitute (mv, m);
347    }
348
349}}}
350
351#endif
Note: See TracBrowser for help on using the repository browser.