ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

updated all code in LAMMPS to use the open-source "math_eigen.h" instead of "math_extra.h". Code in "lib" now uses its own abbreviated local version of the "math_eigen.h" file (which is named "jacobi_pd.h"), since it is not clear how code from "lib/" can access the code in "src/"

Browse Source
This commit is contained in:
Andrew Jewett
2020-09-05 01:39:27 -07:00
parent 7f6fc8a003
commit 57f82abae3
20 changed files with 1299 additions and 768 deletions
Download Patch File Download Diff File Expand all files Collapse all files

197
lib/colvars/colvartypes.cpp
View File

@ -10,6 +10,7 @@
#include <cstdlib>
#include <cstring>
#include "jacobi_pd.h"
#include "colvarmodule.h"
#include "colvartypes.h"
#include "colvarparse.h"
@ -19,20 +20,6 @@ bool colvarmodule::rotation::monitor_crossings = false;
cvm::real colvarmodule::rotation::crossing_threshold = 1.0E-02;
namespace {
/// Numerical recipes diagonalization
static int jacobi(cvm::real **a, cvm::real *d, cvm::real **v, int *nrot);
/// Eigenvector sort
static int eigsrt(cvm::real *d, cvm::real **v);
/// Transpose the matrix
static int transpose(cvm::real **v);
}
std::string cvm::rvector::to_simple_string() const
{
std::ostringstream os;
@ -291,43 +278,6 @@ void colvarmodule::rotation::compute_overlap_matrix()
}
void colvarmodule::rotation::diagonalize_matrix(
cvm::matrix2d<cvm::real> &m,
cvm::vector1d<cvm::real> &eigval,
cvm::matrix2d<cvm::real> &eigvec)
{
eigval.resize(4);
eigval.reset();
eigvec.resize(4, 4);
eigvec.reset();
// diagonalize
int jac_nrot = 0;
if (jacobi(m.c_array(), eigval.c_array(), eigvec.c_array(), &jac_nrot) !=
COLVARS_OK) {
cvm::error("Too many iterations in routine jacobi.\n"
"This is usually the result of an ill-defined set of atoms for "
"rotational alignment (RMSD, rotateReference, etc).\n");
}
eigsrt(eigval.c_array(), eigvec.c_array());
// jacobi saves eigenvectors by columns
transpose(eigvec.c_array());
// normalize eigenvectors
for (size_t ie = 0; ie < 4; ie++) {
cvm::real norm2 = 0.0;
size_t i;
for (i = 0; i < 4; i++) {
norm2 += eigvec[ie][i] * eigvec[ie][i];
}
cvm::real const norm = cvm::sqrt(norm2);
for (i = 0; i < 4; i++) {
eigvec[ie][i] /= norm;
}
}
}
// Calculate the rotation, plus its derivatives
void colvarmodule::rotation::calc_optimal_rotation(
@ -349,7 +299,13 @@ void colvarmodule::rotation::calc_optimal_rotation(
cvm::log("S = "+cvm::to_str(S_backup, cvm::cv_width, cvm::cv_prec)+"\n");
}
diagonalize_matrix(S, S_eigval, S_eigvec);
int ierror = ecalc.Diagonalize(S, S_eigval, S_eigvec);
if (ierror) {
cvm::error("Too many iterations in routine jacobi.\n"
"This is usually the result of an ill-defined set of atoms for "
"rotational alignment (RMSD, rotateReference, etc).\n");
}
// eigenvalues and eigenvectors
cvm::real const L0 = S_eigval[0];
cvm::real const L1 = S_eigval[1];
@ -521,7 +477,7 @@ void colvarmodule::rotation::calc_optimal_rotation(
// cvm::log("S_new = "+cvm::to_str(cvm::to_str (S_new), cvm::cv_width, cvm::cv_prec)+"\n");
diagonalize_matrix(S_new, S_new_eigval, S_new_eigvec);
ecalc.Diagonalize(S_new, S_new_eigval, S_new_eigvec);
cvm::real const &L0_new = S_new_eigval[0];
cvm::quaternion const Q0_new(S_new_eigvec[0]);
@ -544,138 +500,3 @@ void colvarmodule::rotation::calc_optimal_rotation(
// Numerical Recipes routine for diagonalization
#define ROTATE(a,i,j,k,l) g=a[i][j]; \
h=a[k][l]; \
a[i][j]=g-s*(h+g*tau); \
a[k][l]=h+s*(g-h*tau);
#define n 4
namespace {
int jacobi(cvm::real **a, cvm::real *d, cvm::real **v, int *nrot)
{
int j,iq,ip,i;
cvm::real tresh,theta,tau,t,sm,s,h,g,c;
cvm::vector1d<cvm::real> b(n);
cvm::vector1d<cvm::real> z(n);
for (ip=0;ip<n;ip++) {
for (iq=0;iq<n;iq++) {
v[ip][iq]=0.0;
}
v[ip][ip]=1.0;
}
for (ip=0;ip<n;ip++) {
b[ip]=d[ip]=a[ip][ip];
z[ip]=0.0;
}
*nrot=0;
for (i=0;i<=50;i++) {
sm=0.0;
for (ip=0;ip<n-1;ip++) {
for (iq=ip+1;iq<n;iq++)
sm += cvm::fabs(a[ip][iq]);
}
if (sm == 0.0) {
return COLVARS_OK;
}
if (i < 4)
tresh=0.2*sm/(n*n);
else
tresh=0.0;
for (ip=0;ip<n-1;ip++) {
for (iq=ip+1;iq<n;iq++) {
g=100.0*cvm::fabs(a[ip][iq]);
if (i > 4 && (cvm::real)(cvm::fabs(d[ip])+g) == (cvm::real)cvm::fabs(d[ip])
&& (cvm::real)(cvm::fabs(d[iq])+g) == (cvm::real)cvm::fabs(d[iq]))
a[ip][iq]=0.0;
else if (cvm::fabs(a[ip][iq]) > tresh) {
h=d[iq]-d[ip];
if ((cvm::real)(cvm::fabs(h)+g) == (cvm::real)cvm::fabs(h))
t=(a[ip][iq])/h;
else {
theta=0.5*h/(a[ip][iq]);
t=1.0/(cvm::fabs(theta)+cvm::sqrt(1.0+theta*theta));
if (theta < 0.0) t = -t;
}
c=1.0/cvm::sqrt(1+t*t);
s=t*c;
tau=s/(1.0+c);
h=t*a[ip][iq];
z[ip] -= h;
z[iq] += h;
d[ip] -= h;
d[iq] += h;
a[ip][iq]=0.0;
for (j=0;j<=ip-1;j++) {
ROTATE(a,j,ip,j,iq)
}
for (j=ip+1;j<=iq-1;j++) {
ROTATE(a,ip,j,j,iq)
}
for (j=iq+1;j<n;j++) {
ROTATE(a,ip,j,iq,j)
}
for (j=0;j<n;j++) {
ROTATE(v,j,ip,j,iq)
}
++(*nrot);
}
}
}
for (ip=0;ip<n;ip++) {
b[ip] += z[ip];
d[ip]=b[ip];
z[ip]=0.0;
}
}
return COLVARS_ERROR;
}
int eigsrt(cvm::real *d, cvm::real **v)
{
int k,j,i;
cvm::real p;
for (i=0;i<n;i++) {
p=d[k=i];
for (j=i+1;j<n;j++)
if (d[j] >= p) p=d[k=j];
if (k != i) {
d[k]=d[i];
d[i]=p;
for (j=0;j<n;j++) {
p=v[j][i];
v[j][i]=v[j][k];
v[j][k]=p;
}
}
}
return COLVARS_OK;
}
int transpose(cvm::real **v)
{
cvm::real p;
int i,j;
for (i=0;i<n;i++) {
for (j=i+1;j<n;j++) {
p=v[i][j];
v[i][j]=v[j][i];
v[j][i]=p;
}
}
return COLVARS_OK;
}
}
#undef n
#undef ROTATE

14
lib/colvars/colvartypes.h
View File

@ -12,6 +12,7 @@
#include <vector>
#include "jacobi_pd.h"
#include "colvarmodule.h"
#ifndef PI
@ -1390,19 +1391,20 @@ public:
/// Default constructor
inline rotation()
: b_debug_gradients(false)
: b_debug_gradients(false), ecalc(4)
{}
/// Constructor after a quaternion
inline rotation(cvm::quaternion const &qi)
: q(qi),
b_debug_gradients(false)
b_debug_gradients(false),
ecalc(4)
{
}
/// Constructor after an axis of rotation and an angle (in radians)
inline rotation(cvm::real angle, cvm::rvector const &axis)
: b_debug_gradients(false)
: b_debug_gradients(false), ecalc(4)
{
cvm::rvector const axis_n = axis.unit();
cvm::real const sina = cvm::sin(angle/2.0);
@ -1538,9 +1540,9 @@ protected:
void compute_overlap_matrix();
/// Diagonalize a given matrix m (used by calc_optimal_rotation())
static void diagonalize_matrix(cvm::matrix2d<cvm::real> &m,
cvm::vector1d<cvm::real> &eigval,
cvm::matrix2d<cvm::real> &eigvec);
MathEigen::Jacobi<cvm::real,
cvm::vector1d<cvm::real> &,
cvm::matrix2d<cvm::real> &> ecalc;
};

623
lib/jacobi_pd.h Normal file
View File

@ -0,0 +1,623 @@
/// @file jacobi_pd.h
/// @brief Calculate the eigenvalues and eigevectors of a symmetric matrix
/// using the Jacobi eigenvalue algorithm.
/// @author Andrew Jewett
/// @note Benchmarks and tests: https://github.com/jewettaij/jacobi_pd
/// @note A version of this code is included with LAMMPS ("math_eigen.h").
/// @license CC0-1.0 (public domain)
#ifndef _MATH_EIGEN_H
#define _MATH_EIGEN_H
#include <algorithm>
#include <cmath>
//#include <cassert>
namespace MathEigen {
/// @brief Allocate a 2-dimensional array. (Uses row-major order.)
/// @note (Matrix allocation functions can also be found in colvartypes.h)
template<typename Entry>
void Alloc2D(size_t nrows, //!< size of the array (number of rows)
size_t ncols, //!< size of the array (number of columns)
Entry ***paaX //!< pointer to a 2D C-style array
);
/// @brief Deallocate arrays that were created using Alloc2D().
/// @note (Matrix allocation functions can also be found in colvartypes.h)
template<typename Entry>
void Dealloc2D(Entry ***paaX //!< pointer to a 2D C-style array
);
// ---- Eigendecomposition of small dense symmetric matrices ----
/// @class Jacobi
/// @brief Calculate the eigenvalues and eigevectors of a symmetric matrix
/// using the Jacobi eigenvalue algorithm. Code for the Jacobi class
/// (along with tests and benchmarks) is available free of copyright at
/// https://github.com/jewettaij/jacobi_pd
/// @note The "Vector" and "Matrix" type arguments can be any
/// C or C++ object that support indexing, including pointers or vectors.
/// @details
/// Usage example:
/// @code
///
/// int n = 5; // Matrix size
/// double **M; // A symmetric n x n matrix you want to diagonalize
/// double *evals; // Store the eigenvalues here.
/// double **evects; // Store the eigenvectors here.
/// // Allocate space for M, evals, and evects, and load contents of M (omitted)
///
/// // Now create an instance of Jacobi ("eigen_calc"). This will allocate space
/// // for storing intermediate calculations. Once created, it can be reused
/// // multiple times without paying the cost of allocating memory on the heap.
///
/// 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 vectors over C-style pointers, this works also:
/// // Jacobi<double, vector<double>&, vector<vector<double>>&> eigen_calc(n);
///
/// // Now, calculate the eigenvalues and eigenvectors of M
///
/// eigen_calc.Diagonalize(M, evals, evects);
///
/// @endcode
template<typename Scalar,
typename Vector,
typename Matrix,
typename ConstMatrix=Matrix>
class Jacobi
{
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 the eigenvalues and eigevectors of a symmetric matrix
/// using the Jacobi eigenvalue algorithm.
/// @returns 0 if the algorithm converged,
/// 1 if the algorithm failed to converge. (IE, if the number of
/// pivot iterations exceeded max_num_sweeps * iters_per_sweep,
/// where iters_per_sweep = (n*(n-1)/2))
/// @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 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(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);
Jacobi<Scalar, Vector, Matrix, ConstMatrix>& operator = (Jacobi<Scalar, Vector, Matrix, ConstMatrix> source);
}; // class Jacobi
// -------------- IMPLEMENTATION --------------
// Utilities for memory allocation
template<typename Entry>
void Alloc2D(size_t nrows, // size of the array (number of rows)
size_t ncols, // size of the array (number of columns)
Entry ***paaX) // pointer to a 2D C-style array
{
//assert(paaX);
*paaX = new Entry* [nrows]; //conventional 2D C array (pointer-to-pointer)
(*paaX)[0] = new Entry [nrows * ncols]; // 1D C array (contiguous memory)
for(size_t iy=0; iy<nrows; iy++)
(*paaX)[iy] = (*paaX)[0] + iy*ncols;
// The caller can access the contents using (*paaX)[i][j]
}
template<typename Entry>
void Dealloc2D(Entry ***paaX) // pointer to a 2D C-style array
{
if (paaX && *paaX) {
delete [] (*paaX)[0];
delete [] (*paaX);
*paaX = nullptr;
}
}
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
Matrix evec, // store the eigenvectors here (in rows)
SortCriteria sort_criteria, // sort results?
bool calc_evec, // calculate the eigenvectors?
int max_num_sweeps) // limit the number of iterations ("sweeps")
{
// -- Initialization --
for (int i = 0; i < n; i++)
for (int j = i; j < n; j++) //copy mat[][] into M[][]
M[i][j] = mat[i][j]; //(M[][] is a local copy we can modify)
if (calc_evec)
for (int i = 0; i < n; i++)
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
max_idx_row[i] = MaxEntryRow(M, i); //(which is needed by MaxEntry())
// -- Iteration --
int n_iters;
int max_num_iters = max_num_sweeps*n*(n-1)/2; //"sweep" = n*(n-1)/2 iters
for (n_iters=0; n_iters < max_num_iters; n_iters++) {
int i,j;
MaxEntry(M, i, j); // Find the maximum entry in the matrix. Store in i,j
// If M[i][j] is small compared to M[i][i] and M[j][j], set it to 0.
if ((M[i][i] + M[i][j] == M[i][i]) && (M[j][j] + M[i][j] == M[j][j])) {
M[i][j] = 0.0;
max_idx_row[i] = MaxEntryRow(M,i); //must also update max_idx_row[i]
}
if (M[i][j] == 0.0)
break;
// Otherwise, apply a rotation to make M[i][j] = 0
CalcRot(M, i, j); // Calculate the parameters of the rotation matrix.
ApplyRot(M, i, j); // Apply this rotation to the M matrix.
if (calc_evec) // Optional: If the caller wants the eigenvectors, then
ApplyRotLeft(evec,i,j); // apply the rotation to the eigenvector matrix
} //for (int n_iters=0; n_iters < max_num_iters; n_iters++)
// -- Post-processing --
for (int i = 0; i < n; i++)
eval[i] = M[i][i];
// Optional: Sort results by eigenvalue.
SortRows(eval, evec, n, sort_criteria);
return (n_iters == max_num_iters);
}
/// @brief Calculate the components of a rotation matrix which performs a
/// rotation in the i,j plane by an angle (θ) that (when multiplied on
/// 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).
template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>
void Jacobi<Scalar, Vector, Matrix, ConstMatrix>::
CalcRot(Scalar const *const *M, // matrix
int i, // row index
int j) // column index
{
t = 1.0; // = tan(θ)
Scalar M_jj_ii = (M[j][j] - M[i][i]);
if (M_jj_ii != 0.0) {
// kappa = (M[j][j] - M[i][i]) / (2*M[i][j])
Scalar kappa = M_jj_ii;
t = 0.0;
Scalar M_ij = M[i][j];
if (M_ij != 0.0) {
kappa /= (2.0*M_ij);
// t satisfies: t^2 + 2*t*kappa - 1 = 0
// (choose the root which has the smaller absolute value)
t = 1.0 / (std::sqrt(1 + kappa*kappa) + std::abs(kappa));
if (kappa < 0.0)
t = -t;
}
}
//assert(std::abs(t) <= 1.0);
c = 1.0 / std::sqrt(1 + t*t);
s = c*t;
}
/// brief Perform a similarity transformation by multiplying matrix M on both
/// sides by a rotation matrix (and its transpose) to eliminate M[i][j].
/// details This rotation matrix performs a rotation in the i,j plane by
/// angle θ. This function assumes that c=cos(θ). s=som(θ), t=tan(θ)
/// have been calculated in advance (using the CalcRot() function).
/// It also assumes that i<j. The max_idx_row[] array is also updated.
/// To save time, since the matrix is symmetric, the elements
/// below the diagonal (ie. M[u][v] where u>v) are not computed.
/// verbatim
/// M' = R^T * M * R
/// where R the rotation in the i,j plane and ^T denotes the transpose.
/// i j
/// _ _
/// | 1 |
/// | . |
/// | . |
/// | 1 |
/// | c ... s |
/// | . . . |
/// R = | . 1 . |
/// | . . . |
/// | -s ... c |
/// | 1 |
/// | . |
/// | . |
/// |_ 1 _|
/// endverbatim
///
/// Let M' denote the matrix M after multiplication by R^T and R.
/// The components of M' are:
///
/// verbatim
/// M'_uv = Σ_w Σ_z R_wu * M_wz * R_zv
/// endverbatim
///
/// Note that a the rotation at location i,j will modify all of the matrix
/// elements containing at least one index which is either i or j
/// such as: M[w][i], M[i][w], M[w][j], M[j][w].
/// Check and see whether these modified matrix elements exceed the
/// corresponding values in max_idx_row[] array for that row.
/// If so, then update max_idx_row for that row.
/// This is somewhat complicated by the fact that we must only consider
/// matrix elements in the upper-right triangle strictly above the diagonal.
/// (ie. matrix elements whose second index is > the first index).
/// The modified elements we must consider are marked with an "X" below:
///
/// @verbatim
/// i j
/// _ _
/// | . X X |
/// | . X X |
/// | . X X |
/// | . X X |
/// | X X X X X 0 X X X X | i
/// | . X |
/// | . X |
/// M = | . X |
/// | . X |
/// | X X X X X | j
/// | . |
/// | . |
/// | . |
/// |_ . _|
/// @endverbatim
template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>
void Jacobi<Scalar, Vector, Matrix, ConstMatrix>::
ApplyRot(Scalar **M, // matrix
int i, // row index
int j) // column index
{
// Recall that:
// c = cos(θ)
// s = sin(θ)
// t = tan(θ) (which should be <= 1.0)
// Compute the diagonal elements of M which have changed:
M[i][i] -= t * M[i][j];
M[j][j] += t * M[i][j];
// Note: This is algebraically equivalent to:
// M[i][i] = c*c*M[i][i] + s*s*M[j][j] - 2*s*c*M[i][j]
// M[j][j] = s*s*M[i][i] + c*c*M[j][j] + 2*s*c*M[i][j]
//Update the off-diagonal elements of M which will change (above the diagonal)
//assert(i < j);
M[i][j] = 0.0;
//compute M[w][i] and M[i][w] for all w!=i,considering above-diagonal elements
for (int w=0; w < i; w++) { // 0 <= w < i < j < n
M[i][w] = M[w][i]; //backup the previous value. store below diagonal (i>w)
M[w][i] = c*M[w][i] - s*M[w][j]; //M[w][i], M[w][j] from previous iteration
if (i == max_idx_row[w]) max_idx_row[w] = MaxEntryRow(M, w);
else if (std::abs(M[w][i])>std::abs(M[w][max_idx_row[w]])) max_idx_row[w]=i;
//assert(max_idx_row[w] == MaxEntryRow(M, w));
}
for (int w=i+1; w < j; w++) { // 0 <= i < w < j < n
M[w][i] = M[i][w]; //backup the previous value. store below diagonal (w>i)
M[i][w] = c*M[i][w] - s*M[w][j]; //M[i][w], M[w][j] from previous iteration
}
for (int w=j+1; w < n; w++) { // 0 <= i < j+1 <= w < n
M[w][i] = M[i][w]; //backup the previous value. store below diagonal (w>i)
M[i][w] = c*M[i][w] - s*M[j][w]; //M[i][w], M[j][w] from previous iteration
}
// now that we're done modifying row i, we can update max_idx_row[i]
max_idx_row[i] = MaxEntryRow(M, i);
//compute M[w][j] and M[j][w] for all w!=j,considering above-diagonal elements
for (int w=0; w < i; w++) { // 0 <= w < i < j < n
M[w][j] = s*M[i][w] + c*M[w][j]; //M[i][w], M[w][j] from previous iteration
if (j == max_idx_row[w]) max_idx_row[w] = MaxEntryRow(M, w);
else if (std::abs(M[w][j])>std::abs(M[w][max_idx_row[w]])) max_idx_row[w]=j;
//assert(max_idx_row[w] == MaxEntryRow(M, w));
}
for (int w=i+1; w < j; w++) { // 0 <= i+1 <= w < j < n
M[w][j] = s*M[w][i] + c*M[w][j]; //M[w][i], M[w][j] from previous iteration
if (j == max_idx_row[w]) max_idx_row[w] = MaxEntryRow(M, w);
else if (std::abs(M[w][j])>std::abs(M[w][max_idx_row[w]])) max_idx_row[w]=j;
//assert(max_idx_row[w] == MaxEntryRow(M, w));
}
for (int w=j+1; w < n; w++) { // 0 <= i < j < w < n
M[j][w] = s*M[w][i] + c*M[j][w]; //M[w][i], M[j][w] from previous iteration
}
// now that we're done modifying row j, we can update max_idx_row[j]
max_idx_row[j] = MaxEntryRow(M, j);
} //Jacobi::ApplyRot()
///@brief Multiply matrix M on the LEFT side by a transposed rotation matrix R^T
/// This matrix performs a rotation in the i,j plane by angle θ (where
/// the arguments "s" and "c" refer to cos(θ) and sin(θ), respectively).
/// @verbatim
/// E'_uv = Σ_w R_wu * E_wv
/// @endverbatim
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
{
// Recall that c = cos(θ) and s = sin(θ)
for (int v = 0; v < n; v++) {
Scalar Eiv = E[i][v]; //backup E[i][v]
E[i][v] = c*E[i][v] - s*E[j][v];
E[j][v] = s*Eiv + c*E[j][v];
}
}
template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>
int Jacobi<Scalar, Vector, Matrix, ConstMatrix>::
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]))
j_max = j;
return j_max;
}
template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>
void Jacobi<Scalar, Vector, Matrix, ConstMatrix>::
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;
j_max = max_idx_row[i_max];
Scalar max_entry = std::abs(M[i_max][j_max]);
int nm1 = n-1;
for (int i=1; i < nm1; i++) {
int j = max_idx_row[i];
if (std::abs(M[i][j]) > max_entry) {
max_entry = std::abs(M[i][j]);
i_max = i;
j_max = j;
}
}
//#ifndef NDEBUG
//// make sure that the maximum element really is stored at i_max, j_max
//// (comment out the next loop before using. it slows down the code)
//for (int i = 0; i < nm1; i++)
// for (int j = i+1; j < n; j++)
// assert(std::abs(M[i][j]) <= max_entry);
//#endif
}
//Sort the rows of a matrix "evec" by the numbers contained in "eval"
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++) {
int i_max = i;
for (int j = i+1; j < n; j++) {
// find the "maximum" element in the array starting at position i+1
switch (sort_criteria) {
case SORT_DECREASING_EVALS:
if (eval[j] > eval[i_max])
i_max = j;
break;
case SORT_INCREASING_EVALS:
if (eval[j] < eval[i_max])
i_max = j;
break;
case SORT_DECREASING_ABS_EVALS:
if (std::abs(eval[j]) > std::abs(eval[i_max]))
i_max = j;
break;
case SORT_INCREASING_ABS_EVALS:
if (std::abs(eval[j]) < std::abs(eval[i_max]))
i_max = j;
break;
default:
break;
}
}
std::swap(eval[i], eval[i_max]); // sort "eval"
for (int k = 0; k < n; k++)
std::swap(evec[i][k], evec[i_max][k]); // sort "evec"
}
}
template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>
void Jacobi<Scalar, Vector, Matrix, ConstMatrix>::
Init() {
n = 0;
M = nullptr;
max_idx_row = nullptr;
}
template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>
void Jacobi<Scalar, Vector, Matrix, ConstMatrix>::
SetSize(int n) {
Dealloc();
Alloc(n);
}
// memory management:
template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>
void Jacobi<Scalar, Vector, Matrix, ConstMatrix>::
Alloc(int n) {
this->n = n;
if (n > 0) {
max_idx_row = new int[n];
Alloc2D(n, n, &M);
}
}
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();
}
// memory management: copy and move constructor, swap, and assignment operator:
template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>
Jacobi<Scalar, Vector, Matrix, ConstMatrix>::
Jacobi(const Jacobi<Scalar, Vector, Matrix, ConstMatrix>& source)
{
Init();
SetSize(source.n);
//assert(n == source.n);
// The following lines aren't really necessary, because the contents
// of source.M and source.max_idx_row are not needed (since they are
// overwritten every time Jacobi::Diagonalize() is invoked).
std::copy(source.max_idx_row,
source.max_idx_row + n,
max_idx_row);
for (int i = 0; i < n; i++)
std::copy(source.M[i],
source.M[i] + n,
M[i]);
}
template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>
void Jacobi<Scalar, Vector, Matrix, ConstMatrix>::
swap(Jacobi<Scalar, Vector, Matrix, ConstMatrix> &other) {
std::swap(n, other.n);
std::swap(max_idx_row, other.max_idx_row);
std::swap(M, other.M);
}
// Move constructor (C++11)
template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>
Jacobi<Scalar, Vector, Matrix, ConstMatrix>::
Jacobi(Jacobi<Scalar, Vector, Matrix, ConstMatrix>&& other) {
Init();
this->swap(other);
}
// Using the "copy-swap" idiom for the assignment operator
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;
}
} // namespace MathEigen
#endif //#ifndef _MATH_EIGEN_H

3
src/BODY/body_nparticle.cpp
View File

@ -16,6 +16,7 @@
#include <cstdlib>
#include "my_pool_chunk.h"
#include "math_extra.h"
#include "math_eigen.h"
#include "atom_vec_body.h"
#include "atom.h"
#include "force.h"
@ -137,7 +138,7 @@ void BodyNparticle::data_body(int ibonus, int ninteger, int ndouble,
double *inertia = bonus->inertia;
double evectors[3][3];
int ierror = MathExtra::jacobi(tensor,inertia,evectors);
int ierror = MathEigen::jacobi3(tensor,inertia,evectors);
if (ierror) error->one(FLERR,
"Insufficient Jacobi rotations for body nparticle");

3
src/BODY/body_rounded_polygon.cpp
View File

@ -24,6 +24,7 @@
#include "force.h"
#include "domain.h"
#include "math_extra.h"
#include "math_eigen.h"
#include "memory.h"
#include "error.h"
#include "fmt/format.h"
@ -198,7 +199,7 @@ void BodyRoundedPolygon::data_body(int ibonus, int ninteger, int ndouble,
double *inertia = bonus->inertia;
double evectors[3][3];
int ierror = MathExtra::jacobi(tensor,inertia,evectors);
int ierror = MathEigen::jacobi3(tensor,inertia,evectors);
if (ierror) error->one(FLERR,
"Insufficient Jacobi rotations for body nparticle");

3
src/BODY/body_rounded_polyhedron.cpp
View File

@ -24,6 +24,7 @@
#include "atom.h"
#include "force.h"
#include "math_extra.h"
#include "math_eigen.h"
#include "memory.h"
#include "error.h"
#include "fmt/format.h"
@ -243,7 +244,7 @@ void BodyRoundedPolyhedron::data_body(int ibonus, int ninteger, int ndouble,
double *inertia = bonus->inertia;
double evectors[3][3];
int ierror = MathExtra::jacobi(tensor,inertia,evectors);
int ierror = MathEigen::jacobi3(tensor,inertia,evectors);
if (ierror) error->one(FLERR,
"Insufficient Jacobi rotations for body nparticle");

86
src/POEMS/fix_poems.cpp
View File

@ -36,6 +36,7 @@
#include "error.h"
#include "utils.h"
#include "fmt/format.h"
#include "math_eigen.h"
using namespace LAMMPS_NS;
using namespace FixConst;
@ -44,7 +45,6 @@ using namespace FixConst;
#define DELTA 128
#define TOLERANCE 1.0e-6
#define EPSILON 1.0e-7
#define MAXJACOBI 50
static const char cite_fix_poems[] =
"fix poems command:\n\n"
@ -492,7 +492,7 @@ void FixPOEMS::init()
tensor[1][2] = tensor[2][1] = all[ibody][4];
tensor[0][2] = tensor[2][0] = all[ibody][5];
ierror = jacobi(tensor,inertia[ibody],evectors);
ierror = MathEigen::jacobi3(tensor,inertia[ibody],evectors);
if (ierror) error->all(FLERR,"Insufficient Jacobi rotations for POEMS body");
ex_space[ibody][0] = evectors[0][0];
@ -1283,88 +1283,6 @@ int FixPOEMS::loopcheck(int nvert, int nedge, tagint **elist)
return 0;
}
/* ----------------------------------------------------------------------
compute evalues and evectors of 3x3 real symmetric matrix
based on Jacobi rotations
adapted from Numerical Recipes jacobi() function
------------------------------------------------------------------------- */
int FixPOEMS::jacobi(double **matrix, double *evalues, double **evectors)
{
int i,j,k;
double tresh,theta,tau,t,sm,s,h,g,c,b[3],z[3];
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) evectors[i][j] = 0.0;
evectors[i][i] = 1.0;
}
for (i = 0; i < 3; i++) {
b[i] = evalues[i] = matrix[i][i];
z[i] = 0.0;
}
for (int iter = 1; iter <= MAXJACOBI; iter++) {
sm = 0.0;
for (i = 0; i < 2; i++)
for (j = i+1; j < 3; j++)
sm += fabs(matrix[i][j]);
if (sm == 0.0) return 0;
if (iter < 4) tresh = 0.2*sm/(3*3);
else tresh = 0.0;
for (i = 0; i < 2; i++) {
for (j = i+1; j < 3; j++) {
g = 100.0*fabs(matrix[i][j]);
if (iter > 4 && fabs(evalues[i])+g == fabs(evalues[i])
&& fabs(evalues[j])+g == fabs(evalues[j]))
matrix[i][j] = 0.0;
else if (fabs(matrix[i][j]) > tresh) {
h = evalues[j]-evalues[i];
if (fabs(h)+g == fabs(h)) t = (matrix[i][j])/h;
else {
theta = 0.5*h/(matrix[i][j]);
t = 1.0/(fabs(theta)+sqrt(1.0+theta*theta));
if (theta < 0.0) t = -t;
}
c = 1.0/sqrt(1.0+t*t);
s = t*c;
tau = s/(1.0+c);
h = t*matrix[i][j];
z[i] -= h;
z[j] += h;
evalues[i] -= h;
evalues[j] += h;
matrix[i][j] = 0.0;
for (k = 0; k < i; k++) rotate(matrix,k,i,k,j,s,tau);
for (k = i+1; k < j; k++) rotate(matrix,i,k,k,j,s,tau);
for (k = j+1; k < 3; k++) rotate(matrix,i,k,j,k,s,tau);
for (k = 0; k < 3; k++) rotate(evectors,k,i,k,j,s,tau);
}
}
}
for (i = 0; i < 3; i++) {
evalues[i] = b[i] += z[i];
z[i] = 0.0;
}
}
return 1;
}
/* ----------------------------------------------------------------------
perform a single Jacobi rotation
------------------------------------------------------------------------- */
void FixPOEMS::rotate(double **matrix, int i, int j, int k, int l,
double s, double tau)
{
double g = matrix[i][j];
double h = matrix[k][l];
matrix[i][j] = g-s*(h+g*tau);
matrix[k][l] = h+s*(g-h*tau);
}
/* ----------------------------------------------------------------------
compute omega from angular momentum
w = omega = angular velocity in space frame

2
src/POEMS/fix_poems.h
View File

@ -106,8 +106,6 @@ class FixPOEMS : public Fix {
void jointbuild();
void sortlist(int, tagint **);
int loopcheck(int, int, tagint **);
int jacobi(double **, double *, double **);
void rotate(double **, int, int, int, int, double, double);
void omega_from_mq(double *, double *, double *, double *,
double *, double *);
void set_v();

3
src/RIGID/fix_rigid.cpp
View File

@ -17,6 +17,7 @@
#include <cstdlib>
#include <cstring>
#include "math_extra.h"
#include "math_eigen.h"
#include "atom.h"
#include "atom_vec_ellipsoid.h"
#include "atom_vec_line.h"
@ -1948,7 +1949,7 @@ void FixRigid::setup_bodies_static()
tensor[0][2] = tensor[2][0] = all[ibody][4];
tensor[0][1] = tensor[1][0] = all[ibody][5];
ierror = MathExtra::jacobi(tensor,inertia[ibody],evectors);
ierror = MathEigen::jacobi3(tensor,inertia[ibody],evectors);
if (ierror) error->all(FLERR,
"Insufficient Jacobi rotations for rigid body");

3
src/RIGID/fix_rigid_small.cpp
View File

@ -18,6 +18,7 @@
#include <cstring>
#include <utility>
#include "math_extra.h"
#include "math_eigen.h"
#include "atom.h"
#include "atom_vec_ellipsoid.h"
#include "atom_vec_line.h"
@ -2085,7 +2086,7 @@ void FixRigidSmall::setup_bodies_static()
tensor[0][1] = tensor[1][0] = itensor[ibody][5];
inertia = body[ibody].inertia;
ierror = MathExtra::jacobi(tensor,inertia,evectors);
ierror = MathEigen::jacobi3(tensor,inertia,evectors);
if (ierror) error->all(FLERR,
"Insufficient Jacobi rotations for rigid body");

3
src/USER-MISC/compute_gyration_shape.cpp
View File

@ -21,6 +21,7 @@
#include <cstring>
#include "error.h"
#include "math_extra.h"
#include "math_eigen.h"
#include "math_special.h"
#include "modify.h"
#include "update.h"
@ -95,7 +96,7 @@ void ComputeGyrationShape::compute_vector()
ione[1][2] = ione[2][1] = gyration_tensor[4];
ione[0][2] = ione[2][0] = gyration_tensor[5];
int ierror = MathExtra::jacobi(ione,evalues,evectors);
int ierror = MathEigen::jacobi3(ione,evalues,evectors);
if (ierror) error->all(FLERR, "Insufficient Jacobi rotations "
"for gyration/shape");

3
src/USER-MISC/compute_gyration_shape_chunk.cpp
View File

@ -21,6 +21,7 @@
#include <cstring>
#include "error.h"
#include "math_extra.h"
#include "math_eigen.h"
#include "math_special.h"
#include "modify.h"
#include "memory.h"
@ -125,7 +126,7 @@ void ComputeGyrationShapeChunk::compute_array()
ione[0][2] = ione[2][0] = gyration_tensor[ichunk][4];
ione[1][2] = ione[2][1] = gyration_tensor[ichunk][5];
int ierror = MathExtra::jacobi(ione,evalues,evectors);
int ierror = MathEigen::jacobi3(ione,evalues,evectors);
if (ierror) error->all(FLERR, "Insufficient Jacobi rotations "
"for gyration/shape");

3
src/atom_vec_tri.cpp
View File

@ -15,6 +15,7 @@
#include <cmath>
#include <cstring>
#include "math_extra.h"
#include "math_eigen.h"
#include "atom.h"
#include "domain.h"
#include "modify.h"
@ -555,7 +556,7 @@ void AtomVecTri::data_atom_bonus(int m, char **values)
tensor[0][2] = tensor[2][0] = inertia[4];
tensor[0][1] = tensor[1][0] = inertia[5];
int ierror = MathExtra::jacobi(tensor,bonus[nlocal_bonus].inertia,evectors);
int ierror = MathEigen::jacobi3(tensor,bonus[nlocal_bonus].inertia,evectors);
if (ierror) error->one(FLERR,"Insufficient Jacobi rotations for triangle");
double ex_space[3],ey_space[3],ez_space[3];

5
src/compute_omega_chunk.cpp
View File

@ -20,6 +20,7 @@
#include "compute_chunk_atom.h"
#include "domain.h"
#include "math_extra.h"
#include "math_eigen.h"
#include "memory.h"
#include "error.h"
@ -250,10 +251,10 @@ void ComputeOmegaChunk::compute_array()
// handle each (nearly) singular I matrix
// due to 2-atom chunk or linear molecule
// use jacobi() and angmom_to_omega() to calculate valid omega
// use jacobi3() and angmom_to_omega() to calculate valid omega
} else {
int ierror = MathExtra::jacobi(ione,idiag,evectors);
int ierror = MathEigen::jacobi3(ione,idiag,evectors);
if (ierror) error->all(FLERR,
"Insufficient Jacobi rotations for omega/chunk");

5
src/group.cpp
View File

@ -30,6 +30,7 @@
#include "variable.h"
#include "dump.h"
#include "math_extra.h"
#include "math_eigen.h"
#include "memory.h"
#include "error.h"
#include "utils.h"
@ -1729,11 +1730,11 @@ void Group::omega(double *angmom, double inertia[3][3], double *w)
// handle (nearly) singular I matrix
// typically due to 2-atom group or linear molecule
// use jacobi() and angmom_to_omega() to calculate valid omega
// use jacobi3() and angmom_to_omega() to calculate valid omega
// less exact answer than matrix inversion, due to iterative Jacobi method
} else {
int ierror = MathExtra::jacobi(inertia,idiag,evectors);
int ierror = MathEigen::jacobi3(inertia, idiag, evectors);
if (ierror) error->all(FLERR,
"Insufficient Jacobi rotations for group::omega");

140
src/math_eigen.cpp Normal file
View File

@ -0,0 +1,140 @@
#include "math_eigen.h"
#include<array>
using std::vector;
using std::array;
using namespace MathEigen;
// --- Instantiate template class instances ---
// --- to reduce bloat in the compiled binary ---
// When using one of these versions of Jacobi, you can reduce code bloat by
// using an "extern template class" declaration in your cpp file. For example:
//
// #include"math_eigen.h"
// extern template class MathEigen::Jacobi<double, double*, double**, double const*const*>;
//
// ...This should (hopefully) use the version of Jacobi defined in this file
// instead of compiling a new version.
// template instantiations of Jacobi for pointer->pointer arrays
template class MathEigen::Jacobi<double, double*, double**, double const*const*>;
// template instantiations of Jacbi for fixed-size (3x3) arrays
template class MathEigen::Jacobi<double, double*, double (*)[3], double const (*)[3]>;
// template instantiations of Jacobi for vector-of-vectors:
template class MathEigen::Jacobi<double,
vector<double>&,
vector<vector<double> >&,
const vector<vector<double> >& >;
// template instantiations of Jacobi for array-of-arrays:
template class MathEigen::Jacobi<double,
array<double, 3>&,
array<array<double, 3>, 3 >&,
const array<array<double, 3>, 3 >& >;
// template instantiations of LambdaLanczos
template class MathEigen::LambdaLanczos<double>;
// template instantiations of PEidenDense for pointer->pointer arrays
template class MathEigen::PEigenDense<double, double*, double const*const*>;
// template instantiations of PEidenDense for vector<vector<>>
template class MathEigen::PEigenDense<double,
vector<double>&,
const vector<vector<double> >&>;
// ----------------- now create instances for floats ----------------
template class MathEigen::Jacobi<float, float*, float**, float const*const*>;
template class MathEigen::Jacobi<float, float*, float (*)[3], float const (*)[3]>;
template class MathEigen::Jacobi<float,
vector<float>&,
vector<vector<float> >&,
const vector<vector<float> >& >;
template class MathEigen::Jacobi<float,
array<float, 3>&,
array<array<float, 3>, 3 >&,
const array<array<float, 3>, 3 >& >;
template class MathEigen::LambdaLanczos<float>;
template class MathEigen::PEigenDense<float, float*, float const*const*>;
template class MathEigen::PEigenDense<float,
vector<float>&,
const vector<vector<float> >&>;
// If you plan to use other numeric types, add them below...
// Special case: 3x3 matrices
int MathEigen::
jacobi3(double const mat[3][3], // the 3x3 matrix you wish to diagonalize
double *eval, // store the eigenvalues here
double evec[3][3], // store the eigenvectors here...
bool evec_as_columns, // ...as rows or columns?
Jacobi<double,double*,double(*)[3],double const(*)[3]>::SortCriteria
sort_criteria)
{
// Create "ecalc3", an instance of the MathEigen::Jacobi class to calculate
// the eigenvalues oand eigenvectors of matrix "mat". It requires memory
// to be allocated to store a copy of the matrix. To avoid allocating
// this memory on the heap, we create a fixed size 3x3 array on the stack
// (and convert it to type double**).
double mat_cpy[3][3] = { {mat[0][0], mat[0][1], mat[0][2]},
{mat[1][0], mat[1][1], mat[1][2]},
{mat[2][0], mat[2][1], mat[2][2]} };
double *M[3] = { &(mat_cpy[0][0]), &(mat_cpy[1][0]), &(mat_cpy[2][0]) };
// variable "ecalc3" does all the work:
Jacobi<double,double*,double(*)[3],double const(*)[3]> ecalc3(3,M);
int ierror = ecalc3.Diagonalize(mat, eval, evec, sort_criteria);
// This will store the eigenvectors in the rows of "evec".
// Do we want to return the eigenvectors in columns instead of rows?
if (evec_as_columns)
for (int i=0; i<3; i++)
for (int j=0; j<3; j++)
std::swap(evec[i][j], evec[j][i]);
return ierror;
}
int MathEigen::
jacobi3(double const **mat, // the 3x3 matrix you wish to diagonalize
double *eval, // store the eigenvalues here
double **evec, // store the eigenvectors here...
bool evec_as_columns, // ...as rows or columns?
Jacobi<double,double*,double**,double const*const*>::SortCriteria
sort_criteria)
{
// Create "ecalc3", an instance of the MathEigen::Jacobi class to calculate
// the eigenvalues oand eigenvectors of matrix "mat". It requires memory
// to be allocated to store a copy of the matrix. To avoid allocating
// this memory on the heap, we create a fixed size 3x3 array on the stack
// (and convert it to type double**).
double mat_cpy[3][3] = { {mat[0][0], mat[0][1], mat[0][2]},
{mat[1][0], mat[1][1], mat[1][2]},
{mat[2][0], mat[2][1], mat[2][2]} };
double *M[3] = { &(mat_cpy[0][0]), &(mat_cpy[1][0]), &(mat_cpy[2][0]) };
// variable "ecalc3" does all the work:
Jacobi<double,double*,double**,double const*const*> ecalc3(3,M);
int ierror = ecalc3.Diagonalize(mat, eval, evec, sort_criteria);
// This will store the eigenvectors in the rows of "evec".
// Do we want to return the eigenvectors in columns instead of rows?
if (evec_as_columns)
for (int i=0; i<3; i++)
for (int j=0; j<3; j++)
std::swap(evec[i][j], evec[j][i]);
return ierror;
}

597
src/math_eigen.h
View File

@ -16,8 +16,8 @@
Andrew Jewett (Scripps Research, Jacobi algorithm)
------------------------------------------------------------------------- */
#ifndef _MATH_EIGEN_H
#define _MATH_EIGEN_H
#ifndef LMP_MATH_EIGEN_H
#define LMP_MATH_EIGEN_H
/// @file This file contains a library of functions and classes which can
/// efficiently perform eigendecomposition for an extremely broad
@ -40,7 +40,7 @@
/// https://github.com/jewettaij/jacobi_pd (CC0-1.0 license)
/// https://github.com/mrcdr/lambda-lanczos (MIT license)
#include <cassert>
//#include <cassert>
#include <numeric>
#include <complex>
#include <limits>
@ -51,118 +51,118 @@
namespace MathEigen {
// --- Memory allocation for matrices ---
// --- Memory allocation for matrices ---
/// @brief Allocate an arbitrary 2-dimensional array. (Uses row-major order.)
/// @note This function was intended for relatively small matrices (eg 4x4).
/// For large arrays, please use the 2d create() function from "memory.h"
template<typename Entry>
void Alloc2D(size_t nrows, //!< size of the array (number of rows)
/// @brief Allocate an arbitrary 2-dimensional array. (Uses row-major order.)
/// @note This function was intended for relatively small matrices (eg 4x4).
/// For large arrays, please use the 2d create() function from "memory.h"
template<typename Entry>
void Alloc2D(size_t nrows, //!< size of the array (number of rows)
size_t ncols, //!< size of the array (number of columns)
Entry ***paaX); //!< pointer to a 2D C-style array
/// @brief Deallocate arrays that were created using Alloc2D().
template<typename Entry>
void Dealloc2D(Entry ***paaX); //!< pointer to 2D multidimensional array
/// @brief Deallocate arrays that were created using Alloc2D().
template<typename Entry>
void Dealloc2D(Entry ***paaX); //!< pointer to 2D multidimensional array
// --- Complex numbers ---
// --- Complex numbers ---
/// @brief "realTypeMap" struct is used to the define "real_t<T>" type mapper
/// which returns the C++ type corresponding to the real component of T.
/// @details Consider a function ("l2_norm()") that calculates the
/// (Euclidian) length of a vector of numbers (either real or complex):
/// @code
/// template <typename T> real_t<T> l2_norm(const std::vector<T>& vec);
/// @endcode
/// The l2_norm is always real by definition.
/// (See https://en.wikipedia.org/wiki/Norm_(mathematics)#Euclidean_norm)
/// The return type of this function ("real_t<T>") indicates that
/// it returns a real number, even if the entries (of type T)
/// are complex numbers. In other words, by default, real_t<T> returns T.
/// However real_t<std::complex<T>> returns T (not std::complex<T>).
/// We define "real_t<T>" below using C++ template specializations:
/// @brief "realTypeMap" struct is used to the define "real_t<T>" type mapper
/// which returns the C++ type corresponding to the real component of T.
/// @details Consider a function ("l2_norm()") that calculates the
/// (Euclidian) length of a vector of numbers (either real or complex):
/// @code
/// template <typename T> real_t<T> l2_norm(const std::vector<T>& vec);
/// @endcode
/// The l2_norm is always real by definition.
/// (See https://en.wikipedia.org/wiki/Norm_(mathematics)#Euclidean_norm)
/// The return type of this function ("real_t<T>") indicates that
/// it returns a real number, even if the entries (of type T)
/// are complex numbers. In other words, by default, real_t<T> returns T.
/// However real_t<std::complex<T>> returns T (not std::complex<T>).
/// We define "real_t<T>" below using C++ template specializations:
template <typename T>
struct realTypeMap {
template <typename T>
struct realTypeMap {
typedef T type;
};
template <typename T>
struct realTypeMap<std::complex<T>> {
};
template <typename T>
struct realTypeMap<std::complex<T>> {
typedef T type;
};
template <typename T>
using real_t = typename realTypeMap<T>::type;
};
template <typename T>
using real_t = typename realTypeMap<T>::type;
// --- Operations on vectors (of real and complex numbers) ---
// --- Operations on vectors (of real and complex numbers) ---
/// @brief Calculate the inner product of two vectors.
/// (For vectors of complex numbers, std::conj() is used.)
template <typename T>
T inner_prod(const std::vector<T>& v1, const std::vector<T>& v2);
/// @brief Calculate the inner product of two vectors.
/// (For vectors of complex numbers, std::conj() is used.)
template <typename T>
T inner_prod(const std::vector<T>& v1, const std::vector<T>& v2);
/// @brief Compute the sum of the absolute values of the entries in v
/// @returns a real number (of type real_t<T>).
template <typename T>
real_t<T> l1_norm(const std::vector<T>& v);
/// @brief Compute the sum of the absolute values of the entries in v
/// @returns a real number (of type real_t<T>).
template <typename T>
real_t<T> l1_norm(const std::vector<T>& v);
/// @brief Calculate the l2_norm (Euclidian length) of vector v.
/// @returns a real number (of type real_t<T>).
template <typename T>
real_t<T> l2_norm(const std::vector<T>& v);
/// @brief Calculate the l2_norm (Euclidian length) of vector v.
/// @returns a real number (of type real_t<T>).
template <typename T>
real_t<T> l2_norm(const std::vector<T>& v);
/// @brief Multiply a vector (v) by a scalar (c).
template <typename T1, typename T2>
void scalar_mul(T1 c, std::vector<T2>& v);
/// @brief Multiply a vector (v) by a scalar (c).
template <typename T1, typename T2>
void scalar_mul(T1 c, std::vector<T2>& v);
/// @brief Divide vector "v" in-place by it's length (l2_norm(v)).
template <typename T>
void normalize(std::vector<T>& v);
/// @brief Divide vector "v" in-place by it's length (l2_norm(v)).
template <typename T>
void normalize(std::vector<T>& v);
// ---- Eigendecomposition of small dense symmetric matrices ----
// ---- Eigendecomposition of small dense symmetric matrices ----
/// @class Jacobi
/// @brief Calculate the eigenvalues and eigevectors of a symmetric matrix
/// using the Jacobi eigenvalue algorithm. Code for the Jacobi class
/// (along with tests and benchmarks) is available free of copyright at
/// https://github.com/jewettaij/jacobi_pd
/// @note The "Vector" and "Matrix" type arguments can be any
/// C or C++ object that support indexing, including pointers or vectors.
/// @details
/// Usage example:
/// @code
///
/// int n = 5; // Matrix size
/// double **M; // A symmetric n x n matrix you want to diagonalize
/// double *evals; // Store the eigenvalues here.
/// double **evects; // Store the eigenvectors here.
/// // Allocate space for M, evals, and evects, and load contents of M (omitted)
///
/// // Now create an instance of Jacobi ("eigen_calc"). This will allocate space
/// // for storing intermediate calculations. Once created, it can be reused
/// // multiple times without paying the cost of allocating memory on the heap.
///
/// 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 vectors over C-style pointers, this works also:
/// // Jacobi<double, vector<double>&, vector<vector<double>>&> eigen_calc(n);
///
/// // Now, calculate the eigenvalues and eigenvectors of M
///
/// eigen_calc.Diagonalize(M, evals, evects);
///
/// @endcode
/// @class Jacobi
/// @brief Calculate the eigenvalues and eigevectors of a symmetric matrix
/// using the Jacobi eigenvalue algorithm. Code for the Jacobi class
/// (along with tests and benchmarks) is available free of copyright at
/// https://github.com/jewettaij/jacobi_pd
/// @note The "Vector" and "Matrix" type arguments can be any
/// C or C++ object that support indexing, including pointers or vectors.
/// @details
/// -- Example: --
///
/// int n = 5; // Matrix size
/// double **M; // A symmetric n x n matrix you want to diagonalize
/// double *evals; // Store the eigenvalues here.
/// double **evects; // Store the eigenvectors here.
/// // Allocate space for M, evals, and evects, and load contents of M (omitted)
///
/// // Now create an instance of Jacobi ("eigen_calc"). This will allocate space
/// // for storing intermediate calculations. Once created, it can be reused
/// // multiple times without paying the cost of allocating memory on the heap.
///
/// 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 vectors over C-style pointers, this works also:
/// // Jacobi<double, vector<double>&, vector<vector<double>>&> eigen_calc(n);
///
/// // Now, calculate the eigenvalues and eigenvectors of M
///
/// eigen_calc.Diagonalize(M, evals, evects);
///
/// --- end of example ---
template<typename Scalar,
template<typename Scalar,
typename Vector,
typename Matrix,
typename ConstMatrix=Matrix>
class Jacobi
{
class Jacobi
{
int n; //!< the size of the matrices you want to diagonalize
Scalar **M; //!< local copy of the current matrix being analyzed
// Precomputed cosine, sine, and tangent of the most recent rotation angle:
@ -171,16 +171,19 @@ class Jacobi
Scalar t; //!< = tan(θ), (note |t|<=1)
int *max_idx_row; //!< = keep track of the the maximum element in row i (>i)
public:
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.
/// @param workspace optional array for storing intermediate calculations
/// (The vast majority of users can ignore this argument.)
Jacobi(int n = 0, Scalar **workspace=nullptr);
~Jacobi();
/// @brief Change the size of the matrices you want to diagonalize.
/// @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,
@ -192,8 +195,10 @@ public:
/// @brief Calculate the eigenvalues and eigevectors of a symmetric matrix
/// using the Jacobi eigenvalue algorithm.
/// @returns The number_of_sweeps (= number_of_iterations / (n*(n-1)/2)).
/// If this equals max_num_sweeps, the algorithm failed to converge.
/// @returns 0 if the algorithm converged,
/// 1 if the algorithm failed to converge. (IE, if the number of
/// pivot iterations exceeded max_num_sweeps * iters_per_sweep,
/// where iters_per_sweep = (n*(n-1)/2))
/// @note To reduce the computation time further, set calc_evecs=false.
int
Diagonalize(ConstMatrix mat, //!< the matrix you wish to diagonalize (size n)
@ -201,9 +206,11 @@ public:
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
int max_num_sweeps=50 //!< limit the number of iterations
);
private:
private:
bool is_M_preallocated;
// (Descriptions of private functions can be found in their implementation.)
void CalcRot(Scalar const *const *M, int i, int j);
void ApplyRot(Scalar **M, int i, int j);
@ -214,134 +221,182 @@ private:
void Init();
void Alloc(int n);
void Dealloc();
public:
public:
// C++ boilerplate: 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);
Jacobi<Scalar, Vector, Matrix, ConstMatrix>& operator = (Jacobi<Scalar, Vector, Matrix, ConstMatrix> source);
}; // class Jacobi
}; // class Jacobi
/// @brief
/// A specialized function which finds the eigenvalues and eigenvectors
/// of a 3x3 matrix (in double ** format).
/// @param mat the 3x3 matrix you wish to diagonalize
/// @param eval store the eigenvalues here
/// @param evec store the eigenvectors here...
/// @param evec_as_columns ...in the columns or in rows?
/// @param sort_criteria sort the resulting eigenvalues?
/// @note
/// When invoked using default arguments, this function is a stand-in for
/// the the version of "jacobi()" that was defined in "POEMS/fix_poems.cpp".
int jacobi3(double const **mat,
double *eval,
double **evec,
// optional arguments:
bool evec_as_columns = true,
Jacobi<double,double*,double**,double const*const*>::
SortCriteria sort_criteria =
Jacobi<double,double*,double**,double const*const*>::
SORT_DECREASING_EVALS
);
/// @brief
/// This version of jacobi3() finds the eigenvalues and eigenvectors
/// of a 3x3 matrix (which is in double (*)[3] format, instead of double **)
/// @param mat the 3x3 matrix you wish to diagonalize
/// @param eval store the eigenvalues here
/// @param evec store the eigenvectors here...
/// @param evec_as_columns ...in the columns or in rows?
/// @param sort_criteria sort the resulting eigenvalues?
/// @note
/// When invoked using default arguments, this function is a stand-in for
/// the previous version of "jacobi()" that was declared in "math_extra.h".
int jacobi3(double const mat[3][3],//!< the 3x3 matrix you wish to diagonalize
double *eval, //!< store the eigenvalues here
double evec[3][3], //!< store the eigenvectors here...
bool evec_as_columns=true, //!< ...as columns or as rows?
// optional arguments:
Jacobi<double,double*,double(*)[3],double const(*)[3]>::
SortCriteria sort_criteria=
Jacobi<double,double*,double(*)[3],double const(*)[3]>::
SORT_DECREASING_EVALS //!<sort results?
);
// ---- Eigendecomposition of large sparse (and dense) matrices ----
// ---- Eigendecomposition of large sparse (and dense) matrices ----
// The "LambdaLanczos" is a class useful for calculating eigenvalues
// and eigenvectors of large sparse matrices. Unfortunately, before the
// LambdaLanczos class can be declared, several additional expressions,
// classes and functions that it depends on must be declared first.
// The "LambdaLanczos" is a class useful for calculating eigenvalues
// and eigenvectors of large sparse matrices. Unfortunately, before the
// LambdaLanczos class can be declared, several additional expressions,
// classes and functions that it depends on must be declared first.
// @brief Create random vectors used at the beginning of the Lanczos algorithm.
// @note "Partially specialization of function" is not allowed, so
// it is mimicked by wrapping the "init" function with a class template.
template <typename T>
struct VectorRandomInitializer {
public:
// @brief Create random vectors used at the beginning of the Lanczos algorithm.
// @note "Partially specialization of function" is not allowed, so
// it is mimicked by wrapping the "init" function with a class template.
template <typename T>
struct VectorRandomInitializer {
public:
static void init(std::vector<T>&);
};
};
template <typename T>
struct VectorRandomInitializer<std::complex<T>> {
public:
template <typename T>
struct VectorRandomInitializer<std::complex<T>> {
public:
static void init(std::vector<std::complex<T>>&);
};
};
/// @brief Return the number of significant decimal digits of type T.
template <typename T>
inline constexpr int sig_decimal_digit() {
/// @brief Return the number of significant decimal digits of type T.
template <typename T>
inline constexpr int sig_decimal_digit() {
return (int)(std::numeric_limits<T>::digits *
std::log10(std::numeric_limits<T>::radix));
}
}
/// @brief Return 10^-n where n=number of significant decimal digits of type T.
template <typename T>
inline constexpr T minimum_effective_decimal() {
/// @brief Return 10^-n where n=number of significant decimal digits of type T.
template <typename T>
inline constexpr T minimum_effective_decimal() {
return std::pow(10, -sig_decimal_digit<T>());
}
}
/// @brief The LambdaLanczos class provides a general way to calculate
/// the smallest or largest eigenvalue and the corresponding eigenvector
/// of a symmetric (Hermitian) matrix using the Lanczos algorithm.
/// The characteristic feature of LambdaLanczos is that the matrix-vector
/// multiplication routine used in the Lanczos algorithm is adaptable.
/// @details
/// @code
///
/// //Example:
/// const int n = 3;
/// double M[n][n] = { {-1.0, -1.0, 1.0},
/// {-1.0, 1.0, 1.0},
/// { 1.0, 1.0, 1.0} };
/// // (Its eigenvalues are {-2, 1, 2})
///
/// // Specify the matrix-vector multiplication function
/// auto mv_mul = [&](const vector<double>& in, vector<double>& out) {
/// for(int i = 0;i < n;i++) {
/// for(int j = 0;j < n;j++) {
/// out[i] += M[i][j]*in[j];
/// }
/// }
/// };
///
/// LambdaLanczos<double> engine(mv_mul, n, true);
/// // ("true" means to calculate the largest eigenvalue.)
/// engine.eigenvalue_offset = 3.0 # = max_i{Σ_j|Mij|} (see below)
/// double eigenvalue;
/// vector<double> eigenvector(n);
/// int itern = engine.run(eigenvalue, eigenvector);
///
/// cout << "Iteration count: " << itern << endl;
/// cout << "Eigen value: " << setprecision(16) << eigenvalue << endl;
/// cout << "Eigen vector:";
/// for(int i = 0; i < n; i++) {
/// cout << eigenvector[i] << " ";
/// }
/// cout << endl;
///
/// @endcode
/// This feature allows you to use a matrix whose elements are partially given,
/// e.g. a sparse matrix whose non-zero elements are stored as a list of
/// {row-index, column-index, value} tuples. You can also easily combine
/// LambdaLanczos with existing matrix libraries (e.g. Eigen)
///
/// @note
/// If the matrices you want to analyze are ordinary square matrices, (as in
/// the example) it might be easier to use "PEigenDense" instead. (It is a
/// wrapper which takes care of all of the LambdaLanczos details for you.)
///
/// @note
/// IMPORTANT:
/// The Lanczos algorithm finds the largest magnitude eigenvalue, so you
/// MUST ensure that the eigenvalue you are seeking has the largest magnitude
/// (regardless of whether it is the maximum or minimum eigenvalue).
/// To insure that this is so, you can add or subtract a number to all
/// of the eigenvalues of the matrix by specifying the "eigenvalue_offset".
/// This number should exceed the largest magnitude eigenvalue of the matrix.
/// According to the Gershgorin theorem, you can estimate this number using
/// r = max_i{Σ_j|Mij|} = max_j{Σ_i|Mij|}
/// (where Mij are the elements of the matrix and Σ_j denotes the sum over j).
/// If find_maximum == true (if you are seeking the maximum eigenvalue), then
/// eigenvalue_offset = +r
/// If find_maximum == false, then
/// eigenvalue_offset = -r
/// The eigenvalue_offset MUST be specified by the user. LambdaLanczos does
/// not have an efficient and general way to access the elements of the matrix.
///
/// (You can omit this step if you are seeking the maximum eigenvalue,
/// and the matrix is positive definite, or if you are seeking the minimum
/// eigenvalue and the matrix is negative definite.)
///
/// @note
/// LambdaLanczos is available under the MIT license and downloadable at:
/// https://github.com/mrcdr/lambda-lanczos
/// @class LambdaLanczos
/// @brief The LambdaLanczos class provides a general way to calculate
/// the smallest or largest eigenvalue and the corresponding eigenvector
/// of a symmetric (Hermitian) matrix using the Lanczos algorithm.
/// The characteristic feature of LambdaLanczos is that the matrix-vector
/// multiplication routine used in the Lanczos algorithm is adaptable.
/// @details
/// @code
///
/// //Example:
/// const int n = 3;
/// double M[n][n] = { {-1.0, -1.0, 1.0},
/// {-1.0, 1.0, 1.0},
/// { 1.0, 1.0, 1.0} };
/// // (Its eigenvalues are {-2, 1, 2})
///
/// // Specify the matrix-vector multiplication function
/// auto mv_mul = [&](const vector<double>& in, vector<double>& out) {
/// for(int i = 0;i < n;i++) {
/// for(int j = 0;j < n;j++) {
/// out[i] += M[i][j]*in[j];
/// }
/// }
/// };
///
/// LambdaLanczos<double> engine(mv_mul, n, true);
/// // ("true" means to calculate the largest eigenvalue.)
/// engine.eigenvalue_offset = 3.0; // = max_i{Σ_j|Mij|} (see below)
///
/// double eigenvalue; //(must never be a complex number even if M is complex)
/// vector<double> eigenvector(n);
///
/// int itern = engine.run(eigenvalue, eigenvector);
///
/// cout << "Iteration count: " << itern << endl;
/// cout << "Eigen value: " << setprecision(16) << eigenvalue << endl;
/// cout << "Eigen vector:";
/// for(int i = 0; i < n; i++) {
/// cout << eigenvector[i] << " ";
/// }
/// cout << endl;
///
/// @endcode
/// This feature allows you to use a matrix whose elements are partially given,
/// e.g. a sparse matrix whose non-zero elements are stored as a list of
/// {row-index, column-index, value} tuples. You can also easily combine
/// LambdaLanczos with existing matrix libraries (e.g. Eigen)
///
/// @note
/// If the matrices you want to analyze are ordinary square matrices, (as in
/// the example) it might be easier to use "PEigenDense" instead. (It is a
/// wrapper which takes care of all of the LambdaLanczos details for you.)
///
/// @note
/// IMPORTANT:
/// The Lanczos algorithm finds the largest magnitude eigenvalue, so you
/// MUST ensure that the eigenvalue you are seeking has the largest magnitude
/// (regardless of whether it is the maximum or minimum eigenvalue).
/// To insure that this is so, you can add or subtract a number to all
/// of the eigenvalues of the matrix by specifying the "eigenvalue_offset".
/// This number should exceed the largest magnitude eigenvalue of the matrix.
/// According to the Gershgorin theorem, you can estimate this number using
/// r = max_i{Σ_j|Mij|} = max_j{Σ_i|Mij|}
/// (where Mij are the elements of the matrix and Σ_j denotes the sum over j).
/// If find_maximum == true (if you are seeking the maximum eigenvalue), then
/// eigenvalue_offset = +r
/// If find_maximum == false, then
/// eigenvalue_offset = -r
/// The eigenvalue_offset MUST be specified by the user. LambdaLanczos does
/// not have an efficient and general way to access the elements of the matrix.
///
/// (You can omit this step if you are seeking the maximum eigenvalue,
/// and the matrix is positive definite, or if you are seeking the minimum
/// eigenvalue and the matrix is negative definite.)
///
/// @note
/// LambdaLanczos is available under the MIT license and downloadable at:
/// https://github.com/mrcdr/lambda-lanczos
template <typename T>
class LambdaLanczos {
public:
template <typename T>
class LambdaLanczos {
public:
LambdaLanczos();
LambdaLanczos(std::function<void(const std::vector<T>&, std::vector<T>&)> mv_mul, int matrix_size, bool find_maximum);
LambdaLanczos(std::function<void(const std::vector<T>&, std::vector<T>&)> mv_mul, int matrix_size) : LambdaLanczos(mv_mul, matrix_size, true) {}
@ -400,7 +455,7 @@ public:
void SetEpsilon(T eps);
void SetTriEpsRatio(T tridiag_eps_ratio);
private:
private:
static void schmidt_orth(std::vector<T>&, const std::vector<std::vector<T>>&);
real_t<T> find_minimum_eigenvalue(const std::vector<real_t<T>>&,
const std::vector<real_t<T>>&) const;
@ -412,37 +467,47 @@ private:
const std::vector<real_t<T>>&,
const std::vector<real_t<T>>&);
real_t<T> UpperBoundEvals() const;
};
};
/// @brief
/// PEigenDense is a class containing only one useful member function:
/// PrincipalEigen(). This function calculates the principal (largest
/// or smallest) eigenvalue and corresponding eigenvector of a square
/// n x n matrix. This can be faster than diagionalizing the entire matrix.
/// (For example by using the Lanczos algorithm or something similar.)
/// @note
/// This code is a wrapper. Internally, it uses the "LambdaLanczos" class.
/// @note
/// For matrices larger than 13x13, PEigenDense::PrincipleEigen()
/// is usually faster than Jacobi::Diagonalize().)
/// @class PEigenDense
/// @brief
/// PEigenDense is a class containing only one useful member function:
/// PrincipalEigen(). This function calculates the principal (largest
/// or smallest) eigenvalue and corresponding eigenvector of a square
/// n x n matrix. This can be faster than diagionalizing the entire matrix.
/// (For example by using the Lanczos algorithm or something similar.)
/// @note
/// This code is a wrapper. Internally, it uses the "LambdaLanczos" class.
/// @note
/// For dense matrices smaller than 13x13, Jacobi::Diagonalize(),
/// is usually faster than PEigenDense::PrincipleEigen().
/// @details
/// Usage example:
/// @code
///
/// const int n = 100;
///
/// PEigenDense<double, double*, double**, double const*const*> pe(n);
///
/// double **M; // A symmetric n x n matrix you want to diagonalize
/// double evect[n]; // Store the principal eigenvector here.
///
/// // Now, allocate space for M and load it's contents. (omitted)
///
/// double eval = pe.PrincipalEigen(M, evect, true); //Returns the max eval
///
/// @endcode
template<typename Scalar, typename Vector, typename ConstMatrix>
class PEigenDense
{
template<typename Scalar, typename Vector, typename ConstMatrix>
class PEigenDense
{
size_t n; // the size of the matrix
std::vector<Scalar> evec; // preallocated vector
public:
void SetSize(int matrix_size) {
n = matrix_size;
evec.resize(n);
}
PEigenDense(int matrix_size=0):evec(matrix_size) {
SetSize(matrix_size);
}
public:
PEigenDense(int matrix_size=0);
/// @brief Calculate the principal eigenvalue and eigenvector of a matrix.
/// @return Return the principal eigenvalue of the matrix.
@ -452,7 +517,9 @@ public:
Vector evector, //!< the eigenvector is stored here
bool find_max=false); //!< want the max or min eigenvalue?
}; // class PEigenDense
void SetSize(int matrix_size); //change matrix size after instantiation
}; // class PEigenDense
@ -460,16 +527,13 @@ public:
// ----------- IMPLEMENTATION -----------
// --------------------------------------
// --- Implementation: Memory allocation for matrices ---
template<typename Entry>
void Alloc2D(size_t nrows, // size of the array (number of rows)
size_t ncols, // size of the array (number of columns)
Entry ***paaX) // pointer to a 2D C-style array
{
assert(paaX);
//assert(paaX);
*paaX = new Entry* [nrows]; //conventional 2D C array (pointer-to-pointer)
(*paaX)[0] = new Entry [nrows * ncols]; // 1D C array (contiguous memor)
for(size_t iy=0; iy<nrows; iy++)
@ -549,6 +613,28 @@ inline real_t<T> l1_norm(const std::vector<T>& vec) {
// --- Implementation: Eigendecomposition of small dense matrices ---
template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>
Jacobi<Scalar, Vector, Matrix, ConstMatrix>::
Jacobi(int n, Scalar **workspace) {
Init();
if (workspace) {
is_M_preallocated = true;
M = workspace;
this->n = n;
}
else {
is_M_preallocated = false;
SetSize(n);
}
}
template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>
Jacobi<Scalar, Vector, Matrix, ConstMatrix>::
~Jacobi() {
if (! is_M_preallocated)
Dealloc();
}
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)
@ -602,7 +688,7 @@ Diagonalize(ConstMatrix mat, // the matrix you wish to diagonalize (size n)
// Optional: Sort results by eigenvalue.
SortRows(eval, evec, n, sort_criteria);
return n_iters / (n*(n-1)/2); //returns the number of "sweeps" (converged?)
return (n_iters == max_num_iters);
}
@ -634,7 +720,7 @@ CalcRot(Scalar const *const *M, //!< matrix
t = -t;
}
}
assert(std::abs(t) <= 1.0);
//assert(std::abs(t) <= 1.0);
c = 1.0 / std::sqrt(1 + t*t);
s = c*t;
}
@ -699,7 +785,7 @@ ApplyRot(Scalar **M, // matrix
M[j][j] += t * M[i][j];
//Update the off-diagonal elements of M which will change (above the diagonal)
assert(i < j);
//assert(i < j);
M[i][j] = 0.0;
//compute M[w][i] and M[i][w] for all w!=i,considering above-diagonal elements
@ -887,7 +973,7 @@ Jacobi(const Jacobi<Scalar, Vector, Matrix, ConstMatrix>& source)
{
Init();
SetSize(source.n);
assert(n == source.n);
//assert(n == source.n);
// The following lines aren't really necessary, because the contents
// of source.M and source.max_idx_row are not needed (since they are
// overwritten every time Jacobi::Diagonalize() is invoked).
@ -913,7 +999,7 @@ template<typename Scalar,typename Vector,typename Matrix,typename ConstMatrix>
Jacobi<Scalar, Vector, Matrix, ConstMatrix>::
Jacobi(Jacobi<Scalar, Vector, Matrix, ConstMatrix>&& other) {
Init();
swap(*this, other);
this->swap(other);
}
// Using the "copy-swap" idiom for the assignment operator
@ -955,8 +1041,8 @@ template <typename T>
inline int LambdaLanczos<T>::
run(real_t<T>& eigvalue, std::vector<T>& eigvec) const
{
assert(matrix_size > 0);
assert(0 < this->tridiag_eps_ratio && this->tridiag_eps_ratio < 1);
//assert(matrix_size > 0);
//assert(0 < this->tridiag_eps_ratio && this->tridiag_eps_ratio < 1);
std::vector<std::vector<T>> u; // Lanczos vectors
std::vector<real_t<T>> alpha; // Diagonal elements of an approximated tridiagonal matrix
@ -1321,13 +1407,26 @@ init(std::vector<std::complex<T>>& v)
// --- Implementation of PEigenDense
template<typename Scalar, typename Vector, typename ConstMatrix>
void PEigenDense<Scalar, Vector, ConstMatrix>::
SetSize(int matrix_size) {
n = matrix_size;
evec.resize(n);
}
template<typename Scalar, typename Vector, typename ConstMatrix>
PEigenDense<Scalar, Vector, ConstMatrix>::
PEigenDense(int matrix_size):evec(matrix_size) {
SetSize(matrix_size);
}
template<typename Scalar, typename Vector, typename ConstMatrix>
Scalar PEigenDense<Scalar, Vector, ConstMatrix>::
PrincipalEigen(ConstMatrix matrix,
Vector eigenvector,
bool find_max)
{
assert(n > 0);
//assert(n > 0);
auto matmul = [&](const std::vector<Scalar>& in, std::vector<Scalar>& out) {
for(int i = 0; i < n; i++) {
for(int j = 0; j < n; j++) {
@ -1373,8 +1472,12 @@ PrincipalEigen(ConstMatrix matrix,
return eval;
}
} //namespace MathEigen
#endif //#ifndef _MATH_EIGEN_H
#endif //#ifndef LMP_MATH_EIGEN_H

84
src/math_extra.cpp
View File

@ -19,8 +19,6 @@
#include <cstdio>
#include <cstring>
#define MAXJACOBI 50
namespace MathExtra {
/* ----------------------------------------------------------------------
@ -92,88 +90,6 @@ int mldivide3(const double m[3][3], const double *v, double *ans)
return 0;
}
/* ----------------------------------------------------------------------
compute evalues and evectors of 3x3 real symmetric matrix
based on Jacobi rotations
adapted from Numerical Recipes jacobi() function
------------------------------------------------------------------------- */
int jacobi(double matrix[3][3], double *evalues, double evectors[3][3])
{
int i,j,k;
double tresh,theta,tau,t,sm,s,h,g,c,b[3],z[3];
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) evectors[i][j] = 0.0;
evectors[i][i] = 1.0;
}
for (i = 0; i < 3; i++) {
b[i] = evalues[i] = matrix[i][i];
z[i] = 0.0;
}
for (int iter = 1; iter <= MAXJACOBI; iter++) {
sm = 0.0;
for (i = 0; i < 2; i++)
for (j = i+1; j < 3; j++)
sm += fabs(matrix[i][j]);
if (sm == 0.0) return 0;
if (iter < 4) tresh = 0.2*sm/(3*3);
else tresh = 0.0;
for (i = 0; i < 2; i++) {
for (j = i+1; j < 3; j++) {
g = 100.0*fabs(matrix[i][j]);
if (iter > 4 && fabs(evalues[i])+g == fabs(evalues[i])
&& fabs(evalues[j])+g == fabs(evalues[j]))
matrix[i][j] = 0.0;
else if (fabs(matrix[i][j]) > tresh) {
h = evalues[j]-evalues[i];
if (fabs(h)+g == fabs(h)) t = (matrix[i][j])/h;
else {
theta = 0.5*h/(matrix[i][j]);
t = 1.0/(fabs(theta)+sqrt(1.0+theta*theta));
if (theta < 0.0) t = -t;
}
c = 1.0/sqrt(1.0+t*t);
s = t*c;
tau = s/(1.0+c);
h = t*matrix[i][j];
z[i] -= h;
z[j] += h;
evalues[i] -= h;
evalues[j] += h;
matrix[i][j] = 0.0;
for (k = 0; k < i; k++) rotate(matrix,k,i,k,j,s,tau);
for (k = i+1; k < j; k++) rotate(matrix,i,k,k,j,s,tau);
for (k = j+1; k < 3; k++) rotate(matrix,i,k,j,k,s,tau);
for (k = 0; k < 3; k++) rotate(evectors,k,i,k,j,s,tau);
}
}
}
for (i = 0; i < 3; i++) {
evalues[i] = b[i] += z[i];
z[i] = 0.0;
}
}
return 1;
}
/* ----------------------------------------------------------------------
perform a single Jacobi rotation
------------------------------------------------------------------------- */
void rotate(double matrix[3][3], int i, int j, int k, int l,
double s, double tau)
{
double g = matrix[i][j];
double h = matrix[k][l];
matrix[i][j] = g-s*(h+g*tau);
matrix[k][l] = h+s*(g-h*tau);
}
/* ----------------------------------------------------------------------
Richardson iteration to update quaternion from angular momentum
return new normalized quaternion q

1
src/math_extra.h
View File

@ -85,7 +85,6 @@ namespace MathExtra {
void write3(const double mat[3][3]);
int mldivide3(const double mat[3][3], const double *vec, double *ans);
int jacobi(double matrix[3][3], double *evalues, double evectors[3][3]);
void rotate(double matrix[3][3], int i, int j, int k, int l,
double s, double tau);
void richardson(double *q, double *m, double *w, double *moments, double dtq);

3
src/molecule.cpp
View File

@ -23,6 +23,7 @@
#include "comm.h"
#include "domain.h"
#include "math_extra.h"
#include "math_eigen.h"
#include "memory.h"
#include "error.h"
#include "utils.h"
@ -350,7 +351,7 @@ void Molecule::compute_inertia()
tensor[0][2] = tensor[2][0] = itensor[4];
tensor[0][1] = tensor[1][0] = itensor[5];
if (MathExtra::jacobi(tensor,inertia,evectors))
if (MathEigen::jacobi3(tensor,inertia,evectors))
error->all(FLERR,"Insufficient Jacobi rotations for rigid molecule");
ex[0] = evectors[0][0];