{

int n; //!< the size of the matrix

Scalar **M; //!< local copy of the matrix being analyzed

// Precomputed cosine, sine, and tangent of the most recent rotation angle:

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:

/// @brief Specify the size of the matrices you want to diagonalize later.

/// @param n the size (ie. number of rows) of the square matrix

void SetSize(int n);

Jacobi(int n = 0) {

Init();

SetSize(n);

}

~Jacobi() {

Dealloc();

}

// @typedef choose the criteria for sorting eigenvalues and eigenvectors

typedef enum eSortCriteria {

DO_NOT_SORT,

SORT_DECREASING_EVALS,

SORT_INCREASING_EVALS,

SORT_DECREASING_ABS_EVALS,

SORT_INCREASING_ABS_EVALS

} SortCriteria;

/// @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 Jacobi iterations attempted, which should be > 0.

/// If the return value is not strictly > 0 then 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)

Vector eval, //!< store the eigenvalues here

Matrix evec, //!< store the eigenvectors here (in rows)

SortCriteria sort_criteria=SORT_DECREASING_EVALS,//!<sort results?

bool calc_evecs=true, //!< calculate the eigenvectors?

int max_num_sweeps = 50); //!< limit the number of iterations

private:

/// @brief Calculate the components of a rotation matrix which performs a

/// rotation in the i,j plane by an angle (θ) causing M[i][j]=0.

/// The resulting parameters will be stored in this->c, this->s, and

/// this->t (which store cos(θ), sin(θ), and tan(θ), respectively).

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 the elements in the upper-right-triangular

/// region of the matrix are updated. It is assumed that i < j.)

void ApplyRot(Scalar **M, //!< matrix

int i, //!< row index

int j); //!< column index

/// @brief Multiply matrix E on the left by the (previously calculated)

/// rotation matrix.

void ApplyRotLeft(Matrix E, //!< matrix

int i, //!< row index

int j); //!< column index

///@brief Find the off-diagonal index in row i whose absolute value is largest

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 a value. 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(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

Matrix M, //!< matrix whose rows will be sorted according to v

int n, //!< size of the vector and matrix

SortCriteria s=SORT_DECREASING_EVALS //!< sort decreasing order?

) const;

// memory management:

void Alloc(int N);

void Init();

void Dealloc();

public:

// memory management: copy and move constructor, swap, and assignment operator

Jacobi(const Jacobi<Scalar, Vector, Matrix, ConstMatrix>& source);

Jacobi(Jacobi<Scalar, Vector, Matrix, ConstMatrix>&& other);

void swap(Jacobi<Scalar, Vector, Matrix, ConstMatrix> &other) noexcept;

Jacobi<Scalar, Vector, Matrix, ConstMatrix>& operator = (Jacobi<Scalar, Vector, Matrix, ConstMatrix> source);