...one of the most highly
regarded and expertly designed C++ library projects in the
world.
— Herb Sutter and Andrei
Alexandrescu, C++
Coding Standards
00001 // 00002 // Copyright (c) 2000-2002 00003 // Joerg Walter, Mathias Koch 00004 // 00005 // Distributed under the Boost Software License, Version 1.0. (See 00006 // accompanying file LICENSE_1_0.txt or copy at 00007 // http://www.boost.org/LICENSE_1_0.txt) 00008 // 00009 // The authors gratefully acknowledge the support of 00010 // GeNeSys mbH & Co. KG in producing this work. 00011 // 00012 00013 #ifndef _BOOST_UBLAS_OPERATION_SPARSE_ 00014 #define _BOOST_UBLAS_OPERATION_SPARSE_ 00015 00016 #include <boost/numeric/ublas/traits.hpp> 00017 00018 // These scaled additions were borrowed from MTL unashamedly. 00019 // But Alexei Novakov had a lot of ideas to improve these. Thanks. 00020 00021 namespace boost { namespace numeric { namespace ublas { 00022 00023 template<class M, class E1, class E2, class TRI> 00024 BOOST_UBLAS_INLINE 00025 M & 00026 sparse_prod (const matrix_expression<E1> &e1, 00027 const matrix_expression<E2> &e2, 00028 M &m, TRI, 00029 row_major_tag) { 00030 typedef M matrix_type; 00031 typedef TRI triangular_restriction; 00032 typedef const E1 expression1_type; 00033 typedef const E2 expression2_type; 00034 typedef typename M::size_type size_type; 00035 typedef typename M::value_type value_type; 00036 00037 // ISSUE why is there a dense vector here? 00038 vector<value_type> temporary (e2 ().size2 ()); 00039 temporary.clear (); 00040 typename expression1_type::const_iterator1 it1 (e1 ().begin1 ()); 00041 typename expression1_type::const_iterator1 it1_end (e1 ().end1 ()); 00042 while (it1 != it1_end) { 00043 size_type jb (temporary.size ()); 00044 size_type je (0); 00045 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 00046 typename expression1_type::const_iterator2 it2 (it1.begin ()); 00047 typename expression1_type::const_iterator2 it2_end (it1.end ()); 00048 #else 00049 typename expression1_type::const_iterator2 it2 (boost::numeric::ublas::begin (it1, iterator1_tag ())); 00050 typename expression1_type::const_iterator2 it2_end (boost::numeric::ublas::end (it1, iterator1_tag ())); 00051 #endif 00052 while (it2 != it2_end) { 00053 // temporary.plus_assign (*it2 * row (e2 (), it2.index2 ())); 00054 matrix_row<expression2_type> mr (e2 (), it2.index2 ()); 00055 typename matrix_row<expression2_type>::const_iterator itr (mr.begin ()); 00056 typename matrix_row<expression2_type>::const_iterator itr_end (mr.end ()); 00057 while (itr != itr_end) { 00058 size_type j (itr.index ()); 00059 temporary (j) += *it2 * *itr; 00060 jb = (std::min) (jb, j); 00061 je = (std::max) (je, j); 00062 ++ itr; 00063 } 00064 ++ it2; 00065 } 00066 for (size_type j = jb; j < je + 1; ++ j) { 00067 if (temporary (j) != value_type/*zero*/()) { 00068 // FIXME we'll need to extend the container interface! 00069 // m.push_back (it1.index1 (), j, temporary (j)); 00070 // FIXME What to do with adaptors? 00071 // m.insert (it1.index1 (), j, temporary (j)); 00072 if (triangular_restriction::other (it1.index1 (), j)) 00073 m (it1.index1 (), j) = temporary (j); 00074 temporary (j) = value_type/*zero*/(); 00075 } 00076 } 00077 ++ it1; 00078 } 00079 return m; 00080 } 00081 00082 template<class M, class E1, class E2, class TRI> 00083 BOOST_UBLAS_INLINE 00084 M & 00085 sparse_prod (const matrix_expression<E1> &e1, 00086 const matrix_expression<E2> &e2, 00087 M &m, TRI, 00088 column_major_tag) { 00089 typedef M matrix_type; 00090 typedef TRI triangular_restriction; 00091 typedef const E1 expression1_type; 00092 typedef const E2 expression2_type; 00093 typedef typename M::size_type size_type; 00094 typedef typename M::value_type value_type; 00095 00096 // ISSUE why is there a dense vector here? 00097 vector<value_type> temporary (e1 ().size1 ()); 00098 temporary.clear (); 00099 typename expression2_type::const_iterator2 it2 (e2 ().begin2 ()); 00100 typename expression2_type::const_iterator2 it2_end (e2 ().end2 ()); 00101 while (it2 != it2_end) { 00102 size_type ib (temporary.size ()); 00103 size_type ie (0); 00104 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 00105 typename expression2_type::const_iterator1 it1 (it2.begin ()); 00106 typename expression2_type::const_iterator1 it1_end (it2.end ()); 00107 #else 00108 typename expression2_type::const_iterator1 it1 (boost::numeric::ublas::begin (it2, iterator2_tag ())); 00109 typename expression2_type::const_iterator1 it1_end (boost::numeric::ublas::end (it2, iterator2_tag ())); 00110 #endif 00111 while (it1 != it1_end) { 00112 // column (m, it2.index2 ()).plus_assign (*it1 * column (e1 (), it1.index1 ())); 00113 matrix_column<expression1_type> mc (e1 (), it1.index1 ()); 00114 typename matrix_column<expression1_type>::const_iterator itc (mc.begin ()); 00115 typename matrix_column<expression1_type>::const_iterator itc_end (mc.end ()); 00116 while (itc != itc_end) { 00117 size_type i (itc.index ()); 00118 temporary (i) += *it1 * *itc; 00119 ib = (std::min) (ib, i); 00120 ie = (std::max) (ie, i); 00121 ++ itc; 00122 } 00123 ++ it1; 00124 } 00125 for (size_type i = ib; i < ie + 1; ++ i) { 00126 if (temporary (i) != value_type/*zero*/()) { 00127 // FIXME we'll need to extend the container interface! 00128 // m.push_back (i, it2.index2 (), temporary (i)); 00129 // FIXME What to do with adaptors? 00130 // m.insert (i, it2.index2 (), temporary (i)); 00131 if (triangular_restriction::other (i, it2.index2 ())) 00132 m (i, it2.index2 ()) = temporary (i); 00133 temporary (i) = value_type/*zero*/(); 00134 } 00135 } 00136 ++ it2; 00137 } 00138 return m; 00139 } 00140 00141 // Dispatcher 00142 template<class M, class E1, class E2, class TRI> 00143 BOOST_UBLAS_INLINE 00144 M & 00145 sparse_prod (const matrix_expression<E1> &e1, 00146 const matrix_expression<E2> &e2, 00147 M &m, TRI, bool init = true) { 00148 typedef typename M::value_type value_type; 00149 typedef TRI triangular_restriction; 00150 typedef typename M::orientation_category orientation_category; 00151 00152 if (init) 00153 m.assign (zero_matrix<value_type> (e1 ().size1 (), e2 ().size2 ())); 00154 return sparse_prod (e1, e2, m, triangular_restriction (), orientation_category ()); 00155 } 00156 template<class M, class E1, class E2, class TRI> 00157 BOOST_UBLAS_INLINE 00158 M 00159 sparse_prod (const matrix_expression<E1> &e1, 00160 const matrix_expression<E2> &e2, 00161 TRI) { 00162 typedef M matrix_type; 00163 typedef TRI triangular_restriction; 00164 00165 matrix_type m (e1 ().size1 (), e2 ().size2 ()); 00166 // FIXME needed for c_matrix?! 00167 // return sparse_prod (e1, e2, m, triangular_restriction (), false); 00168 return sparse_prod (e1, e2, m, triangular_restriction (), true); 00169 } 00170 template<class M, class E1, class E2> 00171 BOOST_UBLAS_INLINE 00172 M & 00173 sparse_prod (const matrix_expression<E1> &e1, 00174 const matrix_expression<E2> &e2, 00175 M &m, bool init = true) { 00176 typedef typename M::value_type value_type; 00177 typedef typename M::orientation_category orientation_category; 00178 00179 if (init) 00180 m.assign (zero_matrix<value_type> (e1 ().size1 (), e2 ().size2 ())); 00181 return sparse_prod (e1, e2, m, full (), orientation_category ()); 00182 } 00183 template<class M, class E1, class E2> 00184 BOOST_UBLAS_INLINE 00185 M 00186 sparse_prod (const matrix_expression<E1> &e1, 00187 const matrix_expression<E2> &e2) { 00188 typedef M matrix_type; 00189 00190 matrix_type m (e1 ().size1 (), e2 ().size2 ()); 00191 // FIXME needed for c_matrix?! 00192 // return sparse_prod (e1, e2, m, full (), false); 00193 return sparse_prod (e1, e2, m, full (), true); 00194 } 00195 00196 }}} 00197 00198 #endif