diff --git a/Eigen/Sparse b/Eigen/Sparse index 2364c6fe4..54b15b8ae 100644 --- a/Eigen/Sparse +++ b/Eigen/Sparse @@ -8,6 +8,10 @@ #include #include +#ifdef EIGEN_GOOGLEHASH_SUPPORT + #include +#endif + #ifdef EIGEN_CHOLMOD_SUPPORT extern "C" { #include "cholmod.h" @@ -70,10 +74,10 @@ namespace Eigen { #include "src/Sparse/RandomSetter.h" #include "src/Sparse/SparseBlock.h" #include "src/Sparse/SparseMatrix.h" -#include "src/Sparse/HashMatrix.h" -#include "src/Sparse/LinkedVectorMatrix.h" +//#include "src/Sparse/HashMatrix.h" +//#include "src/Sparse/LinkedVectorMatrix.h" #include "src/Sparse/CoreIterators.h" -#include "src/Sparse/SparseSetter.h" +//#include "src/Sparse/SparseSetter.h" #include "src/Sparse/SparseProduct.h" #include "src/Sparse/TriangularSolver.h" #include "src/Sparse/SparseLLT.h" diff --git a/Eigen/src/Sparse/RandomSetter.h b/Eigen/src/Sparse/RandomSetter.h index 8606e03c6..35bc9daee 100644 --- a/Eigen/src/Sparse/RandomSetter.h +++ b/Eigen/src/Sparse/RandomSetter.h @@ -29,6 +29,9 @@ template struct StdMapTraits { typedef int KeyType; typedef std::map Type; + enum { + IsSorted = 1 + }; static void setInvalidKey(Type&, const KeyType&) {} }; @@ -38,6 +41,9 @@ template struct GnuHashMapTraits { typedef int KeyType; typedef __gnu_cxx::hash_map Type; + enum { + IsSorted = 0 + }; static void setInvalidKey(Type&, const KeyType&) {} }; @@ -48,6 +54,9 @@ template struct GoogleDenseHashMapTraits { typedef int KeyType; typedef google::dense_hash_map Type; + enum { + IsSorted = 0 + }; static void setInvalidKey(Type& map, const KeyType& k) { map.set_empty_key(k); } @@ -59,6 +68,9 @@ template struct GoogleSparseHashMapTraits { typedef int KeyType; typedef google::sparse_hash_map Type; + enum { + IsSorted = 0 + }; static void setInvalidKey(Type&, const KeyType&) {} }; @@ -66,10 +78,35 @@ template struct GoogleSparseHashMapTraits /** \class RandomSetter * + * Typical usage: + * \code + * SparseMatrix m(rows,cols); + * { + * RandomSetter > w(m); + * // don't use m but w instead with read/write random access to the coefficients: + * for(;;) + * w(rand(),rand()) = rand; + * } + * // when w is deleted, the data are copied back to m + * // and m is ready to use. + * \endcode + * + * \note for performance and memory consumption reasons it is highly recommended to use + * Google's hash library. To do so you have two options: + * - include yourself \b before Eigen/Sparse header + * - define EIGEN_GOOGLEHASH_SUPPORT + * In the later case the inclusion of is made for you. */ template class HashMapTraits = StdMapTraits, - int OuterPacketBits = 6> + template class MapTraits = +#if defined _DENSE_HASH_MAP_H_ + GoogleDenseHashMapTraits +#elif defined _HASH_MAP + GnuHashMapTraits +#else + StdMapTraits +#endif + ,int OuterPacketBits = 6> class RandomSetter { typedef typename ei_traits::Scalar Scalar; @@ -78,11 +115,13 @@ class RandomSetter ScalarWrapper() : value(0) {} Scalar value; }; - typedef typename HashMapTraits::KeyType KeyType; - typedef typename HashMapTraits::Type HashMapType; + typedef typename MapTraits::KeyType KeyType; + typedef typename MapTraits::Type HashMapType; static const int OuterPacketMask = (1 << OuterPacketBits) - 1; enum { - RowMajor = SparseMatrixType::Flags & RowMajorBit + SwapStorage = 1 - MapTraits::IsSorted, + TargetRowMajor = (SparseMatrixType::Flags & RowMajorBit) ? 1 : 0, + SetterRowMajor = SwapStorage ? 1-TargetRowMajor : TargetRowMajor }; public: @@ -90,31 +129,114 @@ class RandomSetter inline RandomSetter(SparseMatrixType& target) : mp_target(&target) { - m_outerPackets = target.outerSize() >> OuterPacketBits; - if (target.outerSize()&OuterPacketMask) + const int outerSize = SwapStorage ? target.innerSize() : target.outerSize(); + const int innerSize = SwapStorage ? target.outerSize() : target.innerSize(); + m_outerPackets = outerSize >> OuterPacketBits; + if (outerSize&OuterPacketMask) m_outerPackets += 1; m_hashmaps = new HashMapType[m_outerPackets]; - KeyType ik = (1<innerSize()+1; + // compute number of bits needed to store inner indices + int aux = innerSize - 1; + m_keyBitsOffset = 0; + while (aux) + { + m_keyBitsOffset++; + aux = aux >> 1; + } + KeyType ik = (1<<(OuterPacketBits+m_keyBitsOffset)); for (int k=0; k::setInvalidKey(m_hashmaps[k],ik); + MapTraits::setInvalidKey(m_hashmaps[k],ik); + + // insert current coeffs + for (int j=0; jouterSize(); ++j) + for (typename SparseMatrixType::InnerIterator it(*mp_target,j); it; ++it) + (*this)(TargetRowMajor?j:it.index(), TargetRowMajor?it.index():j) = it.value(); } ~RandomSetter() { + KeyType keyBitsMask = (1<startFill(nonZeros()); + for (int k=0; kfirst >> m_keyBitsOffset) + outerOffset; + const int inner = it->first & keyBitsMask; + mp_target->fill(TargetRowMajor ? outer : inner, TargetRowMajor ? inner : outer) = it->second.value; + } + } + mp_target->endFill(); + } + else + { + VectorXi positions(mp_target->outerSize()); + positions.setZero(); + // pass 1 + for (int k=0; kfirst & keyBitsMask; + positions[outer]++; + } + } + // prefix sum + int count = 0; + for (int j=0; jouterSize(); ++j) + { + int tmp = positions[j]; + mp_target->_outerIndexPtr()[j] = count; + positions[j] = count; + count += tmp; + } + mp_target->_outerIndexPtr()[mp_target->outerSize()] = count; + mp_target->resizeNonZeros(count); + // pass 2 + for (int k=0; kfirst >> m_keyBitsOffset) + outerOffset; + const int outer = it->first & keyBitsMask; + // sorted insertion + // Note that we have to deal with at most 2^OuterPacketBits unsorted coefficients, + // moreover those 2^OuterPacketBits coeffs are likely to be sparse, an so only a + // small fraction of them have to be sorted, whence the following simple procedure: + int posStart = mp_target->_outerIndexPtr()[outer]; + int i = (positions[outer]++) - 1; + while ( (i >= posStart) && (mp_target->_innerIndexPtr()[i] > inner) ) + { + mp_target->_valuePtr()[i+1] = mp_target->_valuePtr()[i]; + mp_target->_innerIndexPtr()[i+1] = mp_target->_innerIndexPtr()[i]; + --i; + } + mp_target->_innerIndexPtr()[i+1] = inner; + mp_target->_valuePtr()[i+1] = it->second.value; + } + } + } delete[] m_hashmaps; } Scalar& operator() (int row, int col) { - const int outer = RowMajor ? row : col; - const int inner = RowMajor ? col : row; - const int outerMajor = outer >> OuterPacketBits; - const int outerMinor = outer & OuterPacketMask; - const KeyType key = inner + outerMinor * mp_target->innerSize(); - + const int outer = SetterRowMajor ? row : col; + const int inner = SetterRowMajor ? col : row; + const int outerMajor = outer >> OuterPacketBits; // index of the packet/map + const int outerMinor = outer & OuterPacketMask; // index of the inner vector in the packet + const KeyType key = (KeyType(outerMinor)< { // std::cout << "Derived& operator=(const MatrixBase& other)\n"; const bool transpose = (Flags & RowMajorBit) != (OtherDerived::Flags & RowMajorBit); -// std::cout << "eval transpose = " << transpose << "\n"; + ei_assert((!transpose) && "the transpose operation is supposed to be handled in SparseMatrix::operator="); const int outerSize = other.outerSize(); - typedef typename ei_meta_if, Derived>::ret TempType; - TempType temp(other.rows(), other.cols()); + //typedef typename ei_meta_if, Derived>::ret TempType; + // thanks to shallow copies, we always eval to a tempary + Derived temp(other.rows(), other.cols()); temp.startFill(std::max(this->rows(),this->cols())*2); for (int j=0; j // g++ -DSIZE=10000 -DDENSITY=0.001 sparse_cholesky.cpp -I.. -DDENSEMATRI -O3 -g0 -DNDEBUG -DNBTRIES=1 -I /home/gael/Coding/LinearAlgebra/taucs_full/src/ -I/home/gael/Coding/LinearAlgebra/taucs_full/build/linux/ -L/home/gael/Coding/LinearAlgebra/taucs_full/lib/linux/ -ltaucs /home/gael/Coding/LinearAlgebra/GotoBLAS/libgoto.a -lpthread -I /home/gael/Coding/LinearAlgebra/SuiteSparse/CHOLMOD/Include/ $CHOLLIB -I /home/gael/Coding/LinearAlgebra/SuiteSparse/UFconfig/ /home/gael/Coding/LinearAlgebra/SuiteSparse/CCOLAMD/Lib/libccolamd.a /home/gael/Coding/LinearAlgebra/SuiteSparse/CHOLMOD/Lib/libcholmod.a -lmetis /home/gael/Coding/LinearAlgebra/SuiteSparse/AMD/Lib/libamd.a /home/gael/Coding/LinearAlgebra/SuiteSparse/CAMD/Lib/libcamd.a /home/gael/Coding/LinearAlgebra/SuiteSparse/CCOLAMD/Lib/libccolamd.a /home/gael/Coding/LinearAlgebra/SuiteSparse/COLAMD/Lib/libcolamd.a -llapack && ./a.out @@ -70,7 +70,7 @@ void doEigen(const char* name, const EigenSparseSelfAdjointMatrix& sm1, int flag std::cout << ":\t" << timer.value() << endl; std::cout << " nnz: " << sm1.nonZeros() << " => " << chol.matrixL().nonZeros() << "\n"; - //std::cout << "sparse\n" << chol.matrixL() << "%\n"; + std::cout << "sparse\n" << chol.matrixL() << "%\n"; } int main(int argc, char *argv[]) diff --git a/bench/sparse_randomsetter.cpp b/bench/sparse_randomsetter.cpp index 7868e177a..61753d8c2 100644 --- a/bench/sparse_randomsetter.cpp +++ b/bench/sparse_randomsetter.cpp @@ -37,6 +37,31 @@ X \ } timer.stop(); } + +static double rtime; +static double nentries; + +template +void dostuff(const char* name, EigenSparseMatrix& sm1) +{ + int rows = sm1.rows(); + int cols = sm1.cols(); + sm1.setZero(); + BenchTimer t; + SetterType* set1 = new SetterType(sm1); + t.reset(); t.start(); + for (int k=0; k(0,rows-1),ei_random(0,cols-1)) += 1; + t.stop(); + std::cout << "std::map => \t" << t.value()-rtime + << " nnz=" << set1->nonZeros() << std::flush; + + // getchar(); + + t.reset(); t.start(); delete set1; t.stop(); + std::cout << " back: \t" << t.value() << "\n"; +} + int main(int argc, char *argv[]) { int rows = SIZE; @@ -46,56 +71,52 @@ int main(int argc, char *argv[]) EigenSparseMatrix sm1(rows,cols), sm2(rows,cols); - int n = rows*cols*density; - std::cout << "n = " << n << "\n"; + nentries = rows*cols*density; + std::cout << "n = " << nentries << "\n"; int dummy; BenchTimer t; t.reset(); t.start(); - for (int k=0; k(0,rows-1) + ei_random(0,cols-1); t.stop(); - double rtime = t.value(); + rtime = t.value(); std::cout << "rtime = " << rtime << " (" << dummy << ")\n\n"; const int Bits = 6; for (;;) { - { - RandomSetter set1(sm1); - t.reset(); t.start(); - for (int k=0; k(0,rows-1),ei_random(0,cols-1)) += 1; - t.stop(); - std::cout << "std::map => \t" << t.value()-rtime - << " nnz=" << set1.nonZeros() << "\n";getchar(); - } - { - RandomSetter set1(sm1); - t.reset(); t.start(); - for (int k=0; k(0,rows-1),ei_random(0,cols-1)) += 1; - t.stop(); - std::cout << "gnu::hash_map => \t" << t.value()-rtime - << " nnz=" << set1.nonZeros() << "\n";getchar(); - } - { - RandomSetter set1(sm1); - t.reset(); t.start(); - for (int k=0; k(0,rows-1),ei_random(0,cols-1)) += 1; - t.stop(); - std::cout << "google::dense => \t" << t.value()-rtime - << " nnz=" << set1.nonZeros() << "\n";getchar(); - } - { - RandomSetter set1(sm1); - t.reset(); t.start(); - for (int k=0; k(0,rows-1),ei_random(0,cols-1)) += 1; - t.stop(); - std::cout << "google::sparse => \t" << t.value()-rtime - << " nnz=" << set1.nonZeros() << "\n";getchar(); - } + dostuff >("std::map ", sm1); + dostuff >("gnu::hash_map", sm1); + dostuff >("google::dense", sm1); + dostuff >("google::sparse", sm1); + +// { +// RandomSetter set1(sm1); +// t.reset(); t.start(); +// for (int k=0; k(0,rows-1),ei_random(0,cols-1)) += 1; +// t.stop(); +// std::cout << "gnu::hash_map => \t" << t.value()-rtime +// << " nnz=" << set1.nonZeros() << "\n";getchar(); +// } +// { +// RandomSetter set1(sm1); +// t.reset(); t.start(); +// for (int k=0; k(0,rows-1),ei_random(0,cols-1)) += 1; +// t.stop(); +// std::cout << "google::dense => \t" << t.value()-rtime +// << " nnz=" << set1.nonZeros() << "\n";getchar(); +// } +// { +// RandomSetter set1(sm1); +// t.reset(); t.start(); +// for (int k=0; k(0,rows-1),ei_random(0,cols-1)) += 1; +// t.stop(); +// std::cout << "google::sparse => \t" << t.value()-rtime +// << " nnz=" << set1.nonZeros() << "\n";getchar(); +// } std::cout << "\n\n"; } diff --git a/cmake/FindCholmod.cmake b/cmake/FindCholmod.cmake index 794c14fb5..444279ab9 100644 --- a/cmake/FindCholmod.cmake +++ b/cmake/FindCholmod.cmake @@ -22,9 +22,7 @@ if(CHOLMOD_LIBRARIES) find_library(AMD_LIBRARY amd PATHS ${CHOLMOD_LIBDIR} $ENV{CHOLMODDIR} ${LIB_INSTALL_DIR}) if (AMD_LIBRARY) set(CHOLMOD_LIBRARIES ${CHOLMOD_LIBRARIES} ${AMD_LIBRARY}) - message ("AMD_LIBRARY found") else (AMD_LIBRARY) - message ("AMD_LIBRARY not found") set(CHOLMOD_LIBRARIES FALSE) endif (AMD_LIBRARY) diff --git a/cmake/FindGoogleHash.cmake b/cmake/FindGoogleHash.cmake new file mode 100644 index 000000000..767996b41 --- /dev/null +++ b/cmake/FindGoogleHash.cmake @@ -0,0 +1,16 @@ + +if (GOOGLEHASH_INCLUDES AND GOOGLEHASH_LIBRARIES) + set(GOOGLEHASH_FIND_QUIETLY TRUE) +endif (GOOGLEHASH_INCLUDES AND GOOGLEHASH_LIBRARIES) + +find_path(GOOGLEHASH_INCLUDES + NAMES + google/dense_hash_map + PATHS + ${INCLUDE_INSTALL_DIR} +) + +include(FindPackageHandleStandardArgs) +find_package_handle_standard_args(GOOGLEHASH DEFAULT_MSG GOOGLEHASH_INCLUDES) + +mark_as_advanced(GOOGLEHASH_INCLUDES) diff --git a/Eigen/src/Sparse/HashMatrix.h b/disabled/HashMatrix.h similarity index 100% rename from Eigen/src/Sparse/HashMatrix.h rename to disabled/HashMatrix.h diff --git a/Eigen/src/Sparse/LinkedVectorMatrix.h b/disabled/LinkedVectorMatrix.h similarity index 100% rename from Eigen/src/Sparse/LinkedVectorMatrix.h rename to disabled/LinkedVectorMatrix.h diff --git a/Eigen/src/Sparse/SparseSetter.h b/disabled/SparseSetter.h similarity index 100% rename from Eigen/src/Sparse/SparseSetter.h rename to disabled/SparseSetter.h diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 1320e9b53..39729fe7e 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -13,14 +13,17 @@ if(TAUCS_FOUND) add_definitions("-DEIGEN_TAUCS_SUPPORT") include_directories(${TAUCS_INCLUDES}) set(EXTERNAL_LIBS ${EXTERNAL_LIBS} ${TAUCS_LIBRARIES}) +else(TAUCS_FOUND) + message("TAUCS not found, this optional backend won't be tested") endif(TAUCS_FOUND) find_package(Cholmod) if(CHOLMOD_FOUND) - message("add EIGEN_CHOLMOD_SUPPORT " ${CHOLMOD_LIBRARIES}) add_definitions("-DEIGEN_CHOLMOD_SUPPORT") include_directories(${CHOLMOD_INCLUDES}) set(EXTERNAL_LIBS ${EXTERNAL_LIBS} ${CHOLMOD_LIBRARIES}) +else(CHOLMOD_FOUND) + message("CHOLMOD not found, this optional backend won't be tested") endif(CHOLMOD_FOUND) find_package(Umfpack) @@ -28,6 +31,8 @@ if(UMFPACK_FOUND) add_definitions("-DEIGEN_UMFPACK_SUPPORT") include_directories(${UMFPACK_INCLUDES}) set(EXTERNAL_LIBS ${EXTERNAL_LIBS} ${UMFPACK_LIBRARIES}) +else(UMFPACK_FOUND) + message("UMFPACK not found, this optional backend won't be tested") endif(UMFPACK_FOUND) find_package(SuperLU) @@ -35,8 +40,18 @@ if(SUPERLU_FOUND) add_definitions("-DEIGEN_SUPERLU_SUPPORT") include_directories(${SUPERLU_INCLUDES}) set(EXTERNAL_LIBS ${EXTERNAL_LIBS} ${SUPERLU_LIBRARIES}) +else(SUPERLU_FOUND) + message("SUPERLU not found, this optional backend won't be tested") endif(SUPERLU_FOUND) +find_package(GoogleHash) +if(GOOGLEHASH_FOUND) + add_definitions("-DEIGEN_GOOGLEHASH_SUPPORT") + include_directories(${GOOGLEHASH_INCLUDES}) +else(GOOGLEHASH_FOUND) + message("Google's hash library not found, those map implementations won't be tested") +endif(GOOGLEHASH_FOUND) + if(CMAKE_COMPILER_IS_GNUCXX) if(CMAKE_SYSTEM_NAME MATCHES Linux) set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g2") @@ -152,6 +167,6 @@ ei_add_test(geometry) ei_add_test(hyperplane) ei_add_test(parametrizedline) ei_add_test(regression) -ei_add_test(sparse ${EI_OFLAG}) +ei_add_test(sparse ) endif(BUILD_TESTS) diff --git a/test/sparse.cpp b/test/sparse.cpp index 6470b1bbf..949aed940 100644 --- a/test/sparse.cpp +++ b/test/sparse.cpp @@ -22,6 +22,14 @@ // License and a copy of the GNU General Public License along with // Eigen. If not, see . +#ifdef __GNUC__ +#include +#endif + +#ifdef EIGEN_GOOGLEHASH_SUPPORT + #include +#endif + #include "main.h" #include #include @@ -72,6 +80,24 @@ initSparse(double density, sparseMat.endFill(); } +template +bool test_random_setter(SparseType& sm, const DenseType& ref, const std::vector& nonzeroCoords) +{ + { + sm.setZero(); + SetterType w(sm); + std::vector remaining = nonzeroCoords; + while(!remaining.empty()) + { + int i = ei_random(0,remaining.size()-1); + w(remaining[i].x(),remaining[i].y()) = ref.coeff(remaining[i].x(),remaining[i].y()); + remaining[i] = remaining.back(); + remaining.pop_back(); + } + } + return sm.isApprox(ref); +} + template void sparse(int rows, int cols) { double density = std::max(8./(rows*cols), 0.01); @@ -157,20 +183,48 @@ template void sparse(int rows, int cols) // VERIFY_IS_APPROX(m, refMat); // random setter - { - m.setZero(); - VERIFY_IS_NOT_APPROX(m, refMat); - SparseSetter, RandomAccessPattern> w(m); - std::vector remaining = nonzeroCoords; - while(!remaining.empty()) - { - int i = ei_random(0,remaining.size()-1); - w->coeffRef(remaining[i].x(),remaining[i].y()) = refMat.coeff(remaining[i].x(),remaining[i].y()); - remaining[i] = remaining.back(); - remaining.pop_back(); - } - } - VERIFY_IS_APPROX(m, refMat); +// { +// m.setZero(); +// VERIFY_IS_NOT_APPROX(m, refMat); +// SparseSetter, RandomAccessPattern> w(m); +// std::vector remaining = nonzeroCoords; +// while(!remaining.empty()) +// { +// int i = ei_random(0,remaining.size()-1); +// w->coeffRef(remaining[i].x(),remaining[i].y()) = refMat.coeff(remaining[i].x(),remaining[i].y()); +// remaining[i] = remaining.back(); +// remaining.pop_back(); +// } +// } +// VERIFY_IS_APPROX(m, refMat); + + VERIFY(( test_random_setter, StdMapTraits> >(m,refMat,nonzeroCoords) )); + #ifdef _HASH_MAP + VERIFY(( test_random_setter, GnuHashMapTraits> >(m,refMat,nonzeroCoords) )); + #endif + #ifdef _DENSE_HASH_MAP_H_ + VERIFY(( test_random_setter, GoogleDenseHashMapTraits> >(m,refMat,nonzeroCoords) )); + #endif + #ifdef _SPARSE_HASH_MAP_H_ + VERIFY(( test_random_setter, GoogleSparseHashMapTraits> >(m,refMat,nonzeroCoords) )); + #endif +// { +// m.setZero(); +// VERIFY_IS_NOT_APPROX(m, refMat); +// // RandomSetter > w(m); +// RandomSetter, GoogleDenseHashMapTraits > w(m); +// // RandomSetter, GnuHashMapTraits > w(m); +// std::vector remaining = nonzeroCoords; +// while(!remaining.empty()) +// { +// int i = ei_random(0,remaining.size()-1); +// w(remaining[i].x(),remaining[i].y()) = refMat.coeff(remaining[i].x(),remaining[i].y()); +// remaining[i] = remaining.back(); +// remaining.pop_back(); +// } +// } +// std::cerr << m.transpose() << "\n\n" << refMat.transpose() << "\n\n"; +// VERIFY_IS_APPROX(m, refMat); // test transpose { @@ -180,7 +234,7 @@ template void sparse(int rows, int cols) VERIFY_IS_APPROX(m2.transpose().eval(), refMat2.transpose().eval()); VERIFY_IS_APPROX(m2.transpose(), refMat2.transpose()); } - +#if 0 // test matrix product { DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows); @@ -306,7 +360,7 @@ template void sparse(int rows, int cols) #endif count++; } - +#endif } void test_sparse()