diff --git a/lib/colvars/colvartypes.cpp b/lib/colvars/colvartypes.cpp index da00183323..460eec06ac 100644 --- a/lib/colvars/colvartypes.cpp +++ b/lib/colvars/colvartypes.cpp @@ -10,6 +10,7 @@ #include #include +#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 &m, - cvm::vector1d &eigval, - cvm::matrix2d &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 b(n); - cvm::vector1d z(n); - - for (ip=0;ip 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= p) p=d[k=j]; - if (k != i) { - d[k]=d[i]; - d[i]=p; - for (j=0;j +#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 &m, - cvm::vector1d &eigval, - cvm::matrix2d &eigvec); + MathEigen::Jacobi &, + cvm::matrix2d &> ecalc; }; diff --git a/lib/jacobi_pd.h b/lib/jacobi_pd.h new file mode 100644 index 0000000000..ca68e6ba3a --- /dev/null +++ b/lib/jacobi_pd.h @@ -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 +#include +//#include + + +namespace MathEigen { + +/// @brief Allocate a 2-dimensional array. (Uses row-major order.) +/// @note (Matrix allocation functions can also be found in colvartypes.h) +template +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 +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 eigen_calc(n); +/// +/// // Note: +/// // If the matrix you plan to diagonalize (M) is read-only, use this instead: +/// // Jacobi eigen_calc(n); +/// // If you prefer using vectors over C-style pointers, this works also: +/// // Jacobi&, vector>&> eigen_calc(n); +/// +/// // Now, calculate the eigenvalues and eigenvectors of M +/// +/// eigen_calc.Diagonalize(M, evals, evects); +/// +/// @endcode + +template +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,//!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& source); + Jacobi(Jacobi&& other); + void swap(Jacobi &other); + Jacobi& operator = (Jacobi source); + +}; // class Jacobi + + + + + +// -------------- IMPLEMENTATION -------------- + +// Utilities for memory allocation + +template +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 +void Dealloc2D(Entry ***paaX) // pointer to a 2D C-style array +{ + if (paaX && *paaX) { + delete [] (*paaX)[0]; + delete [] (*paaX); + *paaX = nullptr; + } +} + + + +template +int Jacobi:: +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 +void Jacobi:: +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 iv) 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 +void Jacobi:: +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 +void Jacobi:: +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 +int Jacobi:: +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 +void Jacobi:: +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 +void Jacobi:: +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 +void Jacobi:: +Init() { + n = 0; + M = nullptr; + max_idx_row = nullptr; +} + +template +void Jacobi:: +SetSize(int n) { + Dealloc(); + Alloc(n); +} + +// memory management: + +template +void Jacobi:: +Alloc(int n) { + this->n = n; + if (n > 0) { + max_idx_row = new int[n]; + Alloc2D(n, n, &M); + } +} + +template +void Jacobi:: +Dealloc() { + if (max_idx_row) + delete [] max_idx_row; + Dealloc2D(&M); + Init(); +} + +// memory management: copy and move constructor, swap, and assignment operator: + +template +Jacobi:: +Jacobi(const Jacobi& 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 +void Jacobi:: +swap(Jacobi &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 +Jacobi:: +Jacobi(Jacobi&& other) { + Init(); + this->swap(other); +} + +// Using the "copy-swap" idiom for the assignment operator +template +Jacobi& +Jacobi:: +operator = (Jacobi source) { + this->swap(source); + return *this; +} + +} // namespace MathEigen + + +#endif //#ifndef _MATH_EIGEN_H + diff --git a/src/BODY/body_nparticle.cpp b/src/BODY/body_nparticle.cpp index a99c025e1f..52ae2a2eb8 100644 --- a/src/BODY/body_nparticle.cpp +++ b/src/BODY/body_nparticle.cpp @@ -16,6 +16,7 @@ #include #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"); diff --git a/src/BODY/body_rounded_polygon.cpp b/src/BODY/body_rounded_polygon.cpp index 349ad957c2..375ec7c5d7 100644 --- a/src/BODY/body_rounded_polygon.cpp +++ b/src/BODY/body_rounded_polygon.cpp @@ -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"); diff --git a/src/BODY/body_rounded_polyhedron.cpp b/src/BODY/body_rounded_polyhedron.cpp index cb84517b6f..73a6148b3e 100644 --- a/src/BODY/body_rounded_polyhedron.cpp +++ b/src/BODY/body_rounded_polyhedron.cpp @@ -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"); diff --git a/src/POEMS/fix_poems.cpp b/src/POEMS/fix_poems.cpp index b0b1287717..aeca600a51 100644 --- a/src/POEMS/fix_poems.cpp +++ b/src/POEMS/fix_poems.cpp @@ -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 diff --git a/src/POEMS/fix_poems.h b/src/POEMS/fix_poems.h index 6892c51d7e..7c81889639 100644 --- a/src/POEMS/fix_poems.h +++ b/src/POEMS/fix_poems.h @@ -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(); diff --git a/src/RIGID/fix_rigid.cpp b/src/RIGID/fix_rigid.cpp index 761b957292..293d86f7a9 100644 --- a/src/RIGID/fix_rigid.cpp +++ b/src/RIGID/fix_rigid.cpp @@ -17,6 +17,7 @@ #include #include #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"); diff --git a/src/RIGID/fix_rigid_small.cpp b/src/RIGID/fix_rigid_small.cpp index 01dd9e889e..3af63dc5ec 100644 --- a/src/RIGID/fix_rigid_small.cpp +++ b/src/RIGID/fix_rigid_small.cpp @@ -18,6 +18,7 @@ #include #include #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"); diff --git a/src/USER-MISC/compute_gyration_shape.cpp b/src/USER-MISC/compute_gyration_shape.cpp index aef5ef91a3..7b5d77ea04 100644 --- a/src/USER-MISC/compute_gyration_shape.cpp +++ b/src/USER-MISC/compute_gyration_shape.cpp @@ -21,6 +21,7 @@ #include #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"); diff --git a/src/USER-MISC/compute_gyration_shape_chunk.cpp b/src/USER-MISC/compute_gyration_shape_chunk.cpp index 116d851024..e0fc5c1342 100644 --- a/src/USER-MISC/compute_gyration_shape_chunk.cpp +++ b/src/USER-MISC/compute_gyration_shape_chunk.cpp @@ -21,6 +21,7 @@ #include #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"); diff --git a/src/atom_vec_tri.cpp b/src/atom_vec_tri.cpp index d5b48b9540..a3409ef584 100644 --- a/src/atom_vec_tri.cpp +++ b/src/atom_vec_tri.cpp @@ -15,6 +15,7 @@ #include #include #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]; diff --git a/src/compute_omega_chunk.cpp b/src/compute_omega_chunk.cpp index 327c64493f..5c51d036ac 100644 --- a/src/compute_omega_chunk.cpp +++ b/src/compute_omega_chunk.cpp @@ -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"); diff --git a/src/group.cpp b/src/group.cpp index 2377f028bd..c007ba813e 100644 --- a/src/group.cpp +++ b/src/group.cpp @@ -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"); diff --git a/src/math_eigen.cpp b/src/math_eigen.cpp new file mode 100644 index 0000000000..b86b626309 --- /dev/null +++ b/src/math_eigen.cpp @@ -0,0 +1,140 @@ +#include "math_eigen.h" +#include + +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; +// +// ...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; + +// template instantiations of Jacbi for fixed-size (3x3) arrays +template class MathEigen::Jacobi; + +// template instantiations of Jacobi for vector-of-vectors: +template class MathEigen::Jacobi&, + vector >&, + const vector >& >; + +// template instantiations of Jacobi for array-of-arrays: +template class MathEigen::Jacobi&, + array, 3 >&, + const array, 3 >& >; +// template instantiations of LambdaLanczos +template class MathEigen::LambdaLanczos; + +// template instantiations of PEidenDense for pointer->pointer arrays +template class MathEigen::PEigenDense; + +// template instantiations of PEidenDense for vector> +template class MathEigen::PEigenDense&, + const vector >&>; + +// ----------------- now create instances for floats ---------------- + +template class MathEigen::Jacobi; +template class MathEigen::Jacobi; +template class MathEigen::Jacobi&, + vector >&, + const vector >& >; +template class MathEigen::Jacobi&, + array, 3 >&, + const array, 3 >& >; +template class MathEigen::LambdaLanczos; +template class MathEigen::PEigenDense; +template class MathEigen::PEigenDense&, + const vector >&>; + +// 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::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 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::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 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; +} + diff --git a/src/math_eigen.h b/src/math_eigen.h index f8d2be4e21..3a697584d7 100644 --- a/src/math_eigen.h +++ b/src/math_eigen.h @@ -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 +//#include #include #include #include @@ -51,408 +51,475 @@ 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 -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 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 + 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 -void Dealloc2D(Entry ***paaX); //!< pointer to 2D multidimensional array + /// @brief Deallocate arrays that were created using Alloc2D(). + template + void Dealloc2D(Entry ***paaX); //!< pointer to 2D multidimensional array -// --- Complex numbers --- + // --- Complex numbers --- -/// @brief "realTypeMap" struct is used to the define "real_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 real_t l2_norm(const std::vector& 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") indicates that -/// it returns a real number, even if the entries (of type T) -/// are complex numbers. In other words, by default, real_t returns T. -/// However real_t> returns T (not std::complex). -/// We define "real_t" below using C++ template specializations: + /// @brief "realTypeMap" struct is used to the define "real_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 real_t l2_norm(const std::vector& 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") indicates that + /// it returns a real number, even if the entries (of type T) + /// are complex numbers. In other words, by default, real_t returns T. + /// However real_t> returns T (not std::complex). + /// We define "real_t" below using C++ template specializations: -template -struct realTypeMap { - typedef T type; -}; -template -struct realTypeMap> { - typedef T type; -}; -template -using real_t = typename realTypeMap::type; + template + struct realTypeMap { + typedef T type; + }; + template + struct realTypeMap> { + typedef T type; + }; + template + using real_t = typename realTypeMap::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 -T inner_prod(const std::vector& v1, const std::vector& v2); + /// @brief Calculate the inner product of two vectors. + /// (For vectors of complex numbers, std::conj() is used.) + template + T inner_prod(const std::vector& v1, const std::vector& v2); -/// @brief Compute the sum of the absolute values of the entries in v -/// @returns a real number (of type real_t). -template -real_t l1_norm(const std::vector& v); + /// @brief Compute the sum of the absolute values of the entries in v + /// @returns a real number (of type real_t). + template + real_t l1_norm(const std::vector& v); -/// @brief Calculate the l2_norm (Euclidian length) of vector v. -/// @returns a real number (of type real_t). -template -real_t l2_norm(const std::vector& v); + /// @brief Calculate the l2_norm (Euclidian length) of vector v. + /// @returns a real number (of type real_t). + template + real_t l2_norm(const std::vector& v); -/// @brief Multiply a vector (v) by a scalar (c). -template -void scalar_mul(T1 c, std::vector& v); + /// @brief Multiply a vector (v) by a scalar (c). + template + void scalar_mul(T1 c, std::vector& v); -/// @brief Divide vector "v" in-place by it's length (l2_norm(v)). -template -void normalize(std::vector& v); + /// @brief Divide vector "v" in-place by it's length (l2_norm(v)). + template + void normalize(std::vector& 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 -/// -- 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 eigen_calc(n); -/// -/// // Note: -/// // If the matrix you plan to diagonalize (M) is read-only, use this instead: -/// // Jacobi eigen_calc(n); -/// // If you prefer using vectors over C-style pointers, this works also: -/// // Jacobi&, vector>&> eigen_calc(n); -/// -/// // Now, calculate the eigenvalues and eigenvectors of M -/// -/// eigen_calc.Diagonalize(M, evals, evects); -/// -/// --- end of example --- - -template - -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: - Scalar c; //!< = cos(θ) - Scalar s; //!< = sin(θ) - Scalar t; //!< = tan(θ), (note |t|<=1) - int *max_idx_row; //!< = keep track of the the maximum element in row i (>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; - + /// @class Jacobi /// @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. - /// @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,//! eigen_calc(n); + /// + /// // Note: + /// // If the matrix you plan to diagonalize (M) is read-only, use this instead: + /// // Jacobi eigen_calc(n); + /// // If you prefer using vectors over C-style pointers, this works also: + /// // Jacobi&, vector>&> eigen_calc(n); + /// + /// // Now, calculate the eigenvalues and eigenvectors of M + /// + /// eigen_calc.Diagonalize(M, evals, evects); + /// + /// @endcode -private: - // (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); - void ApplyRotLeft(Matrix E, int i, int j); - int MaxEntryRow(Scalar const *const *M, int i) const; - void MaxEntry(Scalar const *const *M, int& i_max, int& j_max) const; - void SortRows(Vector v, Matrix M, int n, SortCriteria s=SORT_DECREASING_EVALS) const; - void Init(); - void Alloc(int n); - void Dealloc(); -public: - // C++ boilerplate: copy and move constructor, swap, and assignment operator - Jacobi(const Jacobi& source); - Jacobi(Jacobi&& other); - void swap(Jacobi &other); - Jacobi& operator = (Jacobi source); + template -}; // 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: + Scalar c; //!< = cos(θ) + Scalar s; //!< = sin(θ) + Scalar t; //!< = tan(θ), (note |t|<=1) + int *max_idx_row; //!< = keep track of the the maximum element in row i (>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. + /// @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); + + // @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,//!& source); + Jacobi(Jacobi&& other); + void swap(Jacobi &other); + Jacobi& operator = (Jacobi source); + + }; // 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:: + SortCriteria sort_criteria = + Jacobi:: + 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:: + SortCriteria sort_criteria= + Jacobi:: + SORT_DECREASING_EVALS //! + struct VectorRandomInitializer { + public: + static void init(std::vector&); + }; -// @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 -struct VectorRandomInitializer { -public: - static void init(std::vector&); -}; + template + struct VectorRandomInitializer> { + public: + static void init(std::vector>&); + }; -template -struct VectorRandomInitializer> { -public: - static void init(std::vector>&); -}; - -/// @brief Return the number of significant decimal digits of type T. -template -inline constexpr int sig_decimal_digit() { - return (int)(std::numeric_limits::digits * - std::log10(std::numeric_limits::radix)); -} - -/// @brief Return 10^-n where n=number of significant decimal digits of type T. -template -inline constexpr T minimum_effective_decimal() { - return std::pow(10, -sig_decimal_digit()); -} - - -/// @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& in, vector& out) { -/// for(int i = 0;i < n;i++) { -/// for(int j = 0;j < n;j++) { -/// out[i] += M[i][j]*in[j]; -/// } -/// } -/// }; -/// -/// LambdaLanczos 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 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 -class LambdaLanczos { -public: - LambdaLanczos(); - LambdaLanczos(std::function&, std::vector&)> mv_mul, int matrix_size, bool find_maximum); - LambdaLanczos(std::function&, std::vector&)> mv_mul, int matrix_size) : LambdaLanczos(mv_mul, matrix_size, true) {} - - /// @brief Calculate the principal (largest or smallest) eigenvalue - /// of the matrix (and its corresponding eigenvector). - int run(real_t&, std::vector&) const; - - // --- public data members --- - - /// @brief Specify the size of the matrix you will analyze. - /// (This equals the size of the eigenvector which will be returned.) - int matrix_size; - - /// @brief Specify the function used for matrix*vector multiplication - /// used by the Lanczos algorithm. For an ordinary dense matrix, - /// this function is the ordinary matrix*vector product. (See the - /// example above. For a sparse matrix, it will be something else.) - std::function&, std::vector&)> mv_mul; - - /// @brief Are we searching for the maximum or minimum eigenvalue? - /// @note (Usually, you must also specify eigenvalue_offset.) - bool find_maximum = false; - - /// @brief Shift all the eigenvalues by "eigenvalue_offset" during the Lanczos - /// iteration (ie. during LambdaLanczos::run()). The goal is to insure - /// that the correct eigenvalue is selected (the one with the maximum - /// magnitude). - /// @note The eigevalue returned by LambdaLanczos::run() is not effected - /// because after the iteration is finished, it will subtract this - /// number from the eigenvalue before it is returned to the caller. - /// @note Unless your matrix is positive definite or negative definite, - /// you MUST specify eigenvalue_offset. See comment above for details. - real_t eigenvalue_offset = 0.0; - - /// @brief This function sets "eigenvalue_offset" automatically. - /// @note Using this function is not recommended because it is very slow. - /// For efficiency, set the "eigenvalue_offset" yourself. - void ChooseOffset(); - - // The remaining data members usually can be left alone: - int max_iteration; - real_t eps = minimum_effective_decimal>() * 1e3; - real_t tridiag_eps_ratio = 1e-1; - int initial_vector_size = 200; - std::function&)> init_vector = - VectorRandomInitializer::init; - - // (for those who prefer "Set" functions...) - int SetSize(int matrix_size); - void SetMul(std::function&, - std::vector&)> mv_mul); - void SetInitVec(std::function&)> init_vector); - void SetFindMax(bool find_maximum); - void SetEvalOffset(T eigenvalue_offset); - void SetEpsilon(T eps); - void SetTriEpsRatio(T tridiag_eps_ratio); - -private: - static void schmidt_orth(std::vector&, const std::vector>&); - real_t find_minimum_eigenvalue(const std::vector>&, - const std::vector>&) const; - real_t find_maximum_eigenvalue(const std::vector>&, - const std::vector>&) const; - static real_t tridiagonal_eigen_limit(const std::vector>&, - const std::vector>&); - static int num_of_eigs_smaller_than(real_t, - const std::vector>&, - const std::vector>&); - real_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().) - -template -class PEigenDense -{ - size_t n; // the size of the matrix - std::vector evec; // preallocated vector - -public: - void SetSize(int matrix_size) { - n = matrix_size; - evec.resize(n); + /// @brief Return the number of significant decimal digits of type T. + template + inline constexpr int sig_decimal_digit() { + return (int)(std::numeric_limits::digits * + std::log10(std::numeric_limits::radix)); } - PEigenDense(int matrix_size=0):evec(matrix_size) { - SetSize(matrix_size); + /// @brief Return 10^-n where n=number of significant decimal digits of type T. + template + inline constexpr T minimum_effective_decimal() { + return std::pow(10, -sig_decimal_digit()); } - /// @brief Calculate the principal eigenvalue and eigenvector of a matrix. - /// @return Return the principal eigenvalue of the matrix. - /// If you want the eigenvector, pass a non-null "evector" argument. - Scalar - PrincipalEigen(ConstMatrix matrix, //!< the input patrix - Vector evector, //!< the eigenvector is stored here - bool find_max=false); //!< want the max or min eigenvalue? -}; // class PEigenDense + /// @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& in, vector& out) { + /// for(int i = 0;i < n;i++) { + /// for(int j = 0;j < n;j++) { + /// out[i] += M[i][j]*in[j]; + /// } + /// } + /// }; + /// + /// LambdaLanczos 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 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 + class LambdaLanczos { + public: + LambdaLanczos(); + LambdaLanczos(std::function&, std::vector&)> mv_mul, int matrix_size, bool find_maximum); + LambdaLanczos(std::function&, std::vector&)> mv_mul, int matrix_size) : LambdaLanczos(mv_mul, matrix_size, true) {} + + /// @brief Calculate the principal (largest or smallest) eigenvalue + /// of the matrix (and its corresponding eigenvector). + int run(real_t&, std::vector&) const; + + // --- public data members --- + + /// @brief Specify the size of the matrix you will analyze. + /// (This equals the size of the eigenvector which will be returned.) + int matrix_size; + + /// @brief Specify the function used for matrix*vector multiplication + /// used by the Lanczos algorithm. For an ordinary dense matrix, + /// this function is the ordinary matrix*vector product. (See the + /// example above. For a sparse matrix, it will be something else.) + std::function&, std::vector&)> mv_mul; + + /// @brief Are we searching for the maximum or minimum eigenvalue? + /// @note (Usually, you must also specify eigenvalue_offset.) + bool find_maximum = false; + + /// @brief Shift all the eigenvalues by "eigenvalue_offset" during the Lanczos + /// iteration (ie. during LambdaLanczos::run()). The goal is to insure + /// that the correct eigenvalue is selected (the one with the maximum + /// magnitude). + /// @note The eigevalue returned by LambdaLanczos::run() is not effected + /// because after the iteration is finished, it will subtract this + /// number from the eigenvalue before it is returned to the caller. + /// @note Unless your matrix is positive definite or negative definite, + /// you MUST specify eigenvalue_offset. See comment above for details. + real_t eigenvalue_offset = 0.0; + + /// @brief This function sets "eigenvalue_offset" automatically. + /// @note Using this function is not recommended because it is very slow. + /// For efficiency, set the "eigenvalue_offset" yourself. + void ChooseOffset(); + + // The remaining data members usually can be left alone: + int max_iteration; + real_t eps = minimum_effective_decimal>() * 1e3; + real_t tridiag_eps_ratio = 1e-1; + int initial_vector_size = 200; + std::function&)> init_vector = + VectorRandomInitializer::init; + + // (for those who prefer "Set" functions...) + int SetSize(int matrix_size); + void SetMul(std::function&, + std::vector&)> mv_mul); + void SetInitVec(std::function&)> init_vector); + void SetFindMax(bool find_maximum); + void SetEvalOffset(T eigenvalue_offset); + void SetEpsilon(T eps); + void SetTriEpsRatio(T tridiag_eps_ratio); + + private: + static void schmidt_orth(std::vector&, const std::vector>&); + real_t find_minimum_eigenvalue(const std::vector>&, + const std::vector>&) const; + real_t find_maximum_eigenvalue(const std::vector>&, + const std::vector>&) const; + static real_t tridiagonal_eigen_limit(const std::vector>&, + const std::vector>&); + static int num_of_eigs_smaller_than(real_t, + const std::vector>&, + const std::vector>&); + real_t UpperBoundEvals() const; + }; + + + + /// @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 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 + class PEigenDense + { + size_t n; // the size of the matrix + std::vector evec; // preallocated vector + + public: + PEigenDense(int matrix_size=0); + + /// @brief Calculate the principal eigenvalue and eigenvector of a matrix. + /// @return Return the principal eigenvalue of the matrix. + /// If you want the eigenvector, pass a non-null "evector" argument. + Scalar + PrincipalEigen(ConstMatrix matrix, //!< the input patrix + Vector evector, //!< the eigenvector is stored here + bool find_max=false); //!< want the max or min eigenvalue? + + void SetSize(int matrix_size); //change matrix size after instantiation + + }; // class PEigenDense @@ -460,16 +527,13 @@ public: // ----------- IMPLEMENTATION ----------- // -------------------------------------- - - - // --- Implementation: Memory allocation for matrices --- template 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 l1_norm(const std::vector& vec) { // --- Implementation: Eigendecomposition of small dense matrices --- +template +Jacobi:: +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 +Jacobi:: +~Jacobi() { + if (! is_M_preallocated) + Dealloc(); +} + template int Jacobi:: 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& 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). @@ -912,8 +998,8 @@ swap(Jacobi &other) { template Jacobi:: Jacobi(Jacobi&& other) { - Init(); - swap(*this, other); + Init(); + this->swap(other); } // Using the "copy-swap" idiom for the assignment operator @@ -955,8 +1041,8 @@ template inline int LambdaLanczos:: run(real_t& eigvalue, std::vector& 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> u; // Lanczos vectors std::vector> alpha; // Diagonal elements of an approximated tridiagonal matrix @@ -1321,13 +1407,26 @@ init(std::vector>& v) // --- Implementation of PEigenDense +template +void PEigenDense:: +SetSize(int matrix_size) { + n = matrix_size; + evec.resize(n); +} + +template +PEigenDense:: +PEigenDense(int matrix_size):evec(matrix_size) { + SetSize(matrix_size); +} + template Scalar PEigenDense:: PrincipalEigen(ConstMatrix matrix, Vector eigenvector, bool find_max) { - assert(n > 0); + //assert(n > 0); auto matmul = [&](const std::vector& in, std::vector& 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 diff --git a/src/math_extra.cpp b/src/math_extra.cpp index b04414a300..a7a8a98a90 100644 --- a/src/math_extra.cpp +++ b/src/math_extra.cpp @@ -19,8 +19,6 @@ #include #include -#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 diff --git a/src/math_extra.h b/src/math_extra.h index a818bae4f4..6c58bc93c2 100644 --- a/src/math_extra.h +++ b/src/math_extra.h @@ -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); diff --git a/src/molecule.cpp b/src/molecule.cpp index 18f5c38992..9042bf0521 100644 --- a/src/molecule.cpp +++ b/src/molecule.cpp @@ -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];