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
21namespace boost { namespace numeric { namespace ublas {
22
23 template<class M, class E1, class E2, class TRI>
24 BOOST_UBLAS_INLINE
25 M &
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 &
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 &
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
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 &
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
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

source code of include/boost/numeric/ublas/operation_sparse.hpp