diff --git a/Eigen/Sparse b/Eigen/Sparse index bca1c4ceb..864f194f4 100644 --- a/Eigen/Sparse +++ b/Eigen/Sparse @@ -11,64 +11,6 @@ #include #include -#ifdef EIGEN_GOOGLEHASH_SUPPORT - #include -#endif - -#ifdef EIGEN_CHOLMOD_SUPPORT - extern "C" { - #include - } -#endif - -#ifdef EIGEN_TAUCS_SUPPORT - // taucs.h declares a lot of mess - #define isnan - #define finite - #define isinf - extern "C" { - #include - } - #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 - #include - #include - - namespace SuperLU_S { - #include - } - namespace SuperLU_D { - #include - } - namespace SuperLU_C { - #include - } - namespace SuperLU_Z { - #include - } - namespace Eigen { struct SluMatrix; } -#endif - -#ifdef EIGEN_UMFPACK_SUPPORT - #include -#endif - namespace Eigen { /** \defgroup Sparse_Module Sparse module @@ -78,7 +20,7 @@ namespace Eigen { * See the \ref TutorialSparse "Sparse tutorial" * * \code - * #include + * #include * \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 class SparseMatrixBase : public EigenBase // 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) -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 class SparseMatrixBase; template class SparseMatrix; @@ -126,11 +102,6 @@ template struct SparseProductReturnType; -const int CoherentAccessPattern = 0x1; -const int InnerRandomAccessPattern = 0x2 | CoherentAccessPattern; -const int OuterRandomAccessPattern = 0x4 | CoherentAccessPattern; -const int RandomAccessPattern = 0x8 | OuterRandomAccessPattern | InnerRandomAccessPattern; - template struct ei_eval { typedef typename ei_traits::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 -bool test_random_setter(SparseMatrix& sm, const DenseType& ref, const std::vector& nonzeroCoords) -{ - typedef SparseMatrix SparseType; - { - sm.setZero(); - SetterType w(sm); - std::vector remaining = nonzeroCoords; - while(!remaining.empty()) - { - int i = ei_random(0,static_cast(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 -bool test_random_setter(DynamicSparseMatrix& sm, const DenseType& ref, const std::vector& nonzeroCoords) -{ - sm.setZero(); - std::vector remaining = nonzeroCoords; - while(!remaining.empty()) - { - int i = ei_random(0,static_cast(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 void sparse_basic(const SparseMatrixType& ref) { const int rows = ref.rows(); @@ -136,47 +102,6 @@ template void sparse_basic(const SparseMatrixType& re } */ - // test SparseSetters - // coherent setter - // TODO extend the MatrixSetter -// { -// m.setZero(); -// VERIFY_IS_NOT_APPROX(m, refMat); -// SparseSetter w(m); -// for (int i=0; icoeffRef(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 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 >(m,refMat,nonzeroCoords) )); - #ifdef EIGEN_UNORDERED_MAP_SUPPORT - VERIFY(( test_random_setter >(m,refMat,nonzeroCoords) )); - #endif - #ifdef _DENSE_HASH_MAP_H_ - VERIFY(( test_random_setter >(m,refMat,nonzeroCoords) )); - #endif - #ifdef _SPARSE_HASH_MAP_H_ - VERIFY(( test_random_setter >(m,refMat,nonzeroCoords) )); - #endif - // test insert (inner random) { DenseMatrix m1(rows,cols); @@ -213,22 +138,6 @@ template 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(density, refM1, m1); - { - Eigen::RandomSetter setter(m2); - for (int j=0; j 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(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 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(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 +// Copyright (C) 2008-2010 Gael Guennebaud // // 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 void sparse_solvers(int rows, int cols) VERIFY_IS_APPROX(refMat2.template triangularView().solve(vec2), m2.template triangularView().solve(vec3)); } - - // test LLT - { - // TODO fix the issue with complex (see SparseLLT::solveInPlace) - SparseMatrix m2(rows, cols); - DenseMatrix refMat2(rows, cols); - - DenseVector b = DenseVector::Random(cols); - DenseVector refX(cols), x(cols); - - initSparse(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, 0, 0); - for(int i=0; i().llt().solve(b); - if (!NumTraits::IsComplex) - { - x = b; - SparseLLT > (m2).solveInPlace(x); - VERIFY(refX.isApprox(x,test_precision()) && "LLT: default"); - } - #ifdef EIGEN_CHOLMOD_SUPPORT - x = b; - SparseLLT ,Cholmod>(m2).solveInPlace(x); - VERIFY(refX.isApprox(x,test_precision()) && "LLT: cholmod"); - #endif - - #ifdef EIGEN_TAUCS_SUPPORT - // TODO fix TAUCS with complexes - if (!NumTraits::IsComplex) - { - x = b; -// SparseLLT ,Taucs>(m2,IncompleteFactorization).solveInPlace(x); -// VERIFY(refX.isApprox(x,test_precision()) && "LLT: taucs (IncompleteFactorization)"); - - x = b; - SparseLLT ,Taucs>(m2,SupernodalMultifrontal).solveInPlace(x); - VERIFY(refX.isApprox(x,test_precision()) && "LLT: taucs (SupernodalMultifrontal)"); - x = b; - SparseLLT ,Taucs>(m2,SupernodalLeftLooking).solveInPlace(x); - VERIFY(refX.isApprox(x,test_precision()) && "LLT: taucs (SupernodalLeftLooking)"); - } - #endif - } - - // test LDLT - { - SparseMatrix m2(rows, cols); - DenseMatrix refMat2(rows, cols); - - DenseVector b = DenseVector::Random(cols); - DenseVector refX(cols), x(cols); - - initSparse(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, 0, 0); - for(int i=0; i().ldlt().solve(b); - typedef SparseMatrix SparseSelfAdjointMatrix; - x = b; - SparseLDLT ldlt(m2); - if (ldlt.succeeded()) - ldlt.solveInPlace(x); - else - std::cerr << "warning LDLT failed\n"; - - VERIFY_IS_APPROX(refMat2.template selfadjointView() * x, b); - VERIFY(refX.isApprox(x,test_precision()) && "LDLT: default"); - } - - // test LU - { - static int count = 0; - SparseMatrix m2(rows, cols); - DenseMatrix refMat2(rows, cols); - - DenseVector b = DenseVector::Random(cols); - DenseVector refX(cols), x(cols); - - initSparse(density, refMat2, m2, ForceNonZeroDiag, &zeroCoords, &nonzeroCoords); - - FullPivLU refLu(refMat2); - refX = refLu.solve(b); - #if defined(EIGEN_SUPERLU_SUPPORT) || defined(EIGEN_UMFPACK_SUPPORT) - Scalar refDet = refLu.determinant(); - #endif - x.setZero(); - // // SparseLU > (m2).solve(b,&x); - // // VERIFY(refX.isApprox(x,test_precision()) && "LU: default"); - #ifdef EIGEN_SUPERLU_SUPPORT - { - x.setZero(); - SparseLU,SuperLU> slu(m2); - if (slu.succeeded()) - { - if (slu.solve(b,&x)) { - VERIFY(refX.isApprox(x,test_precision()) && "LU: SuperLU"); - } - // std::cerr << refDet << " == " << slu.determinant() << "\n"; - if (slu.solve(b, &x, SvTranspose)) { - VERIFY(b.isApprox(m2.transpose() * x, test_precision())); - } - - if (slu.solve(b, &x, SvAdjoint)) { - VERIFY(b.isApprox(m2.adjoint() * x, test_precision())); - } - - 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,UmfPack> slu(m2); - if (slu.succeeded()) { - if (slu.solve(b,&x)) { - if (count==0) { - VERIFY(refX.isApprox(x,test_precision()) && "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 + } +#endif + +namespace Eigen { + +/** \defgroup CholmodSupport_Module Cholmod Support module + * + * \nonstableyet + * + * + * \code + * #include + * \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 +#include +#include +#include +#include + +#ifdef EIGEN_GOOGLEHASH_SUPPORT + #include +#endif + +namespace Eigen { + +/** \defgroup SparseExtra_Module SparseExtra module + * + * This module contains some experimental features extending the sparse module. + * + * \code + * #include + * \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 +#include +#include + +namespace SuperLU_S { +#include +} +namespace SuperLU_D { +#include +} +namespace SuperLU_C { +#include +} +namespace SuperLU_Z { +#include +} +namespace Eigen { struct SluMatrix; } + +namespace Eigen { + +/** \defgroup SuperLUSupport_Module Super LU support + * + * \nonstableyet + * + * \code + * #include + * \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 +} +#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 + * \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 + +namespace Eigen { + +/** \defgroup Sparse_Module Sparse module + * + * \nonstableyet + * + * + * \code + * #include + * \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 similarity index 85% rename from Eigen/src/Sparse/CholmodSupport.h rename to 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 -cholmod_sparse SparseMatrixBase::asCholmodMatrix() +template +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::asCholmodMatrix() ei_cholmod_configure_matrix(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& mat) return res; } -template -MappedSparseMatrix::MappedSparseMatrix(cholmod_sparse& cm) +template +MappedSparseMatrix ei_map_cholmod(cholmod_sparse& cm) { - m_innerSize = cm.nrow; - m_outerSize = cm.ncol; - m_outerIndex = reinterpret_cast(cm.p); - m_innerIndices = reinterpret_cast(cm.i); - m_values = reinterpret_cast(cm.x); - m_nnz = m_outerIndex[cm.ncol]; + return MappedSparseMatrix + (cm.nrow, cm.ncol, reinterpret_cast(cm.p)[cm.ncol], + reinterpret_cast(cm.p), reinterpret_cast(cm.i),reinterpret_cast(cm.x) ); } template @@ -178,7 +175,7 @@ void SparseLLT::compute(const MatrixType& a) m_cholmodFactor = 0; } - cholmod_sparse A = const_cast(a).asCholmodMatrix(); + cholmod_sparse A = ei_asCholmodMatrix(const_cast(a)); // m_cholmod.supernodal = CHOLMOD_AUTO; // TODO // if (m_flags&IncompleteFactorization) @@ -209,7 +206,7 @@ SparseLLT::matrixL() const ei_assert(!(m_status & SupernodalFactorIsDirty)); cholmod_sparse* cmRes = cholmod_factor_to_sparse(m_cholmodFactor, &m_cholmod); - const_cast(m_matrix) = MappedSparseMatrix(*cmRes); + const_cast(m_matrix) = ei_map_cholmod(*cmRes); free(cmRes); m_status = (m_status & ~MatrixLIsDirty); @@ -235,7 +232,7 @@ bool SparseLLT::solveInPlace(MatrixBase &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::Map(reinterpret_cast(x->x),b.rows()); diff --git a/Eigen/src/Sparse/RandomSetter.h b/unsupported/Eigen/src/SparseExtra/RandomSetter.h similarity index 100% rename from Eigen/src/Sparse/RandomSetter.h rename to unsupported/Eigen/src/SparseExtra/RandomSetter.h diff --git a/Eigen/src/Sparse/SparseLDLT.h b/unsupported/Eigen/src/SparseExtra/SparseLDLT.h similarity index 100% rename from Eigen/src/Sparse/SparseLDLT.h rename to unsupported/Eigen/src/SparseExtra/SparseLDLT.h diff --git a/Eigen/src/Sparse/SparseLLT.h b/unsupported/Eigen/src/SparseExtra/SparseLLT.h similarity index 100% rename from Eigen/src/Sparse/SparseLLT.h rename to unsupported/Eigen/src/SparseExtra/SparseLLT.h diff --git a/Eigen/src/Sparse/SparseLU.h b/unsupported/Eigen/src/SparseExtra/SparseLU.h similarity index 100% rename from Eigen/src/Sparse/SparseLU.h rename to unsupported/Eigen/src/SparseExtra/SparseLU.h diff --git a/Eigen/src/Sparse/SuperLUSupport.h b/unsupported/Eigen/src/SparseExtra/SuperLUSupport.h similarity index 96% rename from Eigen/src/Sparse/SuperLUSupport.h rename to 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 > } }; -template -SluMatrix SparseMatrixBase::asSluMatrix() +template +SluMatrix ei_asSluMatrix(MatrixType& mat) { - return SluMatrix::Map(derived()); + return SluMatrix::Map(mat); } /** View a Super LU matrix as an Eigen expression */ -template -MappedSparseMatrix::MappedSparseMatrix(SluMatrix& sluMat) +template +MappedSparseMatrix 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(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( + sluMat.nrow, sluMat.ncol, sluMat.storage.outerInd[outerSize], + sluMat.storage.outerInd, sluMat.storage.innerInd, reinterpret_cast(sluMat.storage.values) ); } template @@ -401,7 +393,7 @@ void SparseLU::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 similarity index 84% rename from Eigen/src/Sparse/TaucsSupport.h rename to 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 -taucs_ccs_matrix SparseMatrixBase::asTaucsMatrix() +template +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::ret) res.flags |= TAUCS_INT; else if (ei_is_same_type::ret) @@ -63,15 +65,12 @@ taucs_ccs_matrix SparseMatrixBase::asTaucsMatrix() return res; } -template -MappedSparseMatrix::MappedSparseMatrix(taucs_ccs_matrix& taucsMat) +template +MappedSparseMatrix 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(taucsMat.values.v); - m_nnz = taucsMat.colptr[taucsMat.n]; + return MappedSparseMatrix + (taucsMat.m, taucsMat.n, taucsMat.colptr[taucsMat.n], + taucsMat.colptr, taucsMat.rowind, reinterpret_cast(taucsMat.values.v)); } template @@ -79,6 +78,7 @@ class SparseLLT : public SparseLLT { protected: typedef SparseLLT 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::compute(const MatrixType& a) m_taucsSupernodalFactor = 0; } - taucs_ccs_matrix taucsMatA = const_cast(a).asTaucsMatrix(); + taucs_ccs_matrix taucsMatA = ei_asTaucsMatrix(const_cast(a)); if (m_flags & IncompleteFactorization) { @@ -140,7 +140,7 @@ void SparseLLT::compute(const MatrixType& a) } // the matrix returned by Taucs is not necessarily sorted, // so let's copy it in two steps - DynamicSparseMatrix tmp = MappedSparseMatrix(*taucsRes); + DynamicSparseMatrix tmp = ei_map_taucs(*taucsRes); m_matrix = tmp; free(taucsRes); m_status = (m_status & ~(CompleteFactorization|MatrixLIsDirty)) @@ -176,7 +176,7 @@ SparseLLT::matrixL() const // the matrix returned by Taucs is not necessarily sorted, // so let's copy it in two steps - DynamicSparseMatrix tmp = MappedSparseMatrix(*taucsL); + DynamicSparseMatrix tmp = ei_map_taucs(*taucsL); const_cast(m_matrix) = tmp; free(taucsL); m_status = (m_status & ~MatrixLIsDirty); @@ -197,7 +197,7 @@ void SparseLLT::solveInPlace(MatrixBase &b) const } else if (m_flags & IncompleteFactorization) { - taucs_ccs_matrix taucsLLT = const_cast(m_matrix).asTaucsMatrix(); + taucs_ccs_matrix taucsLLT = ei_asTaucsMatrix(const_cast(m_matrix)); typename ei_plain_matrix_type::type x(b.rows()); for (int j=0; j