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 &
sparse_prod(const matrix_expression<E1> & e1,const matrix_expression<E2> & e2,M & m,TRI,row_major_tag)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 &
sparse_prod(const matrix_expression<E1> & e1,const matrix_expression<E2> & e2,M & m,TRI,column_major_tag)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 &
sparse_prod(const matrix_expression<E1> & e1,const matrix_expression<E2> & e2,M & m,TRI,bool init=true)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
sparse_prod(const matrix_expression<E1> & e1,const matrix_expression<E2> & e2,TRI)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 &
sparse_prod(const matrix_expression<E1> & e1,const matrix_expression<E2> & e2,M & m,bool init=true)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
sparse_prod(const matrix_expression<E1> & e1,const matrix_expression<E2> & e2)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
199