Skip to content

Commit

Permalink
fixed all syntax errors. a linker error remains
Browse files Browse the repository at this point in the history
  • Loading branch information
@jewettaij
jewettaij committed Jan 29, 2020
1 parent 25eabd5 commit 678cb68
Show file tree
Hide file tree
Showing 3 changed files with 75 additions and 77 deletions.
83 changes: 43 additions & 40 deletions include/jacobi.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,9 +6,12 @@

#include <cmath>
#include <cassert>
#include "matrix_alloc.hpp"

namespace jacobi_public_domain {

using namespace matrix_alloc;

/// @class Jacobi
/// @brief Calculate the eigenvalues and eigevectors of a symmetric matrix
/// using the Jacobi eigenvalue algorithm.
Expand All @@ -22,13 +25,13 @@ template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>

class Jacobi
{
int n; //!< the size of the matrix
int *max_idx_row; //!< for row i, the index j of the maximum element where j>i
Scalar **M; //!< local copy of the matrix being analyzed
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:
Scalar c; //!< = cos(θ)
Scalar s; //!< = sin(θ)
Scalar t; //!< = tan(θ), (note |t|<=1)
Scalar c; //!< = cos(θ)
Scalar s; //!< = sin(θ)
Scalar t; //!< = tan(θ), (note |t|<=1)
int *max_idx_row; //!< for row i, the index j of the maximum element where j>i

public:

Expand Down Expand Up @@ -112,9 +115,9 @@ class Jacobi
void Dealloc();

// memory management: copy constructor, swap, and assignment operator
Jacobi(const Jacobi<Scalar, Vector, Matrix>& source);
void swap(Jacobi<Scalar, Vector, Matrix> &other);
Jacobi<Scalar, Vector, Matrix>& operator = (Jacobi<Scalar, Vector, Matrix> source);
Jacobi(const Jacobi<Scalar, Vector, Matrix, ConstMatrix>& source);
void swap(Jacobi<Scalar, Vector, Matrix, ConstMatrix> &other);
Jacobi<Scalar, Vector, Matrix, ConstMatrix>& operator = (Jacobi<Scalar, Vector, Matrix, ConstMatrix> source);

}; // class Jacobi

Expand All @@ -126,7 +129,7 @@ class Jacobi



template<typename Scalar, typename Vector, typename Matrix>
template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>
int Jacobi<Scalar, Vector, Matrix, ConstMatrix>::
Diagonalize(ConstMatrix mat, // the matrix you wish to diagonalize (size n)
Vector eval, // store the eigenvalues here
Expand Down Expand Up @@ -191,8 +194,8 @@ Diagonalize(ConstMatrix mat, // the matrix you wish to diagonalize (size n)
/// M[i][j] = 0. The results will be stored in c, s, and t
/// (which store cos(θ), sin(θ), and tan(θ), respectively).

template<typename Scalar, typename Vector, typename Matrix>
void Jacobi<Scalar, Vector, Matrix>::
template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>
void Jacobi<Scalar, Vector, Matrix, ConstMatrix>::
CalcRot(Matrix M, // matrix
int i, // row index
int j) // column index
Expand Down Expand Up @@ -281,8 +284,8 @@ CalcRot(Matrix M, // matrix
/// |_ . _|
/// @endcode

template<typename Scalar, typename Vector, typename Matrix>
void Jacobi<Scalar, Vector, Matrix>::
template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>
void Jacobi<Scalar, Vector, Matrix, ConstMatrix>::
ApplyRot(Matrix M, // matrix
int i, // row index
int j) // column index
Expand Down Expand Up @@ -357,8 +360,8 @@ ApplyRot(Matrix M, // matrix
/// E'_uv = Σ_w R_wu * E_wv
/// @endcode

template<typename Scalar, typename Vector, typename Matrix>
void Jacobi<Scalar, Vector, Matrix>::
template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>
void Jacobi<Scalar, Vector, Matrix, ConstMatrix>::
ApplyRotLeft(Matrix E, // matrix
int i, // row index
int j) // column index
Expand All @@ -373,8 +376,8 @@ ApplyRotLeft(Matrix E, // matrix



template<typename Scalar, typename Vector, typename Matrix>
int Jacobi<Scalar, Vector, Matrix>::
template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>
int Jacobi<Scalar, Vector, Matrix, ConstMatrix>::
MaxEntryRow(Matrix M, int i) const {
int j_max = i+1;
for(int j = i+2; j < n; j++)
Expand All @@ -385,8 +388,8 @@ MaxEntryRow(Matrix M, int i) const {



template<typename Scalar, typename Vector, typename Matrix>
void Jacobi<Scalar, Vector, Matrix>::
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 {
// find the maximum entry in the matrix M in O(n) time
i_max = 0; // (start with an arbitrary
Expand Down Expand Up @@ -414,8 +417,8 @@ MaxEntry(Matrix M, int& i_max, int& j_max) const {


//Sort the rows of a matrix "evec" by the numbers contained in "eval"
template<typename Scalar, typename Vector, typename Matrix>
void Jacobi<Scalar, Vector, Matrix>::
template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>
void Jacobi<Scalar, Vector, Matrix, ConstMatrix>::
SortRows(Vector eval, Matrix evec, int n, SortCriteria sort_criteria) const
{
for (int i = 0; i < n-1; i++) {
Expand Down Expand Up @@ -452,40 +455,40 @@ SortRows(Vector eval, Matrix evec, int n, SortCriteria sort_criteria) const



template<typename Scalar, typename Vector, typename Matrix>
void Jacobi<Scalar, Vector, Matrix>::
template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>
void Jacobi<Scalar, Vector, Matrix, ConstMatrix>::
Init() {
n = 0;
max_idx_row = nullptr;
}

template<typename Scalar, typename Vector, typename Matrix>
void Jacobi<Scalar, Vector, Matrix>::
template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>
void Jacobi<Scalar, Vector, Matrix, ConstMatrix>::
SetSize(int n) {
Dealloc();
Alloc(n);
}

template<typename Scalar, typename Vector, typename Matrix>
void Jacobi<Scalar, Vector, Matrix>::
template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>
void Jacobi<Scalar, Vector, Matrix, ConstMatrix>::
Alloc(int n) {
this->n = n;
max_idx_row = new int[n];
Alloc2D(n, n, &M);
}

template<typename Scalar, typename Vector, typename Matrix>
void Jacobi<Scalar, Vector, Matrix>::
template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>
void Jacobi<Scalar, Vector, Matrix, ConstMatrix>::
Dealloc() {
if (max_idx_row)
delete [] max_idx_row;
Dealloc2D(&M);
Init();
}

template<typename Scalar, typename Vector, typename Matrix>
Jacobi<Scalar, Vector, Matrix>::
Jacobi(const Jacobi<Scalar, Vector, Matrix>& source)
template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>
Jacobi<Scalar, Vector, Matrix, ConstMatrix>::
Jacobi(const Jacobi<Scalar, Vector, Matrix, ConstMatrix>& source)
{
Init();
Alloc(source.n);
Expand All @@ -499,19 +502,19 @@ Jacobi(const Jacobi<Scalar, Vector, Matrix>& source)
M[i]);
}

template<typename Scalar, typename Vector, typename Matrix>
void Jacobi<Scalar, Vector, Matrix>::
swap(Jacobi<Scalar, Vector, Matrix> &other) {
template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>
void Jacobi<Scalar, Vector, Matrix, ConstMatrix>::
swap(Jacobi<Scalar, Vector, Matrix, ConstMatrix> &other) {
std::swap(max_idx_row, other.max_idx_row);
for (int i = 0; i < n; i++)
std::swap(M[i], other.M[i]);
std::swap(n, other.n);
}

template<typename Scalar, typename Vector, typename Matrix>
Jacobi<Scalar, Vector, Matrix>&
Jacobi<Scalar, Vector, Matrix>::
operator = (Jacobi<Scalar, Vector, Matrix> source) {
template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>
Jacobi<Scalar, Vector, Matrix, ConstMatrix>&
Jacobi<Scalar, Vector, Matrix, ConstMatrix>::
operator = (Jacobi<Scalar, Vector, Matrix, ConstMatrix> source) {
this->swap(source);
return *this;
}
Expand Down
40 changes: 17 additions & 23 deletions include/matrix_alloc.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,17 +9,19 @@

namespace matrix_alloc {

/// @brief Allocate a 2-dimensional table row-major order
template<typename Entry, typename Integer>
void Alloc2D(Integer const size[2], //!< size of the array in x,y directions
Entry ***paaX); //!< pointer to 2-D multidimensional array

/// @brief
/// Slightly different version of Alloc2D()
/// Allocate a 2-dimensional table row-major order.
/// In this version, the the size of the array specified by 2 integer arguments.
template<typename Entry>
void Alloc2D(size_t M, //!< size of the array (number of rows)
size_t N, //!< size of the array (number of columns)
void Alloc2D(int M, //!< size of the array (number of rows)
int N, //!< size of the array (number of columns)
Entry ***paaX); //!< pointer to 2-D multidimensional array


/// @brief Allocate a 2-dimensional table. After allocation, the entries in
/// the table (*paaX)[i][j] exist for all 0<=i<size[0] and 0<=j<size[1].
template<typename Entry, typename Integer>
void Alloc2D(Integer const size[2], //!< size of the array in x,y directions
Entry ***paaX); //!< pointer to 2-D multidimensional array

/// @brief
Expand All @@ -29,27 +31,19 @@ template<typename Entry>
void Dealloc2D(Entry ***paaX); //!< pointer to 2-D multidimensional array



// -------------- IMPLEMENTATION --------------
// --- IMPLEMENTATION ---

template<typename Entry>
void Alloc2D(size_t nrows, //!< size of the array (number of rows)
size_t ncolumns, //!< size of the array (number of columns)
void Alloc2D(int nrows, //!< size of the array (outer)
int ncolumns, //!< size of the array (inner)
Entry ***paaX) //!< pointer to 2-D multidimensional array
{
assert(paaX);

Entry *aX = new Entry [size[0] * size[1]];

// Allocate a conventional 2-dimensional
// pointer-to-a-pointer data structure.
*paaX = new Entry* [size[1]];
for(Integer iy=0; iy<size[1]; iy++)
(*paaX)[iy] = &(aX[ iy*size[0] ]);
// The caller can access the contents of *paX using (*paaX)[i][j] notation.
int size[2];
size[0] = nrows;
size[1] = ncolumns;
Alloc2D<Entry, int>(size, paaX);
}


template<typename Entry>
void Dealloc2D(Entry ***paaX) //!< pointer to 2-D multidimensional array
{
Expand Down
29 changes: 15 additions & 14 deletions tests/test.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,7 @@ inline static bool SimilarVecUnsigned(Vector a, Vector b, int n, Scalar eps=1.0e

/// @brief Multiply two matrices A and B, store the result in C. (C = AB).

template<typename Scalar, typename Matrix, typename ConstMatrix>
template<typename Matrix, typename ConstMatrix>
void mmult(ConstMatrix A, //<! input array
ConstMatrix B, //<! input array
Matrix C, //<! store result here
Expand Down Expand Up @@ -174,13 +174,14 @@ void GenRandOrth(Matrix R,
template <typename Scalar, typename Vector, typename Matrix>
void GenRandSymm(Matrix M, //<! store the matrix here
int n, //<! matrix size
Scalar eval_magnitude_range=2.0,//<! range of wanted eigevalues
int n_degeneracy=1,//<!number of repeated eigevalues(1disables)
std::default_random_engine &generator,//<! makes random numbers
Vector evals, //<! store the eigenvalues of here
Matrix evects //<! store the eigenvectors here
Matrix evects, //<! store the eigenvectors here
std::default_random_engine &generator,//<! makes random numbers
Scalar eval_magnitude_range=2.0,//<! range of wanted eigevalues
int n_degeneracy=1//<!number of repeated eigevalues(1disables)
)
{
std::uniform_real_distribution<Scalar> random_real01;
Scalar **D, **tmp;
Alloc2D(n, n, &D);
Alloc2D(n, n, &tmp);
Expand Down Expand Up @@ -224,6 +225,7 @@ void TestJacobi(int n, //<! matrix size
int n_matrices=100, //<! number of matrices to test
Scalar eval_magnitude_range=2.0, //<! range of eigevalues
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)
)
{
Expand All @@ -239,7 +241,6 @@ void TestJacobi(int n, //<! matrix size
if (seed == 0) // if the caller did not specify a seed, use the system clock
seed = std::chrono::system_clock::now().time_since_epoch().count();
std::default_random_engine generator(seed);
std::uniform_real_distribution<Scalar> random_real01;

Scalar **M, **evects, **evects_known;
Alloc2D(n, n, &M);
Expand All @@ -257,11 +258,11 @@ void TestJacobi(int n, //<! matrix size

GenRandSymm<Scalar, Scalar*, Scalar**>(M,
n,
eval_magnitude_range,
n_degeneracy,
generator,
evals_known,
evects_known);
evects_known,
generator,
eval_magnitude_range,
n_degeneracy);

// Sort the matrix evals and eigenvector rows
SortRows<Scalar, Scalar*, Scalar**> (evals_known, evects_known, n);
Expand Down Expand Up @@ -305,8 +306,8 @@ void TestJacobi(int n, //<! matrix size
for (int a = 0; a < n; a++) {
test_evec[i] = 0.0;
for (int b = 0; b < n; b++)
test_evec[a] += M[a][b] * evec[i][b];
assert(Similar(test_evec[a], eval[i] * evec[i][a], 1.0e-06));
test_evec[a] += M[a][b] * evects[i][b];
assert(Similar(test_evec[a], evals[i] * evects[i][a], 1.0e-06));
}
}

Expand All @@ -315,10 +316,10 @@ void TestJacobi(int n, //<! matrix size
} //for(int imat = 0; imat < n_matrices; imat++) {

Dealloc2D(&M);
Dealloc2D(&evects_known);
Dealloc2D(&evects);
delete [] evals_known;
Dealloc2D(&evects_known);
delete [] evals;
delete [] evals_known;
delete [] test_evec;

}
Expand Down

0 comments on commit 678cb68

Please sign in to comment.