diff options
author | Chen-Pang He <jdh8@ms63.hinet.net> | 2013-07-15 00:43:14 +0800 |
---|---|---|
committer | Chen-Pang He <jdh8@ms63.hinet.net> | 2013-07-15 00:43:14 +0800 |
commit | 9be658f7015161989b7ecccd70fd050ce563cad9 (patch) | |
tree | 070c28b551973bffe3452d40d8a6b727559dfc71 /unsupported/test/matrix_functions.h | |
parent | b8f0364a1c56784fa666c783e21c4bd1b218b9b0 (diff) |
generateTestMatrix can use processTriangularMatrix
Diffstat (limited to 'unsupported/test/matrix_functions.h')
-rw-r--r-- | unsupported/test/matrix_functions.h | 41 |
1 files changed, 31 insertions, 10 deletions
diff --git a/unsupported/test/matrix_functions.h b/unsupported/test/matrix_functions.h index 5817caef6..295da16b6 100644 --- a/unsupported/test/matrix_functions.h +++ b/unsupported/test/matrix_functions.h @@ -10,27 +10,48 @@ #include "main.h" #include <unsupported/Eigen/MatrixFunctions> +// For complex matrices, any matrix is fine. +template<typename MatrixType, int IsComplex = NumTraits<typename internal::traits<MatrixType>::Scalar>::IsComplex> +struct processTriangularMatrix +{ + static void run(MatrixType&, MatrixType&, const MatrixType&) + { } +}; + +// For real matrices, make sure none of the eigenvalues are negative. +template<typename MatrixType> +struct processTriangularMatrix<MatrixType,0> +{ + static void run(MatrixType& m, MatrixType& T, const MatrixType& U) + { + typedef typename MatrixType::Index Index; + const Index size = m.cols(); + + for (Index i=0; i < size; ++i) { + if (i == size - 1 || T.coeff(i+1,i) == 0) + T.coeffRef(i,i) = std::abs(T.coeff(i,i)); + else + ++i; + } + m = U * T * U.transpose(); + } +}; + template <typename MatrixType, int IsComplex = NumTraits<typename internal::traits<MatrixType>::Scalar>::IsComplex> struct generateTestMatrix; -// for real matrices, make sure none of the eigenvalues are negative template <typename MatrixType> struct generateTestMatrix<MatrixType,0> { static void run(MatrixType& result, typename MatrixType::Index size) { - MatrixType mat = MatrixType::Random(size, size); - EigenSolver<MatrixType> es(mat); - typename EigenSolver<MatrixType>::EigenvalueType eivals = es.eigenvalues(); - for (typename MatrixType::Index i = 0; i < size; ++i) { - if (eivals(i).imag() == 0 && eivals(i).real() < 0) - eivals(i) = -eivals(i); - } - result = (es.eigenvectors() * eivals.asDiagonal() * es.eigenvectors().inverse()).real(); + result = MatrixType::Random(size, size); + RealSchur<MatrixType> schur(result); + MatrixType T = schur.matrixT(); + processTriangularMatrix<MatrixType>::run(result, T, schur.matrixU()); } }; -// for complex matrices, any matrix is fine template <typename MatrixType> struct generateTestMatrix<MatrixType,1> { |