added test for new SGVector/SGMatrix based eigenvalue routine and a

check for correct results for a simple example.
shogun-toolbox · May 29, 2012 · bea615a · bea615a
1 parent 80e75c8
commit bea615a
Showing 1 changed file with 64 additions and 7 deletions.
diff --git a/examples/undocumented/libshogun/mathematics_lapack.cpp b/examples/undocumented/libshogun/mathematics_lapack.cpp
@@ -5,19 +5,67 @@
  * (at your option) any later version.
  *
  * Written (W) 2011 Sergey Lisitsyn
+ * Written (W) 2012 Heiko Strathmann
  * Copyright (C) 2011 Berlin Institute of Technology and Max-Planck-Society
  */
 
 #include <shogun/base/init.h>
 #include <shogun/lib/config.h>
+#include <shogun/mathematics/Math.h>
 #include <shogun/mathematics/lapack.h>
 
 using namespace shogun;
 
-int main(int argc, char** argv)
+bool is_equal(float64_t a, float64_t b, float64_t eps)
 {
-	init_shogun();
+	return CMath::abs(a-b)<=eps;
+}
+
+void test_ev()
+{
+#ifdef HAVE_LAPACK
+
+	SGMatrix<float64_t> A(3,3);
+	A(0,0)=0;
+	A(0,1)=1;
+	A(0,2)=0;
+	A(1,0)=1;
+	A(1,1)=0;
+	A(1,2)=1;
+	A(1,0)=0;
+	A(2,1)=1;
+	A(2,2)=0;
+
+	SGVector<float64_t> ev=CMath::compute_eigenvectors(A);
+	CMath::display_matrix(A.matrix, A.num_rows, A.num_cols, "A");
+	CMath::display_vector(ev.vector, ev.vlen, "eigenvalues");
+
+	float64_t sqrt22=CMath::sqrt(2.0)/2.0;
+	float64_t eps=10E-16;
+
+	/* check for correct eigenvectors */
+	ASSERT(is_equal(A(0,0), 0.5, eps));
+	ASSERT(is_equal(A(0,1), -sqrt22, eps));
+	ASSERT(is_equal(A(0,2), 0.5, eps));
+
+	ASSERT(is_equal(A(1,0), -sqrt22, eps));
+	ASSERT(is_equal(A(1,1), 0, eps));
+	ASSERT(is_equal(A(1,2), sqrt22, eps));
+
+	ASSERT(is_equal(A(2,0), 0.5, eps));
+	ASSERT(is_equal(A(2,1), sqrt22, eps));
+	ASSERT(is_equal(A(2,2), 0.5, eps));
 
+	/* check for correct eigenvalues */
+	ASSERT(is_equal(ev[0], -sqrt22*2, eps));
+	ASSERT(is_equal(ev[1], 0, eps));
+	ASSERT(is_equal(ev[2], sqrt22*2, eps));
+
+#endif /* HAVE_LAPACK */
+}
+
+void test_lapack()
+{
 #ifdef HAVE_LAPACK
 	// size of square matrix
 	int N = 100;
@@ -41,13 +89,13 @@ int main(int argc, char** argv)
 			double_matrix[i*N+j] = ((double)(i-j))/(i+j+1);
 
 		double_matrix[i*N+i] += 100;
-	}	
+	}
 	status = 0;
 	wrap_dsygvx(1,'V','U',N,double_matrix,N,double_matrix,N,1,3,double_eigenvalues,double_eigenvectors,&status);
 	if (status!=0)
 	{
 		printf("DSYGVX/SSYGVX failed with code %d\n",status);
-		return -1;
+		ASSERT(false);
 	}
 	delete[] double_eigenvectors;
 
@@ -59,7 +107,7 @@ int main(int argc, char** argv)
 	if (status!=0)
 	{
 		printf("DGEQRF/DORGQR failed with code %d\n",status);
-		return -1;
+		ASSERT(false);
 	}
 	delete[] double_tau;
 
@@ -74,7 +122,7 @@ int main(int argc, char** argv)
 	if (status!=0)
 	{
 		printf("DGESVD failed with code %d\n",status);
-		return -1;
+		ASSERT(false);
 	}
 	delete[] double_s;
 	delete[] double_U;
@@ -86,13 +134,22 @@ int main(int argc, char** argv)
 	if (status!=0)
 	{
 		printf("DSYEV failed with code %d\n",status);
-		return -1;
+		ASSERT(false);
 	}
 	delete[] double_eigenvalues;
 	delete[] double_matrix;
 
 	#endif // HAVE_LAPACK
+}
+
+int main(int argc, char** argv)
+{
+	init_shogun();
+
+	test_lapack();
+	test_ev();
 
 	exit_shogun();
 	return 0;
 }
+