Updated ARPACK and LAPACK libshogun examples

shogun-toolbox · Oct 1, 2011 · 3f082b3 · 3f082b3
1 parent a945b61
commit 3f082b3
Show file tree

Hide file tree

Showing 2 changed files with 90 additions and 45 deletions.
diff --git a/examples/undocumented/libshogun/mathematics_arpack.cpp b/examples/undocumented/libshogun/mathematics_arpack.cpp
@@ -22,32 +22,59 @@ int main(int argc, char** argv)
 	int N = 100;
 	int nev = 2;
 
-	double* matrix = new double[N*N];
+	double* double_matrix = new double[N*N];
+	float* float_matrix = new float[N*N];
 
-	double* eigenvalues = new double[nev];
-	double* eigenvectors = new double[nev*N];
+	double* rhs_double_diag = new double[N];
+	float* rhs_float_diag = new float[N];
+
+	double* double_eigenvalues = new double[nev];
+	float* float_eigenvalues = new float[nev];
+
+	double* double_eigenvectors = new double[nev*N];
+	float* float_eigenvectors = new float[nev*N];
 
 	for (int i=0; i<N; i++)
 	{
+		rhs_double_diag[i] = 1.0;
+		rhs_float_diag[i] = 1.0;
 		for (int j=0; j<N; j++)
-			matrix[i*N+j] = i*i+j*j;
+		{
+			double_matrix[i*N+j] = i*i+j*j;
+			float_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);
+	arpack_xsxupd<double>(double_matrix, NULL, N, 2, "LM", 1, false, 0.0, 0.0,
+	                      double_eigenvalues, double_eigenvectors, status);
+	arpack_xsxupd<float>(float_matrix, NULL, N, 2, "LM", 1, false, 0.0, 0.0,
+	                     float_eigenvalues, float_eigenvectors, status);
 	if (status!=0)
 		return -1;
 
-	arpack_dsaeupd_wrap(matrix, NULL, N, 2, "BE", 3, false, 1.0, 0.0,
-	                    eigenvalues, eigenvectors, status);
+	arpack_xsxupd<double>(double_matrix, NULL, N, 2, "BE", 3, false, 1.0, 0.0,
+	                      double_eigenvalues, double_eigenvectors, status);
+	arpack_xsxupd<float>(float_matrix, NULL, N, 2, "BE", 3, false, 1.0, 0.0,
+	                     float_eigenvalues, float_eigenvectors, status);
 	if (status!=0)
 		return -1;
 
+	arpack_xsxupd<double>(double_matrix, rhs_double_diag, N, 2, "SM", 3, false, 0.0, 0.0,
+	                      double_eigenvalues, double_eigenvectors, status);
+	arpack_xsxupd<float>(float_matrix, rhs_float_diag, N, 2, "SM", 3, false, 0.0, 0.0,
+	                     float_eigenvalues, float_eigenvectors, status);
+	if (status!=0)
+		return -1;
 
-	delete[] eigenvalues;
-	delete[] eigenvectors;
-	delete[] matrix;
+	delete[] double_eigenvalues;
+	delete[] double_eigenvectors;
+	delete[] double_matrix;
+	delete[] rhs_double_diag;
+	delete[] float_eigenvalues;
+	delete[] float_eigenvectors;
+	delete[] float_matrix;
+	delete[] rhs_float_diag;
 #endif // HAVE_ARPACK
 
 	exit_shogun();

diff --git a/examples/undocumented/libshogun/mathematics_lapack.cpp b/examples/undocumented/libshogun/mathematics_lapack.cpp
@@ -23,80 +23,98 @@ int main(int argc, char** argv)
 	int N = 100;
 
 	// square matrix
-	double* matrix = new double[N*N];
+	double* double_matrix = new double[N*N];
+	float* float_matrix = new float[N*N];
 	// for storing eigenpairs
-	double* eigenvalues = new double[N];
-	double* eigenvectors = new double[N*N];
+	double* double_eigenvalues = new double[N];
+	float* float_eigenvalues = new float[N];
+	double* double_eigenvectors = new double[N*N];
+	float* float_eigenvectors = new float[N*N];
 	// for SVD
-	double* U = new double[N*N];
-	double* s = new double[N*N];
-	double* Vt = new double[N*N];
+	double* double_U = new double[N*N];
+	double* double_s = new double[N];
+	double* double_Vt = new double[N*N];
+	float* float_U = new float[N*N];
+	float* float_s = new float[N];
+	float* float_Vt = new float[N*N];
 	// status (should be zero)
 	int status;
 
 	// DSYGVX
 	for (int i=0; i<N; i++)
 	{
 		for (int j=0; j<N; j++)
-			matrix[i*N+j] = (i-j)/(i+j+1);
-		matrix[i*N+i] += 10;
+		{
+			double_matrix[i*N+j] = ((double)(i-j))/(i+j+1);
+			float_matrix[i*N+j] = ((float)(i-j))/(i+j+1);
+		}
+		double_matrix[i*N+i] += 100;
+		float_matrix[i*N+i] += 100;
 	}	
 	status = 0;
-	wrap_dsygvx(1,'V','U',N,matrix,N,matrix,N,1,3,eigenvalues,eigenvectors,&status);
+	wrap_dsygvx(1,'V','U',N,double_matrix,N,double_matrix,N,1,3,double_eigenvalues,double_eigenvectors,&status);
+	wrap_ssygvx(1,'V','U',N,float_matrix,N,float_matrix,N,1,3,float_eigenvalues,float_eigenvectors,&status);
 	if (status!=0)
 	{
-		printf("DSYGVX failed with code %d\n",status);
+		printf("DSYGVX/SSYGVX failed with code %d\n",status);
 		return -1;
 	}
-	delete[] eigenvectors;
-
+	delete[] double_eigenvectors;
+	delete[] float_eigenvectors;
 
 	// DGEQRF+DORGQR
 	status = 0;
-	double* tau = new double[N];
-	wrap_dgeqrf(N,N,matrix,N,tau,&status);
-	wrap_dorgqr(N,N,N,matrix,N,tau,&status);
+	double* double_tau = new double[N];
+	float* float_tau = new float[N];
+	wrap_dgeqrf(N,N,double_matrix,N,double_tau,&status);
+	wrap_dorgqr(N,N,N,double_matrix,N,double_tau,&status);
+	wrap_sgeqrf(N,N,float_matrix,N,float_tau,&status);
+	wrap_sorgqr(N,N,N,float_matrix,N,float_tau,&status);
 	if (status!=0)
 	{
-		printf("DGEQRF/DORGQR failed with code %d\n",status);
+		printf("DGEQRF/DORGQR/SGEQRF/SORGQR failed with code %d\n",status);
 		return -1;
 	}
-	delete[] tau;
-
+	delete[] double_tau;
+	delete[] float_tau;
 
 	// DGESVD
 	for (int i=0; i<N; i++)
 	{
 		for (int j=0; j<N; j++)
-			matrix[i*N+j] = i*i+j*j;
+		{
+			double_matrix[i*N+j] = i*i+j*j;
+			float_matrix[i*N+j] = i*i+j*j;
+		}
 	}
 	status = 0;
-	wrap_dgesvd('A','A',N,N,matrix,N,s,U,N,Vt,N,&status);
+	wrap_dgesvd('A','A',N,N,double_matrix,N,double_s,double_U,N,double_Vt,N,&status);
+	wrap_sgesvd('A','A',N,N,float_matrix,N,float_s,float_U,N,float_Vt,N,&status);
 	if (status!=0)
 	{
-		printf("DGESVD failed with code %d\n",status);
+		printf("DGESVD/SGESVD failed with code %d\n",status);
 		return -1;
 	}
-	delete[] s;
-	delete[] U;
-	delete[] Vt;
-
+	delete[] double_s;
+	delete[] double_U;
+	delete[] double_Vt;
+	delete[] float_s;
+	delete[] float_U;
+	delete[] float_Vt;
 
 	// DSYEV
-	for (int i=0; i<N; i++)
-	{
-		for (int j=0; j<N; j++)
-			matrix[i*N+j] = i*i+j*j;
-	}
 	status = 0;
-	wrap_dsyev('V','U',N,matrix,N,eigenvalues,&status);
+	wrap_dsyev('V','U',N,double_matrix,N,double_eigenvalues,&status);
+	wrap_ssyev('V','U',N,float_matrix,N,float_eigenvalues,&status);
 	if (status!=0)
 	{
-		printf("DSYEV failed with code %d\n",status);
+		printf("DSYEV/SSYEV failed with code %d\n",status);
 		return -1;
 	}
-	delete[] eigenvalues;
-	delete[] matrix;
+	delete[] double_eigenvalues;
+	delete[] float_eigenvalues;
+	delete[] double_matrix;
+	delete[] float_matrix;
 
 	#endif // HAVE_LAPACK