jewettaij / jacobi_pd

a short, simple public-domain header-only C++ library for calculating eigenvalues and eigenvectors of real symmetric matrices

jacobi_pd

Description

This repository contains a small C++ header file that implements the Jacobi eigenvalue algorithm. It is free of copyright.

The Jacobi algorithm remains one of the oldest and most popular methods for diagonalizing dense, square, real, symmetric matrices.

The matrices passed to to the "Diagonalize()" function can be any C or C++ object which supports [i][j] indexing, including X** (pointer-to-pointer), vector<vector<X>>&, or fixed-size arrays. (Here X is any real numeric type. Complex numbers are not supported.)

(Memory allocation on the heap is avoided except during instantiation.)

The main feature of this repository is it's license.

As of late 2020, no simple public domain C++11 code yet exists for matrix diagonalization. Other C++ libraries such as Eigen or GSL are typically much larger and use more restrictive licenses. (On several occasions, this has prevented me from including their code in other open-source projects with incompatible licenses.) Some repositories may unwittingly contain code snippets from other sources, such as numerical recipes. This short repository was written from scratch. No lines of code were borrowed or adapted from other sources.

Caveats: The code in this repository does not run in parallel, and only works on dense square real symmetric matrices. However it is reasonably short, simple, fast and reliable. You can do anything you like with this code.

Example usage

#include "jacobi_pd.hpp"
using namespace jacobi_pd;

// ...
int n = 3;       // Matrix size
double **M;      // A symmetric n x n matrix you want to diagonalize
double *evals;   // Store the eigenvalues here.
double **evecs;  // Store the eigenvectors here.
// Allocate space for M, evals, and evecs (omitted)...
M[0][0] = 2.0; M[0][1] = 1.0; M[0][2] = 1.0;
M[1][0] = 1.0; M[1][1] = 2.0; M[1][2] =-1.0;  //Note: The matrix
M[2][0] = 1.0; M[2][1] =-1.0; M[2][2] = 2.0;  //must be symmetric.

// Now create an instance of Jacobi ("eigen_calc").

Jacobi<double, double*, double**> eigen_calc(n);

// Note:
// If the matrix you plan to diagonalize (M) is read-only, use this instead:
//   Jacobi<double, double*, double**, double const*const*> eigen_calc(n);
// If you prefer using C++ vectors over C-style pointers, this works also:
//   Jacobi<double, vector<double>&, vector<vector<double>>&,
//          const vector<vector<double>>&>  eigen_calc(n);

// Now, calculate the eigenvalues and eigenvectors of M

eigen_calc.Diagonalize(M, evals, evecs);  //(successful if return value is > 0)

// If you have many matrices to diagonalize, you can re-use "eigen_calc". (This
// is more efficient than creating a new "Jacobi" class instance for each use.)

std::cout << "eigenvalues:  ";
for (int i=0; i < n; i++)
  cout << evals[i] << " ";
cout << endl;
for (int i=0; i < n; i++) {
  cout << "eigenvector" <<i+1<< ": ";
  for (int j=0; j < n; j++)
    cout << evecs[i][j] << " ";
  cout << endl;
}

Benchmarks

(details here...)

Installation

Copy the file(s) in the include subdirectory, to a location in your include path. No linking is necessary. This is a header-only library.

Development Status: stable

jacobi_pd has been tested for accuracy and memory safety over a wide range of array types, matrix sizes, eigenvalue magnitudes and degeneracies. jacobi_pd code is currently used in the popular LAMMPS and colvars MD simulation tools.

Requirements

A C++11 compatible compiler.

License

jacobi_pd is available under the terms of the Creative-Commons-Zero license.

Please send me corrections or suggestions.

jacobi_pd's Issues

Function fails for example matrix

If given the matrix:

[1, 2, 3]
[2, 5, 6]
[3, 6, 9]

Diagonalize will return the eigen values: 0.0, 0.699265, 14.300735, when they should be 13.0, 0, 0.

I'm not sure why the code is failing on this matrix specifically, as this matrix meets all the criteria for a Jacobi eigensolver (square, symmetric, dense).

The issue appears to be an issue with the application of the rotation matrix for this problem. I've compared against a python implementation of the Jacobi method (in addition to numpys np.linalg.eig). The values for cos(θ), sin(θ), and tan(θ) are the same on the first rotation ( 0.8112421851755608 0.584710284663765 0.7207592200561265, respectively) , but the outcome of the rotation is different. The python code gets:

[ 1,  0,  3],
[ 2,  0,  0],
[ 3,  6, 13]

while this code gets:

[1.000000, -0.131646, 3.603147]
[2.000000, 0.675445, 0.000000]
[0.000000, 0.000000, 13.324555]

(ignoring the lower triangular portion obviously).

For running this C++ code, I simply added:

  double sym_mat[3][3] = {
  {1, 2, 3},
  {2, 5, 6},
  {3, 6, 9} };
  double eigvals[3], eigvecs[3][3];
  const int NF = 3;
        int n = 3;

  Jacobi<double,
         double*,
         double (*)[NF],
         double const (*)[NF]> ecalc(n);


  ecalc.Diagonalize(sym_mat,
                    eigvals,
                    eigvecs,
                    Jacobi<double,
                            double*,
                            double (*)[NF],
                            double const (*)[NF]>::DO_NOT_SORT);

  // TestJacobi(n_size, n_matr, emin, emax, n_tests, n_degeneracy, seed, eps);

  printf("\n");
  printf("eigvals: [ %f, %f, %f]\n", eigvals[0], eigvals[1], eigvals[2]);
  printf("eigvecs: [[%f, %f, %f]\n[%f, %f, %f]\n[%f, %f, %f]]\n", eigvecs[0][0], eigvecs[0][1], eigvecs[0][2], eigvecs[1][0], eigvecs[1][1], eigvecs[1][2], eigvecs[2][0], eigvecs[2][1], eigvecs[2][2]);

to the end of test_jacobi.cpp:main. I was able to reproduce the results when manually translating the C++ to C as well (which was my intended application to begin with).