ceph/src/boost/libs/numeric/ublas/include/boost/numeric/ublas/operation_blocked.hpp

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