Chris@16: //
Chris@16: //  Copyright (c) 2000-2002
Chris@16: //  Joerg Walter, Mathias Koch
Chris@16: //
Chris@16: //  Distributed under the Boost Software License, Version 1.0. (See
Chris@16: //  accompanying file LICENSE_1_0.txt or copy at
Chris@16: //  http://www.boost.org/LICENSE_1_0.txt)
Chris@16: //
Chris@16: //  The authors gratefully acknowledge the support of
Chris@16: //  GeNeSys mbH & Co. KG in producing this work.
Chris@16: //
Chris@16: 
Chris@16: #ifndef _BOOST_UBLAS_OPERATION_SPARSE_
Chris@16: #define _BOOST_UBLAS_OPERATION_SPARSE_
Chris@16: 
Chris@16: #include <boost/numeric/ublas/traits.hpp>
Chris@16: 
Chris@16: // These scaled additions were borrowed from MTL unashamedly.
Chris@16: // But Alexei Novakov had a lot of ideas to improve these. Thanks.
Chris@16: 
Chris@16: namespace boost { namespace numeric { namespace ublas {
Chris@16: 
Chris@16:     template<class M, class E1, class E2, class TRI>
Chris@16:     BOOST_UBLAS_INLINE
Chris@16:     M &
Chris@16:     sparse_prod (const matrix_expression<E1> &e1,
Chris@16:                  const matrix_expression<E2> &e2,
Chris@16:                  M &m, TRI,
Chris@16:                  row_major_tag) {
Chris@16:         typedef M matrix_type;
Chris@16:         typedef TRI triangular_restriction;
Chris@16:         typedef const E1 expression1_type;
Chris@16:         typedef const E2 expression2_type;
Chris@16:         typedef typename M::size_type size_type;
Chris@16:         typedef typename M::value_type value_type;
Chris@16: 
Chris@16:         // ISSUE why is there a dense vector here?
Chris@16:         vector<value_type> temporary (e2 ().size2 ());
Chris@16:         temporary.clear ();
Chris@16:         typename expression1_type::const_iterator1 it1 (e1 ().begin1 ());
Chris@16:         typename expression1_type::const_iterator1 it1_end (e1 ().end1 ());
Chris@16:         while (it1 != it1_end) {
Chris@16:             size_type jb (temporary.size ());
Chris@16:             size_type je (0);
Chris@16: #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
Chris@16:             typename expression1_type::const_iterator2 it2 (it1.begin ());
Chris@16:             typename expression1_type::const_iterator2 it2_end (it1.end ());
Chris@16: #else
Chris@16:             typename expression1_type::const_iterator2 it2 (boost::numeric::ublas::begin (it1, iterator1_tag ()));
Chris@16:             typename expression1_type::const_iterator2 it2_end (boost::numeric::ublas::end (it1, iterator1_tag ()));
Chris@16: #endif
Chris@16:             while (it2 != it2_end) {
Chris@16:                 // temporary.plus_assign (*it2 * row (e2 (), it2.index2 ()));
Chris@16:                 matrix_row<expression2_type> mr (e2 (), it2.index2 ());
Chris@16:                 typename matrix_row<expression2_type>::const_iterator itr (mr.begin ());
Chris@16:                 typename matrix_row<expression2_type>::const_iterator itr_end (mr.end ());
Chris@16:                 while (itr != itr_end) {
Chris@16:                     size_type j (itr.index ());
Chris@16:                     temporary (j) += *it2 * *itr;
Chris@16:                     jb = (std::min) (jb, j);
Chris@16:                     je = (std::max) (je, j);
Chris@16:                     ++ itr;
Chris@16:                 }
Chris@16:                 ++ it2;
Chris@16:             }
Chris@16:             for (size_type j = jb; j < je + 1; ++ j) {
Chris@16:                 if (temporary (j) != value_type/*zero*/()) {
Chris@16:                     // FIXME we'll need to extend the container interface!
Chris@16:                     // m.push_back (it1.index1 (), j, temporary (j));
Chris@16:                     // FIXME What to do with adaptors?
Chris@16:                     // m.insert (it1.index1 (), j, temporary (j));
Chris@16:                     if (triangular_restriction::other (it1.index1 (), j))
Chris@16:                         m (it1.index1 (), j) = temporary (j);
Chris@16:                     temporary (j) = value_type/*zero*/();
Chris@16:                 }
Chris@16:             }
Chris@16:             ++ it1;
Chris@16:         }
Chris@16:         return m;
Chris@16:     }
Chris@16: 
Chris@16:     template<class M, class E1, class E2, class TRI>
Chris@16:     BOOST_UBLAS_INLINE
Chris@16:     M &
Chris@16:     sparse_prod (const matrix_expression<E1> &e1,
Chris@16:                  const matrix_expression<E2> &e2,
Chris@16:                  M &m, TRI,
Chris@16:                  column_major_tag) {
Chris@16:         typedef M matrix_type;
Chris@16:         typedef TRI triangular_restriction;
Chris@16:         typedef const E1 expression1_type;
Chris@16:         typedef const E2 expression2_type;
Chris@16:         typedef typename M::size_type size_type;
Chris@16:         typedef typename M::value_type value_type;
Chris@16: 
Chris@16:         // ISSUE why is there a dense vector here?
Chris@16:         vector<value_type> temporary (e1 ().size1 ());
Chris@16:         temporary.clear ();
Chris@16:         typename expression2_type::const_iterator2 it2 (e2 ().begin2 ());
Chris@16:         typename expression2_type::const_iterator2 it2_end (e2 ().end2 ());
Chris@16:         while (it2 != it2_end) {
Chris@16:             size_type ib (temporary.size ());
Chris@16:             size_type ie (0);
Chris@16: #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
Chris@16:             typename expression2_type::const_iterator1 it1 (it2.begin ());
Chris@16:             typename expression2_type::const_iterator1 it1_end (it2.end ());
Chris@16: #else
Chris@16:             typename expression2_type::const_iterator1 it1 (boost::numeric::ublas::begin (it2, iterator2_tag ()));
Chris@16:             typename expression2_type::const_iterator1 it1_end (boost::numeric::ublas::end (it2, iterator2_tag ()));
Chris@16: #endif
Chris@16:             while (it1 != it1_end) {
Chris@16:                 // column (m, it2.index2 ()).plus_assign (*it1 * column (e1 (), it1.index1 ()));
Chris@16:                 matrix_column<expression1_type> mc (e1 (), it1.index1 ());
Chris@16:                 typename matrix_column<expression1_type>::const_iterator itc (mc.begin ());
Chris@16:                 typename matrix_column<expression1_type>::const_iterator itc_end (mc.end ());
Chris@16:                 while (itc != itc_end) {
Chris@16:                     size_type i (itc.index ());
Chris@16:                     temporary (i) += *it1 * *itc;
Chris@16:                     ib = (std::min) (ib, i);
Chris@16:                     ie = (std::max) (ie, i);
Chris@16:                     ++ itc;
Chris@16:                 }
Chris@16:                 ++ it1;
Chris@16:             }
Chris@16:             for (size_type i = ib; i < ie + 1; ++ i) {
Chris@16:                 if (temporary (i) != value_type/*zero*/()) {
Chris@16:                     // FIXME we'll need to extend the container interface!
Chris@16:                     // m.push_back (i, it2.index2 (), temporary (i));
Chris@16:                     // FIXME What to do with adaptors?
Chris@16:                     // m.insert (i, it2.index2 (), temporary (i));
Chris@16:                     if (triangular_restriction::other (i, it2.index2 ()))
Chris@16:                         m (i, it2.index2 ()) = temporary (i);
Chris@16:                     temporary (i) = value_type/*zero*/();
Chris@16:                 }
Chris@16:             }
Chris@16:             ++ it2;
Chris@16:         }
Chris@16:         return m;
Chris@16:     }
Chris@16: 
Chris@16:     // Dispatcher
Chris@16:     template<class M, class E1, class E2, class TRI>
Chris@16:     BOOST_UBLAS_INLINE
Chris@16:     M &
Chris@16:     sparse_prod (const matrix_expression<E1> &e1,
Chris@16:                  const matrix_expression<E2> &e2,
Chris@16:                  M &m, TRI, bool init = true) {
Chris@16:         typedef typename M::value_type value_type;
Chris@16:         typedef TRI triangular_restriction;
Chris@16:         typedef typename M::orientation_category orientation_category;
Chris@16: 
Chris@16:         if (init)
Chris@16:             m.assign (zero_matrix<value_type> (e1 ().size1 (), e2 ().size2 ()));
Chris@16:         return sparse_prod (e1, e2, m, triangular_restriction (), orientation_category ());
Chris@16:     }
Chris@16:     template<class M, class E1, class E2, class TRI>
Chris@16:     BOOST_UBLAS_INLINE
Chris@16:     M
Chris@16:     sparse_prod (const matrix_expression<E1> &e1,
Chris@16:                  const matrix_expression<E2> &e2,
Chris@16:                  TRI) {
Chris@16:         typedef M matrix_type;
Chris@16:         typedef TRI triangular_restriction;
Chris@16: 
Chris@16:         matrix_type m (e1 ().size1 (), e2 ().size2 ());
Chris@16:         // FIXME needed for c_matrix?!
Chris@16:         // return sparse_prod (e1, e2, m, triangular_restriction (), false);
Chris@16:         return sparse_prod (e1, e2, m, triangular_restriction (), true);
Chris@16:     }
Chris@16:     template<class M, class E1, class E2>
Chris@16:     BOOST_UBLAS_INLINE
Chris@16:     M &
Chris@16:     sparse_prod (const matrix_expression<E1> &e1,
Chris@16:                  const matrix_expression<E2> &e2,
Chris@16:                  M &m, bool init = true) {
Chris@16:         typedef typename M::value_type value_type;
Chris@16:         typedef typename M::orientation_category orientation_category;
Chris@16: 
Chris@16:         if (init)
Chris@16:             m.assign (zero_matrix<value_type> (e1 ().size1 (), e2 ().size2 ()));
Chris@16:         return sparse_prod (e1, e2, m, full (), orientation_category ());
Chris@16:     }
Chris@16:     template<class M, class E1, class E2>
Chris@16:     BOOST_UBLAS_INLINE
Chris@16:     M
Chris@16:     sparse_prod (const matrix_expression<E1> &e1,
Chris@16:                  const matrix_expression<E2> &e2) {
Chris@16:         typedef M matrix_type;
Chris@16: 
Chris@16:         matrix_type m (e1 ().size1 (), e2 ().size2 ());
Chris@16:         // FIXME needed for c_matrix?!
Chris@16:         // return sparse_prod (e1, e2, m, full (), false);
Chris@16:         return sparse_prod (e1, e2, m, full (), true);
Chris@16:     }
Chris@16: 
Chris@16: }}}
Chris@16: 
Chris@16: #endif