aboutsummaryrefslogtreecommitdiffhomepage
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Eigen/src/Core/Matrix.h2
-rw-r--r--Eigen/src/Core/ProductWIP.h4
-rw-r--r--Eigen/src/QR/EigenSolver.h25
-rw-r--r--Eigen/src/QR/QR.h8
-rw-r--r--test/CMakeLists.txt2
-rw-r--r--test/eigensolver.cpp61
-rw-r--r--test/qr.cpp55
7 files changed, 138 insertions, 19 deletions
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<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);
}
diff --git a/Eigen/src/Core/ProductWIP.h b/Eigen/src/Core/ProductWIP.h
index 9e7dc91f1..d1bc86a13 100644
--- a/Eigen/src/Core/ProductWIP.h
+++ b/Eigen/src/Core/ProductWIP.h
@@ -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
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<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.,
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<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;
}
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 <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)) );
+ }
+}
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 <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
+ }
+}
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-29 11:27:58 +0000