SVDLIBC support module for Eigen

SVDLIBCSupport.h

#ifndef EIGEN_SVDLIBCSUPPORT_H
#define EIGEN_SVDLIBCSUPPORT_H

#include <Eigen/Sparse>

namespace Eigen
{
	namespace svdlib_h
	{
		extern "C"
		{
#include <svdlib.h>
		}
	}
	
	class SVDLIBC
	{
	public:
		typedef double Scalar;
		typedef long Index;
		typedef SparseMatrix<Scalar, ColMajor, Index> MatrixType;
		
		typedef Matrix<Scalar, Dynamic, 1> VectorType;
		typedef Matrix<Scalar, Dynamic, Dynamic, ColMajor> DenseMatrixType;
		
		typedef Map<VectorType> SingularValuesType;
		typedef Map<DenseMatrixType> MatrixUType;
		typedef Map<DenseMatrixType> MatrixVType;
		
		SVDLIBC()
		: m_svdRec(0),
		  m_singularValues(0, 0),
		  m_matrixU(0, 0, 0),
		  m_matrixV(0, 0, 0)
		{
			_setDefault();
		}
		
		SVDLIBC(const MatrixType& A)
		: m_svdRec(0),
		  m_singularValues(0, 0),
		  m_matrixU(0, 0, 0),
		  m_matrixV(0, 0, 0)
		{
			_setDefault();
			compute(A);
		}
		
		~SVDLIBC()
		{
			if(m_svdRec)
				svdlib_h::svdFreeSVDRec(m_svdRec);
		}
		
		ComputationInfo info() const
		{
			return m_info;
		}
		
		SVDLIBC& compute(const MatrixType& A);
		
		int rank() const
		{
			eigen_assert(m_svdRec && "SVDLIBC is not initialized.");
			return m_svdRec->d;	
		}
		
		const SingularValuesType& singularValues() const
		{
			eigen_assert(m_svdRec && "SVDLIBC is not initialized.");
			return m_singularValues;
		}
		
		const MatrixUType& matrixU()
		{
			eigen_assert(m_svdRec && "SVDLIBC is not initialized.");
			return m_matrixU;
		}
		
		const MatrixVType& matrixV()
		{
			eigen_assert(m_svdRec && "SVDLIBC is not initialized.");
			return m_matrixV;
		}
		
		long maxIterations() const
		{
			return m_iterations;
		}
		
		SVDLIBC& setMaxIterations(long iterations)
		{
			m_iterations = iterations;
			return *this;
		}
		
		SVDLIBC& setMaxIterations(Default_t)
		{
			m_iterations = 0;
			return *this;
		}
		
		long dimensions() const
		{
			return m_dimensions;
		}
		
		SVDLIBC& setDimensions(long dimensions)
		{
			m_dimensions = dimensions;
			return *this;
		}
		
		SVDLIBC& setDimensions(Default_t)
		{
			m_dimensions = 0;
			return *this;
		}
		
		double threshold() const
		{
			return m_end[1];
		}
		
		SVDLIBC& setThreshold(double threshold)
		{
			m_end[0] = -threshold;
			m_end[1] = threshold;
			return *this;
		}
		
		SVDLIBC& setThreshold(Default_t)
		{
			m_end[0] = -1E-30;
			m_end[1] = +1E-30;
			return *this;
		}
		
		double kappa() const
		{
			return m_kappa;
		}
		
		SVDLIBC& setKappa(double kappa)
		{
			m_kappa = kappa;
			return *this;
		}
		
		SVDLIBC& setKappa(Default_t)
		{
			m_kappa = 1E-6;
			return *this;
		}
		
		static int verbosity()
		{
			return svdlib_h::SVDVerbosity;
		}
		
		static int setVerbosity(int verbosity)
		{
			return svdlib_h::SVDVerbosity = verbosity;
		}
		
	private:
		void _setDefault()
		{
			setDimensions(Default);
			setMaxIterations(Default);
			setThreshold(Default);
			setKappa(Default);
		}
		
		long m_dimensions;
		long m_iterations;
		double m_end[2];
		double m_kappa;
		
		ComputationInfo m_info;
		svdlib_h::SVDRec m_svdRec;
		
		SingularValuesType m_singularValues;
		MatrixUType m_matrixU;
		MatrixVType m_matrixV;
	};
	
	inline SVDLIBC& SVDLIBC::compute(const SVDLIBC::MatrixType& A)
	{
		eigen_assert(A.isCompressed() && "sparse matrix must be compressed");
		
		svdlib_h::smat mA;
		mA.rows = A.rows();
		mA.cols = A.cols();
		mA.vals = A.nonZeros();
		mA.pointr = const_cast<long*>(A.outerIndexPtr());
		mA.rowind = const_cast<long*>(A.innerIndexPtr());
		mA.value = const_cast<double*>(A.valuePtr());
		
		if(m_svdRec)
			svdlib_h::svdFreeSVDRec(m_svdRec);
		
		m_svdRec = svdlib_h::svdLAS2(&mA, m_dimensions, m_iterations, m_end, m_kappa);
		
		if(!m_svdRec)
		{
			m_info = NumericalIssue;
			return *this;
		}
		
		new (&m_singularValues) SingularValuesType(m_svdRec->S, m_svdRec->d);
		new (&m_matrixU) MatrixUType(m_svdRec->Ut->value[0], m_svdRec->Ut->cols, m_svdRec->Ut->rows);
		new (&m_matrixV) MatrixVType(m_svdRec->Vt->value[0], m_svdRec->Vt->cols, m_svdRec->Ut->rows);
		
		m_info = Success;
		
		return *this;
	}
}

#endif