numeric/ublas/operation_blocked.hpp

0001 //
0002 //  Copyright (c) 2000-2002
0003 //  Joerg Walter, Mathias Koch
0004 //
0005 //  Distributed under the Boost Software License, Version 1.0. (See
0006 //  accompanying file LICENSE_1_0.txt or copy at
0007 //  http://www.boost.org/LICENSE_1_0.txt)
0008 //
0009 //  The authors gratefully acknowledge the support of
0010 //  GeNeSys mbH & Co. KG in producing this work.
0011 //
0012
0013 #ifndef _BOOST_UBLAS_OPERATION_BLOCKED_
0014 #define _BOOST_UBLAS_OPERATION_BLOCKED_
0015
0016 #include <boost/numeric/ublas/traits.hpp>
0017 #include <boost/numeric/ublas/detail/vector_assign.hpp> // indexing_vector_assign
0018 #include <boost/numeric/ublas/detail/matrix_assign.hpp> // indexing_matrix_assign
0019
0020
0021 namespace boost { namespace numeric { namespace ublas {
0022
0023     template<class V, typename V::size_type BS, class E1, class E2>
0024     BOOST_UBLAS_INLINE
0025     V
0026     block_prod (const matrix_expression<E1> &e1,
0027                 const vector_expression<E2> &e2) {
0028         typedef V vector_type;
0029         typedef const E1 expression1_type;
0030         typedef const E2 expression2_type;
0031         typedef typename V::size_type size_type;
0032         typedef typename V::value_type value_type;
0033         const size_type block_size = BS;
0034
0035         V v (e1 ().size1 ());
0036 #if BOOST_UBLAS_TYPE_CHECK
0037         vector<value_type> cv (v.size ());
0038         typedef typename type_traits<value_type>::real_type real_type;
0039         real_type verrorbound (norm_1 (v) + norm_1 (e1) * norm_1 (e2));
0040         indexing_vector_assign<scalar_assign> (cv, prod (e1, e2));
0041 #endif
0042         size_type i_size = e1 ().size1 ();
0043         size_type j_size = BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size ());
0044         for (size_type i_begin = 0; i_begin < i_size; i_begin += block_size) {
0045             size_type i_end = i_begin + (std::min) (i_size - i_begin, block_size);
0046             // FIX: never ignore Martin Weiser's advice ;-(
0047 #ifdef BOOST_UBLAS_NO_CACHE
0048             vector_range<vector_type> v_range (v, range (i_begin, i_end));
0049 #else
0050             // vector<value_type, bounded_array<value_type, block_size> > v_range (i_end - i_begin);
0051             vector<value_type> v_range (i_end - i_begin);
0052 #endif
0053             v_range.assign (zero_vector<value_type> (i_end - i_begin));
0054             for (size_type j_begin = 0; j_begin < j_size; j_begin += block_size) {
0055                 size_type j_end = j_begin + (std::min) (j_size - j_begin, block_size);
0056 #ifdef BOOST_UBLAS_NO_CACHE
0057                 const matrix_range<expression1_type> e1_range (e1 (), range (i_begin, i_end), range (j_begin, j_end));
0058                 const vector_range<expression2_type> e2_range (e2 (), range (j_begin, j_end));
0059                 v_range.plus_assign (prod (e1_range, e2_range));
0060 #else
0061                 // 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)));
0062                 // const vector<value_type, bounded_array<value_type, block_size> > e2_range (project (e2 (), range (j_begin, j_end)));
0063                 const matrix<value_type, row_major> e1_range (project (e1 (), range (i_begin, i_end), range (j_begin, j_end)));
0064                 const vector<value_type> e2_range (project (e2 (), range (j_begin, j_end)));
0065                 v_range.plus_assign (prod (e1_range, e2_range));
0066 #endif
0067             }
0068 #ifndef BOOST_UBLAS_NO_CACHE
0069             project (v, range (i_begin, i_end)).assign (v_range);
0070 #endif
0071         }
0072 #if BOOST_UBLAS_TYPE_CHECK
0073         BOOST_UBLAS_CHECK (norm_1 (v - cv) <= 2 * std::numeric_limits<real_type>::epsilon () * verrorbound, internal_logic ());
0074 #endif
0075         return v;
0076     }
0077
0078     template<class V, typename V::size_type BS, class E1, class E2>
0079     BOOST_UBLAS_INLINE
0080     V
0081     block_prod (const vector_expression<E1> &e1,
0082                 const matrix_expression<E2> &e2) {
0083         typedef V vector_type;
0084         typedef const E1 expression1_type;
0085         typedef const E2 expression2_type;
0086         typedef typename V::size_type size_type;
0087         typedef typename V::value_type value_type;
0088         const size_type block_size = BS;
0089
0090         V v (e2 ().size2 ());
0091 #if BOOST_UBLAS_TYPE_CHECK
0092         vector<value_type> cv (v.size ());
0093         typedef typename type_traits<value_type>::real_type real_type;
0094         real_type verrorbound (norm_1 (v) + norm_1 (e1) * norm_1 (e2));
0095         indexing_vector_assign<scalar_assign> (cv, prod (e1, e2));
0096 #endif
0097         size_type i_size = BOOST_UBLAS_SAME (e1 ().size (), e2 ().size1 ());
0098         size_type j_size = e2 ().size2 ();
0099         for (size_type j_begin = 0; j_begin < j_size; j_begin += block_size) {
0100             size_type j_end = j_begin + (std::min) (j_size - j_begin, block_size);
0101             // FIX: never ignore Martin Weiser's advice ;-(
0102 #ifdef BOOST_UBLAS_NO_CACHE
0103             vector_range<vector_type> v_range (v, range (j_begin, j_end));
0104 #else
0105             // vector<value_type, bounded_array<value_type, block_size> > v_range (j_end - j_begin);
0106             vector<value_type> v_range (j_end - j_begin);
0107 #endif
0108             v_range.assign (zero_vector<value_type> (j_end - j_begin));
0109             for (size_type i_begin = 0; i_begin < i_size; i_begin += block_size) {
0110                 size_type i_end = i_begin + (std::min) (i_size - i_begin, block_size);
0111 #ifdef BOOST_UBLAS_NO_CACHE
0112                 const vector_range<expression1_type> e1_range (e1 (), range (i_begin, i_end));
0113                 const matrix_range<expression2_type> e2_range (e2 (), range (i_begin, i_end), range (j_begin, j_end));
0114 #else
0115                 // const vector<value_type, bounded_array<value_type, block_size> > e1_range (project (e1 (), range (i_begin, i_end)));
0116                 // 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)));
0117                 const vector<value_type> e1_range (project (e1 (), range (i_begin, i_end)));
0118                 const matrix<value_type, column_major> e2_range (project (e2 (), range (i_begin, i_end), range (j_begin, j_end)));
0119 #endif
0120                 v_range.plus_assign (prod (e1_range, e2_range));
0121             }
0122 #ifndef BOOST_UBLAS_NO_CACHE
0123             project (v, range (j_begin, j_end)).assign (v_range);
0124 #endif
0125         }
0126 #if BOOST_UBLAS_TYPE_CHECK
0127         BOOST_UBLAS_CHECK (norm_1 (v - cv) <= 2 * std::numeric_limits<real_type>::epsilon () * verrorbound, internal_logic ());
0128 #endif
0129         return v;
0130     }
0131
0132     template<class M, typename M::size_type BS, class E1, class E2>
0133     BOOST_UBLAS_INLINE
0134     M
0135     block_prod (const matrix_expression<E1> &e1,
0136                 const matrix_expression<E2> &e2,
0137                 row_major_tag) {
0138         typedef M matrix_type;
0139         typedef const E1 expression1_type;
0140         typedef const E2 expression2_type;
0141         typedef typename M::size_type size_type;
0142         typedef typename M::value_type value_type;
0143         const size_type block_size = BS;
0144
0145         M m (e1 ().size1 (), e2 ().size2 ());
0146 #if BOOST_UBLAS_TYPE_CHECK
0147         matrix<value_type, row_major> cm (m.size1 (), m.size2 ());
0148         typedef typename type_traits<value_type>::real_type real_type;
0149         real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));
0150         indexing_matrix_assign<scalar_assign> (cm, prod (e1, e2), row_major_tag ());
0151         disable_type_check<bool>::value = true;
0152 #endif
0153         size_type i_size = e1 ().size1 ();
0154         size_type j_size = e2 ().size2 ();
0155         size_type k_size = BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size1 ());
0156         for (size_type i_begin = 0; i_begin < i_size; i_begin += block_size) {
0157             size_type i_end = i_begin + (std::min) (i_size - i_begin, block_size);
0158             for (size_type j_begin = 0; j_begin < j_size; j_begin += block_size) {
0159                 size_type j_end = j_begin + (std::min) (j_size - j_begin, block_size);
0160                 // FIX: never ignore Martin Weiser's advice ;-(
0161 #ifdef BOOST_UBLAS_NO_CACHE
0162                 matrix_range<matrix_type> m_range (m, range (i_begin, i_end), range (j_begin, j_end));
0163 #else
0164                 // matrix<value_type, row_major, bounded_array<value_type, block_size * block_size> > m_range (i_end - i_begin, j_end - j_begin);
0165                 matrix<value_type, row_major> m_range (i_end - i_begin, j_end - j_begin);
0166 #endif
0167                 m_range.assign (zero_matrix<value_type> (i_end - i_begin, j_end - j_begin));
0168                 for (size_type k_begin = 0; k_begin < k_size; k_begin += block_size) {
0169                     size_type k_end = k_begin + (std::min) (k_size - k_begin, block_size);
0170 #ifdef BOOST_UBLAS_NO_CACHE
0171                     const matrix_range<expression1_type> e1_range (e1 (), range (i_begin, i_end), range (k_begin, k_end));
0172                     const matrix_range<expression2_type> e2_range (e2 (), range (k_begin, k_end), range (j_begin, j_end));
0173 #else
0174                     // 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)));
0175                     // 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)));
0176                     const matrix<value_type, row_major> e1_range (project (e1 (), range (i_begin, i_end), range (k_begin, k_end)));
0177                     const matrix<value_type, column_major> e2_range (project (e2 (), range (k_begin, k_end), range (j_begin, j_end)));
0178 #endif
0179                     m_range.plus_assign (prod (e1_range, e2_range));
0180                 }
0181 #ifndef BOOST_UBLAS_NO_CACHE
0182                 project (m, range (i_begin, i_end), range (j_begin, j_end)).assign (m_range);
0183 #endif
0184             }
0185         }
0186 #if BOOST_UBLAS_TYPE_CHECK
0187         disable_type_check<bool>::value = false;
0188         BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());
0189 #endif
0190         return m;
0191     }
0192
0193     template<class M, typename M::size_type BS, class E1, class E2>
0194     BOOST_UBLAS_INLINE
0195     M
0196     block_prod (const matrix_expression<E1> &e1,
0197                 const matrix_expression<E2> &e2,
0198                 column_major_tag) {
0199         typedef M matrix_type;
0200         typedef const E1 expression1_type;
0201         typedef const E2 expression2_type;
0202         typedef typename M::size_type size_type;
0203         typedef typename M::value_type value_type;
0204         const size_type block_size = BS;
0205
0206         M m (e1 ().size1 (), e2 ().size2 ());
0207 #if BOOST_UBLAS_TYPE_CHECK
0208         matrix<value_type, column_major> cm (m.size1 (), m.size2 ());
0209         typedef typename type_traits<value_type>::real_type real_type;
0210         real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));
0211         indexing_matrix_assign<scalar_assign> (cm, prod (e1, e2), column_major_tag ());
0212         disable_type_check<bool>::value = true;
0213 #endif
0214         size_type i_size = e1 ().size1 ();
0215         size_type j_size = e2 ().size2 ();
0216         size_type k_size = BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size1 ());
0217         for (size_type j_begin = 0; j_begin < j_size; j_begin += block_size) {
0218             size_type j_end = j_begin + (std::min) (j_size - j_begin, block_size);
0219             for (size_type i_begin = 0; i_begin < i_size; i_begin += block_size) {
0220                 size_type i_end = i_begin + (std::min) (i_size - i_begin, block_size);
0221                 // FIX: never ignore Martin Weiser's advice ;-(
0222 #ifdef BOOST_UBLAS_NO_CACHE
0223                 matrix_range<matrix_type> m_range (m, range (i_begin, i_end), range (j_begin, j_end));
0224 #else
0225                 // matrix<value_type, column_major, bounded_array<value_type, block_size * block_size> > m_range (i_end - i_begin, j_end - j_begin);
0226                 matrix<value_type, column_major> m_range (i_end - i_begin, j_end - j_begin);
0227 #endif
0228                 m_range.assign (zero_matrix<value_type> (i_end - i_begin, j_end - j_begin));
0229                 for (size_type k_begin = 0; k_begin < k_size; k_begin += block_size) {
0230                     size_type k_end = k_begin + (std::min) (k_size - k_begin, block_size);
0231 #ifdef BOOST_UBLAS_NO_CACHE
0232                     const matrix_range<expression1_type> e1_range (e1 (), range (i_begin, i_end), range (k_begin, k_end));
0233                     const matrix_range<expression2_type> e2_range (e2 (), range (k_begin, k_end), range (j_begin, j_end));
0234 #else
0235                     // 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)));
0236                     // 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)));
0237                     const matrix<value_type, row_major> e1_range (project (e1 (), range (i_begin, i_end), range (k_begin, k_end)));
0238                     const matrix<value_type, column_major> e2_range (project (e2 (), range (k_begin, k_end), range (j_begin, j_end)));
0239 #endif
0240                     m_range.plus_assign (prod (e1_range, e2_range));
0241                 }
0242 #ifndef BOOST_UBLAS_NO_CACHE
0243                 project (m, range (i_begin, i_end), range (j_begin, j_end)).assign (m_range);
0244 #endif
0245             }
0246         }
0247 #if BOOST_UBLAS_TYPE_CHECK
0248         disable_type_check<bool>::value = false;
0249         BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());
0250 #endif
0251         return m;
0252     }
0253
0254     // Dispatcher
0255     template<class M, typename M::size_type BS, class E1, class E2>
0256     BOOST_UBLAS_INLINE
0257     M
0258     block_prod (const matrix_expression<E1> &e1,
0259                 const matrix_expression<E2> &e2) {
0260         typedef typename M::orientation_category orientation_category;
0261         return block_prod<M, BS> (e1, e2, orientation_category ());
0262     }
0263
0264 }}}
0265
0266 #endif