Repository files navigation

This repository contains a small C++header filethat implements theJacobi eigenvalue algorithm.It isfree of copyright.

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

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

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

As of late 2020, no simple public domain C++11 codeyet exists for matrix diagonalization.Other C++ libraries such as Eigen or GSL are typicallymuch larger and use more restrictive licenses.(On several occasions, this has prevented me from includingtheir code in other open-source projects with incompatible licenses.)Some repositories may unwittingly contain codesnippets from other sources, such asnumerical 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 reasonablyshort, simple,fast andreliable.You can do anything you like with this code.

#include"jacobi_pd.hpp"// ...int n =3;// Matrix sizedouble **M;// A symmetric n x n matrix you want to diagonalizedouble *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 matrixM[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_pd::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 Meigen_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;}

(details here...)

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

jacobi_pd has beentestedfor accuracy and memory safetyover a wide range of array types, matrix sizes,eigenvalue magnitudes and degeneracies.jacobi_pd code is currently used in the popularLAMMPSandcolvarsMD simulation tools.

A C++11 compatible compiler.

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

Please send me corrections or suggestions.