mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-01-24 14:45:14 +08:00
Code Issues Packages Projects Releases Wiki Activity

* fix the QR module to use extract/part instead of the previous triangular stuff

Browse Source 
* added qr and eigensolver tests
* fix a compilation warning in Matrix copy constructor
This commit is contained in:
Gael Guennebaud 2008-05-27 09:16:27 +00:00
parent 5aa00f6870
commit 559233c73e
7 changed files with 138 additions and 19 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
Eigen/src/Core/Matrix.h
View File

@ -318,7 +318,7 @@ class Matrix : public MatrixBase<Matrix<_Scalar, _Rows, _Cols, _Flags, _MaxRows,
}
/** Copy constructor */
inline Matrix(const Matrix& other)
: m_storage(other.rows() * other.cols(), other.rows(), other.cols())
: Base(), m_storage(other.rows() * other.cols(), other.rows(), other.cols())
{
Base::lazyAssign(other);
}

4
Eigen/src/Core/ProductWIP.h
View File

@ -167,7 +167,7 @@ template<typename T> class ei_product_eval_to_column_major
template<typename T, int n=1> struct ei_product_nested_rhs
{
typedef typename ei_meta_if<
(ei_traits<T>::Flags & NestByValueBit) && !(ei_traits<T>::Flags & RowMajorBit),
(ei_traits<T>::Flags & NestByValueBit) && (!(ei_traits<T>::Flags & RowMajorBit)) && (int(ei_traits<T>::Flags) & DirectAccessBit),
T,
typename ei_meta_if<
((ei_traits<T>::Flags & EvalBeforeNestingBit)
@ -183,7 +183,7 @@ template<typename T, int n=1> struct ei_product_nested_rhs
template<typename T, int n=1> struct ei_product_nested_lhs
{
typedef typename ei_meta_if<
ei_traits<T>::Flags & NestByValueBit,
ei_traits<T>::Flags & NestByValueBit && (int(ei_traits<T>::Flags) & DirectAccessBit),
T,
typename ei_meta_if<
int(ei_traits<T>::Flags) & EvalBeforeNestingBit

25
Eigen/src/QR/EigenSolver.h
View File

@ -32,7 +32,7 @@
* \param MatrixType the type of the matrix of which we are computing the eigen decomposition
* \param IsSelfadjoint tells the input matrix is guaranteed to be selfadjoint (hermitian). In that case the
* return type of eigenvalues() is a real vector.
*
*
* Currently it only support real matrices.
*
* \note this code was adapted from JAMA (public domain)
@ -49,6 +49,7 @@ template<typename _MatrixType, bool IsSelfadjoint=false> class EigenSolver
typedef std::complex<RealScalar> Complex;
typedef Matrix<typename ei_meta_if<IsSelfadjoint, Scalar, Complex>::ret, MatrixType::ColsAtCompileTime, 1> EigenvalueType;
typedef Matrix<RealScalar, MatrixType::ColsAtCompileTime, 1> RealVectorType;
typedef Matrix<RealScalar, Dynamic, 1> RealVectorTypeX;
EigenSolver(const MatrixType& matrix)
: m_eivec(matrix.rows(), matrix.cols()),
@ -74,7 +75,7 @@ template<typename _MatrixType, bool IsSelfadjoint=false> class EigenSolver
void tql2(RealVectorType& eivalr, RealVectorType& eivali);
void orthes(MatrixType& matH, RealVectorType& ort);
void hqr2(MatrixType& matH, RealVectorType& ort);
void hqr2(MatrixType& matH);
protected:
MatrixType m_eivec;
@ -87,7 +88,7 @@ void EigenSolver<MatrixType,IsSelfadjoint>::computeImpl(const MatrixType& matrix
assert(matrix.cols() == matrix.rows());
int n = matrix.cols();
m_eivalues.resize(n,1);
RealVectorType eivali(n);
m_eivec = matrix;
@ -115,25 +116,25 @@ void EigenSolver<MatrixType,IsSelfadjoint>::computeImpl(const MatrixType& matrix
RealVectorType eivalr(n);
RealVectorType eivali(n);
m_eivec = matrix;
// Tridiagonalize.
tridiagonalization(eivalr, eivali);
// Diagonalize.
tql2(eivalr, eivali);
m_eivalues = eivalr.template cast<Complex>();
}
else
{
MatrixType matH = matrix;
RealVectorType ort(n);
// Reduce to Hessenberg form.
orthes(matH, ort);
// Reduce Hessenberg to real Schur form.
hqr2(matH, ort);
hqr2(matH);
}
}
@ -198,7 +199,7 @@ void EigenSolver<MatrixType,IsSelfadjoint>::tridiagonalization(RealVectorType& e
f = (eivali.start(i).transpose() * eivalr.start(i))(0,0);
eivali.start(i) = (eivali.start(i) - (f / (h + h)) * eivalr.start(i))/h;
m_eivec.corner(TopLeft, i, i).lower() -=
m_eivec.corner(TopLeft, i, i).template part<Lower>() -=
( (eivali.start(i) * eivalr.start(i).transpose()).lazy()
+ (eivalr.start(i) * eivali.start(i).transpose()).lazy());
@ -279,7 +280,7 @@ void EigenSolver<MatrixType,IsSelfadjoint>::tql2(RealVectorType& eivalr, RealVec
Scalar dl1 = eivalr[l+1];
Scalar h = g - eivalr[l];
if (l+2<n)
eivalr.end(n-l-2) -= RealVectorType::constant(n-l-2, h);
eivalr.end(n-l-2) -= RealVectorTypeX::constant(n-l-2, h);
f = f + h;
// Implicit QL transformation.
@ -432,7 +433,7 @@ std::complex<Scalar> cdiv(Scalar xr, Scalar xi, Scalar yr, Scalar yi)
// Nonsymmetric reduction from Hessenberg to real Schur form.
template<typename MatrixType, bool IsSelfadjoint>
void EigenSolver<MatrixType,IsSelfadjoint>::hqr2(MatrixType& matH, RealVectorType& ort)
void EigenSolver<MatrixType,IsSelfadjoint>::hqr2(MatrixType& matH)
{
// This is derived from the Algol procedure hqr2,
// by Martin and Wilkinson, Handbook for Auto. Comp.,

8
Eigen/src/QR/QR.h
View File

@ -46,7 +46,7 @@ template<typename MatrixType> class QR
public:
typedef typename MatrixType::Scalar Scalar;
typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime> RMatrixType;
typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime> MatrixTypeR;
typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType;
QR(const MatrixType& matrix)
@ -59,7 +59,7 @@ template<typename MatrixType> class QR
/** \returns whether or not the matrix is of full rank */
bool isFullRank() const { return ei_isMuchSmallerThan(m_norms.cwiseAbs().minCoeff(), Scalar(1)); }
RMatrixType matrixR(void) const;
MatrixTypeR matrixR(void) const;
MatrixType matrixQ(void) const;
@ -108,10 +108,10 @@ void QR<MatrixType>::_compute(const MatrixType& matrix)
/** \returns the matrix R */
template<typename MatrixType>
typename QR<MatrixType>::RMatrixType QR<MatrixType>::matrixR(void) const
typename QR<MatrixType>::MatrixTypeR QR<MatrixType>::matrixR(void) const
{
int cols = m_qr.cols();
RMatrixType res = m_qr.block(0,0,cols,cols).strictlyUpper();
MatrixTypeR res = m_qr.block(0,0,cols,cols).template extract<StrictlyUpper>();
res.diagonal() = m_norms;
return res;
}

2
test/CMakeLists.txt
View File

@ -75,5 +75,7 @@ EI_ADD_TEST(product)
EI_ADD_TEST(triangular)
EI_ADD_TEST(cholesky)
# EI_ADD_TEST(determinant)
EI_ADD_TEST(qr)
EI_ADD_TEST(eigensolver)
ENDIF(BUILD_TESTS)

61
test/eigensolver.cpp Normal file
View File

@ -0,0 +1,61 @@
// 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 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/>.
#include "main.h"
#include <Eigen/QR>
template<typename MatrixType> void eigensolver(const MatrixType& m)
{
/* this test covers the following files:
EigenSolver.h
*/
int rows = m.rows();
int cols = m.cols();
typedef typename std::complex<typename NumTraits<typename MatrixType::Scalar>::Real> Complex;
MatrixType a = MatrixType::random(rows,cols);
MatrixType covMat = a.adjoint() * a;
EigenSolver<MatrixType,true> eiSymm(covMat);
VERIFY_IS_APPROX(covMat * eiSymm.eigenvectors(), eiSymm.eigenvectors() * eiSymm.eigenvalues().asDiagonal());
EigenSolver<MatrixType,false> eiNotSymmButSymm(covMat);
VERIFY_IS_APPROX((covMat.template cast<Complex>()) * (eiNotSymmButSymm.eigenvectors().template cast<Complex>()),
(eiNotSymmButSymm.eigenvectors().template cast<Complex>()) * (eiNotSymmButSymm.eigenvalues().asDiagonal()));
EigenSolver<MatrixType,false> eiNotSymm(a);
// VERIFY_IS_APPROX(a.template cast<Complex>() * eiNotSymm.eigenvectors().template cast<Complex>(),
// eiNotSymm.eigenvectors().template cast<Complex>() * eiNotSymm.eigenvalues().asDiagonal());
}
void test_eigensolver()
{
for(int i = 0; i < 1; i++) {
CALL_SUBTEST( eigensolver(Matrix3f()) );
CALL_SUBTEST( eigensolver(Matrix4d()) );
CALL_SUBTEST( eigensolver(MatrixXd(7,7)) );
}
}

55
test/qr.cpp Normal file
View File

@ -0,0 +1,55 @@
// 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 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/>.
#include "main.h"
#include <Eigen/QR>
template<typename MatrixType> void qr(const MatrixType& m)
{
/* this test covers the following files:
QR.h
*/
int rows = m.rows();
int cols = m.cols();
typedef typename MatrixType::Scalar Scalar;
typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime> SquareMatrixType;
typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType;
MatrixType a = MatrixType::random(rows,cols);
QR<MatrixType> qrOfA(a);
VERIFY_IS_APPROX(a, qrOfA.matrixQ() * qrOfA.matrixR());
VERIFY_IS_NOT_APPROX(a+MatrixType::identity(rows, cols), qrOfA.matrixQ() * qrOfA.matrixR());
}
void test_qr()
{
for(int i = 0; i < 1; i++) {
CALL_SUBTEST( qr(Matrix2f()) );
CALL_SUBTEST( qr(Matrix3d()) );
CALL_SUBTEST( qr(MatrixXf(12,8)) );
// CALL_SUBTEST( qr(MatrixXcd(17,7)) ); // complex numbers are not supported yet
}
}