aboutsummaryrefslogtreecommitdiffhomepage
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Eigen/Sparse1
-rw-r--r--Eigen/src/Sparse/AmbiVector.h35
-rw-r--r--Eigen/src/Sparse/CholmodSupport.h1
-rw-r--r--Eigen/src/Sparse/SparseLDLT.h346
-rw-r--r--Eigen/src/Sparse/SparseMatrix.h2
-rw-r--r--Eigen/src/Sparse/SparseMatrixBase.h2
-rw-r--r--Eigen/src/Sparse/SparseProduct.h1
-rw-r--r--Eigen/src/Sparse/TriangularSolver.h8
-rw-r--r--bench/sparse_cholesky.cpp4
-rw-r--r--cmake/FindCholmod.cmake4
-rw-r--r--cmake/FindSuperLU.cmake4
-rw-r--r--cmake/FindUmfpack.cmake11
-rw-r--r--test/CMakeLists.txt7
-rw-r--r--test/sparse.h92
-rw-r--r--test/sparse_basic.cpp (renamed from test/sparse.cpp)187
-rw-r--r--test/sparse_solvers.cpp180
16 files changed, 693 insertions, 192 deletions
diff --git a/Eigen/Sparse b/Eigen/Sparse
index 54b15b8ae..b8a10c9e0 100644
--- a/Eigen/Sparse
+++ b/Eigen/Sparse
@@ -81,6 +81,7 @@ namespace Eigen {
#include "src/Sparse/SparseProduct.h"
#include "src/Sparse/TriangularSolver.h"
#include "src/Sparse/SparseLLT.h"
+#include "src/Sparse/SparseLDLT.h"
#include "src/Sparse/SparseLU.h"
#ifdef EIGEN_CHOLMOD_SUPPORT
diff --git a/Eigen/src/Sparse/AmbiVector.h b/Eigen/src/Sparse/AmbiVector.h
index 8d5d88983..4834b42b7 100644
--- a/Eigen/src/Sparse/AmbiVector.h
+++ b/Eigen/src/Sparse/AmbiVector.h
@@ -36,7 +36,7 @@ template<typename _Scalar> class AmbiVector
typedef _Scalar Scalar;
typedef typename NumTraits<Scalar>::Real RealScalar;
AmbiVector(int size)
- : m_buffer(0), m_size(0), m_allocatedSize(0), m_mode(-1)
+ : m_buffer(0), m_size(0), m_allocatedSize(0), m_allocatedElements(0), m_mode(-1)
{
resize(size);
}
@@ -72,18 +72,35 @@ template<typename _Scalar> class AmbiVector
void reallocate(int size)
{
- Scalar* newBuffer = new Scalar[size/* *4 + (size * sizeof(int)*2)/sizeof(Scalar)+1 */];
- int copySize = std::min(size, m_size);
- memcpy(newBuffer, m_buffer, copySize * sizeof(Scalar));
+ // if the size of the matrix is not too large, let's allocate a bit more than needed such
+ // that we can handle dense vector even in sparse mode.
delete[] m_buffer;
- m_buffer = newBuffer;
- m_allocatedSize = size;
-
+ if (size<1000)
+ {
+ int allocSize = (size * sizeof(ListEl))/sizeof(Scalar);
+ m_allocatedElements = (allocSize*sizeof(Scalar))/sizeof(ListEl);
+ m_buffer = new Scalar[allocSize];
+ }
+ else
+ {
+ m_allocatedElements = (size*sizeof(Scalar))/sizeof(ListEl);
+ m_buffer = new Scalar[size];
+ }
m_size = size;
m_start = 0;
m_end = m_size;
}
+ void reallocateSparse()
+ {
+ int copyElements = m_allocatedElements;
+ m_allocatedElements = std::min(int(m_allocatedElements*1.5),m_size);
+ int allocSize = m_allocatedElements * sizeof(ListEl);
+ allocSize = allocSize/sizeof(Scalar) + (allocSize%sizeof(Scalar)>0?1:0);
+ Scalar* newBuffer = new Scalar[allocSize];
+ memcpy(newBuffer, m_buffer, copyElements * sizeof(ListEl));
+ }
+
protected:
// element type of the linked list
struct ListEl
@@ -99,6 +116,7 @@ template<typename _Scalar> class AmbiVector
int m_start;
int m_end;
int m_allocatedSize;
+ int m_allocatedElements;
int m_mode;
// linked list mode
@@ -219,6 +237,9 @@ Scalar& AmbiVector<Scalar>::coeffRef(int i)
}
else
{
+ if (m_llSize>=m_allocatedElements)
+ reallocateSparse();
+ ei_internal_assert(m_llSize<m_size && "internal error: overflow in sparse mode");
// let's insert a new coefficient
ListEl& el = llElements[m_llSize];
el.value = Scalar(0);
diff --git a/Eigen/src/Sparse/CholmodSupport.h b/Eigen/src/Sparse/CholmodSupport.h
index b6c13350c..c98ca2125 100644
--- a/Eigen/src/Sparse/CholmodSupport.h
+++ b/Eigen/src/Sparse/CholmodSupport.h
@@ -155,6 +155,7 @@ void SparseLLT<MatrixType,Cholmod>::compute(const MatrixType& a)
}
cholmod_sparse A = const_cast<MatrixType&>(a).asCholmodMatrix();
+ // TODO
if (m_flags&IncompleteFactorization)
{
m_cholmod.nmethods = 1;
diff --git a/Eigen/src/Sparse/SparseLDLT.h b/Eigen/src/Sparse/SparseLDLT.h
new file mode 100644
index 000000000..66e31aa64
--- /dev/null
+++ b/Eigen/src/Sparse/SparseLDLT.h
@@ -0,0 +1,346 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2008 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
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+/*
+
+NOTE: the _symbolic, and _numeric functions has been adapted from
+ the LDL library:
+
+LDL Copyright (c) 2005 by Timothy A. Davis. All Rights Reserved.
+
+LDL License:
+
+ Your use or distribution of LDL or any modified version of
+ LDL implies that you agree to this License.
+
+ This library is free software; you can redistribute it and/or
+ modify it under the terms of the GNU Lesser General Public
+ License as published by the Free Software Foundation; either
+ version 2.1 of the License, or (at your option) any later version.
+
+ This library is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public
+ License along with this library; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
+ USA
+
+ Permission is hereby granted to use or copy this program under the
+ terms of the GNU LGPL, provided that the Copyright, this License,
+ and the Availability of the original version is retained on all copies.
+ User documentation of any code that uses this code or any modified
+ version of this code must cite the Copyright, this License, the
+ Availability note, and "Used by permission." Permission to modify
+ the code and to distribute modified code is granted, provided the
+ Copyright, this License, and the Availability note are retained,
+ and a notice that the code was modified is included.
+ */
+
+#ifndef EIGEN_SPARSELDLT_H
+#define EIGEN_SPARSELDLT_H
+
+/** \ingroup Sparse_Module
+ *
+ * \class SparseLDLT
+ *
+ * \brief LDLT Cholesky decomposition of a sparse matrix and associated features
+ *
+ * \param MatrixType the type of the matrix of which we are computing the LDLT Cholesky decomposition
+ *
+ * \sa class LDLT, class LDLT
+ */
+template<typename MatrixType, int Backend = DefaultBackend>
+class SparseLDLT
+{
+ protected:
+ typedef typename MatrixType::Scalar Scalar;
+ typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
+ typedef SparseMatrix<Scalar,Lower|UnitDiagBit> CholMatrixType;
+ typedef Matrix<Scalar,MatrixType::ColsAtCompileTime,1> VectorType;
+
+ enum {
+ SupernodalFactorIsDirty = 0x10000,
+ MatrixLIsDirty = 0x20000
+ };
+
+ public:
+
+ /** Creates a dummy LDLT factorization object with flags \a flags. */
+ SparseLDLT(int flags = 0)
+ : m_flags(flags), m_status(0)
+ {
+ ei_assert((MatrixType::Flags&RowMajorBit)==0);
+ m_precision = RealScalar(0.1) * Eigen::precision<RealScalar>();
+ }
+
+ /** Creates a LDLT object and compute the respective factorization of \a matrix using
+ * flags \a flags. */
+ SparseLDLT(const MatrixType& matrix, int flags = 0)
+ : m_matrix(matrix.rows(), matrix.cols()), m_flags(flags), m_status(0)
+ {
+ ei_assert((MatrixType::Flags&RowMajorBit)==0);
+ m_precision = RealScalar(0.1) * Eigen::precision<RealScalar>();
+ compute(matrix);
+ }
+
+ /** Sets the relative threshold value used to prune zero coefficients during the decomposition.
+ *
+ * Setting a value greater than zero speeds up computation, and yields to an imcomplete
+ * factorization with fewer non zero coefficients. Such approximate factors are especially
+ * useful to initialize an iterative solver.
+ *
+ * \warning if precision is greater that zero, the LDLT factorization is not guaranteed to succeed
+ * even if the matrix is positive definite.
+ *
+ * Note that the exact meaning of this parameter might depends on the actual
+ * backend. Moreover, not all backends support this feature.
+ *
+ * \sa precision() */
+ void setPrecision(RealScalar v) { m_precision = v; }
+
+ /** \returns the current precision.
+ *
+ * \sa setPrecision() */
+ RealScalar precision() const { return m_precision; }
+
+ /** Sets the flags. Possible values are:
+ * - CompleteFactorization
+ * - IncompleteFactorization
+ * - MemoryEfficient (hint to use the memory most efficient method offered by the backend)
+ * - SupernodalMultifrontal (implies a complete factorization if supported by the backend,
+ * overloads the MemoryEfficient flags)
+ * - SupernodalLeftLooking (implies a complete factorization if supported by the backend,
+ * overloads the MemoryEfficient flags)
+ *
+ * \sa flags() */
+ void settagss(int f) { m_flags = f; }
+ /** \returns the current flags */
+ int flags() const { return m_flags; }
+
+ /** Computes/re-computes the LDLT factorization */
+ void compute(const MatrixType& matrix);
+
+ /** Perform a symbolic factorization */
+ void _symbolic(const MatrixType& matrix);
+ /** Perform the actual factorization using the previously
+ * computed symbolic factorization */
+ bool _numeric(const MatrixType& matrix);
+
+ /** \returns the lower triangular matrix L */
+ inline const CholMatrixType& matrixL(void) const { return m_matrix; }
+
+ /** \returns the coefficients of the diagonal matrix D */
+ inline VectorType vectorD(void) const { return m_diag; }
+
+ template<typename Derived>
+ bool solveInPlace(MatrixBase<Derived> &b) const;
+
+ /** \returns true if the factorization succeeded */
+ inline bool succeeded(void) const { return m_succeeded; }
+
+ protected:
+ CholMatrixType m_matrix;
+ VectorType m_diag;
+ VectorXi m_parent; // elimination tree
+ VectorXi m_nonZerosPerCol;
+// VectorXi m_w; // workspace
+ RealScalar m_precision;
+ int m_flags;
+ mutable int m_status;
+ bool m_succeeded;
+};
+
+/** Computes / recomputes the LDLT decomposition of matrix \a a
+ * using the default algorithm.
+ */
+template<typename MatrixType, int Backend>
+void SparseLDLT<MatrixType,Backend>::compute(const MatrixType& a)
+{
+ _symbolic(a);
+ m_succeeded = _numeric(a);
+}
+
+template<typename MatrixType, int Backend>
+void SparseLDLT<MatrixType,Backend>::_symbolic(const MatrixType& a)
+{
+ assert(a.rows()==a.cols());
+ const int size = a.rows();
+ m_matrix.resize(size, size);
+ m_parent.resize(size);
+ m_nonZerosPerCol.resize(size);
+ int * tags = ei_alloc_stack(int, size);
+
+ const int* Ap = a._outerIndexPtr();
+ const int* Ai = a._innerIndexPtr();
+ int* Lp = m_matrix._outerIndexPtr();
+ const int* P = 0;
+ int* Pinv = 0;
+
+ if (P)
+ {
+ /* If P is present then compute Pinv, the inverse of P */
+ for (int k = 0; k < size; k++)
+ Pinv[P[k]] = k;
+ }
+ for (int k = 0; k < size; k++)
+ {
+ /* L(k,:) pattern: all nodes reachable in etree from nz in A(0:k-1,k) */
+ m_parent[k] = -1; /* parent of k is not yet known */
+ tags[k] = k; /* mark node k as visited */
+ m_nonZerosPerCol[k] = 0; /* count of nonzeros in column k of L */
+ int kk = P ? P[k] : k; /* kth original, or permuted, column */
+ int p2 = Ap[kk+1];
+ for (int p = Ap[kk]; p < p2; p++)
+ {
+ /* A (i,k) is nonzero (original or permuted A) */
+ int i = Pinv ? Pinv[Ai[p]] : Ai[p];
+ if (i < k)
+ {
+ /* follow path from i to root of etree, stop at flagged node */
+ for (; tags[i] != k; i = m_parent[i])
+ {
+ /* find parent of i if not yet determined */
+ if (m_parent[i] == -1)
+ m_parent[i] = k;
+ m_nonZerosPerCol[i]++; /* L (k,i) is nonzero */
+ tags[i] = k; /* mark i as visited */
+ }
+ }
+ }
+ }
+ /* construct Lp index array from m_nonZerosPerCol column counts */
+ Lp[0] = 0;
+ for (int k = 0; k < size; k++)
+ Lp[k+1] = Lp[k] + m_nonZerosPerCol[k];
+
+ m_matrix.resizeNonZeros(Lp[size]);
+ ei_free_stack(tags, int, size);
+}
+
+template<typename MatrixType, int Backend>
+bool SparseLDLT<MatrixType,Backend>::_numeric(const MatrixType& a)
+{
+ assert(a.rows()==a.cols());
+ const int size = a.rows();
+ assert(m_parent.size()==size);
+ assert(m_nonZerosPerCol.size()==size);
+
+ const int* Ap = a._outerIndexPtr();
+ const int* Ai = a._innerIndexPtr();
+ const Scalar* Ax = a._valuePtr();
+ const int* Lp = m_matrix._outerIndexPtr();
+ int* Li = m_matrix._innerIndexPtr();
+ Scalar* Lx = m_matrix._valuePtr();
+ m_diag.resize(size);
+
+ Scalar * y = ei_alloc_stack(Scalar, size);
+ int * pattern = ei_alloc_stack(int, size);
+ int * tags = ei_alloc_stack(int, size);
+
+ const int* P = 0;
+ const int* Pinv = 0;
+ bool ok = true;
+
+ for (int k = 0; k < size; k++)
+ {
+ /* compute nonzero pattern of kth row of L, in topological order */
+ y[k] = 0.0; /* Y(0:k) is now all zero */
+ int top = size; /* stack for pattern is empty */
+ tags[k] = k; /* mark node k as visited */
+ m_nonZerosPerCol[k] = 0; /* count of nonzeros in column k of L */
+ int kk = (P) ? (P[k]) : (k); /* kth original, or permuted, column */
+ int p2 = Ap[kk+1];
+ for (int p = Ap[kk]; p < p2; p++)
+ {
+ int i = Pinv ? Pinv[Ai[p]] : Ai[p]; /* get A(i,k) */
+ if (i <= k)
+ {
+ y[i] += Ax[p]; /* scatter A(i,k) into Y (sum duplicates) */
+ int len;
+ for (len = 0; tags[i] != k; i = m_parent[i])
+ {
+ pattern[len++] = i; /* L(k,i) is nonzero */
+ tags[i] = k; /* mark i as visited */
+ }
+ while (len > 0)
+ pattern[--top] = pattern[--len];
+ }
+ }
+ /* compute numerical values kth row of L (a sparse triangular solve) */
+ m_diag[k] = y[k]; /* get D(k,k) and clear Y(k) */
+ y[k] = 0.0;
+ for (; top < size; top++)
+ {
+ int i = pattern[top]; /* pattern[top:n-1] is pattern of L(:,k) */
+ Scalar yi = y[i]; /* get and clear Y(i) */
+ y[i] = 0.0;
+ int p2 = Lp[i] + m_nonZerosPerCol[i];
+ int p;
+ for (p = Lp[i]; p < p2; p++)
+ y[Li[p]] -= Lx[p] * yi;
+ Scalar l_ki = yi / m_diag[i]; /* the nonzero entry L(k,i) */
+ m_diag[k] -= l_ki * yi;
+ Li[p] = k; /* store L(k,i) in column form of L */
+ Lx[p] = l_ki;
+ m_nonZerosPerCol[i]++; /* increment count of nonzeros in col i */
+ }
+ if (m_diag[k] == 0.0)
+ {
+ ok = false; /* failure, D(k,k) is zero */
+ break;
+ }
+ }
+
+ ei_free_stack(y, Scalar, size);
+ ei_free_stack(pattern, int, size);
+ ei_free_stack(tags, int, size);
+
+ return ok; /* success, diagonal of D is all nonzero */
+}
+
+/** Computes b = L^-T L^-1 b */
+template<typename MatrixType, int Backend>
+template<typename Derived>
+bool SparseLDLT<MatrixType, Backend>::solveInPlace(MatrixBase<Derived> &b) const
+{
+ const int size = m_matrix.rows();
+ ei_assert(size==b.rows());
+ if (!m_succeeded)
+ return false;
+
+ if (m_matrix.nonZeros()>0) // otherwise L==I
+ m_matrix.solveTriangularInPlace(b);
+ b = b.cwise() / m_diag;
+ // FIXME should be .adjoint() but it fails to compile...
+
+ if (m_matrix.nonZeros()>0) // otherwise L==I
+ m_matrix.transpose().solveTriangularInPlace(b);
+
+ return true;
+}
+
+#endif // EIGEN_SPARSELDLT_H
diff --git a/Eigen/src/Sparse/SparseMatrix.h b/Eigen/src/Sparse/SparseMatrix.h
index 1346b310a..f9540751e 100644
--- a/Eigen/src/Sparse/SparseMatrix.h
+++ b/Eigen/src/Sparse/SparseMatrix.h
@@ -213,7 +213,7 @@ class SparseMatrix
inline void swap(SparseMatrix& other)
{
- EIGEN_DBG_SPARSE(std::cout << "SparseMatrix:: swap\n");
+ //EIGEN_DBG_SPARSE(std::cout << "SparseMatrix:: swap\n");
std::swap(m_outerIndex, other.m_outerIndex);
std::swap(m_innerSize, other.m_innerSize);
std::swap(m_outerSize, other.m_outerSize);
diff --git a/Eigen/src/Sparse/SparseMatrixBase.h b/Eigen/src/Sparse/SparseMatrixBase.h
index 2695f23b1..9168d7cae 100644
--- a/Eigen/src/Sparse/SparseMatrixBase.h
+++ b/Eigen/src/Sparse/SparseMatrixBase.h
@@ -159,7 +159,7 @@ class SparseMatrixBase : public MatrixBase<Derived>
}
else
{
- LinkedVectorMatrix<Scalar, RowMajorBit> trans = m.derived();
+ SparseMatrix<Scalar, RowMajorBit> trans = m.derived();
s << trans;
}
}
diff --git a/Eigen/src/Sparse/SparseProduct.h b/Eigen/src/Sparse/SparseProduct.h
index c430ed928..f7e50542e 100644
--- a/Eigen/src/Sparse/SparseProduct.h
+++ b/Eigen/src/Sparse/SparseProduct.h
@@ -162,6 +162,7 @@ struct ei_sparse_product_selector<Lhs,Rhs,ResultType,ColMajor,ColMajor,ColMajor>
for (typename Rhs::InnerIterator rhsIt(rhs, j); rhsIt; ++rhsIt)
{
// FIXME should be written like this: tmp += rhsIt.value() * lhs.col(rhsIt.index())
+ tempVector.restart();
Scalar x = rhsIt.value();
for (typename Lhs::InnerIterator lhsIt(lhs, rhsIt.index()); lhsIt; ++lhsIt)
{
diff --git a/Eigen/src/Sparse/TriangularSolver.h b/Eigen/src/Sparse/TriangularSolver.h
index 8e71c936d..c32af328d 100644
--- a/Eigen/src/Sparse/TriangularSolver.h
+++ b/Eigen/src/Sparse/TriangularSolver.h
@@ -70,7 +70,9 @@ struct ei_solve_triangular_selector<Lhs,Rhs,Upper,RowMajor|IsSparse>
{
Scalar tmp = other.coeff(i,col);
typename Lhs::InnerIterator it(lhs, i);
- for(++it; it; ++it)
+ if (it.index() == i)
+ ++it;
+ for(; it; ++it)
{
tmp -= it.value() * other.coeff(it.index(),col);
}
@@ -107,7 +109,9 @@ struct ei_solve_triangular_selector<Lhs,Rhs,Lower,ColMajor|IsSparse>
other.coeffRef(i,col) /= it.value();
}
Scalar tmp = other.coeffRef(i,col);
- for(++it; it; ++it)
+ if (it.index()==i)
+ ++it;
+ for(; it; ++it)
other.coeffRef(it.index(), col) -= tmp * it.value();
}
}
diff --git a/bench/sparse_cholesky.cpp b/bench/sparse_cholesky.cpp
index 0a17848a0..5e07fac9b 100644
--- a/bench/sparse_cholesky.cpp
+++ b/bench/sparse_cholesky.cpp
@@ -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[])
@@ -118,7 +118,7 @@ int main(int argc, char *argv[])
if (!ei_isMuchSmallerThan(ei_abs(chol.matrixL()(i,j)), 0.1))
count++;
std::cout << "dense: " << "nnz = " << count << "\n";
- std::cout << "dense:\n" << m1 << "\n\n" << chol.matrixL() << endl;
+// std::cout << "dense:\n" << m1 << "\n\n" << chol.matrixL() << endl;
}
#endif
diff --git a/cmake/FindCholmod.cmake b/cmake/FindCholmod.cmake
index 444279ab9..3b9d5d3dc 100644
--- a/cmake/FindCholmod.cmake
+++ b/cmake/FindCholmod.cmake
@@ -61,6 +61,10 @@ if(CHOLMOD_LIBRARIES)
endif(CHOLMOD_LIBRARIES)
+if(CHOLMOD_LIBRARIES AND CMAKE_COMPILER_IS_GNUCXX)
+ set(CHOLMOD_LIBRARIES ${CHOLMOD_LIBRARIES} -lgfortran)
+endif(CHOLMOD_LIBRARIES AND CMAKE_COMPILER_IS_GNUCXX)
+
include(FindPackageHandleStandardArgs)
find_package_handle_standard_args(CHOLMOD DEFAULT_MSG
CHOLMOD_INCLUDES CHOLMOD_LIBRARIES)
diff --git a/cmake/FindSuperLU.cmake b/cmake/FindSuperLU.cmake
index 9d12ebd77..7383754e2 100644
--- a/cmake/FindSuperLU.cmake
+++ b/cmake/FindSuperLU.cmake
@@ -13,6 +13,10 @@ find_path(SUPERLU_INCLUDES
find_library(SUPERLU_LIBRARIES superlu PATHS $ENV{SUPERLUDIR} ${LIB_INSTALL_DIR})
+if(SUPERLU_LIBRARIES AND CMAKE_COMPILER_IS_GNUCXX)
+ set(SUPERLU_LIBRARIES ${SUPERLU_LIBRARIES} -lgfortran)
+endif(SUPERLU_LIBRARIES AND CMAKE_COMPILER_IS_GNUCXX)
+
include(FindPackageHandleStandardArgs)
find_package_handle_standard_args(SUPERLU DEFAULT_MSG
SUPERLU_INCLUDES SUPERLU_LIBRARIES)
diff --git a/cmake/FindUmfpack.cmake b/cmake/FindUmfpack.cmake
index ae7883778..4f6875a32 100644
--- a/cmake/FindUmfpack.cmake
+++ b/cmake/FindUmfpack.cmake
@@ -3,6 +3,11 @@ if (UMFPACK_INCLUDES AND UMFPACK_LIBRARIES)
set(UMFPACK_FIND_QUIETLY TRUE)
endif (UMFPACK_INCLUDES AND UMFPACK_LIBRARIES)
+enable_language(Fortran)
+find_package(BLAS)
+
+if(BLAS_FOUND)
+
find_path(UMFPACK_INCLUDES
NAMES
umfpack.h
@@ -39,6 +44,12 @@ if(UMFPACK_LIBRARIES)
endif(UMFPACK_LIBRARIES)
+if(UMFPACK_LIBRARIES)
+ set(UMFPACK_LIBRARIES ${UMFPACK_LIBRARIES} ${BLAS_LIBRARIES})
+endif(UMFPACK_LIBRARIES)
+
+endif(BLAS_FOUND)
+
include(FindPackageHandleStandardArgs)
find_package_handle_standard_args(UMFPACK DEFAULT_MSG
UMFPACK_INCLUDES UMFPACK_LIBRARIES)
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index b902125c8..bec17e621 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -158,18 +158,19 @@ ei_add_test(smallvectors)
ei_add_test(map)
ei_add_test(array)
ei_add_test(triangular)
-ei_add_test(cholesky " " ${GSL_LIBRARIES})
+ei_add_test(cholesky " " "${GSL_LIBRARIES}")
ei_add_test(lu ${EI_OFLAG})
ei_add_test(determinant)
ei_add_test(inverse)
ei_add_test(qr)
-ei_add_test(eigensolver " " ${GSL_LIBRARIES})
+ei_add_test(eigensolver " " "${GSL_LIBRARIES}")
ei_add_test(svd)
ei_add_test(geometry)
ei_add_test(hyperplane)
ei_add_test(parametrizedline)
ei_add_test(alignedbox)
ei_add_test(regression)
-ei_add_test(sparse " " ${SPARSE_LIBS})
+ei_add_test(sparse_basic " " "${SPARSE_LIBS}")
+ei_add_test(sparse_solvers " " "${SPARSE_LIBS}")
endif(BUILD_TESTS)
diff --git a/test/sparse.h b/test/sparse.h
new file mode 100644
index 000000000..bf4a2ac04
--- /dev/null
+++ b/test/sparse.h
@@ -0,0 +1,92 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2008 Daniel Gomez Ferro <dgomezferro@gmail.com>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_TESTSPARSE_H
+
+#ifdef __GNUC__
+#include <ext/hash_map>
+#endif
+
+#ifdef EIGEN_GOOGLEHASH_SUPPORT
+ #include <google/sparse_hash_map>
+#endif
+
+#include "main.h"
+#include <Eigen/Cholesky>
+#include <Eigen/LU>
+#include <Eigen/Sparse>
+
+enum {
+ ForceNonZeroDiag = 1,
+ MakeLowerTriangular = 2,
+ MakeUpperTriangular = 4
+};
+
+/* Initializes both a sparse and dense matrix with same random values,
+ * and a ratio of \a density non zero entries.
+ * \param flags is a union of ForceNonZeroDiag, MakeLowerTriangular and MakeUpperTriangular
+ * allowing to control the shape of the matrix.
+ * \param zeroCoords and nonzeroCoords allows to get the coordinate lists of the non zero,
+ * and zero coefficients respectively.
+ */
+template<typename Scalar> void
+initSparse(double density,
+ Matrix<Scalar,Dynamic,Dynamic>& refMat,
+ SparseMatrix<Scalar>& sparseMat,
+ int flags = 0,
+ std::vector<Vector2i>* zeroCoords = 0,
+ std::vector<Vector2i>* nonzeroCoords = 0)
+{
+ sparseMat.startFill(refMat.rows()*refMat.cols()*density);
+ for(int j=0; j<refMat.cols(); j++)
+ {
+ for(int i=0; i<refMat.rows(); i++)
+ {
+ Scalar v = (ei_random<double>(0,1) < density) ? ei_random<Scalar>() : Scalar(0);
+ if ((flags&ForceNonZeroDiag) && (i==j))
+ {
+ v = ei_random<Scalar>()*Scalar(3.);
+ v = v*v + Scalar(5.);
+ }
+ if ((flags & MakeLowerTriangular) && j>i)
+ v = Scalar(0);
+ else if ((flags & MakeUpperTriangular) && j<i)
+ v = Scalar(0);
+ if (v!=Scalar(0))
+ {
+ sparseMat.fill(i,j) = v;
+ if (nonzeroCoords)
+ nonzeroCoords->push_back(Vector2i(i,j));
+ }
+ else if (zeroCoords)
+ {
+ zeroCoords->push_back(Vector2i(i,j));
+ }
+ refMat(i,j) = v;
+ }
+ }
+ sparseMat.endFill();
+}
+
+#endif // EIGEN_TESTSPARSE_H
diff --git a/test/sparse.cpp b/test/sparse_basic.cpp
index 949aed940..c890d9528 100644
--- a/test/sparse.cpp
+++ b/test/sparse_basic.cpp
@@ -22,63 +22,7 @@
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
-#ifdef __GNUC__
-#include <ext/hash_map>
-#endif
-
-#ifdef EIGEN_GOOGLEHASH_SUPPORT
- #include <google/sparse_hash_map>
-#endif
-
-#include "main.h"
-#include <Eigen/Cholesky>
-#include <Eigen/LU>
-#include <Eigen/Sparse>
-
-enum {
- ForceNonZeroDiag = 1,
- MakeLowerTriangular = 2,
- MakeUpperTriangular = 4
-};
-
-template<typename Scalar> void
-initSparse(double density,
- Matrix<Scalar,Dynamic,Dynamic>& refMat,
- SparseMatrix<Scalar>& sparseMat,
- int flags = 0,
- std::vector<Vector2i>* zeroCoords = 0,
- std::vector<Vector2i>* nonzeroCoords = 0)
-{
- sparseMat.startFill(refMat.rows()*refMat.cols()*density);
- for(int j=0; j<refMat.cols(); j++)
- {
- for(int i=0; i<refMat.rows(); i++)
- {
- Scalar v = (ei_random<double>(0,1) < density) ? ei_random<Scalar>() : Scalar(0);
- if ((flags&ForceNonZeroDiag) && (i==j))
- {
- v = ei_random<Scalar>()*Scalar(3.);
- v = v*v + Scalar(5.);
- }
- if ((flags & MakeLowerTriangular) && j>i)
- v = Scalar(0);
- else if ((flags & MakeUpperTriangular) && j<i)
- v = Scalar(0);
- if (v!=Scalar(0))
- {
- sparseMat.fill(i,j) = v;
- if (nonzeroCoords)
- nonzeroCoords->push_back(Vector2i(i,j));
- }
- else if (zeroCoords)
- {
- zeroCoords->push_back(Vector2i(i,j));
- }
- refMat(i,j) = v;
- }
- }
- sparseMat.endFill();
-}
+#include "sparse.h"
template<typename SetterType,typename DenseType, typename SparseType>
bool test_random_setter(SparseType& sm, const DenseType& ref, const std::vector<Vector2i>& nonzeroCoords)
@@ -98,7 +42,7 @@ bool test_random_setter(SparseType& sm, const DenseType& ref, const std::vector<
return sm.isApprox(ref);
}
-template<typename Scalar> void sparse(int rows, int cols)
+template<typename Scalar> void sparse_basic(int rows, int cols)
{
double density = std::max(8./(rows*cols), 0.01);
typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
@@ -113,8 +57,8 @@ template<typename Scalar> void sparse(int rows, int cols)
std::vector<Vector2i> nonzeroCoords;
initSparse<Scalar>(density, refMat, m, 0, &zeroCoords, &nonzeroCoords);
- VERIFY(zeroCoords.size()>0 && "re-run the test");
- VERIFY(nonzeroCoords.size()>0 && "re-run the test");
+ if (zeroCoords.size()==0 || nonzeroCoords.size()==0)
+ return;
// test coeff and coeffRef
for (int i=0; i<(int)zeroCoords.size(); ++i)
@@ -128,7 +72,7 @@ template<typename Scalar> void sparse(int rows, int cols)
refMat.coeffRef(nonzeroCoords[0].x(), nonzeroCoords[0].y()) = Scalar(5);
VERIFY_IS_APPROX(m, refMat);
-
+/*
// test InnerIterators and Block expressions
for (int t=0; t<10; ++t)
{
@@ -167,6 +111,7 @@ template<typename Scalar> void sparse(int rows, int cols)
VERIFY_IS_APPROX(m.row(r) + m.row(r), (m + m).row(r));
VERIFY_IS_APPROX(m.row(r) + m.row(r), refMat.row(r) + refMat.row(r));
}
+ */
// test SparseSetters
// coherent setter
@@ -234,7 +179,7 @@ template<typename Scalar> 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);
@@ -251,123 +196,13 @@ template<typename Scalar> void sparse(int rows, int cols)
VERIFY_IS_APPROX(m4=m2.transpose()*m3.transpose(), refMat4=refMat2.transpose()*refMat3.transpose());
VERIFY_IS_APPROX(m4=m2*m3.transpose(), refMat4=refMat2*refMat3.transpose());
}
-
- // test triangular solver
- {
- DenseVector vec2 = vec1, vec3 = vec1;
- SparseMatrix<Scalar> m2(rows, cols);
- DenseMatrix refMat2 = DenseMatrix::Zero(rows, cols);
-
- // lower
- initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords);
- VERIFY_IS_APPROX(refMat2.template marked<Lower>().solveTriangular(vec2),
- m2.template marked<Lower>().solveTriangular(vec3));
-
- // upper
- initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, &zeroCoords, &nonzeroCoords);
- VERIFY_IS_APPROX(refMat2.template marked<Upper>().solveTriangular(vec2),
- m2.template marked<Upper>().solveTriangular(vec3));
-
- // TODO test row major
- }
-
- // test LLT
- if (!NumTraits<Scalar>::IsComplex)
- {
- // 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, &zeroCoords, &nonzeroCoords);
- refMat2 += refMat2.adjoint();
- refMat2.diagonal() *= 0.5;
-
- refMat2.llt().solve(b, &refX);
- typedef SparseMatrix<Scalar,Lower|SelfAdjoint> SparseSelfAdjointMatrix;
- x = b;
- SparseLLT<SparseSelfAdjointMatrix> (m2).solveInPlace(x);
- //VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: default");
- #ifdef EIGEN_CHOLMOD_SUPPORT
- x = b;
- SparseLLT<SparseSelfAdjointMatrix,Cholmod>(m2).solveInPlace(x);
- VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: cholmod");
- #endif
- #ifdef EIGEN_TAUCS_SUPPORT
- x = b;
- SparseLLT<SparseSelfAdjointMatrix,Taucs>(m2,IncompleteFactorization).solveInPlace(x);
- VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (IncompleteFactorization)");
- x = b;
- SparseLLT<SparseSelfAdjointMatrix,Taucs>(m2,SupernodalMultifrontal).solveInPlace(x);
- VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (SupernodalMultifrontal)");
- x = b;
- SparseLLT<SparseSelfAdjointMatrix,Taucs>(m2,SupernodalLeftLooking).solveInPlace(x);
- VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (SupernodalLeftLooking)");
- #endif
- }
-
- // 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);
-
- LU<DenseMatrix> refLu(refMat2);
- refLu.solve(b, &refX);
- Scalar refDet = refLu.determinant();
- 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 (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++;
- }
-#endif
}
-void test_sparse()
+void test_sparse_basic()
{
for(int i = 0; i < g_repeat; i++) {
- CALL_SUBTEST( sparse<double>(8, 8) );
- CALL_SUBTEST( sparse<std::complex<double> >(16, 16) );
- CALL_SUBTEST( sparse<double>(33, 33) );
+ CALL_SUBTEST( sparse_basic<double>(8, 8) );
+ CALL_SUBTEST( sparse_basic<std::complex<double> >(16, 16) );
+ CALL_SUBTEST( sparse_basic<double>(33, 33) );
}
}
diff --git a/test/sparse_solvers.cpp b/test/sparse_solvers.cpp
new file mode 100644
index 000000000..b73494879
--- /dev/null
+++ b/test/sparse_solvers.cpp
@@ -0,0 +1,180 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2008 Daniel Gomez Ferro <dgomezferro@gmail.com>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#include "sparse.h"
+
+template<typename Scalar> void sparse_solvers(int rows, int cols)
+{
+ double density = std::max(8./(rows*cols), 0.01);
+ typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
+ typedef Matrix<Scalar,Dynamic,1> DenseVector;
+ Scalar eps = 1e-6;
+
+ DenseVector vec1 = DenseVector::Random(rows);
+
+ std::vector<Vector2i> zeroCoords;
+ std::vector<Vector2i> nonzeroCoords;
+
+ // test triangular solver
+ {
+ DenseVector vec2 = vec1, vec3 = vec1;
+ SparseMatrix<Scalar> m2(rows, cols);
+ DenseMatrix refMat2 = DenseMatrix::Zero(rows, cols);
+
+ // lower
+ initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords);
+ VERIFY_IS_APPROX(refMat2.template marked<Lower>().solveTriangular(vec2),
+ m2.template marked<Lower>().solveTriangular(vec3));
+
+ // upper
+ initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, &zeroCoords, &nonzeroCoords);
+ VERIFY_IS_APPROX(refMat2.template marked<Upper>().solveTriangular(vec2),
+ m2.template marked<Upper>().solveTriangular(vec3));
+
+ // TODO test row major
+ }
+
+ // test LLT
+ if (!NumTraits<Scalar>::IsComplex)
+ {
+ // 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, &zeroCoords, &nonzeroCoords);
+ refMat2 += refMat2.adjoint();
+ refMat2.diagonal() *= 0.5;
+
+ refMat2.llt().solve(b, &refX);
+ typedef SparseMatrix<Scalar,Lower|SelfAdjoint> SparseSelfAdjointMatrix;
+ x = b;
+ SparseLLT<SparseSelfAdjointMatrix> (m2).solveInPlace(x);
+ //VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: default");
+ #ifdef EIGEN_CHOLMOD_SUPPORT
+ x = b;
+ SparseLLT<SparseSelfAdjointMatrix,Cholmod>(m2).solveInPlace(x);
+ VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: cholmod");
+ #endif
+ #ifdef EIGEN_TAUCS_SUPPORT
+ x = b;
+ SparseLLT<SparseSelfAdjointMatrix,Taucs>(m2,IncompleteFactorization).solveInPlace(x);
+ VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (IncompleteFactorization)");
+ x = b;
+ SparseLLT<SparseSelfAdjointMatrix,Taucs>(m2,SupernodalMultifrontal).solveInPlace(x);
+ VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (SupernodalMultifrontal)");
+ x = b;
+ SparseLLT<SparseSelfAdjointMatrix,Taucs>(m2,SupernodalLeftLooking).solveInPlace(x);
+ VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (SupernodalLeftLooking)");
+ #endif
+ }
+
+ // test LDLT
+ if (!NumTraits<Scalar>::IsComplex)
+ {
+ // TODO fix the issue with complex (see SparseLDT::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|MakeUpperTriangular, &zeroCoords, &nonzeroCoords);
+ refMat2 += refMat2.adjoint();
+ refMat2.diagonal() *= 0.5;
+
+ refMat2.ldlt().solve(b, &refX);
+ typedef SparseMatrix<Scalar,Lower|SelfAdjoint> SparseSelfAdjointMatrix;
+ x = b;
+ SparseLDLT<SparseSelfAdjointMatrix> ldlt(m2);
+ if (ldlt.succeeded())
+ ldlt.solveInPlace(x);
+ 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);
+
+ LU<DenseMatrix> refLu(refMat2);
+ refLu.solve(b, &refX);
+ Scalar refDet = refLu.determinant();
+ 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 (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()
+{
+ for(int i = 0; i < g_repeat; i++) {
+ CALL_SUBTEST( sparse_solvers<double>(8, 8) );
+ CALL_SUBTEST( sparse_solvers<std::complex<double> >(16, 16) );
+ CALL_SUBTEST( sparse_solvers<double>(33, 33) );
+ }
+}
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-23 04:15:16 +0000