1
2 #include <Eigen/Sparse>
3 #include <bench/BenchTimer.h>
4 #include <set>
5
6 using namespace std;
7 using namespace Eigen;
8 using namespace Eigen;
9
10 #ifndef SIZE
11 #define SIZE 1024
12 #endif
13
14 #ifndef DENSITY
15 #define DENSITY 0.01
16 #endif
17
18 #ifndef SCALAR
19 #define SCALAR double
20 #endif
21
22 typedef SCALAR Scalar;
23 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
24 typedef Matrix<Scalar,Dynamic,1> DenseVector;
25 typedef SparseMatrix<Scalar> EigenSparseMatrix;
26
fillMatrix(float density,int rows,int cols,EigenSparseMatrix & dst)27 void fillMatrix(float density, int rows, int cols, EigenSparseMatrix& dst)
28 {
29 dst.reserve(double(rows)*cols*density);
30 for(int j = 0; j < cols; j++)
31 {
32 for(int i = 0; i < rows; i++)
33 {
34 Scalar v = (internal::random<float>(0,1) < density) ? internal::random<Scalar>() : 0;
35 if (v!=0)
36 dst.insert(i,j) = v;
37 }
38 }
39 dst.finalize();
40 }
41
fillMatrix2(int nnzPerCol,int rows,int cols,EigenSparseMatrix & dst)42 void fillMatrix2(int nnzPerCol, int rows, int cols, EigenSparseMatrix& dst)
43 {
44 // std::cout << "alloc " << nnzPerCol*cols << "\n";
45 dst.reserve(nnzPerCol*cols);
46 for(int j = 0; j < cols; j++)
47 {
48 std::set<int> aux;
49 for(int i = 0; i < nnzPerCol; i++)
50 {
51 int k = internal::random<int>(0,rows-1);
52 while (aux.find(k)!=aux.end())
53 k = internal::random<int>(0,rows-1);
54 aux.insert(k);
55
56 dst.insert(k,j) = internal::random<Scalar>();
57 }
58 }
59 dst.finalize();
60 }
61
eiToDense(const EigenSparseMatrix & src,DenseMatrix & dst)62 void eiToDense(const EigenSparseMatrix& src, DenseMatrix& dst)
63 {
64 dst.setZero();
65 for (int j=0; j<src.cols(); ++j)
66 for (EigenSparseMatrix::InnerIterator it(src.derived(), j); it; ++it)
67 dst(it.index(),j) = it.value();
68 }
69
70 #ifndef NOGMM
71 #include "gmm/gmm.h"
72 typedef gmm::csc_matrix<Scalar> GmmSparse;
73 typedef gmm::col_matrix< gmm::wsvector<Scalar> > GmmDynSparse;
eiToGmm(const EigenSparseMatrix & src,GmmSparse & dst)74 void eiToGmm(const EigenSparseMatrix& src, GmmSparse& dst)
75 {
76 GmmDynSparse tmp(src.rows(), src.cols());
77 for (int j=0; j<src.cols(); ++j)
78 for (EigenSparseMatrix::InnerIterator it(src.derived(), j); it; ++it)
79 tmp(it.index(),j) = it.value();
80 gmm::copy(tmp, dst);
81 }
82 #endif
83
84 #ifndef NOMTL
85 #include <boost/numeric/mtl/mtl.hpp>
86 typedef mtl::compressed2D<Scalar, mtl::matrix::parameters<mtl::tag::col_major> > MtlSparse;
87 typedef mtl::compressed2D<Scalar, mtl::matrix::parameters<mtl::tag::row_major> > MtlSparseRowMajor;
eiToMtl(const EigenSparseMatrix & src,MtlSparse & dst)88 void eiToMtl(const EigenSparseMatrix& src, MtlSparse& dst)
89 {
90 mtl::matrix::inserter<MtlSparse> ins(dst);
91 for (int j=0; j<src.cols(); ++j)
92 for (EigenSparseMatrix::InnerIterator it(src.derived(), j); it; ++it)
93 ins[it.index()][j] = it.value();
94 }
95 #endif
96
97 #ifdef CSPARSE
98 extern "C" {
99 #include "cs.h"
100 }
eiToCSparse(const EigenSparseMatrix & src,cs * & dst)101 void eiToCSparse(const EigenSparseMatrix& src, cs* &dst)
102 {
103 cs* aux = cs_spalloc (0, 0, 1, 1, 1);
104 for (int j=0; j<src.cols(); ++j)
105 for (EigenSparseMatrix::InnerIterator it(src.derived(), j); it; ++it)
106 if (!cs_entry(aux, it.index(), j, it.value()))
107 {
108 std::cout << "cs_entry error\n";
109 exit(2);
110 }
111 dst = cs_compress(aux);
112 // cs_spfree(aux);
113 }
114 #endif // CSPARSE
115
116 #ifndef NOUBLAS
117 #include <boost/numeric/ublas/vector.hpp>
118 #include <boost/numeric/ublas/matrix.hpp>
119 #include <boost/numeric/ublas/io.hpp>
120 #include <boost/numeric/ublas/triangular.hpp>
121 #include <boost/numeric/ublas/vector_sparse.hpp>
122 #include <boost/numeric/ublas/matrix_sparse.hpp>
123 #include <boost/numeric/ublas/vector_of_vector.hpp>
124 #include <boost/numeric/ublas/operation.hpp>
125
126 typedef boost::numeric::ublas::compressed_matrix<Scalar,boost::numeric::ublas::column_major> UBlasSparse;
127
eiToUblas(const EigenSparseMatrix & src,UBlasSparse & dst)128 void eiToUblas(const EigenSparseMatrix& src, UBlasSparse& dst)
129 {
130 dst.resize(src.rows(), src.cols(), false);
131 for (int j=0; j<src.cols(); ++j)
132 for (EigenSparseMatrix::InnerIterator it(src.derived(), j); it; ++it)
133 dst(it.index(),j) = it.value();
134 }
135
136 template <typename EigenType, typename UblasType>
eiToUblasVec(const EigenType & src,UblasType & dst)137 void eiToUblasVec(const EigenType& src, UblasType& dst)
138 {
139 dst.resize(src.size());
140 for (int j=0; j<src.size(); ++j)
141 dst[j] = src.coeff(j);
142 }
143 #endif
144
145 #ifdef OSKI
146 extern "C" {
147 #include <oski/oski.h>
148 }
149 #endif
150