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 & |
---|
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 |
---|