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 25 namespace boost { namespace numeric { namespace ublas { 26 27 template<class M, class E1, class E2, class TRI> 28 BOOST_UBLAS_INLINE 29 M & sparse_prod(const matrix_expression<E1> & e1,const matrix_expression<E2> & e2,M & m,TRI,row_major_tag)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 & sparse_prod(const matrix_expression<E1> & e1,const matrix_expression<E2> & e2,M & m,TRI,column_major_tag)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 & sparse_prod(const matrix_expression<E1> & e1,const matrix_expression<E2> & e2,M & m,TRI,bool init=true)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 sparse_prod(const matrix_expression<E1> & e1,const matrix_expression<E2> & e2,TRI)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 & sparse_prod(const matrix_expression<E1> & e1,const matrix_expression<E2> & e2,M & m,bool init=true)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 sparse_prod(const matrix_expression<E1> & e1,const matrix_expression<E2> & e2)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 221