Merge pull request #1321 from lambday/feature/log_determinant

lapack wrapper for tridiagonal eigensolver added

src/shogun/lib/computation/job/DenseExactLogJob.cpp

@@ -31,7 +31,8 @@ CDenseExactLogJob::CDenseExactLogJob()
 }
 
 CDenseExactLogJob::CDenseExactLogJob(CJobResultAggregator* aggregator,
-	CDenseMatrixOperator<float64_t>* log_operator, SGVector<float64_t> vector)
+	CDenseMatrixOperator<float64_t, float64_t>* log_operator,
+	SGVector<float64_t> vector)
 	: CIndependentJob(aggregator)
 {
 	init();

src/shogun/lib/computation/job/IndividualJobResultAggregator.cpp

@@ -31,7 +31,7 @@ CIndividualJobResultAggregator::CIndividualJobResultAggregator()
 }
 
 CIndividualJobResultAggregator::CIndividualJobResultAggregator(
-	CLinearOperator<float64_t>* linear_operator,
+	CLinearOperator<float64_t, float64_t>* linear_operator,
 	SGVector<float64_t> vector,
 	const float64_t& const_multiplier)
 	: CStoreVectorAggregator<complex64_t>(vector.vlen),

src/shogun/lib/computation/job/RationalApproximationIndividualJob.cpp

@@ -35,7 +35,7 @@ CRationalApproximationIndividualJob::CRationalApproximationIndividualJob()
 CRationalApproximationIndividualJob::CRationalApproximationIndividualJob(
 	CJobResultAggregator* aggregator,
 	CLinearSolver<complex64_t, float64_t>* linear_solver,
-	CLinearOperator<complex64_t>* linear_operator,
+	CLinearOperator<complex64_t, complex64_t>* linear_operator,
 	SGVector<float64_t> vector,
 	complex64_t weight)
 	: CIndependentJob(aggregator)

src/shogun/mathematics/lapack.cpp

@@ -38,6 +38,7 @@ using namespace shogun;
 #define DPOTRS dpotrs
 #define DGETRS dgetrs
 #define DSYGVX dsygvx
+#define DSTEMR dstemr
 #else
 #define DSYEV dsyev_
 #define DGESVD dgesvd_
@@ -52,6 +53,7 @@ using namespace shogun;
 #define DGETRS dgetrs_
 #define DPOTRS dpotrs_
 #define DSYGVX dsygvx_
+#define DSTEMR dstemr_
 #endif
 
 #ifndef HAVE_ATLAS
@@ -367,5 +369,34 @@ void wrap_dsygvx(int itype, char jobz, char uplo, int n, double *a, int lda, dou
 #endif
 }
 #undef DSYGVX
+
+void wrap_dstemr(char jobz, char range, int n, double* diag, double *subdiag,
+		double vl, double vu, int il, int iu, int* m, double* w, double* z__,
+		int ldz, int nzc, int *isuppz, int tryrac, int *info)
+{
+#ifdef HAVE_ACML
+	SG_SNOTIMPLEMENTED
+#else
+	int lwork=-1;
+	int liwork=-1;
+	double work1=0;
+	int iwork1=0;
+	DSTEMR(&jobz, &range, &n, diag, subdiag, &vl, &vu,
+		&il, &iu, m, w, z__, &ldz, &nzc, isuppz, &tryrac,
+		&work1, &lwork, &iwork1, &liwork, info);
+	lwork=(int)work1;
+	liwork=iwork1;
+	double* work=SG_MALLOC(double, lwork);
+	int* iwork=SG_MALLOC(int, liwork);
+	DSTEMR(&jobz, &range, &n, diag, subdiag, &vl, &vu,
+		&il, &iu, m, w, z__, &ldz, &nzc, isuppz, &tryrac,
+		work, &lwork, iwork, &liwork, info);
+
+	SG_FREE(work);
+	SG_FREE(iwork);
+#endif
+}
+#undef DSTEMR
+
 }
 #endif //HAVE_LAPACK

src/shogun/mathematics/lapack.h

@@ -74,6 +74,9 @@ void wrap_dsyevr(char jobz, char uplo, int n, double *a, int lda, int il, int iu
                  double *eigenvalues, double *eigenvectors, int *info);
 void wrap_dsygvx(int itype, char jobz, char uplo, int n, double *a, int lda, double *b,
                  int ldb, int il, int iu, double *eigenvalues, double *eigenvectors, int *info);
+void wrap_dstemr(char jobz, char range, int n, double* d__, double *e, double vl, double vu,
+		int il,	int iu, int* m, double* w, double* z__, int ldz, int nzc, int *isuppz,
+		int tryrac, int *info);
 #endif
 }
 
@@ -99,6 +102,9 @@ int dpotrs_(const char*, int*, int*, double*, int*, double*, int*, int*);
 int dsygvx_(int*, const char*, const char*, const char*, int*, double*, int*,
             double*, int*, double* , double*, int*, int*, double*,
             int*, double*, double*, int*, double*, int*, int*, int*, int*);
+int dstemr_(char*, char*, int*, double*, double*, double*, double*, int*, 
+	int*, int*, double*, double*, int*, int*, int*, int*, double*, 
+	int*, int*, int*, int*);
 #endif
 }
 

src/shogun/mathematics/logdet/DenseMatrixExactLog.cpp

@@ -36,7 +36,7 @@ CDenseMatrixExactLog::CDenseMatrixExactLog()
 }
 
 CDenseMatrixExactLog::CDenseMatrixExactLog(
-	CDenseMatrixOperator<float64_t>* op,
+	CDenseMatrixOperator<float64_t, float64_t>* op,
 	CIndependentComputationEngine* engine)
 	: COperatorFunction<float64_t>(
 		(CLinearOperator<float64_t>*)op, engine, OF_LOG)

src/shogun/mathematics/logdet/DirectEigenSolver.cpp

@@ -26,8 +26,8 @@ CDirectEigenSolver::CDirectEigenSolver()
 }
 
 CDirectEigenSolver::CDirectEigenSolver(
-	CDenseMatrixOperator<float64_t>* linear_operator)
-	: CEigenSolver((CLinearOperator<float64_t>*)linear_operator)
+	CDenseMatrixOperator<float64_t, float64_t>* linear_operator)
+	: CEigenSolver((CLinearOperator<float64_t, float64_t>*)linear_operator)
 {
 	SG_GCDEBUG("%s created (%p)\n", this->get_name(), this)
 }

src/shogun/mathematics/logdet/IterativeShiftedLinearFamilySolver.h

@@ -21,10 +21,10 @@ template <class T, class ST> class CLinearOperator;
 /** 
  * @brief abstract template base for CG based solvers to the solution of 
  * shifted linear systems of the form \f$(A+\sigma)x=b\f$ for several values
- * of \$f\sigma\f$ simultaneously, using only as many matrix-vector operations
+ * of \f$\sigma\f$ simultaneously, using only as many matrix-vector operations
  * as the solution of a single system requires. This class adds another
  * interface to the basic iterative linear solver that takes the shifts,
- * \$f\sigma\f$, and also weights, \f$\alpha\f$, and returns the summation
+ * \f$\sigma\f$, and also weights, \f$\alpha\f$, and returns the summation
  * \f$\sum_{i} \alpha_{i}x_{i}\f$, where \f$x_{i}\f$ is the solution of the
  * system \f$(A+\sigma_{i})x_{i}=b\f$.
  *
@@ -78,8 +78,8 @@ template<class T, class ST=T> class CIterativeShiftedLinearFamilySolver : public
 	 * -\zeta^{\sigma}_{n}+\zeta^{\sigma}_{n-1}\beta_{n-1}(1-\sigma\beta_{n}}\f$
 	 * [see Jergerlehner, eq 2.44]
 	 *
-	 * @param zeta_sh_old \f$\zeta^{\sigma}_{n-1} shifted params
-	 * @param zeta_sh_cur \f$\zeta^{\sigma}_{n} shifted params
+	 * @param zeta_sh_old \f$\zeta^{\sigma}_{n-1}\f$ shifted params
+	 * @param zeta_sh_cur \f$\zeta^{\sigma}_{n}\f$ shifted params
 	 * @param shifts \f$\sigma\f$ shifts
 	 * @param beta_old \f$\beta_{n-1}\f$, non-shifted
 	 * @param beta_cur \f$\beta_{n}\f$, non-shifted
@@ -94,8 +94,8 @@ template<class T, class ST=T> class CIterativeShiftedLinearFamilySolver : public
 	 * compute \f$\beta^{\sigma}_{n}\f$ as \f$\beta_{n}\frac{\zeta^{\sigma}_{n+1}}
 	 * {\zeta^{\sigma}_{n}}\f$
 	 *
-	 * @param zeta_sh_new \f$\zeta^{\sigma}_{n+1} shifted params
-	 * @param zeta_sh_cur \f$\zeta^{\sigma}_{n} shifted params
+	 * @param zeta_sh_new \f$\zeta^{\sigma}_{n+1}\f$ shifted params
+	 * @param zeta_sh_cur \f$\zeta^{\sigma}_{n}\f$ shifted params
 	 * @param beta_cur \f$\beta_{n}\f$, non-shifted
 	 * @param beta_sh \f$\beta^{\sigma}_{n}\f$, to be computed
 	 */
@@ -106,8 +106,8 @@ template<class T, class ST=T> class CIterativeShiftedLinearFamilySolver : public
 	 * compute \f$alpha^{\sigma}_{n}\f$ as \f$\alpha_{n}\frac{\zeta^{\sigma}
 	 * _{n}\beta^{\sigma}_{n-1}}{\zeta^{\sigma}_{n-1}\beta_{n-1}}\f$
 	 *
-	 * @param zeta_sh_cur \f$\zeta^{\sigma}_{n} shifted params
-	 * @param zeta_sh_old \f$\zeta^{\sigma}_{n-1} shifted params
+	 * @param zeta_sh_cur \f$\zeta^{\sigma}_{n}\f$ shifted params
+	 * @param zeta_sh_old \f$\zeta^{\sigma}_{n-1}\f$ shifted params
 	 * @param beta_sh \f$\beta^{\sigma}_{n}\f$, shifted params
 	 * @param beta \f$\beta_{n}\f$, non-shifted
 	 * @param alpha \f$\alpha_{n}\f$, non-shifted

src/shogun/mathematics/logdet/IterativeSolverIterator.h

@@ -113,7 +113,7 @@ typedef Matrix<T, Dynamic, 1> VectorXt;
 	/** maximum iteration limit */
 	const int64_t m_max_iteration_limit;
 
-	/* tolerence of the iterative solver */
+	/** tolerence of the iterative solver */
 	const float64_t m_tolerence;
 
 	/** true if converged successfully, false otherwise */

src/shogun/mathematics/logdet/LogRationalApproximationIndividual.cpp

@@ -32,7 +32,7 @@ CLogRationalApproximationIndividual::CLogRationalApproximationIndividual()
 }
 
 CLogRationalApproximationIndividual::CLogRationalApproximationIndividual(
-	CDenseMatrixOperator<float64_t>* linear_operator,
+	CDenseMatrixOperator<float64_t, float64_t>* linear_operator,
 	CIndependentComputationEngine* computation_engine,
 	CEigenSolver* eigen_solver, 
 	CLinearSolver<complex64_t, float64_t>* linear_solver,

tests/unit/mathematics/Lapack_wrapper_dstemr.cc

@@ -0,0 +1,64 @@
+/*
+ * This program is free software; 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 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Written (W) 2013 Soumyajit De
+ */
+
+#include <shogun/lib/common.h>
+
+#ifdef HAVE_LAPACK
+#ifdef HAVE_EIGEN3
+#ifdef HAVE_ATLAS
+#include <shogun/lib/SGVector.h>
+#include <shogun/mathematics/lapack.h>
+#include <shogun/mathematics/eigen3.h>
+#include <gtest/gtest.h>
+
+using namespace shogun;
+using namespace Eigen;
+
+TEST(Lapack_wrapper, dstemr)
+{
+	int32_t size=10;
+
+	// main diagonal of a fixed tridiagonal matrix
+	SGVector<float64_t> diag(size);
+	for (index_t i=0; i<size; ++i)
+		diag[i]=i+1;
+
+	// subdiagonal, fixed valued
+	SGVector<float64_t> subdiag(size);
+	subdiag.set_const(0.5);
+
+	int32_t M=0;
+	SGVector<float64_t> w(size);
+	int32_t tryrac=0.0;
+	int32_t info=0;
+
+	// computing all eigenvalues
+	wrap_dstemr('N', 'I', size, diag.vector, subdiag.vector,
+		0.0, 0.0, 1, size, &M, w.vector, NULL, 1, 1, NULL, tryrac, &info);
+
+	ASSERT(info==0);
+
+	// checking with eigen3 dense eigen solver
+	MatrixXd m=MatrixXd::Zero(size, size);
+	// filling the diagonal
+	for (index_t i=0; i<size; ++i)
+		m(i,i)=i+1;
+
+	// filling the subdiagonals
+	for (index_t i=0; i<size-1; ++i)
+		m(i,i+1)=m(i+1,i)=0.5;
+
+	VectorXcd eigenvals=m.eigenvalues();
+	Map<VectorXd> map(w.vector, w.vlen);
+
+	EXPECT_NEAR((map.cast<complex64_t>()-eigenvals).norm(), 0.0, 1E-10);
+}
+#endif // HAVE_ATLAS
+#endif // HAVE_EIGEN3
+#endif // HAVE_LAPACK