diff --git a/Eigen/src/Core/Matrix.h b/Eigen/src/Core/Matrix.h index 2fa4b16cf..f7489491f 100644 --- a/Eigen/src/Core/Matrix.h +++ b/Eigen/src/Core/Matrix.h @@ -318,7 +318,7 @@ class Matrix : public MatrixBase class ei_product_eval_to_column_major template struct ei_product_nested_rhs { typedef typename ei_meta_if< - (ei_traits::Flags & NestByValueBit) && !(ei_traits::Flags & RowMajorBit), + (ei_traits::Flags & NestByValueBit) && (!(ei_traits::Flags & RowMajorBit)) && (int(ei_traits::Flags) & DirectAccessBit), T, typename ei_meta_if< ((ei_traits::Flags & EvalBeforeNestingBit) @@ -183,7 +183,7 @@ template struct ei_product_nested_rhs template struct ei_product_nested_lhs { typedef typename ei_meta_if< - ei_traits::Flags & NestByValueBit, + ei_traits::Flags & NestByValueBit && (int(ei_traits::Flags) & DirectAccessBit), T, typename ei_meta_if< int(ei_traits::Flags) & EvalBeforeNestingBit diff --git a/Eigen/src/QR/EigenSolver.h b/Eigen/src/QR/EigenSolver.h index 55584fe0d..cf3ea9c94 100644 --- a/Eigen/src/QR/EigenSolver.h +++ b/Eigen/src/QR/EigenSolver.h @@ -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 class EigenSolver typedef std::complex Complex; typedef Matrix::ret, MatrixType::ColsAtCompileTime, 1> EigenvalueType; typedef Matrix RealVectorType; + typedef Matrix RealVectorTypeX; EigenSolver(const MatrixType& matrix) : m_eivec(matrix.rows(), matrix.cols()), @@ -74,7 +75,7 @@ template 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::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::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(); } 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::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() -= ( (eivali.start(i) * eivalr.start(i).transpose()).lazy() + (eivalr.start(i) * eivali.start(i).transpose()).lazy()); @@ -279,7 +280,7 @@ void EigenSolver::tql2(RealVectorType& eivalr, RealVec Scalar dl1 = eivalr[l+1]; Scalar h = g - eivalr[l]; if (l+2 cdiv(Scalar xr, Scalar xi, Scalar yr, Scalar yi) // Nonsymmetric reduction from Hessenberg to real Schur form. template -void EigenSolver::hqr2(MatrixType& matH, RealVectorType& ort) +void EigenSolver::hqr2(MatrixType& matH) { // This is derived from the Algol procedure hqr2, // by Martin and Wilkinson, Handbook for Auto. Comp., diff --git a/Eigen/src/QR/QR.h b/Eigen/src/QR/QR.h index 782f43a04..bd01cf71e 100644 --- a/Eigen/src/QR/QR.h +++ b/Eigen/src/QR/QR.h @@ -46,7 +46,7 @@ template class QR public: typedef typename MatrixType::Scalar Scalar; - typedef Matrix RMatrixType; + typedef Matrix MatrixTypeR; typedef Matrix VectorType; QR(const MatrixType& matrix) @@ -59,7 +59,7 @@ template 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::_compute(const MatrixType& matrix) /** \returns the matrix R */ template -typename QR::RMatrixType QR::matrixR(void) const +typename QR::MatrixTypeR QR::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(); res.diagonal() = m_norms; return res; } diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index c96578331..6cd98800a 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -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) diff --git a/test/eigensolver.cpp b/test/eigensolver.cpp new file mode 100644 index 000000000..f63d7a535 --- /dev/null +++ b/test/eigensolver.cpp @@ -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 +// +// 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 + +template 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::Real> Complex; + + MatrixType a = MatrixType::random(rows,cols); + MatrixType covMat = a.adjoint() * a; + + EigenSolver eiSymm(covMat); + VERIFY_IS_APPROX(covMat * eiSymm.eigenvectors(), eiSymm.eigenvectors() * eiSymm.eigenvalues().asDiagonal()); + + EigenSolver eiNotSymmButSymm(covMat); + VERIFY_IS_APPROX((covMat.template cast()) * (eiNotSymmButSymm.eigenvectors().template cast()), + (eiNotSymmButSymm.eigenvectors().template cast()) * (eiNotSymmButSymm.eigenvalues().asDiagonal())); + + EigenSolver eiNotSymm(a); +// VERIFY_IS_APPROX(a.template cast() * eiNotSymm.eigenvectors().template cast(), +// eiNotSymm.eigenvectors().template cast() * 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)) ); + } +} diff --git a/test/qr.cpp b/test/qr.cpp new file mode 100644 index 000000000..8facbac96 --- /dev/null +++ b/test/qr.cpp @@ -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 +// +// 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 + +template 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 SquareMatrixType; + typedef Matrix VectorType; + + MatrixType a = MatrixType::random(rows,cols); + QR 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 + } +}