diff --git a/doc/src/pg_developer.rst b/doc/src/pg_developer.rst index d67b69a80c..4f58eeec6b 100644 --- a/doc/src/pg_developer.rst +++ b/doc/src/pg_developer.rst @@ -1286,6 +1286,3 @@ elements of a sparse matrix.) :project: progguide :members: -.. doxygenstruct:: MathEigen::realTypeMap - :project: progguide - :members: diff --git a/src/math_eigen.cpp b/src/math_eigen.cpp index 9255338e09..9e98a161b6 100644 --- a/src/math_eigen.cpp +++ b/src/math_eigen.cpp @@ -90,16 +90,18 @@ jacobi3(double const mat[3][3], // the 3x3 matrix you wish to diagonalize {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]; // (another array which is preallocated to avoid using the heap) + // variable "ecalc3" does all the work: - Jacobi ecalc3(3,M); + Jacobi ecalc3(3, M, midx); 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]); + for (int j=i+1; j<3; j++) + std::swap(evec[i][j], evec[j][i]); // transpose the evec matrix return ierror; } @@ -124,16 +126,18 @@ jacobi3(double const* const* mat, // the 3x3 matrix you wish to diagonalize {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]; // (another array which is preallocated to avoid using the heap) + // variable "ecalc3" does all the work: - Jacobi ecalc3(3,M); + Jacobi ecalc3(3, M, midx); 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]); + for (int j=i+1; j<3; j++) + std::swap(evec[i][j], evec[j][i]); // transpose the evec matrix return ierror; } diff --git a/src/math_eigen.h b/src/math_eigen.h index a0bd292263..8bb61111f3 100644 --- a/src/math_eigen.h +++ b/src/math_eigen.h @@ -150,9 +150,7 @@ namespace MathEigen { 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(int n); ~Jacobi(); @@ -185,9 +183,25 @@ namespace MathEigen { int max_num_sweeps=50 //!< limit the number of iterations ); + // An alternative constructor is provided, if, for some reason, + // you want to avoid allocating memory on the heap during + // initialization. In that case, you can allocate space for the "M" + // and "max_idx_row" arrays on the stack in advance, and pass them + // to the constructor. (The vast majority of users probably + // do not need this feature. Unless you do, I encourage you + // to use the regular constructor instead.) + // n the size (ie. number of rows) of the (square) matrix. + // M optional preallocated n x n array + // max_idx_row optional preallocated array of size n + // note If either the "M" or "max_idx_row" arguments are specified, + // they both must be specified. + Jacobi(int n, Scalar **M, int *max_idx_row); + private: - bool is_M_preallocated; + bool is_preallocated; + // (Descriptions of private functions can be found in their implementation.) + void _Jacobi(int n, Scalar **M, int *max_idx_row); 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); @@ -240,9 +254,9 @@ namespace MathEigen { double *eval, double evec[3][3], // optional arguments: - bool evec_as_columns=true, + bool evec_as_columns = true, Jacobi:: - SortCriteria sort_criteria= + SortCriteria sort_criteria = Jacobi:: SORT_DECREASING_EVALS ); @@ -505,18 +519,40 @@ inline real_t l1_norm(const std::vector& vec) { // --- Implementation: Eigendecomposition of small dense matrices --- + template Jacobi:: -Jacobi(int n, Scalar **workspace) { +Jacobi(int n) { + _Jacobi(n, nullptr, nullptr); +} + + +template +Jacobi:: +Jacobi(int n, Scalar **M, int *max_idx_row) { + _Jacobi(n, M, max_idx_row); +} + + +// _Jacobi() is a function which is invoked by the two constructors and it +// does all of the work. (Splitting the constructor into multiple functions +// was technically unnecessary, but it makes documenting the code easier.) + +template +void +Jacobi:: +_Jacobi(int n, Scalar **M, int *max_idx_row) { Init(); - if (workspace) { - is_M_preallocated = true; - M = workspace; + if (M) { // if caller supplies their own "M" array, don't allocate M + is_preallocated = true; this->n = n; + this->M = M; + this->max_idx_row = max_idx_row; + //assert(this->max_idx_row); } else { - is_M_preallocated = false; - SetSize(n); + is_preallocated = false; + SetSize(n); // allocate the "M" and "max_int_row" arrays } } @@ -524,7 +560,7 @@ Jacobi(int n, Scalar **workspace) { template Jacobi:: ~Jacobi() { - if (! is_M_preallocated) + if (! is_preallocated) Dealloc(); } @@ -829,6 +865,7 @@ Init() { n = 0; M = nullptr; max_idx_row = nullptr; + is_preallocated = false; } template @@ -843,6 +880,7 @@ SetSize(int n) { template void Jacobi:: Alloc(int n) { + //assert(! is_preallocated); this->n = n; if (n > 0) { max_idx_row = new int[n]; @@ -853,8 +891,7 @@ Alloc(int n) { template void Jacobi:: Dealloc() { - if (max_idx_row) - delete [] max_idx_row; + //assert(! is_preallocated); Dealloc2D(&M); Init(); } @@ -867,6 +904,7 @@ 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 @@ -884,6 +922,7 @@ template void Jacobi:: swap(Jacobi &other) { std::swap(n, other.n); + std::swap(is_preallocated, other.is_preallocated); std::swap(max_idx_row, other.max_idx_row); std::swap(M, other.M); }