ublas: operation_blocked.hpp Source File

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