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 & sparse_prod(const matrix_expression<E1> & e1,const matrix_expression<E2> & e2,M & m,TRI,row_major_tag)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 & sparse_prod(const matrix_expression<E1> & e1,const matrix_expression<E2> & e2,M & m,TRI,column_major_tag)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 & sparse_prod(const matrix_expression<E1> & e1,const matrix_expression<E2> & e2,M & m,TRI,bool init=true)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 sparse_prod(const matrix_expression<E1> & e1,const matrix_expression<E2> & e2,TRI)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 & sparse_prod(const matrix_expression<E1> & e1,const matrix_expression<E2> & e2,M & m,bool init=true)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 sparse_prod(const matrix_expression<E1> & e1,const matrix_expression<E2> & e2)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 199