new implementation of diagonal matrices and diagonal matrix expressions

This commit is contained in:
Benoit Jacob 2009-06-28 21:27:37 +02:00
parent fc9000f23e
commit 6809f7b1cd
21 changed files with 385 additions and 376 deletions

View File

@ -159,7 +159,6 @@ namespace Eigen {
#include "src/Core/CwiseNullaryOp.h"
#include "src/Core/CwiseUnaryView.h"
#include "src/Core/Dot.h"
#include "src/Core/DiagonalProduct.h"
#include "src/Core/SolveTriangular.h"
#include "src/Core/MapBase.h"
#include "src/Core/Map.h"
@ -168,6 +167,7 @@ namespace Eigen {
#include "src/Core/Transpose.h"
#include "src/Core/DiagonalMatrix.h"
#include "src/Core/Diagonal.h"
#include "src/Core/DiagonalProduct.h"
#include "src/Core/Redux.h"
#include "src/Core/Visitor.h"
#include "src/Core/Fuzzy.h"

View File

@ -436,10 +436,16 @@ struct ei_assign_selector<Derived,OtherDerived,true,true> {
template<typename Derived>
template<typename OtherDerived>
EIGEN_STRONG_INLINE Derived& MatrixBase<Derived>
::operator=(const MatrixBase<OtherDerived>& other)
EIGEN_STRONG_INLINE Derived& MatrixBase<Derived>::operator=(const MatrixBase<OtherDerived>& other)
{
return ei_assign_selector<Derived,OtherDerived>::run(derived(), other.derived());
}
template<typename Derived>
EIGEN_STRONG_INLINE Derived& MatrixBase<Derived>::operator=(const MatrixBase& other)
{
return ei_assign_selector<Derived,Derived>::run(derived(), other.derived());
}
#endif // EIGEN_ASSIGN_H

View File

@ -2,7 +2,7 @@
// for linear algebra.
//
// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr>
// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com>
// Copyright (C) 2007-2009 Benoit Jacob <jacob.benoit.1@gmail.com>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
@ -26,250 +26,164 @@
#ifndef EIGEN_DIAGONALMATRIX_H
#define EIGEN_DIAGONALMATRIX_H
template<typename CoeffsVectorType, typename Derived>
class DiagonalMatrixBase : ei_no_assignment_operator,
public MatrixBase<Derived>
template<typename Derived>
class DiagonalBase : public MultiplierBase<Derived>
{
public:
typedef MatrixBase<Derived> Base;
typedef typename ei_traits<Derived>::Scalar Scalar;
typedef typename Base::PacketScalar PacketScalar;
using Base::derived;
typedef typename ei_cleantype<CoeffsVectorType>::type _CoeffsVectorType;
protected:
// MSVC gets crazy if we define default parameters
template<typename OtherDerived, bool IsVector, bool IsDiagonal> struct construct_from_expression;
// = vector
template<typename OtherDerived>
struct construct_from_expression<OtherDerived,true,false>
{
static void run(Derived& dst, const OtherDerived& src)
{ dst.diagonal() = src; }
};
// = diagonal expression
template<typename OtherDerived, bool IsVector>
struct construct_from_expression<OtherDerived,IsVector,true>
{
static void run(Derived& dst, const OtherDerived& src)
{ dst.diagonal() = src.diagonal(); }
};
/** Default constructor without initialization */
inline DiagonalMatrixBase() {}
/** Constructs a diagonal matrix with given dimension */
inline DiagonalMatrixBase(int dim) : m_coeffs(dim) {}
/** Generic constructor from an expression */
template<typename OtherDerived>
inline DiagonalMatrixBase(const MatrixBase<OtherDerived>& other)
{
construct_from_expression<OtherDerived,OtherDerived::IsVectorAtCompileTime,ei_is_diagonal<OtherDerived>::ret>
::run(derived(),other.derived());
}
typedef typename ei_traits<Derived>::DiagonalVectorType DiagonalVectorType;
typedef typename DiagonalVectorType::Scalar Scalar;
template<typename NewType,int dummy=0>
struct cast_selector {
typedef const DiagonalMatrixWrapper<NestByValue<CwiseUnaryOp<ei_scalar_cast_op<Scalar, NewType>, _CoeffsVectorType> > > return_type;
inline static return_type run(const DiagonalMatrixBase& d) {
return d.m_coeffs.template cast<NewType>().nestByValue().asDiagonal();
}
enum {
RowsAtCompileTime = DiagonalVectorType::SizeAtCompileTime,
ColsAtCompileTime = DiagonalVectorType::SizeAtCompileTime,
MaxRowsAtCompileTime = DiagonalVectorType::MaxSizeAtCompileTime,
MaxColsAtCompileTime = DiagonalVectorType::MaxSizeAtCompileTime,
IsVectorAtCompileTime = 0
};
template<int dummy>
struct cast_selector<Scalar,dummy> {
typedef const Derived& return_type;
inline static return_type run(const DiagonalMatrixBase& d) {
return d.derived();
}
};
typedef Matrix<Scalar, RowsAtCompileTime, ColsAtCompileTime, 0, MaxRowsAtCompileTime, MaxColsAtCompileTime> DenseMatrixType;
#ifndef EIGEN_PARSED_BY_DOXYGEN
inline const Derived& derived() const { return *static_cast<const Derived*>(this); }
inline Derived& derived() { return *static_cast<Derived*>(this); }
#endif // not EIGEN_PARSED_BY_DOXYGEN
DenseMatrixType toDenseMatrix() const { return derived(); }
public:
inline DiagonalMatrixBase(const _CoeffsVectorType& coeffs) : m_coeffs(coeffs)
{
EIGEN_STATIC_ASSERT_VECTOR_ONLY(_CoeffsVectorType);
ei_assert(coeffs.size() > 0);
}
template<typename NewType>
inline typename cast_selector<NewType,0>::return_type
cast() const
{
return cast_selector<NewType,0>::run(*this);
}
/** Assignment operator.
* The right-hand-side \a other must be either a vector representing the diagonal
* coefficients or a diagonal matrix expression.
*/
template<typename OtherDerived>
inline Derived& operator=(const MatrixBase<OtherDerived>& other)
{
construct_from_expression<OtherDerived,OtherDerived::IsVectorAtCompileTime,ei_is_diagonal<OtherDerived>::ret>
::run(derived(),other);
return derived();
}
inline int rows() const { return m_coeffs.size(); }
inline int cols() const { return m_coeffs.size(); }
inline const Scalar coeff(int row, int col) const
{
return row == col ? m_coeffs.coeff(row) : static_cast<Scalar>(0);
}
inline Scalar& coeffRef(int row, int col)
{
ei_assert(row==col);
return m_coeffs.coeffRef(row);
}
inline _CoeffsVectorType& diagonal() { return m_coeffs; }
inline const _CoeffsVectorType& diagonal() const { return m_coeffs; }
protected:
CoeffsVectorType m_coeffs;
inline const DiagonalVectorType& diagonal() const { return derived().diagonal(); }
inline DiagonalVectorType& diagonal() { return derived().diagonal(); }
template<typename MatrixDerived>
const DiagonalProduct<MatrixDerived, Derived, DiagonalOnTheLeft>
operator*(const MatrixBase<MatrixDerived> &matrix) const;
};
template<typename Derived>
template<typename DiagonalDerived>
Derived& MatrixBase<Derived>::operator=(const DiagonalBase<DiagonalDerived> &other)
{
setZero();
diagonal() = other.diagonal();
return derived();
}
/** \class DiagonalMatrix
* \nonstableyet
*
* \brief Represent a diagonal matrix with its storage
* \brief Represents a diagonal matrix with its storage
*
* \param _Scalar the type of coefficients
* \param _Size the dimension of the matrix
* \param _Size the dimension of the matrix, or Dynamic
*
* \sa class Matrix
*/
template<typename _Scalar,int _Size>
struct ei_traits<DiagonalMatrix<_Scalar,_Size> > : ei_traits<Matrix<_Scalar,_Size,_Size> >
template<typename _Scalar, int _Size>
struct ei_traits<DiagonalMatrix<_Scalar,_Size> >
{
enum {
Flags = (ei_traits<Matrix<_Scalar,_Size,_Size> >::Flags & HereditaryBits) | DiagonalBits
};
typedef Matrix<_Scalar,_Size,1> DiagonalVectorType;
};
template<typename _Scalar, int _Size>
class DiagonalMatrix
: public DiagonalMatrixBase<Matrix<_Scalar,_Size,1>, DiagonalMatrix<_Scalar,_Size> >
: public DiagonalBase<DiagonalMatrix<_Scalar,_Size> >
{
public:
EIGEN_GENERIC_PUBLIC_INTERFACE(DiagonalMatrix)
typedef DiagonalMatrixBase<Matrix<_Scalar,_Size,1>, DiagonalMatrix<_Scalar,_Size> > DiagonalBase;
typedef typename ei_traits<DiagonalMatrix>::DiagonalVectorType DiagonalVectorType;
typedef const DiagonalMatrix& Nested;
typedef _Scalar Scalar;
protected:
typedef Matrix<_Scalar,_Size,1> CoeffVectorType;
using DiagonalBase::m_coeffs;
DiagonalVectorType m_diagonal;
public:
inline const DiagonalVectorType& diagonal() const { return m_diagonal; }
inline DiagonalVectorType& diagonal() { return m_diagonal; }
/** Default constructor without initialization */
inline DiagonalMatrix() : DiagonalBase()
{}
inline DiagonalMatrix() {}
/** Constructs a diagonal matrix with given dimension */
inline DiagonalMatrix(int dim) : DiagonalBase(dim)
{}
inline DiagonalMatrix(int dim) : m_diagonal(dim) {}
/** 2D only */
inline DiagonalMatrix(const Scalar& sx, const Scalar& sy)
{
EIGEN_STATIC_ASSERT_MATRIX_SPECIFIC_SIZE(DiagonalMatrix,2,2);
m_coeffs.x() = sx;
m_coeffs.y() = sy;
}
inline DiagonalMatrix(const Scalar& x, const Scalar& y) : m_diagonal(x,y) {}
/** 3D only */
inline DiagonalMatrix(const Scalar& sx, const Scalar& sy, const Scalar& sz)
inline DiagonalMatrix(const Scalar& x, const Scalar& y, const Scalar& z) : m_diagonal(x,y,z) {}
template<typename OtherDerived>
inline DiagonalMatrix(const DiagonalBase<OtherDerived>& other) : m_diagonal(other.diagonal()) {}
/** copy constructor. prevent a default copy constructor from hiding the other templated constructor */
inline DiagonalMatrix(const DiagonalMatrix& other) : m_diagonal(other.diagonal()) {}
/** generic constructor from expression of the diagonal coefficients */
template<typename OtherDerived>
explicit inline DiagonalMatrix(const MatrixBase<OtherDerived>& other) : m_diagonal(other)
{}
template<typename OtherDerived>
DiagonalMatrix& operator=(const DiagonalBase<OtherDerived>& other)
{
EIGEN_STATIC_ASSERT_MATRIX_SPECIFIC_SIZE(DiagonalMatrix,3,3);
m_coeffs.x() = sx;
m_coeffs.y() = sy;
m_coeffs.z() = sz;
m_diagonal = other.diagonal();
return *this;
}
/** copy constructor */
inline DiagonalMatrix(const DiagonalMatrix& other) : DiagonalBase(other.m_coeffs)
{}
/** generic constructor from expression */
template<typename OtherDerived>
explicit inline DiagonalMatrix(const MatrixBase<OtherDerived>& other) : DiagonalBase(other)
{}
/** This is a special case of the templated operator=. Its purpose is to
* prevent a default operator= from hiding the templated operator=.
*/
DiagonalMatrix& operator=(const DiagonalMatrix& other)
{
m_coeffs = other.m_coeffs;
return *this;
}
template<typename OtherDerived>
DiagonalMatrix& operator=(const MatrixBase<OtherDerived>& other)
{
EIGEN_STATIC_ASSERT(ei_is_diagonal<OtherDerived>::ret, THIS_METHOD_IS_ONLY_FOR_DIAGONAL_MATRIX);
m_coeffs = other.diagonal();
m_diagonal = other.m_diagonal();
return *this;
}
inline void resize(int size)
{
m_coeffs.resize(size);
}
inline void resize(int rows, int cols)
{
ei_assert(rows==cols && "a diagonal matrix must be square");
m_coeffs.resize(rows);
}
inline void setZero() { m_coeffs.setZero(); }
inline void resize(int size) { m_diagonal.resize(size); }
inline void setZero() { m_diagonal.setZero(); }
inline void setZero(int size) { m_diagonal.setZero(size); }
inline void setIdentity() { m_diagonal.setIdentity(); }
inline void setIdentity(int size) { m_diagonal.setIdentity(size); }
};
/** \class DiagonalMatrixWrapper
/** \class DiagonalWrapper
* \nonstableyet
*
* \brief Expression of a diagonal matrix
*
* \param CoeffsVectorType the type of the vector of diagonal coefficients
* \param _DiagonalVectorType the type of the vector of diagonal coefficients
*
* This class is an expression of a diagonal matrix with given vector of diagonal
* coefficients. It is the return type of MatrixBase::diagonal(const OtherDerived&)
* and most of the time this is the only way it is used.
* coefficients. It is the return type of MatrixBase::asDiagonal()
* and most of the time this is the only way that it is used.
*
* \sa class DiagonalMatrixBase, class DiagonalMatrix, MatrixBase::asDiagonal()
* \sa class DiagonalMatrix, class DiagonalBase, MatrixBase::asDiagonal()
*/
template<typename CoeffsVectorType>
struct ei_traits<DiagonalMatrixWrapper<CoeffsVectorType> >
template<typename _DiagonalVectorType>
struct ei_traits<DiagonalWrapper<_DiagonalVectorType> >
{
typedef typename CoeffsVectorType::Scalar Scalar;
typedef typename ei_nested<CoeffsVectorType>::type CoeffsVectorTypeNested;
typedef typename ei_unref<CoeffsVectorTypeNested>::type _CoeffsVectorTypeNested;
enum {
RowsAtCompileTime = CoeffsVectorType::SizeAtCompileTime,
ColsAtCompileTime = CoeffsVectorType::SizeAtCompileTime,
MaxRowsAtCompileTime = CoeffsVectorType::MaxSizeAtCompileTime,
MaxColsAtCompileTime = CoeffsVectorType::MaxSizeAtCompileTime,
Flags = (_CoeffsVectorTypeNested::Flags & HereditaryBits) | DiagonalBits,
CoeffReadCost = _CoeffsVectorTypeNested::CoeffReadCost
};
typedef _DiagonalVectorType DiagonalVectorType;
};
template<typename CoeffsVectorType>
class DiagonalMatrixWrapper
: public DiagonalMatrixBase<typename CoeffsVectorType::Nested, DiagonalMatrixWrapper<CoeffsVectorType> >
template<typename _DiagonalVectorType>
class DiagonalWrapper
: public DiagonalBase<DiagonalWrapper<_DiagonalVectorType> >, ei_no_assignment_operator
{
typedef typename CoeffsVectorType::Nested CoeffsVectorTypeNested;
typedef DiagonalMatrixBase<CoeffsVectorTypeNested, DiagonalMatrixWrapper<CoeffsVectorType> > DiagonalBase;
public:
EIGEN_GENERIC_PUBLIC_INTERFACE(DiagonalMatrixWrapper)
inline DiagonalMatrixWrapper(const CoeffsVectorType& coeffs) : DiagonalBase(coeffs)
{}
typedef _DiagonalVectorType DiagonalVectorType;
typedef DiagonalWrapper Nested;
inline DiagonalWrapper(const DiagonalVectorType& diagonal) : m_diagonal(diagonal) {}
const DiagonalVectorType& diagonal() const { return m_diagonal; }
protected:
const typename DiagonalVectorType::Nested m_diagonal;
};
/** \nonstableyet
* \returns an expression of a diagonal matrix with *this as vector of diagonal coefficients
* \returns a pseudo-expression of a diagonal matrix with *this as vector of diagonal coefficients
*
* \only_for_vectors
*
@ -278,10 +192,10 @@ class DiagonalMatrixWrapper
* Example: \include MatrixBase_asDiagonal.cpp
* Output: \verbinclude MatrixBase_asDiagonal.out
*
* \sa class DiagonalMatrix, isDiagonal()
* \sa class DiagonalWrapper, class DiagonalMatrix, diagonal(), isDiagonal()
**/
template<typename Derived>
inline const DiagonalMatrixWrapper<Derived>
inline const DiagonalWrapper<Derived>
MatrixBase<Derived>::asDiagonal() const
{
return derived();

View File

@ -2,7 +2,7 @@
// for linear algebra.
//
// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com>
// Copyright (C) 2007-2009 Benoit Jacob <jacob.benoit.1@gmail.com>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
@ -26,105 +26,66 @@
#ifndef EIGEN_DIAGONALPRODUCT_H
#define EIGEN_DIAGONALPRODUCT_H
/** \internal Specialization of ei_nested for DiagonalMatrix.
* Unlike ei_nested, if the argument is a DiagonalMatrix and if it must be evaluated,
* then it evaluated to a DiagonalMatrix having its own argument evaluated.
*/
template<typename T, int N, bool IsDiagonal = ei_is_diagonal<T>::ret> struct ei_nested_diagonal : ei_nested<T,N> {};
template<typename T, int N> struct ei_nested_diagonal<T,N,true>
: ei_nested<T, N, DiagonalMatrix<typename T::Scalar, EIGEN_ENUM_MIN(T::RowsAtCompileTime,T::ColsAtCompileTime)> >
{};
// specialization of ProductReturnType
template<typename Lhs, typename Rhs>
struct ProductReturnType<Lhs,Rhs,DiagonalProduct>
template<typename MatrixType, typename DiagonalType, int Order>
struct ei_traits<DiagonalProduct<MatrixType, DiagonalType, Order> >
{
typedef typename ei_nested_diagonal<Lhs,Rhs::ColsAtCompileTime>::type LhsNested;
typedef typename ei_nested_diagonal<Rhs,Lhs::RowsAtCompileTime>::type RhsNested;
typedef Product<LhsNested, RhsNested, DiagonalProduct> Type;
};
template<typename LhsNested, typename RhsNested>
struct ei_traits<Product<LhsNested, RhsNested, DiagonalProduct> >
{
// clean the nested types:
typedef typename ei_cleantype<LhsNested>::type _LhsNested;
typedef typename ei_cleantype<RhsNested>::type _RhsNested;
typedef typename _LhsNested::Scalar Scalar;
typedef typename MatrixType::Scalar Scalar;
enum {
LhsFlags = _LhsNested::Flags,
RhsFlags = _RhsNested::Flags,
RowsAtCompileTime = _LhsNested::RowsAtCompileTime,
ColsAtCompileTime = _RhsNested::ColsAtCompileTime,
MaxRowsAtCompileTime = _LhsNested::MaxRowsAtCompileTime,
MaxColsAtCompileTime = _RhsNested::MaxColsAtCompileTime,
LhsIsDiagonal = ei_is_diagonal<_LhsNested>::ret,
RhsIsDiagonal = ei_is_diagonal<_RhsNested>::ret,
CanVectorizeRhs = (!RhsIsDiagonal) && (RhsFlags & RowMajorBit) && (RhsFlags & PacketAccessBit)
&& (ColsAtCompileTime % ei_packet_traits<Scalar>::size == 0),
CanVectorizeLhs = (!LhsIsDiagonal) && (!(LhsFlags & RowMajorBit)) && (LhsFlags & PacketAccessBit)
&& (RowsAtCompileTime % ei_packet_traits<Scalar>::size == 0),
RemovedBits = ~((RhsFlags & RowMajorBit) && (!CanVectorizeLhs) ? 0 : RowMajorBit),
Flags = ((unsigned int)(LhsFlags | RhsFlags) & HereditaryBits & RemovedBits)
| (((CanVectorizeLhs&&RhsIsDiagonal) || (CanVectorizeRhs&&LhsIsDiagonal)) ? PacketAccessBit : 0),
CoeffReadCost = NumTraits<Scalar>::MulCost + _LhsNested::CoeffReadCost + _RhsNested::CoeffReadCost
RowsAtCompileTime = MatrixType::RowsAtCompileTime,
ColsAtCompileTime = MatrixType::ColsAtCompileTime,
MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
Flags = (unsigned int)(MatrixType::Flags) & HereditaryBits,
CoeffReadCost = NumTraits<Scalar>::MulCost + MatrixType::CoeffReadCost + DiagonalType::DiagonalVectorType::CoeffReadCost
};
};
template<typename LhsNested, typename RhsNested> class Product<LhsNested, RhsNested, DiagonalProduct> : ei_no_assignment_operator,
public MatrixBase<Product<LhsNested, RhsNested, DiagonalProduct> >
template<typename MatrixType, typename DiagonalType, int Order>
class DiagonalProduct : ei_no_assignment_operator,
public MatrixBase<DiagonalProduct<MatrixType, DiagonalType, Order> >
{
typedef typename ei_traits<Product>::_LhsNested _LhsNested;
typedef typename ei_traits<Product>::_RhsNested _RhsNested;
enum {
RhsIsDiagonal = ei_is_diagonal<_RhsNested>::ret
};
public:
EIGEN_GENERIC_PUBLIC_INTERFACE(Product)
EIGEN_GENERIC_PUBLIC_INTERFACE(DiagonalProduct)
template<typename Lhs, typename Rhs>
inline Product(const Lhs& lhs, const Rhs& rhs)
: m_lhs(lhs), m_rhs(rhs)
inline DiagonalProduct(const MatrixType& matrix, const DiagonalType& diagonal)
: m_matrix(matrix), m_diagonal(diagonal)
{
ei_assert(lhs.cols() == rhs.rows());
ei_assert(diagonal.diagonal().size() == (Order == DiagonalOnTheLeft ? matrix.rows() : matrix.cols()));
}
inline int rows() const { return m_lhs.rows(); }
inline int cols() const { return m_rhs.cols(); }
inline int rows() const { return m_matrix.rows(); }
inline int cols() const { return m_matrix.cols(); }
const Scalar coeff(int row, int col) const
{
const int unique = RhsIsDiagonal ? col : row;
return m_lhs.coeff(row, unique) * m_rhs.coeff(unique, col);
}
template<int LoadMode>
const PacketScalar packet(int row, int col) const
{
if (RhsIsDiagonal)
{
return ei_pmul(m_lhs.template packet<LoadMode>(row, col), ei_pset1(m_rhs.coeff(col, col)));
}
else
{
return ei_pmul(ei_pset1(m_lhs.coeff(row, row)), m_rhs.template packet<LoadMode>(row, col));
}
return m_diagonal.diagonal().coeff(Order == DiagonalOnTheLeft ? row : col) * m_matrix.coeff(row, col);
}
protected:
const LhsNested m_lhs;
const RhsNested m_rhs;
const typename MatrixType::Nested m_matrix;
const typename DiagonalType::Nested m_diagonal;
};
/** \returns the diagonal matrix product of \c *this by the diagonal matrix \a diagonal.
*/
template<typename Derived>
template<typename DiagonalDerived>
inline const DiagonalProduct<Derived, DiagonalDerived, DiagonalOnTheRight>
MatrixBase<Derived>::operator*(const DiagonalBase<DiagonalDerived> &diagonal) const
{
return DiagonalProduct<Derived, DiagonalDerived, DiagonalOnTheRight>(derived(), diagonal.derived());
}
/** \returns the diagonal matrix product of \c *this by the matrix \a matrix.
*/
template<typename DiagonalDerived>
template<typename MatrixDerived>
inline const DiagonalProduct<MatrixDerived, DiagonalDerived, DiagonalOnTheLeft>
DiagonalBase<DiagonalDerived>::operator*(const MatrixBase<MatrixDerived> &matrix) const
{
return DiagonalProduct<MatrixDerived, DiagonalDerived, DiagonalOnTheLeft>(matrix.derived(), derived());
}
#endif // EIGEN_DIAGONALPRODUCT_H

View File

@ -171,11 +171,7 @@ template<typename Derived> class MapBase
return Base::operator=(other);
}
template<typename OtherDerived>
Derived& operator=(const MatrixBase<OtherDerived>& other)
{
return Base::operator=(other);
}
using Base::operator=;
using Base::operator*=;

View File

@ -124,6 +124,7 @@ class Matrix
{
public:
EIGEN_GENERIC_PUBLIC_INTERFACE(Matrix)
enum { Options = _Options };
friend class Eigen::Map<Matrix, Unaligned>;
typedef class Eigen::Map<Matrix, Unaligned> UnalignedMapType;
@ -335,11 +336,11 @@ class Matrix
EIGEN_STRONG_INLINE Matrix& operator=(const ReturnByValue<OtherDerived,OtherEvalType>& func)
{ return Base::operator=(func); }
EIGEN_INHERIT_ASSIGNMENT_OPERATOR(Matrix, +=)
EIGEN_INHERIT_ASSIGNMENT_OPERATOR(Matrix, -=)
EIGEN_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(Matrix, *=)
EIGEN_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(Matrix, /=)
using Base::operator +=;
using Base::operator -=;
using Base::operator *=;
using Base::operator /=;
/** Default constructor.
*
* For fixed-size matrices, does nothing.
@ -438,6 +439,22 @@ class Matrix
{ other.evalTo(*this); }
/** Destructor */
inline ~Matrix() {}
template<typename DiagonalDerived>
EIGEN_STRONG_INLINE Matrix& operator=(const DiagonalBase<DiagonalDerived> &other)
{
resize(other.diagonal().size(), other.diagonal().size());
Base::operator=(other);
return *this;
}
template<typename DiagonalDerived>
EIGEN_STRONG_INLINE Matrix(const DiagonalBase<DiagonalDerived> &other)
: m_storage(other.diagonal().size() * other.diagonal().size(), other.diagonal().size(), other.diagonal().size())
{
*this = other;
}
/** Override MatrixBase::swap() since for dynamic-sized matrices of same type it is enough to swap the
* data pointers.

View File

@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com>
// Copyright (C) 2006-2009 Benoit Jacob <jacob.benoit.1@gmail.com>
// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
//
// Eigen is free software; you can redistribute it and/or
@ -26,6 +26,32 @@
#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
{
Derived& derived() { return *static_cast<Derived*>(this); }
const Derived& derived() const { return *static_cast<const Derived*>(this); }
};
/** Common base class for all classes T such that there are overloaded operator* allowing to
* multiply a MatrixBase by a T on both sides.
*
* In other words, an AnyMatrixBase object is an object that can be multiplied a MatrixBase, the result being again a MatrixBase.
*
* Besides MatrixBase-derived classes, this also includes certain special matrix classes, such as diagonal matrices.
*/
template<typename Derived> struct MultiplierBase : public AnyMatrixBase<Derived>
{
using AnyMatrixBase<Derived>::derived;
};
/** \class MatrixBase
*
* \brief Base class for all matrices, vectors, and expressions
@ -50,11 +76,11 @@
cout << x.row(0) << endl;
}
* \endcode
*
*/
template<typename Derived> class MatrixBase
#ifndef EIGEN_PARSED_BY_DOXYGEN
: public ei_special_scalar_op_base<Derived,typename ei_traits<Derived>::Scalar,
: public MultiplierBase<Derived>,
public ei_special_scalar_op_base<Derived,typename ei_traits<Derived>::Scalar,
typename NumTraits<typename ei_traits<Derived>::Scalar>::Real>
#endif // not EIGEN_PARSED_BY_DOXYGEN
{
@ -256,10 +282,7 @@ template<typename Derived> class MatrixBase
/** Special case of the template operator=, in order to prevent the compiler
* from generating a default operator= (issue hit with g++ 4.1)
*/
inline Derived& operator=(const MatrixBase& other)
{
return this->operator=<Derived>(other);
}
Derived& operator=(const MatrixBase& other);
#ifndef EIGEN_PARSED_BY_DOXYGEN
/** Copies \a other into *this without evaluating other. \returns a reference to *this. */
@ -358,17 +381,30 @@ template<typename Derived> class MatrixBase
operator*(const Scalar& scalar, const MatrixBase& matrix)
{ return matrix*scalar; }
template<typename OtherDerived>
const typename ProductReturnType<Derived,OtherDerived>::Type
operator*(const MatrixBase<OtherDerived> &other) const;
/** replaces \c *this by \c *this * \a other.
*
* \returns a reference to \c *this
*/
template<typename OtherDerived>
Derived& operator*=(const MatrixBase<OtherDerived>& other);
Derived& operator*=(const MultiplierBase<OtherDerived>& other)
{
return *this = *this * other.derived();
}
template<typename DiagonalDerived>
const DiagonalProduct<Derived, DiagonalDerived, DiagonalOnTheRight>
operator*(const DiagonalBase<DiagonalDerived> &diagonal) const;
template<typename DiagonalDerived>
Derived& operator=(const DiagonalBase<DiagonalDerived> &other);
template<typename OtherDerived>
typename ei_plain_matrix_type_column_major<OtherDerived>::type
solveTriangular(const MatrixBase<OtherDerived>& other) const;
solveTriangular(const MatrixBase<OtherDerived>& other) const;
template<typename OtherDerived>
void solveTriangularInPlace(const MatrixBase<OtherDerived>& other) const;
@ -477,7 +513,7 @@ template<typename Derived> class MatrixBase
static const BasisReturnType UnitZ();
static const BasisReturnType UnitW();
const DiagonalMatrixWrapper<Derived> asDiagonal() const;
const DiagonalWrapper<Derived> asDiagonal() const;
void fill(const Scalar& value);
Derived& setConstant(const Scalar& value);
@ -588,8 +624,7 @@ template<typename Derived> class MatrixBase
void visit(Visitor& func) const;
#ifndef EIGEN_PARSED_BY_DOXYGEN
inline const Derived& derived() const { return *static_cast<const Derived*>(this); }
inline Derived& derived() { return *static_cast<Derived*>(this); }
using MultiplierBase<Derived>::derived;
inline Derived& const_cast_derived() const
{ return *static_cast<Derived*>(const_cast<MatrixBase*>(this)); }
#endif // not EIGEN_PARSED_BY_DOXYGEN
@ -685,16 +720,16 @@ template<typename Derived> class MatrixBase
/////////// Sparse module ///////////
// dense = spasre * dense
// dense = sparse * dense
template<typename Derived1, typename Derived2>
Derived& lazyAssign(const SparseProduct<Derived1,Derived2,SparseTimeDenseProduct>& product);
// dense = dense * spasre
// dense = dense * sparse
template<typename Derived1, typename Derived2>
Derived& lazyAssign(const SparseProduct<Derived1,Derived2,DenseTimeSparseProduct>& product);
template<typename OtherDerived,typename OtherEvalType>
Derived& operator=(const ReturnByValue<OtherDerived,OtherEvalType>& func);
#ifdef EIGEN_MATRIXBASE_PLUGIN
#include EIGEN_MATRIXBASE_PLUGIN
#endif

View File

@ -47,10 +47,10 @@ struct ei_product_packet_impl;
* This class defines the typename Type representing the optimized product expression
* between two matrix expressions. In practice, using ProductReturnType<Lhs,Rhs>::Type
* is the recommended way to define the result type of a function returning an expression
* which involve a matrix product. The class Product or DiagonalProduct should never be
* which involve a matrix product. The class Product should never be
* used directly.
*
* \sa class Product, class DiagonalProduct, MatrixBase::operator*(const MatrixBase<OtherDerived>&)
* \sa class Product, MatrixBase::operator*(const MatrixBase<OtherDerived>&)
*/
template<typename Lhs, typename Rhs, int ProductMode>
struct ProductReturnType
@ -62,7 +62,6 @@ struct ProductReturnType
};
// cache friendly specialization
// note that there is a DiagonalProduct specialization in DiagonalProduct.h
template<typename Lhs, typename Rhs>
struct ProductReturnType<Lhs,Rhs,CacheFriendlyProduct>
{
@ -78,15 +77,12 @@ struct ProductReturnType<Lhs,Rhs,CacheFriendlyProduct>
/* Helper class to determine the type of the product, can be either:
* - NormalProduct
* - CacheFriendlyProduct
* - DiagonalProduct
*/
template<typename Lhs, typename Rhs> struct ei_product_mode
{
enum{
value = ei_is_diagonal<Rhs>::ret || ei_is_diagonal<Lhs>::ret
? DiagonalProduct
: Lhs::MaxColsAtCompileTime == Dynamic
value = Lhs::MaxColsAtCompileTime == Dynamic
&& ( Lhs::MaxRowsAtCompileTime == Dynamic
|| Rhs::MaxColsAtCompileTime == Dynamic )
&& (!(Rhs::IsVectorAtCompileTime && (Lhs::Flags&RowMajorBit) && (!(Lhs::Flags&DirectAccessBit))))
@ -290,18 +286,6 @@ 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 MatrixBase<OtherDerived> &other)
{
return *this = *this * other;
}
/***************************************************************************
* Normal product .coeff() implementation (with meta-unrolling)
***************************************************************************/

View File

@ -205,12 +205,14 @@ template<typename T> struct ei_is_diagonal
};
};
enum { DiagonalOnTheLeft, DiagonalOnTheRight };
enum { Aligned, Unaligned };
enum { ForceAligned, AsRequested };
enum { ConditionalJumpCost = 5 };
enum CornerType { TopLeft, TopRight, BottomLeft, BottomRight };
enum DirectionType { Vertical, Horizontal, BothDirections };
enum ProductEvaluationMode { NormalProduct, CacheFriendlyProduct, DiagonalProduct, SparseTimeSparseProduct, SparseTimeDenseProduct, DenseTimeSparseProduct };
enum ProductEvaluationMode { NormalProduct, CacheFriendlyProduct, SparseTimeSparseProduct, SparseTimeDenseProduct, DenseTimeSparseProduct };
enum {
/** \internal Equivalent to a slice vectorization for fixed-size matrices having good alignment

View File

@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com>
// Copyright (C) 2007-2009 Benoit Jacob <jacob.benoit.1@gmail.com>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
@ -28,6 +28,9 @@
template<typename T> struct ei_traits;
template<typename T> struct NumTraits;
template<typename Derived> struct AnyMatrixBase;
template<typename Derived> struct MultiplierBase;
template<typename _Scalar, int _Rows, int _Cols,
int _Options = EIGEN_DEFAULT_MATRIX_STORAGE_ORDER_OPTION | AutoAlign,
int _MaxRows = _Rows, int _MaxCols = _Cols> class Matrix;
@ -46,10 +49,13 @@ template<typename UnaryOp, typename MatrixType> class CwiseUnaryOp;
template<typename ViewOp, typename MatrixType> class CwiseUnaryView;
template<typename BinaryOp, typename Lhs, typename Rhs> class CwiseBinaryOp;
template<typename Lhs, typename Rhs, int ProductMode> class Product;
template<typename CoeffsVectorType, typename Derived> class DiagonalMatrixBase;
template<typename CoeffsVectorType> class DiagonalMatrixWrapper;
template<typename Derived> class DiagonalBase;
template<typename _DiagonalVectorType> class DiagonalWrapper;
template<typename _Scalar, int _Size> class DiagonalMatrix;
template<typename MatrixType, typename DiagonalType, int Order> class DiagonalProduct;
template<typename MatrixType, int Index> class Diagonal;
template<typename MatrixType, int PacketAccess = AsRequested> class Map;
template<typename MatrixType, unsigned int Mode> class Part;
template<typename MatrixType, unsigned int Mode> class Extract;

View File

@ -233,30 +233,16 @@ using Eigen::ei_cos;
// needed to define it here as escaping characters in CMake add_definition's argument seems very problematic.
#define EIGEN_DOCS_IO_FORMAT IOFormat(3, AlignCols, " ", "\n", "", "")
#define EIGEN_INHERIT_ASSIGNMENT_OPERATOR(Derived, Op) \
template<typename OtherDerived> \
EIGEN_STRONG_INLINE Derived& operator Op(const Eigen::MatrixBase<OtherDerived>& other) \
{ \
return Base::operator Op(other.derived()); \
} \
EIGEN_STRONG_INLINE Derived& operator Op(const Derived& other) \
{ \
return Base::operator Op(other); \
}
#define EIGEN_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(Derived, Op) \
template<typename Other> \
EIGEN_STRONG_INLINE Derived& operator Op(const Other& scalar) \
{ \
return Base::operator Op(scalar); \
}
#define EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Derived) \
EIGEN_INHERIT_ASSIGNMENT_OPERATOR(Derived, =) \
EIGEN_INHERIT_ASSIGNMENT_OPERATOR(Derived, +=) \
EIGEN_INHERIT_ASSIGNMENT_OPERATOR(Derived, -=) \
EIGEN_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(Derived, *=) \
EIGEN_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(Derived, /=)
using Base::operator =; \
using Base::operator +=; \
using Base::operator -=; \
using Base::operator *=; \
using Base::operator /=; \
EIGEN_STRONG_INLINE Derived& operator=(const Derived& other) \
{ \
return Base::operator=(other); \
}
#define _EIGEN_GENERIC_PUBLIC_INTERFACE(Derived, BaseClass) \
typedef BaseClass Base; \

View File

@ -91,12 +91,12 @@ class RotationBase
*/
template<typename OtherDerived>
inline typename generic_product_selector<OtherDerived,OtherDerived::IsVectorAtCompileTime>::ReturnType
operator*(const MatrixBase<OtherDerived>& e) const
operator*(const MultiplierBase<OtherDerived>& e) const
{ return generic_product_selector<OtherDerived>::run(derived(), e.derived()); }
/** \returns the concatenation of a linear transformation \a l with the rotation \a r */
template<typename OtherDerived> friend
inline RotationMatrixType operator*(const MatrixBase<OtherDerived>& l, const Derived& r)
inline RotationMatrixType operator*(const MultiplierBase<OtherDerived>& l, const Derived& r)
{ return l.derived() * r.toRotationMatrix(); }
/** \returns the concatenation of the rotation \c *this with a transformation \a t */

View File

@ -141,7 +141,7 @@ static inline DiagonalMatrix<Scalar,3> Scaling(Scalar sx, Scalar sy, Scalar sz)
* This is an alias for coeffs.asDiagonal()
*/
template<typename Derived>
static inline const DiagonalMatrixWrapper<Derived> Scaling(const MatrixBase<Derived>& coeffs)
static inline const DiagonalWrapper<Derived> Scaling(const MatrixBase<Derived>& coeffs)
{ return coeffs.asDiagonal(); }
/** \addtogroup Geometry_Module */

View File

@ -226,14 +226,14 @@ public:
/** Constructs and initializes a transformation from a Dim^2 or a (Dim+1)^2 matrix. */
template<typename OtherDerived>
inline explicit Transform(const MatrixBase<OtherDerived>& other)
inline explicit Transform(const AnyMatrixBase<OtherDerived>& other)
{
ei_transform_construct_from_matrix<OtherDerived,Mode,Dim,HDim>::run(this, other.derived());
}
/** Set \c *this from a Dim^2 or (Dim+1)^2 matrix. */
template<typename OtherDerived>
inline Transform& operator=(const MatrixBase<OtherDerived>& other)
inline Transform& operator=(const AnyMatrixBase<OtherDerived>& other)
{
ei_transform_construct_from_matrix<OtherDerived,Mode,Dim,HDim>::run(this, other.derived());
return *this;
@ -310,7 +310,7 @@ public:
// note: this function is defined here because some compilers cannot find the respective declaration
template<typename OtherDerived>
inline const typename ei_transform_right_product_impl<OtherDerived,Mode,_Dim,_Dim+1>::ResultType
operator * (const MatrixBase<OtherDerived> &other) const
operator * (const MultiplierBase<OtherDerived> &other) const
{ return ei_transform_right_product_impl<OtherDerived,Mode,Dim,HDim>::run(*this,other.derived()); }
/** \returns the product expression of a transformation matrix \a a times a transform \a b
@ -322,11 +322,11 @@ public:
*/
template<typename OtherDerived> friend
inline const typename ei_transform_left_product_impl<OtherDerived,Mode,_Dim,_Dim+1>::ResultType
operator * (const MatrixBase<OtherDerived> &a, const Transform &b)
operator * (const MultiplierBase<OtherDerived> &a, const Transform &b)
{ return ei_transform_left_product_impl<OtherDerived,Mode,Dim,HDim>::run(a.derived(),b); }
template<typename OtherDerived>
inline Transform& operator*=(const MatrixBase<OtherDerived>& other) { return *this = *this * other; }
inline Transform& operator*=(const MultiplierBase<OtherDerived>& other) { return *this = *this * other; }
/** Contatenates two transformations */
inline const Transform operator * (const Transform& other) const
@ -944,7 +944,7 @@ struct ei_transform_take_affine_part<Transform<Scalar,Dim,AffineCompact> > {
template<typename Other, int Mode, int Dim, int HDim>
struct ei_transform_construct_from_matrix<Other, Mode,Dim,HDim, Dim,Dim>
{
static inline void run(Transform<typename ei_traits<Other>::Scalar,Dim,Mode> *transform, const Other& other)
static inline void run(Transform<typename Other::Scalar,Dim,Mode> *transform, const Other& other)
{
transform->linear() = other;
transform->translation().setZero();
@ -955,7 +955,7 @@ struct ei_transform_construct_from_matrix<Other, Mode,Dim,HDim, Dim,Dim>
template<typename Other, int Mode, int Dim, int HDim>
struct ei_transform_construct_from_matrix<Other, Mode,Dim,HDim, Dim,HDim>
{
static inline void run(Transform<typename ei_traits<Other>::Scalar,Dim,Mode> *transform, const Other& other)
static inline void run(Transform<typename Other::Scalar,Dim,Mode> *transform, const Other& other)
{
transform->affine() = other;
transform->makeAffine();
@ -965,20 +965,32 @@ struct ei_transform_construct_from_matrix<Other, Mode,Dim,HDim, Dim,HDim>
template<typename Other, int Mode, int Dim, int HDim>
struct ei_transform_construct_from_matrix<Other, Mode,Dim,HDim, HDim,HDim>
{
static inline void run(Transform<typename ei_traits<Other>::Scalar,Dim,Mode> *transform, const Other& other)
static inline void run(Transform<typename Other::Scalar,Dim,Mode> *transform, const Other& other)
{ transform->matrix() = other; }
};
template<typename Other, int Dim, int HDim>
struct ei_transform_construct_from_matrix<Other, AffineCompact,Dim,HDim, HDim,HDim>
{
static inline void run(Transform<typename ei_traits<Other>::Scalar,Dim,AffineCompact> *transform, const Other& other)
static inline void run(Transform<typename Other::Scalar,Dim,AffineCompact> *transform, const Other& other)
{ transform->matrix() = other.template block<Dim,HDim>(0,0); }
};
/*****************************************************
*** Specializations of operator* with a MatrixBase ***
*****************************************************/
/*********************************************************
*** Specializations of operator* with a MultiplierBase ***
*********************************************************/
// ei_general_product_return_type is a generalization of ProductReturnType, for all types (including e.g. DiagonalBase...),
// instead of being restricted to MatrixBase.
template<typename Lhs, typename Rhs> struct ei_general_product_return_type;
template<typename D1, typename D2> struct ei_general_product_return_type<MatrixBase<D1>, MatrixBase<D2> >
: ProductReturnType<D1,D2> {};
template<typename Lhs, typename D2> struct ei_general_product_return_type<Lhs, MatrixBase<D2> >
{ typedef D2 Type; };
template<typename D1, typename Rhs> struct ei_general_product_return_type<MatrixBase<D1>, Rhs >
{ typedef D1 Type; };
// Projective * set of homogeneous column vectors
template<typename Other, int Dim, int HDim>

View File

@ -93,7 +93,7 @@ public:
/** Concatenates a translation and a linear transformation */
template<typename OtherDerived>
inline AffineTransformType operator* (const MatrixBase<OtherDerived>& linear) const;
inline AffineTransformType operator* (const MultiplierBase<OtherDerived>& linear) const;
/** Concatenates a translation and a rotation */
template<typename Derived>
@ -103,7 +103,7 @@ public:
/** \returns the concatenation of a linear transformation \a l with the translation \a t */
// its a nightmare to define a templated friend function outside its declaration
template<typename OtherDerived> friend
inline AffineTransformType operator*(const MatrixBase<OtherDerived>& linear, const Translation& t)
inline AffineTransformType operator*(const MultiplierBase<OtherDerived>& linear, const Translation& t)
{
AffineTransformType res;
res.matrix().setZero();
@ -182,7 +182,7 @@ Translation<Scalar,Dim>::operator* (const UniformScaling<Scalar>& other) const
template<typename Scalar, int Dim>
template<typename OtherDerived>
inline typename Translation<Scalar,Dim>::AffineTransformType
Translation<Scalar,Dim>::operator* (const MatrixBase<OtherDerived>& linear) const
Translation<Scalar,Dim>::operator* (const MultiplierBase<OtherDerived>& linear) const
{
AffineTransformType res;
res.matrix().setZero();

View File

@ -98,6 +98,7 @@ ei_add_test(redux)
ei_add_test(product_small)
ei_add_test(product_large ${EI_OFLAG})
ei_add_test(product_selfadjoint)
ei_add_test(diagonalmatrices)
ei_add_test(adjoint)
ei_add_test(submatrices)
ei_add_test(miscmatrices)

95
test/diagonalmatrices.cpp Normal file
View File

@ -0,0 +1,95 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
//
// 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"
template<typename MatrixType> void diagonalmatrices(const MatrixType& m)
{
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
enum { Rows = MatrixType::RowsAtCompileTime, Cols = MatrixType::ColsAtCompileTime };
typedef Matrix<Scalar, Rows, 1> VectorType;
typedef Matrix<Scalar, 1, Cols> RowVectorType;
typedef Matrix<Scalar, Rows, Rows> SquareMatrixType;
typedef DiagonalMatrix<Scalar, Rows> LeftDiagonalMatrix;
typedef DiagonalMatrix<Scalar, Cols> RightDiagonalMatrix;
int rows = m.rows();
int cols = m.cols();
MatrixType m1 = MatrixType::Random(rows, cols),
m2 = MatrixType::Random(rows, cols);
VectorType v1 = VectorType::Random(rows),
v2 = VectorType::Random(rows);
RowVectorType rv1 = RowVectorType::Random(cols),
rv2 = RowVectorType::Random(cols);
LeftDiagonalMatrix ldm1(v1), ldm2(v2);
RightDiagonalMatrix rdm1(rv1), rdm2(rv2);
int i = ei_random<int>(0, rows-1);
int j = ei_random<int>(0, cols-1);
VERIFY_IS_APPROX( ((ldm1 * m1)(i,j)) , ldm1.diagonal()(i) * m1(i,j) );
VERIFY_IS_APPROX( ((ldm1 * (m1+m2))(i,j)) , ldm1.diagonal()(i) * (m1+m2)(i,j) );
VERIFY_IS_APPROX( ((m1 * rdm1)(i,j)) , rdm1.diagonal()(j) * m1(i,j) );
VERIFY_IS_APPROX( ((v1.asDiagonal() * m1)(i,j)) , v1(i) * m1(i,j) );
VERIFY_IS_APPROX( ((m1 * rv1.asDiagonal())(i,j)) , rv1(j) * m1(i,j) );
VERIFY_IS_APPROX( (((v1+v2).asDiagonal() * m1)(i,j)) , (v1+v2)(i) * m1(i,j) );
VERIFY_IS_APPROX( (((v1+v2).asDiagonal() * (m1+m2))(i,j)) , (v1+v2)(i) * (m1+m2)(i,j) );
VERIFY_IS_APPROX( ((m1 * (rv1+rv2).asDiagonal())(i,j)) , (rv1+rv2)(j) * m1(i,j) );
VERIFY_IS_APPROX( (((m1+m2) * (rv1+rv2).asDiagonal())(i,j)) , (rv1+rv2)(j) * (m1+m2)(i,j) );
SquareMatrixType sq_m1 (v1.asDiagonal());
VERIFY_IS_APPROX(sq_m1, v1.asDiagonal().toDenseMatrix());
sq_m1 = v1.asDiagonal();
VERIFY_IS_APPROX(sq_m1, v1.asDiagonal().toDenseMatrix());
SquareMatrixType sq_m2 = v1.asDiagonal();
VERIFY_IS_APPROX(sq_m1, sq_m2);
ldm1 = v1.asDiagonal();
LeftDiagonalMatrix ldm3(v1);
VERIFY_IS_APPROX(ldm1.diagonal(), ldm3.diagonal());
LeftDiagonalMatrix ldm4 = v1.asDiagonal();
VERIFY_IS_APPROX(ldm1.diagonal(), ldm4.diagonal());
sq_m1.block(0,0,rows,rows) = ldm1;
VERIFY_IS_APPROX(sq_m1, ldm1.toDenseMatrix());
sq_m1.transpose() = ldm1;
VERIFY_IS_APPROX(sq_m1, ldm1.toDenseMatrix());
}
void test_diagonalmatrices()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST( diagonalmatrices(Matrix<float, 1, 1>()) );
CALL_SUBTEST( diagonalmatrices(Matrix3f()) );
CALL_SUBTEST( diagonalmatrices(Matrix<double,3,3,RowMajor>()) );
CALL_SUBTEST( diagonalmatrices(Matrix4d()) );
CALL_SUBTEST( diagonalmatrices(Matrix<float,4,4,RowMajor>()) );
CALL_SUBTEST( diagonalmatrices(MatrixXcf(3, 5)) );
CALL_SUBTEST( diagonalmatrices(MatrixXi(10, 8)) );
CALL_SUBTEST( diagonalmatrices(Matrix<double,Dynamic,Dynamic,RowMajor>(20, 20)) );
CALL_SUBTEST( diagonalmatrices(MatrixXf(21, 24)) );
}
}

View File

@ -57,7 +57,7 @@ template<typename MatrixType> void eigensolver(const MatrixType& m)
EigenSolver<MatrixType> ei1(a);
VERIFY_IS_APPROX(a * ei1.pseudoEigenvectors(), ei1.pseudoEigenvectors() * ei1.pseudoEigenvalueMatrix());
VERIFY_IS_APPROX(a.template cast<Complex>() * ei1.eigenvectors(),
ei1.eigenvectors() * ei1.eigenvalues().asDiagonal().eval());
ei1.eigenvectors() * ei1.eigenvalues().asDiagonal());
}

View File

@ -75,7 +75,7 @@ template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m)
convert(gEvec, _evec);
// test gsl itself !
VERIFY((symmA * _evec).isApprox(_evec * _eval.asDiagonal().eval(), largerEps));
VERIFY((symmA * _evec).isApprox(_evec * _eval.asDiagonal(), largerEps));
// compare with eigen
VERIFY_IS_APPROX(_eval, eiSymm.eigenvalues());
@ -86,7 +86,7 @@ template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m)
convert(gEval, _eval);
convert(gEvec, _evec);
// test GSL itself:
VERIFY((symmA * _evec).isApprox(symmB * (_evec * _eval.asDiagonal().eval()), largerEps));
VERIFY((symmA * _evec).isApprox(symmB * (_evec * _eval.asDiagonal()), largerEps));
// compare with eigen
// std::cerr << _eval.transpose() << "\n" << eiSymmGen.eigenvalues().transpose() << "\n\n";
@ -102,11 +102,11 @@ template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m)
#endif
VERIFY((symmA * eiSymm.eigenvectors()).isApprox(
eiSymm.eigenvectors() * eiSymm.eigenvalues().asDiagonal().eval(), largerEps));
eiSymm.eigenvectors() * eiSymm.eigenvalues().asDiagonal(), largerEps));
// generalized eigen problem Ax = lBx
VERIFY((symmA * eiSymmGen.eigenvectors()).isApprox(
symmB * (eiSymmGen.eigenvectors() * eiSymmGen.eigenvalues().asDiagonal().eval()), largerEps));
symmB * (eiSymmGen.eigenvectors() * eiSymmGen.eigenvalues().asDiagonal()), largerEps));
MatrixType sqrtSymmA = eiSymm.operatorSqrt();
VERIFY_IS_APPROX(symmA, sqrtSymmA*sqrtSymmA);

View File

@ -332,12 +332,6 @@ template<typename Scalar, int Mode> void transformations(void)
Translation<double,3> tr1d = tr1.template cast<double>();
VERIFY_IS_APPROX(tr1d.template cast<Scalar>(),tr1);
AlignedScaling3 sc1(v0);
DiagonalMatrix<float,3> sc1f; sc1f = sc1.template cast<float>();
VERIFY_IS_APPROX(sc1f.template cast<Scalar>(),sc1);
DiagonalMatrix<double,3> sc1d; sc1d = (sc1.template cast<double>());
VERIFY_IS_APPROX(sc1d.template cast<Scalar>(),sc1);
AngleAxis<float> aa1f = aa1.template cast<float>();
VERIFY_IS_APPROX(aa1f.template cast<Scalar>(),aa1);
AngleAxis<double> aa1d = aa1.template cast<double>();

View File

@ -43,7 +43,7 @@ template<typename MatrixType> void miscMatrices(const MatrixType& m)
VectorType v1 = VectorType::Random(rows);
v1[0];
Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime>
square = v1.asDiagonal();
square(v1.asDiagonal());
if(r==r2) VERIFY_IS_APPROX(square(r,r2), v1[r]);
else VERIFY_IS_MUCH_SMALLER_THAN(square(r,r2), static_cast<Scalar>(1));
square = MatrixType::Zero(rows, rows);