]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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 |