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

Revision 857, 13.9 KB checked in by igarcia, 19 years ago (diff)
Line 
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_OPERATION_BLOCKED_
18#define _BOOST_UBLAS_OPERATION_BLOCKED_
19
20#include <boost/numeric/ublas/traits.hpp>
21
22
23namespace boost { namespace numeric { namespace ublas {
24
25    template<class V, typename V::size_type BS, class E1, class E2>
26    BOOST_UBLAS_INLINE
27    V
28    block_prod (const matrix_expression<E1> &e1,
29                const vector_expression<E2> &e2) {
30        typedef V vector_type;
31        typedef const E1 expression1_type;
32        typedef const E2 expression2_type;
33        typedef typename V::size_type size_type;
34        typedef typename V::value_type value_type;
35        const size_type block_size = BS;
36
37        V v (e1 ().size1 ());
38#if BOOST_UBLAS_TYPE_CHECK
39        vector<value_type> cv (v.size ());
40        typedef typename type_traits<value_type>::real_type real_type;
41        real_type verrorbound (norm_1 (v) + norm_1 (e1) * norm_1 (e2));
42        indexing_vector_assign<scalar_assign> (cv, prod (e1, e2));
43#endif
44        size_type i_size = e1 ().size1 ();
45        size_type j_size = BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size ());
46        for (size_type i_begin = 0; i_begin < i_size; i_begin += block_size) {
47            size_type i_end = i_begin + (std::min) (i_size - i_begin, block_size);
48            // FIX: never ignore Martin Weiser's advice ;-(
49#ifdef BOOST_UBLAS_NO_CACHE
50            vector_range<vector_type> v_range (v, range (i_begin, i_end));
51#else
52            // vector<value_type, bounded_array<value_type, block_size> > v_range (i_end - i_begin);
53            vector<value_type> v_range (i_end - i_begin);
54#endif
55            v_range.assign (zero_vector<value_type> (i_end - i_begin));
56            for (size_type j_begin = 0; j_begin < j_size; j_begin += block_size) {
57                size_type j_end = j_begin + (std::min) (j_size - j_begin, block_size);
58#ifdef BOOST_UBLAS_NO_CACHE
59                const matrix_range<expression1_type> e1_range (e1 (), range (i_begin, i_end), range (j_begin, j_end));
60                const vector_range<expression2_type> e2_range (e2 (), range (j_begin, j_end));
61                v_range.plus_assign (prod (e1_range, e2_range));
62#else
63                // const matrix<value_type, row_major, bounded_array<value_type, block_size * block_size> > e1_range (project (e1 (), range (i_begin, i_end), range (j_begin, j_end)));
64                // const vector<value_type, bounded_array<value_type, block_size> > e2_range (project (e2 (), range (j_begin, j_end)));
65                const matrix<value_type, row_major> e1_range (project (e1 (), range (i_begin, i_end), range (j_begin, j_end)));
66                const vector<value_type> e2_range (project (e2 (), range (j_begin, j_end)));
67                v_range.plus_assign (prod (e1_range, e2_range));
68#endif
69            }
70#ifndef BOOST_UBLAS_NO_CACHE
71            project (v, range (i_begin, i_end)).assign (v_range);
72#endif
73        }
74#if BOOST_UBLAS_TYPE_CHECK
75        BOOST_UBLAS_CHECK (norm_1 (v - cv) <= 2 * std::numeric_limits<real_type>::epsilon () * verrorbound, internal_logic ());
76#endif
77        return v;
78    }
79
80    template<class V, typename V::size_type BS, class E1, class E2>
81    BOOST_UBLAS_INLINE
82    V
83    block_prod (const vector_expression<E1> &e1,
84                const matrix_expression<E2> &e2) {
85        typedef V vector_type;
86        typedef const E1 expression1_type;
87        typedef const E2 expression2_type;
88        typedef typename V::size_type size_type;
89        typedef typename V::value_type value_type;
90        const size_type block_size = BS;
91
92        V v (e2 ().size2 ());
93#if BOOST_UBLAS_TYPE_CHECK
94        vector<value_type> cv (v.size ());
95        typedef typename type_traits<value_type>::real_type real_type;
96        real_type verrorbound (norm_1 (v) + norm_1 (e1) * norm_1 (e2));
97        indexing_vector_assign<scalar_assign> (cv, prod (e1, e2));
98#endif
99        size_type i_size = BOOST_UBLAS_SAME (e1 ().size (), e2 ().size1 ());
100        size_type j_size = e2 ().size2 ();
101        for (size_type j_begin = 0; j_begin < j_size; j_begin += block_size) {
102            size_type j_end = j_begin + (std::min) (j_size - j_begin, block_size);
103            // FIX: never ignore Martin Weiser's advice ;-(
104#ifdef BOOST_UBLAS_NO_CACHE
105            vector_range<vector_type> v_range (v, range (j_begin, j_end));
106#else
107            // vector<value_type, bounded_array<value_type, block_size> > v_range (j_end - j_begin);
108            vector<value_type> v_range (j_end - j_begin);
109#endif
110            v_range.assign (zero_vector<value_type> (j_end - j_begin));
111            for (size_type i_begin = 0; i_begin < i_size; i_begin += block_size) {
112                size_type i_end = i_begin + (std::min) (i_size - i_begin, block_size);
113#ifdef BOOST_UBLAS_NO_CACHE
114                const vector_range<expression1_type> e1_range (e1 (), range (i_begin, i_end));
115                const matrix_range<expression2_type> e2_range (e2 (), range (i_begin, i_end), range (j_begin, j_end));
116#else
117                // const vector<value_type, bounded_array<value_type, block_size> > e1_range (project (e1 (), range (i_begin, i_end)));
118                // const matrix<value_type, column_major, bounded_array<value_type, block_size * block_size> > e2_range (project (e2 (), range (i_begin, i_end), range (j_begin, j_end)));
119                const vector<value_type> e1_range (project (e1 (), range (i_begin, i_end)));
120                const matrix<value_type, column_major> e2_range (project (e2 (), range (i_begin, i_end), range (j_begin, j_end)));
121#endif
122                v_range.plus_assign (prod (e1_range, e2_range));
123            }
124#ifndef BOOST_UBLAS_NO_CACHE
125            project (v, range (j_begin, j_end)).assign (v_range);
126#endif
127        }
128#if BOOST_UBLAS_TYPE_CHECK
129        BOOST_UBLAS_CHECK (norm_1 (v - cv) <= 2 * std::numeric_limits<real_type>::epsilon () * verrorbound, internal_logic ());
130#endif
131        return v;
132    }
133
134    template<class M, typename M::size_type BS, class E1, class E2>
135    BOOST_UBLAS_INLINE
136    M
137    block_prod (const matrix_expression<E1> &e1,
138                const matrix_expression<E2> &e2,
139                row_major_tag) {
140        typedef M matrix_type;
141        typedef const E1 expression1_type;
142        typedef const E2 expression2_type;
143        typedef typename M::size_type size_type;
144        typedef typename M::value_type value_type;
145        const size_type block_size = BS;
146
147        M m (e1 ().size1 (), e2 ().size2 ());
148#if BOOST_UBLAS_TYPE_CHECK
149        matrix<value_type, row_major> cm (m.size1 (), m.size2 ());
150        typedef typename type_traits<value_type>::real_type real_type;
151        real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));
152        indexing_matrix_assign<scalar_assign> (cm, prod (e1, e2), row_major_tag ());
153        disable_type_check<bool>::value = true;
154#endif
155        size_type i_size = e1 ().size1 ();
156        size_type j_size = e2 ().size2 ();
157        size_type k_size = BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size1 ());
158        for (size_type i_begin = 0; i_begin < i_size; i_begin += block_size) {
159            size_type i_end = i_begin + (std::min) (i_size - i_begin, block_size);
160            for (size_type j_begin = 0; j_begin < j_size; j_begin += block_size) {
161                size_type j_end = j_begin + (std::min) (j_size - j_begin, block_size);
162                // FIX: never ignore Martin Weiser's advice ;-(
163#ifdef BOOST_UBLAS_NO_CACHE
164                matrix_range<matrix_type> m_range (m, range (i_begin, i_end), range (j_begin, j_end));
165#else
166                // matrix<value_type, row_major, bounded_array<value_type, block_size * block_size> > m_range (i_end - i_begin, j_end - j_begin);
167                matrix<value_type, row_major> m_range (i_end - i_begin, j_end - j_begin);
168#endif
169                m_range.assign (zero_matrix<value_type> (i_end - i_begin, j_end - j_begin));
170                for (size_type k_begin = 0; k_begin < k_size; k_begin += block_size) {
171                    size_type k_end = k_begin + (std::min) (k_size - k_begin, block_size);
172#ifdef BOOST_UBLAS_NO_CACHE
173                    const matrix_range<expression1_type> e1_range (e1 (), range (i_begin, i_end), range (k_begin, k_end));
174                    const matrix_range<expression2_type> e2_range (e2 (), range (k_begin, k_end), range (j_begin, j_end));
175#else
176                    // const matrix<value_type, row_major, bounded_array<value_type, block_size * block_size> > e1_range (project (e1 (), range (i_begin, i_end), range (k_begin, k_end)));
177                    // const matrix<value_type, column_major, bounded_array<value_type, block_size * block_size> > e2_range (project (e2 (), range (k_begin, k_end), range (j_begin, j_end)));
178                    const matrix<value_type, row_major> e1_range (project (e1 (), range (i_begin, i_end), range (k_begin, k_end)));
179                    const matrix<value_type, column_major> e2_range (project (e2 (), range (k_begin, k_end), range (j_begin, j_end)));
180#endif
181                    m_range.plus_assign (prod (e1_range, e2_range));
182                }
183#ifndef BOOST_UBLAS_NO_CACHE
184                project (m, range (i_begin, i_end), range (j_begin, j_end)).assign (m_range);
185#endif
186            }
187        }
188#if BOOST_UBLAS_TYPE_CHECK
189        disable_type_check<bool>::value = false;
190        BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());
191#endif
192        return m;
193    }
194
195    template<class M, typename M::size_type BS, class E1, class E2>
196    BOOST_UBLAS_INLINE
197    M
198    block_prod (const matrix_expression<E1> &e1,
199                const matrix_expression<E2> &e2,
200                column_major_tag) {
201        typedef M matrix_type;
202        typedef const E1 expression1_type;
203        typedef const E2 expression2_type;
204        typedef typename M::size_type size_type;
205        typedef typename M::value_type value_type;
206        const size_type block_size = BS;
207
208        M m (e1 ().size1 (), e2 ().size2 ());
209#if BOOST_UBLAS_TYPE_CHECK
210        matrix<value_type, column_major> cm (m.size1 (), m.size2 ());
211        typedef typename type_traits<value_type>::real_type real_type;
212        real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));
213        indexing_matrix_assign<scalar_assign> (cm, prod (e1, e2), column_major_tag ());
214        disable_type_check<bool>::value = true;
215#endif
216        size_type i_size = e1 ().size1 ();
217        size_type j_size = e2 ().size2 ();
218        size_type k_size = BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size1 ());
219        for (size_type j_begin = 0; j_begin < j_size; j_begin += block_size) {
220            size_type j_end = j_begin + (std::min) (j_size - j_begin, block_size);
221            for (size_type i_begin = 0; i_begin < i_size; i_begin += block_size) {
222                size_type i_end = i_begin + (std::min) (i_size - i_begin, block_size);
223                // FIX: never ignore Martin Weiser's advice ;-(
224#ifdef BOOST_UBLAS_NO_CACHE
225                matrix_range<matrix_type> m_range (m, range (i_begin, i_end), range (j_begin, j_end));
226#else
227                // matrix<value_type, column_major, bounded_array<value_type, block_size * block_size> > m_range (i_end - i_begin, j_end - j_begin);
228                matrix<value_type, column_major> m_range (i_end - i_begin, j_end - j_begin);
229#endif
230                m_range.assign (zero_matrix<value_type> (i_end - i_begin, j_end - j_begin));
231                for (size_type k_begin = 0; k_begin < k_size; k_begin += block_size) {
232                    size_type k_end = k_begin + (std::min) (k_size - k_begin, block_size);
233#ifdef BOOST_UBLAS_NO_CACHE
234                    const matrix_range<expression1_type> e1_range (e1 (), range (i_begin, i_end), range (k_begin, k_end));
235                    const matrix_range<expression2_type> e2_range (e2 (), range (k_begin, k_end), range (j_begin, j_end));
236#else
237                    // const matrix<value_type, row_major, bounded_array<value_type, block_size * block_size> > e1_range (project (e1 (), range (i_begin, i_end), range (k_begin, k_end)));
238                    // const matrix<value_type, column_major, bounded_array<value_type, block_size * block_size> > e2_range (project (e2 (), range (k_begin, k_end), range (j_begin, j_end)));
239                    const matrix<value_type, row_major> e1_range (project (e1 (), range (i_begin, i_end), range (k_begin, k_end)));
240                    const matrix<value_type, column_major> e2_range (project (e2 (), range (k_begin, k_end), range (j_begin, j_end)));
241#endif
242                    m_range.plus_assign (prod (e1_range, e2_range));
243                }
244#ifndef BOOST_UBLAS_NO_CACHE
245                project (m, range (i_begin, i_end), range (j_begin, j_end)).assign (m_range);
246#endif
247            }
248        }
249#if BOOST_UBLAS_TYPE_CHECK
250        disable_type_check<bool>::value = false;
251        BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());
252#endif
253        return m;
254    }
255
256    // Dispatcher
257    template<class M, typename M::size_type BS, class E1, class E2>
258    BOOST_UBLAS_INLINE
259    M
260    block_prod (const matrix_expression<E1> &e1,
261                const matrix_expression<E2> &e2) {
262        typedef typename M::orientation_category orientation_category;
263        return block_prod<M, BS> (e1, e2, orientation_category ());
264    }
265
266}}}
267
268#endif
Note: See TracBrowser for help on using the repository browser.