1) Changed "matrix_alloc" to "matrix_alloc_jpd" to avoid possible fut…

…ure naming collisions. 2) Jacobi::Diagonalize() returns number of iters (>0) not "sweeps", or 0 when convergence failed, 3) Now runs on 1x1 matrices without crashing 4) removed dependence on <cassert> (although the old assert statements remain commented out).

.travis.yml

@@ -16,6 +16,8 @@ matrix:
         - cd tests
         - g++ -g -O0 -coverage -fno-inline -fno-inline-small-functions -fno-default-inline -I../include -o test_jacobi test_jacobi.cpp
         - ./test_jacobi 5 1000 1e-03 1e+03 1 0  # ("0" tests for code-coverage)
+        # Edge case: Check if it runs without crashing on 1x1 matrices:
+        - ./test_jacobi 1 1 0.1 10.0
         #- cppcheck --quiet --error-exitcode=1
         # Now turn off code-coverage testing and try different size matrices.
         # Note: With such a huge spread of eigenvalue magnitudes(1e-09...1e+09),
@@ -63,6 +65,8 @@ matrix:
         - cd tests
         - clang++ -g -O0 -coverage -fno-inline -I../include -o test_jacobi test_jacobi.cpp
         - ./test_jacobi 5 1000 1e-03 1e+03 1 0  # ("0" tests for code-coverage)
+        # Edge case: Check if it runs without crashing on 1x1 matrices:
+        - ./test_jacobi 1 1 0.1 10.0
         #- cppcheck --quiet --error-exitcode=1
         # Now turn off code-coverage testing and try different size matrices.
         # Note: With such a huge spread of eigenvalue magnitudes(1e-09...1e+09),

include/jacobi_pd.hpp

@@ -9,12 +9,12 @@
 
 #include <algorithm>
 #include <cmath>
-#include <cassert>
-#include "matrix_alloc.hpp"
+//#include <cassert>
+#include "matrix_alloc_jpd.hpp"
 
 namespace jacobi_pd {
 
-using namespace matrix_alloc;
+using namespace matrix_alloc_jpd;
 
 /// @class Jacobi
 /// @brief Calculate the eigenvalues and eigevectors of a symmetric matrix
@@ -63,8 +63,8 @@ class Jacobi
   /// @brief Calculate all the eigenvalues and eigevectors of a symmetric matrix
   ///        using the Jacobi eigenvalue algorithm:
   ///        https://en.wikipedia.org/wiki/Jacobi_eigenvalue_algorithm
-  /// @returns The number_of_sweeps (= number_of_iterations / (n*(n-1)/2)).
-  ///          If this equals max_num_sweeps, the algorithm failed to converge.
+  /// @returns The number of Jacobi rotations attempted if successful (>0),
+  ///          and 0 otherwise.  (0 indicates that convergence failed.)
   /// @note  To reduce the computation time further, set calc_evecs=false.
   int
   Diagonalize(ConstMatrix mat, //!< the matrix you wish to diagonalize (size n)
@@ -163,7 +163,7 @@ Diagonalize(ConstMatrix mat,    // the matrix you wish to diagonalize (size n)
   // -- Iteration --
   int n_iters;
   int max_num_iters = max_num_sweeps*n*(n-1)/2; //"sweep" = n*(n-1)/2 iters
-  for (n_iters=0; n_iters < max_num_iters; n_iters++) {
+  for (n_iters=1; n_iters <= max_num_iters; n_iters++) {
     int i,j;
     MaxEntry(M, i, j); // Find the maximum entry in the matrix. Store in i,j
 
@@ -191,7 +191,11 @@ Diagonalize(ConstMatrix mat,    // the matrix you wish to diagonalize (size n)
   // Optional: Sort results by eigenvalue.
   SortRows(eval, evec, n, sort_criteria);
 
-  return n_iters / (n*(n-1)/2); //returns the number of "sweeps" (converged?)
+  if ((n_iters > max_num_iters) && (n>1))   // If we exceeded max_num_iters,
+    return 0;                               // indicate an error occured.
+
+  //assert(n_iters > 0);
+  return n_iters;
 }
 
 
@@ -223,7 +227,7 @@ CalcRot(Scalar const *const *M,    // matrix
         t = -t;
     }
   }
-  assert(std::abs(t) <= 1.0);
+  //assert(std::abs(t) <= 1.0);
   c = 1.0 / std::sqrt(1 + t*t);
   s = c*t;
 }
@@ -313,7 +317,7 @@ ApplyRot(Scalar **M,  // matrix
 
   //Update the off-diagonal elements of M which will change (above the diagonal)
 
-  assert(i < j);
+  //assert(i < j);
 
   M[i][j] = 0.0;
 
@@ -507,7 +511,7 @@ Jacobi(const Jacobi<Scalar, Vector, Matrix, ConstMatrix>& source)
 {
   Init();
   SetSize(source.n);
-  assert(n == source.n);
+  //assert(n == source.n);
   // The following lines aren't really necessary, because the contents
   // of source.M and source.max_idx_row are not needed (since they are
   // overwritten every time Jacobi::Diagonalize() is invoked).

include/matrix_alloc.hpp → include/matrix_alloc_jpd.hpp

@@ -1,15 +1,16 @@
-/// @file    matrix_alloc.hpp
-/// @brief   Because I allocate 2-dimensional arrays frequently, I created a 
-///          few optional functions that make this more convenient.
+/// @file    matrix_alloc_jpd.hpp
+/// @brief   Short functions for allocating 2-dimensional C-style ** arrays.
+///          Perhaps it is useful to put this short code in a separate file 
+///          to make it independently of "jacobi_pd.hpp".
 /// @author  Andrew Jewett
 /// @license CC0-1.0
 
 #include<cassert>
 
-#ifndef _MATRIX_ALLOC_H
-#define _MATRIX_ALLOC_H
+#ifndef _MATRIX_ALLOC_JPD_H
+#define _MATRIX_ALLOC_JPD_H
 
-namespace matrix_alloc {
+namespace matrix_alloc_jpd {
 
 /// @brief  Allocate a 2-dimensional array.  (Uses row-major order.)
 template<typename Entry>
@@ -51,6 +52,6 @@ void Dealloc2D(Entry ***paaX)       // pointer to a 2D C-style array
   }
 }
 
-} // namespace matrix_alloc
+} // namespace matrix_alloc_jpd
 
-#endif //#ifndef _MATRIX_ALLOC_H
+#endif //#ifndef _MATRIX_ALLOC_JPD_H

tests/test_jacobi.cpp

@@ -3,14 +3,15 @@
 
 #include <iostream>
 #include <cmath>
+#include <cassert>
 #include <iomanip>
 #include <cstdlib>
 #include <chrono>
 #include <random>
 #include <algorithm>
 #include <vector>
 #include <array>
-#include "matrix_alloc.hpp"
+#include "matrix_alloc_jpd.hpp"
 #include "jacobi_pd.hpp"
 
 using std::cout;
@@ -19,7 +20,7 @@ using std::endl;
 using std::setprecision;
 using std::vector;
 using std::array;
-using namespace matrix_alloc;
+using namespace matrix_alloc_jpd;
 using namespace jacobi_pd;