1 //
2 //  Copyright (c) 2000-2002
3 //  Joerg Walter, Mathias Koch
4 //
5 //  Permission to use, copy, modify, distribute and sell this software
6 //  and its documentation for any purpose is hereby granted without fee,
7 //  provided that the above copyright notice appear in all copies and
8 //  that both that copyright notice and this permission notice appear
9 //  in supporting documentation.  The authors make no representations
10 //  about the suitability of this software for any purpose.
11 //  It is provided "as is" without express or implied warranty.
12 //
13 //  The authors gratefully acknowledge the support of
14 //  GeNeSys mbH & Co. KG in producing this work.
15 //
16 
17 #ifndef _BOOST_UBLAS_OPERATION_SPARSE_
18 #define _BOOST_UBLAS_OPERATION_SPARSE_
19 
20 #include <boost/numeric/ublas/traits.hpp>
21 
22 // These scaled additions were borrowed from MTL unashamedly.
23 // But Alexei Novakov had a lot of ideas to improve these. Thanks.
24 
25 namespace boost { namespace numeric { namespace ublas {
26 
27     template<class M, class E1, class E2, class TRI>
28     BOOST_UBLAS_INLINE
29     M &
sparse_prod(const matrix_expression<E1> & e1,const matrix_expression<E2> & e2,M & m,TRI,row_major_tag)30     sparse_prod (const matrix_expression<E1> &e1,
31                  const matrix_expression<E2> &e2,
32                  M &m, TRI,
33                  row_major_tag) {
34         typedef M matrix_type;
35         typedef TRI triangular_restriction;
36         typedef const E1 expression1_type;
37         typedef const E2 expression2_type;
38         typedef typename M::size_type size_type;
39         typedef typename M::value_type value_type;
40 
41         // ISSUE why is there a dense vector here?
42         vector<value_type> temporary (e2 ().size2 ());
43         temporary.clear ();
44 #if BOOST_UBLAS_TYPE_CHECK
45         matrix<value_type, row_major> cm (m.size1 (), m.size2 ());
46         typedef typename type_traits<value_type>::real_type real_type;
47         real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));
48         indexing_matrix_assign<scalar_assign> (cm, prod (e1, e2), row_major_tag ());
49 #endif
50         typename expression1_type::const_iterator1 it1 (e1 ().begin1 ());
51         typename expression1_type::const_iterator1 it1_end (e1 ().end1 ());
52         while (it1 != it1_end) {
53             size_type jb (temporary.size ());
54             size_type je (0);
55 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
56             typename expression1_type::const_iterator2 it2 (it1.begin ());
57             typename expression1_type::const_iterator2 it2_end (it1.end ());
58 #else
59             typename expression1_type::const_iterator2 it2 (boost::numeric::ublas::begin (it1, iterator1_tag ()));
60             typename expression1_type::const_iterator2 it2_end (boost::numeric::ublas::end (it1, iterator1_tag ()));
61 #endif
62             while (it2 != it2_end) {
63                 // temporary.plus_assign (*it2 * row (e2 (), it2.index2 ()));
64                 matrix_row<expression2_type> mr (e2 (), it2.index2 ());
65                 typename matrix_row<expression2_type>::const_iterator itr (mr.begin ());
66                 typename matrix_row<expression2_type>::const_iterator itr_end (mr.end ());
67                 while (itr != itr_end) {
68                     size_type j (itr.index ());
69                     temporary (j) += *it2 * *itr;
70                     jb = (std::min) (jb, j);
71                     je = (std::max) (je, j);
72                     ++ itr;
73                 }
74                 ++ it2;
75             }
76             for (size_type j = jb; j < je + 1; ++ j) {
77                 if (temporary (j) != value_type/*zero*/()) {
78                     // FIXME we'll need to extend the container interface!
79                     // m.push_back (it1.index1 (), j, temporary (j));
80                     // FIXME What to do with adaptors?
81                     // m.insert (it1.index1 (), j, temporary (j));
82                     if (triangular_restriction::other (it1.index1 (), j))
83                         m (it1.index1 (), j) = temporary (j);
84                     temporary (j) = value_type/*zero*/();
85                 }
86             }
87             ++ it1;
88         }
89 #if BOOST_UBLAS_TYPE_CHECK
90         BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());
91 #endif
92         return m;
93     }
94 
95     template<class M, class E1, class E2, class TRI>
96     BOOST_UBLAS_INLINE
97     M &
sparse_prod(const matrix_expression<E1> & e1,const matrix_expression<E2> & e2,M & m,TRI,column_major_tag)98     sparse_prod (const matrix_expression<E1> &e1,
99                  const matrix_expression<E2> &e2,
100                  M &m, TRI,
101                  column_major_tag) {
102         typedef M matrix_type;
103         typedef TRI triangular_restriction;
104         typedef const E1 expression1_type;
105         typedef const E2 expression2_type;
106         typedef typename M::size_type size_type;
107         typedef typename M::value_type value_type;
108 
109         // ISSUE why is there a dense vector here?
110         vector<value_type> temporary (e1 ().size1 ());
111         temporary.clear ();
112 #if BOOST_UBLAS_TYPE_CHECK
113         matrix<value_type, column_major> cm (m.size1 (), m.size2 ());
114         typedef typename type_traits<value_type>::real_type real_type;
115         real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));
116         indexing_matrix_assign<scalar_assign> (cm, prod (e1, e2), column_major_tag ());
117 #endif
118         typename expression2_type::const_iterator2 it2 (e2 ().begin2 ());
119         typename expression2_type::const_iterator2 it2_end (e2 ().end2 ());
120         while (it2 != it2_end) {
121             size_type ib (temporary.size ());
122             size_type ie (0);
123 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
124             typename expression2_type::const_iterator1 it1 (it2.begin ());
125             typename expression2_type::const_iterator1 it1_end (it2.end ());
126 #else
127             typename expression2_type::const_iterator1 it1 (boost::numeric::ublas::begin (it2, iterator2_tag ()));
128             typename expression2_type::const_iterator1 it1_end (boost::numeric::ublas::end (it2, iterator2_tag ()));
129 #endif
130             while (it1 != it1_end) {
131                 // column (m, it2.index2 ()).plus_assign (*it1 * column (e1 (), it1.index1 ()));
132                 matrix_column<expression1_type> mc (e1 (), it1.index1 ());
133                 typename matrix_column<expression1_type>::const_iterator itc (mc.begin ());
134                 typename matrix_column<expression1_type>::const_iterator itc_end (mc.end ());
135                 while (itc != itc_end) {
136                     size_type i (itc.index ());
137                     temporary (i) += *it1 * *itc;
138                     ib = (std::min) (ib, i);
139                     ie = (std::max) (ie, i);
140                     ++ itc;
141                 }
142                 ++ it1;
143             }
144             for (size_type i = ib; i < ie + 1; ++ i) {
145                 if (temporary (i) != value_type/*zero*/()) {
146                     // FIXME we'll need to extend the container interface!
147                     // m.push_back (i, it2.index2 (), temporary (i));
148                     // FIXME What to do with adaptors?
149                     // m.insert (i, it2.index2 (), temporary (i));
150                     if (triangular_restriction::other (i, it2.index2 ()))
151                         m (i, it2.index2 ()) = temporary (i);
152                     temporary (i) = value_type/*zero*/();
153                 }
154             }
155             ++ it2;
156         }
157 #if BOOST_UBLAS_TYPE_CHECK
158         BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());
159 #endif
160         return m;
161     }
162 
163     // Dispatcher
164     template<class M, class E1, class E2, class TRI>
165     BOOST_UBLAS_INLINE
166     M &
sparse_prod(const matrix_expression<E1> & e1,const matrix_expression<E2> & e2,M & m,TRI,bool init=true)167     sparse_prod (const matrix_expression<E1> &e1,
168                  const matrix_expression<E2> &e2,
169                  M &m, TRI, bool init = true) {
170         typedef typename M::value_type value_type;
171         typedef TRI triangular_restriction;
172         typedef typename M::orientation_category orientation_category;
173 
174         if (init)
175             m.assign (zero_matrix<value_type> (e1 ().size1 (), e2 ().size2 ()));
176         return sparse_prod (e1, e2, m, triangular_restriction (), orientation_category ());
177     }
178     template<class M, class E1, class E2, class TRI>
179     BOOST_UBLAS_INLINE
180     M
sparse_prod(const matrix_expression<E1> & e1,const matrix_expression<E2> & e2,TRI)181     sparse_prod (const matrix_expression<E1> &e1,
182                  const matrix_expression<E2> &e2,
183                  TRI) {
184         typedef M matrix_type;
185         typedef TRI triangular_restriction;
186 
187         matrix_type m (e1 ().size1 (), e2 ().size2 ());
188         // FIXME needed for c_matrix?!
189         // return sparse_prod (e1, e2, m, triangular_restriction (), false);
190         return sparse_prod (e1, e2, m, triangular_restriction (), true);
191     }
192     template<class M, class E1, class E2>
193     BOOST_UBLAS_INLINE
194     M &
sparse_prod(const matrix_expression<E1> & e1,const matrix_expression<E2> & e2,M & m,bool init=true)195     sparse_prod (const matrix_expression<E1> &e1,
196                  const matrix_expression<E2> &e2,
197                  M &m, bool init = true) {
198         typedef typename M::value_type value_type;
199         typedef typename M::orientation_category orientation_category;
200 
201         if (init)
202             m.assign (zero_matrix<value_type> (e1 ().size1 (), e2 ().size2 ()));
203         return sparse_prod (e1, e2, m, full (), orientation_category ());
204     }
205     template<class M, class E1, class E2>
206     BOOST_UBLAS_INLINE
207     M
sparse_prod(const matrix_expression<E1> & e1,const matrix_expression<E2> & e2)208     sparse_prod (const matrix_expression<E1> &e1,
209                  const matrix_expression<E2> &e2) {
210         typedef M matrix_type;
211 
212         matrix_type m (e1 ().size1 (), e2 ().size2 ());
213         // FIXME needed for c_matrix?!
214         // return sparse_prod (e1, e2, m, full (), false);
215         return sparse_prod (e1, e2, m, full (), true);
216     }
217 
218 }}}
219 
220 #endif
221