mirror of
https://gitlab.com/libeigen/eigen.git
synced 2025-04-19 08:09:36 +08:00
Sparse module:
* enable complex support for the CHOLMOD LLT backend using CHOLMOD's triangular solver * quick fix for complex support in SparseLLT::solve
This commit is contained in:
Gael Guennebaud
parent
361225068d
commit
ce3984844d
@ -25,6 +25,35 @@
|
||||
#ifndef EIGEN_CHOLMODSUPPORT_H
|
||||
#define EIGEN_CHOLMODSUPPORT_H
|
||||
|
||||
template<typename Scalar, typename CholmodType>
|
||||
void ei_cholmod_configure_matrix(CholmodType& mat)
|
||||
{
|
||||
if (ei_is_same_type<Scalar,float>::ret)
|
||||
{
|
||||
mat.xtype = CHOLMOD_REAL;
|
||||
mat.dtype = 1;
|
||||
}
|
||||
else if (ei_is_same_type<Scalar,double>::ret)
|
||||
{
|
||||
mat.xtype = CHOLMOD_REAL;
|
||||
mat.dtype = 0;
|
||||
}
|
||||
else if (ei_is_same_type<Scalar,std::complex<float> >::ret)
|
||||
{
|
||||
mat.xtype = CHOLMOD_COMPLEX;
|
||||
mat.dtype = 1;
|
||||
}
|
||||
else if (ei_is_same_type<Scalar,std::complex<double> >::ret)
|
||||
{
|
||||
mat.xtype = CHOLMOD_COMPLEX;
|
||||
mat.dtype = 0;
|
||||
}
|
||||
else
|
||||
{
|
||||
ei_assert(false && "Scalar type not supported by CHOLMOD");
|
||||
}
|
||||
}
|
||||
|
||||
template<typename Scalar, int Flags>
|
||||
cholmod_sparse SparseMatrix<Scalar,Flags>::asCholmodMatrix()
|
||||
{
|
||||
@ -42,30 +71,7 @@ cholmod_sparse SparseMatrix<Scalar,Flags>::asCholmodMatrix()
|
||||
res.dtype = 0;
|
||||
res.stype = -1;
|
||||
|
||||
if (ei_is_same_type<Scalar,float>::ret)
|
||||
{
|
||||
res.xtype = CHOLMOD_REAL;
|
||||
res.dtype = 1;
|
||||
}
|
||||
else if (ei_is_same_type<Scalar,double>::ret)
|
||||
{
|
||||
res.xtype = CHOLMOD_REAL;
|
||||
res.dtype = 0;
|
||||
}
|
||||
else if (ei_is_same_type<Scalar,std::complex<float> >::ret)
|
||||
{
|
||||
res.xtype = CHOLMOD_COMPLEX;
|
||||
res.dtype = 1;
|
||||
}
|
||||
else if (ei_is_same_type<Scalar,std::complex<double> >::ret)
|
||||
{
|
||||
res.xtype = CHOLMOD_COMPLEX;
|
||||
res.dtype = 0;
|
||||
}
|
||||
else
|
||||
{
|
||||
ei_assert(false && "Scalar type not supported by CHOLMOD");
|
||||
}
|
||||
ei_cholmod_configure_matrix<Scalar>(res);
|
||||
|
||||
if (Flags & SelfAdjoint)
|
||||
{
|
||||
@ -82,6 +88,24 @@ cholmod_sparse SparseMatrix<Scalar,Flags>::asCholmodMatrix()
|
||||
return res;
|
||||
}
|
||||
|
||||
template<typename Derived>
|
||||
cholmod_dense ei_cholmod_map_eigen_to_dense(MatrixBase<Derived>& mat)
|
||||
{
|
||||
typedef typename Derived::Scalar Scalar;
|
||||
|
||||
cholmod_dense res;
|
||||
res.nrow = mat.rows();
|
||||
res.ncol = mat.cols();
|
||||
res.nzmax = res.nrow * res.ncol;
|
||||
res.d = mat.derived().stride();
|
||||
res.x = mat.derived().data();
|
||||
res.z = 0;
|
||||
|
||||
ei_cholmod_configure_matrix<Scalar>(res);
|
||||
|
||||
return res;
|
||||
}
|
||||
|
||||
template<typename Scalar, int Flags>
|
||||
SparseMatrix<Scalar,Flags> SparseMatrix<Scalar,Flags>::Map(cholmod_sparse& cm)
|
||||
{
|
||||
@ -103,7 +127,7 @@ class SparseLLT<MatrixType,Cholmod> : public SparseLLT<MatrixType>
|
||||
{
|
||||
protected:
|
||||
typedef SparseLLT<MatrixType> Base;
|
||||
using Base::Scalar;
|
||||
using typename Base::Scalar;
|
||||
using Base::RealScalar;
|
||||
using Base::MatrixLIsDirty;
|
||||
using Base::SupernodalFactorIsDirty;
|
||||
@ -155,6 +179,7 @@ void SparseLLT<MatrixType,Cholmod>::compute(const MatrixType& a)
|
||||
}
|
||||
|
||||
cholmod_sparse A = const_cast<MatrixType&>(a).asCholmodMatrix();
|
||||
m_cholmod.supernodal = CHOLMOD_AUTO;
|
||||
// TODO
|
||||
if (m_flags&IncompleteFactorization)
|
||||
{
|
||||
@ -196,13 +221,19 @@ template<typename MatrixType>
|
||||
template<typename Derived>
|
||||
void SparseLLT<MatrixType,Cholmod>::solveInPlace(MatrixBase<Derived> &b) const
|
||||
{
|
||||
if (m_status & MatrixLIsDirty)
|
||||
matrixL();
|
||||
|
||||
const int size = m_matrix.rows();
|
||||
const int size = m_cholmodFactor->n;
|
||||
ei_assert(size==b.rows());
|
||||
|
||||
// this uses Eigen's triangular sparse solver
|
||||
if (m_status & MatrixLIsDirty)
|
||||
matrixL();
|
||||
Base::solveInPlace(b);
|
||||
// as long as our own triangular sparse solver is not fully optimal,
|
||||
// let's use CHOLMOD's one:
|
||||
// cholmod_dense cdb = ei_cholmod_map_eigen_to_dense(b);
|
||||
// cholmod_dense* x = cholmod_solve(CHOLMOD_LDLt, m_cholmodFactor, &cdb, &m_cholmod);
|
||||
// b = Matrix<typename Base::Scalar,Dynamic,1>::Map(reinterpret_cast<typename Base::Scalar*>(x->x),b.rows());
|
||||
// cholmod_free_dense(&x, &m_cholmod);
|
||||
}
|
||||
|
||||
#endif // EIGEN_CHOLMODSUPPORT_H
|
||||
|
||||
@ -190,7 +190,13 @@ bool SparseLLT<MatrixType, Backend>::solveInPlace(MatrixBase<Derived> &b) const
|
||||
ei_assert(size==b.rows());
|
||||
|
||||
m_matrix.solveTriangularInPlace(b);
|
||||
// FIXME should be .adjoint() but it fails to compile...
|
||||
// FIXME should be simply .adjoint() but it fails to compile...
|
||||
if (NumTraits<Scalar>::IsComplex)
|
||||
{
|
||||
CholMatrixType aux = m_matrix.conjugate();
|
||||
aux.transpose().solveTriangularInPlace(b);
|
||||
}
|
||||
else
|
||||
m_matrix.transpose().solveTriangularInPlace(b);
|
||||
|
||||
return true;
|
||||
|
||||
@ -148,7 +148,7 @@ class SparseMatrix
|
||||
*/
|
||||
inline void startFill(int reserveSize = 1000)
|
||||
{
|
||||
std::cerr << this << " startFill\n";
|
||||
// std::cerr << this << " startFill\n";
|
||||
setZero();
|
||||
m_data.reserve(reserveSize);
|
||||
}
|
||||
@ -214,7 +214,7 @@ class SparseMatrix
|
||||
m_data.index(id+1) = inner;
|
||||
//return (m_data.value(id+1) = 0);
|
||||
m_data.value(id+1) = 0;
|
||||
std::cerr << m_outerIndex[outer] << " " << m_outerIndex[outer+1] << "\n";
|
||||
// std::cerr << m_outerIndex[outer] << " " << m_outerIndex[outer+1] << "\n";
|
||||
return m_data.value(id+1);
|
||||
}
|
||||
|
||||
@ -222,7 +222,7 @@ class SparseMatrix
|
||||
|
||||
inline void endFill()
|
||||
{
|
||||
std::cerr << this << " endFill\n";
|
||||
// std::cerr << this << " endFill\n";
|
||||
int size = m_data.size();
|
||||
int i = m_outerSize;
|
||||
// find the last filled column
|
||||
@ -238,7 +238,7 @@ class SparseMatrix
|
||||
|
||||
void resize(int rows, int cols)
|
||||
{
|
||||
std::cerr << this << " resize " << rows << "x" << cols << "\n";
|
||||
// std::cerr << this << " resize " << rows << "x" << cols << "\n";
|
||||
const int outerSize = RowMajor ? rows : cols;
|
||||
m_innerSize = RowMajor ? cols : rows;
|
||||
m_data.clear();
|
||||
@ -273,6 +273,12 @@ class SparseMatrix
|
||||
*this = other.derived();
|
||||
}
|
||||
|
||||
inline SparseMatrix(const SparseMatrix& other)
|
||||
: m_outerSize(0), m_innerSize(0), m_outerIndex(0)
|
||||
{
|
||||
*this = other.derived();
|
||||
}
|
||||
|
||||
inline void swap(SparseMatrix& other)
|
||||
{
|
||||
//EIGEN_DBG_SPARSE(std::cout << "SparseMatrix:: swap\n");
|
||||
|
||||
@ -24,6 +24,27 @@
|
||||
|
||||
#include "sparse.h"
|
||||
|
||||
template<typename Scalar> void
|
||||
initSPD(double density,
|
||||
Matrix<Scalar,Dynamic,Dynamic>& refMat,
|
||||
SparseMatrix<Scalar>& sparseMat)
|
||||
{
|
||||
Matrix<Scalar,Dynamic,Dynamic> aux(refMat.rows(),refMat.cols());
|
||||
initSparse(density,refMat,sparseMat);
|
||||
refMat = refMat * refMat.adjoint();
|
||||
for (int k=0; k<2; ++k)
|
||||
{
|
||||
initSparse(density,aux,sparseMat,ForceNonZeroDiag);
|
||||
refMat += aux * aux.adjoint();
|
||||
}
|
||||
sparseMat.startFill();
|
||||
for (int j=0 ; j<sparseMat.cols(); ++j)
|
||||
for (int i=j ; i<sparseMat.rows(); ++i)
|
||||
if (refMat(i,j)!=Scalar(0))
|
||||
sparseMat.fill(i,j) = refMat(i,j);
|
||||
sparseMat.endFill();
|
||||
}
|
||||
|
||||
template<typename Scalar> void sparse_solvers(int rows, int cols)
|
||||
{
|
||||
double density = std::max(8./(rows*cols), 0.01);
|
||||
@ -56,7 +77,6 @@ template<typename Scalar> void sparse_solvers(int rows, int cols)
|
||||
}
|
||||
|
||||
// test LLT
|
||||
if (!NumTraits<Scalar>::IsComplex)
|
||||
{
|
||||
// TODO fix the issue with complex (see SparseLLT::solveInPlace)
|
||||
SparseMatrix<Scalar> m2(rows, cols);
|
||||
@ -65,20 +85,23 @@ template<typename Scalar> void sparse_solvers(int rows, int 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;
|
||||
initSPD(density, refMat2, m2);
|
||||
|
||||
refMat2.llt().solve(b, &refX);
|
||||
typedef SparseMatrix<Scalar,LowerTriangular|SelfAdjoint> SparseSelfAdjointMatrix;
|
||||
if (!NumTraits<Scalar>::IsComplex)
|
||||
{
|
||||
x = b;
|
||||
SparseLLT<SparseSelfAdjointMatrix> (m2).solveInPlace(x);
|
||||
//VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: default");
|
||||
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
|
||||
if (!NumTraits<Scalar>::IsComplex)
|
||||
{
|
||||
#ifdef EIGEN_TAUCS_SUPPORT
|
||||
x = b;
|
||||
SparseLLT<SparseSelfAdjointMatrix,Taucs>(m2,IncompleteFactorization).solveInPlace(x);
|
||||
@ -91,23 +114,25 @@ template<typename Scalar> void sparse_solvers(int rows, int cols)
|
||||
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)
|
||||
// TODO fix the issue with complex (see SparseLDLT::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);
|
||||
//initSPD(density, refMat2, m2);
|
||||
initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, 0, 0);
|
||||
refMat2 += refMat2.adjoint();
|
||||
refMat2.diagonal() *= 0.5;
|
||||
|
||||
refMat2.ldlt().solve(b, &refX);
|
||||
typedef SparseMatrix<Scalar,LowerTriangular|SelfAdjoint> SparseSelfAdjointMatrix;
|
||||
typedef SparseMatrix<Scalar,UpperTriangular|SelfAdjoint> SparseSelfAdjointMatrix;
|
||||
x = b;
|
||||
SparseLDLT<SparseSelfAdjointMatrix> ldlt(m2);
|
||||
if (ldlt.succeeded())
|
||||
@ -177,6 +202,6 @@ 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) );
|
||||
CALL_SUBTEST( sparse_solvers<double>(101, 101) );
|
||||
}
|
||||
}
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user