calculating eigen values

include/cif++/matrix.hpp

@@ -28,8 +28,10 @@
 
 #include <array>
 #include <cassert>
+#include <cmath>
 #include <cstdint>
 #include <ostream>
+#include <tuple>
 #include <type_traits>
 #include <vector>
 
@@ -378,7 +380,7 @@ class matrix_matrix_multiplication : public matrix_expression<matrix_matrix_mult
 	const M2 &m_m2;
 };
 
-template<typename M, typename T>
+template <typename M, typename T>
 class matrix_scalar_multiplication : public matrix_expression<matrix_scalar_multiplication<M, T>>
 {
   public:
@@ -403,7 +405,6 @@ class matrix_scalar_multiplication : public matrix_expression<matrix_scalar_mult
 	value_type m_v;
 };
 
-
 template <typename M1, typename T, std::enable_if_t<std::is_floating_point_v<T>, int> = 0>
 auto operator*(const matrix_expression<M1> &m, T v)
 {
@@ -452,6 +453,121 @@ matrix3x3<F> inverse(const matrix3x3<F> &m)
 	return result;
 }
 
+template <typename M>
+auto eigen(const matrix_expression<M> &mat)
+{
+	using value_type = decltype(mat(0, 0));
+
+	assert(mat.dim_m() == mat.dim_n());
+
+	const size_t N = mat.dim_m();
+
+	matrix<float> m = mat;
+
+	matrix<value_type> em = identity_matrix(N);
+	std::vector<value_type> ev(N);
+	std::vector<value_type> b(N);
+
+	// Set ev & b to diagonal.
+	for (size_t i = 0; i < N; ++i)
+		ev[i] = b[i] = m(i, i);
+
+	for (int cyc = 1; cyc <= 50; ++cyc)
+	{
+		// calc sum of diagonal, off-diagonal
+		value_type spp = 0, spq = 0;
+
+		for (size_t i = 0; i < N - 1; i++)
+		{
+			for (size_t j = i + 1; j < N; j++)
+				spq += std::fabs(m(i, j));
+			spp += std::fabs(m(i, i));
+		}
+
+		if (spq <= 1.0e-12 * spp)
+			break;
+
+		std::vector<value_type> z(N);
+
+		// now try and reduce each off-diagonal element in turn
+		for (size_t i = 0; i < N - 1; i++)
+		{
+			for (size_t j = i + 1; j < N; j++)
+			{
+				value_type a_ij = m(i, j);
+				value_type h = ev[j] - ev[i];
+				value_type t;
+
+				if (std::fabs(a_ij) > 1.0e-12 * std::fabs(h))
+				{
+					value_type theta = 0.5 * h / a_ij;
+					t = 1.0 / (std::fabs(theta) + std::sqrt(1 + theta * theta));
+					if (theta < 0)
+						t = -t;
+				}
+				else
+					t = a_ij / h;
+
+				// calc trig properties
+				value_type c = 1.f / std::sqrt(1 + t * t);
+				value_type s = t * c;
+				value_type tau = s / (1 + c);
+				h = t * a_ij;
+
+				// update eigenvalues
+				z[i] -= h;
+				z[j] += h;
+				ev[i] -= h;
+				ev[j] += h;
+
+				// rotate the upper diagonal of the matrix
+				m(i, j) = 0;
+
+				for (size_t k = 0; k < i; k++)
+				{
+					value_type ai = m(k, i);
+					value_type aj = m(k, j);
+					m(k, i) = ai - s * (aj + ai * tau);
+					m(k, j) = aj + s * (ai - aj * tau);
+				}
+
+				for (size_t k = i + 1; k < j; k++)
+				{
+					value_type ai = m(i, k);
+					value_type aj = m(k, j);
+					m(i, k) = ai - s * (aj + ai * tau);
+					m(k, j) = aj + s * (ai - aj * tau);
+				}
+
+				for (size_t k = j + 1; k < N; k++)
+				{
+					value_type ai = m(i, k);
+					value_type aj = m(j, k);
+					m(i, k) = ai - s * (aj + ai * tau);
+					m(j, k) = aj + s * (ai - aj * tau);
+				}
+
+				// apply corresponding rotation to result
+				for (size_t k = 0; k < N; k++)
+				{
+					value_type ai = em(k, i);
+					value_type aj = em(k, j);
+					em(k, i) = ai - s * (aj + ai * tau);
+					em(k, j) = aj + s * (ai - aj * tau);
+				}
+			}
+		}
+
+		for (size_t p = 0; p < N; p++)
+		{
+			b[p] += z[p];
+			ev[p] = b[p];
+		}
+	}
+
+	return std::make_tuple(ev, em);
+}
+
 // --------------------------------------------------------------------
 
 template <typename M>

src/symmetry.cpp

@@ -404,8 +404,6 @@ std::tuple<float,point,sym_op> closest_symmetry_copy(const spacegroup &sg, const
 	sym_op result_s;
 
 	auto o = offsetToOrigin(c, a);
-	transformation tlo(identity_matrix<float>(3), o);
-	auto itlo = cif::inverse(tlo);
 
 	a += o;
 	b += o;

test/unit-3d-test.cpp

@@ -360,3 +360,78 @@ BOOST_AUTO_TEST_CASE(symm_2bi3_1, *utf::tolerance(0.1f))
 		BOOST_TEST(so.string() == symm2);
 	}
 }
+
+// --------------------------------------------------------------------
+
+BOOST_AUTO_TEST_CASE(eigen_1, *utf::tolerance(0.1f))
+{
+	cif::symmetric_matrix4x4<float> m;
+
+	m(0, 0) = 4;
+	m(0, 1) = -30;
+	m(0, 2) = 60;
+	m(0, 3) = -35;
+	m(1, 1) = 300;
+	m(1, 2) = -675;
+	m(1, 3) = 420;
+	m(2, 2) = 1620;
+	m(2, 3) = -1050;
+	m(3, 3) = 700;
+
+	cif::matrix4x4<float> m2;
+	m2 = m;
+
+	const auto &[ev, em] = cif::eigen(m2);
+
+	BOOST_TEST(ev[0] == 0.1666428611718905f);
+	BOOST_TEST(ev[1] == 1.4780548447781369f);
+	BOOST_TEST(ev[2] == 37.1014913651276582f);
+	BOOST_TEST(ev[3] == 2585.25381092892231f);
+
+
+
+
+
+
+// 	=== Example ===
+
+// Let 
+// <math>
+// 	S = \begin{pmatrix} 4 & -30 & 60 & -35 \\ -30 & 300 & -675 & 420 \\ 60 & -675 & 1620 & -1050 \\ -35 & 420 & -1050 & 700 \end{pmatrix}
+// </math>
+
+// Then ''jacobi'' produces the following eigenvalues and eigenvectors after 3 sweeps (19 iterations) :
+
+// <math>
+// 	e_1 = 2585.25381092892231
+// </math>
+
+// <math>
+// 	E_1 = \begin{pmatrix}0.0291933231647860588\\ -0.328712055763188997\\ 0.791411145833126331\\ -0.514552749997152907\end{pmatrix}
+// </math>
+
+// <math>
+// 	e_2 = 37.1014913651276582
+// </math>
+
+// <math>
+// 	E_2 = \begin{pmatrix}-0.179186290535454826\\ 0.741917790628453435\\ -0.100228136947192199\\ -0.638282528193614892\end{pmatrix}
+// </math>
+
+// <math>
+// 	e_3 = 1.4780548447781369
+// </math>
+
+// <math>
+// 	E_3 = \begin{pmatrix}-0.582075699497237650\\ 0.370502185067093058\\ 0.509578634501799626\\ 0.514048272222164294\end{pmatrix}
+// </math>
+
+// <math>
+// 	e_4 = 0.1666428611718905
+// </math>
+
+// <math>
+// 	E_4 = \begin{pmatrix}0.792608291163763585\\ 0.451923120901599794\\ 0.322416398581824992\\ 0.252161169688241933\end{pmatrix}
+// </math>
+
+}
\ No newline at end of file