00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #ifndef EIGEN_CONSERVATIVESPARSESPARSEPRODUCT_H
00026 #define EIGEN_CONSERVATIVESPARSESPARSEPRODUCT_H
00027
00028 namespace Eigen {
00029
00030 namespace internal {
00031
00032 template<typename Lhs, typename Rhs, typename ResultType>
00033 static void conservative_sparse_sparse_product_impl(const Lhs& lhs, const Rhs& rhs, ResultType& res)
00034 {
00035 typedef typename remove_all<Lhs>::type::Scalar Scalar;
00036 typedef typename remove_all<Lhs>::type::Index Index;
00037
00038
00039 Index rows = lhs.innerSize();
00040 Index cols = rhs.outerSize();
00041 eigen_assert(lhs.outerSize() == rhs.innerSize());
00042
00043 std::vector<bool> mask(rows,false);
00044 Matrix<Scalar,Dynamic,1> values(rows);
00045 Matrix<Index,Dynamic,1> indices(rows);
00046
00047
00048
00049
00050
00051
00052
00053 Index estimated_nnz_prod = lhs.nonZeros() + rhs.nonZeros();
00054
00055 res.setZero();
00056 res.reserve(Index(estimated_nnz_prod));
00057
00058 for (Index j=0; j<cols; ++j)
00059 {
00060
00061 res.startVec(j);
00062 Index nnz = 0;
00063 for (typename Rhs::InnerIterator rhsIt(rhs, j); rhsIt; ++rhsIt)
00064 {
00065 Scalar y = rhsIt.value();
00066 Index k = rhsIt.index();
00067 for (typename Lhs::InnerIterator lhsIt(lhs, k); lhsIt; ++lhsIt)
00068 {
00069 Index i = lhsIt.index();
00070 Scalar x = lhsIt.value();
00071 if(!mask[i])
00072 {
00073 mask[i] = true;
00074 values[i] = x * y;
00075 indices[nnz] = i;
00076 ++nnz;
00077 }
00078 else
00079 values[i] += x * y;
00080 }
00081 }
00082
00083
00084 for(int k=0; k<nnz; ++k)
00085 {
00086 int i = indices[k];
00087 res.insertBackByOuterInnerUnordered(j,i) = values[i];
00088 mask[i] = false;
00089 }
00090
00091 #if 0
00092
00093
00094 int t200 = rows/(log2(200)*1.39);
00095 int t = (rows*100)/139;
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105 if(true)
00106 {
00107 if(nnz>1) std::sort(indices.data(),indices.data()+nnz);
00108 for(int k=0; k<nnz; ++k)
00109 {
00110 int i = indices[k];
00111 res.insertBackByOuterInner(j,i) = values[i];
00112 mask[i] = false;
00113 }
00114 }
00115 else
00116 {
00117
00118 for(int i=0; i<rows; ++i)
00119 {
00120 if(mask[i])
00121 {
00122 mask[i] = false;
00123 res.insertBackByOuterInner(j,i) = values[i];
00124 }
00125 }
00126 }
00127 #endif
00128
00129 }
00130 res.finalize();
00131 }
00132
00133
00134 }
00135
00136 namespace internal {
00137
00138 template<typename Lhs, typename Rhs, typename ResultType,
00139 int LhsStorageOrder = traits<Lhs>::Flags&RowMajorBit,
00140 int RhsStorageOrder = traits<Rhs>::Flags&RowMajorBit,
00141 int ResStorageOrder = traits<ResultType>::Flags&RowMajorBit>
00142 struct conservative_sparse_sparse_product_selector;
00143
00144 template<typename Lhs, typename Rhs, typename ResultType>
00145 struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,ColMajor,ColMajor,ColMajor>
00146 {
00147 typedef typename remove_all<Lhs>::type LhsCleaned;
00148 typedef typename LhsCleaned::Scalar Scalar;
00149
00150 static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
00151 {
00152 typedef SparseMatrix<typename ResultType::Scalar,RowMajor> RowMajorMatrix;
00153 typedef SparseMatrix<typename ResultType::Scalar,ColMajor> ColMajorMatrix;
00154 ColMajorMatrix resCol(lhs.rows(),rhs.cols());
00155 internal::conservative_sparse_sparse_product_impl<Lhs,Rhs,ColMajorMatrix>(lhs, rhs, resCol);
00156
00157 RowMajorMatrix resRow(resCol);
00158 res = resRow;
00159 }
00160 };
00161
00162 template<typename Lhs, typename Rhs, typename ResultType>
00163 struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,ColMajor,ColMajor>
00164 {
00165 static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
00166 {
00167 typedef SparseMatrix<typename ResultType::Scalar,RowMajor> RowMajorMatrix;
00168 RowMajorMatrix rhsRow = rhs;
00169 RowMajorMatrix resRow(lhs.rows(), rhs.cols());
00170 internal::conservative_sparse_sparse_product_impl<RowMajorMatrix,Lhs,RowMajorMatrix>(rhsRow, lhs, resRow);
00171 res = resRow;
00172 }
00173 };
00174
00175 template<typename Lhs, typename Rhs, typename ResultType>
00176 struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,ColMajor,RowMajor,ColMajor>
00177 {
00178 static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
00179 {
00180 typedef SparseMatrix<typename ResultType::Scalar,RowMajor> RowMajorMatrix;
00181 RowMajorMatrix lhsRow = lhs;
00182 RowMajorMatrix resRow(lhs.rows(), rhs.cols());
00183 internal::conservative_sparse_sparse_product_impl<Rhs,RowMajorMatrix,RowMajorMatrix>(rhs, lhsRow, resRow);
00184 res = resRow;
00185 }
00186 };
00187
00188 template<typename Lhs, typename Rhs, typename ResultType>
00189 struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,RowMajor,ColMajor>
00190 {
00191 static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
00192 {
00193 typedef SparseMatrix<typename ResultType::Scalar,RowMajor> RowMajorMatrix;
00194 RowMajorMatrix resRow(lhs.rows(), rhs.cols());
00195 internal::conservative_sparse_sparse_product_impl<Rhs,Lhs,RowMajorMatrix>(rhs, lhs, resRow);
00196 res = resRow;
00197 }
00198 };
00199
00200
00201 template<typename Lhs, typename Rhs, typename ResultType>
00202 struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,ColMajor,ColMajor,RowMajor>
00203 {
00204 typedef typename traits<typename remove_all<Lhs>::type>::Scalar Scalar;
00205
00206 static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
00207 {
00208 typedef SparseMatrix<typename ResultType::Scalar,ColMajor> ColMajorMatrix;
00209 ColMajorMatrix resCol(lhs.rows(), rhs.cols());
00210 internal::conservative_sparse_sparse_product_impl<Lhs,Rhs,ColMajorMatrix>(lhs, rhs, resCol);
00211 res = resCol;
00212 }
00213 };
00214
00215 template<typename Lhs, typename Rhs, typename ResultType>
00216 struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,ColMajor,RowMajor>
00217 {
00218 static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
00219 {
00220 typedef SparseMatrix<typename ResultType::Scalar,ColMajor> ColMajorMatrix;
00221 ColMajorMatrix lhsCol = lhs;
00222 ColMajorMatrix resCol(lhs.rows(), rhs.cols());
00223 internal::conservative_sparse_sparse_product_impl<ColMajorMatrix,Rhs,ColMajorMatrix>(lhsCol, rhs, resCol);
00224 res = resCol;
00225 }
00226 };
00227
00228 template<typename Lhs, typename Rhs, typename ResultType>
00229 struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,ColMajor,RowMajor,RowMajor>
00230 {
00231 static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
00232 {
00233 typedef SparseMatrix<typename ResultType::Scalar,ColMajor> ColMajorMatrix;
00234 ColMajorMatrix rhsCol = rhs;
00235 ColMajorMatrix resCol(lhs.rows(), rhs.cols());
00236 internal::conservative_sparse_sparse_product_impl<Lhs,ColMajorMatrix,ColMajorMatrix>(lhs, rhsCol, resCol);
00237 res = resCol;
00238 }
00239 };
00240
00241 template<typename Lhs, typename Rhs, typename ResultType>
00242 struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,RowMajor,RowMajor>
00243 {
00244 static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
00245 {
00246 typedef SparseMatrix<typename ResultType::Scalar,RowMajor> RowMajorMatrix;
00247 typedef SparseMatrix<typename ResultType::Scalar,ColMajor> ColMajorMatrix;
00248 RowMajorMatrix resRow(lhs.rows(),rhs.cols());
00249 internal::conservative_sparse_sparse_product_impl<Rhs,Lhs,RowMajorMatrix>(rhs, lhs, resRow);
00250
00251 ColMajorMatrix resCol(resRow);
00252 res = resCol;
00253 }
00254 };
00255
00256 }
00257
00258 }
00259
00260 #endif // EIGEN_CONSERVATIVESPARSESPARSEPRODUCT_H