added unit tests for degenerate eigenvalues

.travis.yml

@@ -20,7 +20,14 @@ matrix:
         - ./test_jacobi 4 1000 3.0
         - ./test_jacobi 5 1000 3.0
         - ./test_jacobi 10 100 3.0
-        - ./test_jacobi 50 20 3.0
+        - ./test_jacobi 30 20 3.0
+        #- now test with degenerate eigenvalues
+        - ./test_jacobi 2 1000 3.0 2
+        - ./test_jacobi 3 1000 3.0 2
+        - ./test_jacobi 4 1000 3.0 2
+        - ./test_jacobi 5 1000 3.0 3
+        - ./test_jacobi 10 100 3.0 5
+        - ./test_jacobi 30 20 3.0 20
         - cd ../
 
     - language: cpp
@@ -38,6 +45,13 @@ matrix:
         - ./test_jacobi 4 1000 3.0
         - ./test_jacobi 5 1000 3.0
         - ./test_jacobi 10 100 3.0
-        - ./test_jacobi 50 20 3.0
+        - ./test_jacobi 30 20 3.0
+        #- now test with degenerate eigenvalues
+        - ./test_jacobi 2 1000 3.0 2
+        - ./test_jacobi 3 1000 3.0 2
+        - ./test_jacobi 4 1000 3.0 2
+        - ./test_jacobi 5 1000 3.0 3
+        - ./test_jacobi 10 100 3.0 5
+        - ./test_jacobi 30 20 3.0 20
         - cd ../
 

README.md

@@ -78,13 +78,11 @@ This is a header-only library.
 
 # Development Status: *alpha*
 
-As of 2020-1-23, basic functionality appears to be working.
+As of 2020-1-29, basic functionality appears to be working.
 More testing is needed including tests for copy constructors,
 memory leaks and profiling.
 The API might change slightly, but existing code built using
 it should still work.
-(Later I might add "const" to the first argument of Jacobi::Diagonalize().)
-
 
 ## License
 

include/jacobi.hpp

@@ -27,7 +27,7 @@ class Jacobi
 {
   int n;            //!< the size of the matrix
   Scalar **M;       //!< local copy of the matrix being analyzed
-  // Precomputed cosine, sin, and tangent of the most recent rotation angle:
+  // Precomputed cosine, sine, and tangent of the most recent rotation angle:
   Scalar c;         //!< = cos(θ)
   Scalar s;         //!< = sin(θ)
   Scalar t;         //!< = tan(θ),  (note |t|<=1)
@@ -222,12 +222,12 @@ CalcRot(Matrix M,    // matrix
 }
 
 
-/// brief   Perform a similarity transform by multiplying matrix M on both
-///         sides by a rotation matrix (transposing one of them).
+/// brief   Perform a similarity transformation by multiplying matrix M on both
+///         sides by a rotation matrix (and its transpose).
 ///         This rotation matrix performs a rotation in the i,j plane by
-///         angle θ.  Also updates the max_idx_row[] array.  Assumes i < j.
-/// details This function assumes that i<j and that cos(θ), sin(θ), and tan(θ)
-///         have already been computed (by invoking CalcRot()).
+///         angle θ.  This function assumes that c=cos(θ). s=som(θ), t=tan(θ)
+///         have been calculated in advance (using the CalcRot() function).
+///         It also assumes that i<j.  The max_idx_row[] array is also updated.
 ///         To save time, since the matrix is symmetric, the elements
 ///         below the diagonal (ie. M[u][v] where u>v) are not computed.
 /// @code

include/matrix_alloc.hpp

@@ -65,8 +65,6 @@ void Alloc2D(int nrows,             //!< size of the array (outer)
 template<typename Entry>
 void Dealloc2D(Entry ***paaX)        //!< pointer to 2-D multidimensional array
 {
-  //*paaX = new Entry* [size[0]];
-  //(*paaX)[0] = new Entry [size[0] * size[1]];
   if (paaX && *paaX) {
     delete [] (*paaX)[0];
     delete [] (*paaX);

tests/test.cpp

@@ -4,6 +4,7 @@
 #include <cstdlib>
 #include <chrono>
 #include <random>
+#include <algorithm>
 #include "matrix_alloc.hpp"
 #include "jacobi.hpp"
 
@@ -17,7 +18,9 @@ using namespace jacobi_public_domain;
 // @brief  Are two numbers "similar"?
 template<typename T>
 inline static bool Similar(T a, T b, T eps=1.0e-06) {
-  return std::abs(a - b) <= std::abs(eps)*std::sqrt(std::abs(a)*std::abs(b));
+  return ((std::abs(a-b) <= std::abs(eps))
+          ||
+          (std::abs(a-b) <= std::abs(eps)*std::sqrt(std::abs(a)*std::abs(b))));
 }
 
 /// @brief  Are two vectors (containing n numbers) similar?
@@ -186,11 +189,6 @@ void GenRandSymm(Matrix M, //<! store the matrix here
   Alloc2D(n, n, &D);
   Alloc2D(n, n, &tmp);
 
-  // D is a diagonal matrix whose diagonal elements are the eigenvalues
-  for (int i = 0; i < n; i++)
-    for (int j = 0; j < n; j++)
-      D[i][j] = 0.0;
-
   // Randomly generate the eigenvalues
   for (int i = 0; i < n; i++) {
     // generate some random eigenvalues
@@ -202,11 +200,26 @@ void GenRandSymm(Matrix M, //<! store the matrix here
     // also consider both positive and negative eigenvalues:
     if (random_real01(generator) < 0.5)
       evals[i] = -evals[i];
-    D[i][i] = evals[i];
   }
 
-  // Now randomly generate the R and Rt matrices (ie. the eigenvectors):
-  GenRandOrth<Scalar, Matrix>(evects, n, generator);
+  // Does the user want us to force some of the eigenvalues to be the same?
+  if (n_degeneracy > 1) {
+    unsigned *permutation = new unsigned[n]; //a random permutation from 0...n-1
+    for (int i = 0; i < n; i++)
+      permutation[i] = i;
+    std::shuffle(permutation, permutation+n, generator);
+    for (int i = 1; i < n_degeneracy; i++) //set the first n_degeneracy to same
+      evals[permutation[i]] = evals[permutation[0]];
+    delete [] permutation;
+  }
+
+  // D is a diagonal matrix whose diagonal elements are the eigenvalues
+  for (int i = 0; i < n; i++)
+    for (int j = 0; j < n; j++)
+      D[i][j] = ((i == j) ? evals[i] : 0.0);
+
+  // Now randomly generate the (transpose of) the "evects" matrix
+  GenRandOrth<Scalar, Matrix>(evects, n, generator); //(will transpose it later)
 
   // Construct the test matrix, M, where M = Rt * D * R
   mmult(evects, D, tmp, n);  //tmp = Rt * D
@@ -301,7 +314,8 @@ void TestJacobi(int n, //<! matrix size
       }
 
       assert(SimilarVec(evals, evals_known, n, 1.0e-06));
-      //Check each eigenvector
+      //Check that each eigenvector behaves like an eigenvector
+      //   Σ_b  M[a][b]*evects[i][b] = evals[i]*evects[i][b]   (for all a)
       for (int i = 0; i < n; i++) {
         for (int a = 0; a < n; a++) {
           test_evec[a] = 0.0;
@@ -322,7 +336,7 @@ void TestJacobi(int n, //<! matrix size
   delete [] evals_known;
   delete [] test_evec;
 
-}
+} //TestJacobi()
 
 
 int main(int argc, char **argv) {
@@ -363,10 +377,10 @@ int main(int argc, char **argv) {
     erange = std::stof(argv[3]);
   if (argc > 4)
     n_tests = std::stoi(argv[4]);
+  if (argc > 5)
+    n_degeneracy = std::stoi(argv[5]);
   if (argc > 6)
-    n_degeneracy = std::stoi(argv[6]);
-  if (argc > 7)
-    seed = std::stoi(argv[7]);
+    seed = std::stoi(argv[6]);
 
   TestJacobi(n_size, n_matr, erange, n_tests, n_degeneracy, seed);