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 block_prod(const matrix_expression<E1> & e1,const vector_expression<E2> & e2)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 block_prod(const vector_expression<E1> & e1,const matrix_expression<E2> & e2)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 block_prod(const matrix_expression<E1> & e1,const matrix_expression<E2> & e2,row_major_tag)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 block_prod(const matrix_expression<E1> & e1,const matrix_expression<E2> & e2,column_major_tag)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 block_prod(const matrix_expression<E1> & e1,const matrix_expression<E2> & e2)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 267