diff options
22 files changed, 350 insertions, 489 deletions
diff --git a/Eigen/Sparse b/Eigen/Sparse index bca1c4ceb..864f194f4 100644 --- a/Eigen/Sparse +++ b/Eigen/Sparse @@ -11,64 +11,6 @@ #include <cstring> #include <algorithm> -#ifdef EIGEN_GOOGLEHASH_SUPPORT - #include <google/dense_hash_map> -#endif - -#ifdef EIGEN_CHOLMOD_SUPPORT - extern "C" { - #include <cholmod.h> - } -#endif - -#ifdef EIGEN_TAUCS_SUPPORT - // taucs.h declares a lot of mess - #define isnan - #define finite - #define isinf - extern "C" { - #include <taucs.h> - } - #undef isnan - #undef finite - #undef isinf - - #ifdef min - #undef min - #endif - #ifdef max - #undef max - #endif - #ifdef complex - #undef complex - #endif -#endif - -#ifdef EIGEN_SUPERLU_SUPPORT - typedef int int_t; - #include <slu_Cnames.h> - #include <supermatrix.h> - #include <slu_util.h> - - namespace SuperLU_S { - #include <slu_sdefs.h> - } - namespace SuperLU_D { - #include <slu_ddefs.h> - } - namespace SuperLU_C { - #include <slu_cdefs.h> - } - namespace SuperLU_Z { - #include <slu_zdefs.h> - } - namespace Eigen { struct SluMatrix; } -#endif - -#ifdef EIGEN_UMFPACK_SUPPORT - #include <umfpack.h> -#endif - namespace Eigen { /** \defgroup Sparse_Module Sparse module @@ -78,7 +20,7 @@ namespace Eigen { * See the \ref TutorialSparse "Sparse tutorial" * * \code - * #include <Eigen/QR> + * #include <Eigen/Sparse> * \endcode */ @@ -89,13 +31,12 @@ struct Sparse {}; #include "src/Sparse/SparseMatrixBase.h" #include "src/Sparse/CompressedStorage.h" #include "src/Sparse/AmbiVector.h" -#include "src/Sparse/RandomSetter.h" -#include "src/Sparse/SparseBlock.h" #include "src/Sparse/SparseMatrix.h" #include "src/Sparse/DynamicSparseMatrix.h" #include "src/Sparse/MappedSparseMatrix.h" #include "src/Sparse/SparseVector.h" #include "src/Sparse/CoreIterators.h" +#include "src/Sparse/SparseBlock.h" #include "src/Sparse/SparseTranspose.h" #include "src/Sparse/SparseCwiseUnaryOp.h" #include "src/Sparse/SparseCwiseBinaryOp.h" @@ -108,31 +49,11 @@ struct Sparse {}; #include "src/Sparse/SparseTriangularView.h" #include "src/Sparse/SparseSelfAdjointView.h" #include "src/Sparse/TriangularSolver.h" -#include "src/Sparse/SparseLLT.h" -#include "src/Sparse/SparseLDLT.h" -#include "src/Sparse/SparseLU.h" #include "src/Sparse/SparseView.h" -#ifdef EIGEN_CHOLMOD_SUPPORT -# include "src/Sparse/CholmodSupport.h" -#endif - -#ifdef EIGEN_TAUCS_SUPPORT -# include "src/Sparse/TaucsSupport.h" -#endif - -#ifdef EIGEN_SUPERLU_SUPPORT -# include "src/Sparse/SuperLUSupport.h" -#endif - -#ifdef EIGEN_UMFPACK_SUPPORT -# include "src/Sparse/UmfPackSupport.h" -#endif - } // namespace Eigen #include "src/Core/util/EnableMSVCWarnings.h" #endif // EIGEN_SPARSE_MODULE_H -/* vim: set filetype=cpp et sw=2 ts=2 ai: */ diff --git a/Eigen/src/Sparse/MappedSparseMatrix.h b/Eigen/src/Sparse/MappedSparseMatrix.h index 99aeeb106..6fc8adf32 100644 --- a/Eigen/src/Sparse/MappedSparseMatrix.h +++ b/Eigen/src/Sparse/MappedSparseMatrix.h @@ -119,18 +119,6 @@ class MappedSparseMatrix m_innerIndices(innerIndexPtr), m_values(valuePtr) {} - #ifdef EIGEN_TAUCS_SUPPORT - explicit MappedSparseMatrix(taucs_ccs_matrix& taucsMatrix); - #endif - - #ifdef EIGEN_CHOLMOD_SUPPORT - explicit MappedSparseMatrix(cholmod_sparse& cholmodMatrix); - #endif - - #ifdef EIGEN_SUPERLU_SUPPORT - explicit MappedSparseMatrix(SluMatrix& sluMatrix); - #endif - /** Empty destructor */ inline ~MappedSparseMatrix() {} }; diff --git a/Eigen/src/Sparse/SparseMatrixBase.h b/Eigen/src/Sparse/SparseMatrixBase.h index bde8868d5..3ec893119 100644 --- a/Eigen/src/Sparse/SparseMatrixBase.h +++ b/Eigen/src/Sparse/SparseMatrixBase.h @@ -676,18 +676,6 @@ template<typename Derived> class SparseMatrixBase : public EigenBase<Derived> // return res; // } - #ifdef EIGEN_TAUCS_SUPPORT - taucs_ccs_matrix asTaucsMatrix(); - #endif - - #ifdef EIGEN_CHOLMOD_SUPPORT - cholmod_sparse asCholmodMatrix(); - #endif - - #ifdef EIGEN_SUPERLU_SUPPORT - SluMatrix asSluMatrix(); - #endif - protected: bool m_isRValue; diff --git a/Eigen/src/Sparse/SparseUtil.h b/Eigen/src/Sparse/SparseUtil.h index deaf70bc8..423a5ff40 100644 --- a/Eigen/src/Sparse/SparseUtil.h +++ b/Eigen/src/Sparse/SparseUtil.h @@ -75,34 +75,10 @@ EIGEN_SPARSE_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(Derived, /=) #define EIGEN_SPARSE_PUBLIC_INTERFACE(Derived) \ _EIGEN_SPARSE_PUBLIC_INTERFACE(Derived, Eigen::SparseMatrixBase<Derived>) -enum SparseBackend { - DefaultBackend, - Taucs, - Cholmod, - SuperLU, - UmfPack -}; - -// solver flags -enum { - CompleteFactorization = 0x0000, // the default - IncompleteFactorization = 0x0001, - MemoryEfficient = 0x0002, - - // For LLT Cholesky: - SupernodalMultifrontal = 0x0010, - SupernodalLeftLooking = 0x0020, - - // Ordering methods: - NaturalOrdering = 0x0100, // the default - MinimumDegree_AT_PLUS_A = 0x0200, - MinimumDegree_ATA = 0x0300, - ColApproxMinimumDegree = 0x0400, - Metis = 0x0500, - Scotch = 0x0600, - Chaco = 0x0700, - OrderingMask = 0x0f00 -}; +const int CoherentAccessPattern = 0x1; +const int InnerRandomAccessPattern = 0x2 | CoherentAccessPattern; +const int OuterRandomAccessPattern = 0x4 | CoherentAccessPattern; +const int RandomAccessPattern = 0x8 | OuterRandomAccessPattern | InnerRandomAccessPattern; template<typename Derived> class SparseMatrixBase; template<typename _Scalar, int _Flags = 0, typename _Index = int> class SparseMatrix; @@ -126,11 +102,6 @@ template<typename Lhs, typename Rhs, template<typename Lhs, typename Rhs> struct SparseProductReturnType; -const int CoherentAccessPattern = 0x1; -const int InnerRandomAccessPattern = 0x2 | CoherentAccessPattern; -const int OuterRandomAccessPattern = 0x4 | CoherentAccessPattern; -const int RandomAccessPattern = 0x8 | OuterRandomAccessPattern | InnerRandomAccessPattern; - template<typename T> struct ei_eval<T,Sparse> { typedef typename ei_traits<T>::Scalar _Scalar; diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 4b5ee0d92..75e4e116d 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -16,56 +16,6 @@ else(GSL_FOUND) set(GSL_LIBRARIES " ") endif(GSL_FOUND) -set(SPARSE_LIBS "") - -find_package(Taucs) -if(TAUCS_FOUND) - add_definitions("-DEIGEN_TAUCS_SUPPORT") - include_directories(${TAUCS_INCLUDES}) - set(SPARSE_LIBS ${SPARSE_LIBS} ${TAUCS_LIBRARIES}) - ei_add_property(EIGEN_TESTED_BACKENDS "Taucs, ") -else(TAUCS_FOUND) - ei_add_property(EIGEN_MISSING_BACKENDS "Taucs, ") -endif(TAUCS_FOUND) - -find_package(Cholmod) -if(CHOLMOD_FOUND) - add_definitions("-DEIGEN_CHOLMOD_SUPPORT") - include_directories(${CHOLMOD_INCLUDES}) - set(SPARSE_LIBS ${SPARSE_LIBS} ${CHOLMOD_LIBRARIES}) - ei_add_property(EIGEN_TESTED_BACKENDS "Cholmod, ") -else(CHOLMOD_FOUND) - ei_add_property(EIGEN_MISSING_BACKENDS "Cholmod, ") -endif(CHOLMOD_FOUND) - -find_package(Umfpack) -if(UMFPACK_FOUND) - add_definitions("-DEIGEN_UMFPACK_SUPPORT") - include_directories(${UMFPACK_INCLUDES}) - set(SPARSE_LIBS ${SPARSE_LIBS} ${UMFPACK_LIBRARIES}) - ei_add_property(EIGEN_TESTED_BACKENDS "UmfPack, ") -else(UMFPACK_FOUND) - ei_add_property(EIGEN_MISSING_BACKENDS "UmfPack, ") -endif(UMFPACK_FOUND) - -find_package(SuperLU) -if(SUPERLU_FOUND) - add_definitions("-DEIGEN_SUPERLU_SUPPORT") - include_directories(${SUPERLU_INCLUDES}) - set(SPARSE_LIBS ${SPARSE_LIBS} ${SUPERLU_LIBRARIES}) - ei_add_property(EIGEN_TESTED_BACKENDS "SuperLU, ") -else(SUPERLU_FOUND) - ei_add_property(EIGEN_MISSING_BACKENDS "SuperLU, ") -endif(SUPERLU_FOUND) - -find_package(GoogleHash) -if(GOOGLEHASH_FOUND) - add_definitions("-DEIGEN_GOOGLEHASH_SUPPORT") - include_directories(${GOOGLEHASH_INCLUDES}) - ei_add_property(EIGEN_TESTED_BACKENDS "GoogleHash, ") -else(GOOGLEHASH_FOUND) - ei_add_property(EIGEN_MISSING_BACKENDS "GoogleHash, ") -endif(GOOGLEHASH_FOUND) option(EIGEN_TEST_NOQT "Disable Qt support in unit tests" OFF) if(NOT EIGEN_TEST_NOQT) diff --git a/test/sparse_basic.cpp b/test/sparse_basic.cpp index 43d9c6254..3d2210930 100644 --- a/test/sparse_basic.cpp +++ b/test/sparse_basic.cpp @@ -24,40 +24,6 @@ #include "sparse.h" -template<typename SetterType,typename DenseType, typename Scalar, int Options> -bool test_random_setter(SparseMatrix<Scalar,Options>& sm, const DenseType& ref, const std::vector<Vector2i>& nonzeroCoords) -{ - typedef SparseMatrix<Scalar,Options> SparseType; - { - sm.setZero(); - SetterType w(sm); - std::vector<Vector2i> remaining = nonzeroCoords; - while(!remaining.empty()) - { - int i = ei_random<int>(0,static_cast<int>(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<typename SetterType,typename DenseType, typename T> -bool test_random_setter(DynamicSparseMatrix<T>& sm, const DenseType& ref, const std::vector<Vector2i>& nonzeroCoords) -{ - sm.setZero(); - std::vector<Vector2i> remaining = nonzeroCoords; - while(!remaining.empty()) - { - int i = ei_random<int>(0,static_cast<int>(remaining.size())-1); - sm.coeffRef(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<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& ref) { const int rows = ref.rows(); @@ -136,47 +102,6 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re } */ - // test SparseSetters - // coherent setter - // TODO extend the MatrixSetter -// { -// m.setZero(); -// VERIFY_IS_NOT_APPROX(m, refMat); -// SparseSetter<SparseMatrixType, FullyCoherentAccessPattern> w(m); -// for (int i=0; i<nonzeroCoords.size(); ++i) -// { -// w->coeffRef(nonzeroCoords[i].x(),nonzeroCoords[i].y()) = refMat.coeff(nonzeroCoords[i].x(),nonzeroCoords[i].y()); -// } -// } -// VERIFY_IS_APPROX(m, refMat); - - // random setter -// { -// m.setZero(); -// VERIFY_IS_NOT_APPROX(m, refMat); -// SparseSetter<SparseMatrixType, RandomAccessPattern> w(m); -// std::vector<Vector2i> remaining = nonzeroCoords; -// while(!remaining.empty()) -// { -// int i = ei_random<int>(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<RandomSetter<SparseMatrixType, StdMapTraits> >(m,refMat,nonzeroCoords) )); - #ifdef EIGEN_UNORDERED_MAP_SUPPORT - VERIFY(( test_random_setter<RandomSetter<SparseMatrixType, StdUnorderedMapTraits> >(m,refMat,nonzeroCoords) )); - #endif - #ifdef _DENSE_HASH_MAP_H_ - VERIFY(( test_random_setter<RandomSetter<SparseMatrixType, GoogleDenseHashMapTraits> >(m,refMat,nonzeroCoords) )); - #endif - #ifdef _SPARSE_HASH_MAP_H_ - VERIFY(( test_random_setter<RandomSetter<SparseMatrixType, GoogleSparseHashMapTraits> >(m,refMat,nonzeroCoords) )); - #endif - // test insert (inner random) { DenseMatrix m1(rows,cols); @@ -213,22 +138,6 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re VERIFY_IS_APPROX(m2,m1); } - // test RandomSetter - /*{ - SparseMatrixType m1(rows,cols), m2(rows,cols); - DenseMatrix refM1 = DenseMatrix::Zero(rows, rows); - initSparse<Scalar>(density, refM1, m1); - { - Eigen::RandomSetter<SparseMatrixType > setter(m2); - for (int j=0; j<m1.outerSize(); ++j) - for (typename SparseMatrixType::InnerIterator i(m1,j); i; ++i) - setter(i.index(), j) = i.value(); - } - VERIFY_IS_APPROX(m1, m2); - }*/ -// std::cerr << m.transpose() << "\n\n" << refMat.transpose() << "\n\n"; -// VERIFY_IS_APPROX(m, refMat); - // test basic computations { DenseMatrix refM1 = DenseMatrix::Zero(rows, rows); @@ -263,6 +172,17 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re // VERIFY_IS_APPROX(m3.cwise()/refM4, refM3.cwise()/refM4); } + // test transpose + { + DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows); + SparseMatrixType m2(rows, rows); + initSparse<Scalar>(density, refMat2, m2); + VERIFY_IS_APPROX(m2.transpose().eval(), refMat2.transpose().eval()); + VERIFY_IS_APPROX(m2.transpose(), refMat2.transpose()); + + VERIFY_IS_APPROX(SparseMatrixType(m2.adjoint()), refMat2.adjoint()); + } + // test innerVector() { DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows); @@ -292,17 +212,6 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re //refMat2.block(0,j0,rows,n0) = refMat2.block(0,j0,rows,n0) + refMat2.block(0,j1,rows,n0); } - // test transpose - { - DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows); - SparseMatrixType m2(rows, rows); - initSparse<Scalar>(density, refMat2, m2); - VERIFY_IS_APPROX(m2.transpose().eval(), refMat2.transpose().eval()); - VERIFY_IS_APPROX(m2.transpose(), refMat2.transpose()); - - VERIFY_IS_APPROX(SparseMatrixType(m2.adjoint()), refMat2.adjoint()); - } - // test prune { SparseMatrixType m2(rows, rows); diff --git a/test/sparse_solvers.cpp b/test/sparse_solvers.cpp index 9d30af146..30ee72af0 100644 --- a/test/sparse_solvers.cpp +++ b/test/sparse_solvers.cpp @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2008 Daniel Gomez Ferro <dgomezferro@gmail.com> +// Copyright (C) 2008-2010 Gael Guennebaud <g.gael@free.fr> // // Eigen is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public @@ -105,139 +105,6 @@ template<typename Scalar> void sparse_solvers(int rows, int cols) VERIFY_IS_APPROX(refMat2.template triangularView<Lower>().solve(vec2), m2.template triangularView<Lower>().solve(vec3)); } - - // test LLT - { - // TODO fix the issue with complex (see SparseLLT::solveInPlace) - SparseMatrix<Scalar> m2(rows, cols); - DenseMatrix refMat2(rows, cols); - - DenseVector b = DenseVector::Random(cols); - DenseVector refX(cols), x(cols); - - initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, 0, 0); - for(int i=0; i<rows; ++i) - m2.coeffRef(i,i) = refMat2(i,i) = ei_abs(ei_real(refMat2(i,i))); - - refX = refMat2.template selfadjointView<Lower>().llt().solve(b); - if (!NumTraits<Scalar>::IsComplex) - { - x = b; - SparseLLT<SparseMatrix<Scalar> > (m2).solveInPlace(x); - VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: default"); - } - #ifdef EIGEN_CHOLMOD_SUPPORT - x = b; - SparseLLT<SparseMatrix<Scalar> ,Cholmod>(m2).solveInPlace(x); - VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: cholmod"); - #endif - - #ifdef EIGEN_TAUCS_SUPPORT - // TODO fix TAUCS with complexes - if (!NumTraits<Scalar>::IsComplex) - { - x = b; -// SparseLLT<SparseMatrix<Scalar> ,Taucs>(m2,IncompleteFactorization).solveInPlace(x); -// VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (IncompleteFactorization)"); - - x = b; - SparseLLT<SparseMatrix<Scalar> ,Taucs>(m2,SupernodalMultifrontal).solveInPlace(x); - VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (SupernodalMultifrontal)"); - x = b; - SparseLLT<SparseMatrix<Scalar> ,Taucs>(m2,SupernodalLeftLooking).solveInPlace(x); - VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (SupernodalLeftLooking)"); - } - #endif - } - - // test LDLT - { - SparseMatrix<Scalar> m2(rows, cols); - DenseMatrix refMat2(rows, cols); - - DenseVector b = DenseVector::Random(cols); - DenseVector refX(cols), x(cols); - - initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, 0, 0); - for(int i=0; i<rows; ++i) - m2.coeffRef(i,i) = refMat2(i,i) = ei_abs(ei_real(refMat2(i,i))); - - refX = refMat2.template selfadjointView<Upper>().ldlt().solve(b); - typedef SparseMatrix<Scalar,Upper|SelfAdjoint> SparseSelfAdjointMatrix; - x = b; - SparseLDLT<SparseSelfAdjointMatrix> ldlt(m2); - if (ldlt.succeeded()) - ldlt.solveInPlace(x); - else - std::cerr << "warning LDLT failed\n"; - - VERIFY_IS_APPROX(refMat2.template selfadjointView<Upper>() * x, b); - VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LDLT: default"); - } - - // test LU - { - static int count = 0; - SparseMatrix<Scalar> m2(rows, cols); - DenseMatrix refMat2(rows, cols); - - DenseVector b = DenseVector::Random(cols); - DenseVector refX(cols), x(cols); - - initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag, &zeroCoords, &nonzeroCoords); - - FullPivLU<DenseMatrix> refLu(refMat2); - refX = refLu.solve(b); - #if defined(EIGEN_SUPERLU_SUPPORT) || defined(EIGEN_UMFPACK_SUPPORT) - Scalar refDet = refLu.determinant(); - #endif - x.setZero(); - // // SparseLU<SparseMatrix<Scalar> > (m2).solve(b,&x); - // // VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LU: default"); - #ifdef EIGEN_SUPERLU_SUPPORT - { - x.setZero(); - SparseLU<SparseMatrix<Scalar>,SuperLU> slu(m2); - if (slu.succeeded()) - { - if (slu.solve(b,&x)) { - VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LU: SuperLU"); - } - // std::cerr << refDet << " == " << slu.determinant() << "\n"; - if (slu.solve(b, &x, SvTranspose)) { - VERIFY(b.isApprox(m2.transpose() * x, test_precision<Scalar>())); - } - - if (slu.solve(b, &x, SvAdjoint)) { - VERIFY(b.isApprox(m2.adjoint() * x, test_precision<Scalar>())); - } - - if (count==0) { - VERIFY_IS_APPROX(refDet,slu.determinant()); // FIXME det is not very stable for complex - } - } - } - #endif - #ifdef EIGEN_UMFPACK_SUPPORT - { - // check solve - x.setZero(); - SparseLU<SparseMatrix<Scalar>,UmfPack> slu(m2); - if (slu.succeeded()) { - if (slu.solve(b,&x)) { - if (count==0) { - VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LU: umfpack"); // FIXME solve is not very stable for complex - } - } - VERIFY_IS_APPROX(refDet,slu.determinant()); - // TODO check the extracted data - //std::cerr << slu.matrixL() << "\n"; - } - } - #endif - count++; - } - } void test_sparse_solvers() diff --git a/unsupported/Eigen/CholmodSupport b/unsupported/Eigen/CholmodSupport new file mode 100644 index 000000000..e1c50e536 --- /dev/null +++ b/unsupported/Eigen/CholmodSupport @@ -0,0 +1,33 @@ +#ifndef EIGEN_CHOLMODSUPPORT_MODULE_H +#define EIGEN_CHOLMODSUPPORT_MODULE_H + +#include "SparseExtra" + +#include "src/Core/util/DisableMSVCWarnings.h" + +#ifdef EIGEN_CHOLMOD_SUPPORT + extern "C" { + #include <cholmod.h> + } +#endif + +namespace Eigen { + +/** \defgroup CholmodSupport_Module Cholmod Support module + * + * \nonstableyet + * + * + * \code + * #include <Eigen/CholmodSupport> + * \endcode + */ + +#include "src/SparseExtra/CholmodSupport.h" + +} // namespace Eigen + +#include "src/Core/util/EnableMSVCWarnings.h" + +#endif // EIGEN_CHOLMODSUPPORT_MODULE_H + diff --git a/unsupported/Eigen/SparseExtra b/unsupported/Eigen/SparseExtra new file mode 100644 index 000000000..ba44a9c42 --- /dev/null +++ b/unsupported/Eigen/SparseExtra @@ -0,0 +1,67 @@ +#ifndef EIGEN_SPARSE_EXTRA_MODULE_H +#define EIGEN_SPARSE_EXTRA_MODULE_H + +#include "../Eigen/Sparse" + +#include "src/Core/util/DisableMSVCWarnings.h" + +#include <vector> +#include <map> +#include <cstdlib> +#include <cstring> +#include <algorithm> + +#ifdef EIGEN_GOOGLEHASH_SUPPORT + #include <google/dense_hash_map> +#endif + +namespace Eigen { + +/** \defgroup SparseExtra_Module SparseExtra module + * + * This module contains some experimental features extending the sparse module. + * + * \code + * #include <Eigen/SparseExtra> + * \endcode + */ + +enum SparseBackend { + DefaultBackend, + Taucs, + Cholmod, + SuperLU, + UmfPack +}; + +// solver flags +enum { + CompleteFactorization = 0x0000, // the default + IncompleteFactorization = 0x0001, + MemoryEfficient = 0x0002, + + // For LLT Cholesky: + SupernodalMultifrontal = 0x0010, + SupernodalLeftLooking = 0x0020, + + // Ordering methods: + NaturalOrdering = 0x0100, // the default + MinimumDegree_AT_PLUS_A = 0x0200, + MinimumDegree_ATA = 0x0300, + ColApproxMinimumDegree = 0x0400, + Metis = 0x0500, + Scotch = 0x0600, + Chaco = 0x0700, + OrderingMask = 0x0f00 +}; + +#include "src/SparseExtra/RandomSetter.h" +#include "src/SparseExtra/SparseLLT.h" +#include "src/SparseExtra/SparseLDLT.h" +#include "src/SparseExtra/SparseLU.h" + +} // namespace Eigen + +#include "src/Core/util/EnableMSVCWarnings.h" + +#endif // EIGEN_SPARSE_EXTRA_MODULE_H diff --git a/unsupported/Eigen/SuperLUSupport b/unsupported/Eigen/SuperLUSupport new file mode 100644 index 000000000..6a2a3b15a --- /dev/null +++ b/unsupported/Eigen/SuperLUSupport @@ -0,0 +1,44 @@ +#ifndef EIGEN_SUPERLUSUPPORT_MODULE_H +#define EIGEN_SUPERLUSUPPORT_MODULE_H + +#include "SparseExtra" + +#include "src/Core/util/DisableMSVCWarnings.h" + +typedef int int_t; +#include <slu_Cnames.h> +#include <supermatrix.h> +#include <slu_util.h> + +namespace SuperLU_S { +#include <slu_sdefs.h> +} +namespace SuperLU_D { +#include <slu_ddefs.h> +} +namespace SuperLU_C { +#include <slu_cdefs.h> +} +namespace SuperLU_Z { +#include <slu_zdefs.h> +} +namespace Eigen { struct SluMatrix; } + +namespace Eigen { + +/** \defgroup SuperLUSupport_Module Super LU support + * + * \nonstableyet + * + * \code + * #include <Eigen/SuperLUSupport> + * \endcode + */ + +#include "src/SparseExtra/SuperLUSupport.h" + +} // namespace Eigen + +#include "src/Core/util/EnableMSVCWarnings.h" + +#endif // EIGEN_SUPERLUSUPPORT_MODULE_H diff --git a/unsupported/Eigen/TaucsSupport b/unsupported/Eigen/TaucsSupport new file mode 100644 index 000000000..93e4bd7aa --- /dev/null +++ b/unsupported/Eigen/TaucsSupport @@ -0,0 +1,49 @@ +#ifndef EIGEN_TAUCSSUPPORT_MODULE_H +#define EIGEN_TAUCSSUPPORT_MODULE_H + +#include "SparseExtra" + +#include "src/Core/util/DisableMSVCWarnings.h" + +// taucs.h declares a lot of mess +#define isnan +#define finite +#define isinf +extern "C" { + #include <taucs.h> +} +#undef isnan +#undef finite +#undef isinf + +#ifdef min + #undef min +#endif +#ifdef max + #undef max +#endif +#ifdef complex + #undef complex +#endif + +namespace Eigen { + +/** \defgroup TaucsSupport_Module Taucs support module + * + * \nonstableyet + * + * + * \code + * #include <Eigen/TaucsSupport> + * \endcode + */ + +#ifdef EIGEN_TAUCS_SUPPORT +# include "src/SparseExtra/TaucsSupport.h" +#endif + +} // namespace Eigen + +#include "src/Core/util/EnableMSVCWarnings.h" + +#endif // EIGEN_TAUCSSUPPORT_MODULE_H diff --git a/unsupported/Eigen/UmfPackSupport b/unsupported/Eigen/UmfPackSupport new file mode 100644 index 000000000..6fcf40ec0 --- /dev/null +++ b/unsupported/Eigen/UmfPackSupport @@ -0,0 +1,28 @@ +#ifndef EIGEN_UMFPACKSUPPORT_MODULE_H +#define EIGEN_UMFPACKSUPPORT_MODULE_H + +#include "SparseExtra" + +#include "src/Core/util/DisableMSVCWarnings.h" + +#include <umfpack.h> + +namespace Eigen { + +/** \defgroup Sparse_Module Sparse module + * + * \nonstableyet + * + * + * \code + * #include <Eigen/UmfPackSupport> + * \endcode + */ + +#include "src/SparseExtra/UmfPackSupport.h" + +} // namespace Eigen + +#include "src/Core/util/EnableMSVCWarnings.h" + +#endif // EIGEN_UMFPACKSUPPORT_MODULE_H diff --git a/unsupported/Eigen/src/CMakeLists.txt b/unsupported/Eigen/src/CMakeLists.txt index 035a84ca8..4578c2e85 100644 --- a/unsupported/Eigen/src/CMakeLists.txt +++ b/unsupported/Eigen/src/CMakeLists.txt @@ -6,3 +6,4 @@ ADD_SUBDIRECTORY(MoreVectorization) # ADD_SUBDIRECTORY(Skyline) ADD_SUBDIRECTORY(MatrixFunctions) ADD_SUBDIRECTORY(Polynomials) +ADD_SUBDIRECTORY(SparseExtra) diff --git a/Eigen/src/Sparse/CholmodSupport.h b/unsupported/Eigen/src/SparseExtra/CholmodSupport.h index a8d7a8fec..92e79f2a6 100644 --- a/Eigen/src/Sparse/CholmodSupport.h +++ b/unsupported/Eigen/src/SparseExtra/CholmodSupport.h @@ -54,17 +54,17 @@ void ei_cholmod_configure_matrix(CholmodType& mat) } } -template<typename Derived> -cholmod_sparse SparseMatrixBase<Derived>::asCholmodMatrix() +template<typename MatrixType> +cholmod_sparse ei_asCholmodMatrix(MatrixType& mat) { - typedef typename Derived::Scalar Scalar; + typedef typename MatrixType::Scalar Scalar; cholmod_sparse res; - res.nzmax = nonZeros(); - res.nrow = rows();; - res.ncol = cols(); - res.p = derived()._outerIndexPtr(); - res.i = derived()._innerIndexPtr(); - res.x = derived()._valuePtr(); + res.nzmax = mat.nonZeros(); + res.nrow = mat.rows();; + res.ncol = mat.cols(); + res.p = mat._outerIndexPtr(); + res.i = mat._innerIndexPtr(); + res.x = mat._valuePtr(); res.xtype = CHOLMOD_REAL; res.itype = CHOLMOD_INT; res.sorted = 1; @@ -75,11 +75,11 @@ cholmod_sparse SparseMatrixBase<Derived>::asCholmodMatrix() ei_cholmod_configure_matrix<Scalar>(res); - if (Derived::Flags & SelfAdjoint) + if (MatrixType::Flags & SelfAdjoint) { - if (Derived::Flags & Upper) + if (MatrixType::Flags & Upper) res.stype = 1; - else if (Derived::Flags & Lower) + else if (MatrixType::Flags & Lower) res.stype = -1; else res.stype = 0; @@ -109,15 +109,12 @@ cholmod_dense ei_cholmod_map_eigen_to_dense(MatrixBase<Derived>& mat) return res; } -template<typename Scalar, int Flags, typename _Index> -MappedSparseMatrix<Scalar,Flags,_Index>::MappedSparseMatrix(cholmod_sparse& cm) +template<typename Scalar, int Flags, typename Index> +MappedSparseMatrix<Scalar,Flags,Index> ei_map_cholmod(cholmod_sparse& cm) { - m_innerSize = cm.nrow; - m_outerSize = cm.ncol; - m_outerIndex = reinterpret_cast<Index*>(cm.p); - m_innerIndices = reinterpret_cast<Index*>(cm.i); - m_values = reinterpret_cast<Scalar*>(cm.x); - m_nnz = m_outerIndex[cm.ncol]; + return MappedSparseMatrix<Scalar,Flags,Index> + (cm.nrow, cm.ncol, reinterpret_cast<Index*>(cm.p)[cm.ncol], + reinterpret_cast<Index*>(cm.p), reinterpret_cast<Index*>(cm.i),reinterpret_cast<Scalar*>(cm.x) ); } template<typename MatrixType> @@ -178,7 +175,7 @@ void SparseLLT<MatrixType,Cholmod>::compute(const MatrixType& a) m_cholmodFactor = 0; } - cholmod_sparse A = const_cast<MatrixType&>(a).asCholmodMatrix(); + cholmod_sparse A = ei_asCholmodMatrix(const_cast<MatrixType&>(a)); // m_cholmod.supernodal = CHOLMOD_AUTO; // TODO // if (m_flags&IncompleteFactorization) @@ -209,7 +206,7 @@ SparseLLT<MatrixType,Cholmod>::matrixL() const ei_assert(!(m_status & SupernodalFactorIsDirty)); cholmod_sparse* cmRes = cholmod_factor_to_sparse(m_cholmodFactor, &m_cholmod); - const_cast<typename Base::CholMatrixType&>(m_matrix) = MappedSparseMatrix<Scalar>(*cmRes); + const_cast<typename Base::CholMatrixType&>(m_matrix) = ei_map_cholmod<Scalar,ColMajor,Index>(*cmRes); free(cmRes); m_status = (m_status & ~MatrixLIsDirty); @@ -235,7 +232,7 @@ bool SparseLLT<MatrixType,Cholmod>::solveInPlace(MatrixBase<Derived> &b) const cholmod_dense* x = cholmod_solve(CHOLMOD_A, m_cholmodFactor, &cdb, &m_cholmod); if(!x) { - //std::cerr << "Eigen: cholmod_solve failed\n"; + std::cerr << "Eigen: cholmod_solve failed\n"; return false; } b = Matrix<typename Base::Scalar,Dynamic,1>::Map(reinterpret_cast<typename Base::Scalar*>(x->x),b.rows()); diff --git a/Eigen/src/Sparse/RandomSetter.h b/unsupported/Eigen/src/SparseExtra/RandomSetter.h index 18777e23d..18777e23d 100644 --- a/Eigen/src/Sparse/RandomSetter.h +++ b/unsupported/Eigen/src/SparseExtra/RandomSetter.h diff --git a/Eigen/src/Sparse/SparseLDLT.h b/unsupported/Eigen/src/SparseExtra/SparseLDLT.h index a6785d0af..a6785d0af 100644 --- a/Eigen/src/Sparse/SparseLDLT.h +++ b/unsupported/Eigen/src/SparseExtra/SparseLDLT.h diff --git a/Eigen/src/Sparse/SparseLLT.h b/unsupported/Eigen/src/SparseExtra/SparseLLT.h index 47d58f8e6..47d58f8e6 100644 --- a/Eigen/src/Sparse/SparseLLT.h +++ b/unsupported/Eigen/src/SparseExtra/SparseLLT.h diff --git a/Eigen/src/Sparse/SparseLU.h b/unsupported/Eigen/src/SparseExtra/SparseLU.h index 0211b78e2..0211b78e2 100644 --- a/Eigen/src/Sparse/SparseLU.h +++ b/unsupported/Eigen/src/SparseExtra/SparseLU.h diff --git a/Eigen/src/Sparse/SuperLUSupport.h b/unsupported/Eigen/src/SparseExtra/SuperLUSupport.h index d93f69df8..9766d8fa2 100644 --- a/Eigen/src/Sparse/SuperLUSupport.h +++ b/unsupported/Eigen/src/SparseExtra/SuperLUSupport.h @@ -260,32 +260,24 @@ struct SluMatrixMapHelper<SparseMatrixBase<Derived> > } }; -template<typename Derived> -SluMatrix SparseMatrixBase<Derived>::asSluMatrix() +template<typename MatrixType> +SluMatrix ei_asSluMatrix(MatrixType& mat) { - return SluMatrix::Map(derived()); + return SluMatrix::Map(mat); } /** View a Super LU matrix as an Eigen expression */ -template<typename Scalar, int Flags, typename _Index> -MappedSparseMatrix<Scalar,Flags,_Index>::MappedSparseMatrix(SluMatrix& sluMat) +template<typename Scalar, int Flags, typename Index> +MappedSparseMatrix<Scalar,Flags,Index> ei_map_superlu(SluMatrix& sluMat) { - if ((Flags&RowMajorBit)==RowMajorBit) - { - assert(sluMat.Stype == SLU_NR); - m_innerSize = sluMat.ncol; - m_outerSize = sluMat.nrow; - } - else - { - assert(sluMat.Stype == SLU_NC); - m_innerSize = sluMat.nrow; - m_outerSize = sluMat.ncol; - } - m_outerIndex = sluMat.storage.outerInd; - m_innerIndices = sluMat.storage.innerInd; - m_values = reinterpret_cast<Scalar*>(sluMat.storage.values); - m_nnz = sluMat.storage.outerInd[m_outerSize]; + ei_assert((Flags&RowMajor)==RowMajor && sluMat.Stype == SLU_NR + || (Flags&ColMajor)==ColMajor && sluMat.Stype == SLU_NC); + + Index outerSize = (Flags&RowMajor)==RowMajor ? sluMat.ncol : sluMat.nrow; + + return MappedSparseMatrix<Scalar,Flags,Index>( + sluMat.nrow, sluMat.ncol, sluMat.storage.outerInd[outerSize], + sluMat.storage.outerInd, sluMat.storage.innerInd, reinterpret_cast<Scalar*>(sluMat.storage.values) ); } template<typename MatrixType> @@ -401,7 +393,7 @@ void SparseLU<MatrixType,SuperLU>::compute(const MatrixType& a) m_sluOptions.ColPerm = NATURAL; }; - m_sluA = m_matrix.asSluMatrix(); + m_sluA = ei_asSluMatrix(m_matrix); memset(&m_sluL,0,sizeof m_sluL); memset(&m_sluU,0,sizeof m_sluU); //m_sluEqued = 'B'; diff --git a/Eigen/src/Sparse/TaucsSupport.h b/unsupported/Eigen/src/SparseExtra/TaucsSupport.h index c189e0127..1e7adfef4 100644 --- a/Eigen/src/Sparse/TaucsSupport.h +++ b/unsupported/Eigen/src/SparseExtra/TaucsSupport.h @@ -25,16 +25,18 @@ #ifndef EIGEN_TAUCSSUPPORT_H #define EIGEN_TAUCSSUPPORT_H -template<typename Derived> -taucs_ccs_matrix SparseMatrixBase<Derived>::asTaucsMatrix() +template<typename MatrixType> +taucs_ccs_matrix ei_asTaucsMatrix(MatrixType& mat) { + typedef typename MatrixType::Scalar Scalar; +enum { Flags = MatrixType::Flags }; taucs_ccs_matrix res; - res.n = cols(); - res.m = rows(); + res.n = mat.cols(); + res.m = mat.rows(); res.flags = 0; - res.colptr = derived()._outerIndexPtr(); - res.rowind = derived()._innerIndexPtr(); - res.values.v = derived()._valuePtr(); + res.colptr = mat._outerIndexPtr(); + res.rowind = mat._innerIndexPtr(); + res.values.v = mat._valuePtr(); if (ei_is_same_type<Scalar,int>::ret) res.flags |= TAUCS_INT; else if (ei_is_same_type<Scalar,float>::ret) @@ -63,15 +65,12 @@ taucs_ccs_matrix SparseMatrixBase<Derived>::asTaucsMatrix() return res; } -template<typename Scalar, int Flags, typename _Index> -MappedSparseMatrix<Scalar,Flags,_Index>::MappedSparseMatrix(taucs_ccs_matrix& taucsMat) +template<typename Scalar, int Flags, typename Index> +MappedSparseMatrix<Scalar,Flags,Index> ei_map_taucs(taucs_ccs_matrix& taucsMat) { - m_innerSize = taucsMat.m; - m_outerSize = taucsMat.n; - m_outerIndex = taucsMat.colptr; - m_innerIndices = taucsMat.rowind; - m_values = reinterpret_cast<Scalar*>(taucsMat.values.v); - m_nnz = taucsMat.colptr[taucsMat.n]; + return MappedSparseMatrix<Scalar,Flags,Index> + (taucsMat.m, taucsMat.n, taucsMat.colptr[taucsMat.n], + taucsMat.colptr, taucsMat.rowind, reinterpret_cast<Scalar*>(taucsMat.values.v)); } template<typename MatrixType> @@ -79,6 +78,7 @@ class SparseLLT<MatrixType,Taucs> : public SparseLLT<MatrixType> { protected: typedef SparseLLT<MatrixType> Base; + typedef typename Base::Index Index; typedef typename Base::Scalar Scalar; typedef typename Base::RealScalar RealScalar; typedef typename Base::CholMatrixType CholMatrixType; @@ -128,7 +128,7 @@ void SparseLLT<MatrixType,Taucs>::compute(const MatrixType& a) m_taucsSupernodalFactor = 0; } - taucs_ccs_matrix taucsMatA = const_cast<MatrixType&>(a).asTaucsMatrix(); + taucs_ccs_matrix taucsMatA = ei_asTaucsMatrix(const_cast<MatrixType&>(a)); if (m_flags & IncompleteFactorization) { @@ -140,7 +140,7 @@ void SparseLLT<MatrixType,Taucs>::compute(const MatrixType& a) } // the matrix returned by Taucs is not necessarily sorted, // so let's copy it in two steps - DynamicSparseMatrix<Scalar,RowMajor> tmp = MappedSparseMatrix<Scalar>(*taucsRes); + DynamicSparseMatrix<Scalar,RowMajor> tmp = ei_map_taucs<Scalar,ColMajor,Index>(*taucsRes); m_matrix = tmp; free(taucsRes); m_status = (m_status & ~(CompleteFactorization|MatrixLIsDirty)) @@ -176,7 +176,7 @@ SparseLLT<MatrixType,Taucs>::matrixL() const // the matrix returned by Taucs is not necessarily sorted, // so let's copy it in two steps - DynamicSparseMatrix<Scalar,RowMajor> tmp = MappedSparseMatrix<Scalar>(*taucsL); + DynamicSparseMatrix<Scalar,RowMajor> tmp = ei_map_taucs<Scalar,ColMajor,Index>(*taucsL); const_cast<typename Base::CholMatrixType&>(m_matrix) = tmp; free(taucsL); m_status = (m_status & ~MatrixLIsDirty); @@ -197,7 +197,7 @@ void SparseLLT<MatrixType,Taucs>::solveInPlace(MatrixBase<Derived> &b) const } else if (m_flags & IncompleteFactorization) { - taucs_ccs_matrix taucsLLT = const_cast<typename Base::CholMatrixType&>(m_matrix).asTaucsMatrix(); + taucs_ccs_matrix taucsLLT = ei_asTaucsMatrix(const_cast<typename Base::CholMatrixType&>(m_matrix)); typename ei_plain_matrix_type<Derived>::type x(b.rows()); for (int j=0; j<b.cols(); ++j) { diff --git a/Eigen/src/Sparse/UmfPackSupport.h b/unsupported/Eigen/src/SparseExtra/UmfPackSupport.h index 950624758..950624758 100644 --- a/Eigen/src/Sparse/UmfPackSupport.h +++ b/unsupported/Eigen/src/SparseExtra/UmfPackSupport.h diff --git a/unsupported/test/CMakeLists.txt b/unsupported/test/CMakeLists.txt index f01d99c36..c7cfc8881 100644 --- a/unsupported/test/CMakeLists.txt +++ b/unsupported/test/CMakeLists.txt @@ -1,6 +1,57 @@ find_package(Adolc) -include_directories(../../test) +set(SPARSE_LIBS "") + +find_package(Taucs) +if(TAUCS_FOUND) + add_definitions("-DEIGEN_TAUCS_SUPPORT") + include_directories(${TAUCS_INCLUDES}) + set(SPARSE_LIBS ${SPARSE_LIBS} ${TAUCS_LIBRARIES}) + ei_add_property(EIGEN_TESTED_BACKENDS "Taucs, ") +else(TAUCS_FOUND) + ei_add_property(EIGEN_MISSING_BACKENDS "Taucs, ") +endif(TAUCS_FOUND) + +find_package(Cholmod) +if(CHOLMOD_FOUND) + add_definitions("-DEIGEN_CHOLMOD_SUPPORT") + include_directories(${CHOLMOD_INCLUDES}) + set(SPARSE_LIBS ${SPARSE_LIBS} ${CHOLMOD_LIBRARIES}) + ei_add_property(EIGEN_TESTED_BACKENDS "Cholmod, ") +else(CHOLMOD_FOUND) + ei_add_property(EIGEN_MISSING_BACKENDS "Cholmod, ") +endif(CHOLMOD_FOUND) + +find_package(Umfpack) +if(UMFPACK_FOUND) + add_definitions("-DEIGEN_UMFPACK_SUPPORT") + include_directories(${UMFPACK_INCLUDES}) + set(SPARSE_LIBS ${SPARSE_LIBS} ${UMFPACK_LIBRARIES}) + ei_add_property(EIGEN_TESTED_BACKENDS "UmfPack, ") +else(UMFPACK_FOUND) + ei_add_property(EIGEN_MISSING_BACKENDS "UmfPack, ") +endif(UMFPACK_FOUND) + +find_package(SuperLU) +if(SUPERLU_FOUND) + add_definitions("-DEIGEN_SUPERLU_SUPPORT") + include_directories(${SUPERLU_INCLUDES}) + set(SPARSE_LIBS ${SPARSE_LIBS} ${SUPERLU_LIBRARIES}) + ei_add_property(EIGEN_TESTED_BACKENDS "SuperLU, ") +else(SUPERLU_FOUND) + ei_add_property(EIGEN_MISSING_BACKENDS "SuperLU, ") +endif(SUPERLU_FOUND) + +find_package(GoogleHash) +if(GOOGLEHASH_FOUND) + add_definitions("-DEIGEN_GOOGLEHASH_SUPPORT") + include_directories(${GOOGLEHASH_INCLUDES}) + ei_add_property(EIGEN_TESTED_BACKENDS "GoogleHash, ") +else(GOOGLEHASH_FOUND) + ei_add_property(EIGEN_MISSING_BACKENDS "GoogleHash, ") +endif(GOOGLEHASH_FOUND) + +include_directories(../../test ../../unsupported ../../Eigen) if(ADOLC_FOUND) include_directories(${ADOLC_INCLUDES}) @@ -19,6 +70,11 @@ ei_add_test(matrix_function) ei_add_test(alignedvector3) ei_add_test(FFT) +ei_add_test(sparse_llt " " "${SPARSE_LIBS}") +ei_add_test(sparse_ldlt " " "${SPARSE_LIBS}") +ei_add_test(sparse_lu " " "${SPARSE_LIBS}") +ei_add_test(sparse_extra " " " ") + find_package(FFTW) if(FFTW_FOUND) ei_add_test(FFTW "-DEIGEN_FFTW_DEFAULT " "-lfftw3 -lfftw3f -lfftw3l" ) |