diff options
author | Gael Guennebaud <g.gael@free.fr> | 2016-06-09 17:11:03 +0200 |
---|---|---|
committer | Gael Guennebaud <g.gael@free.fr> | 2016-06-09 17:11:03 +0200 |
commit | c1f9ca925405c8fad126f327b4bdca7f983fc96e (patch) | |
tree | 06f6cb2be66d89cc553426d335a8c3ee6b0c4689 | |
parent | 15890c304edbccedc8a989468ed3fc475f428059 (diff) |
Update RealQZ to reduce 2x2 diagonal block of T corresponding to non reduced diagonal block of S to positive diagonal form.
This step involve a real 2x2 SVD problem. The respective routine is thus in src/misc/ to be shared by both EVD and AVD modules.
-rw-r--r-- | Eigen/Eigenvalues | 1 | ||||
-rw-r--r-- | Eigen/SVD | 1 | ||||
-rw-r--r-- | Eigen/src/Eigenvalues/RealQZ.h | 32 | ||||
-rw-r--r-- | Eigen/src/SVD/JacobiSVD.h | 32 | ||||
-rw-r--r-- | Eigen/src/misc/RealSvd2x2.h | 54 |
5 files changed, 87 insertions, 33 deletions
diff --git a/Eigen/Eigenvalues b/Eigen/Eigenvalues index ea93eb303..216a6d51d 100644 --- a/Eigen/Eigenvalues +++ b/Eigen/Eigenvalues @@ -32,6 +32,7 @@ * \endcode */ +#include "src/misc/RealSvd2x2.h" #include "src/Eigenvalues/Tridiagonalization.h" #include "src/Eigenvalues/RealSchur.h" #include "src/Eigenvalues/EigenSolver.h" @@ -31,6 +31,7 @@ * \endcode */ +#include "src/misc/RealSvd2x2.h" #include "src/SVD/UpperBidiagonalization.h" #include "src/SVD/SVDBase.h" #include "src/SVD/JacobiSVD.h" diff --git a/Eigen/src/Eigenvalues/RealQZ.h b/Eigen/src/Eigenvalues/RealQZ.h index a62071d42..c4715b954 100644 --- a/Eigen/src/Eigenvalues/RealQZ.h +++ b/Eigen/src/Eigenvalues/RealQZ.h @@ -552,7 +552,6 @@ namespace Eigen { m_T.coeffRef(l,l-1) = Scalar(0.0); } - template<typename MatrixType> RealQZ<MatrixType>& RealQZ<MatrixType>::compute(const MatrixType& A_in, const MatrixType& B_in, bool computeQZ) { @@ -616,6 +615,37 @@ namespace Eigen { } // check if we converged before reaching iterations limit m_info = (local_iter<m_maxIters) ? Success : NoConvergence; + + // For each non triangular 2x2 diagonal block of S, + // reduce the respective 2x2 diagonal block of T to positive diagonal form using 2x2 SVD. + // This step is not mandatory for QZ, but it does help further extraction of eigenvalues/eigenvectors, + // and is in par with Lapack/Matlab QZ. + if(m_info==Success) + { + for(Index i=0; i<dim-1; ++i) + { + if(m_S.coeff(i+1, i) != Scalar(0)) + { + JacobiRotation<Scalar> j_left, j_right; + internal::real_2x2_jacobi_svd(m_T, i, i+1, &j_left, &j_right); + + // Apply resulting Jacobi rotations + m_T.applyOnTheLeft(i,i+1,j_left); + m_T.applyOnTheRight(i,i+1,j_right); + m_S.applyOnTheLeft(i,i+1,j_left); + m_S.applyOnTheRight(i,i+1,j_right); + m_T(i,i+1) = Scalar(0); + + if(m_computeQZ) { + m_Q.applyOnTheRight(i,i+1,j_left.transpose()); + m_Z.applyOnTheLeft(i,i+1,j_right.transpose()); + } + + i++; + } + } + } + return *this; } // end compute diff --git a/Eigen/src/SVD/JacobiSVD.h b/Eigen/src/SVD/JacobiSVD.h index 1940c8294..b83fd7a4d 100644 --- a/Eigen/src/SVD/JacobiSVD.h +++ b/Eigen/src/SVD/JacobiSVD.h @@ -419,38 +419,6 @@ struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, true> } }; -template<typename MatrixType, typename RealScalar, typename Index> -void real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q, - JacobiRotation<RealScalar> *j_left, - JacobiRotation<RealScalar> *j_right) -{ - using std::sqrt; - using std::abs; - Matrix<RealScalar,2,2> m; - m << numext::real(matrix.coeff(p,p)), numext::real(matrix.coeff(p,q)), - numext::real(matrix.coeff(q,p)), numext::real(matrix.coeff(q,q)); - JacobiRotation<RealScalar> rot1; - RealScalar t = m.coeff(0,0) + m.coeff(1,1); - RealScalar d = m.coeff(1,0) - m.coeff(0,1); - if(d == RealScalar(0)) - { - rot1.s() = RealScalar(0); - rot1.c() = RealScalar(1); - } - else - { - // If d!=0, then t/d cannot overflow because the magnitude of the - // entries forming d are not too small compared to the ones forming t. - RealScalar u = t / d; - RealScalar tmp = sqrt(RealScalar(1) + numext::abs2(u)); - rot1.s() = RealScalar(1) / tmp; - rot1.c() = u / tmp; - } - m.applyOnTheLeft(0,1,rot1); - j_right->makeJacobi(m,0,1); - *j_left = rot1 * j_right->transpose(); -} - template<typename _MatrixType, int QRPreconditioner> struct traits<JacobiSVD<_MatrixType,QRPreconditioner> > { diff --git a/Eigen/src/misc/RealSvd2x2.h b/Eigen/src/misc/RealSvd2x2.h new file mode 100644 index 000000000..cdd7777d2 --- /dev/null +++ b/Eigen/src/misc/RealSvd2x2.h @@ -0,0 +1,54 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009-2010 Benoit Jacob <jacob.benoit.1@gmail.com> +// Copyright (C) 2013-2016 Gael Guennebaud <gael.guennebaud@inria.fr> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#ifndef EIGEN_REALSVD2X2_H +#define EIGEN_REALSVD2X2_H + +namespace Eigen { + +namespace internal { + +template<typename MatrixType, typename RealScalar, typename Index> +void real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q, + JacobiRotation<RealScalar> *j_left, + JacobiRotation<RealScalar> *j_right) +{ + using std::sqrt; + using std::abs; + Matrix<RealScalar,2,2> m; + m << numext::real(matrix.coeff(p,p)), numext::real(matrix.coeff(p,q)), + numext::real(matrix.coeff(q,p)), numext::real(matrix.coeff(q,q)); + JacobiRotation<RealScalar> rot1; + RealScalar t = m.coeff(0,0) + m.coeff(1,1); + RealScalar d = m.coeff(1,0) - m.coeff(0,1); + if(d == RealScalar(0)) + { + rot1.s() = RealScalar(0); + rot1.c() = RealScalar(1); + } + else + { + // If d!=0, then t/d cannot overflow because the magnitude of the + // entries forming d are not too small compared to the ones forming t. + RealScalar u = t / d; + RealScalar tmp = sqrt(RealScalar(1) + numext::abs2(u)); + rot1.s() = RealScalar(1) / tmp; + rot1.c() = u / tmp; + } + m.applyOnTheLeft(0,1,rot1); + j_right->makeJacobi(m,0,1); + *j_left = rot1 * j_right->transpose(); +} + +} // end namespace internal + +} // end namespace Eigen + +#endif // EIGEN_REALSVD2X2_H
\ No newline at end of file |