diff options
author | Hauke Heibel <hauke.heibel@gmail.com> | 2009-10-01 07:20:09 +0200 |
---|---|---|
committer | Hauke Heibel <hauke.heibel@gmail.com> | 2009-10-01 07:20:09 +0200 |
commit | 5409ce1625431f9f8149e66028db17e8515ab166 (patch) | |
tree | 7c3dabfd1cce1f73977a8783c514572f42ed48fd /Eigen/src/Eigenvalues/ComplexEigenSolver.h | |
parent | d7a2a37a4c1728dbddbfae9d577e2e49f5edda65 (diff) |
Fixed wrong line endings.
Diffstat (limited to 'Eigen/src/Eigenvalues/ComplexEigenSolver.h')
-rw-r--r-- | Eigen/src/Eigenvalues/ComplexEigenSolver.h | 298 |
1 files changed, 149 insertions, 149 deletions
diff --git a/Eigen/src/Eigenvalues/ComplexEigenSolver.h b/Eigen/src/Eigenvalues/ComplexEigenSolver.h index 224918795..86206ce79 100644 --- a/Eigen/src/Eigenvalues/ComplexEigenSolver.h +++ b/Eigen/src/Eigenvalues/ComplexEigenSolver.h @@ -1,149 +1,149 @@ -// This file is part of Eigen, a lightweight C++ template library
-// for linear algebra.
-//
-// Copyright (C) 2009 Claire Maurice
-// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr>
-//
-// 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 <http://www.gnu.org/licenses/>.
-
-#ifndef EIGEN_COMPLEX_EIGEN_SOLVER_H
-#define EIGEN_COMPLEX_EIGEN_SOLVER_H
-
-/** \eigenvalues_module \ingroup Eigenvalues_Module
- * \nonstableyet
- *
- * \class ComplexEigenSolver
- *
- * \brief Eigen values/vectors solver for general complex matrices
- *
- * \param MatrixType the type of the matrix of which we are computing the eigen decomposition
- *
- * \sa class EigenSolver, class SelfAdjointEigenSolver
- */
-template<typename _MatrixType> class ComplexEigenSolver
-{
- public:
- typedef _MatrixType MatrixType;
- typedef typename MatrixType::Scalar Scalar;
- typedef typename NumTraits<Scalar>::Real RealScalar;
- typedef std::complex<RealScalar> Complex;
- typedef Matrix<Complex, MatrixType::ColsAtCompileTime,1> EigenvalueType;
- typedef Matrix<Complex, MatrixType::RowsAtCompileTime,MatrixType::ColsAtCompileTime> 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<typename MatrixType>
-void ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix)
-{
- // this code is inspired from Jampack
- assert(matrix.cols() == matrix.rows());
- int n = matrix.cols();
- m_eivalues.resize(n,1);
- m_eivec.resize(n,n);
-
- RealScalar eps = epsilon<RealScalar>();
-
- // Reduce to complex Schur form
- ComplexSchur<MatrixType> 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))
- ei_real_ref(z) = 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<n; i++)
- {
- int k;
- m_eivalues.cwise().abs().end(n-i).minCoeff(&k);
- if (k != 0)
- {
- k += i;
- std::swap(m_eivalues[k],m_eivalues[i]);
- m_eivec.col(i).swap(m_eivec.col(k));
- }
- }
- }
-}
-
-
-
-#endif // EIGEN_COMPLEX_EIGEN_SOLVER_H
+// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 Claire Maurice +// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr> +// +// 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 <http://www.gnu.org/licenses/>. + +#ifndef EIGEN_COMPLEX_EIGEN_SOLVER_H +#define EIGEN_COMPLEX_EIGEN_SOLVER_H + +/** \eigenvalues_module \ingroup Eigenvalues_Module + * \nonstableyet + * + * \class ComplexEigenSolver + * + * \brief Eigen values/vectors solver for general complex matrices + * + * \param MatrixType the type of the matrix of which we are computing the eigen decomposition + * + * \sa class EigenSolver, class SelfAdjointEigenSolver + */ +template<typename _MatrixType> class ComplexEigenSolver +{ + public: + typedef _MatrixType MatrixType; + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; + typedef std::complex<RealScalar> Complex; + typedef Matrix<Complex, MatrixType::ColsAtCompileTime,1> EigenvalueType; + typedef Matrix<Complex, MatrixType::RowsAtCompileTime,MatrixType::ColsAtCompileTime> 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<typename MatrixType> +void ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix) +{ + // this code is inspired from Jampack + assert(matrix.cols() == matrix.rows()); + int n = matrix.cols(); + m_eivalues.resize(n,1); + m_eivec.resize(n,n); + + RealScalar eps = epsilon<RealScalar>(); + + // Reduce to complex Schur form + ComplexSchur<MatrixType> 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)) + ei_real_ref(z) = 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<n; i++) + { + int k; + m_eivalues.cwise().abs().end(n-i).minCoeff(&k); + if (k != 0) + { + k += i; + std::swap(m_eivalues[k],m_eivalues[i]); + m_eivec.col(i).swap(m_eivec.col(k)); + } + } + } +} + + + +#endif // EIGEN_COMPLEX_EIGEN_SOLVER_H |