boost/numeric/ublas/operation_sparse.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_SPARSE_
  14 #define _BOOST_UBLAS_OPERATION_SPARSE_
  15
  16 #include <boost/numeric/ublas/traits.hpp>
  17
  18 // These scaled additions were borrowed from MTL unashamedly.
  19 // But Alexei Novakov had a lot of ideas to improve these. Thanks.
  20
  21 namespace boost { namespace numeric { namespace ublas {
  22
  23     template<class M, class E1, class E2, class TRI>
  24     BOOST_UBLAS_INLINE
  25     M &
  26     sparse_prod (const matrix_expression<E1> &e1,
  27                  const matrix_expression<E2> &e2,
  28                  M &m, TRI,
  29                  row_major_tag) {
  30         typedef M matrix_type;
  31         typedef TRI triangular_restriction;
  32         typedef const E1 expression1_type;
  33         typedef const E2 expression2_type;
  34         typedef typename M::size_type size_type;
  35         typedef typename M::value_type value_type;
  36
  37         // ISSUE why is there a dense vector here?
  38         vector<value_type> temporary (e2 ().size2 ());
  39         temporary.clear ();
  40         typename expression1_type::const_iterator1 it1 (e1 ().begin1 ());
  41         typename expression1_type::const_iterator1 it1_end (e1 ().end1 ());
  42         while (it1 != it1_end) {
  43             size_type jb (temporary.size ());
  44             size_type je (0);
  45 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
  46             typename expression1_type::const_iterator2 it2 (it1.begin ());
  47             typename expression1_type::const_iterator2 it2_end (it1.end ());
  48 #else
  49             typename expression1_type::const_iterator2 it2 (boost::numeric::ublas::begin (it1, iterator1_tag ()));
  50             typename expression1_type::const_iterator2 it2_end (boost::numeric::ublas::end (it1, iterator1_tag ()));
  51 #endif
  52             while (it2 != it2_end) {
  53                 // temporary.plus_assign (*it2 * row (e2 (), it2.index2 ()));
  54                 matrix_row<expression2_type> mr (e2 (), it2.index2 ());
  55                 typename matrix_row<expression2_type>::const_iterator itr (mr.begin ());
  56                 typename matrix_row<expression2_type>::const_iterator itr_end (mr.end ());
  57                 while (itr != itr_end) {
  58                     size_type j (itr.index ());
  59                     temporary (j) += *it2 * *itr;
  60                     jb = (std::min) (jb, j);
  61                     je = (std::max) (je, j);
  62                     ++ itr;
  63                 }
  64                 ++ it2;
  65             }
  66             for (size_type j = jb; j < je + 1; ++ j) {
  67                 if (temporary (j) != value_type/*zero*/()) {
  68                     // FIXME we'll need to extend the container interface!
  69                     // m.push_back (it1.index1 (), j, temporary (j));
  70                     // FIXME What to do with adaptors?
  71                     // m.insert (it1.index1 (), j, temporary (j));
  72                     if (triangular_restriction::other (it1.index1 (), j))
  73                         m (it1.index1 (), j) = temporary (j);
  74                     temporary (j) = value_type/*zero*/();
  75                 }
  76             }
  77             ++ it1;
  78         }
  79         return m;
  80     }
  81
  82     template<class M, class E1, class E2, class TRI>
  83     BOOST_UBLAS_INLINE
  84     M &
  85     sparse_prod (const matrix_expression<E1> &e1,
  86                  const matrix_expression<E2> &e2,
  87                  M &m, TRI,
  88                  column_major_tag) {
  89         typedef M matrix_type;
  90         typedef TRI triangular_restriction;
  91         typedef const E1 expression1_type;
  92         typedef const E2 expression2_type;
  93         typedef typename M::size_type size_type;
  94         typedef typename M::value_type value_type;
  95
  96         // ISSUE why is there a dense vector here?
  97         vector<value_type> temporary (e1 ().size1 ());
  98         temporary.clear ();
  99         typename expression2_type::const_iterator2 it2 (e2 ().begin2 ());
 100         typename expression2_type::const_iterator2 it2_end (e2 ().end2 ());
 101         while (it2 != it2_end) {
 102             size_type ib (temporary.size ());
 103             size_type ie (0);
 104 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
 105             typename expression2_type::const_iterator1 it1 (it2.begin ());
 106             typename expression2_type::const_iterator1 it1_end (it2.end ());
 107 #else
 108             typename expression2_type::const_iterator1 it1 (boost::numeric::ublas::begin (it2, iterator2_tag ()));
 109             typename expression2_type::const_iterator1 it1_end (boost::numeric::ublas::end (it2, iterator2_tag ()));
 110 #endif
 111             while (it1 != it1_end) {
 112                 // column (m, it2.index2 ()).plus_assign (*it1 * column (e1 (), it1.index1 ()));
 113                 matrix_column<expression1_type> mc (e1 (), it1.index1 ());
 114                 typename matrix_column<expression1_type>::const_iterator itc (mc.begin ());
 115                 typename matrix_column<expression1_type>::const_iterator itc_end (mc.end ());
 116                 while (itc != itc_end) {
 117                     size_type i (itc.index ());
 118                     temporary (i) += *it1 * *itc;
 119                     ib = (std::min) (ib, i);
 120                     ie = (std::max) (ie, i);
 121                     ++ itc;
 122                 }
 123                 ++ it1;
 124             }
 125             for (size_type i = ib; i < ie + 1; ++ i) {
 126                 if (temporary (i) != value_type/*zero*/()) {
 127                     // FIXME we'll need to extend the container interface!
 128                     // m.push_back (i, it2.index2 (), temporary (i));
 129                     // FIXME What to do with adaptors?
 130                     // m.insert (i, it2.index2 (), temporary (i));
 131                     if (triangular_restriction::other (i, it2.index2 ()))
 132                         m (i, it2.index2 ()) = temporary (i);
 133                     temporary (i) = value_type/*zero*/();
 134                 }
 135             }
 136             ++ it2;
 137         }
 138         return m;
 139     }
 140
 141     // Dispatcher
 142     template<class M, class E1, class E2, class TRI>
 143     BOOST_UBLAS_INLINE
 144     M &
 145     sparse_prod (const matrix_expression<E1> &e1,
 146                  const matrix_expression<E2> &e2,
 147                  M &m, TRI, bool init = true) {
 148         typedef typename M::value_type value_type;
 149         typedef TRI triangular_restriction;
 150         typedef typename M::orientation_category orientation_category;
 151
 152         if (init)
 153             m.assign (zero_matrix<value_type> (e1 ().size1 (), e2 ().size2 ()));
 154         return sparse_prod (e1, e2, m, triangular_restriction (), orientation_category ());
 155     }
 156     template<class M, class E1, class E2, class TRI>
 157     BOOST_UBLAS_INLINE
 158     M
 159     sparse_prod (const matrix_expression<E1> &e1,
 160                  const matrix_expression<E2> &e2,
 161                  TRI) {
 162         typedef M matrix_type;
 163         typedef TRI triangular_restriction;
 164
 165         matrix_type m (e1 ().size1 (), e2 ().size2 ());
 166         // FIXME needed for c_matrix?!
 167         // return sparse_prod (e1, e2, m, triangular_restriction (), false);
 168         return sparse_prod (e1, e2, m, triangular_restriction (), true);
 169     }
 170     template<class M, class E1, class E2>
 171     BOOST_UBLAS_INLINE
 172     M &
 173     sparse_prod (const matrix_expression<E1> &e1,
 174                  const matrix_expression<E2> &e2,
 175                  M &m, bool init = true) {
 176         typedef typename M::value_type value_type;
 177         typedef typename M::orientation_category orientation_category;
 178
 179         if (init)
 180             m.assign (zero_matrix<value_type> (e1 ().size1 (), e2 ().size2 ()));
 181         return sparse_prod (e1, e2, m, full (), orientation_category ());
 182     }
 183     template<class M, class E1, class E2>
 184     BOOST_UBLAS_INLINE
 185     M
 186     sparse_prod (const matrix_expression<E1> &e1,
 187                  const matrix_expression<E2> &e2) {
 188         typedef M matrix_type;
 189
 190         matrix_type m (e1 ().size1 (), e2 ().size2 ());
 191         // FIXME needed for c_matrix?!
 192         // return sparse_prod (e1, e2, m, full (), false);
 193         return sparse_prod (e1, e2, m, full (), true);
 194     }
 195
 196 }}}
 197
 198 #endif