jacobi_pd finally works with vector<vector<X>> and array<array<X,n>,n…

…>. (The tests I wrote in the previous commit were failing to compile the way I thought they were, causing me to erroneously believe that this feature was already working.) I have not come up with a working example which uses fixed-size C arrays yet. Working on that next.

.travis.yml

@@ -32,7 +32,12 @@ matrix:
         #- now use valgrind to find memory leaks and other errors:
         - g++ -g -O0 -I../include -o test_jacobi test.cpp
         - valgrind --leak-check=yes --error-exitcode=1 ./test_jacobi 5 10 1e-09 1e+09
-        - valgrind --leak-check=yes --error-exitcode=1 ./test_jacobi 5 10 1e-09 1e+09 2
+        #- now test again using vector<vector<double>> instead of double **
+        - g++ -g -O0 -DUSE_VECTOR_OF_VECTORS -I../include -o test_jacobi test.cpp
+        - valgrind --leak-check=yes --error-exitcode=1 ./test_jacobi 5 10 1e-09 1e+09
+        #- now test again using array<array<double,5>,5>
+        - g++ -g -O0 -DUSE_ARRAY_OF_ARRAY -I../include -o test_jacobi test.cpp
+        - valgrind --leak-check=yes --error-exitcode=1 ./test_jacobi 5 10 1e-09 1e+09
         - cd ../
 
     - language: cpp
@@ -61,6 +66,11 @@ matrix:
         - ./test_jacobi 20 100 1e-09 1e+09 5
         - clang++ -g -O0 -I../include -o test_jacobi test.cpp
         - valgrind --leak-check=yes --error-exitcode=1 ./test_jacobi 5 10 1e-09 1e+09
-        - valgrind --leak-check=yes --error-exitcode=1 ./test_jacobi 5 10 1e-09 1e+09 2
+        #- now test again using vector<vector<double>> instead of double **
+        - clang++ -g -O0 -DUSE_VECTOR_OF_VECTORS -I../include -o test_jacobi test.cpp
+        - valgrind --leak-check=yes --error-exitcode=1 ./test_jacobi 5 10 1e-09 1e+09
+        #- now test again using array<array<double,5>,5>
+        - clang++ -g -O0 -DUSE_ARRAY_OF_ARRAY -I../include -o test_jacobi test.cpp
+        - valgrind --leak-check=yes --error-exitcode=1 ./test_jacobi 5 10 1e-09 1e+09
         - cd ../
 

README.md

@@ -16,7 +16,7 @@ This remains one of the oldest and most popular algorithms for
 diagonalizing dense, square, real, symmetric matrices.
 
 The matrices themselves can be implemented as \*\*X (pointer-to-pointer),
-vector\<vector\<X\>\>, array\<array\<X\>\>, fixed-size arrays,
+vector\<vector\<X\>\>&, fixed-size arrays,
 or any other C or C++ object which supports double-indexing.
 (Here **X** is any real numeric type.  Complex numbers are not supported.)
 
@@ -63,10 +63,14 @@ double **evects; // Store the eigenvectors here.
 // multiple times without incurring the cost of allocating memory on the heap.
 Jacobi<double, double*, double**, double const*const*> eigen_calc(n);
 
+// Note: If you prefer using vectors over C-style pointers, this works also:
+// Jacobi<double, vector<double>&, vector<vector<double>>&,
+//        vector<vector<double>>&>  eigen_calc(n);
+
 // Now, calculate the eigenvalues and eigenvectors of M
 eigen_calc.Diagonalize(M, evals, evects);
 ```
-*(A complete working example can be found [here](tests/test.cpp).)*
+*(A working example can be found [here](tests/test.cpp).)*
 
 ## Installation
 

include/jacobi.hpp

@@ -76,15 +76,15 @@ class Jacobi
   ///        both sides) will zero the ij'th element of M, so that afterwards
   ///        M[i][j] = 0.  The results will be stored in c, s, and t
   ///        (which store cos(θ), sin(θ), and tan(θ), respectively).
-  void CalcRot(Matrix M,   //!< matrix
+  void CalcRot(Scalar const *const *M,   //!< matrix
                int i,      //!< row index
                int j);     //!< column index
 
   /// @brief Apply the (previously calculated) rotation matrix to matrix M
   ///        by multiplying it on both sides (a "similarity transform").
   ///        (To save time, only update the elements in the upper-right
   ///         triangular region of the matrix.  It is assumed that i < j.)
-  void ApplyRot(Matrix M,  //!< matrix
+  void ApplyRot(Scalar **M,  //!< matrix
                 int i,     //!< row index
                 int j);    //!< column index
 
@@ -95,15 +95,15 @@ class Jacobi
                     int j);    //!< column index
 
   ///@brief Find the off-diagonal index in row i whose absolute value is largest
-  int MaxEntryRow(Matrix M, int i) const;
+  int MaxEntryRow(Scalar const *const *M, int i) const;
 
   /// @brief Find the indices (i_max, j_max) marking the location of the
   ///        entry in the matrix with the largest absolute value.  This
   ///        uses the max_idx_row[] array to find the answer in O(n) time.
-  /// @returns This function does not return avalue.  However after it is
+  /// @returns This function does not return a avalue.  However after it is
   ///          invoked, the location of the largest matrix element will be
   ///          stored in the i_max and j_max arguments.
-  void MaxEntry(Matrix M, int& i_max, int& j_max) const;
+  void MaxEntry(Scalar const *const *M, int& i_max, int& j_max) const;
 
   // @brief Sort the rows in M according to the numbers in v (also sorted)
   void SortRows(Vector v, //!< vector containing the keys used for sorting
@@ -151,7 +151,7 @@ Diagonalize(ConstMatrix mat,    // the matrix you wish to diagonalize (size n)
 
   if (calc_evec)
     for (int i = 0; i < n; i++)
-      for (int j = 0; j < n; j++)        //Initialize the evec[][] matrix.
+      for (int j = 0; j < n; j++)
         evec[i][j] = (i==j) ? 1.0 : 0.0; //Set evec equal to the identity matrix
 
   for (int i = 0; i < n-1; i++)          //Initialize the "max_idx_row[]" array 
@@ -202,7 +202,7 @@ Diagonalize(ConstMatrix mat,    // the matrix you wish to diagonalize (size n)
 
 template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>
 void Jacobi<Scalar, Vector, Matrix, ConstMatrix>::
-CalcRot(Matrix M,    // matrix
+CalcRot(Scalar const *const *M,    // matrix
         int i,       // row index
         int j)       // column index
 {
@@ -292,7 +292,7 @@ CalcRot(Matrix M,    // matrix
 
 template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>
 void Jacobi<Scalar, Vector, Matrix, ConstMatrix>::
-ApplyRot(Matrix M,  // matrix
+ApplyRot(Scalar **M,  // matrix
          int i,     // row index
          int j)     // column index
 {
@@ -384,7 +384,7 @@ ApplyRotLeft(Matrix E,  // matrix
 
 template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>
 int Jacobi<Scalar, Vector, Matrix, ConstMatrix>::
-MaxEntryRow(Matrix M, int i) const {
+MaxEntryRow(Scalar const *const *M, int i) const {
   int j_max = i+1;
   for(int j = i+2; j < n; j++)
     if (std::abs(M[i][j]) > std::abs(M[i][j_max]))
@@ -396,7 +396,7 @@ MaxEntryRow(Matrix M, int i) const {
 
 template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>
 void Jacobi<Scalar, Vector, Matrix, ConstMatrix>::
-MaxEntry(Matrix M, int& i_max, int& j_max) const {
+MaxEntry(Scalar const *const *M, int& i_max, int& j_max) const {
   // find the maximum entry in the matrix M in O(n) time
   i_max = 0; // (start with an arbitrary
   j_max = 1; //  off-diagonal element: M[0][1])
@@ -531,7 +531,7 @@ Jacobi(Jacobi<Scalar, Vector, Matrix, ConstMatrix>&& other) {
   swap(*this, other);
 }
 
-// Using the "copy-swap" idiom for the assignment operator (=)
+// Using the "copy-swap" idiom for the assignment operator
 template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>
 Jacobi<Scalar, Vector, Matrix, ConstMatrix>&
 Jacobi<Scalar, Vector, Matrix, ConstMatrix>::

tests/test.cpp

@@ -5,16 +5,33 @@
 #include <chrono>
 #include <random>
 #include <algorithm>
+#include <vector>
+#include <array>
 #include "matrix_alloc.hpp"
 #include "jacobi.hpp"
 
 using std::cout;
 using std::cerr;
 using std::endl;
 using std::setprecision;
+using std::vector;
+using std::array;
 using namespace matrix_alloc;
 using namespace jacobi_public_domain;
 
+
+// This code works with various types of C++ matrices (for example,
+// double **, vector<vector<double>> array<array<double,5>,5>).
+// I use "#if defined" statements to test different matrix types.
+// For some of these (eg. array<array<double,5>,5>), the size of the matrix
+// must be known at compile time.  I specify that size now.
+#if defined USE_ARRAY_OF_ARRAYS
+const int NF=5;  //(the array size must be known at compile time)
+#elif defined USE_C_FIXED_SIZE_ARRAYS
+const int NF=5;  //(the array size must be known at compile time)
+#endif
+
+
 // @brief  Are two numbers "similar"?
 template<typename Scalar>
 inline static bool Similar(Scalar a, Scalar b,
@@ -63,7 +80,6 @@ void mmult(ConstMatrix A, //<! input array
            int K=0     //<! optional: number of columns of A = num rows of B (=m by default)
            )
 {
-  assert((C != A) && (C != B));
   if (n == 0) n = m; // if not specified, then assume the matrices are square
   if (K == 0) K = m; // if not specified, then assume the matrices are square
 
@@ -193,9 +209,20 @@ void GenRandSymm(Matrix M,       //<! store the matrix here
 {
   assert(n_degeneracy <= n);
   std::uniform_real_distribution<Scalar> random_real01;
+  #if defined USE_VECTOR_OF_VECTORS
+  vector<vector<Scalar> > D(n, vector<Scalar>(n));
+  vector<vector<Scalar> > tmp(n, vector<Scalar>(n));
+  #elif defined USE_ARRAY_OF_ARRAYS
+  array<array<Scalar, NF>, NF> D;
+  array<array<Scalar, NF>, NF> tmp;
+  #elif defined USE_C_FIXED_SIZE_ARRAYS
+  Scalar D[NF][NF], tmp[NF][NF];
+  #else
+  #define USE_C_POINTER_TO_POINTERS
   Scalar  **D, **tmp;
   Alloc2D(n, n, &D);
   Alloc2D(n, n, &tmp);
+  #endif
 
   // Randomly generate the eigenvalues
   for (int i = 0; i < n; i++) {
@@ -231,13 +258,46 @@ void GenRandSymm(Matrix M,       //<! store the matrix here
   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
+
+  // Original code:
+  //mmult(evects, D, tmp, n);  // <--> tmp = Rt * D
+  // Unfortunately, C++ guesses the types incorrectly.  Must manually specify:
+  // #ifdefs making the code ugly again:
+  #if defined USE_VECTOR_OF_VECTORS
+  mmult<vector<vector<Scalar> >&, vector<vector<Scalar> >&>
+  #elif defined USE_ARRAY_OF_ARRAYS
+  mmult<array<array<Scalar,NF>,NF>&, array<array<Scalar,NF>,NF>&>
+  #elif defined USE_C_FIXED_SIZE_ARRAYS
+  mmult<Scalar *[NF], Scalar *[NF]>
+  #else
+  mmult<Scalar**, Scalar const *const *>
+  #endif
+       (evects, D, tmp, n);
+
   for (int i = 0; i < n-1; i++)
     for (int j = i+1; j < n; j++)
       std::swap(evects[i][j], evects[j][i]); //transpose "evects"
-  mmult(tmp, evects, M, n);  //at this point M = Rt * D * R (where "R"="evects")
+
+  // Original code:
+  //mmult(tmp, evects, M, n);
+  // Unfortunately, C++ guesses the types incorrectly.  Must manually specify:
+  // #ifdefs making the code ugly again:
+  #if defined USE_VECTOR_OF_VECTORS
+  mmult<vector<vector<Scalar> >&, vector<vector<Scalar> >&>
+  #elif defined USE_ARRAY_OF_ARRAYS
+  mmult<array<array<Scalar,NF>,NF>&, array<array<Scalar,NF>,NF>&>
+  #elif defined USE_C_FIXED_SIZE_ARRAYS
+  mmult<Scalar *[NF], Scalar *[NF]>
+  #else
+  mmult<Scalar**, Scalar const *const *>
+  #endif
+       (tmp, evects, M, n);
+  //at this point M = Rt*D*R (where "R"="evects")
+
+  #if defined USE_C_POINTER_TO_POINTERS
   Dealloc2D(&D);
   Dealloc2D(&tmp);
+  #endif
 } // GenRandSymm()
 
 
@@ -248,6 +308,7 @@ void TestJacobi(int n, //<! matrix size
                 Scalar min_eval_size=0.1,  //<! minimum possible eigenvalue sizw
                 Scalar max_eval_size=10.0, //<! maximum possible eigenvalue size
                 int n_tests_per_matrix=1, //<! repeat test for benchmarking?
+
                 int n_degeneracy=1, //<! repeated eigenvalues?
                 unsigned seed=0 //<! random seed (if 0 then use the clock)
                 )
@@ -262,12 +323,20 @@ void TestJacobi(int n, //<! matrix size
 
 
 
-  // --------- The following code 
+  // Create an instance of the Jacobi diagonalizer, and allocate the matrix
+  // we will test it on, as well as the arrays that will store the resulting
+  // eigenvalues and eigenvectors.
+  // The way we do this dependa on what version of the code we are using.
+  // This is controlled by "#if defined" statements.
 
-  #if defined USE_VECTORS_OF_VECTORS
+  #if defined USE_VECTOR_OF_VECTORS
+
+  Jacobi<Scalar,
+         vector<Scalar>&,
+         vector<vector<Scalar> >&,
+         vector<vector<Scalar> >& >
+    ecalc(n);
 
-  Jacobi<Scalar, vector<Scalar>, vector<vector<Scalar> >,
-         const vector<const vector<Scalar> > > ecalc(n);
   // allocate the matrix, eigenvalues, eigenvectors
   vector<vector<Scalar> > M(n, vector<Scalar>(n));
   vector<vector<Scalar> > evects(n, vector<Scalar>(n));
@@ -276,34 +345,37 @@ void TestJacobi(int n, //<! matrix size
   vector<Scalar> evals_known(n);
   vector<Scalar> test_evec(n);
 
-  #elseif defined USE_ARRAYS_OF_ARRAYS
+  #elif defined USE_ARRAY_OF_ARRAYS
 
-  n = 5;
+  n = NF;
   cout << "Testing std::array (fixed size).\n"
     "(Ignoring first argument, and setting matrix size to " << n << ")" << endl;
-  Jacobi<Scalar, array<Scalar, 5>, array<array<Scalar, 5>, 5>,
-         const array<const array<Scalar, 5>, 5> > ecalc(n);
+  Jacobi<Scalar,
+         array<Scalar, NF>&,
+         array<array<Scalar, NF>, NF>&,
+         array<array<Scalar, NF>, NF>& >
+    ecalc(n);
   // allocate the matrix, eigenvalues, eigenvectors
-  array<array<Scalar, 5>, 5> M;
-  array<array<Scalar, 5>, 5> evects;
-  array<array<Scalar, 5>, 5> evects_known;
-  array<Scalar, 5> evals;
-  array<Scalar, 5> evals_known;
-  array<Scalar, 5> test_evec;
+  array<array<Scalar, NF>, NF> M;
+  array<array<Scalar, NF>, NF> evects;
+  array<array<Scalar, NF>, NF> evects_known;
+  array<Scalar, NF> evals;
+  array<Scalar, NF> evals_known;
+  array<Scalar, NF> test_evec;
 
-  #elseif defined USE_C_FIXED_SIZE_ARRAYS
+  #elif defined USE_C_FIXED_SIZE_ARRAYS
 
-  n = 5;
+  n = NF;
   cout << "Testing C fixed size arrays.\n"
     "(Ignoring first argument, and setting matrix size to " << n << ")" << endl;
-  Jacobi<Scalar, Scalar*, Scalar *[5], Scalar *[5]> ecalc(n);
+  Jacobi<Scalar, Scalar*, Scalar *[NF], Scalar *[NF]> ecalc(n);
   // allocate the matrix, eigenvalues, eigenvectors
-  Scalar[5][5] M;
-  Scalar[5][5] evects;
-  Scalar[5][5] evects_known;
-  Scalar[5] evals;
-  Scalar[5] evals_known;
-  Scalar[5] test_evec;
+  Scalar M[NF][NF];
+  Scalar evects[NF][NF];
+  Scalar evects_known[NF][NF];
+  Scalar evals[NF];
+  Scalar evals_known[NF];
+  Scalar test_evec[NF];
 
   #else
 
@@ -333,15 +405,26 @@ void TestJacobi(int n, //<! matrix size
   #endif
 
 
-  // Now, generate random matrices and test Jacobi::Diagonalize() on them
+  // --------------------------------------------------------------------
+  // Now, generate random matrices and test Jacobi::Diagonalize() on them.
+  // --------------------------------------------------------------------
 
   for(int imat = 0; imat < n_matrices; imat++) {
 
     // Create a randomly generated symmetric matrix.
     //This function generates random numbers for the eigenvalues ("evals_known")
     //as well as the eigenvectors ("evects_known"), and uses them to generate M.
 
-    GenRandSymm(M,
+    #if defined USE_VECTOR_OF_VECTORS
+    GenRandSymm<Scalar, vector<Scalar>&, vector<vector<Scalar> >&>
+    #elif defined USE_ARRAY_OF_ARRAYS
+    GenRandSymm<Scalar, array<Scalar,NF>&, array<array<Scalar,NF>,NF>&>
+    #elif defined USE_C_FIXED_SIZE_ARRAYS
+    GenRandSymm<Scalar, Scalar*, Scalar *[NF]>
+    #else
+    GenRandSymm<Scalar, Scalar*, Scalar**>
+    #endif
+               (M,
                 n,
                 evals_known,
                 evects_known,
@@ -350,8 +433,21 @@ void TestJacobi(int n, //<! matrix size
                 max_eval_size,
                 n_degeneracy);
 
-    // Sort the matrix evals and eigenvector rows
-    SortRows<Scalar>(evals_known, evects_known, n);
+    // Sort the matrix evals and eigenvector rows:
+    // Original code:
+    //SortRows<Scalar>(evals_known, evects_known, n);
+    // Unfortunately, C++ guesses the types incorrectly. Must use #ifdefs again:
+    #if defined USE_VECTOR_OF_VECTORS
+    SortRows<Scalar, vector<Scalar>&, vector<vector<Scalar> >&>
+    #elif defined USE_ARRAY_OF_ARRAYS
+    SortRows<Scalar, array<Scalar,NF>&, array<array<Scalar,NF>,NF>&>
+    #elif defined USE_C_FIXED_SIZE_ARRAYS
+    SortRows<Scalar, Scalar*, Scalar *[NF]>
+    #else
+    SortRows<Scalar, Scalar*, Scalar**>
+    #endif
+            (evals_known, evects_known, n);
+
 
     if (n_matrices == 1) {
       cout << "Eigenvalues (after sorting):\n";