ublas: operation_blocked.hpp Source File

- Main Page
- Related Pages
- Classes
- Files

- File List

# operation_blocked.hpp

    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