From 7aa7c061ceb47ddb385a702a076d343ba63e4e96 Mon Sep 17 00:00:00 2001
From: toshinagata1964 <toshinagata1964@a2be9bc6-48de-4e38-9406-05402d4bc13c>
Date: Sun, 12 Aug 2012 13:18:31 +0000
Subject: [PATCH] LAMatrix singular value decomposition was implemented.

git-svn-id: svn+ssh://svn.sourceforge.jp/svnroot/molby/trunk@277 a2be9bc6-48de-4e38-9406-05402d4bc13c
---
 MolLib/Ruby_bind/ruby_types.c | 29 ++++++++++++++++++++++++++++-
 MolLib/Types.c                | 35 +++++++++++++++++++++++++++++++++++
 MolLib/Types.h                |  1 +
 3 files changed, 64 insertions(+), 1 deletion(-)

diff --git a/MolLib/Ruby_bind/ruby_types.c b/MolLib/Ruby_bind/ruby_types.c
index fe267b0..7ffff61 100644
--- a/MolLib/Ruby_bind/ruby_types.c
+++ b/MolLib/Ruby_bind/ruby_types.c
@@ -2234,7 +2234,7 @@ s_LAMatrix_Eigenvalues(VALUE self)
 	int info;
 	Data_Get_Struct(self, LAMatrix, mp1);
 	if (mp1->column != mp1->row)
-		rb_raise(rb_eArgError, "cannot get eigenvectors for non-symmetric matrix");
+		rb_raise(rb_eArgError, "cannot get eigenvectors for non-square matrix");
 	mp2 = LAMatrixNew(mp1->row, 1);
 	mp3 = LAMatrixNew(mp1->row, mp1->column);
 	if ((info = LAMatrixSymDiagonalize(mp2, mp3, mp1)) != 0)
@@ -2244,6 +2244,32 @@ s_LAMatrix_Eigenvalues(VALUE self)
 
 /*
  *  call-seq:
+ *     svd -> [left_matrix, singular_values, right_matrix]
+ *
+ *  Decompose the given (m,n) matrix to a product of three matrices, U, S, V.
+ *  U is a (m,m) orthogonal matrix, S is a (m,n) matrix which is zero except for
+ *  min(m,n) diagonal elements, and V is a (n,n) orthogonal matrix.
+ *  (Usually SVD is defined as M = U*S*transpose(V), but this methods returns
+ *  transpose(V) rather than V.) The singular_values is a min(m,n) dimension
+ *  column vector.
+ */
+static VALUE
+s_LAMatrix_SVD(VALUE self)
+{
+	LAMatrix *mp1, *mp2, *mp3, *mp4;
+	int info, n;
+	Data_Get_Struct(self, LAMatrix, mp1);
+	mp2 = LAMatrixNew(mp1->row, mp1->row);
+	n = (mp1->column > mp1->row ? mp1->row : mp1->column);
+	mp3 = LAMatrixNew(n, 1);
+	mp4 = LAMatrixNew(mp1->column, mp1->column);
+	if ((info = LAMatrixSingularValueDecomposition(mp2, mp3, mp4, mp1)) != 0)
+		rb_raise(rb_eArgError, "cannot perform singular value decomposition");
+	return rb_ary_new3(3, Data_Wrap_Struct(rb_cLAMatrix, 0, -1, mp2), Data_Wrap_Struct(rb_cLAMatrix, 0, -1, mp3), Data_Wrap_Struct(rb_cLAMatrix, 0, -1, mp4));
+}
+
+/*
+ *  call-seq:
  *     LAMatrix[f1, f2, ..., fn] -> (new) LAMatrix
  *     LAMatrix[a1, a2, ..., an] -> (new) LAMatrix
  *     LAMatrix[obj] -> (new) LAMatrix
@@ -2901,6 +2927,7 @@ Init_MolbyTypes(void)
 	rb_define_method(rb_cLAMatrix, "column_size", s_LAMatrix_ColumnSize, 0);
 	rb_define_method(rb_cLAMatrix, "row_size", s_LAMatrix_RowSize, 0);
 	rb_define_method(rb_cLAMatrix, "eigenvalues", s_LAMatrix_Eigenvalues, 0);
+	rb_define_method(rb_cLAMatrix, "svd", s_LAMatrix_SVD, 0);
 	rb_define_method(rb_cLAMatrix, "to_a", s_LAMatrix_ToArray, 0);
 	rb_define_method(rb_cLAMatrix, "inspect", s_LAMatrix_Inspect, 0);
 	rb_define_alias(rb_cLAMatrix, "to_s", "inspect");
diff --git a/MolLib/Types.c b/MolLib/Types.c
index e7b041f..e0e95f8 100755
--- a/MolLib/Types.c
+++ b/MolLib/Types.c
@@ -576,6 +576,41 @@ LAMatrixSymDiagonalize(LAMatrix *eigenValues, LAMatrix *eigenVectors, const LAMa
 	return info;
 }
 
+/*  Singular value decomposition  */
+/*  M = U W Vt, M: (m, n) Matrix, U: (m, m) orthogonal matrix, W: min(m, n) column vector, Vt: (n, n) orthogonal matrix (not V)  */
+int
+LAMatrixSingularValueDecomposition(LAMatrix *matU, LAMatrix *matW, LAMatrix *matV, const LAMatrix *mat)
+{
+	__CLPK_integer m, n, lda, num, *iwork, lwork, info;
+	__CLPK_doublereal workSize, *work, *matData;
+	m = mat->row;
+	n = mat->column;
+	lda = m;
+	num = (m < n ? m : n);
+	LAMatrixResize(matW, num, 1);
+	LAMatrixResize(matU, m, m);
+	LAMatrixResize(matV, n, n);
+	matData = (__CLPK_doublereal *)malloc(sizeof(__CLPK_doublereal) * n * m);
+	memmove(matData, mat->data, sizeof(__CLPK_doublereal) * n * m);
+	iwork = (__CLPK_integer *)malloc(sizeof(__CLPK_integer) * num);
+	lwork = -1;
+	info = 0;
+	dgesdd_("A", &m, &n, matData, &lda, matW->data, matU->data, &m, matV->data, &n, &workSize, &lwork, iwork, &info);
+	
+	if (info != 0) {
+		free(matData);
+		free(iwork);
+		return info;
+	}
+	lwork = workSize;
+	work = (__CLPK_doublereal *)malloc(sizeof(__CLPK_doublereal) * lwork);
+	dgesdd_("A", &m, &n, matData, &lda, matW->data, matU->data, &m, matV->data, &n, work, &lwork, iwork, &info);
+	free(work);
+	free(iwork);
+	free(matData);
+	return info;
+}
+
 #pragma mark ==== Array ====
 
 /*  Assign a value to an array. An array is represented by two fields; count and base,
diff --git a/MolLib/Types.h b/MolLib/Types.h
index 3c36863..0462e83 100755
--- a/MolLib/Types.h
+++ b/MolLib/Types.h
@@ -145,6 +145,7 @@ int LAMatrixInvert(LAMatrix *mat1, const LAMatrix *mat2);
 Double LAMatrixDeterminant(const LAMatrix *mat);
 void LAMatrixTranspose(LAMatrix *mat1, const LAMatrix *mat2);
 int LAMatrixSymDiagonalize(LAMatrix *vec, LAMatrix *mat1, const LAMatrix *mat2);
+int LAMatrixSingularValueDecomposition(LAMatrix *matU, LAMatrix *matW, LAMatrix *matV, const LAMatrix *mat);
 
 /*  Utility functions  */
 void SetPanicFunc(void (*func)(const char *, ...));
-- 
2.11.0