- Main Page
- Related Pages
- Classes
- Files

- File List

# operation_sparse.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_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