numeric/ublas/operation_sparse.hpp

0001 //
0002 //  Copyright (c) 2000-2002
0003 //  Joerg Walter, Mathias Koch
0004 //
0005 //  Distributed under the Boost Software License, Version 1.0. (See
0006 //  accompanying file LICENSE_1_0.txt or copy at
0007 //  http://www.boost.org/LICENSE_1_0.txt)
0008 //
0009 //  The authors gratefully acknowledge the support of
0010 //  GeNeSys mbH & Co. KG in producing this work.
0011 //
0012
0013 #ifndef _BOOST_UBLAS_OPERATION_SPARSE_
0014 #define _BOOST_UBLAS_OPERATION_SPARSE_
0015
0016 #include <boost/numeric/ublas/traits.hpp>
0017
0018 // These scaled additions were borrowed from MTL unashamedly.
0019 // But Alexei Novakov had a lot of ideas to improve these. Thanks.
0020
0021 namespace boost { namespace numeric { namespace ublas {
0022
0023     template<class M, class E1, class E2, class TRI>
0024     BOOST_UBLAS_INLINE
0025     M &
0026     sparse_prod (const matrix_expression<E1> &e1,
0027                  const matrix_expression<E2> &e2,
0028                  M &m, TRI,
0029                  row_major_tag) {
0030         typedef M matrix_type;
0031         typedef TRI triangular_restriction;
0032         typedef const E1 expression1_type;
0033         typedef const E2 expression2_type;
0034         typedef typename M::size_type size_type;
0035         typedef typename M::value_type value_type;
0036
0037         // ISSUE why is there a dense vector here?
0038         vector<value_type> temporary (e2 ().size2 ());
0039         temporary.clear ();
0040         typename expression1_type::const_iterator1 it1 (e1 ().begin1 ());
0041         typename expression1_type::const_iterator1 it1_end (e1 ().end1 ());
0042         while (it1 != it1_end) {
0043             size_type jb (temporary.size ());
0044             size_type je (0);
0045 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
0046             typename expression1_type::const_iterator2 it2 (it1.begin ());
0047             typename expression1_type::const_iterator2 it2_end (it1.end ());
0048 #else
0049             typename expression1_type::const_iterator2 it2 (boost::numeric::ublas::begin (it1, iterator1_tag ()));
0050             typename expression1_type::const_iterator2 it2_end (boost::numeric::ublas::end (it1, iterator1_tag ()));
0051 #endif
0052             while (it2 != it2_end) {
0053                 // temporary.plus_assign (*it2 * row (e2 (), it2.index2 ()));
0054                 matrix_row<expression2_type> mr (e2 (), it2.index2 ());
0055                 typename matrix_row<expression2_type>::const_iterator itr (mr.begin ());
0056                 typename matrix_row<expression2_type>::const_iterator itr_end (mr.end ());
0057                 while (itr != itr_end) {
0058                     size_type j (itr.index ());
0059                     temporary (j) += *it2 * *itr;
0060                     jb = (std::min) (jb, j);
0061                     je = (std::max) (je, j);
0062                     ++ itr;
0063                 }
0064                 ++ it2;
0065             }
0066             for (size_type j = jb; j < je + 1; ++ j) {
0067                 if (temporary (j) != value_type/*zero*/()) {
0068                     // FIXME we'll need to extend the container interface!
0069                     // m.push_back (it1.index1 (), j, temporary (j));
0070                     // FIXME What to do with adaptors?
0071                     // m.insert (it1.index1 (), j, temporary (j));
0072                     if (triangular_restriction::other (it1.index1 (), j))
0073                         m (it1.index1 (), j) = temporary (j);
0074                     temporary (j) = value_type/*zero*/();
0075                 }
0076             }
0077             ++ it1;
0078         }
0079         return m;
0080     }
0081
0082     template<class M, class E1, class E2, class TRI>
0083     BOOST_UBLAS_INLINE
0084     M &
0085     sparse_prod (const matrix_expression<E1> &e1,
0086                  const matrix_expression<E2> &e2,
0087                  M &m, TRI,
0088                  column_major_tag) {
0089         typedef M matrix_type;
0090         typedef TRI triangular_restriction;
0091         typedef const E1 expression1_type;
0092         typedef const E2 expression2_type;
0093         typedef typename M::size_type size_type;
0094         typedef typename M::value_type value_type;
0095
0096         // ISSUE why is there a dense vector here?
0097         vector<value_type> temporary (e1 ().size1 ());
0098         temporary.clear ();
0099         typename expression2_type::const_iterator2 it2 (e2 ().begin2 ());
0100         typename expression2_type::const_iterator2 it2_end (e2 ().end2 ());
0101         while (it2 != it2_end) {
0102             size_type ib (temporary.size ());
0103             size_type ie (0);
0104 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
0105             typename expression2_type::const_iterator1 it1 (it2.begin ());
0106             typename expression2_type::const_iterator1 it1_end (it2.end ());
0107 #else
0108             typename expression2_type::const_iterator1 it1 (boost::numeric::ublas::begin (it2, iterator2_tag ()));
0109             typename expression2_type::const_iterator1 it1_end (boost::numeric::ublas::end (it2, iterator2_tag ()));
0110 #endif
0111             while (it1 != it1_end) {
0112                 // column (m, it2.index2 ()).plus_assign (*it1 * column (e1 (), it1.index1 ()));
0113                 matrix_column<expression1_type> mc (e1 (), it1.index1 ());
0114                 typename matrix_column<expression1_type>::const_iterator itc (mc.begin ());
0115                 typename matrix_column<expression1_type>::const_iterator itc_end (mc.end ());
0116                 while (itc != itc_end) {
0117                     size_type i (itc.index ());
0118                     temporary (i) += *it1 * *itc;
0119                     ib = (std::min) (ib, i);
0120                     ie = (std::max) (ie, i);
0121                     ++ itc;
0122                 }
0123                 ++ it1;
0124             }
0125             for (size_type i = ib; i < ie + 1; ++ i) {
0126                 if (temporary (i) != value_type/*zero*/()) {
0127                     // FIXME we'll need to extend the container interface!
0128                     // m.push_back (i, it2.index2 (), temporary (i));
0129                     // FIXME What to do with adaptors?
0130                     // m.insert (i, it2.index2 (), temporary (i));
0131                     if (triangular_restriction::other (i, it2.index2 ()))
0132                         m (i, it2.index2 ()) = temporary (i);
0133                     temporary (i) = value_type/*zero*/();
0134                 }
0135             }
0136             ++ it2;
0137         }
0138         return m;
0139     }
0140
0141     // Dispatcher
0142     template<class M, class E1, class E2, class TRI>
0143     BOOST_UBLAS_INLINE
0144     M &
0145     sparse_prod (const matrix_expression<E1> &e1,
0146                  const matrix_expression<E2> &e2,
0147                  M &m, TRI, bool init = true) {
0148         typedef typename M::value_type value_type;
0149         typedef TRI triangular_restriction;
0150         typedef typename M::orientation_category orientation_category;
0151
0152         if (init)
0153             m.assign (zero_matrix<value_type> (e1 ().size1 (), e2 ().size2 ()));
0154         return sparse_prod (e1, e2, m, triangular_restriction (), orientation_category ());
0155     }
0156     template<class M, class E1, class E2, class TRI>
0157     BOOST_UBLAS_INLINE
0158     M
0159     sparse_prod (const matrix_expression<E1> &e1,
0160                  const matrix_expression<E2> &e2,
0161                  TRI) {
0162         typedef M matrix_type;
0163         typedef TRI triangular_restriction;
0164
0165         matrix_type m (e1 ().size1 (), e2 ().size2 ());
0166         // FIXME needed for c_matrix?!
0167         // return sparse_prod (e1, e2, m, triangular_restriction (), false);
0168         return sparse_prod (e1, e2, m, triangular_restriction (), true);
0169     }
0170     template<class M, class E1, class E2>
0171     BOOST_UBLAS_INLINE
0172     M &
0173     sparse_prod (const matrix_expression<E1> &e1,
0174                  const matrix_expression<E2> &e2,
0175                  M &m, bool init = true) {
0176         typedef typename M::value_type value_type;
0177         typedef typename M::orientation_category orientation_category;
0178
0179         if (init)
0180             m.assign (zero_matrix<value_type> (e1 ().size1 (), e2 ().size2 ()));
0181         return sparse_prod (e1, e2, m, full (), orientation_category ());
0182     }
0183     template<class M, class E1, class E2>
0184     BOOST_UBLAS_INLINE
0185     M
0186     sparse_prod (const matrix_expression<E1> &e1,
0187                  const matrix_expression<E2> &e2) {
0188         typedef M matrix_type;
0189
0190         matrix_type m (e1 ().size1 (), e2 ().size2 ());
0191         // FIXME needed for c_matrix?!
0192         // return sparse_prod (e1, e2, m, full (), false);
0193         return sparse_prod (e1, e2, m, full (), true);
0194     }
0195
0196 }}}
0197
0198 #endif