diff --git a/lib/jacobi_pd.h b/lib/jacobi_pd.h deleted file mode 100644 index ca68e6ba3a..0000000000 --- a/lib/jacobi_pd.h +++ /dev/null @@ -1,623 +0,0 @@ -/// @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/USER-REACTION/superpose3d.h b/src/USER-REACTION/superpose3d.h index 63c5bf5de8..b0b290ddad 100644 --- a/src/USER-REACTION/superpose3d.h +++ b/src/USER-REACTION/superpose3d.h @@ -22,10 +22,10 @@ /// @author Andrew Jewett /// @license MIT -#ifndef _SUPERPOSE3D_H -#define _SUPERPOSE3D_H +#ifndef LMP_SUPERPOSE3D_H +#define LMP_SUPERPOSE3D_H -#include "math_eigen.h" //functions to calculate eigenvalues and eigenvectors +#include "math_eigen_impl.h" //functions to calculate eigenvalues and eigenvectors // ----------------------------------------------------------- // ------------------------ INTERFACE ------------------------ @@ -462,4 +462,4 @@ operator = (Superpose3D source) { } -#endif //#ifndef _SUPERPOSE3D_H +#endif //#ifndef LMP_SUPERPOSE3D_H diff --git a/src/math_eigen.cpp b/src/math_eigen.cpp index 2caa6414d1..e2a811bdab 100644 --- a/src/math_eigen.cpp +++ b/src/math_eigen.cpp @@ -1,125 +1,62 @@ + #include "math_eigen.h" +#include "math_eigen_impl.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 >&>; - -// If you plan to use other numeric types (eg floats), add them to this list. - - - // 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) +typedef Jacobi Jacobi_v1; +typedef Jacobi Jacobi_v2; + +int MathEigen::jacobi3(double const mat[3][3], double *eval, double evec[3][3]) { - // 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**). + // make copy of const matrix + 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]) }; + int midx[3]; - int midx[3]; // (another array which is preallocated to avoid using the heap) + // create instance of generic Jacobi class and get eigenvalues and -vectors - // variable "ecalc3" does all the work: - Jacobi ecalc3(3, M, midx); - int ierror = ecalc3.Diagonalize(mat, eval, evec, sort_criteria); + Jacobi_v1 ecalc3(3, M, midx); + int ierror = ecalc3.Diagonalize(mat, eval, evec, Jacobi_v1::SORT_DECREASING_EVALS); - // 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=i+1; j<3; j++) - std::swap(evec[i][j], evec[j][i]); // transpose the evec matrix + // transpose the evec matrix + + for (int i=0; i<3; i++) + for (int j=i+1; j<3; j++) + std::swap(evec[i][j], evec[j][i]); return ierror; } - - -int MathEigen:: -jacobi3(double const* 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) +int MathEigen::jacobi3(double const* const* mat, double *eval, double **evec) { - // 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**). + // make copy of const matrix + 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]) }; + int midx[3]; - int midx[3]; // (another array which is preallocated to avoid using the heap) + // create instance of generic Jacobi class and get eigenvalues and -vectors - // variable "ecalc3" does all the work: - Jacobi ecalc3(3, M, midx); - int ierror = ecalc3.Diagonalize(mat, eval, evec, sort_criteria); + Jacobi_v2 ecalc3(3, M, midx); + int ierror = ecalc3.Diagonalize(mat, eval, evec, Jacobi_v2::SORT_DECREASING_EVALS); - // 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=i+1; j<3; j++) - std::swap(evec[i][j], evec[j][i]); // transpose the evec matrix + // transpose the evec matrix + + for (int i=0; i<3; i++) + for (int j=i+1; j<3; j++) + std::swap(evec[i][j], evec[j][i]); return ierror; } diff --git a/src/math_eigen.h b/src/math_eigen_impl.h similarity index 96% rename from src/math_eigen.h rename to src/math_eigen_impl.h index 8bb61111f3..5f7dc0a255 100644 --- a/src/math_eigen.h +++ b/src/math_eigen_impl.h @@ -16,8 +16,8 @@ Andrew Jewett (Scripps Research, Jacobi algorithm) ------------------------------------------------------------------------- */ -#ifndef LMP_MATH_EIGEN_H -#define LMP_MATH_EIGEN_H +#ifndef LMP_MATH_EIGEN_IMPL_H +#define LMP_MATH_EIGEN_IMPL_H // This file contains a library of functions and classes which can // efficiently perform eigendecomposition for an extremely broad @@ -220,50 +220,6 @@ namespace MathEigen { }; // 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* 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 **) - /// @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], - double *eval, - double evec[3][3], - // optional arguments: - bool evec_as_columns = true, - Jacobi:: - SortCriteria sort_criteria = - Jacobi:: - SORT_DECREASING_EVALS - ); - - - - // ---- Eigendecomposition of large sparse (or dense) matrices ---- // The "LambdaLanczos" is a class useful for calculating eigenvalues @@ -1413,4 +1369,4 @@ PrincipalEigen(ConstMatrix matrix, -#endif //#ifndef LMP_MATH_EIGEN_H +#endif //#ifndef LMP_MATH_EIGEN_IMPL_H