]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // |
2 | // Copyright (c) 2009 | |
3 | // Gunter Winkler | |
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 | // | |
10 | ||
11 | #ifndef _BOOST_UBLAS_SPARSE_VIEW_ | |
12 | #define _BOOST_UBLAS_SPARSE_VIEW_ | |
13 | ||
14 | #include <boost/numeric/ublas/matrix_expression.hpp> | |
15 | #include <boost/numeric/ublas/detail/matrix_assign.hpp> | |
16 | #if BOOST_UBLAS_TYPE_CHECK | |
17 | #include <boost/numeric/ublas/matrix.hpp> | |
18 | #endif | |
19 | ||
20 | #include <boost/next_prior.hpp> | |
21 | #include <boost/type_traits/remove_cv.hpp> | |
22 | #include <boost/numeric/ublas/storage.hpp> | |
23 | ||
24 | namespace boost { namespace numeric { namespace ublas { | |
25 | ||
26 | // view a chunk of memory as ublas array | |
27 | ||
28 | template < class T > | |
29 | class c_array_view | |
30 | : public storage_array< c_array_view<T> > { | |
31 | private: | |
32 | typedef c_array_view<T> self_type; | |
33 | typedef T * pointer; | |
34 | ||
35 | public: | |
36 | // TODO: think about a const pointer | |
37 | typedef const pointer array_type; | |
38 | ||
39 | typedef std::size_t size_type; | |
40 | typedef std::ptrdiff_t difference_type; | |
41 | ||
42 | typedef T value_type; | |
43 | typedef const T &const_reference; | |
44 | typedef const T *const_pointer; | |
45 | ||
46 | typedef const_pointer const_iterator; | |
47 | typedef std::reverse_iterator<const_iterator> const_reverse_iterator; | |
48 | ||
49 | // | |
50 | // typedefs required by vector concept | |
51 | // | |
52 | ||
53 | typedef dense_tag storage_category; | |
54 | typedef const vector_reference<const self_type> const_closure_type; | |
55 | ||
56 | c_array_view(size_type size, array_type data) : | |
57 | size_(size), data_(data) | |
58 | {} | |
59 | ||
60 | ~c_array_view() | |
61 | {} | |
62 | ||
63 | // | |
64 | // immutable methods of container concept | |
65 | // | |
66 | ||
67 | BOOST_UBLAS_INLINE | |
68 | size_type size () const { | |
69 | return size_; | |
70 | } | |
71 | ||
72 | BOOST_UBLAS_INLINE | |
73 | const_reference operator [] (size_type i) const { | |
74 | BOOST_UBLAS_CHECK (i < size_, bad_index ()); | |
75 | return data_ [i]; | |
76 | } | |
77 | ||
78 | BOOST_UBLAS_INLINE | |
79 | const_iterator begin () const { | |
80 | return data_; | |
81 | } | |
82 | BOOST_UBLAS_INLINE | |
83 | const_iterator end () const { | |
84 | return data_ + size_; | |
85 | } | |
86 | ||
87 | BOOST_UBLAS_INLINE | |
88 | const_reverse_iterator rbegin () const { | |
89 | return const_reverse_iterator (end ()); | |
90 | } | |
91 | BOOST_UBLAS_INLINE | |
92 | const_reverse_iterator rend () const { | |
93 | return const_reverse_iterator (begin ()); | |
94 | } | |
95 | ||
96 | private: | |
97 | size_type size_; | |
98 | array_type data_; | |
99 | }; | |
100 | ||
101 | ||
102 | /** \brief Present existing arrays as compressed array based | |
103 | * sparse matrix. | |
104 | * This class provides CRS / CCS storage layout. | |
105 | * | |
106 | * see also http://www.netlib.org/utk/papers/templates/node90.html | |
107 | * | |
108 | * \param L layout type, either row_major or column_major | |
109 | * \param IB index base, use 0 for C indexing and 1 for | |
110 | * FORTRAN indexing of the internal index arrays. This | |
111 | * does not affect the operator()(int,int) where the first | |
112 | * row/column has always index 0. | |
113 | * \param IA index array type, e.g., int[] | |
114 | * \param TA value array type, e.g., double[] | |
115 | */ | |
116 | template<class L, std::size_t IB, class IA, class JA, class TA> | |
117 | class compressed_matrix_view: | |
118 | public matrix_expression<compressed_matrix_view<L, IB, IA, JA, TA> > { | |
119 | ||
120 | public: | |
121 | typedef typename vector_view_traits<TA>::value_type value_type; | |
122 | ||
123 | private: | |
124 | typedef value_type &true_reference; | |
125 | typedef value_type *pointer; | |
126 | typedef const value_type *const_pointer; | |
127 | typedef L layout_type; | |
128 | typedef compressed_matrix_view<L, IB, IA, JA, TA> self_type; | |
129 | ||
130 | public: | |
131 | #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS | |
132 | using matrix_expression<self_type>::operator (); | |
133 | #endif | |
134 | // ISSUE require type consistency check | |
135 | // is_convertable (IA::size_type, TA::size_type) | |
136 | typedef typename boost::remove_cv<typename vector_view_traits<JA>::value_type>::type index_type; | |
137 | // for compatibility, should be removed some day ... | |
138 | typedef index_type size_type; | |
139 | // size_type for the data arrays. | |
140 | typedef typename vector_view_traits<JA>::size_type array_size_type; | |
141 | typedef typename vector_view_traits<JA>::difference_type difference_type; | |
142 | typedef const value_type & const_reference; | |
143 | ||
144 | // do NOT define reference type, because class is read only | |
145 | // typedef value_type & reference; | |
146 | ||
147 | typedef IA rowptr_array_type; | |
148 | typedef JA index_array_type; | |
149 | typedef TA value_array_type; | |
150 | typedef const matrix_reference<const self_type> const_closure_type; | |
151 | typedef matrix_reference<self_type> closure_type; | |
152 | ||
153 | // FIXME: define a corresponding temporary type | |
154 | // typedef compressed_vector<T, IB, IA, TA> vector_temporary_type; | |
155 | ||
156 | // FIXME: define a corresponding temporary type | |
157 | // typedef self_type matrix_temporary_type; | |
158 | ||
159 | typedef sparse_tag storage_category; | |
160 | typedef typename L::orientation_category orientation_category; | |
161 | ||
162 | // | |
163 | // private types for internal use | |
164 | // | |
165 | ||
166 | private: | |
167 | typedef typename vector_view_traits<index_array_type>::const_iterator const_subiterator_type; | |
168 | ||
169 | // | |
170 | // Construction and destruction | |
171 | // | |
172 | private: | |
173 | /// private default constructor because data must be filled by caller | |
174 | BOOST_UBLAS_INLINE | |
175 | compressed_matrix_view () { } | |
176 | ||
177 | public: | |
178 | BOOST_UBLAS_INLINE | |
179 | compressed_matrix_view (index_type n_rows, index_type n_cols, array_size_type nnz | |
180 | , const rowptr_array_type & iptr | |
181 | , const index_array_type & jptr | |
182 | , const value_array_type & values): | |
183 | matrix_expression<self_type> (), | |
184 | size1_ (n_rows), size2_ (n_cols), | |
185 | nnz_ (nnz), | |
186 | index1_data_ (iptr), | |
187 | index2_data_ (jptr), | |
188 | value_data_ (values) { | |
189 | storage_invariants (); | |
190 | } | |
191 | ||
192 | BOOST_UBLAS_INLINE | |
193 | compressed_matrix_view(const compressed_matrix_view& o) : | |
194 | size1_(o.size1_), size2_(o.size2_), | |
195 | nnz_(o.nnz_), | |
196 | index1_data_(o.index1_data_), | |
197 | index2_data_(o.index2_data_), | |
198 | value_data_(o.value_data_) | |
199 | {} | |
200 | ||
201 | // | |
202 | // implement immutable iterator types | |
203 | // | |
204 | ||
205 | class const_iterator1 {}; | |
206 | class const_iterator2 {}; | |
207 | ||
208 | typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1; | |
209 | typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2; | |
210 | ||
211 | // | |
212 | // implement all read only methods for the matrix expression concept | |
213 | // | |
214 | ||
215 | //! return the number of rows | |
216 | index_type size1() const { | |
217 | return size1_; | |
218 | } | |
219 | ||
220 | //! return the number of columns | |
221 | index_type size2() const { | |
222 | return size2_; | |
223 | } | |
224 | ||
225 | //! return value at position (i,j) | |
226 | value_type operator()(index_type i, index_type j) const { | |
227 | const_pointer p = find_element(i,j); | |
228 | if (!p) { | |
229 | return zero_; | |
230 | } else { | |
231 | return *p; | |
232 | } | |
233 | } | |
234 | ||
235 | ||
236 | private: | |
237 | // | |
238 | // private helper functions | |
239 | // | |
240 | ||
241 | const_pointer find_element (index_type i, index_type j) const { | |
242 | index_type element1 (layout_type::index_M (i, j)); | |
243 | index_type element2 (layout_type::index_m (i, j)); | |
244 | ||
245 | const array_size_type itv = zero_based( index1_data_[element1] ); | |
246 | const array_size_type itv_next = zero_based( index1_data_[element1+1] ); | |
247 | ||
248 | const_subiterator_type it_start = boost::next(vector_view_traits<index_array_type>::begin(index2_data_),itv); | |
249 | const_subiterator_type it_end = boost::next(vector_view_traits<index_array_type>::begin(index2_data_),itv_next); | |
250 | const_subiterator_type it = find_index_in_row(it_start, it_end, element2) ; | |
251 | ||
252 | if (it == it_end || *it != k_based (element2)) | |
253 | return 0; | |
254 | return &value_data_ [it - vector_view_traits<index_array_type>::begin(index2_data_)]; | |
255 | } | |
256 | ||
257 | const_subiterator_type find_index_in_row(const_subiterator_type it_start | |
258 | , const_subiterator_type it_end | |
259 | , index_type index) const { | |
260 | return std::lower_bound( it_start | |
261 | , it_end | |
262 | , k_based (index) ); | |
263 | } | |
264 | ||
265 | ||
266 | private: | |
267 | void storage_invariants () const { | |
268 | BOOST_UBLAS_CHECK (index1_data_ [layout_type::size_M (size1_, size2_)] == k_based (nnz_), external_logic ()); | |
269 | } | |
270 | ||
271 | index_type size1_; | |
272 | index_type size2_; | |
273 | ||
274 | array_size_type nnz_; | |
275 | ||
276 | const rowptr_array_type & index1_data_; | |
277 | const index_array_type & index2_data_; | |
278 | const value_array_type & value_data_; | |
279 | ||
280 | static const value_type zero_; | |
281 | ||
282 | BOOST_UBLAS_INLINE | |
283 | static index_type zero_based (index_type k_based_index) { | |
284 | return k_based_index - IB; | |
285 | } | |
286 | BOOST_UBLAS_INLINE | |
287 | static index_type k_based (index_type zero_based_index) { | |
288 | return zero_based_index + IB; | |
289 | } | |
290 | ||
291 | friend class iterator1; | |
292 | friend class iterator2; | |
293 | friend class const_iterator1; | |
294 | friend class const_iterator2; | |
295 | }; | |
296 | ||
297 | template<class L, std::size_t IB, class IA, class JA, class TA > | |
298 | const typename compressed_matrix_view<L,IB,IA,JA,TA>::value_type | |
299 | compressed_matrix_view<L,IB,IA,JA,TA>::zero_ = value_type/*zero*/(); | |
300 | ||
301 | ||
302 | template<class L, std::size_t IB, class IA, class JA, class TA > | |
303 | compressed_matrix_view<L,IB,IA,JA,TA> | |
304 | make_compressed_matrix_view(typename vector_view_traits<JA>::value_type n_rows | |
305 | , typename vector_view_traits<JA>::value_type n_cols | |
306 | , typename vector_view_traits<JA>::size_type nnz | |
307 | , const IA & ia | |
308 | , const JA & ja | |
309 | , const TA & ta) { | |
310 | ||
311 | return compressed_matrix_view<L,IB,IA,JA,TA>(n_rows, n_cols, nnz, ia, ja, ta); | |
312 | ||
313 | } | |
314 | ||
315 | }}} | |
316 | ||
317 | #endif |