ublas: operation_sparse.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_SPARSE_
00014 #define _BOOST_UBLAS_OPERATION_SPARSE_
00015 
00016 #include <boost/numeric/ublas/traits.hpp>
00017 
00018 // These scaled additions were borrowed from MTL unashamedly.
00019 // But Alexei Novakov had a lot of ideas to improve these. Thanks.
00020 
00021 namespace boost { namespace numeric { namespace ublas {
00022 
00023     template<class M, class E1, class E2, class TRI>
00024     BOOST_UBLAS_INLINE
00025     M &
00026     sparse_prod (const matrix_expression<E1> &e1,
00027                  const matrix_expression<E2> &e2,
00028                  M &m, TRI,
00029                  row_major_tag) {
00030         typedef M matrix_type;
00031         typedef TRI triangular_restriction;
00032         typedef const E1 expression1_type;
00033         typedef const E2 expression2_type;
00034         typedef typename M::size_type size_type;
00035         typedef typename M::value_type value_type;
00036 
00037         // ISSUE why is there a dense vector here?
00038         vector<value_type> temporary (e2 ().size2 ());
00039         temporary.clear ();
00040         typename expression1_type::const_iterator1 it1 (e1 ().begin1 ());
00041         typename expression1_type::const_iterator1 it1_end (e1 ().end1 ());
00042         while (it1 != it1_end) {
00043             size_type jb (temporary.size ());
00044             size_type je (0);
00045 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
00046             typename expression1_type::const_iterator2 it2 (it1.begin ());
00047             typename expression1_type::const_iterator2 it2_end (it1.end ());
00048 #else
00049             typename expression1_type::const_iterator2 it2 (boost::numeric::ublas::begin (it1, iterator1_tag ()));
00050             typename expression1_type::const_iterator2 it2_end (boost::numeric::ublas::end (it1, iterator1_tag ()));
00051 #endif
00052             while (it2 != it2_end) {
00053                 // temporary.plus_assign (*it2 * row (e2 (), it2.index2 ()));
00054                 matrix_row<expression2_type> mr (e2 (), it2.index2 ());
00055                 typename matrix_row<expression2_type>::const_iterator itr (mr.begin ());
00056                 typename matrix_row<expression2_type>::const_iterator itr_end (mr.end ());
00057                 while (itr != itr_end) {
00058                     size_type j (itr.index ());
00059                     temporary (j) += *it2 * *itr;
00060                     jb = (std::min) (jb, j);
00061                     je = (std::max) (je, j);
00062                     ++ itr;
00063                 }
00064                 ++ it2;
00065             }
00066             for (size_type j = jb; j < je + 1; ++ j) {
00067                 if (temporary (j) != value_type/*zero*/()) {
00068                     // FIXME we'll need to extend the container interface!
00069                     // m.push_back (it1.index1 (), j, temporary (j));
00070                     // FIXME What to do with adaptors?
00071                     // m.insert (it1.index1 (), j, temporary (j));
00072                     if (triangular_restriction::other (it1.index1 (), j))
00073                         m (it1.index1 (), j) = temporary (j);
00074                     temporary (j) = value_type/*zero*/();
00075                 }
00076             }
00077             ++ it1;
00078         }
00079         return m;
00080     }
00081 
00082     template<class M, class E1, class E2, class TRI>
00083     BOOST_UBLAS_INLINE
00084     M &
00085     sparse_prod (const matrix_expression<E1> &e1,
00086                  const matrix_expression<E2> &e2,
00087                  M &m, TRI,
00088                  column_major_tag) {
00089         typedef M matrix_type;
00090         typedef TRI triangular_restriction;
00091         typedef const E1 expression1_type;
00092         typedef const E2 expression2_type;
00093         typedef typename M::size_type size_type;
00094         typedef typename M::value_type value_type;
00095 
00096         // ISSUE why is there a dense vector here?
00097         vector<value_type> temporary (e1 ().size1 ());
00098         temporary.clear ();
00099         typename expression2_type::const_iterator2 it2 (e2 ().begin2 ());
00100         typename expression2_type::const_iterator2 it2_end (e2 ().end2 ());
00101         while (it2 != it2_end) {
00102             size_type ib (temporary.size ());
00103             size_type ie (0);
00104 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
00105             typename expression2_type::const_iterator1 it1 (it2.begin ());
00106             typename expression2_type::const_iterator1 it1_end (it2.end ());
00107 #else
00108             typename expression2_type::const_iterator1 it1 (boost::numeric::ublas::begin (it2, iterator2_tag ()));
00109             typename expression2_type::const_iterator1 it1_end (boost::numeric::ublas::end (it2, iterator2_tag ()));
00110 #endif
00111             while (it1 != it1_end) {
00112                 // column (m, it2.index2 ()).plus_assign (*it1 * column (e1 (), it1.index1 ()));
00113                 matrix_column<expression1_type> mc (e1 (), it1.index1 ());
00114                 typename matrix_column<expression1_type>::const_iterator itc (mc.begin ());
00115                 typename matrix_column<expression1_type>::const_iterator itc_end (mc.end ());
00116                 while (itc != itc_end) {
00117                     size_type i (itc.index ());
00118                     temporary (i) += *it1 * *itc;
00119                     ib = (std::min) (ib, i);
00120                     ie = (std::max) (ie, i);
00121                     ++ itc;
00122                 }
00123                 ++ it1;
00124             }
00125             for (size_type i = ib; i < ie + 1; ++ i) {
00126                 if (temporary (i) != value_type/*zero*/()) {
00127                     // FIXME we'll need to extend the container interface!
00128                     // m.push_back (i, it2.index2 (), temporary (i));
00129                     // FIXME What to do with adaptors?
00130                     // m.insert (i, it2.index2 (), temporary (i));
00131                     if (triangular_restriction::other (i, it2.index2 ()))
00132                         m (i, it2.index2 ()) = temporary (i);
00133                     temporary (i) = value_type/*zero*/();
00134                 }
00135             }
00136             ++ it2;
00137         }
00138         return m;
00139     }
00140 
00141     // Dispatcher
00142     template<class M, class E1, class E2, class TRI>
00143     BOOST_UBLAS_INLINE
00144     M &
00145     sparse_prod (const matrix_expression<E1> &e1,
00146                  const matrix_expression<E2> &e2,
00147                  M &m, TRI, bool init = true) {
00148         typedef typename M::value_type value_type;
00149         typedef TRI triangular_restriction;
00150         typedef typename M::orientation_category orientation_category;
00151 
00152         if (init)
00153             m.assign (zero_matrix<value_type> (e1 ().size1 (), e2 ().size2 ()));
00154         return sparse_prod (e1, e2, m, triangular_restriction (), orientation_category ());
00155     }
00156     template<class M, class E1, class E2, class TRI>
00157     BOOST_UBLAS_INLINE
00158     M
00159     sparse_prod (const matrix_expression<E1> &e1,
00160                  const matrix_expression<E2> &e2,
00161                  TRI) {
00162         typedef M matrix_type;
00163         typedef TRI triangular_restriction;
00164 
00165         matrix_type m (e1 ().size1 (), e2 ().size2 ());
00166         // FIXME needed for c_matrix?!
00167         // return sparse_prod (e1, e2, m, triangular_restriction (), false);
00168         return sparse_prod (e1, e2, m, triangular_restriction (), true);
00169     }
00170     template<class M, class E1, class E2>
00171     BOOST_UBLAS_INLINE
00172     M &
00173     sparse_prod (const matrix_expression<E1> &e1,
00174                  const matrix_expression<E2> &e2,
00175                  M &m, bool init = true) {
00176         typedef typename M::value_type value_type;
00177         typedef typename M::orientation_category orientation_category;
00178 
00179         if (init)
00180             m.assign (zero_matrix<value_type> (e1 ().size1 (), e2 ().size2 ()));
00181         return sparse_prod (e1, e2, m, full (), orientation_category ());
00182     }
00183     template<class M, class E1, class E2>
00184     BOOST_UBLAS_INLINE
00185     M
00186     sparse_prod (const matrix_expression<E1> &e1,
00187                  const matrix_expression<E2> &e2) {
00188         typedef M matrix_type;
00189 
00190         matrix_type m (e1 ().size1 (), e2 ().size2 ());
00191         // FIXME needed for c_matrix?!
00192         // return sparse_prod (e1, e2, m, full (), false);
00193         return sparse_prod (e1, e2, m, full (), true);
00194     }
00195 
00196 }}}
00197 
00198 #endif