diff --git a/.travis.yml b/.travis.yml index 433c6e1..92a8101 100644 --- a/.travis.yml +++ b/.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> 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,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> 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,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 ../ diff --git a/README.md b/README.md index 4424edc..1afbcfa 100644 --- a/README.md +++ b/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\\>, array\\>, fixed-size arrays, +vector\\>&, 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 eigen_calc(n); +// Note: If you prefer using vectors over C-style pointers, this works also: +// Jacobi&, vector>&, +// vector>&> 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 diff --git a/include/jacobi.hpp b/include/jacobi.hpp index 7e9c0b7..175a9dc 100644 --- a/include/jacobi.hpp +++ b/include/jacobi.hpp @@ -76,7 +76,7 @@ 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 @@ -84,7 +84,7 @@ class Jacobi /// 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 void Jacobi:: -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 void Jacobi:: -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 int Jacobi:: -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 void Jacobi:: -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&& other) { swap(*this, other); } -// Using the "copy-swap" idiom for the assignment operator (=) +// Using the "copy-swap" idiom for the assignment operator template Jacobi& Jacobi:: diff --git a/tests/test.cpp b/tests/test.cpp index 1a84dba..29282bf 100644 --- a/tests/test.cpp +++ b/tests/test.cpp @@ -5,6 +5,8 @@ #include #include #include +#include +#include #include "matrix_alloc.hpp" #include "jacobi.hpp" @@ -12,9 +14,24 @@ 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> array,5>). +// I use "#if defined" statements to test different matrix types. +// For some of these (eg. array,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 inline static bool Similar(Scalar a, Scalar b, @@ -63,7 +80,6 @@ void mmult(ConstMatrix A, // random_real01; + #if defined USE_VECTOR_OF_VECTORS + vector > D(n, vector(n)); + vector > tmp(n, vector(n)); + #elif defined USE_ARRAY_OF_ARRAYS + array, NF> D; + array, 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, //(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 >&> + #elif defined USE_ARRAY_OF_ARRAYS + mmult,NF>&, array,NF>&> + #elif defined USE_C_FIXED_SIZE_ARRAYS + mmult + #else + mmult + #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 >&> + #elif defined USE_ARRAY_OF_ARRAYS + mmult,NF>&, array,NF>&> + #elif defined USE_C_FIXED_SIZE_ARRAYS + mmult + #else + mmult + #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, //&, + vector >&, + vector >& > + ecalc(n); - Jacobi, vector >, - const vector > > ecalc(n); // allocate the matrix, eigenvalues, eigenvectors vector > M(n, vector(n)); vector > evects(n, vector(n)); @@ -276,34 +345,37 @@ void TestJacobi(int n, // evals_known(n); vector 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, array, 5>, - const array, 5> > ecalc(n); + Jacobi&, + array, NF>&, + array, NF>& > + ecalc(n); // allocate the matrix, eigenvalues, eigenvectors - array, 5> M; - array, 5> evects; - array, 5> evects_known; - array evals; - array evals_known; - array test_evec; + array, NF> M; + array, NF> evects; + array, NF> evects_known; + array evals; + array evals_known; + array 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 ecalc(n); + Jacobi 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,7 +405,9 @@ void TestJacobi(int n, //&, vector >&> + #elif defined USE_ARRAY_OF_ARRAYS + GenRandSymm&, array,NF>&> + #elif defined USE_C_FIXED_SIZE_ARRAYS + GenRandSymm + #else + GenRandSymm + #endif + (M, n, evals_known, evects_known, @@ -350,8 +433,21 @@ void TestJacobi(int n, //(evals_known, evects_known, n); + // Sort the matrix evals and eigenvector rows: + // Original code: + //SortRows(evals_known, evects_known, n); + // Unfortunately, C++ guesses the types incorrectly. Must use #ifdefs again: + #if defined USE_VECTOR_OF_VECTORS + SortRows&, vector >&> + #elif defined USE_ARRAY_OF_ARRAYS + SortRows&, array,NF>&> + #elif defined USE_C_FIXED_SIZE_ARRAYS + SortRows + #else + SortRows + #endif + (evals_known, evects_known, n); + if (n_matrices == 1) { cout << "Eigenvalues (after sorting):\n";