...one of the most highly
regarded and expertly designed C++ library projects in the
world.
— Herb Sutter and Andrei
Alexandrescu, C++
Coding Standards
00001 // 00002 // Copyright (c) 2000-2002 00003 // Joerg Walter, Mathias Koch 00004 // 00005 // Distributed under the Boost Software License, Version 1.0. (See 00006 // accompanying file LICENSE_1_0.txt or copy at 00007 // http://www.boost.org/LICENSE_1_0.txt) 00008 // 00009 // The authors gratefully acknowledge the support of 00010 // GeNeSys mbH & Co. KG in producing this work. 00011 // 00012 00013 #ifndef _BOOST_UBLAS_OPERATION_BLOCKED_ 00014 #define _BOOST_UBLAS_OPERATION_BLOCKED_ 00015 00016 #include <boost/numeric/ublas/traits.hpp> 00017 #include <boost/numeric/ublas/detail/vector_assign.hpp> // indexing_vector_assign 00018 #include <boost/numeric/ublas/detail/matrix_assign.hpp> // indexing_matrix_assign 00019 00020 00021 namespace boost { namespace numeric { namespace ublas { 00022 00023 template<class V, typename V::size_type BS, class E1, class E2> 00024 BOOST_UBLAS_INLINE 00025 V 00026 block_prod (const matrix_expression<E1> &e1, 00027 const vector_expression<E2> &e2) { 00028 typedef V vector_type; 00029 typedef const E1 expression1_type; 00030 typedef const E2 expression2_type; 00031 typedef typename V::size_type size_type; 00032 typedef typename V::value_type value_type; 00033 const size_type block_size = BS; 00034 00035 V v (e1 ().size1 ()); 00036 #if BOOST_UBLAS_TYPE_CHECK 00037 vector<value_type> cv (v.size ()); 00038 typedef typename type_traits<value_type>::real_type real_type; 00039 real_type verrorbound (norm_1 (v) + norm_1 (e1) * norm_1 (e2)); 00040 indexing_vector_assign<scalar_assign> (cv, prod (e1, e2)); 00041 #endif 00042 size_type i_size = e1 ().size1 (); 00043 size_type j_size = BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size ()); 00044 for (size_type i_begin = 0; i_begin < i_size; i_begin += block_size) { 00045 size_type i_end = i_begin + (std::min) (i_size - i_begin, block_size); 00046 // FIX: never ignore Martin Weiser's advice ;-( 00047 #ifdef BOOST_UBLAS_NO_CACHE 00048 vector_range<vector_type> v_range (v, range (i_begin, i_end)); 00049 #else 00050 // vector<value_type, bounded_array<value_type, block_size> > v_range (i_end - i_begin); 00051 vector<value_type> v_range (i_end - i_begin); 00052 #endif 00053 v_range.assign (zero_vector<value_type> (i_end - i_begin)); 00054 for (size_type j_begin = 0; j_begin < j_size; j_begin += block_size) { 00055 size_type j_end = j_begin + (std::min) (j_size - j_begin, block_size); 00056 #ifdef BOOST_UBLAS_NO_CACHE 00057 const matrix_range<expression1_type> e1_range (e1 (), range (i_begin, i_end), range (j_begin, j_end)); 00058 const vector_range<expression2_type> e2_range (e2 (), range (j_begin, j_end)); 00059 v_range.plus_assign (prod (e1_range, e2_range)); 00060 #else 00061 // 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))); 00062 // const vector<value_type, bounded_array<value_type, block_size> > e2_range (project (e2 (), range (j_begin, j_end))); 00063 const matrix<value_type, row_major> e1_range (project (e1 (), range (i_begin, i_end), range (j_begin, j_end))); 00064 const vector<value_type> e2_range (project (e2 (), range (j_begin, j_end))); 00065 v_range.plus_assign (prod (e1_range, e2_range)); 00066 #endif 00067 } 00068 #ifndef BOOST_UBLAS_NO_CACHE 00069 project (v, range (i_begin, i_end)).assign (v_range); 00070 #endif 00071 } 00072 #if BOOST_UBLAS_TYPE_CHECK 00073 BOOST_UBLAS_CHECK (norm_1 (v - cv) <= 2 * std::numeric_limits<real_type>::epsilon () * verrorbound, internal_logic ()); 00074 #endif 00075 return v; 00076 } 00077 00078 template<class V, typename V::size_type BS, class E1, class E2> 00079 BOOST_UBLAS_INLINE 00080 V 00081 block_prod (const vector_expression<E1> &e1, 00082 const matrix_expression<E2> &e2) { 00083 typedef V vector_type; 00084 typedef const E1 expression1_type; 00085 typedef const E2 expression2_type; 00086 typedef typename V::size_type size_type; 00087 typedef typename V::value_type value_type; 00088 const size_type block_size = BS; 00089 00090 V v (e2 ().size2 ()); 00091 #if BOOST_UBLAS_TYPE_CHECK 00092 vector<value_type> cv (v.size ()); 00093 typedef typename type_traits<value_type>::real_type real_type; 00094 real_type verrorbound (norm_1 (v) + norm_1 (e1) * norm_1 (e2)); 00095 indexing_vector_assign<scalar_assign> (cv, prod (e1, e2)); 00096 #endif 00097 size_type i_size = BOOST_UBLAS_SAME (e1 ().size (), e2 ().size1 ()); 00098 size_type j_size = e2 ().size2 (); 00099 for (size_type j_begin = 0; j_begin < j_size; j_begin += block_size) { 00100 size_type j_end = j_begin + (std::min) (j_size - j_begin, block_size); 00101 // FIX: never ignore Martin Weiser's advice ;-( 00102 #ifdef BOOST_UBLAS_NO_CACHE 00103 vector_range<vector_type> v_range (v, range (j_begin, j_end)); 00104 #else 00105 // vector<value_type, bounded_array<value_type, block_size> > v_range (j_end - j_begin); 00106 vector<value_type> v_range (j_end - j_begin); 00107 #endif 00108 v_range.assign (zero_vector<value_type> (j_end - j_begin)); 00109 for (size_type i_begin = 0; i_begin < i_size; i_begin += block_size) { 00110 size_type i_end = i_begin + (std::min) (i_size - i_begin, block_size); 00111 #ifdef BOOST_UBLAS_NO_CACHE 00112 const vector_range<expression1_type> e1_range (e1 (), range (i_begin, i_end)); 00113 const matrix_range<expression2_type> e2_range (e2 (), range (i_begin, i_end), range (j_begin, j_end)); 00114 #else 00115 // const vector<value_type, bounded_array<value_type, block_size> > e1_range (project (e1 (), range (i_begin, i_end))); 00116 // 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))); 00117 const vector<value_type> e1_range (project (e1 (), range (i_begin, i_end))); 00118 const matrix<value_type, column_major> e2_range (project (e2 (), range (i_begin, i_end), range (j_begin, j_end))); 00119 #endif 00120 v_range.plus_assign (prod (e1_range, e2_range)); 00121 } 00122 #ifndef BOOST_UBLAS_NO_CACHE 00123 project (v, range (j_begin, j_end)).assign (v_range); 00124 #endif 00125 } 00126 #if BOOST_UBLAS_TYPE_CHECK 00127 BOOST_UBLAS_CHECK (norm_1 (v - cv) <= 2 * std::numeric_limits<real_type>::epsilon () * verrorbound, internal_logic ()); 00128 #endif 00129 return v; 00130 } 00131 00132 template<class M, typename M::size_type BS, class E1, class E2> 00133 BOOST_UBLAS_INLINE 00134 M 00135 block_prod (const matrix_expression<E1> &e1, 00136 const matrix_expression<E2> &e2, 00137 row_major_tag) { 00138 typedef M matrix_type; 00139 typedef const E1 expression1_type; 00140 typedef const E2 expression2_type; 00141 typedef typename M::size_type size_type; 00142 typedef typename M::value_type value_type; 00143 const size_type block_size = BS; 00144 00145 M m (e1 ().size1 (), e2 ().size2 ()); 00146 #if BOOST_UBLAS_TYPE_CHECK 00147 matrix<value_type, row_major> cm (m.size1 (), m.size2 ()); 00148 typedef typename type_traits<value_type>::real_type real_type; 00149 real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2)); 00150 indexing_matrix_assign<scalar_assign> (cm, prod (e1, e2), row_major_tag ()); 00151 disable_type_check<bool>::value = true; 00152 #endif 00153 size_type i_size = e1 ().size1 (); 00154 size_type j_size = e2 ().size2 (); 00155 size_type k_size = BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size1 ()); 00156 for (size_type i_begin = 0; i_begin < i_size; i_begin += block_size) { 00157 size_type i_end = i_begin + (std::min) (i_size - i_begin, block_size); 00158 for (size_type j_begin = 0; j_begin < j_size; j_begin += block_size) { 00159 size_type j_end = j_begin + (std::min) (j_size - j_begin, block_size); 00160 // FIX: never ignore Martin Weiser's advice ;-( 00161 #ifdef BOOST_UBLAS_NO_CACHE 00162 matrix_range<matrix_type> m_range (m, range (i_begin, i_end), range (j_begin, j_end)); 00163 #else 00164 // matrix<value_type, row_major, bounded_array<value_type, block_size * block_size> > m_range (i_end - i_begin, j_end - j_begin); 00165 matrix<value_type, row_major> m_range (i_end - i_begin, j_end - j_begin); 00166 #endif 00167 m_range.assign (zero_matrix<value_type> (i_end - i_begin, j_end - j_begin)); 00168 for (size_type k_begin = 0; k_begin < k_size; k_begin += block_size) { 00169 size_type k_end = k_begin + (std::min) (k_size - k_begin, block_size); 00170 #ifdef BOOST_UBLAS_NO_CACHE 00171 const matrix_range<expression1_type> e1_range (e1 (), range (i_begin, i_end), range (k_begin, k_end)); 00172 const matrix_range<expression2_type> e2_range (e2 (), range (k_begin, k_end), range (j_begin, j_end)); 00173 #else 00174 // 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))); 00175 // 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))); 00176 const matrix<value_type, row_major> e1_range (project (e1 (), range (i_begin, i_end), range (k_begin, k_end))); 00177 const matrix<value_type, column_major> e2_range (project (e2 (), range (k_begin, k_end), range (j_begin, j_end))); 00178 #endif 00179 m_range.plus_assign (prod (e1_range, e2_range)); 00180 } 00181 #ifndef BOOST_UBLAS_NO_CACHE 00182 project (m, range (i_begin, i_end), range (j_begin, j_end)).assign (m_range); 00183 #endif 00184 } 00185 } 00186 #if BOOST_UBLAS_TYPE_CHECK 00187 disable_type_check<bool>::value = false; 00188 BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ()); 00189 #endif 00190 return m; 00191 } 00192 00193 template<class M, typename M::size_type BS, class E1, class E2> 00194 BOOST_UBLAS_INLINE 00195 M 00196 block_prod (const matrix_expression<E1> &e1, 00197 const matrix_expression<E2> &e2, 00198 column_major_tag) { 00199 typedef M matrix_type; 00200 typedef const E1 expression1_type; 00201 typedef const E2 expression2_type; 00202 typedef typename M::size_type size_type; 00203 typedef typename M::value_type value_type; 00204 const size_type block_size = BS; 00205 00206 M m (e1 ().size1 (), e2 ().size2 ()); 00207 #if BOOST_UBLAS_TYPE_CHECK 00208 matrix<value_type, column_major> cm (m.size1 (), m.size2 ()); 00209 typedef typename type_traits<value_type>::real_type real_type; 00210 real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2)); 00211 indexing_matrix_assign<scalar_assign> (cm, prod (e1, e2), column_major_tag ()); 00212 disable_type_check<bool>::value = true; 00213 #endif 00214 size_type i_size = e1 ().size1 (); 00215 size_type j_size = e2 ().size2 (); 00216 size_type k_size = BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size1 ()); 00217 for (size_type j_begin = 0; j_begin < j_size; j_begin += block_size) { 00218 size_type j_end = j_begin + (std::min) (j_size - j_begin, block_size); 00219 for (size_type i_begin = 0; i_begin < i_size; i_begin += block_size) { 00220 size_type i_end = i_begin + (std::min) (i_size - i_begin, block_size); 00221 // FIX: never ignore Martin Weiser's advice ;-( 00222 #ifdef BOOST_UBLAS_NO_CACHE 00223 matrix_range<matrix_type> m_range (m, range (i_begin, i_end), range (j_begin, j_end)); 00224 #else 00225 // matrix<value_type, column_major, bounded_array<value_type, block_size * block_size> > m_range (i_end - i_begin, j_end - j_begin); 00226 matrix<value_type, column_major> m_range (i_end - i_begin, j_end - j_begin); 00227 #endif 00228 m_range.assign (zero_matrix<value_type> (i_end - i_begin, j_end - j_begin)); 00229 for (size_type k_begin = 0; k_begin < k_size; k_begin += block_size) { 00230 size_type k_end = k_begin + (std::min) (k_size - k_begin, block_size); 00231 #ifdef BOOST_UBLAS_NO_CACHE 00232 const matrix_range<expression1_type> e1_range (e1 (), range (i_begin, i_end), range (k_begin, k_end)); 00233 const matrix_range<expression2_type> e2_range (e2 (), range (k_begin, k_end), range (j_begin, j_end)); 00234 #else 00235 // 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))); 00236 // 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))); 00237 const matrix<value_type, row_major> e1_range (project (e1 (), range (i_begin, i_end), range (k_begin, k_end))); 00238 const matrix<value_type, column_major> e2_range (project (e2 (), range (k_begin, k_end), range (j_begin, j_end))); 00239 #endif 00240 m_range.plus_assign (prod (e1_range, e2_range)); 00241 } 00242 #ifndef BOOST_UBLAS_NO_CACHE 00243 project (m, range (i_begin, i_end), range (j_begin, j_end)).assign (m_range); 00244 #endif 00245 } 00246 } 00247 #if BOOST_UBLAS_TYPE_CHECK 00248 disable_type_check<bool>::value = false; 00249 BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ()); 00250 #endif 00251 return m; 00252 } 00253 00254 // Dispatcher 00255 template<class M, typename M::size_type BS, class E1, class E2> 00256 BOOST_UBLAS_INLINE 00257 M 00258 block_prod (const matrix_expression<E1> &e1, 00259 const matrix_expression<E2> &e2) { 00260 typedef typename M::orientation_category orientation_category; 00261 return block_prod<M, BS> (e1, e2, orientation_category ()); 00262 } 00263 00264 }}} 00265 00266 #endif