From cbdfd8fe515e0b0edffb869f763317f89c1ec3f9 Mon Sep 17 00:00:00 2001 From: rajaditya-m Date: Wed, 3 Jun 2015 20:33:46 -0400 Subject: [PATCH 01/20] Hessenberg The Hessenberg decomposition module for square matrix is completed... --- include/boost/numeric/ublas/hessenberg.hpp | 23 +++++++++++++++++++--- 1 file changed, 20 insertions(+), 3 deletions(-) diff --git a/include/boost/numeric/ublas/hessenberg.hpp b/include/boost/numeric/ublas/hessenberg.hpp index 05c73ad9b..ca95831e2 100644 --- a/include/boost/numeric/ublas/hessenberg.hpp +++ b/include/boost/numeric/ublas/hessenberg.hpp @@ -25,9 +25,26 @@ namespace boost { namespace numeric { namespace ublas { 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; - vector temp(n); - for (size_type i = 0; i < n; ++i) { - //Insert householder modules here + 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); + + //@TODO: Ask mentor why uncommented lines are more accurate than commented lines... + //project(m, range(i + 1, num_rows), range(i, num_rows)).assign(prod(imvvt, project(m, range(i + 1, num_rows), range(i, num_rows)))); + //project(m, range(0, num_rows), range(i+1, num_rows)).assign(prod(project(m, range(0, num_rows), range(i+1, num_rows)),imvvt)); + } } From 060f43706eb231eabb87ebb20bc4e341f32c6fad Mon Sep 17 00:00:00 2001 From: rajaditya-m Date: Wed, 10 Jun 2015 13:19:53 +0530 Subject: [PATCH 02/20] Francis QR Step --- .../numeric/ublas/schur_decomposition.hpp | 94 +++++++++++++++++++ 1 file changed, 94 insertions(+) create mode 100644 include/boost/numeric/ublas/schur_decomposition.hpp diff --git a/include/boost/numeric/ublas/schur_decomposition.hpp b/include/boost/numeric/ublas/schur_decomposition.hpp new file mode 100644 index 000000000..2bf220e09 --- /dev/null +++ b/include/boost/numeric/ublas/schur_decomposition.hpp @@ -0,0 +1,94 @@ +// 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) + +#ifndef _BOOST_UBLAS_SCHURDECOMPOSITION_ +#define _BOOST_UBLAS_SCHURDECOMPOSITION_ + +#include +#include +#include +#include +#include + +namespace boost { + namespace numeric { + namespace ublas { + + 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 m = n - size_type(1); + + value_type s = h(m, m) + h(n, n); + value_type t = h(m, m) * h(n, n) - h(m, n) * h(n, m); + 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 <= n - 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 q = (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(q - size_type(1),n - size_type(1)))); + project(h, range(k,k + size_type(3)), range(q - size_type(1), n - size_type(1))).assign(t1); + + size_type r = (std::min)(k + size_type(4), n); + + 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 < (n - size_type(3))) { + z = h(k + size_type(3), k); + } + + } + + vector tx(2); + tx(0) = x; tx(1) = y; + vector v; + value_type beta; + householder >(tx, 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(h, range(n - size_type(2),n), range(n - size_type(3), n))); + project(h, range(n - size_type(2) ,n), range(n - size_type(3), n)).assign(t1); + + matrix t2 = prod(project(h, range(0,n), range(n - size_type(2), n)),imvvt); + project(h, range(0,n), range(n - size_type(2), n)).assign(t2); + + + } + + + + +}}} + + +#endif From 6c8137ad496d179907611e8736b79ae20627af99 Mon Sep 17 00:00:00 2001 From: rajaditya-m Date: Mon, 15 Jun 2015 10:56:10 +0530 Subject: [PATCH 03/20] FRancis Double Step QR This finished the recovery of eigen values --- include/boost/numeric/ublas/hessenberg.hpp | 4 +- .../numeric/ublas/schur_decomposition.hpp | 112 ++++++++++-------- 2 files changed, 66 insertions(+), 50 deletions(-) diff --git a/include/boost/numeric/ublas/hessenberg.hpp b/include/boost/numeric/ublas/hessenberg.hpp index ca95831e2..40031ff17 100644 --- a/include/boost/numeric/ublas/hessenberg.hpp +++ b/include/boost/numeric/ublas/hessenberg.hpp @@ -42,8 +42,8 @@ namespace boost { namespace numeric { namespace ublas { project(m, range(0, n), range(i + 1, n)).assign(t2); //@TODO: Ask mentor why uncommented lines are more accurate than commented lines... - //project(m, range(i + 1, num_rows), range(i, num_rows)).assign(prod(imvvt, project(m, range(i + 1, num_rows), range(i, num_rows)))); - //project(m, range(0, num_rows), range(i+1, num_rows)).assign(prod(project(m, range(0, num_rows), range(i+1, num_rows)),imvvt)); + //project(m, range(i + 1, n), range(i, n)).assign(prod(imvvt, project(m, range(i + 1, n), range(i, n)))); + //project(m, range(0, n), range(i+1, n)).assign(prod(project(m, range(0, n), range(i+1, n)),imvvt)); } diff --git a/include/boost/numeric/ublas/schur_decomposition.hpp b/include/boost/numeric/ublas/schur_decomposition.hpp index 2bf220e09..7bdfd97c1 100644 --- a/include/boost/numeric/ublas/schur_decomposition.hpp +++ b/include/boost/numeric/ublas/schur_decomposition.hpp @@ -5,11 +5,14 @@ #ifndef _BOOST_UBLAS_SCHURDECOMPOSITION_ #define _BOOST_UBLAS_SCHURDECOMPOSITION_ +#define EPSA 1.0e-20 + #include #include #include #include #include +#include namespace boost { namespace numeric { @@ -22,67 +25,80 @@ namespace boost { typedef typename M::value_type value_type; size_type n = h.size1(); - size_type m = n - size_type(1); - - value_type s = h(m, m) + h(n, n); - value_type t = h(m, m) * h(n, n) - h(m, n) * h(n, m); - 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 <= n - 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 p = n; + while (p > size_type(2)) { + size_type q = p - size_type(1); - size_type q = (std::max)(size_type(1), k); + 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)); - matrix vvt = outer_prod(v, v); - vvt *= beta; - size_type n_vvt = vvt.size1(); - matrix imvvt = identity_matrix(n_vvt) - vvt; + //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); - matrix t1 = prod(imvvt, project(h, range(k,k + size_type(3)), range(q - size_type(1),n - size_type(1)))); - project(h, range(k,k + size_type(3)), range(q - size_type(1), n - size_type(1))).assign(t1); + size_type r = (std::max)(size_type(1), k); - size_type r = (std::min)(k + size_type(4), n); + matrix vvt = outer_prod(v, v); + vvt *= beta; + size_type n_vvt = vvt.size1(); + matrix imvvt = identity_matrix(n_vvt) -vvt; - 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 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); - x = h(k + size_type(1), k); - y = h(k + size_type(2), k); + r = (std::min)(k + size_type(4), p); - if (k < (n - size_type(3))) { - z = h(k + size_type(3), k); - } + 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); - vector tx(2); - tx(0) = x; tx(1) = y; - vector v; - value_type beta; - householder >(tx, v, beta); + if (k < (p - size_type(3))) { + z = h(k + size_type(3), 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(n - size_type(2),n), range(n - size_type(3), n))); - project(h, range(n - size_type(2) ,n), range(n - size_type(3), n)).assign(t1); - - matrix t2 = prod(project(h, range(0,n), range(n - size_type(2), n)),imvvt); - project(h, range(0,n), range(n - size_type(2), n)).assign(t2); + vector tx(2); + tx(0) = x; tx(1) = y; + vector v; + value_type beta; + householder >(tx, 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(h, range(q - size_type(1), p), range(p - size_type(3), n))); + project(h, range(q - size_type(1), p), range(p - size_type(3), n)).assign(t1); + matrix t2 = prod(project(h, range(0, p), range(p - size_type(2), p)), imvvt); + 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))) < (EPSA)*((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))) < (EPSA)*((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); + } + std::cout << p << " " << q << "\n"; + } } From b3dbc13bc3e576f103257bb9cbfe6c9bc92fe9f8 Mon Sep 17 00:00:00 2001 From: rajaditya-m Date: Thu, 18 Jun 2015 23:28:59 +0530 Subject: [PATCH 04/20] Transform Accumulations This is needed for the computation of Eigen Vectors --- include/boost/numeric/ublas/hessenberg.hpp | 48 ++++++ .../numeric/ublas/schur_decomposition.hpp | 143 +++++++++++++++++- 2 files changed, 188 insertions(+), 3 deletions(-) diff --git a/include/boost/numeric/ublas/hessenberg.hpp b/include/boost/numeric/ublas/hessenberg.hpp index 40031ff17..d6ff29911 100644 --- a/include/boost/numeric/ublas/hessenberg.hpp +++ b/include/boost/numeric/ublas/hessenberg.hpp @@ -10,6 +10,9 @@ #include #include #include +#include + +#include namespace boost { namespace numeric { namespace ublas { @@ -49,6 +52,51 @@ namespace boost { namespace numeric { namespace ublas { } + template + typename M::size_type 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); + + //std::cout << i << " " << x << " " << v << "\n"; + + 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); + + //@TODO: Ask mentor why uncommented lines are more accurate than commented lines... + //project(m, range(i + 1, n), range(i, n)).assign(prod(imvvt, project(m, range(i + 1, n), range(i, n)))); + //project(m, range(0, n), range(i+1, n)).assign(prod(project(m, range(0, n), range(i+1, n)),imvvt)); + + } + } }}} diff --git a/include/boost/numeric/ublas/schur_decomposition.hpp b/include/boost/numeric/ublas/schur_decomposition.hpp index 7bdfd97c1..a9c39a731 100644 --- a/include/boost/numeric/ublas/schur_decomposition.hpp +++ b/include/boost/numeric/ublas/schur_decomposition.hpp @@ -100,10 +100,147 @@ namespace boost { std::cout << p << " " << q << "\n"; } } + + 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; + 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; + vector v; + value_type beta; + householder >(tx, 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(h, range(q - size_type(1), p), range(p - size_type(3), n))); + project(h, range(q - size_type(1), p), range(p - size_type(3), n)).assign(t1); + + matrix t2 = prod(project(h, range(0, p), range(p - size_type(2), p)), imvvt); + project(h, range(0, p), range(p - size_type(2), p)).assign(t2); + + matrix t3 = prod(project(qv, range(0,n), range(p - size_type(2),p)), imvvt); + project(qv, range(0, p), range(p - size_type(2), p)).assign(t3); + + if ((std::abs)(h(p - size_type(1), q - size_type(1))) < (EPSA)*((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))) < (EPSA)*((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); + } + //std::cout << p << " " << q << "\n"; + } + } + - - - + template + void inverse_iterations(M &h, M &u0, M &v) { + typedef typename M::size_type size_type; + typedef typename M::value_type value_type; + + size_type n = h.size1(); + + v = identity_matrix (n); + matrix in = identity_matrix(n); + //This is the main loop + for (size_type i = size_type(0); i < n; i++) { + value_type mu = h(i, i); + matrix mu_i = mu * in; + matrix mat_lhs = h - mu_i; + size_type max_iters = size_type(1000); //This should be adjustable + vector q = column(v, i); + vector z; + int counter = 0; + //This is the iteration loop + bool convergence_reached = false; + while (!convergence_reached) { + //Start solve of mat_lhs * z = q via column oriented back-substitution + //vector q_copy = q; + for (size_type j = n; j >= size_type(2); j--) { + q(j - size_type(1)) = q(j - size_type(1)) / mat_lhs(j - size_type(1), j - size_type(1)); + vector u_ranged = project(column(mat_lhs, j - size_type(1)), range(size_type(0), j - size_type(1))); + u_ranged *= q(j - size_type(1)); + vector b_ranged = project(q, range(size_type(0), j - size_type(1))); + b_ranged -= u_ranged; + project(q, range(size_type(0), j - size_type(1))).assign(b_ranged); + } + q(size_type(0)) = q(size_type(0)) / mat_lhs(size_type(0), size_type(0)); + + //Normalize the q + value_type q_norm = norm_2(q); + //q = q / q_norm; + + counter++; + std::cout << counter << " " << q << "\n"; + if (counter==10) + convergence_reached = true; + } + + column(v, i) = q; + + + } + + + } }}} From 6ac85ee4a5a095a2bbba1462968153c8d205e74c Mon Sep 17 00:00:00 2001 From: rajaditya-m Date: Tue, 23 Jun 2015 13:19:53 +0530 Subject: [PATCH 05/20] Eigen Solver Wrapper eigen_solver.hpp which can compute real eigenvalyes of real matrices.. --- include/boost/numeric/ublas/eigen_solver.hpp | 97 +++++++++++++++++++ .../numeric/ublas/schur_decomposition.hpp | 32 ++++++ 2 files changed, 129 insertions(+) create mode 100644 include/boost/numeric/ublas/eigen_solver.hpp diff --git a/include/boost/numeric/ublas/eigen_solver.hpp b/include/boost/numeric/ublas/eigen_solver.hpp new file mode 100644 index 000000000..30278643c --- /dev/null +++ b/include/boost/numeric/ublas/eigen_solver.hpp @@ -0,0 +1,97 @@ +// 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) + +#ifndef _BOOST_UBLAS_EIGENSOLVER_ +#define _BOOST_UBLAS_EIGENSOLVER_ + +#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 +typedef enum eig_solver_params {EIGVAL,EIGVEC}; + +template +class eigen_solver { +private: + M matrix; + M hessenberg_form; + M real_schur_form; + M transform_accumulations; + bool has_complex; + bool has_eigenvalues; + bool has_eigenvectors; + M eigenvalues; + M eigenvectors; + eig_solver_params solver_params; +public: + //This is the constructor + BOOST_UBLAS_INLINE + explicit eigen_solver(M &m, eig_solver_params params = EIGVAL) : + matrix(m),has_complex(false),solver_params(params) + { + has_eigenvalues = has_eigenvectors = false; + hessenberg_form = M(m); + compute(solver_params); + } + + //The compute function + 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; + } + else { + to_hessenberg(hessenberg_form, transform_accumulations); + real_schur_form = M(hessenberg_form); + schur_decomposition(real_schur_form, transform_accumulations); + extract_eigenvalues_from_schur(); + eigenvectors = M(transform_accumulations); + M temp_q_storage = M(real_schur_form); + block_diag(temp_q_storage, eigenvectors); + has_eigenvectors = true; + } + + + } + + //This is a tricky function because it will do many things... + //Right now it doesn't + BOOST_UBLAS_INLINE + void extract_eigenvalues_from_schur() { + eigenvalues = diagonal_adaptor(real_schur_form); + has_eigenvalues = true; + } + + //Extraction function + BOOST_UBLAS_INLINE + M& get_real_eigenvalues() { + return eigenvalues; + } + + BOOST_UBLAS_INLINE + M& get_real_eigenvectors() { + return eigenvectors; + } + + +}; + + +}}} + + +#endif diff --git a/include/boost/numeric/ublas/schur_decomposition.hpp b/include/boost/numeric/ublas/schur_decomposition.hpp index a9c39a731..9f1a250d4 100644 --- a/include/boost/numeric/ublas/schur_decomposition.hpp +++ b/include/boost/numeric/ublas/schur_decomposition.hpp @@ -241,6 +241,38 @@ namespace boost { } + + template + void block_diag(M &t, M &q) { + typedef typename M::size_type size_type; + typedef typename M::value_type value_type; + + size_type n = t.size1(); + + for (size_type j = size_type(2); j <= n; j++) { + for (size_type i = size_type(1); i < j; i++) { + value_type tii = t(i - size_type(1), i - size_type(1)); + value_type tjj = t(j - size_type(1), j - size_type(1)); + value_type tij = t(i - size_type(1), j - size_type(1)); + value_type z = -tij / (tii - tjj); + for (size_type k = j + size_type(1); k <= n; k++) { + t(i - size_type(1), k - size_type(1)) -= (z*t(j - size_type(1), k - size_type(1))); + } + for (size_type k = size_type(1); k <= n; k++) { + q(k - size_type(1), j - size_type(1)) += (z*q(k - size_type(1), i - size_type(1))); + } + } + } + + //We need to normalize the columns also + for (size_type i = 0; i < n; i++) { + vector vi = column(q, i); + value_type norm_vi = norm_2(vi); + vector normalized_vi = vi / norm_vi; + column(q, i) = normalized_vi; + } + + } }}} From 62872a857b003fd9b14d33831a5cf0b90b0de9e1 Mon Sep 17 00:00:00 2001 From: rajaditya-m Date: Tue, 7 Jul 2015 13:17:23 +0530 Subject: [PATCH 06/20] Givens Rotation and Modified Schur Decomposition Implemented Givens rotation and modified schur decomposition for better convergence in complex case --- include/boost/numeric/ublas/eigen_solver.hpp | 12 ++ include/boost/numeric/ublas/householder.hpp | 28 +++++ .../numeric/ublas/schur_decomposition.hpp | 103 +++++++++++++++--- 3 files changed, 129 insertions(+), 14 deletions(-) diff --git a/include/boost/numeric/ublas/eigen_solver.hpp b/include/boost/numeric/ublas/eigen_solver.hpp index 30278643c..fb4a80437 100644 --- a/include/boost/numeric/ublas/eigen_solver.hpp +++ b/include/boost/numeric/ublas/eigen_solver.hpp @@ -88,6 +88,18 @@ class eigen_solver { } + //THis is my functions - just kept for debugging will be deleted (or maybe not) + BOOST_UBLAS_INLINE + M& get_real_schur_form() { + return real_schur_form; + } + + BOOST_UBLAS_INLINE + M& get_hessenberg_form() { + return hessenberg_form; + } + + }; diff --git a/include/boost/numeric/ublas/householder.hpp b/include/boost/numeric/ublas/householder.hpp index 619b9bc29..c6e9a64ad 100644 --- a/include/boost/numeric/ublas/householder.hpp +++ b/include/boost/numeric/ublas/householder.hpp @@ -57,6 +57,34 @@ namespace boost { return size_type(0); } + + 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/schur_decomposition.hpp b/include/boost/numeric/ublas/schur_decomposition.hpp index 9f1a250d4..d24e1ad96 100644 --- a/include/boost/numeric/ublas/schur_decomposition.hpp +++ b/include/boost/numeric/ublas/schur_decomposition.hpp @@ -13,6 +13,8 @@ #include #include #include +#include + namespace boost { namespace numeric { @@ -120,15 +122,62 @@ namespace boost { 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)); + /*value_type h_nn = h(p - size_type(1), p - size_type(1)); //gsl_matrix_get(H, N - 1, N - 1); + value_type h_nm1nm1 = h(q - size_type(1), q - size_type(1)); //gsl_matrix_get(H, N - 2, N - 2); + value_type h_cross = h(p - size_type(1), q - size_type(1)) * h(q - size_type(1), p - size_type(1)); //gsl_matrix_get(H, N - 1, N - 2) * gsl_matrix_get(H, N - 2, N - 1); + + value_type disc = value_type(0.5) * (h_nm1nm1 - h_nn); + disc = disc * disc + h_cross; + if (disc > value_type(0)) + { + value_type ave; + + + disc = sqrt(disc); + ave = value_type(0.5) * (h_nm1nm1 + h_nn); + value_type sign_ave = value_type(1.0); + if (ave < value_type(0)) { + sign_ave = value_type(-1.0); + } + if (((std::abs)(h_nm1nm1) - (std::abs)(h_nn) )> value_type(0)) + { + h_nm1nm1 = h_nm1nm1 * h_nn - h_cross; + h_nn = h_nm1nm1 / (disc * sign_ave + ave); + } + else + { + h_nn = disc * sign_ave + ave; + } + + h_nm1nm1 = h_nn; + h_cross = 0.0; + } + + value_type h_tmp1 = h_nm1nm1 - h(0, 0); + value_type h_tmp2 = h_nn - h(0, 0); + + value_type x = (h_tmp1*h_tmp2 - h_cross) / h( 1, 0) + h(0, 1); + value_type y = h(1, 1) - h(0, 0) - h_tmp1 - h_tmp2; + value_type z = h(2, 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); @@ -158,30 +207,56 @@ namespace boost { vector tx(2); tx(0) = x; tx(1) = y; - vector v; - value_type beta; - householder >(tx, v, beta); - matrix vvt = outer_prod(v, v); - vvt *= beta; - size_type n_vvt = vvt.size1(); - matrix imvvt = identity_matrix(n_vvt) -vvt; + value_type scale = (std::abs)(tx(0)) + (std::abs)(tx(1)) ; + if (scale != value_type(0)) { + tx(0) /= scale; + tx(1) /= scale; + } - matrix t1 = prod(imvvt, project(h, range(q - size_type(1), p), range(p - size_type(3), n))); + 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 = prod(project(h, range(0, p), range(p - size_type(2), p)), imvvt); + 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 = prod(project(qv, range(0,n), range(p - size_type(2),p)), imvvt); - project(qv, range(0, p), range(p - size_type(2), p)).assign(t3); + 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); - if ((std::abs)(h(p - size_type(1), q - size_type(1))) < (EPSA)*((std::abs)(h(p - size_type(1), p - size_type(1))) + (std::abs)(h(q - size_type(1), q - size_type(1))))) { + 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))) < (EPSA)*((std::abs)(h(q - size_type(2), q - size_type(2))) + (std::abs)(h(q - size_type(1), q - 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); From ae0467243e5dff9d4f9a1c8d8019cde86db9c11d Mon Sep 17 00:00:00 2001 From: rajaditya-m Date: Tue, 7 Jul 2015 22:01:53 +0530 Subject: [PATCH 07/20] Complex Eigen Values Complex Eigen Values are supported --- include/boost/numeric/ublas/eigen_solver.hpp | 41 +++++++++++++++++--- 1 file changed, 36 insertions(+), 5 deletions(-) diff --git a/include/boost/numeric/ublas/eigen_solver.hpp b/include/boost/numeric/ublas/eigen_solver.hpp index fb4a80437..939218a70 100644 --- a/include/boost/numeric/ublas/eigen_solver.hpp +++ b/include/boost/numeric/ublas/eigen_solver.hpp @@ -22,17 +22,18 @@ class eigen_solver { M hessenberg_form; M real_schur_form; M transform_accumulations; - bool has_complex; + bool has_complex_part; bool has_eigenvalues; bool has_eigenvectors; - M eigenvalues; + M eigenvalues_real; + M eigenvalues_complex; M eigenvectors; eig_solver_params solver_params; public: //This is the constructor BOOST_UBLAS_INLINE explicit eigen_solver(M &m, eig_solver_params params = EIGVAL) : - matrix(m),has_complex(false),solver_params(params) + matrix(m), has_complex_part(false), solver_params(params) { has_eigenvalues = has_eigenvectors = false; hessenberg_form = M(m); @@ -72,14 +73,44 @@ class eigen_solver { //Right now it doesn't BOOST_UBLAS_INLINE void extract_eigenvalues_from_schur() { - eigenvalues = diagonal_adaptor(real_schur_form); + 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))) { + 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; } //Extraction function BOOST_UBLAS_INLINE M& get_real_eigenvalues() { - return eigenvalues; + return eigenvalues_real; + } + + BOOST_UBLAS_INLINE + M& get_complex_eigenvalues() { + return eigenvalues_complex; } BOOST_UBLAS_INLINE From 3d48bb281b1c4852a6af592a717fd5115019ea60 Mon Sep 17 00:00:00 2001 From: rajaditya-m Date: Wed, 15 Jul 2015 19:11:59 +0530 Subject: [PATCH 08/20] eigenvalue new routine --- include/boost/numeric/ublas/eigen_solver.hpp | 106 +++++++++++++++++-- include/boost/numeric/ublas/hessenberg.hpp | 6 +- 2 files changed, 103 insertions(+), 9 deletions(-) diff --git a/include/boost/numeric/ublas/eigen_solver.hpp b/include/boost/numeric/ublas/eigen_solver.hpp index 939218a70..3d285da91 100644 --- a/include/boost/numeric/ublas/eigen_solver.hpp +++ b/include/boost/numeric/ublas/eigen_solver.hpp @@ -5,6 +5,8 @@ #ifndef _BOOST_UBLAS_EIGENSOLVER_ #define _BOOST_UBLAS_EIGENSOLVER_ +#include +#include #include #include #include @@ -54,23 +56,21 @@ class eigen_solver { 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(); - eigenvectors = M(transform_accumulations); - M temp_q_storage = M(real_schur_form); - block_diag(temp_q_storage, eigenvectors); + extract_eigenvectors(); + has_eigenvalues = true; has_eigenvectors = true; } } - //This is a tricky function because it will do many things... - //Right now it doesn't BOOST_UBLAS_INLINE void extract_eigenvalues_from_schur() { typedef typename M::size_type size_type; @@ -83,7 +83,7 @@ class eigen_solver { eigenvalues_complex = zero_matrix(n, n); while (i < n) { - if ((i == n - size_type(1)) || (real_schur_form(i + 1, i) == value_type(0))) { + 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); } @@ -119,6 +119,95 @@ class eigen_solver { } + 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); + } + + + for (size_type k = n; k--!=size_type(0); ) { + + value_type p = eigenvalues_real(k, k); + value_type q = eigenvalues_complex(k, k); + + 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; + } + } + } + } + + 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; + } + + for (size_type i = 0; i < n; i++) { + vector vi = column(eigenvectors, i); + value_type norm_vi = norm_2(vi); + vector normalized_vi = vi / norm_vi; + column(eigenvectors, i) = normalized_vi; + } + } + //THis is my functions - just kept for debugging will be deleted (or maybe not) BOOST_UBLAS_INLINE M& get_real_schur_form() { @@ -130,6 +219,11 @@ class eigen_solver { return hessenberg_form; } + BOOST_UBLAS_INLINE + M& get_transform_accumulations() { + return transform_accumulations; + } + }; diff --git a/include/boost/numeric/ublas/hessenberg.hpp b/include/boost/numeric/ublas/hessenberg.hpp index d6ff29911..915deb8c7 100644 --- a/include/boost/numeric/ublas/hessenberg.hpp +++ b/include/boost/numeric/ublas/hessenberg.hpp @@ -17,7 +17,7 @@ namespace boost { namespace numeric { namespace ublas { template - typename M::size_type to_hessenberg(M &m) { + void to_hessenberg(M &m) { typedef typename M::size_type size_type; typedef typename M::value_type value_type; @@ -53,7 +53,7 @@ namespace boost { namespace numeric { namespace ublas { } template - typename M::size_type to_hessenberg(M &m, M &u0) { + void to_hessenberg(M &m, M &u0) { typedef typename M::size_type size_type; typedef typename M::value_type value_type; @@ -63,7 +63,7 @@ namespace boost { namespace numeric { namespace ublas { //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 + 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); From 58e0bb89418aeb5eab153a49af634d87d2586be1 Mon Sep 17 00:00:00 2001 From: rajaditya-m Date: Sat, 25 Jul 2015 00:22:33 -0400 Subject: [PATCH 09/20] Complex Eigenvector --- include/boost/numeric/ublas/eigen_solver.hpp | 153 ++++++++++++++++-- include/boost/numeric/ublas/hessenberg.hpp | 2 +- .../numeric/ublas/schur_decomposition.hpp | 94 +++++------ 3 files changed, 183 insertions(+), 66 deletions(-) diff --git a/include/boost/numeric/ublas/eigen_solver.hpp b/include/boost/numeric/ublas/eigen_solver.hpp index 3d285da91..436f3523f 100644 --- a/include/boost/numeric/ublas/eigen_solver.hpp +++ b/include/boost/numeric/ublas/eigen_solver.hpp @@ -5,6 +5,8 @@ #ifndef _BOOST_UBLAS_EIGENSOLVER_ #define _BOOST_UBLAS_EIGENSOLVER_ +#include + #include #include #include @@ -30,6 +32,8 @@ class eigen_solver { M eigenvalues_real; M eigenvalues_complex; M eigenvectors; + M eigenvectors_real; + M eigenvectors_complex; eig_solver_params solver_params; public: //This is the constructor @@ -115,7 +119,12 @@ class eigen_solver { BOOST_UBLAS_INLINE M& get_real_eigenvectors() { - return eigenvectors; + return eigenvectors_real; + } + + BOOST_UBLAS_INLINE + M& get_complex_eigenvectors() { + return eigenvectors_complex; } @@ -137,19 +146,23 @@ class eigen_solver { norm += norm_1(v); } + if (norm == value_type(0)) + return; + + for (size_type k = n; k-- != size_type(0);) { - for (size_type k = n; k--!=size_type(0); ) { - value_type p = eigenvalues_real(k, k); value_type q = eigenvalues_complex(k, k); - value_type lastr, lastw; - lastr = lastw = value_type(0); - size_type l = 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); + T(k, k) = value_type(1); - for (size_type i = k; i--!=size_type(0);) { + 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))); @@ -185,11 +198,94 @@ class eigen_solver { // 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; + 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, n - 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(); + } + + } + + } + + } + k--; + } + } for (size_type j = n; j--!=size_type(0);) @@ -200,11 +296,38 @@ class eigen_solver { column(eigenvectors,j) = v; } - for (size_type i = 0; i < n; i++) { - vector vi = column(eigenvectors, i); - value_type norm_vi = norm_2(vi); - vector normalized_vi = vi / norm_vi; - column(eigenvectors, i) = normalized_vi; + 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++; + } } } diff --git a/include/boost/numeric/ublas/hessenberg.hpp b/include/boost/numeric/ublas/hessenberg.hpp index 915deb8c7..6832dd246 100644 --- a/include/boost/numeric/ublas/hessenberg.hpp +++ b/include/boost/numeric/ublas/hessenberg.hpp @@ -25,7 +25,7 @@ namespace boost { namespace numeric { namespace ublas { //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 + 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) { diff --git a/include/boost/numeric/ublas/schur_decomposition.hpp b/include/boost/numeric/ublas/schur_decomposition.hpp index d24e1ad96..ceebfa5e8 100644 --- a/include/boost/numeric/ublas/schur_decomposition.hpp +++ b/include/boost/numeric/ublas/schur_decomposition.hpp @@ -74,32 +74,51 @@ namespace boost { vector tx(2); tx(0) = x; tx(1) = y; - vector v; - value_type beta; - householder >(tx, v, beta); - matrix vvt = outer_prod(v, v); - vvt *= beta; - size_type n_vvt = vvt.size1(); - matrix imvvt = identity_matrix(n_vvt) -vvt; + 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); - matrix t1 = prod(imvvt, project(h, range(q - size_type(1), p), range(p - size_type(3), n))); + //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 = prod(project(h, range(0, p), range(p - size_type(2), p)), imvvt); + 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))) < (EPSA)*((std::abs)(h(p - size_type(1), p - size_type(1))) + (std::abs)(h(q - size_type(1), q - size_type(1))))) { + 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))) < (EPSA)*((std::abs)(h(q - size_type(2), q - size_type(2))) + (std::abs)(h(q - size_type(1), q - 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); } - std::cout << p << " " << q << "\n"; + } } @@ -122,44 +141,6 @@ namespace boost { 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)); - /*value_type h_nn = h(p - size_type(1), p - size_type(1)); //gsl_matrix_get(H, N - 1, N - 1); - value_type h_nm1nm1 = h(q - size_type(1), q - size_type(1)); //gsl_matrix_get(H, N - 2, N - 2); - value_type h_cross = h(p - size_type(1), q - size_type(1)) * h(q - size_type(1), p - size_type(1)); //gsl_matrix_get(H, N - 1, N - 2) * gsl_matrix_get(H, N - 2, N - 1); - - value_type disc = value_type(0.5) * (h_nm1nm1 - h_nn); - disc = disc * disc + h_cross; - if (disc > value_type(0)) - { - value_type ave; - - - disc = sqrt(disc); - ave = value_type(0.5) * (h_nm1nm1 + h_nn); - value_type sign_ave = value_type(1.0); - if (ave < value_type(0)) { - sign_ave = value_type(-1.0); - } - if (((std::abs)(h_nm1nm1) - (std::abs)(h_nn) )> value_type(0)) - { - h_nm1nm1 = h_nm1nm1 * h_nn - h_cross; - h_nn = h_nm1nm1 / (disc * sign_ave + ave); - } - else - { - h_nn = disc * sign_ave + ave; - } - - h_nm1nm1 = h_nn; - h_cross = 0.0; - } - - value_type h_tmp1 = h_nm1nm1 - h(0, 0); - value_type h_tmp2 = h_nn - h(0, 0); - - value_type x = (h_tmp1*h_tmp2 - h_cross) / h( 1, 0) + h(0, 1); - value_type y = h(1, 1) - h(0, 0) - h_tmp1 - h_tmp2; - value_type z = h(2, 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++) { @@ -251,6 +232,19 @@ namespace boost { } 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); From 2e733521190ff24875c63b833be45f4e7c7f22ee Mon Sep 17 00:00:00 2001 From: rajaditya-m Date: Mon, 27 Jul 2015 12:24:59 -0400 Subject: [PATCH 10/20] Small Numerical Fix --- include/boost/numeric/ublas/eigen_solver.hpp | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/include/boost/numeric/ublas/eigen_solver.hpp b/include/boost/numeric/ublas/eigen_solver.hpp index 436f3523f..2e09064bd 100644 --- a/include/boost/numeric/ublas/eigen_solver.hpp +++ b/include/boost/numeric/ublas/eigen_solver.hpp @@ -242,7 +242,6 @@ class eigen_solver { } else { l = i; - if (eigenvalues_complex(i, i) == value_type(0)) { std::complex x(-ra, -sa); std::complex y(w, q); @@ -268,7 +267,7 @@ class eigen_solver { 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, n - size_type(1))) / 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)); @@ -282,12 +281,21 @@ class eigen_solver { } + 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--; } - } + std::cout << eigenvectors << "\n"; + for (size_type j = n; j--!=size_type(0);) { M matrix_left = project(eigenvectors, range(0, n),range(0, j + size_type(1))); @@ -296,6 +304,8 @@ class eigen_solver { column(eigenvectors,j) = v; } + //std::cout << eigenvectors << "\n"; + eigenvectors_real = zero_matrix(n, n); eigenvectors_complex = zero_matrix(n, n); From 803d24d38b33219a0c03a374a55a5be0eccc01d3 Mon Sep 17 00:00:00 2001 From: rajaditya-m Date: Tue, 28 Jul 2015 16:30:11 -0400 Subject: [PATCH 11/20] Matrix Balancing --- include/boost/numeric/ublas/eigen_solver.hpp | 1 + .../boost/numeric/ublas/matrix_balancing.hpp | 103 ++++++++++++++++++ 2 files changed, 104 insertions(+) create mode 100644 include/boost/numeric/ublas/matrix_balancing.hpp diff --git a/include/boost/numeric/ublas/eigen_solver.hpp b/include/boost/numeric/ublas/eigen_solver.hpp index 2e09064bd..7f2f7548c 100644 --- a/include/boost/numeric/ublas/eigen_solver.hpp +++ b/include/boost/numeric/ublas/eigen_solver.hpp @@ -12,6 +12,7 @@ #include #include #include +#include namespace boost { namespace numeric { namespace ublas { diff --git a/include/boost/numeric/ublas/matrix_balancing.hpp b/include/boost/numeric/ublas/matrix_balancing.hpp new file mode 100644 index 000000000..3596b04b6 --- /dev/null +++ b/include/boost/numeric/ublas/matrix_balancing.hpp @@ -0,0 +1,103 @@ +// 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) + + +#ifndef _BOOST_UBLAS_MATRIXBALANCING_ +#define _BOOST_UBLAS_MATRIXBALANCING_ + +#include +#include +#include +#include +#include + +namespace boost { + namespace numeric { + namespace ublas { + + 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; + + } + } + } + } + + 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 From 147766ed708f633a86c02e25b50d7cd72473fb85 Mon Sep 17 00:00:00 2001 From: rajaditya-m Date: Sat, 1 Aug 2015 01:06:17 -0400 Subject: [PATCH 12/20] Doxygen Comments Documentations --- include/boost/numeric/ublas/eigen_solver.hpp | 37 ++++++- include/boost/numeric/ublas/hessenberg.hpp | 35 +++--- include/boost/numeric/ublas/householder.hpp | 37 +++++-- .../boost/numeric/ublas/matrix_balancing.hpp | 12 +++ .../numeric/ublas/schur_decomposition.hpp | 100 +++--------------- 5 files changed, 109 insertions(+), 112 deletions(-) diff --git a/include/boost/numeric/ublas/eigen_solver.hpp b/include/boost/numeric/ublas/eigen_solver.hpp index 7f2f7548c..e6b3125c5 100644 --- a/include/boost/numeric/ublas/eigen_solver.hpp +++ b/include/boost/numeric/ublas/eigen_solver.hpp @@ -1,7 +1,11 @@ +// 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_ @@ -20,6 +24,12 @@ namespace boost { namespace numeric { namespace ublas { //Ask mentor about how to "formally" add them to ublas typedef 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: @@ -37,7 +47,12 @@ class eigen_solver { M eigenvectors_complex; eig_solver_params solver_params; public: - //This is the constructor + /// \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) @@ -47,7 +62,11 @@ class eigen_solver { compute(solver_params); } - //The compute function + /// \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) { @@ -107,22 +126,25 @@ class eigen_solver { has_eigenvalues = true; } - //Extraction function + /// \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; @@ -342,22 +364,29 @@ class eigen_solver { } } - //THis is my functions - just kept for debugging will be deleted (or maybe not) + /// \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; + } }; diff --git a/include/boost/numeric/ublas/hessenberg.hpp b/include/boost/numeric/ublas/hessenberg.hpp index 6832dd246..abc171a3a 100644 --- a/include/boost/numeric/ublas/hessenberg.hpp +++ b/include/boost/numeric/ublas/hessenberg.hpp @@ -1,6 +1,10 @@ -// 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) +// 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_ @@ -11,11 +15,18 @@ #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) { @@ -43,15 +54,17 @@ namespace boost { namespace numeric { namespace ublas { 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); - - //@TODO: Ask mentor why uncommented lines are more accurate than commented lines... - //project(m, range(i + 1, n), range(i, n)).assign(prod(imvvt, project(m, range(i + 1, n), range(i, n)))); - //project(m, range(0, n), range(i+1, n)).assign(prod(project(m, range(0, n), range(i+1, n)),imvvt)); - } } + /// \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) { @@ -73,8 +86,6 @@ namespace boost { namespace numeric { namespace ublas { value_type beta; householder(x, v, beta); - //std::cout << i << " " << x << " " << v << "\n"; - matrix vvt = outer_prod(v, v); vvt *= beta; size_type n_vvt = vvt.size1(); @@ -91,10 +102,6 @@ namespace boost { namespace numeric { namespace ublas { matrix t2 = prod(project(m, range(0, n), range(i + 1, n)), imvvt); project(m, range(0, n), range(i + 1, n)).assign(t2); - //@TODO: Ask mentor why uncommented lines are more accurate than commented lines... - //project(m, range(i + 1, n), range(i, n)).assign(prod(imvvt, project(m, range(i + 1, n), range(i, n)))); - //project(m, range(0, n), range(i+1, n)).assign(prod(project(m, range(0, n), range(i+1, n)),imvvt)); - } } diff --git a/include/boost/numeric/ublas/householder.hpp b/include/boost/numeric/ublas/householder.hpp index c6e9a64ad..b14668b51 100644 --- a/include/boost/numeric/ublas/householder.hpp +++ b/include/boost/numeric/ublas/householder.hpp @@ -1,6 +1,10 @@ -// 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) +// 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_ @@ -15,8 +19,16 @@ 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 - typename M::size_type householder(M &x, M &v, typename M::value_type &beta) { + 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; @@ -26,7 +38,7 @@ namespace boost { size_type n = x.size(); value_type sigma = inner_prod(project(x, range(1, n)), project(x, range(1, n))); - v = M(x); //Ask Mentor how to evoke copy constructor + v = M(x); v(size_type(0)) = value_type(1); @@ -38,7 +50,7 @@ namespace boost { v(size_type(0)) = x(size_type(0)) - mu; } else { - v(size_type(0)) = (-sigma) / (x(size_type(0)) + mu); // Ask mentor is -sigma is valid + 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)); @@ -49,15 +61,20 @@ namespace boost { beta = value_type(0); } else { - beta = value_type(-2); // Ask mentor if this is Valid + beta = value_type(-2); } } - - - return size_type(0); } + /// \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)) { diff --git a/include/boost/numeric/ublas/matrix_balancing.hpp b/include/boost/numeric/ublas/matrix_balancing.hpp index 3596b04b6..2bfd19126 100644 --- a/include/boost/numeric/ublas/matrix_balancing.hpp +++ b/include/boost/numeric/ublas/matrix_balancing.hpp @@ -1,7 +1,10 @@ +// 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_ @@ -16,6 +19,11 @@ 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; @@ -80,6 +88,10 @@ namespace boost { } } + /// \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; diff --git a/include/boost/numeric/ublas/schur_decomposition.hpp b/include/boost/numeric/ublas/schur_decomposition.hpp index ceebfa5e8..7a3900623 100644 --- a/include/boost/numeric/ublas/schur_decomposition.hpp +++ b/include/boost/numeric/ublas/schur_decomposition.hpp @@ -1,7 +1,12 @@ +// 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_ @@ -20,6 +25,11 @@ namespace boost { namespace numeric { namespace ublas { + /// \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). + /// \param m matrix type (like matrix) - input hessenberg form output - schur form template void schur_decomposition(M &h) { @@ -122,6 +132,12 @@ namespace boost { } } + /// \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) { @@ -255,93 +271,9 @@ namespace boost { p = p - size_type(2); q = p - size_type(1); } - //std::cout << p << " " << q << "\n"; - } - } - - - template - void inverse_iterations(M &h, M &u0, M &v) { - typedef typename M::size_type size_type; - typedef typename M::value_type value_type; - - size_type n = h.size1(); - - v = identity_matrix (n); - matrix in = identity_matrix(n); - //This is the main loop - for (size_type i = size_type(0); i < n; i++) { - value_type mu = h(i, i); - matrix mu_i = mu * in; - matrix mat_lhs = h - mu_i; - size_type max_iters = size_type(1000); //This should be adjustable - vector q = column(v, i); - vector z; - int counter = 0; - //This is the iteration loop - bool convergence_reached = false; - while (!convergence_reached) { - //Start solve of mat_lhs * z = q via column oriented back-substitution - //vector q_copy = q; - for (size_type j = n; j >= size_type(2); j--) { - q(j - size_type(1)) = q(j - size_type(1)) / mat_lhs(j - size_type(1), j - size_type(1)); - vector u_ranged = project(column(mat_lhs, j - size_type(1)), range(size_type(0), j - size_type(1))); - u_ranged *= q(j - size_type(1)); - vector b_ranged = project(q, range(size_type(0), j - size_type(1))); - b_ranged -= u_ranged; - project(q, range(size_type(0), j - size_type(1))).assign(b_ranged); - } - q(size_type(0)) = q(size_type(0)) / mat_lhs(size_type(0), size_type(0)); - - //Normalize the q - value_type q_norm = norm_2(q); - //q = q / q_norm; - - counter++; - std::cout << counter << " " << q << "\n"; - if (counter==10) - convergence_reached = true; - } - - column(v, i) = q; - - } - - } - template - void block_diag(M &t, M &q) { - typedef typename M::size_type size_type; - typedef typename M::value_type value_type; - - size_type n = t.size1(); - - for (size_type j = size_type(2); j <= n; j++) { - for (size_type i = size_type(1); i < j; i++) { - value_type tii = t(i - size_type(1), i - size_type(1)); - value_type tjj = t(j - size_type(1), j - size_type(1)); - value_type tij = t(i - size_type(1), j - size_type(1)); - value_type z = -tij / (tii - tjj); - for (size_type k = j + size_type(1); k <= n; k++) { - t(i - size_type(1), k - size_type(1)) -= (z*t(j - size_type(1), k - size_type(1))); - } - for (size_type k = size_type(1); k <= n; k++) { - q(k - size_type(1), j - size_type(1)) += (z*q(k - size_type(1), i - size_type(1))); - } - } - } - - //We need to normalize the columns also - for (size_type i = 0; i < n; i++) { - vector vi = column(q, i); - value_type norm_vi = norm_2(vi); - vector normalized_vi = vi / norm_vi; - column(q, i) = normalized_vi; - } - - } }}} From 24531ff62f89372e280b514c1b55ddf9bcf132af Mon Sep 17 00:00:00 2001 From: rajaditya-m Date: Wed, 5 Aug 2015 19:28:57 -0400 Subject: [PATCH 13/20] SingleTest and sample --- doc/samples/eigen_solver.cpp | 54 +++++++++++ .../numeric/ublas/schur_decomposition.hpp | 12 ++- test/common/testhelper.hpp | 25 +++++ test/test_eigensolver.cpp | 92 +++++++++++++++++++ 4 files changed, 178 insertions(+), 5 deletions(-) create mode 100644 doc/samples/eigen_solver.cpp create mode 100644 test/test_eigensolver.cpp 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/schur_decomposition.hpp b/include/boost/numeric/ublas/schur_decomposition.hpp index 7a3900623..cc8dd4760 100644 --- a/include/boost/numeric/ublas/schur_decomposition.hpp +++ b/include/boost/numeric/ublas/schur_decomposition.hpp @@ -25,11 +25,13 @@ namespace boost { namespace numeric { namespace ublas { - /// \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). - /// \param m matrix type (like matrix) - input hessenberg form output - schur form + /*! @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) { diff --git a/test/common/testhelper.hpp b/test/common/testhelper.hpp index b554511fa..36182ad14 100644 --- a/test/common/testhelper.hpp +++ b/test/common/testhelper.hpp @@ -69,6 +69,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..e90e85b9e --- /dev/null +++ b/test/test_eigensolver.cpp @@ -0,0 +1,92 @@ +// 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"; +static const string matrix_EVALR = "[3,3]((1,0.0,0.0),(0.0,2,0.0),(0.0,0.0,3))\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"; + +int main() { + + typedef double TYPE; + + 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); + + + { + std::istringstream is(matrix_IN); + is >> A; + } + + { + std::istringstream is(matrix_EVALR); + is >> GT_EVALR; + } + + { + std::istringstream is(matrix_EVALC); + 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)); + + return (getResults().second > 0) ? boost::exit_failure : boost::exit_success; +} From 929ec7062003a6bb2c614a7cdb87a4ce72fa5d19 Mon Sep 17 00:00:00 2001 From: rajaditya-m Date: Wed, 12 Aug 2015 16:06:04 -0400 Subject: [PATCH 14/20] Test Cases --- include/boost/numeric/ublas/eigen_solver.hpp | 2 - .../numeric/ublas/schur_decomposition.hpp | 4 +- test/test_eigensolver.cpp | 113 +++++++++++------- 3 files changed, 71 insertions(+), 48 deletions(-) diff --git a/include/boost/numeric/ublas/eigen_solver.hpp b/include/boost/numeric/ublas/eigen_solver.hpp index e6b3125c5..1890329e8 100644 --- a/include/boost/numeric/ublas/eigen_solver.hpp +++ b/include/boost/numeric/ublas/eigen_solver.hpp @@ -317,8 +317,6 @@ class eigen_solver { } } - std::cout << eigenvectors << "\n"; - for (size_type j = n; j--!=size_type(0);) { M matrix_left = project(eigenvectors, range(0, n),range(0, j + size_type(1))); diff --git a/include/boost/numeric/ublas/schur_decomposition.hpp b/include/boost/numeric/ublas/schur_decomposition.hpp index cc8dd4760..3c2caa478 100644 --- a/include/boost/numeric/ublas/schur_decomposition.hpp +++ b/include/boost/numeric/ublas/schur_decomposition.hpp @@ -251,7 +251,7 @@ namespace boost { 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 ci = 0; ci < p; ++ci) { for (size_type ri = p; ri < n; ri++) { h(ri, ci) = value_type(0); @@ -261,7 +261,7 @@ namespace boost { 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); diff --git a/test/test_eigensolver.cpp b/test/test_eigensolver.cpp index e90e85b9e..af706bb51 100644 --- a/test/test_eigensolver.cpp +++ b/test/test_eigensolver.cpp @@ -18,9 +18,19 @@ 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"; -static const string matrix_EVALR = "[3,3]((1,0.0,0.0),(0.0,2,0.0),(0.0,0.0,3))\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"; +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"}; +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" }; +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" }; + int main() { @@ -37,56 +47,71 @@ int main() { MATRIX GT_EVALC; MATRIX Zero_Matrix = zero_matrix(3,3); + int numTestCases = 4; - { - std::istringstream is(matrix_IN); - is >> A; - } - { - std::istringstream is(matrix_EVALR); - is >> GT_EVALR; - } + for (int i = 0; i < numTestCases; i++){ - { - std::istringstream is(matrix_EVALC); - is >> GT_EVALC; - } + std::cout << "Running test case# " << (i + 1) << "\n"; - 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)); + { + std::istringstream is(matrix_IN[i]); + is >> A; + } - matrix > V(3, 3); - matrix > D(3, 3); - matrix > M(3, 3); - matrix > Lambda; + { + std::istringstream is(matrix_EVALR[i]); + is >> GT_EVALR; + } - MATRIX Lambda_Real(3, 3); - MATRIX Lambda_Complex(3, 3); + { + std::istringstream is(matrix_EVALC[i]); + is >> GT_EVALC; + } - 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); + 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(); + 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)); + } - 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)); - - return (getResults().second > 0) ? boost::exit_failure : boost::exit_success; + if (!_fail_counter) { + std::cout << "\nEigen Solver regression suite passed.\n"; + return boost::exit_failure; + } + else { + std::cout << "\nEigen Solver regression suite failed.\n"; + return boost::exit_success; + } } From ac89993eea0079a3470ffeb4221abd2bf9e43bda Mon Sep 17 00:00:00 2001 From: rajaditya-m Date: Wed, 12 Aug 2015 22:24:51 -0400 Subject: [PATCH 15/20] HTML Docs --- doc/eigensolver.html | 146 +++++++++++++++++++++++++++++++++++++++++++ doc/index.html | 6 ++ 2 files changed, 152 insertions(+) create mode 100644 doc/eigensolver.html 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 + + +

Boost uBLAS EigenSolver

+
+

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

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
MemberDescription
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
+ Use, modification and distribution are subject to 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 + ). +

+ + + diff --git a/doc/index.html b/doc/index.html index 6a3a15dfd..249a60c07 100644 --- a/doc/index.html +++ b/doc/index.html @@ -297,6 +297,12 @@

Documentation

+ +
  • + Numerical Algebra + +
    • + From 1cbafefb4bccadc970655ca982a91747e663a4db Mon Sep 17 00:00:00 2001 From: rajaditya-m Date: Thu, 15 Oct 2015 17:52:47 -0400 Subject: [PATCH 16/20] Matlab's Magic Shift Function --- .../numeric/ublas/schur_decomposition.hpp | 251 ++++++++++++++++++ 1 file changed, 251 insertions(+) diff --git a/include/boost/numeric/ublas/schur_decomposition.hpp b/include/boost/numeric/ublas/schur_decomposition.hpp index 3c2caa478..6ff7ad4a6 100644 --- a/include/boost/numeric/ublas/schur_decomposition.hpp +++ b/include/boost/numeric/ublas/schur_decomposition.hpp @@ -276,6 +276,257 @@ namespace boost { } } + 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); + } + + } + }}} From aaee56c87ef4ed29f02f17696e8304cfec9288cb Mon Sep 17 00:00:00 2001 From: Thomas Yang Date: Sat, 27 Apr 2019 15:27:08 -0500 Subject: [PATCH 17/20] Added test cases for eigensolver; run test for different data type --- test/test_eigensolver.cpp | 94 +++++++++++++++++++++++++++++++++++---- 1 file changed, 86 insertions(+), 8 deletions(-) diff --git a/test/test_eigensolver.cpp b/test/test_eigensolver.cpp index af706bb51..c0c22e179 100644 --- a/test/test_eigensolver.cpp +++ b/test/test_eigensolver.cpp @@ -20,21 +20,43 @@ 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"}; + "[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" }; + "[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.000,00),(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" }; + "[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; @@ -47,7 +69,7 @@ int main() { MATRIX GT_EVALC; MATRIX Zero_Matrix = zero_matrix(3,3); - int numTestCases = 4; + int numTestCases = 9; for (int i = 0; i < numTestCases; i++){ @@ -59,7 +81,7 @@ int main() { is >> A; } - { + { std::istringstream is(matrix_EVALR[i]); is >> GT_EVALR; } @@ -91,8 +113,64 @@ int main() { 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)); + + } + + for (int i = 0; i < numTestCases; i++){ + + std::cout << "Running test case# " << (i + 1) << "\n"; + + { + std::istringstream is(matrix_IN[i]); + is >> A; + } + + { + std::istringstream is(matrix_EVALR[i]); + is >> GT_EVALR; + } + + { + std::istringstream is(matrix_EVALC[i]); + 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++){ From 0b808761ad24037ec849fb83a98f28b9801b9462 Mon Sep 17 00:00:00 2001 From: Thomas Yang Date: Sat, 4 May 2019 16:34:14 -0500 Subject: [PATCH 18/20] Modified test Jamfile --- test/Jamfile | 5 +++++ 1 file changed, 5 insertions(+) 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 : : From 85c3867b436378899a9e21ab90678f01718c21e8 Mon Sep 17 00:00:00 2001 From: Thomas Yang Date: Sun, 5 May 2019 01:35:55 -0500 Subject: [PATCH 19/20] Fixed backslash error in eigen_solver.hpp and matrix_balancing.hpp --- include/boost/numeric/ublas/eigen_solver.hpp | 12 +++---- .../boost/numeric/ublas/matrix_balancing.hpp | 8 ++--- test/test_eigensolver.cpp | 31 ++++++++++--------- 3 files changed, 26 insertions(+), 25 deletions(-) diff --git a/include/boost/numeric/ublas/eigen_solver.hpp b/include/boost/numeric/ublas/eigen_solver.hpp index 1890329e8..cc5c2688a 100644 --- a/include/boost/numeric/ublas/eigen_solver.hpp +++ b/include/boost/numeric/ublas/eigen_solver.hpp @@ -11,12 +11,12 @@ #include -#include -#include -#include -#include -#include -#include +#include +#include +#include +#include +#include +#include namespace boost { namespace numeric { namespace ublas { diff --git a/include/boost/numeric/ublas/matrix_balancing.hpp b/include/boost/numeric/ublas/matrix_balancing.hpp index 2bfd19126..f755ffc8c 100644 --- a/include/boost/numeric/ublas/matrix_balancing.hpp +++ b/include/boost/numeric/ublas/matrix_balancing.hpp @@ -10,10 +10,10 @@ #define _BOOST_UBLAS_MATRIXBALANCING_ #include -#include -#include -#include -#include +#include +#include +#include +#include namespace boost { namespace numeric { diff --git a/test/test_eigensolver.cpp b/test/test_eigensolver.cpp index c0c22e179..e000d12a8 100644 --- a/test/test_eigensolver.cpp +++ b/test/test_eigensolver.cpp @@ -27,18 +27,18 @@ static const string matrix_IN[] = { "[3,3]((-149.0,-50.0,-154.0),(537.0,180.0,54 "[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" + "[10,10]((0, 20000, 30000, 40000, 30000, 60000, 70000, 80000, 90000, 100000),(10000, 20000, 30000, 30000, 50000, 0, 70000, 80000, 90000, 100000),(30000, 20000, 0, 42000, 50000, 61000, 0, 80 000, 90000, 0),(10000, 20000, 30000, 40000, 50000, 60000, 0, 80000, 90000, 10000),(10000, 20000, 30000, 45000, 50000, 60000, 70000, 80000, 90000, 0),(50000, 60000, 0, 40000, 50000, 60000, 70000, 80000, 90000, 10000),(0, 20000, 30000, 0, 50000, 60000, 70000, 10000, 90000, 40000),(10000, 20000, 0, 40000, 50000, 60000, 70000, 80000, 93000, 0),(0, 20000, 30000, 40000, 0, 60000, 70000, 80000, 90000, 1000 00),(0, 20000, 30000, 40000, 50000, 60000, 70000, 80000, 90000, 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.000,00),(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", + "[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" + "[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]((445477, 0, 0, 0, 0, 0, 0, 0, 0, 0),(0, -16305.1, 0, 0, 0, 0, 0, 0, 0, 0),(0, 0, -16305.1, 0, 0, 0, 0, 0, 0, 0),(0, 0, 0, -39098, 0, 0, 0, 0, 0, 0),(0, 0, 0, 0, -15122.7, 0, 0 , 0, 0, 0),(0, 0, 0, 0, 0, -15122.7, 0, 0, 0, 0),(0, 0, 0, 0, 0, 0, 25925.2, 0, 0, 0),(0, 0, 0, 0, 0, 0, 0, 25925.2, 0, 0),(0, 0, 0, 0, 0, 0, 0, 0, 7313.04, 0),(0, 0, 0, 0, 0, 0, 0, 0, 0, 7313.04))\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", @@ -49,7 +49,7 @@ static const string matrix_EVALC[] = { "[3,3]((0.0,0.0,0.0),(0.0,0.0,0.0),(0.0,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" + "[10,10]((0, 0, 0, 0, 0, 0, 0, 0, 0, 0),(0, 53455.2, 0, 0, 0, 0, 0, 0, 0, 0),(0, 0, -53455.2, 0, 0, 0, 0, 0, 0, 0),(0, 0, 0, 0, 0, 0, 0, 0, 0, 0),(0, 0, 0, 0, 25827.8, 0, 0, 0, 0, 0),( 0, 0, 0, 0, 0, -25827.8, 0, 0, 0, 0),(0, 0, 0, 0, 0, 0, 15857.9, 0, 0, 0),(0, 0, 0, 0, 0, 0, 0, -15857.9, 0, 0),(0, 0, 0, 0, 0, 0, 0, 0, 5775.81, 0),(0, 0, 0, 0, 0, 0, 0, 0, 0, -5775.81))\0" }; @@ -72,22 +72,22 @@ int main() { int numTestCases = 9; - for (int i = 0; i < numTestCases; i++){ + for (int k = 0; k < numTestCases; k++){ - std::cout << "Running test case# " << (i + 1) << "\n"; + std::cout << "Running test case# " << (k + 1) << "\n"; { - std::istringstream is(matrix_IN[i]); + std::istringstream is(matrix_IN[k]); is >> A; } { - std::istringstream is(matrix_EVALR[i]); + std::istringstream is(matrix_EVALR[k]); is >> GT_EVALR; } { - std::istringstream is(matrix_EVALC[i]); + std::istringstream is(matrix_EVALC[k]); is >> GT_EVALC; } @@ -97,6 +97,7 @@ int main() { 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)); @@ -128,22 +129,22 @@ int main() { } - for (int i = 0; i < numTestCases; i++){ + for (int k = 0; k < numTestCases; k++){ - std::cout << "Running test case# " << (i + 1) << "\n"; + std::cout << "Running test case# " << (k + 1) << "\n"; { - std::istringstream is(matrix_IN[i]); + std::istringstream is(matrix_IN[k]); is >> A; } { - std::istringstream is(matrix_EVALR[i]); + std::istringstream is(matrix_EVALR[k]); is >> GT_EVALR; } { - std::istringstream is(matrix_EVALC[i]); + std::istringstream is(matrix_EVALC[k]); is >> GT_EVALC; } From 6b89a4a4453c160a42b431fa1d3da538ec2f884b Mon Sep 17 00:00:00 2001 From: Thomas Yang Date: Sun, 5 May 2019 13:46:40 -0500 Subject: [PATCH 20/20] Fixed test cases; removed "typedef" in Eigensolver file --- include/boost/numeric/ublas/eigen_solver.hpp | 2 +- test/test_eigensolver.cpp | 69 ++------------------ 2 files changed, 8 insertions(+), 63 deletions(-) diff --git a/include/boost/numeric/ublas/eigen_solver.hpp b/include/boost/numeric/ublas/eigen_solver.hpp index cc5c2688a..fd50c19d0 100644 --- a/include/boost/numeric/ublas/eigen_solver.hpp +++ b/include/boost/numeric/ublas/eigen_solver.hpp @@ -22,7 +22,7 @@ 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 -typedef enum eig_solver_params {EIGVAL,EIGVEC}; +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. diff --git a/test/test_eigensolver.cpp b/test/test_eigensolver.cpp index e000d12a8..a4a33ce4a 100644 --- a/test/test_eigensolver.cpp +++ b/test/test_eigensolver.cpp @@ -27,7 +27,7 @@ static const string matrix_IN[] = { "[3,3]((-149.0,-50.0,-154.0),(537.0,180.0,54 "[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, 20000, 30000, 40000, 30000, 60000, 70000, 80000, 90000, 100000),(10000, 20000, 30000, 30000, 50000, 0, 70000, 80000, 90000, 100000),(30000, 20000, 0, 42000, 50000, 61000, 0, 80 000, 90000, 0),(10000, 20000, 30000, 40000, 50000, 60000, 0, 80000, 90000, 10000),(10000, 20000, 30000, 45000, 50000, 60000, 70000, 80000, 90000, 0),(50000, 60000, 0, 40000, 50000, 60000, 70000, 80000, 90000, 10000),(0, 20000, 30000, 0, 50000, 60000, 70000, 10000, 90000, 40000),(10000, 20000, 0, 40000, 50000, 60000, 70000, 80000, 93000, 0),(0, 20000, 30000, 40000, 0, 60000, 70000, 80000, 90000, 1000 00),(0, 20000, 30000, 40000, 50000, 60000, 70000, 80000, 90000, 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", @@ -38,7 +38,7 @@ static const string matrix_EVALR[] = { "[3,3]((1,0.0,0.0),(0.0,2,0.0),(0.0,0.0,3 "[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]((445477, 0, 0, 0, 0, 0, 0, 0, 0, 0),(0, -16305.1, 0, 0, 0, 0, 0, 0, 0, 0),(0, 0, -16305.1, 0, 0, 0, 0, 0, 0, 0),(0, 0, 0, -39098, 0, 0, 0, 0, 0, 0),(0, 0, 0, 0, -15122.7, 0, 0 , 0, 0, 0),(0, 0, 0, 0, 0, -15122.7, 0, 0, 0, 0),(0, 0, 0, 0, 0, 0, 25925.2, 0, 0, 0),(0, 0, 0, 0, 0, 0, 0, 25925.2, 0, 0),(0, 0, 0, 0, 0, 0, 0, 0, 7313.04, 0),(0, 0, 0, 0, 0, 0, 0, 0, 0, 7313.04))\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", @@ -49,7 +49,7 @@ static const string matrix_EVALC[] = { "[3,3]((0.0,0.0,0.0),(0.0,0.0,0.0),(0.0,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, 53455.2, 0, 0, 0, 0, 0, 0, 0, 0),(0, 0, -53455.2, 0, 0, 0, 0, 0, 0, 0),(0, 0, 0, 0, 0, 0, 0, 0, 0, 0),(0, 0, 0, 0, 25827.8, 0, 0, 0, 0, 0),( 0, 0, 0, 0, 0, -25827.8, 0, 0, 0, 0),(0, 0, 0, 0, 0, 0, 15857.9, 0, 0, 0),(0, 0, 0, 0, 0, 0, 0, -15857.9, 0, 0),(0, 0, 0, 0, 0, 0, 0, 0, 5775.81, 0),(0, 0, 0, 0, 0, 0, 0, 0, 0, -5775.81))\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" }; @@ -69,7 +69,7 @@ int main() { MATRIX GT_EVALC; MATRIX Zero_Matrix = zero_matrix(3,3); - int numTestCases = 9; + int numTestCases = 8; for (int k = 0; k < numTestCases; k++){ @@ -91,6 +91,7 @@ int main() { is >> GT_EVALC; } + eigen_solver es(A, EIGVEC); matrix evals_r = es.get_real_eigenvalues(); matrix evals_c = es.get_complex_eigenvalues(); @@ -129,68 +130,12 @@ int main() { } - 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_failure; + return boost::exit_success; } else { std::cout << "\nEigen Solver regression suite failed.\n"; - return boost::exit_success; + return boost::exit_failure; } }