From a5a9352372a2063193126cd371359725eeff80fb Mon Sep 17 00:00:00 2001 From: Andrew Jewett Date: Fri, 27 Nov 2020 21:38:44 -0800 Subject: [PATCH] 1) Changed "matrix_alloc" to "matrix_alloc_jpd" to avoid possible future 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 (although the old assert statements remain commented out). --- .travis.yml | 4 ++++ include/jacobi_pd.hpp | 24 +++++++++++-------- ...{matrix_alloc.hpp => matrix_alloc_jpd.hpp} | 17 ++++++------- tests/test_jacobi.cpp | 5 ++-- 4 files changed, 30 insertions(+), 20 deletions(-) rename include/{matrix_alloc.hpp => matrix_alloc_jpd.hpp} (77%) diff --git a/.travis.yml b/.travis.yml index 917a5f5..d64a6a2 100644 --- a/.travis.yml +++ b/.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), diff --git a/include/jacobi_pd.hpp b/include/jacobi_pd.hpp index a88d544..2561c08 100644 --- a/include/jacobi_pd.hpp +++ b/include/jacobi_pd.hpp @@ -9,12 +9,12 @@ #include #include -#include -#include "matrix_alloc.hpp" +//#include +#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& 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). diff --git a/include/matrix_alloc.hpp b/include/matrix_alloc_jpd.hpp similarity index 77% rename from include/matrix_alloc.hpp rename to include/matrix_alloc_jpd.hpp index 62326bb..2e4c742 100644 --- a/include/matrix_alloc.hpp +++ b/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 -#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 @@ -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 diff --git a/tests/test_jacobi.cpp b/tests/test_jacobi.cpp index 004ffc8..3788fd3 100644 --- a/tests/test_jacobi.cpp +++ b/tests/test_jacobi.cpp @@ -3,6 +3,7 @@ #include #include +#include #include #include #include @@ -10,7 +11,7 @@ #include #include #include -#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;