diff --git a/doc/eigensolver.html b/doc/eigensolver.html
new file mode 100644
index 000000000..63988173b
--- /dev/null
+++ b/doc/eigensolver.html
@@ -0,0 +1,146 @@
+
+
+
+Eigen Solver
+
+
+
+Description
+The class eigen_solver is used to compute the eigenvalues and eigenvectors of real matrices.
+Example
+
+#include <boost/numeric/ublas/storage.hpp>
+
+#include <boost/numeric/ublas/matrix.hpp>
+#include <boost/numeric/ublas/io.hpp>
+#include <boost/numeric/ublas/eigen_solver.hpp>
+#include <complex>
+
+int main () {
+ using namespace boost::numeric::ublas;
+ matrix m(3, 3);
+ m <<= 2.0, 3.0, 5.0,
+ 2.0, -3.0, 7.0,
+ 4.0, 1.0, 1.0;
+ eigen_solver > es(m,EIGVEC);
+
+ matrix evals_r = es.get_real_eigenvalues();
+ matrix evals_c = es.get_complex_eigenvalues();
+ matrix evecs_r = es.get_real_eigenvectors();
+ matrix evecs_c = es.get_complex_eigenvectors();
+
+ std::cout << "Eigenvalues (Real Part)\n";
+ std::cout << evals_r << std::endl;
+
+ std::cout << "Eigenvalues (Imag. Part)\n";
+ std::cout << evals_c << std::endl;
+
+ std::cout << "Eigenvectors (Real Part)\n";
+ std::cout << evecs_r << std::endl;
+
+ std::cout << "Eigenvectors (Imag. Part)\n";
+ std::cout << evecs_c << std::endl;
+
+ std::cout << "Verification\n";
+ matrix > V(3,3);
+ matrix > D(3,3);
+ matrix > A(3,3);
+ matrix > Lambda;
+ for (int i = 0; i < 3; i++){
+ for (int j = 0; j < 3; j++){
+ V(i, j) = std::complex(evecs_r(i, j), evecs_c(i, j));
+ D(i, j) = std::complex(evals_r(i, j), evals_c(i, j));
+ A(i, j) = std::complex(m(i, j), 0.0);
+ }
+ }
+ Lambda = prod(A, V) - prod(V, D);
+ std::cout << Lambda << std::endl;
+}
+
+
+Definition
+Defined in the header eigen_solver.hpp.
+Model of
+None.
+Type requirements
+Expected to be of type double or float. int may give unexpected results.
+Public base classes
+None.
+Members
+
+
+
+Member
+Description
+
+
+eigen_solver(M &m, eig_solver_params params = EIGVAL)Default constructor that takes a object of type M and returns either its eigenvalues or both eigenvalues and eigenvectors.
+
+
+void compute(eig_solver_params params = EIGVAL)Mostly internal method that starts the computation of eigenvalues. To be used in case the user wanted EigenValues option in the constructor and then later on wants the Eigenvectors as well, then we should call the compute with the EIGVEC option.
+
+
+M& get_real_eigenvalues()Returns the real portion of the Eigenvalues.
+
+
+M& get_complex_eigenvalues()Returns the complex portion of the Eigenvalues.
+
+
+M& get_real_eigenvectors()Returns the real portion of the Eigenvectors(if computed).
+
+
+M& get_complex_eigenvectors()Returns the complex portion of the Eigenvectors(if computed).
+
+
+bool has_complex_eigenvalues()Returns if there is any complex portion to the eigen values.
+
+
+M& get_real_schur_form()Returns the real schur decomposition of the hessenberg form of the original matrix.
+
+
+M& get_hessenberg_form()Returns the hessenberg decomposition of the original matrix.
+
+
+
+eig_solver_params Param Values
+EIGVAL - Just computes the eigenvalues of the matrix.EIGVEC - Computes the eigenvalues as well as eigenvectors of the matrix.
+
+
+
+ Copyright (©) 2000-2004 Michael Stevens, Mathias Koch,
+ Joerg Walter, Gunter Winkler
+ http://www.boost.org/LICENSE_1_0.txt
+ ).
+
+
+
+
diff --git a/doc/index.html b/doc/index.html
index bedf58133..976fd2b2e 100644
--- a/doc/index.html
+++ b/doc/index.html
@@ -347,6 +347,12 @@ Documentation
+
+
+ Numerical Algebra
+
+
+
Release notes
diff --git a/doc/samples/eigen_solver.cpp b/doc/samples/eigen_solver.cpp
new file mode 100644
index 000000000..f21812116
--- /dev/null
+++ b/doc/samples/eigen_solver.cpp
@@ -0,0 +1,54 @@
+//
+// Rajaditya Mukherjee
+//
+// Distributed under the Boost Software License, Version 1.0. (See
+// accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+//
+
+#include
+#include
+#include
+#include
+
+int main () {
+ using namespace boost::numeric::ublas;
+ matrix m(3, 3);
+ m <<= 2.0, 3.0, 5.0,
+ 2.0, -3.0, 7.0,
+ 4.0, 1.0, 1.0;
+ eigen_solver > es(m,EIGVEC);
+
+ matrix evals_r = es.get_real_eigenvalues();
+ matrix evals_c = es.get_complex_eigenvalues();
+ matrix evecs_r = es.get_real_eigenvectors();
+ matrix evecs_c = es.get_complex_eigenvectors();
+
+ std::cout << "Eigenvalues (Real Part)\n";
+ std::cout << evals_r << std::endl;
+
+ std::cout << "Eigenvalues (Imag. Part)\n";
+ std::cout << evals_c << std::endl;
+
+ std::cout << "Eigenvectors (Real Part)\n";
+ std::cout << evecs_r << std::endl;
+
+ std::cout << "Eigenvectors (Imag. Part)\n";
+ std::cout << evecs_c << std::endl;
+
+ std::cout << "Verification\n";
+ matrix > V(3,3);
+ matrix > D(3,3);
+ matrix > A(3,3);
+ matrix > Lambda;
+ for (int i = 0; i < 3; i++){
+ for (int j = 0; j < 3; j++){
+ V(i, j) = std::complex(evecs_r(i, j), evecs_c(i, j));
+ D(i, j) = std::complex(evals_r(i, j), evals_c(i, j));
+ A(i, j) = std::complex(m(i, j), 0.0);
+ }
+ }
+ Lambda = prod(A, V) - prod(V, D);
+ std::cout << Lambda << std::endl;
+}
+
diff --git a/include/boost/numeric/ublas/eigen_solver.hpp b/include/boost/numeric/ublas/eigen_solver.hpp
new file mode 100644
index 000000000..fd50c19d0
--- /dev/null
+++ b/include/boost/numeric/ublas/eigen_solver.hpp
@@ -0,0 +1,395 @@
+// Rajaditya Mukherjee
+//
+// Distributed under the Boost Software License, Version 1.0. (See
+// accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+
+/// \file eigen_solver.hpp Contains method for computing the eigenvalues of real matrices
+
+#ifndef _BOOST_UBLAS_EIGENSOLVER_
+#define _BOOST_UBLAS_EIGENSOLVER_
+
+#include
+
+#include
+#include
+#include
+#include
+#include
+#include
+
+namespace boost { namespace numeric { namespace ublas {
+
+//Add some sloppily defined enums here : later on ask mentor how we can use this to
+//Ask mentor about how to "formally" add them to ublas
+enum eig_solver_params {EIGVAL,EIGVEC};
+
+/** \brief Computes Eigenvalues and Eigenvectors for matrix type M \c T.
+* Given a matrix type \c M, this class computes the EigenValues and (optionally) Eigenvectors. The input matrix is expected to be real.
+* The output Eigenvalues as well as Eigenvectors can be both real or complex.
+* \tparam M type of the objects stored in the vector. Expected to be a matrix (like matrix)
+*/
+
+template
+class eigen_solver {
+private:
+ M matrix;
+ M hessenberg_form;
+ M real_schur_form;
+ M transform_accumulations;
+ bool has_complex_part;
+ bool has_eigenvalues;
+ bool has_eigenvectors;
+ M eigenvalues_real;
+ M eigenvalues_complex;
+ M eigenvectors;
+ M eigenvectors_real;
+ M eigenvectors_complex;
+ eig_solver_params solver_params;
+public:
+ /// \brief Constructor that takes a matrix and an (optional) argument whether to extract eigenvectors also.
+ /// The only argument that it needs is the matrix for which we wish to compute the Eigenvalues (and optionally Eigenvectors). The second option
+ /// specifies whether we wish to compute only the Eigenvalues (EIGVAL) or both Eigenvalues as well as Eigenvectors(EIGVEC). The default option
+ /// is that it only computes the Eigenvalues.
+ /// \param m Matrix for which we need to extract the Eigenvalues
+ /// \param params Can have two values EIGVAL or EIGVEC.
+ BOOST_UBLAS_INLINE
+ explicit eigen_solver(M &m, eig_solver_params params = EIGVAL) :
+ matrix(m), has_complex_part(false), solver_params(params)
+ {
+ has_eigenvalues = has_eigenvectors = false;
+ hessenberg_form = M(m);
+ compute(solver_params);
+ }
+
+ /// \brief Method to order the class to (re)compute the Eigenvalues and Eigenvectors.
+ /// Normally the user doesn't need to call this method explicitly as it called with the constructor. However, if initially the second argument
+ /// specified EIGVAL option and then the user wants to compute the Eigenvectors also, he will need to call this to recompute the
+ /// Eigenvectors with the EIGVEC option.
+ /// \param params Can have two values EIGVAL or EIGVEC.
+ BOOST_UBLAS_INLINE
+ void compute(eig_solver_params params = EIGVAL) {
+ if (params != solver_params) {
+ solver_params = params;
+ }
+
+ //If only eigen values are desired
+ if (params == EIGVAL) {
+ to_hessenberg(hessenberg_form);
+ real_schur_form = M(hessenberg_form);
+ schur_decomposition(real_schur_form);
+ extract_eigenvalues_from_schur();
+ has_eigenvalues = true;
+ has_eigenvectors = false;
+ }
+ else {
+ to_hessenberg(hessenberg_form, transform_accumulations);
+ real_schur_form = M(hessenberg_form);
+ schur_decomposition(real_schur_form, transform_accumulations);
+ extract_eigenvalues_from_schur();
+ extract_eigenvectors();
+ has_eigenvalues = true;
+ has_eigenvectors = true;
+ }
+
+
+ }
+
+ BOOST_UBLAS_INLINE
+ void extract_eigenvalues_from_schur() {
+ typedef typename M::size_type size_type;
+ typedef typename M::value_type value_type;
+
+ size_type n = real_schur_form.size1();
+ size_type i = size_type(0);
+
+ eigenvalues_real = zero_matrix(n,n);
+ eigenvalues_complex = zero_matrix(n, n);
+
+ while (i < n) {
+ if ((i == n - size_type(1)) || (real_schur_form(i + 1, i) == value_type(0)) || ((std::abs)(real_schur_form(i + 1, i)) <= 1.0e-5)) {
+ eigenvalues_real(i, i) = real_schur_form(i, i);
+ i += size_type(1);
+ }
+ else {
+ value_type p = value_type(0.5) * (real_schur_form(i, i) - real_schur_form(i + size_type(1), i + size_type(1)));
+ value_type z = (std::sqrt)((std::abs)(p * p + real_schur_form(i + size_type(1), i) * real_schur_form(i, i + size_type(1))));
+ eigenvalues_real(i, i) = real_schur_form(i + size_type(1), i + size_type(1)) + p;
+ eigenvalues_complex(i, i) = z;
+ eigenvalues_real(i + size_type(1), i + size_type(1)) = real_schur_form(i + size_type(1), i + size_type(1)) + p;
+ eigenvalues_complex(i + size_type(1), i + size_type(1)) = -z;
+ i += size_type(2);
+ has_complex_part = true;
+ }
+ }
+
+ has_eigenvalues = true;
+ }
+
+ /// \brief Method to get the real portion of the Eigenvalues. Output type same as input type.
+ BOOST_UBLAS_INLINE
+ M& get_real_eigenvalues() {
+ return eigenvalues_real;
+ }
+
+ /// \brief Method to get the imaginary portion of the Eigenvalues. Output type same as input type.
+ BOOST_UBLAS_INLINE
+ M& get_complex_eigenvalues() {
+ return eigenvalues_complex;
+ }
+
+ /// \brief Method to get the real portion of the Eigenvectors. Output type same as input type.
+ BOOST_UBLAS_INLINE
+ M& get_real_eigenvectors() {
+ return eigenvectors_real;
+ }
+
+ /// \brief Method to get the imaginary portion of the Eigenvectors. Output type same as input type.
+ BOOST_UBLAS_INLINE
+ M& get_complex_eigenvectors() {
+ return eigenvectors_complex;
+ }
+
+
+ BOOST_UBLAS_INLINE
+ void extract_eigenvectors() {
+
+ typedef typename M::size_type size_type;
+ typedef typename M::value_type value_type;
+
+ M T(real_schur_form);
+ eigenvectors = M(transform_accumulations);
+
+ size_type n = eigenvalues_real.size1();
+
+ value_type norm = value_type(0);
+ for (size_type j = size_type(0); j < n; j++) {
+ size_type idx = (j==size_type(0))?size_type(0):(j - size_type(1));
+ vector v = project(row(T,j), range(idx, n));
+ norm += norm_1(v);
+ }
+
+ if (norm == value_type(0))
+ return;
+
+ for (size_type k = n; k-- != size_type(0);) {
+
+ value_type p = eigenvalues_real(k, k);
+ value_type q = eigenvalues_complex(k, k);
+
+ if (q == value_type(0)) {
+
+ value_type lastr, lastw;
+ lastr = lastw = value_type(0);
+ size_type l = k;
+
+ T(k, k) = value_type(1);
+
+ for (size_type i = k; i-- != size_type(0);) {
+ value_type w = T(i, i) - p;
+ vector r_left = project(row(T, i), range(l, k + size_type(1)));
+ vector r_right = project(column(T, k), range(l, k + size_type(1)));
+ value_type r = inner_prod(r_left, r_right);
+ if (eigenvalues_complex(i, i) < value_type(0)) {
+ lastw = w;
+ lastr = r;
+ }
+ else {
+ l = i;
+ if (eigenvalues_complex(i, i) == value_type(0)) {
+ if (w != value_type(0)) {
+ T(i, k) = -r / w;
+ }
+ else {
+ T(i, k) = -r / (norm * std::numeric_limits::epsilon());
+ }
+ }
+ else {
+ value_type x = T(i, i + size_type(1));
+ value_type y = T(i + size_type(1), i);
+ value_type denom = (eigenvalues_real(i, i) - p)*(eigenvalues_real(i, i) - p) + (eigenvalues_complex(i, i))*(eigenvalues_complex(i, i));
+ value_type t = (x * lastr - lastw * r) / denom;
+ T(i, k) = t;
+ if ((std::abs)(x) > (std::abs)(lastw)) {
+ T(i + size_type(1), k) = (-r - w * t) / x;
+ }
+ else {
+ T(i + size_type(1), k) = (-lastr - y * t) / lastw;
+ }
+ }
+
+ // Overflow control
+ value_type t = (std::abs)(T(i, k));
+ if ((std::numeric_limits::epsilon() * t) * t > value_type(1)) {
+ for (size_type j = n - k - i; j < n; j++)
+ T(k, j) /= t;
+ }
+ }
+ }
+ }
+ else if (q < value_type(0) && k>0) {
+
+ value_type lastra, lastsa, lastw;
+ lastra = lastsa = lastw = value_type(0);
+ size_type l = k - 1;
+
+ if ((std::abs)(T(k, k - size_type(1))) > (std::abs)(T(k - size_type(1), k))) {
+ T(k - size_type(1), k - size_type(1)) = q / T(k, k - size_type(1));
+ T(k - size_type(1), k) = -(T(k, k) - p) / T(k, k - size_type(1));
+ }
+ else {
+ std::complex x(value_type(0), -T(k - size_type(1), k));
+ std::complex y(T(k - size_type(1), k - size_type(1)) - p, q);
+ std::complex cc = x / y;
+ T(k - size_type(1), k - size_type(1)) = cc.real();
+ T(k - size_type(1), k) = cc.imag();
+ }
+ T(k ,k - size_type(1)) = value_type(0);
+ T(k, k) = value_type(1);
+
+ for (size_type i = k - size_type(1); i-- != size_type(0);) {
+
+ vector ra_left = project(row(T, i), range(l, k + size_type(1)));
+ vector ra_right = project(column(T, k - size_type(1)), range(l, k + size_type(1)));
+ value_type ra = inner_prod(ra_left, ra_right);
+
+ vector sa_right = project(column(T, k), range(l, k + size_type(1)));
+ value_type sa = inner_prod(ra_left, sa_right);
+
+ value_type w = T(i, i) - p;
+
+ if (eigenvalues_complex(i, i) < value_type(0)) {
+ lastw = w;
+ lastra = ra;
+ lastsa = sa;
+ }
+ else {
+ l = i;
+ if (eigenvalues_complex(i, i) == value_type(0)) {
+ std::complex x(-ra, -sa);
+ std::complex y(w, q);
+ std::complex cc = x / y;
+ T(i, k - size_type(1)) = cc.real();
+ T(i, k) = cc.imag();
+ }
+ else {
+ value_type x = T(i, i + size_type(1));
+ value_type y = T(i + size_type(1), i);
+ value_type vr = (eigenvalues_real(i, i) - p) * (eigenvalues_real(i, i) - p) + eigenvalues_complex(i, i) * eigenvalues_complex(i, i) - q * q;
+ value_type vi = (eigenvalues_real(i, i) - p) * value_type(2) * q;
+
+ if (vr == value_type(0) && vi == value_type(0)) {
+ vr = (std::numeric_limits::epsilon()) * norm * ((std::abs)(w)+(std::abs)(q)+(std::abs)(x)+(std::abs)(y)+(std::abs)(lastw));
+ }
+
+ std::complex x1(x*lastra - lastw*ra + q*sa, x*lastsa - lastw*sa - q*ra);
+ std::complex y1(vr, vi);
+ std::complex cc = x1 / y1;
+ T(i, k - size_type(1)) = cc.real();
+ T(i, k) = cc.imag();
+
+ if ((std::abs)(x) > ((std::abs)(lastw)+(std::abs)(q))) {
+ T(i + size_type(1), k - size_type(1)) = (-ra - w * T(i, k - size_type(1)) + q * T(i, k)) / x;
+ T(i + size_type(1), k) = (-sa - w * T(i, k) - q * T(i, k - size_type(1))) / x;
+ }
+ else {
+ x1 = std::complex(-lastra - y*T(i, k - size_type(1)), -lastsa - y*T(i, k));
+ y1 = std::complex(lastw, q);
+ cc = x1 / y1;
+ T(i + size_type(1), k - size_type(1)) = cc.real();
+ T(i + size_type(1), k) = cc.imag();
+ }
+
+ }
+
+ }
+
+ value_type t = (std::max)((std::abs)(T(i, k - size_type(1))), (std::abs)(T(i, k)));
+ if ((std::numeric_limits::epsilon() * t)*t > value_type(1)) {
+ for (size_type j = i; j < k + i - size_type(1); j++){
+ T(j, n - i) /= t;
+ T(j, n - i + size_type(1)) /= t;
+ }
+ }
+
+ }
+ k--;
+ }
+ }
+
+ for (size_type j = n; j--!=size_type(0);)
+ {
+ M matrix_left = project(eigenvectors, range(0, n),range(0, j + size_type(1)));
+ vector vector_right = project(column(T,j), range(0, j + size_type(1)));
+ vector v = prod(matrix_left,vector_right);
+ column(eigenvectors,j) = v;
+ }
+
+ //std::cout << eigenvectors << "\n";
+
+ eigenvectors_real = zero_matrix(n, n);
+ eigenvectors_complex = zero_matrix(n, n);
+
+ for (size_type j = 0; j < n; j++) {
+ if (j + size_type(1) == n || (eigenvalues_complex(j, j) == value_type(0))){
+ vector vj = column(eigenvectors, j);
+ value_type norm_vj = norm_2(vj);
+ vector normalized_vj = vj / norm_vj;
+ column(eigenvectors_real, j) = normalized_vj;
+ }
+ else {
+ vector > col_ij(n);
+ vector > col_ijp1(n);
+ for (size_type i = 0; i < n; ++i) {
+ col_ij(i) = std::complex(eigenvectors(i, j), eigenvectors(i, j + 1));
+ col_ijp1(i) = std::complex(eigenvectors(i, j), -eigenvectors(i, j + 1));
+ }
+
+ std::complex norm_ij = norm_2(col_ij);
+ vector > normalized_ij = col_ij / norm_ij;
+ std::complex norm_ijp1 = norm_2(col_ijp1);
+ vector > normalized_ijp1 = col_ijp1 / norm_ijp1;
+
+ for (size_type i = 0; i < n; ++i) {
+ eigenvectors_real(i, j) = normalized_ij(i).real();
+ eigenvectors_complex(i, j) = normalized_ij(i).imag();
+
+ eigenvectors_real(i, j+1) = normalized_ijp1(i).real();
+ eigenvectors_complex(i, j+1) = normalized_ijp1(i).imag();
+ }
+ j++;
+ }
+ }
+ }
+
+ /// \brief Method to get the real schur decomposition of input matrix. Output type same as input type.
+ BOOST_UBLAS_INLINE
+ M& get_real_schur_form() {
+ return real_schur_form;
+ }
+
+ /// \brief Method to get the Hessenberg Form of input matrix. Output type same as input type.
+ BOOST_UBLAS_INLINE
+ M& get_hessenberg_form() {
+ return hessenberg_form;
+ }
+
+ /// \brief Method to get the accumulated transformations for Hessenberg Form of input matrix. Output type same as input type.
+ BOOST_UBLAS_INLINE
+ M& get_transform_accumulations() {
+ return transform_accumulations;
+ }
+
+ /// \brief Method to check if the given matrix has real or complex Eigenvalues (and hence Eigenvectors).
+ BOOST_UBLAS_INLINE
+ bool has_complex_eigenvalues() {
+ return has_complex_part;
+ }
+
+};
+
+
+}}}
+
+
+#endif
diff --git a/include/boost/numeric/ublas/hessenberg.hpp b/include/boost/numeric/ublas/hessenberg.hpp
new file mode 100644
index 000000000..abc171a3a
--- /dev/null
+++ b/include/boost/numeric/ublas/hessenberg.hpp
@@ -0,0 +1,111 @@
+// Rajaditya Mukherjee
+//
+// Distributed under the Boost Software License, Version 1.0. (See
+// accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+
+/// \file hessenberg.hpp Definition for the methods performing Hessenberg Transformation
+
+#ifndef _BOOST_UBLAS_HESSENBERG_
+#define _BOOST_UBLAS_HESSENBERG_
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+#include
+
+namespace boost { namespace numeric { namespace ublas {
+
+ /// \brief Performs Hessenberg for Matrix \c m without accumulating transformations.
+ /// Replaces the input matrix by its Hessenberg Form which is a quasi-upper triangular matrix
+ /// preserved by QR Algorithms. It is easier and faster for the QR Algorithm to work on the Hessenberg form instead of the dense form.
+ /// We use a series of householder reflections to generate the Hessenberg Forms.
+ /// We use this form when we need only the eigen-values.
+ /// \param m matrix type (like matrix)
+ template
+ void to_hessenberg(M &m) {
+
+ typedef typename M::size_type size_type;
+ typedef typename M::value_type value_type;
+
+ //At this point, we should probably check if the number of rows and columns of M are the same or not
+ size_type num_rows = m.size1();
+ size_type num_cols = m.size2();
+ BOOST_UBLAS_CHECK(num_rows == num_cols, singular()); //Throw some kind of assertion error saying that eigen solver works only for sqaure matrices
+
+ size_type n = num_rows;
+ for (size_type i = 0; i < n-2; ++i) {
+ vector x = project(column(m,i), range(i + 1, n));
+ vector v;
+ value_type beta;
+ householder(x, v, beta);
+
+ matrix vvt = outer_prod(v, v);
+ vvt *= beta;
+ size_type n_vvt = vvt.size1();
+ matrix imvvt = identity_matrix(n_vvt) - vvt;
+
+ matrix t1 = prod(imvvt, project(m, range(i + 1, n), range(i, n)));
+ project(m, range(i + 1, n), range(i, n)).assign(t1);
+ matrix t2 = prod(project(m, range(0, n), range(i + 1, n)), imvvt);
+ project(m, range(0, n), range(i + 1, n)).assign(t2);
+ }
+
+ }
+
+ /// \brief Performs Hessenberg for Matrix \c m with accumulating transformations.
+ /// Replaces the input matrix by its Hessenberg Form which is a quasi-upper triangular matrix
+ /// preserved by QR Algorithms. It is easier and faster for the QR Algorithm to work on the Hessenberg form instead of the dense form.
+ /// We use a series of householder reflections to generate the Hessenberg Forms. In this format, we also accumulate the transformations which is needed
+ /// to calculate the eigen-vectors. We use this form when we need both the eigen-values as well as eigen-vectors.
+ /// \param m matrix type (like matrix)
+ /// \param u0 output with the transformations accumulated. (same type as M)
+ template
+ void to_hessenberg(M &m, M &u0) {
+
+ typedef typename M::size_type size_type;
+ typedef typename M::value_type value_type;
+
+ M m_copy = M(m);
+
+ //At this point, we should probably check if the number of rows and columns of M are the same or not
+ size_type num_rows = m.size1();
+ size_type num_cols = m.size2();
+ BOOST_UBLAS_CHECK(num_rows == num_cols, singular()); //Throw some kind of assertion error saying that eigen solver works only for sqaure matrices
+
+ size_type n = num_rows;
+ u0 = identity_matrix(n);
+ for (size_type i = 0; i < n - 2; ++i) {
+ vector x = project(column(m, i), range(i + 1, n));
+ vector v;
+ value_type beta;
+ householder(x, v, beta);
+
+ matrix vvt = outer_prod(v, v);
+ vvt *= beta;
+ size_type n_vvt = vvt.size1();
+ matrix imvvt = identity_matrix(n_vvt) -vvt;
+
+ matrix pk = identity_matrix(n);
+ project(pk, range(i+1, n), range(i+1, n)).assign(imvvt);
+
+ matrix t0 = prod(u0, pk);
+ u0 = t0;
+
+ matrix t1 = prod(imvvt, project(m, range(i + 1, n), range(i, n)));
+ project(m, range(i + 1, n), range(i, n)).assign(t1);
+ matrix t2 = prod(project(m, range(0, n), range(i + 1, n)), imvvt);
+ project(m, range(0, n), range(i + 1, n)).assign(t2);
+
+ }
+ }
+
+}}}
+
+
+#endif
\ No newline at end of file
diff --git a/include/boost/numeric/ublas/householder.hpp b/include/boost/numeric/ublas/householder.hpp
new file mode 100644
index 000000000..b14668b51
--- /dev/null
+++ b/include/boost/numeric/ublas/householder.hpp
@@ -0,0 +1,107 @@
+// Rajaditya Mukherjee
+//
+// Distributed under the Boost Software License, Version 1.0. (See
+// accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+
+/// \file householder.hpp Definition for the methods performing Householder and Givens Rotation
+
+#ifndef _BOOST_UBLAS_HOUSEHOLDER_
+#define _BOOST_UBLAS_HOUSEHOLDER_
+
+#include
+#include
+#include
+#include
+#include
+
+namespace boost {
+ namespace numeric {
+ namespace ublas {
+
+ /// \brief Performs Householder reflections on the vector x.
+ /// Given a vector \c x, returns a vector \c v and a scalar \c beta
+ /// which is used to introduce zeros in the vector.
+ /// More specifically we have (I - beta * v & v^T)x = ||x||e_1 where e_1 is a
+ /// unit vector.
+ /// \param x input vector type (like vector)
+ /// \param v output vector type (like vector)
+ /// \param beta scalar output (same as M::value_type) containing the multiplication factor
+ template
+ void householder(M &x, M &v, typename M::value_type &beta) {
+
+ typedef typename M::size_type size_type;
+ typedef typename M::value_type value_type;
+
+ // x must be a vector
+ // output will be a vector in R^n
+ size_type n = x.size();
+ value_type sigma = inner_prod(project(x, range(1, n)), project(x, range(1, n)));
+
+ v = M(x);
+ v(size_type(0)) = value_type(1);
+
+
+ if (sigma != value_type()) {
+ value_type mu = ((x(size_type(0)))*(x(size_type(0)))) + sigma;
+ mu = value_type(boost_numeric_ublas_sqrt(mu));
+
+ if (x(size_type(0)) <= value_type(0)) {
+ v(size_type(0)) = x(size_type(0)) - mu;
+ }
+ else {
+ v(size_type(0)) = (-sigma) / (x(size_type(0)) + mu);
+ }
+ beta = (value_type(2) * v(size_type(0)) * v(size_type(0))) / (v(size_type(0)) * v(size_type(0)) + sigma);
+ value_type inv_v_0 = value_type(1) / v(size_type(0));
+ v *= inv_v_0;
+ }
+ else {
+ if (x(size_type(0)) >= value_type(0)) {
+ beta = value_type(0);
+ }
+ else {
+ beta = value_type(-2);
+ }
+ }
+ }
+
+
+ /// \brief Determines the Givens Rotation Coefficient for two scalars [a,b].
+ /// Given a vector \c [a,b], returns two scalars \c c and \c d such that the 2x2 rotation matrix formed by the said values of
+ /// cosine and sine will rotate the vector [a,b] to form [r,0].
+ /// \param a input scalar type (like double)
+ /// \param b input scalar type (like double)
+ /// \param c output scalar type (like double) containing cosine for Givens Rotation
+ /// \param s output scalar type (like double) containing sine for Givens Rotation
+
+ template
+ void givens_rotation(M &a, M &b, M&c, M&s) {
+ if (b == M(0)) {
+ c = M(1);
+ s = M(0);
+ }
+ else {
+ if ((std::abs)(b) > (std::abs)(a)) {
+ M tau = -a / b;
+ M s_temp = M(1.0) + (tau*tau);
+ s_temp = (std::sqrt)(s_temp);
+
+ s = M(1.0) / s_temp;
+ c = s * tau;
+ }
+ else {
+ M tau = -b / a;
+ M c_temp = M(1.0) + (tau*tau);
+ c_temp = (std::sqrt)(c_temp);
+
+ c = M(1.0) / c_temp;
+ s = c * tau;
+ }
+ }
+ }
+
+
+}}}
+
+#endif
\ No newline at end of file
diff --git a/include/boost/numeric/ublas/matrix_balancing.hpp b/include/boost/numeric/ublas/matrix_balancing.hpp
new file mode 100644
index 000000000..f755ffc8c
--- /dev/null
+++ b/include/boost/numeric/ublas/matrix_balancing.hpp
@@ -0,0 +1,115 @@
+// Rajaditya Mukherjee
+//
+// Distributed under the Boost Software License, Version 1.0. (See
+// accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+
+/// \file matrix_balancing.hpp Contains methods for balancing a real-nonsymmetric matrix.
+
+#ifndef _BOOST_UBLAS_MATRIXBALANCING_
+#define _BOOST_UBLAS_MATRIXBALANCING_
+
+#include
+#include
+#include
+#include
+#include
+
+namespace boost {
+ namespace numeric {
+ namespace ublas {
+
+ /// \brief Method that computes the balancing multiplier vector.
+ /// The only argument that it needs is the matrix which we wish to balance. Note that if the matrix is symmetric, then the balancing vector will be all ones.
+ /// The output of this method is the diagonal matrix (we represent the matrix as a vector for space constaints).
+ /// \param m Matrix of type \c T for which we need to balance.
+ /// \param d Vector of type \c T which is the diagonal balancing matrix.
+ template
+ void matrix_balance(matrix &m, vector &d) {
+ typedef typename matrix::size_type size_type;
+
+ size_type n = m.size1();
+
+ T row_norm,col_norm;
+ bool converged = false;
+ d = scalar_vector(n, T(1));
+
+ while (!converged) {
+ T g, f, s;
+
+ converged = true;
+
+ for (size_type i = size_type(0); i < n; i++) {
+ row_norm = col_norm = T(0);
+ for (size_type j = size_type(0); j < n; j++){
+ if (j != i){
+ col_norm += (std::abs)(m(j, i));
+ row_norm += (std::abs)(m(i, j));
+ }
+ }
+
+ if (col_norm == T(0) || row_norm == T(0)) {
+ continue;
+ }
+
+ g = row_norm / T(2);
+ f = T(1);
+ s = col_norm + row_norm;
+
+ while (col_norm < g) {
+ f *= T(2);
+ col_norm *= T(4);
+ }
+
+ g = row_norm * T(2);
+
+ while (col_norm > g) {
+ f /= T(2);
+ col_norm /= T(4);
+ }
+
+ if ((row_norm + col_norm) < T(0.95)*s*f){
+ converged = false;
+
+ g = T(1) / f;
+
+ vector v = row(m, i);
+ v *= g;
+ row(m, i) = v;
+
+ v = column(m, i);
+ v *= f;
+ column(m, i) = v;
+
+ d(i) *= f;
+
+ }
+ }
+ }
+ }
+
+ /// \brief Method that applies the balancing transform.
+ /// Essentially given a matrix m, this returns m = md where d is a diagonal matrix whose diagonal elements are represented as a vector.
+ /// \param m Matrix of type \c T for which we need to balance.
+ /// \param d Vector of type \c T which is the diagonal balancing matrix.
+ template
+ void apply_transformation(matrix &m, vector &d) {
+ typedef typename matrix::size_type size_type;
+ size_type n = m.size1();
+ //Personal note: Ask Mentor is this is right???
+ assert(n == d.size());
+ T s;
+ for (size_type i = size_type(0); i < n; i++){
+ s = d(i);
+ vector r = row(m, i);
+ r *= s;
+ row(m, i) = r;
+ }
+ }
+
+
+}}}
+
+
+
+#endif
\ No newline at end of file
diff --git a/include/boost/numeric/ublas/schur_decomposition.hpp b/include/boost/numeric/ublas/schur_decomposition.hpp
new file mode 100644
index 000000000..6ff7ad4a6
--- /dev/null
+++ b/include/boost/numeric/ublas/schur_decomposition.hpp
@@ -0,0 +1,533 @@
+// Rajaditya Mukherjee
+//
+// Distributed under the Boost Software License, Version 1.0. (See
+// accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+
+
+/// \file schur_decomposition.hpp Contains methods for real schur decomposition
+
+#ifndef _BOOST_UBLAS_SCHURDECOMPOSITION_
+#define _BOOST_UBLAS_SCHURDECOMPOSITION_
+
+#define EPSA 1.0e-20
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+
+namespace boost {
+ namespace numeric {
+ namespace ublas {
+
+ /*! @fn schur_decomposition(M &h)
+ * @brief Performs Real Schur Decomposition for Matrix \c m (Version without EigenVectors).
+ * Replaces the input matrix which is assumed to be in Hessenberg Form by the
+ real schur form of the matrix. This uses the Francis Double Shift QR Algorithm to compute the real Schur Form.
+ * The diagonals of the schur form are either 1x1 blocks (real eigenvalues) or 2x2 blocks (complex eigenvalues).
+ * @tparam[in] m matrix type (like matrix) - input hessenberg form output - schur form
+ */
+ template
+ void schur_decomposition(M &h) {
+
+ typedef typename M::size_type size_type;
+ typedef typename M::value_type value_type;
+
+ size_type n = h.size1();
+ size_type p = n;
+ while (p > size_type(2)) {
+ size_type q = p - size_type(1);
+
+ value_type s = h(q - size_type(1), q - size_type(1)) + h(p - size_type(1), p - size_type(1));
+ value_type t = h(q - size_type(1), q - size_type(1)) * h(p - size_type(1), p - size_type(1)) - h(q - size_type(1), p - size_type(1)) * h(p - size_type(1), q - size_type(1));
+ value_type x = h(size_type(0), size_type(0)) * h(size_type(0), size_type(0)) +
+ h(size_type(0), size_type(1)) * h(size_type(1), size_type(0)) -
+ s * h(size_type(0), size_type(0)) + t;
+ value_type y = h(size_type(1), size_type(0)) * (h(size_type(0), size_type(0)) + h(size_type(1), size_type(1)) - s);
+ value_type z = h(size_type(1), size_type(0)) * h(size_type(2), size_type(1));
+
+ //Ask mentor if we can assume that size_type will support signed type (otherwise we need a different mechanism
+ // On second thought lets make no such assumptions
+ for (size_type k = size_type(0); k <= p - size_type(3); k++) {
+ vector tx(3);
+ tx(0) = x; tx(1) = y; tx(2) = z;
+ vector v;
+ value_type beta;
+ householder >(tx, v, beta);
+
+ size_type r = (std::max)(size_type(1), k);
+
+ matrix vvt = outer_prod(v, v);
+ vvt *= beta;
+ size_type n_vvt = vvt.size1();
+ matrix imvvt = identity_matrix(n_vvt) -vvt;
+
+ matrix t1 = prod(imvvt, project(h, range(k, k + size_type(3)), range(r - size_type(1), n)));
+ project(h, range(k, k + size_type(3)), range(r - size_type(1), n)).assign(t1);
+
+ r = (std::min)(k + size_type(4), p);
+
+ matrix t2 = prod(project(h, range(0, r), range(k, k + size_type(3))), imvvt);
+ project(h, range(0, r), range(k, k + size_type(3))).assign(t2);
+
+ x = h(k + size_type(1), k);
+ y = h(k + size_type(2), k);
+
+ if (k < (p - size_type(3))) {
+ z = h(k + size_type(3), k);
+ }
+
+ }
+
+ vector tx(2);
+ tx(0) = x; tx(1) = y;
+
+ value_type scale = (std::abs)(tx(0)) + (std::abs)(tx(1));
+ if (scale != value_type(0)) {
+ tx(0) /= scale;
+ tx(1) /= scale;
+ }
+
+ value_type cg, sg;
+ givens_rotation(tx(0), tx(1), cg, sg);
+
+ //Apply rotation to matrices as needed
+ matrix t1 = project(h, range(q - size_type(1), p), range(p - size_type(3), n));
+ size_type t1_cols = t1.size2();
+ size_type t1_rows = t1.size1();
+ for (size_type j = size_type(0); j < t1_cols; j++) {
+ value_type tau1 = t1(0, j);
+ value_type tau2 = t1(t1_rows - 1, j);
+ t1(0, j) = cg*tau1 - sg*tau2;
+ t1(t1_rows - 1, j) = sg*tau1 + cg*tau2;
+ }
+ //std::cout << t1 << "\n";
+ project(h, range(q - size_type(1), p), range(p - size_type(3), n)).assign(t1);
+
+ matrix t2 = project(h, range(0, p), range(p - size_type(2), p));
+ size_type t2_cols = t2.size2();
+ size_type t2_rows = t2.size1();
+ for (size_type j = size_type(0); j < t2_rows; j++) {
+ value_type tau1 = t2(j, 0);
+ value_type tau2 = t2(j, t2_cols - 1);
+ t2(j, 0) = cg*tau1 - sg*tau2;
+ t2(j, t2_cols - 1) = sg*tau1 + cg*tau2;
+ }
+ project(h, range(0, p), range(p - size_type(2), p)).assign(t2);
+
+ if ((std::abs)(h(p - size_type(1), q - size_type(1))) < (std::numeric_limits::epsilon())*((std::abs)(h(p - size_type(1), p - size_type(1))) + (std::abs)(h(q - size_type(1), q - size_type(1))))) {
+ h(p - size_type(1), q - size_type(1)) = value_type(0);
+ p = p - size_type(1);
+ q = p - size_type(1);
+ }
+ else if ((std::abs)(h(p - size_type(2), q - size_type(2))) < (std::numeric_limits::epsilon())*((std::abs)(h(q - size_type(2), q - size_type(2))) + (std::abs)(h(q - size_type(1), q - size_type(1))))) {
+ h(p - size_type(2), q - size_type(2)) = value_type(0);
+ p = p - size_type(2);
+ q = p - size_type(1);
+ }
+
+ }
+ }
+
+ /// \brief Performs Real Schur Decomposition for Matrix \c m (Version with EigenVectors).
+ /// Replaces the input matrix which is assumed to be in Hessenberg Form by the
+ /// real schur form of the matrix. This uses the Francis Double Shift QR Algorithm to compute the real Schur Form.
+ /// The diagonals of the schur form are either 1x1 blocks (real eigenvalues) or 2x2 blocks (complex eigenvalues).
+ /// \param m matrix type (like matrix) - input hessenberg form output - schur form
+ /// \param qv matrix type (like matrix) - input accumulated transforms
+ template
+ void schur_decomposition(M &h, M &qv) {
+
+ typedef typename M::size_type size_type;
+ typedef typename M::value_type value_type;
+
+ size_type n = h.size1();
+ size_type p = n;
+ while (p > size_type(2)) {
+ size_type q = p - size_type(1);
+
+ value_type s = h(q - size_type(1), q - size_type(1)) + h(p - size_type(1), p - size_type(1));
+ value_type t = h(q - size_type(1), q - size_type(1)) * h(p - size_type(1), p - size_type(1)) - h(q - size_type(1), p - size_type(1)) * h(p - size_type(1), q - size_type(1));
+ value_type x = h(size_type(0), size_type(0)) * h(size_type(0), size_type(0)) +
+ h(size_type(0), size_type(1)) * h(size_type(1), size_type(0)) -
+ s * h(size_type(0), size_type(0)) + t;
+ value_type y = h(size_type(1), size_type(0)) * (h(size_type(0), size_type(0)) + h(size_type(1), size_type(1)) - s);
+ value_type z = h(size_type(1), size_type(0)) * h(size_type(2), size_type(1));
+
+ //Ask mentor if we can assume that size_type will support signed type (otherwise we need a different mechanism
+ // On second thought lets make no such assumptions
+ for (size_type k = size_type(0); k <= p - size_type(3); k++) {
+ vector tx(3);
+ tx(0) = x; tx(1) = y; tx(2) = z;
+
+ //Add some scaling factor here
+ value_type scale = (std::abs)(tx(0)) + (std::abs)(tx(1)) + (std::abs)(tx(2));
+ if (scale != value_type(0)) {
+ tx(0) /= scale;
+ tx(1) /= scale;
+ tx(2) /= scale;
+ }
+
+ vector v;
+ value_type beta;
+ householder >(tx, v, beta);
+
+ size_type r = (std::max)(size_type(1), k);
+
+ matrix vvt = outer_prod(v, v);
+ vvt *= beta;
+ size_type n_vvt = vvt.size1();
+ matrix imvvt = identity_matrix(n_vvt) -vvt;
+
+ matrix t1 = prod(imvvt, project(h, range(k, k + size_type(3)), range(r - size_type(1), n)));
+ project(h, range(k, k + size_type(3)), range(r - size_type(1), n)).assign(t1);
+
+ r = (std::min)(k + size_type(4), p);
+
+ matrix t2 = prod(project(h, range(0, r), range(k, k + size_type(3))), imvvt);
+ project(h, range(0, r), range(k, k + size_type(3))).assign(t2);
+
+ matrix t3 = prod(project(qv, range(0, n), range(k, k + size_type(3))), imvvt);
+ project(qv, range(0, n), range(k, k + size_type(3))).assign(t3);
+
+ x = h(k + size_type(1), k);
+ y = h(k + size_type(2), k);
+
+ if (k < (p - size_type(3))) {
+ z = h(k + size_type(3), k);
+ }
+
+ }
+
+ vector tx(2);
+ tx(0) = x; tx(1) = y;
+
+ value_type scale = (std::abs)(tx(0)) + (std::abs)(tx(1)) ;
+ if (scale != value_type(0)) {
+ tx(0) /= scale;
+ tx(1) /= scale;
+ }
+
+ value_type cg, sg;
+ givens_rotation(tx(0), tx(1), cg, sg);
+
+ //Apply rotation to matrices as needed
+ matrix t1 = project(h, range(q - size_type(1), p), range(p - size_type(3), n));
+ size_type t1_cols = t1.size2();
+ size_type t1_rows = t1.size1();
+ for (size_type j = size_type(0); j < t1_cols; j++) {
+ value_type tau1 = t1(0, j);
+ value_type tau2 = t1(t1_rows - 1, j);
+ t1(0, j) = cg*tau1 - sg*tau2;
+ t1(t1_rows - 1, j) = sg*tau1 + cg*tau2;
+ }
+ project(h, range(q - size_type(1), p), range(p - size_type(3), n)).assign(t1);
+
+ matrix t2 = project(h, range(0, p), range(p - size_type(2), p));
+ size_type t2_cols = t2.size2();
+ size_type t2_rows = t2.size1();
+ for (size_type j = size_type(0); j < t2_rows; j++) {
+ value_type tau1 = t2(j, 0);
+ value_type tau2 = t2(j,t2_cols - 1);
+ t2(j, 0) = cg*tau1 - sg*tau2;
+ t2(j, t2_cols - 1) = sg*tau1 + cg*tau2;
+ }
+ project(h, range(0, p), range(p - size_type(2), p)).assign(t2);
+
+ matrix t3 = project(qv, range(0, n), range(p - size_type(2), p));
+ size_type t3_cols = t3.size2();
+ size_type t3_rows = t3.size1();
+ for (size_type j = size_type(0); j < t3_rows; j++) {
+ value_type tau1 = t3(j, 0);
+ value_type tau2 = t3(j, t3_cols - 1);
+ t3(j, 0) = cg*tau1 - sg*tau2;
+ t3(j, t3_cols - 1) = sg*tau1 + cg*tau2;
+ }
+ project(qv, range(0, n), range(p - size_type(2), p)).assign(t3);
+
+ //Pollution Cleanup
+ /*for (size_type ci = 0; ci < p; ++ci)
+ {
+ for (size_type ri = p; ri < n; ri++) {
+ h(ri, ci) = value_type(0);
+ }
+ }
+ for (size_type ci = p; ci < n - q; ci++) {
+ for (size_type ri = n - q; ri < n; ri++){
+ h(ri, ci) = value_type(0);
+ }
+ }*/
+
+ if ((std::abs)(h(p - size_type(1), q - size_type(1))) < (std::numeric_limits::epsilon())*((std::abs)(h(p - size_type(1), p - size_type(1))) + (std::abs)(h(q - size_type(1), q - size_type(1))))) {
+ h(p - size_type(1), q - size_type(1)) = value_type(0);
+ p = p - size_type(1);
+ q = p - size_type(1);
+ }
+ else if ((std::abs)(h(p - size_type(2), q - size_type(2))) < (std::numeric_limits::epsilon())*((std::abs)(h(q - size_type(2), q - size_type(2))) + (std::abs)(h(q - size_type(1), q - size_type(1))))) {
+ h(p - size_type(2), q - size_type(2)) = value_type(0);
+ p = p - size_type(2);
+ q = p - size_type(1);
+ }
+ }
+ }
+
+ template
+ void find_small_diag_entry(M &h, typename M::size_type end, typename M::value_type l1_norm, typename M::size_type &small_index)
+ {
+ typedef typename M::size_type size_type;
+ typedef typename M::value_type value_type;
+
+ size_type k = end;
+ while (k > size_type(0))
+ {
+ value_type s = (std::abs)(h(k - size_type(1), k - size_type(1))) + (std::abs)(h(k, k));
+ if (s == value_type(0))
+ s = l1_norm;
+ if ((std::abs)(h(k, k - size_type(1))) < (std::numeric_limits::epsilon() * s))
+ break;
+ --k;
+ }
+ small_index = k;
+ }
+
+ template
+ void row_split(M &h, M &qv, typename M::size_type end, typename M::value_type exceptional_shift_sum)
+ {
+ typedef typename M::size_type size_type;
+ typedef typename M::value_type value_type;
+
+ size_type n = h.size1();
+
+ value_type p = value_type(0.5) * (h(end - size_type(1), end - size_type(1)) - h(end, end));
+ value_type q = p * p + h(end, end - size_type(1)) *h(end - size_type(1), end);
+ h(end, end) += exceptional_shift_sum;
+ h(end - size_type(1), end - size_type(1)) += exceptional_shift_sum;
+
+ //Two real eigenvalues are present so separate them
+ if (q >= value_type(0))
+ {
+ value_type z = (std::sqrt)((std::abs)(q));
+ vector tx(2);
+ if (p >= value_type(0))
+ {
+ tx(0) = p + z;
+ }
+ else
+ {
+ tx(0) = p - z;
+ }
+ tx(1) = h(end, end - size_type(1));
+ value_type cg, sg;
+ givens_rotation(tx(0), tx(1), cg, sg);
+
+ //Apply rotation to matrices as needed
+ matrix t1 = project(h, range(end - size_type(1), end + size_type(1)), range(end - size_type(1), n));
+ size_type t1_cols = t1.size2();
+ size_type t1_rows = t1.size1();
+ for (size_type j = size_type(0); j < t1_cols; ++j) {
+ value_type tau1 = t1(0, j);
+ value_type tau2 = t1(t1_rows - 1, j);
+ t1(0, j) = cg*tau1 - sg*tau2;
+ t1(t1_rows - 1, j) = sg*tau1 + cg*tau2;
+ }
+ project(h, range(end - size_type(1), end + size_type(1)), range(end - size_type(1), n)).assign(t1);
+
+ matrix t2 = project(h, range(0, end + size_type(1)), range(end - size_type(1), end + size_type(1)));
+ size_type t2_cols = t2.size2();
+ size_type t2_rows = t2.size1();
+ for (size_type j = size_type(0); j < t2_rows; ++j) {
+ value_type tau1 = t2(j, 0);
+ value_type tau2 = t2(j, t2_cols - 1);
+ t2(j, 0) = cg*tau1 - sg*tau2;
+ t2(j, t2_cols - 1) = sg*tau1 + cg*tau2;
+ }
+ project(h, range(0, end + size_type(1)), range(end - size_type(1), end + size_type(1))).assign(t2);
+
+ matrix t3 = project(qv, range(0, n), range(end - size_type(1), end + size_type(1)));
+ size_type t3_cols = t3.size2();
+ size_type t3_rows = t3.size1();
+ for (size_type j = size_type(0); j < t3_rows; ++j) {
+ value_type tau1 = t3(j, 0);
+ value_type tau2 = t3(j, t3_cols - 1);
+ t3(j, 0) = cg*tau1 - sg*tau2;
+ t3(j, t3_cols - 1) = sg*tau1 + cg*tau2;
+ }
+ project(qv, range(0, n), range(end - size_type(1), end + size_type(1))).assign(t3);
+ }
+
+ if (end > size_type(1))
+ h(end - size_type(1), end - size_type(2)) = value_type(0);
+ }
+
+ template
+ void infer_shifts(M &h, typename M::size_type end, typename M::size_type iter_nos, typename M::value_type exceptional_shift_sum, vector &shift_vector)
+ {
+ typedef typename M::size_type size_type;
+ typedef typename M::value_type value_type;
+
+ shift_vector = vector(3);
+
+ shift_vector(size_type(0)) = h(end, end);
+ shift_vector(size_type(1)) = h(end - size_type(1), end - size_type(1));
+ shift_vector(size_type(2)) = h(end, end - size_type(1)) * h(end - size_type(1), end);
+
+ //Original Shift
+ if (iter_nos == size_type(10))
+ {
+ exceptional_shift_sum += shift_vector(0);
+ for (size_type i = size_type(0); i <= end; ++i)
+ {
+ h(i, i) -= shift_vector(0);
+ }
+ value_type s = (std::abs)(h(end, end - size_type(1))) + (std::abs)(h(end - size_type(1), end - size_type(2)));
+ shift_vector(0) = value_type(0.75) * s;
+ shift_vector(1) = value_type(0.75) * s;
+ shift_vector(2) = value_type(-0.4375) * s * s;
+ }
+
+ // Matlabs ad hoc shift (somehow Eigen people got this and its not in public domain)
+ // Sometimes I just think Numerical Computing is a big insider job
+ if (iter_nos == size_type(30))
+ {
+ value_type s = (shift_vector(1) - shift_vector(0))*value_type(0.5);
+ s = s * s + shift_vector(2);
+ if (s > value_type(0))
+ {
+ s = (std::sqrt)(s);
+ if (shift_vector(1) < shift_vector(0))
+ {
+ s = -s;
+ }
+ s = s + ((shift_vector(1) - shift_vector(0))*value_type(0.5));
+ s = (shift_vector(0) - shift_vector(2)) / s;
+ exceptional_shift_sum += s;
+ for (size_type i = size_type(0); i <= end; ++i)
+ h(i, i) -= s;
+ shift_vector(0) = shift_vector(1) = shift_vector(2) = value_type(0.964);
+ }
+ }
+
+ }
+
+ template
+ void francis_qr_step(M &h, M &qv, typename M::size_type rowStart, typename M::size_type end, vector shifts)
+ {
+ typedef typename M::size_type size_type;
+ typedef typename M::value_type value_type;
+
+ size_type colStart;
+ vector householderVec(3);
+
+ for (colStart = end - size_type(2); colStart >= rowStart; --colStart)
+ {
+ value_type tmm = h(colStart, colStart);
+ value_type r = shifts(0) - tmm;
+ value_type s = shifts(1) - tmm;
+ householderVec(0) = (r * s - shifts(2)) / h(colStart + size_type(1), colStart) + h(colStart, colStart + size_type(1));
+ householderVec(1) = h(colStart + size_type(1), colStart + size_type(1)) - tmm - r - s;
+ householderVec(2) = h(colStart + size_type(2), colStart + size_type(1));
+ if (colStart == rowStart) {
+ break;
+ }
+ value_type lhs = h(colStart, colStart - size_type(1)) * ((std::abs)(householderVec(1)) + (std::abs)(householderVec(2)));
+ value_type rhs = householderVec(0) * ((std::abs)(h(colStart - size_type(1), colStart - size_type(1))) + (std::abs)(tmm)+(std::abs)(h(colStart + size_type(1), colStart + size_type(1))));
+ if ((std::abs)(lhs) < (std::numeric_limits::epsilon()) * rhs)
+ {
+ break;
+ }
+ }
+
+ size_type n = h.size1();
+ bool firstTime = true;
+ for (size_type k = colStart; k <= end - size_type(2); ++k)
+ {
+ vector tx(3);
+ if (firstTime)
+ {
+ tx = householderVec;
+ firstTime = false;
+ }
+ else
+ {
+ tx(0) = h(k, k - size_type(1));
+ tx(1) = h(k + size_type(1), k - size_type(1));
+ tx(2) = h(k + size_type(1), k - size_type(1));
+ }
+ vector v;
+ value_type beta;
+ householder >(tx, v, beta);
+
+ if (beta != value_type(0))
+ {
+ if (firstTime && k > rowStart)
+ {
+ h(k, k - size_type(1)) = -h(k, k - size_type(1));
+ }
+ else if (!firstTime)
+ {
+ h(k, k - size_type(1)) = beta;
+ }
+
+ matrix vvt = outer_prod(v, v);
+ vvt *= beta;
+ size_type n_vvt = vvt.size1();
+ matrix imvvt = identity_matrix(n_vvt) - vvt;
+
+ matrix t1 = prod(imvvt, project(h, range(k, k + size_type(3)), range(k, n)));
+ project(h, range(k, k + size_type(3)), range(k, n)).assign(t1);
+
+ value_type r = (std::min)(end, k + size_type(3)) + size_type(1);
+ matrix t2 = prod(project(h, range(0, r), range(k, k + size_type(3))), imvvt);
+ project(h, range(0, r), range(k, k + size_type(3))).assign(t2);
+
+ matrix t3 = prod(project(qv, range(0, n), range(k, k + size_type(3))), imvvt);
+ project(qv, range(0, n), range(k, k + size_type(3))).assign(t3);
+
+ }
+ }
+
+ vector tx(2);
+ vector v;
+ tx(0) = h(end - size_type(1), end - size_type(2));
+ tx(1) = h(end , end - size_type(2));
+ value_type beta;
+ householder >(tx, v, beta);
+
+ if (beta != value_type(0))
+ {
+ h(end - size_type(1), end - size_type(2)) = beta;
+
+ matrix vvt = outer_prod(v, v);
+ vvt *= beta;
+ size_type n_vvt = vvt.size1();
+ matrix imvvt = identity_matrix(n_vvt) -vvt;
+
+ matrix t1 = prod(imvvt, project(h, range(end - size_type(1), end + size_type(1)), range(end - size_type(1), n)));
+ project(h, range(end - size_type(1), end + size_type(1)), range(end - size_type(1), n)).assign(t1);
+
+ matrix t2 = prod(project(h, range(0, end + size_type(1)), range(end - size_type(1), end + size_type(1))), imvvt);
+ project(h, range(0, end + size_type(1)), range(end - size_type(1), end + size_type(1))).assign(t2);
+
+ matrix t3 = prod(project(qv, range(0, n), range(end - size_type(1), end + size_type(1))), imvvt);
+ project(qv, range(0, n), range(end - size_type(1), end + size_type(1))).assign(t3);
+ }
+
+ // Round off errors (this creates a lot of issues)
+ for (size_type i = colStart + size_type(2); i <= end; ++i)
+ {
+ h(i, i - size_type(2)) = value_type(0);
+ if (i > colStart + size_type(2))
+ h(i, i - size_type(3)) = value_type(0);
+ }
+
+ }
+
+}}}
+
+
+#endif
diff --git a/test/Jamfile b/test/Jamfile
index 26b67b045..d706ad1e0 100644
--- a/test/Jamfile
+++ b/test/Jamfile
@@ -245,6 +245,11 @@ test-suite numeric/uBLAS
:
:
]
+ [ run test_eigensolver.cpp
+ :
+ :
+ :
+ ]
[ run test_fixed_containers.cpp
:
:
diff --git a/test/common/testhelper.hpp b/test/common/testhelper.hpp
index 4bc152ca3..2976d31cd 100644
--- a/test/common/testhelper.hpp
+++ b/test/common/testhelper.hpp
@@ -76,6 +76,31 @@ bool compare( const boost::numeric::ublas::matrix_expression & m1,
return true;
}
+template < class M1, class M2 >
+bool compare_on_tolerance(const boost::numeric::ublas::matrix_expression & m1,
+ const boost::numeric::ublas::matrix_expression & m2, double tolerance = 1.0e7) {
+ if ((m1().size1() != m2().size1()) ||
+ (m1().size2() != m2().size2())) {
+ return false;
+ }
+
+ size_t size1 = m1().size1();
+ size_t size2 = m1().size2();
+ for (size_t i = 0; i < size1; ++i) {
+ for (size_t j = 0; j < size2; ++j) {
+ if ((std::abs)(m1()(i, j) - m2()(i, j)) < tolerance)
+ {
+ continue;
+ }
+ else
+ {
+ return false;
+ }
+ }
+ }
+ return true;
+}
+
template < class M1, class M2 >
bool compare( const boost::numeric::ublas::vector_expression & m1,
const boost::numeric::ublas::vector_expression & m2 ) {
diff --git a/test/test_eigensolver.cpp b/test/test_eigensolver.cpp
new file mode 100644
index 000000000..a4a33ce4a
--- /dev/null
+++ b/test/test_eigensolver.cpp
@@ -0,0 +1,141 @@
+// Copyright 2008 Gunter Winkler
+// Distributed under the Boost Software License, Version 1.0. (See
+// accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+
+// switch automatic singular check off
+#define BOOST_UBLAS_TYPE_CHECK 0
+
+#include
+#include
+#include
+
+#include "common/testhelper.hpp"
+
+#include
+#include
+
+using namespace boost::numeric::ublas;
+using std::string;
+
+static const string matrix_IN[] = { "[3,3]((-149.0,-50.0,-154.0),(537.0,180.0,546.0),(-27.0,-9.0,-25.0))\0",
+ "[3,3]((2.0,3.0,5.0),(2.0,-3.0,7.0),(4.0,1.0,1.0))\0",
+ "[6,6]((7.0, 3.0, 4.0, -11.0, -9.0, -2.0),(-6.0, 4.0, -5.0, 7.0, 1.0, 12.0),(-1.0, -9.0, 2.0, 2.0, 9.0, 1.0),(-8.0, 0.0, -1.0, 5.0, 0.0, 8.0),(-4.0, 3.0, -5.0, 7.0, 2.0, 10.0),(6.0, 1.0, 4.0, -11.0, -7.0, -1.0))\0",
+ "[4,4]((0.421761282626275,0.655740699156587,0.678735154857774,0.655477890177557),(0.915735525189067,0.0357116785741896,0.757740130578333,0.171186687811562),(0.792207329559554,0.849129305868777,0.743132468124916,0.706046088019609),(0.959492426392903, 0.933993247757551, 0.392227019534168, 0.0318328463774207))\0",
+ "[10,10]((1.0, 2.0, 3.0, 4.0, 3.0, 6.0, 7.0, 8.0, 9.0, 10.0),(1.0, 2.0, 3.0, 3.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0),(3.0, 2.0, 3.0, 4.2, 5.0, 6.1, 7.0, 8.0, 9.0, 14.0),(1.0, 2.1, 3.0, 4.0, 5.0, 6.0, 7.1, 8.0, 9.0, 10.0),(1.0, 2.0, 3.0, 4.5, 5.0, 6.0, 7.0, 8.0, 9.0, 30.0),(5.0, 6.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0),(1.0, 2.0, 3.0, 41.0, 5.0, 6.0, 7.0, 1.0, 9.0, 10.0),(1.0, 2.0, 32.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.3, 10.0),(9.0, 2.0, 3.0, 4.0, 0.0, 6.0, 7.0, 8.0, 9.0, 10.0),(12.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0))\0",
+ "[10,10]((0, 2, 3, 4, 3, 6, 7, 8, 9, 10),(1, 2, 3, 3, 5, 0, 7, 8, 9, 10),(3, 2, 0, 4.2, 5, 6.1, 0, 8, 9, 0),(1, 2, 3, 4, 5, 6, 0, 8, 9, 1),(1, 2, 3, 4.5, 5, 6, 7, 8, 9, 0),(5, 6, 0, 4, 5, 6, 7, 8, 9, 1),(0, 2, 3, 0, 5, 6, 7, 1, 9, 4),(1, 2, 0, 4, 5, 6, 7, 8, 9.3, 0),(0, 2, 3, 4, 0, 6, 7, 8, 9, 10),(0, 2, 3, 4, 5, 6, 7, 8, 9, 0))\0",
+
+ "[10,10]((0, 2e-05, 3e-05, 4e-05, 3e-05, 6e-05, 7e-05, 8e-05, 9e-05, 0.0001),(1e-05, 2e-05, 3e-05, 3e-05, 5e-05, 0, 7e-05, 8e-05, 9e-05, 0.0001),(3e-05, 2e-05, 0, 4.2e-05, 5e-05, 6.1e-05, 0, 8e-05, 9e-05, 0),(1e-05, 2e-05, 3e-05, 4e-05, 5e-05, 6e-05, 0, 8e-05, 9e-05, 1e-05),(1e-05, 2e-05, 3e-05, 4.5e-05, 5e-05, 6e-05, 7e-05, 8e-05, 9e-05, 0),(5e-05, 6e-05, 0, 4e-05, 5e-05, 6e-05, 7e-05, 8e-05, 9e-05, 1e-05),(0, 2e-05, 3e-05, 0, 5e-05, 6e-05, 7e-05, 1e-05, 9e-05, 4e-05),(1e-05, 2e-05, 0, 4e-05, 5e-05, 6e-05, 7e-05, 8e-05, 9.3e-05, 0),(0, 2e-05, 3e-05, 4e-05, 0, 6e-05, 7e-05, 8e-05, 9e-05, 0.0001),(0, 2e-05, 3e-05, 4e-05, 5e-05, 6e-05, 7e-05, 8e-05, 9e-05, 0))\0",
+ "[10,10]((4, 2, 3, 4, 3, 6, 7, 8, 9, 10),(0, 2, 3, 3, 5, 0, 7, 8, 9, 10),(3, 2, 0, 4.2, 5, 6.1, 0, 8, 9, 0),(1, 2, 3, 4, 5, 6, 0, 8, 9, 1),(3, 2, 3, 4.5, 5, 6, 0, 8, 9, 0),(5, 2, 0, 4, 5, 6, 1, 8, 9, 1),(0, 2, 3, 0, 5, 6, 7, 0, 9, 4),(1, 13, 0, 4, 5, 6, 2, 8, 9.3, 0),(100, 2, 3, 0, 0, 6, 7, 8, 9, 10),(0, 2, 3, 4, 5, 6, 0, 8, 9, 0),)\0",
+ "[10,10]((0, 2e+10, 3e+10, 4e+10, 3e+10, 6e+10, 7e+10, 8e+10, 9e+10, 1e+11),(1e+10, 2e+10, 3e+10, 3e+10, 5e+10, 0, 7e+10, 8e+10, 9e+10, 1e+11),(3e+10, 2e+10, 0, 4.2e+10, 5e+10, 6.1e+10, 0, 8e+10, 9e+10, 0),(1e+10, 2e+10, 3e+10, 4e+10, 5e+10, 6e+10, 0, 8e+10, 9e+10, 1e+10),(1e+10, 2e+10, 3e+10, 4.5e+10, 5e+10, 6e+10, 7e+10, 8e+10, 9e+10, 0),(5e+10, 6e+10, 0, 4e+10, 5e+10, 6e+10, 7e+10, 8e+10, 9e+10, 1e+10),(0, 2e+10, 3e+10, 0, 5e+10, 6e+10, 7e+10, 1e+10, 9e+10, 4e+10),(1e+10, 2e+10, 0, 4e+10, 5e+10, 6e+10, 7e+10, 8e+10, 9.3e+10, 0),(0, 2e+10, 3e+10, 4e+10, 0, 6e+10, 7e+10, 8e+10, 9e+10, 1e+11),(0, 2e+10, 3e+10, 4e+10, 5e+10, 6e+10, 7e+10, 8e+10, 9e+10, 0))\0"
+ };
+
+static const string matrix_EVALR[] = { "[3,3]((1,0.0,0.0),(0.0,2,0.0),(0.0,0.0,3))\0",
+ "[3,3]((7.547183,0.0,0.0),(0.0,-3.773591,0.0),(0.0,0.0,-3.773591))\0",
+ "[6,6]((5.0,0.0,0.0,0.0,0.0,0.0),(0.0,5.0,0.0,0.0,0.0,0.0),(0.0,0.0,1.0,0.0,0.0,0.0),(0.0,0.0,0.0,1.0,0.0,0.0),(0.0,0.0,0.0,0.0,4.0,0.0),(0.0,0.0,0.0,0.0,0.0,3.0))\0",
+ "[4,4]((2.44784,0.00000,0.00000,0.00000),(0.00000,-0.56040,0.00000,0.00000),(0.00000,0.00000,-0.56040,0.00000),(0.00000,0.00000,0.00000,-0.09461))\0",
+ "[10,10]((66.1927, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000), (0.000000, -10.9351, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000), (0.000000, 0.000000, -7.02399, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000), (0.000000, 0.000000, 0.000000, 2.03077, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000), (0.000000, 0.000000, 0.000000, 0.000000, 2.03077, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000), (0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 3.56282, 0.000000, 0.000000, 0.000000, 0.000000), (0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, -0.593575, 0.000000, 0.000000, 0.000000), (0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, -0.593575, 0.000000, 0.000000), (0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, -0.237288, 0.000000), (0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.566443))\0",
+ "[10,10]((44.5477, 0, 0, 0, 0, 0, 0, 0, 0, 0),(0, -1.63051, 0, 0, 0, 0, 0, 0, 0, 0),(0, 0, -1.63051, 0, 0, 0, 0, 0, 0, 0),(0, 0, 0, -3.9098, 0, 0, 0, 0, 0, 0),(0, 0, 0, 0, -1.51227, 0, 0, 0, 0, 0),(0, 0, 0, 0, 0, -1.51227, 0, 0, 0, 0),(0, 0, 0, 0, 0, 0, 2.59252, 0, 0, 0),(0, 0, 0, 0, 0, 0, 0, 2.59252, 0, 0),(0, 0, 0, 0, 0, 0, 0, 0, 0.731304, 0),(0, 0, 0, 0, 0, 0, 0, 0, 0, 0.731304))\0",
+ "[10,10]((0.000445477, 0, 0, 0, 0, 0, 0, 0, 0, 0),(0, -1.63051e-05, 0, 0, 0, 0, 0, 0, 0, 0),(0, 0, -1.63051e-05, 0, 0, 0, 0, 0, 0, 0),(0, 0, 0, -3.9098e-05, 0, 0, 0, 0, 0, 0),(0, 0, 0, 0, -1.51227e-05, 0, 0, 0, 0, 0),(0, 0, 0, 0, 0, -1.51227e-05, 0, 0, 0, 0),(0, 0, 0, 0, 0, 0, 2.59252e-05, 0, 0, 0),(0, 0, 0, 0, 0, 0, 0, 2.59252e-05, 0, 0),(0, 0, 0, 0, 0, 0, 0, 0, 7.31304e-06, 0),(0, 0, 0, 0, 0, 0, 0, 0, 0, 7.31304e-06))\0",
+ "[10,10]((60.8688, 0, 0, 0, 0, 0, 0, 0, 0, 0),(0, -9.82873, 0, 0, 0, 0, 0, 0, 0, 0),(0, 0, -9.82873, 0, 0, 0, 0, 0, 0, 0),(0, 0, 0, -5.64923, 0, 0, 0, 0, 0, 0),(0, 0, 0, 0, 4.19257, 0, 0, 0, 0, 0),(0, 0, 0, 0, 0, 4.19257, 0, 0, 0, 0),(0, 0, 0, 0, 0, 0, -0.578372, 0, 0, 0),(0, 0, 0, 0, 0, 0, 0, -0.578372, 0, 0),(0, 0, 0, 0, 0, 0, 0, 0, 1.10473, 0),(0, 0, 0, 0, 0, 0, 0, 0, 0, 1.10473))\0",
+ "[10,10]((4.45477e+11, 0, 0, 0, 0, 0, 0, 0, 0, 0),(0, -1.63051e+10, 0, 0, 0, 0, 0, 0, 0, 0),(0, 0, -1.63051e+10, 0, 0, 0, 0, 0, 0, 0),(0, 0, 0, -3.9098e+10, 0, 0, 0, 0, 0, 0),(0, 0, 0, 0, -1.51227e+10, 0, 0, 0, 0, 0),(0, 0, 0, 0, 0, -1.51227e+10, 0, 0, 0, 0),(0, 0, 0, 0, 0, 0, 2.59252e+10, 0, 0, 0),(0, 0, 0, 0, 0, 0, 0, 2.59252e+10, 0, 0),(0, 0, 0, 0, 0, 0, 0, 0, 7.31304e+09, 0),(0, 0, 0, 0, 0, 0, 0, 0, 0, 7.31304e+09))\0"
+ };
+
+static const string matrix_EVALC[] = { "[3,3]((0.0,0.0,0.0),(0.0,0.0,0.0),(0.0,0.0,0.0))\0",
+ "[3,3]((0.0,0.0,0.0),(0.0,1.649236,0.0),(0.0,0.0,-1.649236))\0",
+ "[6,6]((6.0,0.0,0.0,0.0,0.0,0.0),(0.0,-6.0,0.0,0.0,0.0,0.0),(0.0,0.0,2.0,0.0,0.0,0.0),(0.0,0.0,0.0,-2.0,0.0,0.0),(0.0,0.0,0.0,0.0,0.0,0.0),(0.0,0.0,0.0,0.0,0.0,0.0))\0",
+ "[4,4]((0.00000,0.00000,0.00000,0.00000),(0.00000,0.31770,0.00000,0.00000),(0.00000,0.00000,-0.31770,0.00000),(0.00000,0.00000,0.00000,0.00000))\0",
+ "[10,10]((0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000), (0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000), (0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000), (0.000000, 0.000000, 0.000000, 4.66122, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000), (0.000000, 0.000000, 0.000000, 0.000000, -4.66122, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000), (0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000), (0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.884104, 0.000000, 0.000000, 0.000000), (0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, -0.884104, 0.000000, 0.000000), (0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000), (0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000))\0",
+ "[10,10]((0, 0, 0, 0, 0, 0, 0, 0, 0, 0),(0, 5.34552, 0, 0, 0, 0, 0, 0, 0, 0),(0, 0, -5.34552, 0, 0, 0, 0, 0, 0, 0),(0, 0, 0, 0, 0, 0, 0, 0, 0, 0),(0, 0, 0, 0, 2.58278, 0, 0, 0, 0, 0),(0, 0, 0, 0, 0, -2.58278, 0, 0, 0, 0),(0, 0, 0, 0, 0, 0, 1.58579, 0, 0, 0),(0, 0, 0, 0, 0, 0, 0, -1.58579, 0, 0),(0, 0, 0, 0, 0, 0, 0, 0, 0.577581, 0),(0, 0, 0, 0, 0, 0, 0, 0, 0, -0.577581))\0",
+ "[10,10]((0, 0, 0, 0, 0, 0, 0, 0, 0, 0),(0, 5.34552e-05, 0, 0, 0, 0, 0, 0, 0, 0),(0, 0, -5.34552e-05, 0, 0, 0, 0, 0, 0, 0),(0, 0, 0, 0, 0, 0, 0, 0, 0, 0),(0, 0, 0, 0, 2.58278e-05, 0, 0, 0, 0, 0),(0, 0, 0, 0, 0, -2.58278e-05, 0, 0, 0, 0),(0, 0, 0, 0, 0, 0, 1.58579e-05, 0, 0, 0),(0, 0, 0, 0, 0, 0, 0, -1.58579e-05, 0, 0),(0, 0, 0, 0, 0, 0, 0, 0, 5.77581e-06, 0),(0, 0, 0, 0, 0, 0, 0, 0, 0, -5.77581e-06))\0",
+ "[10,10]((0, 0, 0, 0, 0, 0, 0, 0, 0, 0),(0, 11.8224, 0, 0, 0, 0, 0, 0, 0, 0),(0, 0, -11.8224, 0, 0, 0, 0, 0, 0, 0),(0, 0, 0, 0, 0, 0, 0, 0, 0, 0),(0, 0, 0, 0, 3.47865, 0, 0, 0, 0, 0),(0, 0, 0, 0, 0, -3.47865, 0, 0, 0, 0),(0, 0, 0, 0, 0, 0, 0.5707, 0, 0, 0),(0, 0, 0, 0, 0, 0, 0, -0.5707, 0, 0),(0, 0, 0, 0, 0, 0, 0, 0, 0.0717088, 0),(0, 0, 0, 0, 0, 0, 0, 0, 0, -0.0717088))\0",
+ "[10,10]((0, 0, 0, 0, 0, 0, 0, 0, 0, 0),(0, 5.34552e+10, 0, 0, 0, 0, 0, 0, 0, 0),(0, 0, -5.34552e+10, 0, 0, 0, 0, 0, 0, 0),(0, 0, 0, 0, 0, 0, 0, 0, 0, 0),(0, 0, 0, 0, 2.58278e+10, 0, 0, 0, 0, 0),(0, 0, 0, 0, 0, -2.58278e+10, 0, 0, 0, 0),(0, 0, 0, 0, 0, 0, 1.58579e+10, 0, 0, 0),(0, 0, 0, 0, 0, 0, 0, -1.58579e+10, 0, 0),(0, 0, 0, 0, 0, 0, 0, 0, 5.77581e+09, 0),(0, 0, 0, 0, 0, 0, 0, 0, 0, -5.77581e+09))\0"
+ };
+
+
+int main() {
+
+ typedef double TYPE;
+ typedef float TYPE2;
+
+ typedef matrix MATRIX;
+
+ MATRIX A;
+ MATRIX EVALR;
+ MATRIX EVALC;
+ MATRIX EVECR;
+ MATRIX EVECC;
+ MATRIX GT_EVALR;
+ MATRIX GT_EVALC;
+ MATRIX Zero_Matrix = zero_matrix(3,3);
+
+ int numTestCases = 8;
+
+
+ for (int k = 0; k < numTestCases; k++){
+
+ std::cout << "Running test case# " << (k + 1) << "\n";
+
+ {
+ std::istringstream is(matrix_IN[k]);
+ is >> A;
+ }
+
+ {
+ std::istringstream is(matrix_EVALR[k]);
+ is >> GT_EVALR;
+ }
+
+ {
+ std::istringstream is(matrix_EVALC[k]);
+ is >> GT_EVALC;
+ }
+
+
+ eigen_solver es(A, EIGVEC);
+ matrix evals_r = es.get_real_eigenvalues();
+ matrix evals_c = es.get_complex_eigenvalues();
+ matrix evecs_r = es.get_real_eigenvectors();
+ matrix evecs_c = es.get_complex_eigenvectors();
+
+
+ assertTrue("Real portion of Eigen Values Match:", compare_on_tolerance(evals_r, GT_EVALR));
+ assertTrue("Imag. portion of Eigen Values Match:", compare_on_tolerance(evals_c, GT_EVALC));
+
+ matrix > V(3, 3);
+ matrix > D(3, 3);
+ matrix > M(3, 3);
+ matrix > Lambda;
+
+ MATRIX Lambda_Real(3, 3);
+ MATRIX Lambda_Complex(3, 3);
+
+ for (int i = 0; i < 3; i++){
+ for (int j = 0; j < 3; j++){
+ V(i, j) = std::complex(evecs_r(i, j), evecs_c(i, j));
+ D(i, j) = std::complex(evals_r(i, j), evals_c(i, j));
+ M(i, j) = std::complex(A(i, j), 0.0);
+ }
+ }
+ Lambda = prod(M, V) - prod(V, D);
+ for (int i = 0; i < 3; i++){
+ for (int j = 0; j < 3; j++){
+ Lambda_Real(i, j) = Lambda(i, j).real();
+ Lambda_Complex(i, j) = Lambda(i, j).imag();
+ }
+ }
+
+ assertTrue("Real portion of verifier is zero:", compare_on_tolerance(Lambda_Real, Zero_Matrix));
+ assertTrue("Imag. portion of verifier is zero:", compare_on_tolerance(Lambda_Complex, Zero_Matrix));
+
+ }
+
+ if (!_fail_counter) {
+ std::cout << "\nEigen Solver regression suite passed.\n";
+ return boost::exit_success;
+ }
+ else {
+ std::cout << "\nEigen Solver regression suite failed.\n";
+ return boost::exit_failure;
+ }
+}
Boost uBLAS EigenSolver