start work on eigen decomposition

3 files changed, 146 insertions, 1 deletions

diff --git a/src/msvmmaj_kernel.c b/src/msvmmaj_kernel.c
index 9b421ac..3a92ff5 100644
--- a/src/msvmmaj_kernel.c
+++ b/src/msvmmaj_kernel.c
@@ -62,6 +62,12 @@ void msvmmaj_make_kernel(struct MajModel *model, struct MajData *data)
 		}
 	}
 
+	print_matrix(K, n, n);
+
+	double *P = Malloc(double, n*n);
+	double *Lambda = Malloc(double, n);
+	msvmmaj_make_eigen(K, n, P, Lambda);
+
 	// copy kernel to data
 	data->Z = realloc(data->Z, n*(n+1)*(sizeof(double)));
 	for (i=0; i<n; i++) {
@@ -98,6 +104,106 @@ void msvmmaj_make_kernel(struct MajModel *model, struct MajData *data)
 }
 
 /**
+ * @brief Find the eigendecomposition of a kernel matrix.
+ *
+ * @details.
+ * tbd
+ *
+ *
+ */
+void msvmmaj_make_eigen(double *K, long n, double *P, double *Lambda)
+{
+	int M, status, LWORK, *IWORK, *IFAIL;
+	double abstol, *WORK;
+
+	double *tempP = Malloc(double, n*n);
+
+	IWORK = Malloc(int, 5*n);
+	IFAIL = Malloc(int, n);
+	
+	// highest precision eigenvalues, may reduce for speed	
+	abstol = 2.0*dlamch('S');
+
+	// first perform a workspace query to determine optimal size of the 
+	// WORK array.
+	WORK = Malloc(double, 1);
+	status = dsyevx(
+			'V',
+			'A',
+			'U',
+			n,
+			K,
+			n,
+			0,
+			0,
+			0,
+			0,
+			abstol,
+			&M,
+			Lambda,
+			tempP,
+			n,
+			WORK,
+			-1,
+			IWORK,
+			IFAIL);
+	LWORK = WORK[0];
+
+	// allocate the requested memory for the eigendecomposition 
+	WORK = (double *)realloc(WORK, LWORK*sizeof(double));
+	status = dsyevx(
+			'V',
+			'A',
+			'U',
+			n,
+			K,
+			n,
+			0,
+			0,
+			0,
+			0,
+			abstol,
+			&M,
+			Lambda,
+			tempP,
+			n,
+			WORK,
+			LWORK,
+			IWORK,
+			IFAIL);
+
+	printf("status = %i\n", status);
+	printf("Number of eigenvalues found: %i\n", M);
+
+	if (status != 0) {
+		fprintf(stderr, "Nonzero exit status from dsyevx. Exiting...");
+		exit(1);
+	}
+	// Here you can put code to select the necessary eigenvectors, 
+	// depending on the size of the eigenvalues.
+	// For now, let's just print the eigenvalues and exit
+	long i;
+	for (i=0; i<n; i++) {
+		printf("%+16.16f\n", Lambda[i]);
+	}
+	printf("\n");
+
+	long j;
+	for (i=0; i<n; i++)
+		for (j=0; j<n; j++)
+			P[i*n+j] = tempP[j*n+i];
+
+	print_matrix(P, n, n);
+
+	printf("\n");
+	print_matrix(K, n, n);
+
+	free(tempP);	
+
+	exit(1);
+}
+
+/**
  * @brief Compute the RBF kernel between two vectors
  * 
  * @details
diff --git a/src/msvmmaj_lapack.c b/src/msvmmaj_lapack.c
index e8b5d6b..34bd132 100644
--- a/src/msvmmaj_lapack.c
+++ b/src/msvmmaj_lapack.c
@@ -94,3 +94,42 @@ int dsysv(char UPLO, int N, int NRHS, double *A, int LDA, int *IPIV,
 	dsysv_(&UPLO, &N, &NRHS, A, &LDA, IPIV, B, &LDB, WORK, &LWORK, &INFO);
 	return INFO;
 }
+
+/**
+ * @brief Compute the eigenvalues and optionally the eigenvectors of a 
+ * symmetric matrix.
+ *
+ * @details
+ * See the LAPACK documentation at:
+ * http://www.netlib.org/lapack/explore-html/d2/d97/dsyevx_8f.html
+ *
+ *
+ */
+int dsyevx(char JOBZ, char RANGE, char UPLO, int N, double *A, int LDA,
+	       	double VL, double VU, int IL, int IU, double ABSTOL, int *M,
+		double *W, double *Z, int LDZ, double *WORK, int LWORK,
+		int *IWORK, int *IFAIL)
+{
+	extern void dsyevx_(char *JOBZ, char *RANGE, char *UPLO, int *Np,
+			double *A, int *LDAp, double *VLp, double *VUp,
+		       	int *ILp, int *IUp, double *ABSTOLp, int *M,
+			double *W, double *Z, int *LDZp, double *WORK,
+		       	int *LWORKp, int *IWORK, int *IFAIL, int *INFOp);
+	int INFO;
+	dsyevx_(&JOBZ, &RANGE, &UPLO, &N, A, &LDA, &VL, &VU, &IL, &IU, &ABSTOL,
+			M, W, Z, &LDZ, WORK, &LWORK, IWORK, IFAIL, &INFO);
+	return INFO;
+}
+
+/**
+ * @brief Determine double precision machine parameters.
+ *
+ * @details
+ * See the LAPACK documentation at:
+ * http://www.netlib.org/lapack/explore-html/d5/dd4/dlamch_8f.html
+ */
+double dlamch(char CMACH)
+{
+	extern double dlamch_(char *CMACH);
+	return dlamch_(&CMACH);
+}
diff --git a/src/msvmmaj_matrix.c b/src/msvmmaj_matrix.c
index ffa0c21..3f5bf4a 100644
--- a/src/msvmmaj_matrix.c
+++ b/src/msvmmaj_matrix.c
@@ -146,7 +146,7 @@ void print_matrix(double *M, long rows, long cols)
 
 	for (i=0; i<rows; i++) {
 		for (j=0; j<cols; j++)
-			note("%8.8f ", matrix_get(M, cols, i, j));
+			note("%+6.6f ", matrix_get(M, cols, i, j));
 		note("\n");
 	}
 	note("\n");