mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-03-19 18:40:38 +08:00
Code Issues Packages Projects Releases Wiki Activity

* add a HouseholderSequence class (not good enough yet for Triadiagonalization and HessenbergDecomposition)

Browse Source 
* rework a bit AnyMatrixBase, and mobe it to a separate file
This commit is contained in:
Gael Guennebaud 2009-09-16 14:35:42 +02:00
parent 77f858f6ab
commit 49dd5d7847
12 changed files with 339 additions and 103 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
Eigen/Core
View File

@ -143,6 +143,7 @@ namespace Eigen {
#include "src/Core/Functors.h"
#include "src/Core/MatrixBase.h"
#include "src/Core/AnyMatrixBase.h"
#include "src/Core/Coeffs.h"
#ifndef EIGEN_PARSED_BY_DOXYGEN // work around Doxygen bug triggered by Assign.h r814874

1
Eigen/Householder
View File

@ -16,6 +16,7 @@ namespace Eigen {
*/
#include "src/Householder/Householder.h"
#include "src/Householder/HouseholderSequence.h"
} // namespace Eigen

153
Eigen/src/Core/AnyMatrixBase.h Normal file
View File

@ -0,0 +1,153 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
// 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_ANYMATRIXBASE_H
#define EIGEN_ANYMATRIXBASE_H
/** Common base class for all classes T such that MatrixBase has an operator=(T) and a constructor MatrixBase(T).
*
* In other words, an AnyMatrixBase object is an object that can be copied into a MatrixBase.
*
* Besides MatrixBase-derived classes, this also includes special matrix classes such as diagonal matrices, etc.
*
* Notice that this class is trivial, it is only used to disambiguate overloaded functions.
*/
template<typename Derived> struct AnyMatrixBase
{
typedef typename ei_plain_matrix_type<Derived>::type PlainMatrixType;
Derived& derived() { return *static_cast<Derived*>(this); }
const Derived& derived() const { return *static_cast<const Derived*>(this); }
/** \returns the number of rows. \sa cols(), RowsAtCompileTime */
inline int rows() const { return derived().rows(); }
/** \returns the number of columns. \sa rows(), ColsAtCompileTime*/
inline int cols() const { return derived().cols(); }
/** \internal Don't use it, but do the equivalent: \code dst = *this; \endcode */
template<typename Dest> inline void evalTo(Dest& dst) const
{ derived().evalTo(dst); }
/** \internal Don't use it, but do the equivalent: \code dst += *this; \endcode */
template<typename Dest> inline void addToDense(Dest& dst) const
{
// This is the default implementation,
// derived class can reimplement it in a more optimized way.
typename Dest::PlainMatrixType res(rows(),cols());
evalTo(res);
dst += res;
}
/** \internal Don't use it, but do the equivalent: \code dst -= *this; \endcode */
template<typename Dest> inline void subToDense(Dest& dst) const
{
// This is the default implementation,
// derived class can reimplement it in a more optimized way.
typename Dest::PlainMatrixType res(rows(),cols());
evalTo(res);
dst -= res;
}
/** \internal Don't use it, but do the equivalent: \code dst.applyOnTheRight(*this); \endcode */
template<typename Dest> inline void applyThisOnTheRight(Dest& dst) const
{
// This is the default implementation,
// derived class can reimplement it in a more optimized way.
dst = dst * this->derived();
}
/** \internal Don't use it, but do the equivalent: \code dst.applyOnTheLeft(*this); \endcode */
template<typename Dest> inline void applyThisOnTheLeft(Dest& dst) const
{
// This is the default implementation,
// derived class can reimplement it in a more optimized way.
dst = this->derived() * dst;
}
};
/***************************************************************************
* Implementation of matrix base methods
***************************************************************************/
/** Copies the generic expression \a other into *this. \returns a reference to *this.
* The expression must provide a (templated) evalToDense(Derived& dst) const function
* which does the actual job. In practice, this allows any user to write its own
* special matrix without having to modify MatrixBase */
template<typename Derived>
template<typename OtherDerived>
Derived& MatrixBase<Derived>::operator=(const AnyMatrixBase<OtherDerived> &other)
{
other.derived().evalTo(derived());
return derived();
}
template<typename Derived>
template<typename OtherDerived>
Derived& MatrixBase<Derived>::operator+=(const AnyMatrixBase<OtherDerived> &other)
{
other.derived().addToDense(derived());
return derived();
}
template<typename Derived>
template<typename OtherDerived>
Derived& MatrixBase<Derived>::operator-=(const AnyMatrixBase<OtherDerived> &other)
{
other.derived().subToDense(derived());
return derived();
}
/** replaces \c *this by \c *this * \a other.
*
* \returns a reference to \c *this
*/
template<typename Derived>
template<typename OtherDerived>
inline Derived&
MatrixBase<Derived>::operator*=(const AnyMatrixBase<OtherDerived> &other)
{
other.derived().applyThisOnTheRight(derived());
return derived();
}
/** replaces \c *this by \c *this * \a other. It is equivalent to MatrixBase::operator*=() */
template<typename Derived>
template<typename OtherDerived>
inline void MatrixBase<Derived>::applyOnTheRight(const AnyMatrixBase<OtherDerived> &other)
{
other.derived().applyThisOnTheRight(derived());
}
/** replaces \c *this by \c *this * \a other. */
template<typename Derived>
template<typename OtherDerived>
inline void MatrixBase<Derived>::applyOnTheLeft(const AnyMatrixBase<OtherDerived> &other)
{
other.derived().applyThisOnTheLeft(derived());
}
#endif // EIGEN_ANYMATRIXBASE_H

58
Eigen/src/Core/MatrixBase.h
View File

@ -26,44 +26,6 @@
#ifndef EIGEN_MATRIXBASE_H
#define EIGEN_MATRIXBASE_H
/** Common base class for all classes T such that MatrixBase has an operator=(T) and a constructor MatrixBase(T).
*
* In other words, an AnyMatrixBase object is an object that can be copied into a MatrixBase.
*
* Besides MatrixBase-derived classes, this also includes special matrix classes such as diagonal matrices, etc.
*
* Notice that this class is trivial, it is only used to disambiguate overloaded functions.
*/
template<typename Derived> struct AnyMatrixBase
{
typedef typename ei_plain_matrix_type<Derived>::type PlainMatrixType;
Derived& derived() { return *static_cast<Derived*>(this); }
const Derived& derived() const { return *static_cast<const Derived*>(this); }
/** \returns the number of rows. \sa cols(), RowsAtCompileTime */
inline int rows() const { return derived().rows(); }
/** \returns the number of columns. \sa rows(), ColsAtCompileTime*/
inline int cols() const { return derived().cols(); }
template<typename Dest> inline void evalTo(Dest& dst) const
{ derived().evalTo(dst); }
template<typename Dest> inline void addToDense(Dest& dst) const
{
typename Dest::PlainMatrixType res(rows(),cols());
evalToDense(res);
dst += res;
}
template<typename Dest> inline void subToDense(Dest& dst) const
{
typename Dest::PlainMatrixType res(rows(),cols());
evalToDense(res);
dst -= res;
}
};
/** \class MatrixBase
*
* \brief Base class for all matrices, vectors, and expressions
@ -96,7 +58,6 @@ template<typename Derived> class MatrixBase
#endif // not EIGEN_PARSED_BY_DOXYGEN
{
public:
#ifndef EIGEN_PARSED_BY_DOXYGEN
using ei_special_scalar_op_base<Derived,typename ei_traits<Derived>::Scalar,
typename NumTraits<typename ei_traits<Derived>::Scalar>::Real>::operator*;
@ -301,21 +262,14 @@ template<typename Derived> class MatrixBase
*/
Derived& operator=(const MatrixBase& other);
/** Copies the generic expression \a other into *this. \returns a reference to *this.
* The expression must provide a (templated) evalToDense(Derived& dst) const function
* which does the actual job. In practice, this allows any user to write its own
* special matrix without having to modify MatrixBase */
template<typename OtherDerived>
Derived& operator=(const AnyMatrixBase<OtherDerived> &other)
{ other.derived().evalToDense(derived()); return derived(); }
Derived& operator=(const AnyMatrixBase<OtherDerived> &other);
template<typename OtherDerived>
Derived& operator+=(const AnyMatrixBase<OtherDerived> &other)
{ other.derived().addToDense(derived()); return derived(); }
Derived& operator+=(const AnyMatrixBase<OtherDerived> &other);
template<typename OtherDerived>
Derived& operator-=(const AnyMatrixBase<OtherDerived> &other)
{ other.derived().subToDense(derived()); return derived(); }
Derived& operator-=(const AnyMatrixBase<OtherDerived> &other);
template<typename OtherDerived,typename OtherEvalType>
Derived& operator=(const ReturnByValue<OtherDerived,OtherEvalType>& func);
@ -436,6 +390,12 @@ template<typename Derived> class MatrixBase
template<typename OtherDerived>
Derived& operator*=(const AnyMatrixBase<OtherDerived>& other);
template<typename OtherDerived>
void applyOnTheLeft(const AnyMatrixBase<OtherDerived>& other);
template<typename OtherDerived>
void applyOnTheRight(const AnyMatrixBase<OtherDerived>& other);
template<typename DiagonalDerived>
const DiagonalProduct<Derived, DiagonalDerived, DiagonalOnTheRight>
operator*(const DiagonalBase<DiagonalDerived> &diagonal) const;

14
Eigen/src/Core/Product.h
View File

@ -434,18 +434,4 @@ MatrixBase<Derived>::operator*(const MatrixBase<OtherDerived> &other) const
return typename ProductReturnType<Derived,OtherDerived>::Type(derived(), other.derived());
}
/** replaces \c *this by \c *this * \a other.
*
* \returns a reference to \c *this
*/
template<typename Derived>
template<typename OtherDerived>
inline Derived &
MatrixBase<Derived>::operator*=(const AnyMatrixBase<OtherDerived> &other)
{
return derived() = derived() * other.derived();
}
#endif // EIGEN_PRODUCT_H

12
Eigen/src/Core/TriangularMatrix.h
View File

@ -91,9 +91,9 @@ template<typename Derived> class TriangularBase : public AnyMatrixBase<Derived>
#endif // not EIGEN_PARSED_BY_DOXYGEN
template<typename DenseDerived>
void evalToDense(MatrixBase<DenseDerived> &other) const;
void evalTo(MatrixBase<DenseDerived> &other) const;
template<typename DenseDerived>
void evalToDenseLazy(MatrixBase<DenseDerived> &other) const;
void evalToLazy(MatrixBase<DenseDerived> &other) const;
protected:
@ -546,23 +546,23 @@ void TriangularView<MatrixType, Mode>::lazyAssign(const TriangularBase<OtherDeri
* If the matrix is triangular, the opposite part is set to zero. */
template<typename Derived>
template<typename DenseDerived>
void TriangularBase<Derived>::evalToDense(MatrixBase<DenseDerived> &other) const
void TriangularBase<Derived>::evalTo(MatrixBase<DenseDerived> &other) const
{
if(ei_traits<Derived>::Flags & EvalBeforeAssigningBit)
{
typename Derived::PlainMatrixType other_evaluated(rows(), cols());
evalToDenseLazy(other_evaluated);
evalToLazy(other_evaluated);
other.derived().swap(other_evaluated);
}
else
evalToDenseLazy(other.derived());
evalToLazy(other.derived());
}
/** Assigns a triangular or selfadjoint matrix to a dense matrix.
* If the matrix is triangular, the opposite part is set to zero. */
template<typename Derived>
template<typename DenseDerived>
void TriangularBase<Derived>::evalToDenseLazy(MatrixBase<DenseDerived> &other) const
void TriangularBase<Derived>::evalToLazy(MatrixBase<DenseDerived> &other) const
{
const bool unroll = DenseDerived::SizeAtCompileTime * Derived::CoeffReadCost / 2
<= EIGEN_UNROLLING_LIMIT;

1
Eigen/src/Core/util/ForwardDeclarations.h
View File

@ -123,6 +123,7 @@ template<typename MatrixType> class SVD;
template<typename MatrixType, unsigned int Options = 0> class JacobiSVD;
template<typename MatrixType, int UpLo = LowerTriangular> class LLT;
template<typename MatrixType> class LDLT;
template<typename VectorsType, typename CoeffsType> class HouseholderSequence;
template<typename Scalar> class PlanarRotation;
// Geometry module:

146
Eigen/src/Householder/HouseholderSequence.h Normal file
View File

@ -0,0 +1,146 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// 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_HOUSEHOLDER_SEQUENCE_H
#define EIGEN_HOUSEHOLDER_SEQUENCE_H
/** \ingroup Householder_Module
* \householder_module
* \class HouseholderSequence
* \brief Represents a sequence of householder reflections with decreasing size
*
* This class represents a product sequence of householder reflections \f$ H = \Pi_0^{n-1} H_i \f$
* where \f$ H_i \f$ is the i-th householder transformation \f$ I - h_i v_i v_i^* \f$,
* \f$ v_i \f$ is the i-th householder vector \f$ [ 1, m_vectors(i+1,i), m_vectors(i+2,i), ...] \f$
* and \f$ h_i \f$ is the i-th householder coefficient \c m_coeffs[i].
*
* Typical usages are listed below, where H is a HouseholderSequence:
* \code
* A.applyOnTheRight(H); // A = A * H
* A.applyOnTheLeft(H); // A = H * A
* A.applyOnTheRight(H.adjoint()); // A = A * H^*
* A.applyOnTheLeft(H.adjoint()); // A = H^* * A
* MatrixXd Q = H; // conversion to a dense matrix
* \endcode
* In addition to the adjoint, you can also apply the inverse (=adjoint), the transpose, and the conjugate.
*
* \sa MatrixBase::applyOnTheLeft(), MatrixBase::applyOnTheRight()
*/
template<typename VectorsType, typename CoeffsType>
struct ei_traits<HouseholderSequence<VectorsType,CoeffsType> >
{
typedef typename VectorsType::Scalar Scalar;
enum {
RowsAtCompileTime = ei_traits<VectorsType>::RowsAtCompileTime,
ColsAtCompileTime = ei_traits<VectorsType>::RowsAtCompileTime,
MaxRowsAtCompileTime = ei_traits<VectorsType>::MaxRowsAtCompileTime,
MaxColsAtCompileTime = ei_traits<VectorsType>::MaxRowsAtCompileTime,
Flags = 0
};
};
template<typename VectorsType, typename CoeffsType> class HouseholderSequence
: public AnyMatrixBase<HouseholderSequence<VectorsType,CoeffsType> >
{
typedef typename VectorsType::Scalar Scalar;
public:
typedef HouseholderSequence<VectorsType,
typename ei_meta_if<NumTraits<Scalar>::IsComplex,
NestByValue<typename ei_cleantype<typename CoeffsType::ConjugateReturnType>::type >,
CoeffsType>::ret> ConjugateReturnType;
HouseholderSequence(const VectorsType& v, const CoeffsType& h, bool trans = false)
: m_vectors(v), m_coeffs(h), m_trans(trans)
{}
int rows() const { return m_vectors.rows(); }
int cols() const { return m_vectors.rows(); }
HouseholderSequence transpose() const
{ return HouseholderSequence(m_vectors, m_coeffs, !m_trans); }
ConjugateReturnType conjugate() const
{ return ConjugateReturnType(m_vectors, m_coeffs.conjugate(), m_trans); }
ConjugateReturnType adjoint() const
{ return ConjugateReturnType(m_vectors, m_coeffs.conjugate(), !m_trans); }
ConjugateReturnType inverse() const { return adjoint(); }
/** \internal */
template<typename DestType> void evalTo(DestType& dst) const
{
int vecs = std::min(m_vectors.cols(),m_vectors.rows());
int length = m_vectors.rows();
dst.setIdentity();
Matrix<Scalar,1,DestType::RowsAtCompileTime> temp(dst.rows());
for(int k = vecs-1; k >= 0; --k)
{
if(m_trans)
dst.corner(BottomRight, length-k, length-k)
.applyHouseholderOnTheRight(m_vectors.col(k).end(length-k-1), m_coeffs.coeff(k), &temp.coeffRef(0));
else
dst.corner(BottomRight, length-k, length-k)
.applyHouseholderOnTheLeft(m_vectors.col(k).end(length-k-1), m_coeffs.coeff(k), &temp.coeffRef(k));
}
}
/** \internal */
template<typename Dest> inline void applyThisOnTheRight(Dest& dst) const
{
int vecs = std::min(m_vectors.cols(),m_vectors.rows()); // number of householder vectors
int length = m_vectors.rows(); // size of the largest householder vector
Matrix<Scalar,1,Dest::ColsAtCompileTime> temp(dst.rows());
for(int k = 0; k < vecs; ++k)
{
int actual_k = m_trans ? vecs-k-1 : k;
dst.corner(BottomRight, dst.rows(), length-k)
.applyHouseholderOnTheRight(m_vectors.col(k).end(length-k-1), m_coeffs.coeff(k), &temp.coeffRef(0));
}
}
/** \internal */
template<typename Dest> inline void applyThisOnTheLeft(Dest& dst) const
{
int vecs = std::min(m_vectors.cols(),m_vectors.rows()); // number of householder vectors
int length = m_vectors.rows(); // size of the largest householder vector
Matrix<Scalar,1,Dest::ColsAtCompileTime> temp(dst.cols());
for(int k = 0; k < vecs; ++k)
{
int actual_k = m_trans ? k : vecs-k-1;
dst.corner(BottomRight, length-actual_k, dst.cols())
.applyHouseholderOnTheLeft(m_vectors.col(actual_k).end(length-actual_k-1), m_coeffs.coeff(actual_k), &temp.coeffRef(0));
}
}
protected:
typename VectorsType::Nested m_vectors;
typename CoeffsType::Nested m_coeffs;
bool m_trans;
};
#endif // EIGEN_HOUSEHOLDER_SEQUENCE_H

37
Eigen/src/QR/HouseholderQR.h
View File

@ -56,12 +56,13 @@ template<typename MatrixType> class HouseholderQR
Options = MatrixType::Options,
DiagSizeAtCompileTime = EIGEN_ENUM_MIN(ColsAtCompileTime,RowsAtCompileTime)
};
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixQType;
typedef Matrix<Scalar, DiagSizeAtCompileTime, 1> HCoeffsType;
typedef Matrix<Scalar, 1, ColsAtCompileTime> RowVectorType;
typedef typename HouseholderSequence<MatrixQType,HCoeffsType>::ConjugateReturnType HouseholderSequenceType;
/**
* \brief Default Constructor.
@ -97,7 +98,12 @@ template<typename MatrixType> class HouseholderQR
template<typename OtherDerived, typename ResultType>
void solve(const MatrixBase<OtherDerived>& b, ResultType *result) const;
MatrixQType matrixQ(void) const;
MatrixQType matrixQ() const;
HouseholderSequenceType matrixQAsHouseholderSequence() const
{
return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate());
}
/** \returns a reference to the matrix where the Householder QR decomposition is stored
* in a LAPACK-compatible way.
@ -169,7 +175,7 @@ HouseholderQR<MatrixType>& HouseholderQR<MatrixType>::compute(const MatrixType&
int rows = matrix.rows();
int cols = matrix.cols();
int size = std::min(rows,cols);
m_qr = matrix;
m_hCoeffs.resize(size);
@ -206,15 +212,7 @@ void HouseholderQR<MatrixType>::solve(
result->resize(rows, cols);
*result = b;
Matrix<Scalar,1,MatrixType::ColsAtCompileTime> temp(cols);
for (int k = 0; k < cols; ++k)
{
int remainingSize = rows-k;
result->corner(BottomRight, remainingSize, cols)
.applyHouseholderOnTheLeft(m_qr.col(k).end(remainingSize-1), m_hCoeffs.coeff(k), &temp.coeffRef(0));
}
result->applyOnTheLeft(matrixQAsHouseholderSequence().inverse());
const int rank = std::min(result->rows(), result->cols());
m_qr.corner(TopLeft, rank, rank)
@ -227,20 +225,7 @@ template<typename MatrixType>
typename HouseholderQR<MatrixType>::MatrixQType HouseholderQR<MatrixType>::matrixQ() const
{
ei_assert(m_isInitialized && "HouseholderQR is not initialized.");
// compute the product H'_0 H'_1 ... H'_n-1,
// where H_k is the k-th Householder transformation I - h_k v_k v_k'
// and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...]
int rows = m_qr.rows();
int cols = m_qr.cols();
int size = std::min(rows,cols);
MatrixQType res = MatrixQType::Identity(rows, rows);
Matrix<Scalar,1,MatrixType::RowsAtCompileTime> temp(rows);
for (int k = size-1; k >= 0; k--)
{
res.block(k, k, rows-k, rows-k)
.applyHouseholderOnTheLeft(m_qr.col(k).end(rows-k-1), ei_conj(m_hCoeffs.coeff(k)), &temp.coeffRef(k));
}
return res;
return matrixQAsHouseholderSequence();
}
#endif // EIGEN_HIDE_HEAVY_CODE

9
test/householder.cpp
View File

@ -43,7 +43,7 @@ template<typename MatrixType> void householder(const MatrixType& m)
Matrix<Scalar, EIGEN_ENUM_MAX(MatrixType::RowsAtCompileTime,MatrixType::ColsAtCompileTime), 1> _tmp(std::max(rows,cols));
Scalar* tmp = &_tmp.coeffRef(0,0);
Scalar beta;
RealScalar alpha;
EssentialVectorType essential;
@ -58,7 +58,7 @@ template<typename MatrixType> void householder(const MatrixType& m)
v2 = v1;
v1.applyHouseholderOnTheLeft(essential,beta,tmp);
VERIFY_IS_APPROX(v1.norm(), v2.norm());
MatrixType m1(rows, cols),
m2(rows, cols);
@ -72,7 +72,7 @@ template<typename MatrixType> void householder(const MatrixType& m)
VERIFY_IS_MUCH_SMALLER_THAN(m1.block(1,0,rows-1,cols).norm(), m1.norm());
VERIFY_IS_MUCH_SMALLER_THAN(ei_imag(m1(0,0)), ei_real(m1(0,0)));
VERIFY_IS_APPROX(ei_real(m1(0,0)), alpha);
v1 = VectorType::Random(rows);
if(even) v1.end(rows-1).setZero();
SquareMatrixType m3(rows,rows), m4(rows,rows);
@ -84,6 +84,9 @@ template<typename MatrixType> void householder(const MatrixType& m)
VERIFY_IS_MUCH_SMALLER_THAN(m3.block(0,1,rows,rows-1).norm(), m3.norm());
VERIFY_IS_MUCH_SMALLER_THAN(ei_imag(m3(0,0)), ei_real(m3(0,0)));
VERIFY_IS_APPROX(ei_real(m3(0,0)), alpha);
// test householder sequence
// TODO test HouseholderSequence
}
void test_householder()

8
test/jacobisvd.cpp
View File

@ -36,14 +36,14 @@ template<typename MatrixType, unsigned int Options> void svd(const MatrixType& m
RowsAtCompileTime = MatrixType::RowsAtCompileTime,
ColsAtCompileTime = MatrixType::ColsAtCompileTime
};
typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<Scalar>::Real RealScalar;
typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixUType;
typedef Matrix<Scalar, ColsAtCompileTime, ColsAtCompileTime> MatrixVType;
typedef Matrix<Scalar, RowsAtCompileTime, 1> ColVectorType;
typedef Matrix<Scalar, ColsAtCompileTime, 1> InputVectorType;
MatrixType a;
if(pickrandom) a = MatrixType::Random(rows,cols);
else a = m;
@ -53,7 +53,7 @@ template<typename MatrixType, unsigned int Options> void svd(const MatrixType& m
sigma.diagonal() = svd.singularValues().template cast<Scalar>();
MatrixUType u = svd.matrixU();
MatrixVType v = svd.matrixV();
VERIFY_IS_APPROX(a, u * sigma * v.adjoint());
VERIFY_IS_UNITARY(u);
VERIFY_IS_UNITARY(v);
@ -98,7 +98,7 @@ void test_jacobisvd()
}
CALL_SUBTEST(( svd<MatrixXf,0>(MatrixXf(300,200)) ));
CALL_SUBTEST(( svd<MatrixXcd,AtLeastAsManyColsAsRows>(MatrixXcd(100,150)) ));
CALL_SUBTEST(( svd_verify_assert<Matrix3f>() ));
CALL_SUBTEST(( svd_verify_assert<Matrix3d>() ));
CALL_SUBTEST(( svd_verify_assert<MatrixXf>() ));

2
test/qr.cpp
View File

@ -78,7 +78,7 @@ template<typename MatrixType> void qr_invertible()
m3 = MatrixType::Random(size,size);
qr.solve(m3, &m2);
VERIFY_IS_APPROX(m3, m1*m2);
// now construct a matrix with prescribed determinant
m1.setZero();
for(int i = 0; i < size; i++) m1(i,i) = ei_random<Scalar>();