Improved and renamed arpack wrapper, added libshogun arpack example

shogun-toolbox · Aug 22, 2011 · 80149f4 · 80149f4
1 parent 876909b
commit 80149f4
Show file tree

Hide file tree

Showing 7 changed files with 85 additions and 22 deletions.
diff --git a/examples/undocumented/libshogun/Makefile b/examples/undocumented/libshogun/Makefile
@@ -34,7 +34,7 @@ TARGETS =	  basic_minimal classifier_libsvm classifier_minimal_svm \
 		  mathematics_confidence_intervals \
 		  clustering_kmeans base_parameter_map \
 		  splitting_stratified_crossvalidation \
-
+		  mathematics_arpack \
 
 all: $(TARGETS)
 

diff --git a/examples/undocumented/libshogun/mathematics_arpack.cpp b/examples/undocumented/libshogun/mathematics_arpack.cpp
@@ -0,0 +1,52 @@
+/*
+ * 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) 2011 Sergey Lisitsyn
+ * Copyright (C) 2011 Berlin Institute of Technology and Max-Planck-Society
+ */
+
+#include <shogun/base/init.h>
+#include <shogun/mathematics/arpack.h>
+
+using namespace shogun;
+
+int main(int argc, char** argv)
+{
+	init_shogun();
+
+	int N = 100;
+	int nev = 2;
+
+	double* matrix = new double[N*N];
+
+	double* eigenvalues = new double[nev];
+	double* eigenvectors = new double[nev*N];
+
+	for (int i=0; i<N; i++)
+	{
+		for (int j=0; j<N; j++)
+			matrix[i*N+j] = i*i+j*j;
+	}
+
+	int status = 0;
+	arpack_dsaeupd_wrap(matrix, NULL, N, 2, "LM", 1, false, 0.0, 0.0,
+	                    eigenvalues, eigenvectors, status);
+	if (status!=0)
+		return -1;
+
+	arpack_dsaeupd_wrap(matrix, NULL, N, 2, "BE", 3, false, 1.0, 0.0,
+	                    eigenvalues, eigenvectors, status);
+	if (status!=0)
+		return -1;
+
+
+	delete[] eigenvalues;
+	delete[] eigenvectors;
+	delete[] matrix;
+
+	exit_shogun();
+	return 0;
+}
diff --git a/src/shogun/mathematics/arpack.cpp b/src/shogun/mathematics/arpack.cpp
@@ -22,9 +22,10 @@ using namespace shogun;
 
 namespace shogun
 {
-void arpack_dsaupd(double* matrix, double* rhs_diag, int n, int nev, const char* which, 
-                   int mode, bool pos, double shift, double* eigenvalues, 
-                   double* eigenvectors, int& status)
+
+void arpack_dsaeupd_wrap(double* matrix, double* rhs_diag, int n, int nev, const char* which, 
+                         int mode, bool pos, double shift, double tolerance, double* eigenvalues, 
+                         double* eigenvectors, int& status)
 {
 	// loop vars
 	int i,j;
@@ -48,7 +49,7 @@ void arpack_dsaupd(double* matrix, double* rhs_diag, int n, int nev, const char*
 		bmat[0] = 'G';
 
 	// init tolerance (zero means machine precision)
-	double tol = 0.0;
+	double tol = tolerance;
 
 	// allocate array to hold residuals
 	double* resid = SG_MALLOC(double, n);
@@ -97,7 +98,7 @@ void arpack_dsaupd(double* matrix, double* rhs_diag, int n, int nev, const char*
 		// subtract shift from main diagonal if necessary
 		if (shift!=0.0)
 		{
-			SG_SDEBUG("Subtracting shift\n");
+			SG_SDEBUG("ARPACK: Subtracting shift\n");
 			// if right hand side diagonal matrix is provided
 			if (rhs_diag)
 				// subtract I*diag(rhs_diag)
@@ -113,18 +114,19 @@ void arpack_dsaupd(double* matrix, double* rhs_diag, int n, int nev, const char*
 		if (pos)
 		{
 			// with Cholesky
-			SG_SDEBUG("ARPACK: Using Cholesky factorization\n");
+			SG_SDEBUG("ARPACK: Using Cholesky factorization.\n");
 			clapack_dpotrf(CblasColMajor,CblasUpper,n,matrix,n);
 		}
 		else
 		{
 			// with LUP
-			SG_SDEBUG("ARPACK: Using LUP factorization\n");
+			SG_SDEBUG("ARPACK: Using LUP factorization.\n");
 			ipiv = SG_CALLOC(int, n);
 			clapack_dgetrf(CblasColMajor,n,n,matrix,n,ipiv);
 		}
 	}
 	// main computation loop
+	SG_SDEBUG("ARPACK: Starting main computation DSAUPD loop.\n");
 	do	 
 	{
 		dsaupd_(&ido, bmat, &n, which_, &nev, &tol, resid,
@@ -191,18 +193,18 @@ void arpack_dsaupd(double* matrix, double* rhs_diag, int n, int nev, const char*
 	if (info<0) 
 	{
 		if ((info<=-1)&&(info>=-6))
-			SG_SWARNING("DSAUPD failed. Wrong parameter passed.");
+			SG_SWARNING("ARPACK: DSAUPD failed. Wrong parameter passed.\n");
 		else if (info==-7)
-			SG_SWARNING("DSAUPD failed. Workaround array size is not sufficient.");
+			SG_SWARNING("ARPACK: DSAUPD failed. Workaround array size is not sufficient.\n");
 		else 
-			SG_SWARNING("DSAUPD failed. Error code: %d.", info);
+			SG_SWARNING("ARPACK: DSAUPD failed. Error code: %d.\n", info);
 
 		status = -1;
 	}
 	else 
 	{
 		if (info==1)
-			SG_SWARNING("Maximum number of iterations reached.\n");
+			SG_SWARNING("ARPACK: Maximum number of iterations reached.\n");
 
 		// allocate select for dseupd
 		int* select = SG_MALLOC(int, ncv);
@@ -217,6 +219,8 @@ void arpack_dsaupd(double* matrix, double* rhs_diag, int n, int nev, const char*
 		// specify that eigenvectors are going to be computed too		
 		int rvec = 1;
 
+		SG_SDEBUG("APRACK: Starting DSEUPD.\n");
+
 		// call dseupd_ routine
 		dseupd_(&rvec, all_, select, d, v, &ldv, &sigma, bmat,
 		        &n, which_, &nev, &tol, resid, &ncv, v, &ldv,
@@ -225,11 +229,13 @@ void arpack_dsaupd(double* matrix, double* rhs_diag, int n, int nev, const char*
 		// check for errors
 		if (ierr!=0)
 		{
-			SG_SWARNING("DSEUPD failed with status=%d", ierr);
+			SG_SWARNING("ARPACK: DSEUPD failed with status %d.\n", ierr);
 			status = -1;
 		}
 		else
 		{
+			SG_SDEBUG("ARPACK: Storing eigenpairs.\n");
+
 			// store eigenpairs to specified arrays
 			for (i=0; i<nev; i++)
 			{	

diff --git a/src/shogun/mathematics/arpack.h b/src/shogun/mathematics/arpack.h
@@ -10,6 +10,7 @@
 
 #ifndef ARPACK_H_
 #define ARPACK_H_
+#include <shogun/lib/config.h>
 #ifdef HAVE_ARPACK
 #ifdef HAVE_LAPACK
 #include <cblas.h>
@@ -59,16 +60,18 @@ namespace shogun
  * @param mode shift-mode of IRLM. Possible values:
  *        - 1: regular mode
  *        - 3: shift-invert mode
- * @param pos true if matrix is positive definite (Cholesky factorization)
+ * @param pos true if matrix is positive definite (Cholesky factorization is used in
+ *        this case instead of LUP factorization))
  * @param shift shift for shift-invert (3) mode of IRLM. In this mode
  *        routine will compute eigenvalues near provided shift
+ * @param tolerance tolerance with eigenvalues should be computed (zero means machine precision)
  * @param eigenvalues array of size nev to hold computed eigenvalues
  * @param eigenvectors array of size nev*n to hold computed eigenvectors
  * @param status on output -1 if computation failed
  */
-void arpack_dsaupd(double* matrix, double* rhs_diag, int n, int nev, 
-                   const char* which, int mode, bool pos, double shift,
-                   double* eigenvalues, double* eigenvectors, int& status);
+void arpack_dsaeupd_wrap(double* matrix, double* rhs_diag, int n, int nev, const char* which,
+                         int mode, bool pos, double shift, double tolerance,
+                         double* eigenvalues, double* eigenvectors, int& status);
 }
 #endif /* HAVE_LAPACK */
 #endif /* HAVE_ARPACK */

diff --git a/src/shogun/preprocessor/LaplacianEigenmaps.cpp b/src/shogun/preprocessor/LaplacianEigenmaps.cpp
@@ -149,7 +149,8 @@ SGMatrix<float64_t> CLaplacianEigenmaps::apply_to_feature_matrix(CFeatures* feat
 		// using ARPACK DS{E,A}UPD
 		int eigenproblem_status = 0;
 		float64_t* eigenvalues_vector = SG_MALLOC(float64_t,m_target_dim+1);
-		arpack_dsaupd(W_matrix,D_diag_vector,N,m_target_dim+1,"LA",3,false,0.0,eigenvalues_vector,W_matrix,eigenproblem_status);
+		arpack_dsaeupd_wrap(W_matrix,D_diag_vector,N,m_target_dim+1,"LA",3,false,0.0,0.0,
+		                    eigenvalues_vector,W_matrix,eigenproblem_status);
 		ASSERT(eigenproblem_status==0);
 		SG_FREE(eigenvalues_vector);
 	#else

diff --git a/src/shogun/preprocessor/LocallyLinearEmbedding.cpp b/src/shogun/preprocessor/LocallyLinearEmbedding.cpp
@@ -260,7 +260,8 @@ SGMatrix<float64_t> CLocallyLinearEmbedding::find_null_space(SGMatrix<float64_t>
 		// using ARPACK (faster)
 		eigenvalues_vector = SG_MALLOC(float64_t, dimension+1);
 		#ifdef HAVE_ARPACK
-		arpack_dsaupd(matrix.matrix,NULL,N,dimension+1,"LA",3,m_posdef,-1e-6,eigenvalues_vector,matrix.matrix,eigenproblem_status);
+		arpack_dsaeupd_wrap(matrix.matrix,NULL,N,dimension+1,"LA",3,m_posdef,-1e-6, 0.0,
+		                    eigenvalues_vector,matrix.matrix,eigenproblem_status);
 		#endif
 	}
 	else

diff --git a/src/shogun/preprocessor/MultidimensionalScaling.cpp b/src/shogun/preprocessor/MultidimensionalScaling.cpp
@@ -185,9 +185,9 @@ SGMatrix<float64_t> CMultidimensionalScaling::classic_embedding(SGMatrix<float64
 	// using ARPACK
 	float64_t* eigenvalues_vector = SG_MALLOC(float64_t, m_target_dim);
 	// solve eigenproblem with ARPACK (faster)
-	arpack_dsaupd(distance_matrix.matrix, NULL, N, m_target_dim, "LM", 1, false, 0.0,
-	              eigenvalues_vector, replace_feature_matrix,
-	              eigenproblem_status);
+	arpack_dsaeupd_wrap(distance_matrix.matrix, NULL, N, m_target_dim, "LM", 1, false, 0.0, 0.0,
+	                    eigenvalues_vector, replace_feature_matrix,
+	                    eigenproblem_status);
 	// check for failure
 	ASSERT(eigenproblem_status == 0);
 	// reverse eigenvectors order