From 4d91229bdce3ab91ef25d853126989e4cd235b9f Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Tue, 1 Sep 2009 16:20:56 +0200 Subject: [PATCH] [mq]: eigensolver --- Eigen/QR | 2 + Eigen/src/QR/ComplexEigenSolver.h | 138 +++++++++++++++ Eigen/src/QR/ComplexSchur.h | 273 ++++++++++++++++++++++++++++++ test/CMakeLists.txt | 1 + test/eigensolver_complex.cpp | 70 ++++++++ 5 files changed, 484 insertions(+) create mode 100644 Eigen/src/QR/ComplexEigenSolver.h create mode 100644 Eigen/src/QR/ComplexSchur.h create mode 100644 test/eigensolver_complex.cpp diff --git a/Eigen/QR b/Eigen/QR index 1cc94d8eb..4b49004c3 100644 --- a/Eigen/QR +++ b/Eigen/QR @@ -42,6 +42,8 @@ namespace Eigen { #include "src/QR/EigenSolver.h" #include "src/QR/SelfAdjointEigenSolver.h" #include "src/QR/HessenbergDecomposition.h" +#include "src/QR/ComplexSchur.h" +#include "src/QR/ComplexEigenSolver.h" // declare all classes for a given matrix type #define EIGEN_QR_MODULE_INSTANTIATE_TYPE(MATRIXTYPE,PREFIX) \ diff --git a/Eigen/src/QR/ComplexEigenSolver.h b/Eigen/src/QR/ComplexEigenSolver.h new file mode 100644 index 000000000..af47c2195 --- /dev/null +++ b/Eigen/src/QR/ComplexEigenSolver.h @@ -0,0 +1,138 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 Claire Maurice +// Copyright (C) 2009 Gael Guennebaud +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see . + +#ifndef EIGEN_COMPLEX_EIGEN_SOLVER_H +#define EIGEN_COMPLEX_EIGEN_SOLVER_H + +#define MAXITER 30 + +template class ComplexEigenSolver +{ + public: + typedef _MatrixType MatrixType; + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits::Real RealScalar; + typedef std::complex Complex; + typedef Matrix EigenvalueType; + typedef Matrix EigenvectorType; + + /** + * \brief Default Constructor. + * + * The default constructor is useful in cases in which the user intends to + * perform decompositions via ComplexEigenSolver::compute(const MatrixType&). + */ + ComplexEigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false) + {} + + ComplexEigenSolver(const MatrixType& matrix) + : m_eivec(matrix.rows(),matrix.cols()), + m_eivalues(matrix.cols()), + m_isInitialized(false) + { + compute(matrix); + } + + EigenvectorType eigenvectors(void) const + { + ei_assert(m_isInitialized && "ComplexEigenSolver is not initialized."); + return m_eivec; + } + + EigenvalueType eigenvalues() const + { + ei_assert(m_isInitialized && "ComplexEigenSolver is not initialized."); + return m_eivalues; + } + + void compute(const MatrixType& matrix); + + protected: + MatrixType m_eivec; + EigenvalueType m_eivalues; + bool m_isInitialized; +}; + + +template +void ComplexEigenSolver::compute(const MatrixType& matrix) +{ + assert(matrix.cols() == matrix.rows()); + int n = matrix.cols(); + m_eivalues.resize(n,1); + + RealScalar eps = epsilon(); + + // Reduce to complex Schur form + ComplexSchur schur(matrix); + + m_eivalues = schur.matrixT().diagonal(); + + m_eivec.setZero(); + + Scalar d2, z; + RealScalar norm = matrix.norm(); + + // compute the (normalized) eigenvectors + for(int k=n-1 ; k>=0 ; k--) + { + d2 = schur.matrixT().coeff(k,k); + m_eivec.coeffRef(k,k) = Scalar(1.0,0.0); + for(int i=k-1 ; i>=0 ; i--) + { + m_eivec.coeffRef(i,k) = -schur.matrixT().coeff(i,k); + if(k-i-1>0) + m_eivec.coeffRef(i,k) -= (schur.matrixT().row(i).segment(i+1,k-i-1) * m_eivec.col(k).segment(i+1,k-i-1)).value(); + z = schur.matrixT().coeff(i,i) - d2; + if(z==Scalar(0)) + z.real() = eps * norm; + m_eivec.coeffRef(i,k) = m_eivec.coeff(i,k) / z; + + } + m_eivec.col(k).normalize(); + } + + m_eivec = schur.matrixU() * m_eivec; + m_isInitialized = true; + + // sort the eigenvalues + { + for (int i=0; i +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see . + +#ifndef EIGEN_COMPLEX_SCHUR_H +#define EIGEN_COMPLEX_SCHUR_H + +#define MAXITER 30 + +/** \ingroup QR + * + * \class ComplexShur + * + * \brief Performs a complex Shur decomposition of a real or complex square matrix + * + */ +template class ComplexSchur +{ + public: + typedef _MatrixType MatrixType; + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits::Real RealScalar; + typedef std::complex Complex; + typedef Matrix ComplexMatrixType; + + /** + * \brief Default Constructor. + * + * The default constructor is useful in cases in which the user intends to + * perform decompositions via ComplexSchur::compute(const MatrixType&). + */ + ComplexSchur() : m_matT(), m_matU(), m_isInitialized(false) + {} + + ComplexSchur(const MatrixType& matrix) + : m_matT(matrix.rows(),matrix.cols()), + m_matU(matrix.rows(),matrix.cols()), + m_isInitialized(false) + { + compute(matrix); + } + + ComplexMatrixType matrixU() const + { + ei_assert(m_isInitialized && "ComplexSchur is not initialized."); + return m_matU; + } + + ComplexMatrixType matrixT() const + { + ei_assert(m_isInitialized && "ComplexShur is not initialized."); + return m_matT; + } + + void compute(const MatrixType& matrix); + + protected: + ComplexMatrixType m_matT, m_matU; + bool m_isInitialized; +}; + +// computes the plane rotation G such that G' x |p| = | c s' |' |p| = |z| +// |q| |-s c' | |q| |0| +// and returns z if requested. Note that the returned c is real. +template void ei_make_givens(const std::complex& p, const std::complex& q, + JacobiRotation >& rot, std::complex* z=0) +{ + typedef std::complex Complex; + T scale, absx, absxy; + if(p==Complex(0)) + { + // return identity + rot.c() = Complex(1,0); + rot.s() = Complex(0,0); + if(z) *z = p; + } + else + { + scale = cnorm1(p); + absx = scale * ei_sqrt(ei_abs2(p/scale)); + scale = ei_abs(scale) + cnorm1(q); + absxy = scale * ei_sqrt((absx/scale)*(absx/scale) + ei_abs2(q/scale)); + rot.c() = Complex(absx / absxy); + Complex np = p/absx; + rot.s() = -ei_conj(np) * q / absxy; + if(z) *z = np * absxy; + } +} + +template +void ComplexSchur::compute(const MatrixType& matrix) +{ + assert(matrix.cols() == matrix.rows()); + int n = matrix.cols(); + + // Reduce to Hessenberg form + HessenbergDecomposition hess(matrix); + + m_matT = hess.matrixH(); + m_matU = hess.matrixQ(); + + int iu = m_matT.cols() - 1; + int il; + RealScalar d,sd,sf; + Complex c,b,disc,r1,r2,kappa; + + RealScalar eps = epsilon(); + + int iter = 0; + while(true) + { + //locate the range in which to iterate + while(iu > 0) + { + d = cnorm1(m_matT.coeffRef(iu,iu)) + cnorm1(m_matT.coeffRef(iu-1,iu-1)); + sd = cnorm1(m_matT.coeffRef(iu,iu-1)); + + if(sd >= eps * d) break; // FIXME : precision criterion ?? + + m_matT.coeffRef(iu,iu-1) = Complex(0); + iter = 0; + --iu; + } + if(iu==0) break; + iter++; + + if(iter >= MAXITER) + { + // FIXME : what to do when iter==MAXITER ?? + std::cerr << "MAXITER" << std::endl; + return; + } + + il = iu-1; + while( il > 0 ) + { + // check if the current 2x2 block on the diagonal is upper triangular + d = cnorm1(m_matT.coeffRef(il,il)) + cnorm1(m_matT.coeffRef(il-1,il-1)); + sd = cnorm1(m_matT.coeffRef(il,il-1)); + + if(sd < eps * d) break; // FIXME : precision criterion ?? + + --il; + } + + if( il != 0 ) m_matT.coeffRef(il,il-1) = Complex(0); + + // compute the shift (the normalization by sf is to avoid under/overflow) + Matrix t = m_matT.template block<2,2>(iu-1,iu-1); + sf = t.cwise().abs().sum(); + t /= sf; + + c = t.determinant(); + b = t.diagonal().sum(); + + disc = csqrt(b*b - RealScalar(4)*c); + + r1 = (b+disc)/RealScalar(2); + r2 = (b-disc)/RealScalar(2); + + if(cnorm1(r1) > cnorm1(r2)) + r2 = c/r1; + else + r1 = c/r2; + + if(cnorm1(r1-t.coeff(1,1)) < cnorm1(r2-t.coeff(1,1))) + kappa = sf * r1; + else + kappa = sf * r2; + + // perform the QR step using Givens rotations + JacobiRotation rot; + ei_make_givens(m_matT.coeff(il,il) - kappa, m_matT.coeff(il+1,il), rot); + + for(int i=il ; i +T cnorm1(const std::complex &Z) +{ + return(ei_abs(Z.real()) + ei_abs(Z.imag())); +} + +/** + * Computes the principal value of the square root of the complex \a z. + */ +template +std::complex csqrt(const std::complex &z) +{ + RealScalar t, tre, tim; + + t = ei_abs(z); + + if (ei_abs(z.real()) <= ei_abs(z.imag())) + { + // No cancellation in these formulas + tre = ei_sqrt(0.5*(t+z.real())); + tim = ei_sqrt(0.5*(t-z.real())); + } + else + { + // Stable computation of the above formulas + if (z.real() > 0) + { + tre = t + z.real(); + tim = ei_abs(z.imag())*ei_sqrt(0.5/tre); + tre = ei_sqrt(0.5*tre); + } + else + { + tim = t - z.real(); + tre = ei_abs(z.imag())*ei_sqrt(0.5/tim); + tim = ei_sqrt(0.5*tim); + } + } + if(z.imag() < 0) + tim = -tim; + + return (std::complex(tre,tim)); + +} + + +#endif // EIGEN_COMPLEX_SCHUR_H diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 896392188..56d36aa70 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -125,6 +125,7 @@ ei_add_test(qr_colpivoting) ei_add_test(qr_fullpivoting) ei_add_test(eigensolver_selfadjoint " " "${GSL_LIBRARIES}") ei_add_test(eigensolver_generic " " "${GSL_LIBRARIES}") +ei_add_test(eigensolver_complex) ei_add_test(svd) ei_add_test(jacobisvd ${EI_OFLAG}) ei_add_test(geo_orthomethods) diff --git a/test/eigensolver_complex.cpp b/test/eigensolver_complex.cpp new file mode 100644 index 000000000..32de327dd --- /dev/null +++ b/test/eigensolver_complex.cpp @@ -0,0 +1,70 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2008-2009 Gael Guennebaud +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see . + +#include "main.h" +#include +#include + +template void eigensolver(const MatrixType& m) +{ + /* this test covers the following files: + ComplexEigenSolver.h + */ + int rows = m.rows(); + int cols = m.cols(); + + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits::Real RealScalar; + typedef Matrix VectorType; + typedef Matrix RealVectorType; + typedef typename std::complex::Real> Complex; + + // RealScalar largerEps = 10*test_precision(); + + MatrixType a = MatrixType::Random(rows,cols); + MatrixType a1 = MatrixType::Random(rows,cols); + MatrixType symmA = a.adjoint() * a + a1.adjoint() * a1; + +// ComplexEigenSolver ei0(symmA); + +// VERIFY_IS_APPROX(symmA * ei0.eigenvectors(), ei0.eigenvectors() * ei0.eigenvalues().asDiagonal()); + +// a.imag().setZero(); +// std::cerr << a << "\n\n"; + ComplexEigenSolver ei1(a); +// exit(1); + VERIFY_IS_APPROX(a * ei1.eigenvectors(), ei1.eigenvectors() * ei1.eigenvalues().asDiagonal()); + +} + +void test_eigensolver_complex() +{ + for(int i = 0; i < g_repeat; i++) { +// CALL_SUBTEST( eigensolver(Matrix4cf()) ); +// CALL_SUBTEST( eigensolver(MatrixXcd(4,4)) ); + CALL_SUBTEST( eigensolver(MatrixXcd(6,6)) ); +// CALL_SUBTEST( eigensolver(MatrixXd(14,14)) ); + } +} +