GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-09-17 20:03:17 +08:00
Code Issues Packages Projects Releases Wiki Activity

Fixed wrong line endings.

Browse Source
This commit is contained in:
Hauke Heibel 2009-10-01 07:20:09 +02:00
parent d7a2a37a4c
commit 5409ce1625
1 changed files with 149 additions and 149 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

298
Eigen/src/Eigenvalues/ComplexEigenSolver.h
View File

@ -1,149 +1,149 @@
// This file is part of Eigen, a lightweight C++ template library // This file is part of Eigen, a lightweight C++ template library
// for linear algebra. // for linear algebra.
// //
// Copyright (C) 2009 Claire Maurice // Copyright (C) 2009 Claire Maurice
// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr> // Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr>
// //
// Eigen is free software; you can redistribute it and/or // Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public // modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either // License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version. // version 3 of the License, or (at your option) any later version.
// //
// Alternatively, you can redistribute it and/or // Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as // modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of // published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version. // the License, or (at your option) any later version.
// //
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details. // GNU General Public License for more details.
// //
// You should have received a copy of the GNU Lesser General Public // You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with // License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>. // Eigen. If not, see <http://www.gnu.org/licenses/>.
#ifndef EIGEN_COMPLEX_EIGEN_SOLVER_H #ifndef EIGEN_COMPLEX_EIGEN_SOLVER_H
#define EIGEN_COMPLEX_EIGEN_SOLVER_H #define EIGEN_COMPLEX_EIGEN_SOLVER_H
/** \eigenvalues_module \ingroup Eigenvalues_Module /** \eigenvalues_module \ingroup Eigenvalues_Module
* \nonstableyet * \nonstableyet
* *
* \class ComplexEigenSolver * \class ComplexEigenSolver
* *
* \brief Eigen values/vectors solver for general complex matrices * \brief Eigen values/vectors solver for general complex matrices
* *
* \param MatrixType the type of the matrix of which we are computing the eigen decomposition * \param MatrixType the type of the matrix of which we are computing the eigen decomposition
* *
* \sa class EigenSolver, class SelfAdjointEigenSolver * \sa class EigenSolver, class SelfAdjointEigenSolver
*/ */
template<typename _MatrixType> class ComplexEigenSolver template<typename _MatrixType> class ComplexEigenSolver
{ {
public: public:
typedef _MatrixType MatrixType; typedef _MatrixType MatrixType;
typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<Scalar>::Real RealScalar; typedef typename NumTraits<Scalar>::Real RealScalar;
typedef std::complex<RealScalar> Complex; typedef std::complex<RealScalar> Complex;
typedef Matrix<Complex, MatrixType::ColsAtCompileTime,1> EigenvalueType; typedef Matrix<Complex, MatrixType::ColsAtCompileTime,1> EigenvalueType;
typedef Matrix<Complex, MatrixType::RowsAtCompileTime,MatrixType::ColsAtCompileTime> EigenvectorType; typedef Matrix<Complex, MatrixType::RowsAtCompileTime,MatrixType::ColsAtCompileTime> EigenvectorType;
/** /**
* \brief Default Constructor. * \brief Default Constructor.
* *
* The default constructor is useful in cases in which the user intends to * The default constructor is useful in cases in which the user intends to
* perform decompositions via ComplexEigenSolver::compute(const MatrixType&). * perform decompositions via ComplexEigenSolver::compute(const MatrixType&).
*/ */
ComplexEigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false) ComplexEigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false)
{} {}
ComplexEigenSolver(const MatrixType& matrix) ComplexEigenSolver(const MatrixType& matrix)
: m_eivec(matrix.rows(),matrix.cols()), : m_eivec(matrix.rows(),matrix.cols()),
m_eivalues(matrix.cols()), m_eivalues(matrix.cols()),
m_isInitialized(false) m_isInitialized(false)
{ {
compute(matrix); compute(matrix);
} }
EigenvectorType eigenvectors(void) const EigenvectorType eigenvectors(void) const
{ {
ei_assert(m_isInitialized && "ComplexEigenSolver is not initialized."); ei_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
return m_eivec; return m_eivec;
} }
EigenvalueType eigenvalues() const EigenvalueType eigenvalues() const
{ {
ei_assert(m_isInitialized && "ComplexEigenSolver is not initialized."); ei_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
return m_eivalues; return m_eivalues;
} }
void compute(const MatrixType& matrix); void compute(const MatrixType& matrix);
protected: protected:
MatrixType m_eivec; MatrixType m_eivec;
EigenvalueType m_eivalues; EigenvalueType m_eivalues;
bool m_isInitialized; bool m_isInitialized;
}; };
template<typename MatrixType> template<typename MatrixType>
void ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix) void ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix)
{ {
// this code is inspired from Jampack // this code is inspired from Jampack
assert(matrix.cols() == matrix.rows()); assert(matrix.cols() == matrix.rows());
int n = matrix.cols(); int n = matrix.cols();
m_eivalues.resize(n,1); m_eivalues.resize(n,1);
m_eivec.resize(n,n); m_eivec.resize(n,n);
RealScalar eps = epsilon<RealScalar>(); RealScalar eps = epsilon<RealScalar>();
// Reduce to complex Schur form // Reduce to complex Schur form
ComplexSchur<MatrixType> schur(matrix); ComplexSchur<MatrixType> schur(matrix);
m_eivalues = schur.matrixT().diagonal(); m_eivalues = schur.matrixT().diagonal();
m_eivec.setZero(); m_eivec.setZero();
Scalar d2, z; Scalar d2, z;
RealScalar norm = matrix.norm(); RealScalar norm = matrix.norm();
// compute the (normalized) eigenvectors // compute the (normalized) eigenvectors
for(int k=n-1 ; k>=0 ; k--) for(int k=n-1 ; k>=0 ; k--)
{ {
d2 = schur.matrixT().coeff(k,k); d2 = schur.matrixT().coeff(k,k);
m_eivec.coeffRef(k,k) = Scalar(1.0,0.0); m_eivec.coeffRef(k,k) = Scalar(1.0,0.0);
for(int i=k-1 ; i>=0 ; i--) for(int i=k-1 ; i>=0 ; i--)
{ {
m_eivec.coeffRef(i,k) = -schur.matrixT().coeff(i,k); m_eivec.coeffRef(i,k) = -schur.matrixT().coeff(i,k);
if(k-i-1>0) 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(); 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; z = schur.matrixT().coeff(i,i) - d2;
if(z==Scalar(0)) if(z==Scalar(0))
ei_real_ref(z) = eps * norm; ei_real_ref(z) = eps * norm;
m_eivec.coeffRef(i,k) = m_eivec.coeff(i,k) / z; m_eivec.coeffRef(i,k) = m_eivec.coeff(i,k) / z;
} }
m_eivec.col(k).normalize(); m_eivec.col(k).normalize();
} }
m_eivec = schur.matrixU() * m_eivec; m_eivec = schur.matrixU() * m_eivec;
m_isInitialized = true; m_isInitialized = true;
// sort the eigenvalues // sort the eigenvalues
{ {
for (int i=0; i<n; i++) for (int i=0; i<n; i++)
{ {
int k; int k;
m_eivalues.cwise().abs().end(n-i).minCoeff(&k); m_eivalues.cwise().abs().end(n-i).minCoeff(&k);
if (k != 0) if (k != 0)
{ {
k += i; k += i;
std::swap(m_eivalues[k],m_eivalues[i]); std::swap(m_eivalues[k],m_eivalues[i]);
m_eivec.col(i).swap(m_eivec.col(k)); m_eivec.col(i).swap(m_eivec.col(k));
} }
} }
} }
} }
#endif // EIGEN_COMPLEX_EIGEN_SOLVER_H #endif // EIGEN_COMPLEX_EIGEN_SOLVER_H