← Back to branch summary

~ubuntu-branches/debian/sid/boost1.49/sid

« back to all changes in this revision

Viewing changes to boost/numeric/ublas/operation_sparse.hpp

  • Committer: Package Import Robot
  • Author(s): Steve M. Robbins
  • Date: 2012-02-26 00:31:44 UTC
  • Revision ID: package-import@ubuntu.com-20120226003144-eaytp12cbf6ubpms
Tags: upstream-1.49.0
Import upstream version 1.49.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
boost/numeric/ublas/operation_sparse.hpp
 
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