From e1c50a3cb17a0375761aab55065aab85596d8407 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 20 Oct 2008 11:37:45 +0000 Subject: [PATCH] add unit tests for sparse LU and fix a couple of warnings --- Eigen/src/Sparse/AmbiVector.h | 2 +- Eigen/src/Sparse/LinkedVectorMatrix.h | 8 +++--- Eigen/src/Sparse/SuperLUSupport.h | 35 ++++++++++++--------------- test/CMakeLists.txt | 22 ++++++++++------- test/sparse.cpp | 31 +++++++++++++++++++++--- 5 files changed, 61 insertions(+), 37 deletions(-) diff --git a/Eigen/src/Sparse/AmbiVector.h b/Eigen/src/Sparse/AmbiVector.h index 44697bceb..8d5d88983 100644 --- a/Eigen/src/Sparse/AmbiVector.h +++ b/Eigen/src/Sparse/AmbiVector.h @@ -124,7 +124,7 @@ void AmbiVector::nonZeros() const template void AmbiVector::init(RealScalar estimatedDensity) { - if (m_mode = estimatedDensity>0.1) + if (estimatedDensity>0.1) init(IsDense); else init(IsSparse); diff --git a/Eigen/src/Sparse/LinkedVectorMatrix.h b/Eigen/src/Sparse/LinkedVectorMatrix.h index cb7d2120c..4e5a4862e 100644 --- a/Eigen/src/Sparse/LinkedVectorMatrix.h +++ b/Eigen/src/Sparse/LinkedVectorMatrix.h @@ -43,7 +43,7 @@ struct ei_traits > template struct LinkedVectorChunk { - LinkedVectorChunk() : size(0), next(0), prev(0) {} + LinkedVectorChunk() : next(0), prev(0), size(0) {} Element data[ChunkSize]; LinkedVectorChunk* next; LinkedVectorChunk* prev; @@ -141,7 +141,7 @@ class LinkedVectorMatrix inline void startFill(int reserveSize = 1000) { clear(); - for (int i=0; i > res.nrow = mat.rows(); res.ncol = mat.cols(); - res.lda = mat.stride(); - res.values = mat.data(); + res.storage.lda = mat.stride(); + res.storage.values = mat.data(); } }; @@ -157,10 +154,10 @@ struct SluMatrixMapHelper > res.Mtype = SLU_GE; - res.nnz = mat.nonZeros(); - res.values = mat._valuePtr(); - res.innerInd = mat._innerIndexPtr(); - res.outerInd = mat._outerIndexPtr(); + res.storage.nnz = mat.nonZeros(); + res.storage.values = mat._valuePtr(); + res.storage.innerInd = mat._innerIndexPtr(); + res.storage.outerInd = mat._outerIndexPtr(); res.setScalarType(); @@ -196,11 +193,11 @@ SparseMatrix SparseMatrix::Map(SluMatrix& sluMat) res.m_innerSize = sluMat.nrow; res.m_outerSize = sluMat.ncol; } - res.m_outerIndex = sluMat.outerInd; + res.m_outerIndex = sluMat.storage.outerInd; SparseArray data = SparseArray::Map( - sluMat.innerInd, - reinterpret_cast(sluMat.values), - sluMat.outerInd[res.m_outerSize]); + sluMat.storage.innerInd, + reinterpret_cast(sluMat.storage.values), + sluMat.storage.outerInd[res.m_outerSize]); res.m_data.swap(data); res.markAsRValue(); return res; @@ -295,9 +292,9 @@ void SparseLU::compute(const MatrixType& a) m_sluB.setStorageType(SLU_DN); m_sluB.setScalarType(); m_sluB.Mtype = SLU_GE; - m_sluB.values = 0; + m_sluB.storage.values = 0; m_sluB.nrow = m_sluB.ncol = 0; - m_sluB.lda = size; + m_sluB.storage.lda = size; m_sluX = m_sluB; StatInit(&m_sluStat); diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 146c560f3..6c6deee4c 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -23,15 +23,19 @@ if(CHOLMOD_FOUND) set(EXTERNAL_LIBS ${EXTERNAL_LIBS} ${CHOLMOD_LIBRARIES}) endif(CHOLMOD_FOUND) -# find_package(Umfpack) -# if(UMFPACK_FOUND) -# add_definitions("-DEIGEN_UMFPACK_SUPPORT") -# endif(UMFPACK_FOUND) -# -# find_package(SuperLU) -# if(SUPERLU_FOUND) -# add_definitions("-DEIGEN_SUPERLU_SUPPORT") -# endif(SUPERLU_FOUND) +find_package(Umfpack) +if(UMFPACK_FOUND) + add_definitions("-DEIGEN_UMFPACK_SUPPORT") + include_directories(${UMFPACK_INCLUDES}) + set(EXTERNAL_LIBS ${EXTERNAL_LIBS} ${UMFPACK_LIBRARIES}) +endif(UMFPACK_FOUND) + +find_package(SuperLU) +if(SUPERLU_FOUND) + add_definitions("-DEIGEN_SUPERLU_SUPPORT") + include_directories(${SUPERLU_INCLUDES}) + set(EXTERNAL_LIBS ${EXTERNAL_LIBS} ${SUPERLU_LIBRARIES}) +endif(SUPERLU_FOUND) if(CMAKE_COMPILER_IS_GNUCXX) if(CMAKE_SYSTEM_NAME MATCHES Linux) diff --git a/test/sparse.cpp b/test/sparse.cpp index ca80e6362..048a5e5cb 100644 --- a/test/sparse.cpp +++ b/test/sparse.cpp @@ -99,7 +99,7 @@ template void sparse(int rows, int cols) refMat.coeffRef(nonzeroCoords[0].x(), nonzeroCoords[0].y()) = Scalar(5); VERIFY_IS_APPROX(m, refMat); - #if 0 + // test InnerIterators and Block expressions for(int j=0; j void sparse(int rows, int cols) // TODO test row major } - #endif // test LLT { @@ -232,8 +231,6 @@ template void sparse(int rows, int cols) refMat2.diagonal() *= 0.5; refMat2.llt().solve(b, &refX); -// std::cerr << refMat2 << "\n\n" << refMat2.llt().matrixL() << "\n\n"; -// std::cerr << m2 << "\n\n"; typedef SparseMatrix SparseSelfAdjointMatrix; x = b; SparseLLT (m2).solveInPlace(x); @@ -256,6 +253,32 @@ template void sparse(int rows, int cols) #endif } + // test LU + { + 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); + + refMat2.lu().solve(b, &refX); +// x.setZero(); +// SparseLU > (m2).solve(b,&x); +// VERIFY(refX.isApprox(x,test_precision()) && "LU: default"); + #ifdef EIGEN_SUPERLU_SUPPORT + x.setZero(); + SparseLU,SuperLU>(m2).solve(b,&x); + VERIFY(refX.isApprox(x,test_precision()) && "LU: SuperLU"); + #endif + #ifdef EIGEN_UMFPACK_SUPPORT + x.setZero(); + SparseLU,UmfPack>(m2).solve(b,&x); + VERIFY(refX.isApprox(x,test_precision()) && "LU: umfpack"); + #endif + } + } void test_sparse()