source: NonGTP/Boost/boost/numeric/ublas/operation_sparse.hpp @ 857

Revision 857, 9.4 KB checked in by igarcia, 19 years ago (diff)
Line 
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
25namespace boost { namespace numeric { namespace ublas {
26
27    template<class M, class E1, class E2, class TRI>
28    BOOST_UBLAS_INLINE
29    M &
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 &
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 &
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
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 &
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
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
Note: See TracBrowser for help on using the repository browser.