mirror of
https://gitlab.com/libeigen/eigen.git
synced 2025-04-19 16:19:37 +08:00
- introduce Part and Extract classes, splitting and extending the former
Triangular class - full meta-unrolling in Part - move inverseProduct() to MatrixBase - compilation fix in ProductWIP: introduce a meta-selector to only do direct access on types that support it. - phase out the old Product, remove the WIP_DIRTY stuff. - misc renaming and fixes
This commit is contained in:
Benoit Jacob
parent
8f1fc80a77
commit
953efdbfe7
@ -7,7 +7,7 @@ SET(CMAKE_INCLUDE_CURRENT_DIR ON)
|
||||
|
||||
IF(CMAKE_COMPILER_IS_GNUCXX)
|
||||
IF(CMAKE_SYSTEM_NAME MATCHES Linux)
|
||||
SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -pedantic -Wnon-virtual-dtor -Wno-long-long -ansi -Wundef -Wcast-align -Wchar-subscripts -Wall -W -Wpointer-arith -Wwrite-strings -Wformat-security -fno-exceptions -fno-check-new -fno-common -fstrict-aliasing")
|
||||
SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O1 -g -pedantic -Wnon-virtual-dtor -Wno-long-long -ansi -Wundef -Wcast-align -Wchar-subscripts -Wall -W -Wpointer-arith -Wwrite-strings -Wformat-security -fno-exceptions -fno-check-new -fno-common -fstrict-aliasing")
|
||||
IF(TEST_OPENMP)
|
||||
SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fopenmp")
|
||||
MESSAGE("Enabling OpenMP in tests/examples")
|
||||
|
||||
@ -35,11 +35,7 @@ namespace Eigen {
|
||||
#include "src/Core/CwiseBinaryOp.h"
|
||||
#include "src/Core/CwiseUnaryOp.h"
|
||||
#include "src/Core/CwiseNullaryOp.h"
|
||||
#if (defined EIGEN_WIP_PRODUCT)
|
||||
#include "src/Core/ProductWIP.h"
|
||||
#else
|
||||
#include "src/Core/Product.h"
|
||||
#endif
|
||||
#include "src/Core/Block.h"
|
||||
#include "src/Core/Minor.h"
|
||||
#include "src/Core/Transpose.h"
|
||||
@ -54,7 +50,8 @@ namespace Eigen {
|
||||
#include "src/Core/Swap.h"
|
||||
#include "src/Core/CommaInitializer.h"
|
||||
#include "src/Core/Triangular.h"
|
||||
#include "src/Core/TriangularAssign.h"
|
||||
#include "src/Core/Part.h"
|
||||
#include "src/Core/InverseProduct.h"
|
||||
|
||||
} // namespace Eigen
|
||||
|
||||
|
||||
@ -1,10 +1,8 @@
|
||||
#ifndef EIGEN_CORE_DECLARATIONS_H
|
||||
#define EIGEN_CORE_DECLARATIONS_H
|
||||
|
||||
#define EIGEN_WIP_PRODUCT
|
||||
|
||||
#ifdef __GNUC__
|
||||
#define EIGEN_GNUC_AT_LEAST(x,y) (__GNUC__>=x && __GNUC_MINOR__>=y) || __GNUC__>x
|
||||
#define EIGEN_GNUC_AT_LEAST(x,y) ((__GNUC__>=x && __GNUC_MINOR__>=y) || __GNUC__>x)
|
||||
#else
|
||||
#define EIGEN_GNUC_AT_LEAST(x,y) 0
|
||||
#endif
|
||||
|
||||
@ -58,9 +58,9 @@ template<typename MatrixType> class Cholesky
|
||||
compute(matrix);
|
||||
}
|
||||
|
||||
Triangular<Lower, MatrixType> matrixL(void) const
|
||||
Extract<MatrixType, Lower> matrixL(void) const
|
||||
{
|
||||
return m_matrix.lower();
|
||||
return m_matrix;
|
||||
}
|
||||
|
||||
bool isPositiveDefinite(void) const { return m_isPositiveDefinite; }
|
||||
@ -118,7 +118,7 @@ typename Derived::Eval Cholesky<MatrixType>::solve(MatrixBase<Derived> &b)
|
||||
const int size = m_matrix.rows();
|
||||
ei_assert(size==b.size());
|
||||
|
||||
return m_matrix.adjoint().upper().inverseProduct(m_matrix.lower().inverseProduct(b));
|
||||
return m_matrix.adjoint().template extract<Upper>().inverseProduct(matrixL().inverseProduct(b));
|
||||
}
|
||||
|
||||
|
||||
|
||||
@ -58,9 +58,9 @@ template<typename MatrixType> class CholeskyWithoutSquareRoot
|
||||
}
|
||||
|
||||
/** \returns the lower triangular matrix L */
|
||||
Triangular<Lower|UnitDiagBit, MatrixType > matrixL(void) const
|
||||
Extract<MatrixType, UnitLower> matrixL(void) const
|
||||
{
|
||||
return m_matrix.lowerWithUnitDiag();
|
||||
return m_matrix;
|
||||
}
|
||||
|
||||
/** \returns the coefficients of the diagonal matrix D */
|
||||
@ -131,9 +131,9 @@ typename Derived::Eval CholeskyWithoutSquareRoot<MatrixType>::solve(MatrixBase<D
|
||||
const int size = m_matrix.rows();
|
||||
ei_assert(size==vecB.size());
|
||||
|
||||
return m_matrix.adjoint().upperWithUnitDiag()
|
||||
return m_matrix.adjoint().template extract<UnitUpper>()
|
||||
.inverseProduct(
|
||||
(m_matrix.lowerWithUnitDiag()
|
||||
(matrixL()
|
||||
.inverseProduct(vecB))
|
||||
.cwiseQuotient(m_matrix.diagonal())
|
||||
);
|
||||
|
||||
@ -104,8 +104,7 @@ bool Vectorize = (int(Derived::Flags) & int(OtherDerived::Flags) & VectorizableB
|
||||
&& ( (int(Derived::Flags) & int(OtherDerived::Flags) & Like1DArrayBit)
|
||||
||((int(Derived::Flags)&RowMajorBit)
|
||||
? int(Derived::ColsAtCompileTime)!=Dynamic && (int(Derived::ColsAtCompileTime)%ei_packet_traits<typename Derived::Scalar>::size==0)
|
||||
: int(Derived::RowsAtCompileTime)!=Dynamic && (int(Derived::RowsAtCompileTime)%ei_packet_traits<typename Derived::Scalar>::size==0)) ),
|
||||
bool TriangularAssign = Derived::Flags & (NullLowerBit | NullUpperBit)>
|
||||
: int(Derived::RowsAtCompileTime)!=Dynamic && (int(Derived::RowsAtCompileTime)%ei_packet_traits<typename Derived::Scalar>::size==0)) )>
|
||||
struct ei_assignment_impl;
|
||||
|
||||
template<typename Derived>
|
||||
@ -146,7 +145,7 @@ inline Derived& MatrixBase<Derived>
|
||||
}
|
||||
|
||||
template <typename Derived, typename OtherDerived>
|
||||
struct ei_assignment_impl<Derived, OtherDerived, false, false>
|
||||
struct ei_assignment_impl<Derived, OtherDerived, false>
|
||||
{
|
||||
static void execute(Derived & dst, const OtherDerived & src)
|
||||
{
|
||||
@ -180,7 +179,7 @@ struct ei_assignment_impl<Derived, OtherDerived, false, false>
|
||||
};
|
||||
|
||||
template <typename Derived, typename OtherDerived>
|
||||
struct ei_assignment_impl<Derived, OtherDerived, true, false>
|
||||
struct ei_assignment_impl<Derived, OtherDerived, true>
|
||||
{
|
||||
static void execute(Derived & dst, const OtherDerived & src)
|
||||
{
|
||||
|
||||
@ -65,11 +65,11 @@ struct ei_traits<CwiseBinaryOp<BinaryOp, Lhs, Rhs> >
|
||||
ColsAtCompileTime = Lhs::ColsAtCompileTime,
|
||||
MaxRowsAtCompileTime = Lhs::MaxRowsAtCompileTime,
|
||||
MaxColsAtCompileTime = Lhs::MaxColsAtCompileTime,
|
||||
Flags = ((int(LhsFlags) | int(RhsFlags)) & (
|
||||
Flags = (int(LhsFlags) | int(RhsFlags)) & (
|
||||
HereditaryBits
|
||||
| int(LhsFlags) & int(RhsFlags) & Like1DArrayBit
|
||||
| (int(LhsFlags) & int(RhsFlags) & Like1DArrayBit)
|
||||
| (ei_functor_traits<BinaryOp>::IsVectorizable && ((int(LhsFlags) & RowMajorBit)==(int(RhsFlags) & RowMajorBit))
|
||||
? int(LhsFlags) & int(RhsFlags) & VectorizableBit : 0))),
|
||||
? int(LhsFlags) & int(RhsFlags) & VectorizableBit : 0)),
|
||||
CoeffReadCost = LhsCoeffReadCost + RhsCoeffReadCost + ei_functor_traits<BinaryOp>::Cost
|
||||
};
|
||||
};
|
||||
|
||||
@ -60,50 +60,70 @@ template<typename ExpressionType, unsigned int Added, unsigned int Removed> clas
|
||||
|
||||
EIGEN_GENERIC_PUBLIC_INTERFACE(Flagged)
|
||||
|
||||
inline Flagged(const ExpressionType& matrix) : m_expression(matrix) {}
|
||||
inline Flagged(const ExpressionType& matrix) : m_matrix(matrix) {}
|
||||
|
||||
/** \internal */
|
||||
inline ExpressionType _expression() const { return m_expression; }
|
||||
inline ExpressionType _expression() const { return m_matrix; }
|
||||
|
||||
private:
|
||||
|
||||
inline int _rows() const { return m_expression.rows(); }
|
||||
inline int _cols() const { return m_expression.cols(); }
|
||||
inline int _rows() const { return m_matrix.rows(); }
|
||||
inline int _cols() const { return m_matrix.cols(); }
|
||||
inline int _stride() const { return m_matrix.stride(); }
|
||||
|
||||
inline const Scalar _coeff(int row, int col) const
|
||||
{
|
||||
return m_expression.coeff(row, col);
|
||||
return m_matrix.coeff(row, col);
|
||||
}
|
||||
|
||||
inline Scalar& _coeffRef(int row, int col)
|
||||
{
|
||||
return m_expression.const_cast_derived().coeffRef(row, col);
|
||||
return m_matrix.const_cast_derived().coeffRef(row, col);
|
||||
}
|
||||
|
||||
template<int LoadMode>
|
||||
inline const PacketScalar _packetCoeff(int row, int col) const
|
||||
{
|
||||
return m_expression.template packetCoeff<LoadMode>(row, col);
|
||||
return m_matrix.template packetCoeff<LoadMode>(row, col);
|
||||
}
|
||||
|
||||
template<int LoadMode>
|
||||
inline void _writePacketCoeff(int row, int col, const PacketScalar& x)
|
||||
{
|
||||
m_expression.const_cast_derived().template writePacketCoeff<LoadMode>(row, col, x);
|
||||
m_matrix.const_cast_derived().template writePacketCoeff<LoadMode>(row, col, x);
|
||||
}
|
||||
|
||||
protected:
|
||||
const ExpressionType m_expression;
|
||||
const ExpressionType m_matrix;
|
||||
};
|
||||
|
||||
/** \returns an expression of the temporary version of *this.
|
||||
/** \returns an expression of *this with added flags
|
||||
*/
|
||||
template<typename Derived>
|
||||
template<unsigned int Added, unsigned int Removed>
|
||||
inline const Flagged<Derived, Added, Removed>
|
||||
MatrixBase<Derived>::flagged() const
|
||||
template<unsigned int Added>
|
||||
inline const Flagged<Derived, Added, 0>
|
||||
MatrixBase<Derived>::marked() const
|
||||
{
|
||||
return Flagged<Derived, Added, Removed>(derived());
|
||||
return derived();
|
||||
}
|
||||
|
||||
/** \returns an expression of *this with the following flags removed:
|
||||
* EvalBeforeNestingBit and EvalBeforeAssigningBit.
|
||||
*/
|
||||
template<typename Derived>
|
||||
inline const Flagged<Derived, 0, EvalBeforeNestingBit | EvalBeforeAssigningBit>
|
||||
MatrixBase<Derived>::lazy() const
|
||||
{
|
||||
return derived();
|
||||
}
|
||||
|
||||
/** \returns an expression of *this with the NestByValueBit flag added.
|
||||
*/
|
||||
template<typename Derived>
|
||||
inline const Flagged<Derived, NestByValueBit, 0>
|
||||
MatrixBase<Derived>::temporary() const
|
||||
{
|
||||
return derived();
|
||||
}
|
||||
|
||||
#endif // EIGEN_FLAGGED_H
|
||||
|
||||
83
Eigen/src/Core/InverseProduct.h
Executable file
83
Eigen/src/Core/InverseProduct.h
Executable file
@ -0,0 +1,83 @@
|
||||
// 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/>.
|
||||
|
||||
#ifndef EIGEN_INVERSEPRODUCT_H
|
||||
#define EIGEN_INVERSEPRODUCT_H
|
||||
|
||||
/** \returns the product of the inverse of *this with \a other.
|
||||
*
|
||||
* This function computes the inverse-matrix matrix product inverse(*this) * \a other
|
||||
* It works as a forward (resp. backward) substitution if *this is an upper (resp. lower)
|
||||
* triangular matrix.
|
||||
*/
|
||||
template<typename Derived>
|
||||
template<typename OtherDerived>
|
||||
typename OtherDerived::Eval MatrixBase<Derived>::inverseProduct(const MatrixBase<OtherDerived>& other) const
|
||||
{
|
||||
assert(cols() == other.rows());
|
||||
assert(!(Flags & ZeroDiagBit));
|
||||
assert(Flags & (UpperTriangularBit|LowerTriangularBit));
|
||||
|
||||
typename OtherDerived::Eval res(other.rows(), other.cols());
|
||||
|
||||
for(int c=0 ; c<other.cols() ; ++c)
|
||||
{
|
||||
if(Flags & LowerTriangularBit)
|
||||
{
|
||||
// forward substitution
|
||||
if(Flags & UnitDiagBit)
|
||||
res.coeffRef(0,c) = other.coeff(0,c);
|
||||
else
|
||||
res.coeffRef(0,c) = other.coeff(0,c)/coeff(0, 0);
|
||||
for(int i=1; i<rows(); ++i)
|
||||
{
|
||||
Scalar tmp = other.coeff(i,c) - ((this->row(i).start(i)) * res.col(c).start(i)).coeff(0,0);
|
||||
if (Flags & UnitDiagBit)
|
||||
res.coeffRef(i,c) = tmp;
|
||||
else
|
||||
res.coeffRef(i,c) = tmp/coeff(i,i);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
// backward substitution
|
||||
if(Flags & UnitDiagBit)
|
||||
res.coeffRef(cols()-1,c) = other.coeff(cols()-1,c);
|
||||
else
|
||||
res.coeffRef(cols()-1,c) = other.coeff(cols()-1, c)/coeff(rows()-1, cols()-1);
|
||||
for(int i=rows()-2 ; i>=0 ; --i)
|
||||
{
|
||||
Scalar tmp = other.coeff(i,c)
|
||||
- ((this->row(i).end(cols()-i-1)) * res.col(c).end(cols()-i-1)).coeff(0,0);
|
||||
if (Flags & UnitDiagBit)
|
||||
res.coeffRef(i,c) = tmp;
|
||||
else
|
||||
res.coeffRef(i,c) = tmp/coeff(i,i);
|
||||
}
|
||||
}
|
||||
}
|
||||
return res;
|
||||
}
|
||||
|
||||
#endif // EIGEN_INVERSEPRODUCT_H
|
||||
@ -136,6 +136,12 @@ template<typename Derived> class MatrixBase
|
||||
CoeffReadCost = ei_traits<Derived>::CoeffReadCost
|
||||
};
|
||||
|
||||
MatrixBase()
|
||||
{
|
||||
assert(!( (Flags&UnitDiagBit && Flags&ZeroDiagBit)
|
||||
|| (Flags&UpperTriangularBit && Flags&LowerTriangularBit) ));
|
||||
}
|
||||
|
||||
/** This is the "real scalar" type; if the \a Scalar type is already real numbers
|
||||
* (e.g. int, float or double) then \a RealScalar is just the same as \a Scalar. If
|
||||
* \a Scalar is \a std::complex<T> then RealScalar is \a T.
|
||||
@ -280,6 +286,10 @@ template<typename Derived> class MatrixBase
|
||||
|
||||
template<typename OtherDerived>
|
||||
Derived& operator*=(const MatrixBase<OtherDerived>& other);
|
||||
|
||||
template<typename OtherDerived>
|
||||
typename OtherDerived::Eval inverseProduct(const MatrixBase<OtherDerived>& other) const;
|
||||
|
||||
//@}
|
||||
|
||||
/** \name Dot product and related notions
|
||||
@ -296,7 +306,7 @@ template<typename Derived> class MatrixBase
|
||||
const Transpose<Derived> transpose() const;
|
||||
const Transpose<
|
||||
Flagged<CwiseUnaryOp<ei_scalar_conjugate_op<typename ei_traits<Derived>::Scalar>, Derived>
|
||||
, TemporaryBit, 0> >
|
||||
, NestByValueBit, 0> >
|
||||
adjoint() const;
|
||||
//@}
|
||||
|
||||
@ -347,6 +357,9 @@ template<typename Derived> class MatrixBase
|
||||
|
||||
DiagonalCoeffs<Derived> diagonal();
|
||||
const DiagonalCoeffs<Derived> diagonal() const;
|
||||
|
||||
template<unsigned int PartType> Part<Derived, PartType> part();
|
||||
template<unsigned int ExtractType> const Extract<Derived, ExtractType> extract() const;
|
||||
//@}
|
||||
|
||||
/// \name Generating special matrices
|
||||
@ -406,6 +419,9 @@ template<typename Derived> class MatrixBase
|
||||
bool isIdentity(RealScalar prec = precision<Scalar>()) const;
|
||||
bool isDiagonal(RealScalar prec = precision<Scalar>()) const;
|
||||
|
||||
bool isUpper(RealScalar prec = precision<Scalar>()) const;
|
||||
bool isLower(RealScalar prec = precision<Scalar>()) const;
|
||||
|
||||
template<typename OtherDerived>
|
||||
bool isOrtho(const MatrixBase<OtherDerived>& other,
|
||||
RealScalar prec = precision<Scalar>()) const;
|
||||
@ -433,24 +449,17 @@ template<typename Derived> class MatrixBase
|
||||
template<typename OtherDerived>
|
||||
void swap(const MatrixBase<OtherDerived>& other);
|
||||
|
||||
template<unsigned int Added, unsigned int Removed>
|
||||
const Flagged<Derived, Added, Removed> flagged() const;
|
||||
|
||||
const Flagged<Derived, 0, EvalBeforeNestingBit | EvalBeforeAssigningBit> lazy() const
|
||||
{
|
||||
return derived();
|
||||
}
|
||||
const Flagged<Derived, TemporaryBit, 0> temporary() const
|
||||
{
|
||||
return derived();
|
||||
}
|
||||
template<unsigned int Added>
|
||||
const Flagged<Derived, Added, 0> marked() const;
|
||||
const Flagged<Derived, 0, EvalBeforeNestingBit | EvalBeforeAssigningBit> lazy() const;
|
||||
const Flagged<Derived, NestByValueBit, 0> temporary() const;
|
||||
|
||||
/** \returns number of elements to skip to pass from one row (resp. column) to another
|
||||
* for a row-major (resp. column-major) matrix.
|
||||
* Combined with coeffRef() and the compile times flags, it allows a direct access to the data
|
||||
* of the underlying matrix.
|
||||
*/
|
||||
int stride(void) const { return derived()._stride(); }
|
||||
inline int stride(void) const { return derived()._stride(); }
|
||||
//@}
|
||||
|
||||
/// \name Coefficient-wise operations
|
||||
@ -553,22 +562,6 @@ template<typename Derived> class MatrixBase
|
||||
{ return *static_cast<Derived*>(const_cast<MatrixBase*>(this)); }
|
||||
//@}
|
||||
|
||||
/// \name Triangular matrices
|
||||
//@{
|
||||
Triangular<Upper, Derived> upper(void);
|
||||
const Triangular<Upper, Derived> upper(void) const;
|
||||
const Triangular<Upper|UnitDiagBit, Derived> upperWithUnitDiag(void) const;
|
||||
const Triangular<Upper|NullDiagBit, Derived> upperWithNullDiag(void) const;
|
||||
|
||||
Triangular<Lower, Derived> lower(void);
|
||||
const Triangular<Lower, Derived> lower(void) const;
|
||||
const Triangular<Lower|UnitDiagBit, Derived> lowerWithUnitDiag(void) const;
|
||||
const Triangular<Lower|NullDiagBit, Derived> lowerWithNullDiag(void) const;
|
||||
|
||||
bool isUpper(RealScalar prec = precision<Scalar>()) const;
|
||||
bool isLower(RealScalar prec = precision<Scalar>()) const;
|
||||
//@}
|
||||
|
||||
/** \name LU module
|
||||
*
|
||||
* \code #include <Eigen/LU> \endcode
|
||||
|
||||
220
Eigen/src/Core/Part.h
Normal file
220
Eigen/src/Core/Part.h
Normal file
@ -0,0 +1,220 @@
|
||||
// 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 Benoit Jacob <jacob@math.jussieu.fr>
|
||||
// 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/>.
|
||||
|
||||
#ifndef EIGEN_PART_H
|
||||
#define EIGEN_PART_H
|
||||
|
||||
template<typename MatrixType, unsigned int Mode>
|
||||
class Part
|
||||
{
|
||||
public:
|
||||
Part(MatrixType& matrix) : m_matrix(matrix) {}
|
||||
template<typename Other> void lazyAssign(const Other& other);
|
||||
template<typename Other> void operator=(const Other& other);
|
||||
template<typename Other> void operator+=(const Other& other);
|
||||
template<typename Other> void operator-=(const Other& other);
|
||||
void operator*=(const typename ei_traits<MatrixType>::Scalar& other);
|
||||
void operator/=(const typename ei_traits<MatrixType>::Scalar& other);
|
||||
void setConstant(const typename ei_traits<MatrixType>::Scalar& value);
|
||||
void setZero();
|
||||
void setOnes();
|
||||
void setRandom();
|
||||
void setIdentity();
|
||||
|
||||
private:
|
||||
MatrixType& m_matrix;
|
||||
};
|
||||
|
||||
template<typename MatrixType, unsigned int Mode>
|
||||
template<typename Other>
|
||||
inline void Part<MatrixType, Mode>::operator=(const Other& other)
|
||||
{
|
||||
if(Other::Flags & EvalBeforeAssigningBit)
|
||||
lazyAssign(other.eval());
|
||||
else
|
||||
lazyAssign(other.derived());
|
||||
}
|
||||
|
||||
template<typename Derived1, typename Derived2, unsigned int Mode, int UnrollCount>
|
||||
struct ei_part_assignment_unroller
|
||||
{
|
||||
enum {
|
||||
col = (UnrollCount-1) / Derived1::RowsAtCompileTime,
|
||||
row = (UnrollCount-1) % Derived1::RowsAtCompileTime
|
||||
};
|
||||
|
||||
inline static void run(Derived1 &dst, const Derived2 &src)
|
||||
{
|
||||
ei_part_assignment_unroller<Derived1, Derived2, Mode, UnrollCount-1>::run(dst, src);
|
||||
|
||||
if(Mode == SelfAdjoint)
|
||||
{
|
||||
if(row == col)
|
||||
dst.coeffRef(row, col) = ei_real(src.coeff(row, col));
|
||||
else if(row < col)
|
||||
dst.coeffRef(col, row) = ei_conj(dst.coeffRef(row, col) = src.coeff(row, col));
|
||||
}
|
||||
else
|
||||
{
|
||||
if((Mode == Upper && row <= col)
|
||||
|| (Mode == Lower && row >= col)
|
||||
|| (Mode == StrictlyUpper && row < col)
|
||||
|| (Mode == StrictlyLower && row > col))
|
||||
dst.coeffRef(row, col) = src.coeff(row, col);
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
template<typename Derived1, typename Derived2, unsigned int Mode>
|
||||
struct ei_part_assignment_unroller<Derived1, Derived2, Mode, 1>
|
||||
{
|
||||
inline static void run(Derived1 &dst, const Derived2 &src)
|
||||
{
|
||||
if(!(Mode & ZeroDiagBit))
|
||||
dst.coeffRef(0, 0) = src.coeff(0, 0);
|
||||
}
|
||||
};
|
||||
|
||||
// prevent buggy user code from causing an infinite recursion
|
||||
template<typename Derived1, typename Derived2, unsigned int Mode>
|
||||
struct ei_part_assignment_unroller<Derived1, Derived2, Mode, 0>
|
||||
{
|
||||
inline static void run(Derived1 &, const Derived2 &) {}
|
||||
};
|
||||
|
||||
template<typename Derived1, typename Derived2, unsigned int Mode>
|
||||
struct ei_part_assignment_unroller<Derived1, Derived2, Mode, Dynamic>
|
||||
{
|
||||
inline static void run(Derived1 &, const Derived2 &) {}
|
||||
};
|
||||
|
||||
|
||||
template<typename MatrixType, unsigned int Mode>
|
||||
template<typename Other>
|
||||
void Part<MatrixType, Mode>::lazyAssign(const Other& other)
|
||||
{
|
||||
const bool unroll = MatrixType::SizeAtCompileTime * Other::CoeffReadCost / 2 <= EIGEN_UNROLLING_LIMIT;
|
||||
ei_assert(m_matrix.rows() == other.rows() && m_matrix.cols() == other.cols());
|
||||
if(unroll)
|
||||
{
|
||||
ei_part_assignment_unroller
|
||||
<MatrixType, Other, Mode,
|
||||
unroll ? int(MatrixType::SizeAtCompileTime) : Dynamic
|
||||
>::run(m_matrix, other.derived());
|
||||
}
|
||||
else
|
||||
{
|
||||
switch(Mode)
|
||||
{
|
||||
case Upper:
|
||||
for(int j = 0; j < m_matrix.cols(); j++)
|
||||
for(int i = 0; i <= j; i++)
|
||||
m_matrix.coeffRef(i, j) = other.coeff(i, j);
|
||||
break;
|
||||
case Lower:
|
||||
for(int j = 0; j < m_matrix.cols(); j++)
|
||||
for(int i = j; i < m_matrix.rows(); i++)
|
||||
m_matrix.coeffRef(i, j) = other.coeff(i, j);
|
||||
break;
|
||||
case StrictlyUpper:
|
||||
for(int j = 0; j < m_matrix.cols(); j++)
|
||||
for(int i = 0; i < j; i++)
|
||||
m_matrix.coeffRef(i, j) = other.coeff(i, j);
|
||||
break;
|
||||
case StrictlyLower:
|
||||
for(int j = 0; j < m_matrix.cols(); j++)
|
||||
for(int i = j+1; i < m_matrix.rows(); i++)
|
||||
m_matrix.coeffRef(i, j) = other.coeff(i, j);
|
||||
break;
|
||||
case SelfAdjoint:
|
||||
for(int j = 0; j < m_matrix.cols(); j++)
|
||||
{
|
||||
for(int i = 0; i < j; i++)
|
||||
m_matrix.coeffRef(j, i) = ei_conj(m_matrix.coeffRef(i, j) = other.coeff(i, j));
|
||||
m_matrix.coeffRef(j, j) = ei_real(other.coeff(j, j));
|
||||
}
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template<typename MatrixType, unsigned int Mode>
|
||||
template<typename Other> inline void Part<MatrixType, Mode>::operator+=(const Other& other)
|
||||
{
|
||||
*this = m_matrix + other;
|
||||
}
|
||||
|
||||
template<typename MatrixType, unsigned int Mode>
|
||||
template<typename Other> inline void Part<MatrixType, Mode>::operator-=(const Other& other)
|
||||
{
|
||||
*this = m_matrix - other;
|
||||
}
|
||||
|
||||
template<typename MatrixType, unsigned int Mode>
|
||||
inline void Part<MatrixType, Mode>::operator*=
|
||||
(const typename ei_traits<MatrixType>::Scalar& other)
|
||||
{
|
||||
*this = m_matrix * other;
|
||||
}
|
||||
|
||||
template<typename MatrixType, unsigned int Mode>
|
||||
inline void Part<MatrixType, Mode>::operator/=
|
||||
(const typename ei_traits<MatrixType>::Scalar& other)
|
||||
{
|
||||
*this = m_matrix / other;
|
||||
}
|
||||
|
||||
template<typename MatrixType, unsigned int Mode>
|
||||
inline void Part<MatrixType, Mode>::setConstant(const typename ei_traits<MatrixType>::Scalar& value)
|
||||
{
|
||||
*this = MatrixType::constant(m_matrix.rows(), m_matrix.cols(), value);
|
||||
}
|
||||
|
||||
template<typename MatrixType, unsigned int Mode>
|
||||
inline void Part<MatrixType, Mode>::setZero()
|
||||
{
|
||||
setConstant((typename ei_traits<MatrixType>::Scalar)(0));
|
||||
}
|
||||
|
||||
template<typename MatrixType, unsigned int Mode>
|
||||
inline void Part<MatrixType, Mode>::setOnes()
|
||||
{
|
||||
setConstant((typename ei_traits<MatrixType>::Scalar)(1));
|
||||
}
|
||||
|
||||
template<typename MatrixType, unsigned int Mode>
|
||||
inline void Part<MatrixType, Mode>::setRandom()
|
||||
{
|
||||
*this = MatrixType::random(m_matrix.rows(), m_matrix.cols());
|
||||
}
|
||||
|
||||
template<typename Derived>
|
||||
template<unsigned int Mode>
|
||||
inline Part<Derived, Mode> MatrixBase<Derived>::part()
|
||||
{
|
||||
return Part<Derived, Mode>(derived());
|
||||
}
|
||||
|
||||
#endif // EIGEN_PART_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_is_temporary<T>::ret && !(ei_traits<T>::Flags & RowMajorBit)),
|
||||
(ei_traits<T>::Flags & NestByValueBit) && !(ei_traits<T>::Flags & RowMajorBit),
|
||||
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_is_temporary<T>::ret,
|
||||
ei_traits<T>::Flags & NestByValueBit,
|
||||
T,
|
||||
typename ei_meta_if<
|
||||
int(ei_traits<T>::Flags) & EvalBeforeNestingBit
|
||||
@ -202,7 +202,7 @@ struct ei_traits<Product<Lhs, Rhs, EvalMode> >
|
||||
// the cache friendly product evals lhs once only
|
||||
// FIXME what to do if we chose to dynamically call the normal product from the cache friendly one for small matrices ?
|
||||
typedef typename ei_meta_if<EvalMode==CacheFriendlyProduct,
|
||||
typename ei_product_nested_lhs<Lhs,0>::type,
|
||||
typename ei_product_nested_lhs<Lhs,1>::type,
|
||||
typename ei_nested<Lhs,Rhs::ColsAtCompileTime>::type>::ret LhsNested;
|
||||
|
||||
// NOTE that rhs must be ColumnMajor, so we might need a special nested type calculation
|
||||
@ -231,9 +231,7 @@ struct ei_traits<Product<Lhs, Rhs, EvalMode> >
|
||||
(_RowMajor ? 0 : RowMajorBit)
|
||||
| ((RowsAtCompileTime == Dynamic || ColsAtCompileTime == Dynamic) ? 0 : LargeBit)),
|
||||
Flags = ((unsigned int)(LhsFlags | RhsFlags) & _LostBits)
|
||||
#ifndef EIGEN_WIP_PRODUCT_DIRTY
|
||||
| EvalBeforeAssigningBit //FIXME
|
||||
#endif
|
||||
| EvalBeforeAssigningBit
|
||||
| EvalBeforeNestingBit
|
||||
| (_Vectorizable ? VectorizableBit : 0),
|
||||
CoeffReadCost
|
||||
@ -335,15 +333,15 @@ template<typename Lhs, typename Rhs, int EvalMode> class Product : ei_no_assignm
|
||||
return res;
|
||||
}
|
||||
|
||||
template<typename Lhs_, typename Rhs_, int EvalMode_, typename DestDerived_, bool DirectAccess_>
|
||||
friend struct ei_cache_friendly_selector;
|
||||
|
||||
protected:
|
||||
const LhsNested m_lhs;
|
||||
const RhsNested m_rhs;
|
||||
};
|
||||
|
||||
/** \returns the matrix product of \c *this and \a other.
|
||||
*
|
||||
* \note This function causes an immediate evaluation. If you want to perform a matrix product
|
||||
* without immediate evaluation, call .lazy() on one of the matrices before taking the product.
|
||||
*
|
||||
* \sa lazy(), operator*=(const MatrixBase&)
|
||||
*/
|
||||
@ -385,54 +383,89 @@ inline Derived& MatrixBase<Derived>::lazyAssign(const Product<Lhs,Rhs,CacheFrien
|
||||
return derived();
|
||||
}
|
||||
|
||||
template<typename Lhs, typename Rhs, int EvalMode>
|
||||
template<typename DestDerived>
|
||||
inline void Product<Lhs,Rhs,EvalMode>::_cacheFriendlyEval(DestDerived& res) const
|
||||
template<typename Lhs, typename Rhs, int EvalMode, typename DestDerived, bool DirectAccess>
|
||||
struct ei_cache_friendly_selector
|
||||
{
|
||||
if ( _rows()>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD
|
||||
&& _cols()>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD
|
||||
&& m_lhs.cols()>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD
|
||||
typedef Product<Lhs,Rhs,EvalMode> Prod;
|
||||
typedef typename Prod::_LhsNested _LhsNested;
|
||||
typedef typename Prod::_RhsNested _RhsNested;
|
||||
typedef typename Prod::Scalar Scalar;
|
||||
static inline void eval(const Prod& product, DestDerived& res)
|
||||
{
|
||||
if ( product._rows()>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD
|
||||
&& product._cols()>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD
|
||||
&& product.m_lhs.cols()>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD
|
||||
)
|
||||
{
|
||||
#ifndef EIGEN_WIP_PRODUCT_DIRTY
|
||||
res.setZero();
|
||||
#endif
|
||||
|
||||
ei_cache_friendly_product<Scalar>(
|
||||
_rows(), _cols(), m_lhs.cols(),
|
||||
_LhsNested::Flags&RowMajorBit, &(m_lhs.const_cast_derived().coeffRef(0,0)), m_lhs.stride(),
|
||||
_RhsNested::Flags&RowMajorBit, &(m_rhs.const_cast_derived().coeffRef(0,0)), m_rhs.stride(),
|
||||
Flags&RowMajorBit, &(res.coeffRef(0,0)), res.stride()
|
||||
product._rows(), product._cols(), product.m_lhs.cols(),
|
||||
_LhsNested::Flags&RowMajorBit, &(product.m_lhs.const_cast_derived().coeffRef(0,0)), product.m_lhs.stride(),
|
||||
_RhsNested::Flags&RowMajorBit, &(product.m_rhs.const_cast_derived().coeffRef(0,0)), product.m_rhs.stride(),
|
||||
Prod::Flags&RowMajorBit, &(res.coeffRef(0,0)), res.stride()
|
||||
);
|
||||
}
|
||||
else
|
||||
{
|
||||
// lazy product
|
||||
res = Product<_LhsNested,_RhsNested,NormalProduct>(m_lhs, m_rhs).lazy();
|
||||
res = Product<_LhsNested,_RhsNested,NormalProduct>(product.m_lhs, product.m_rhs).lazy();
|
||||
}
|
||||
}
|
||||
|
||||
static inline void eval_and_add(const Prod& product, DestDerived& res)
|
||||
{
|
||||
if ( product._rows()>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD
|
||||
&& product._cols()>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD
|
||||
&& product.m_lhs.cols()>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD
|
||||
)
|
||||
{
|
||||
ei_cache_friendly_product<Scalar>(
|
||||
product._rows(), product._cols(), product.m_lhs.cols(),
|
||||
_LhsNested::Flags&RowMajorBit, &(product.m_lhs.const_cast_derived().coeffRef(0,0)), product.m_lhs.stride(),
|
||||
_RhsNested::Flags&RowMajorBit, &(product.m_rhs.const_cast_derived().coeffRef(0,0)), product.m_rhs.stride(),
|
||||
Prod::Flags&RowMajorBit, &(res.coeffRef(0,0)), res.stride()
|
||||
);
|
||||
}
|
||||
else
|
||||
{
|
||||
res += Product<_LhsNested,_RhsNested,NormalProduct>(product.m_lhs, product.m_rhs).lazy();
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
template<typename Lhs, typename Rhs, int EvalMode, typename DestDerived>
|
||||
struct ei_cache_friendly_selector<Lhs,Rhs,EvalMode,DestDerived,false>
|
||||
{
|
||||
typedef Product<Lhs,Rhs,EvalMode> Prod;
|
||||
typedef typename Prod::_LhsNested _LhsNested;
|
||||
typedef typename Prod::_RhsNested _RhsNested;
|
||||
typedef typename Prod::Scalar Scalar;
|
||||
static inline void eval(const Prod& product, DestDerived& res)
|
||||
{
|
||||
res = Product<_LhsNested,_RhsNested,NormalProduct>(product.m_lhs, product.m_rhs).lazy();
|
||||
}
|
||||
|
||||
static inline void eval_and_add(const Prod& product, DestDerived& res)
|
||||
{
|
||||
res += Product<_LhsNested,_RhsNested,NormalProduct>(product.m_lhs, product.m_rhs).lazy();
|
||||
}
|
||||
};
|
||||
|
||||
template<typename Lhs, typename Rhs, int EvalMode>
|
||||
template<typename DestDerived>
|
||||
inline void Product<Lhs,Rhs,EvalMode>::_cacheFriendlyEval(DestDerived& res) const
|
||||
{
|
||||
ei_cache_friendly_selector<Lhs,Rhs,EvalMode,DestDerived,
|
||||
_LhsNested::Flags&_RhsNested::Flags&DirectAccessBit>
|
||||
::eval(*this, res);
|
||||
}
|
||||
|
||||
template<typename Lhs, typename Rhs, int EvalMode>
|
||||
template<typename DestDerived>
|
||||
inline void Product<Lhs,Rhs,EvalMode>::_cacheFriendlyEvalAndAdd(DestDerived& res) const
|
||||
{
|
||||
if ( _rows()>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD
|
||||
&& _cols()>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD
|
||||
&& m_lhs.cols()>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD
|
||||
)
|
||||
{
|
||||
ei_cache_friendly_product<Scalar>(
|
||||
_rows(), _cols(), m_lhs.cols(),
|
||||
_LhsNested::Flags&RowMajorBit, &(m_lhs.const_cast_derived().coeffRef(0,0)), m_lhs.stride(),
|
||||
_RhsNested::Flags&RowMajorBit, &(m_rhs.const_cast_derived().coeffRef(0,0)), m_rhs.stride(),
|
||||
Flags&RowMajorBit, &(res.coeffRef(0,0)), res.stride()
|
||||
);
|
||||
}
|
||||
else
|
||||
{
|
||||
// lazy product
|
||||
res += Product<_LhsNested,_RhsNested,NormalProduct>(m_lhs, m_rhs).lazy();
|
||||
}
|
||||
ei_cache_friendly_selector<Lhs,Rhs,EvalMode,DestDerived,
|
||||
_LhsNested::Flags&_RhsNested::Flags&DirectAccessBit>
|
||||
::eval_and_add(*this, res);
|
||||
}
|
||||
|
||||
#endif // EIGEN_PRODUCT_H
|
||||
|
||||
@ -127,7 +127,7 @@ MatrixBase<Derived>::transpose() const
|
||||
template<typename Derived>
|
||||
inline const Transpose<
|
||||
Flagged<CwiseUnaryOp<ei_scalar_conjugate_op<typename ei_traits<Derived>::Scalar>, Derived >
|
||||
, TemporaryBit, 0> >
|
||||
, NestByValueBit, 0> >
|
||||
MatrixBase<Derived>::adjoint() const
|
||||
{
|
||||
return conjugate().temporary().transpose();
|
||||
|
||||
@ -2,6 +2,7 @@
|
||||
// for linear algebra. Eigen itself is part of the KDE project.
|
||||
//
|
||||
// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
|
||||
// Copyright (C) 2008 Benoit Jacob <jacob@math.jussieu.fr>
|
||||
//
|
||||
// Eigen is free software; you can redistribute it and/or
|
||||
// modify it under the terms of the GNU Lesser General Public
|
||||
@ -22,42 +23,27 @@
|
||||
// License and a copy of the GNU General Public License along with
|
||||
// Eigen. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
#ifndef EIGEN_TRIANGULAR_H
|
||||
#define EIGEN_TRIANGULAR_H
|
||||
#ifndef EIGEN_EXTRACT_H
|
||||
#define EIGEN_EXTRACT_H
|
||||
|
||||
/** \class Triangular
|
||||
/** \class Extract
|
||||
*
|
||||
* \brief Expression of a triangular matrix from a square matrix
|
||||
* \brief Expression of a triangular matrix extracted from a given matrix
|
||||
*
|
||||
* \param Mode or-ed bit field indicating the triangular part (Upper or Lower) we are taking,
|
||||
* and the property of the diagonal if any (UnitDiagBit or NullDiagBit).
|
||||
* \param MatrixType the type of the object in which we are taking the triangular part
|
||||
* \param Mode the kind of triangular matrix expression to construct. Can be Upper, StrictlyUpper,
|
||||
* UnitUpper, Lower, StrictlyLower, UnitLower. This is in fact a bit field; it must have either
|
||||
* UpperTriangularBit or LowerTriangularBit, and additionnaly it may have either ZeroDiagBit or
|
||||
* UnitDiagBit.
|
||||
*
|
||||
* This class represents an expression of the upper or lower triangular part of
|
||||
* a squared matrix. It is the return type of MatrixBase::upper(), MatrixBase::lower(),
|
||||
* MatrixBase::upperWithUnitDiagBit(), etc., and used to optimize operations involving
|
||||
* triangular matrices. Most of the time this is the only way it is used.
|
||||
* a square matrix, possibly with a further assumption on the diagonal. It is the return type
|
||||
* of MatrixBase::extract() and most of the time this is the only way it is used.
|
||||
*
|
||||
* Examples of some key features:
|
||||
* \code
|
||||
* m1 = (<any expression>).upper();
|
||||
* \endcode
|
||||
* In this example, the strictly lower part of the expression is not evaluated,
|
||||
* m1 might be resized and the strict lower part of m1 == 0.
|
||||
*
|
||||
* \code
|
||||
* m1.upper() = <any expression>;
|
||||
* \endcode
|
||||
* This example diverge from the previous one in the sense that the strictly
|
||||
* lower part of m1 is left unchanged, and optimal loops are employed. Note that
|
||||
* m1 might also be resized.
|
||||
*
|
||||
* Of course, in both examples \c <any \c expression> has to be a square matrix.
|
||||
*
|
||||
* \sa MatrixBase::upper(), MatrixBase::lower(), class TriangularProduct
|
||||
* \sa MatrixBase::extract()
|
||||
*/
|
||||
template<int Mode, typename MatrixType>
|
||||
struct ei_traits<Triangular<Mode, MatrixType> >
|
||||
template<typename MatrixType, unsigned int Mode>
|
||||
struct ei_traits<Extract<MatrixType, Mode> >
|
||||
{
|
||||
typedef typename MatrixType::Scalar Scalar;
|
||||
typedef typename ei_nested<MatrixType>::type MatrixTypeNested;
|
||||
@ -72,105 +58,30 @@ struct ei_traits<Triangular<Mode, MatrixType> >
|
||||
};
|
||||
};
|
||||
|
||||
template<int Mode, typename MatrixType> class Triangular
|
||||
: public MatrixBase<Triangular<Mode,MatrixType> >
|
||||
template<typename MatrixType, unsigned int Mode> class Extract
|
||||
: public MatrixBase<Extract<MatrixType, Mode> >
|
||||
{
|
||||
public:
|
||||
|
||||
EIGEN_GENERIC_PUBLIC_INTERFACE(Triangular)
|
||||
EIGEN_GENERIC_PUBLIC_INTERFACE(Extract)
|
||||
|
||||
inline Triangular(const MatrixType& matrix)
|
||||
: m_matrix(matrix)
|
||||
{
|
||||
assert(!( (Flags&UnitDiagBit) && (Flags&NullDiagBit)));
|
||||
assert(matrix.rows()==matrix.cols());
|
||||
}
|
||||
inline Extract(const MatrixType& matrix) : m_matrix(matrix) {}
|
||||
|
||||
EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Triangular)
|
||||
|
||||
/** Overloaded to keep a Triangular expression */
|
||||
inline Triangular<(Upper | Lower) ^ Mode, Flagged<Transpose<MatrixType>,TemporaryBit,0> > transpose()
|
||||
{
|
||||
return m_matrix.transpose().temporary();
|
||||
}
|
||||
|
||||
/** Overloaded to keep a Triangular expression */
|
||||
inline const Triangular<(Upper | Lower) ^ Mode, Flagged<Transpose<MatrixType>,TemporaryBit,0> > transpose() const
|
||||
{
|
||||
return m_matrix.transpose().temporary();
|
||||
}
|
||||
|
||||
/** \returns the product of the inverse of *this with \a other.
|
||||
*
|
||||
* This function computes the inverse-matrix matrix product inverse(*this) * \a other
|
||||
* It works as a forward (resp. backward) substitution if *this is an upper (resp. lower)
|
||||
* triangular matrix.
|
||||
*/
|
||||
template<typename OtherDerived>
|
||||
typename OtherDerived::Eval inverseProduct(const MatrixBase<OtherDerived>& other) const
|
||||
{
|
||||
assert(_cols() == other.rows());
|
||||
assert(!(Flags & NullDiagBit));
|
||||
|
||||
typename OtherDerived::Eval res(other.rows(), other.cols());
|
||||
|
||||
for (int c=0 ; c<other.cols() ; ++c)
|
||||
{
|
||||
if (Flags & Lower)
|
||||
{
|
||||
// forward substitution
|
||||
if (Flags & UnitDiagBit)
|
||||
res(0,c) = other(0,c);
|
||||
else
|
||||
res(0,c) = other(0,c)/_coeff(0, 0);
|
||||
for (int i=1 ; i<_rows() ; ++i)
|
||||
{
|
||||
Scalar tmp = other(i,c) - ((this->row(i).start(i)) * res.col(c).start(i))(0,0);
|
||||
if (Flags & UnitDiagBit)
|
||||
res(i,c) = tmp;
|
||||
else
|
||||
res(i,c) = tmp/_coeff(i,i);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
// backward substitution
|
||||
if (Flags & UnitDiagBit)
|
||||
res(_cols()-1,c) = other(_cols()-1,c);
|
||||
else
|
||||
res(_cols()-1,c) = other(_cols()-1, c)/_coeff(_rows()-1, _cols()-1);
|
||||
for (int i=_rows()-2 ; i>=0 ; --i)
|
||||
{
|
||||
Scalar tmp = other(i,c) - ((this->row(i).end(_cols()-i-1)) * res.col(c).end(_cols()-i-1))(0,0);
|
||||
if (Flags & UnitDiagBit)
|
||||
res(i,c) = tmp;
|
||||
else
|
||||
res(i,c) = tmp/_coeff(i,i);
|
||||
}
|
||||
}
|
||||
}
|
||||
return res;
|
||||
}
|
||||
EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Extract)
|
||||
|
||||
private:
|
||||
|
||||
inline int _rows() const { return m_matrix.rows(); }
|
||||
inline int _cols() const { return m_matrix.cols(); }
|
||||
|
||||
inline Scalar& _coeffRef(int row, int col)
|
||||
{
|
||||
ei_assert( ((! (Flags & Lower)) && row<=col) || (Flags & Lower && col<=row));
|
||||
return m_matrix.const_cast_derived().coeffRef(row, col);
|
||||
}
|
||||
|
||||
inline Scalar _coeff(int row, int col) const
|
||||
{
|
||||
if ((Flags & Lower) ? col>row : row>col)
|
||||
return 0;
|
||||
if(Flags & LowerTriangularBit ? col>row : row>col)
|
||||
return (Scalar)0;
|
||||
if(Flags & UnitDiagBit)
|
||||
return col==row ? 1 : m_matrix.coeff(row, col);
|
||||
else if (Flags & NullDiagBit)
|
||||
return col==row ? 0 : m_matrix.coeff(row, col);
|
||||
return col==row ? (Scalar)1 : m_matrix.coeff(row, col);
|
||||
else if(Flags & ZeroDiagBit)
|
||||
return col==row ? (Scalar)0 : m_matrix.coeff(row, col);
|
||||
else
|
||||
return m_matrix.coeff(row, col);
|
||||
}
|
||||
@ -180,86 +91,21 @@ template<int Mode, typename MatrixType> class Triangular
|
||||
const typename MatrixType::Nested m_matrix;
|
||||
};
|
||||
|
||||
/** \returns an expression of a upper triangular matrix
|
||||
/** \returns an expression of a triangular matrix extracted from the current matrix
|
||||
*
|
||||
* \sa isUpper(), upperWithNullDiagBit(), upperWithNullDiagBit(), lower()
|
||||
* \sa part(), marked()
|
||||
*/
|
||||
template<typename Derived>
|
||||
inline Triangular<Upper, Derived> MatrixBase<Derived>::upper(void)
|
||||
template<unsigned int Mode>
|
||||
const Extract<Derived, Mode> MatrixBase<Derived>::extract() const
|
||||
{
|
||||
return Triangular<Upper,Derived>(derived());
|
||||
}
|
||||
|
||||
/** This is the const version of upper(). */
|
||||
template<typename Derived>
|
||||
inline const Triangular<Upper, Derived> MatrixBase<Derived>::upper(void) const
|
||||
{
|
||||
return Triangular<Upper,Derived>(derived());
|
||||
}
|
||||
|
||||
/** \returns an expression of a lower triangular matrix
|
||||
*
|
||||
* \sa isLower(), lowerWithUnitDiag(), lowerWithNullDiag(), upper()
|
||||
*/
|
||||
template<typename Derived>
|
||||
inline Triangular<Lower, Derived> MatrixBase<Derived>::lower(void)
|
||||
{
|
||||
return Triangular<Lower,Derived>(derived());
|
||||
}
|
||||
|
||||
/** This is the const version of lower().*/
|
||||
template<typename Derived>
|
||||
inline const Triangular<Lower, Derived> MatrixBase<Derived>::lower(void) const
|
||||
{
|
||||
return Triangular<Lower,Derived>(derived());
|
||||
}
|
||||
|
||||
/** \returns an expression of a upper triangular matrix with a unit diagonal
|
||||
*
|
||||
* \sa upper(), lowerWithUnitDiagBit()
|
||||
*/
|
||||
template<typename Derived>
|
||||
inline const Triangular<Upper|UnitDiagBit, Derived> MatrixBase<Derived>::upperWithUnitDiag(void) const
|
||||
{
|
||||
return Triangular<Upper|UnitDiagBit, Derived>(derived());
|
||||
}
|
||||
|
||||
/** \returns an expression of a strictly upper triangular matrix (diagonal==zero)
|
||||
* FIXME could also be called strictlyUpper() or upperStrict()
|
||||
*
|
||||
* \sa upper(), lowerWithNullDiag()
|
||||
*/
|
||||
template<typename Derived>
|
||||
inline const Triangular<Upper|NullDiagBit, Derived> MatrixBase<Derived>::upperWithNullDiag(void) const
|
||||
{
|
||||
return Triangular<Upper|NullDiagBit, Derived>(derived());
|
||||
}
|
||||
|
||||
/** \returns an expression of a lower triangular matrix with a unit diagonal
|
||||
*
|
||||
* \sa lower(), upperWithUnitDiag()
|
||||
*/
|
||||
template<typename Derived>
|
||||
inline const Triangular<Lower|UnitDiagBit, Derived> MatrixBase<Derived>::lowerWithUnitDiag(void) const
|
||||
{
|
||||
return Triangular<Lower|UnitDiagBit, Derived>(derived());
|
||||
}
|
||||
|
||||
/** \returns an expression of a strictly lower triangular matrix (diagonal==zero)
|
||||
* FIXME could also be called strictlyLower() or lowerStrict()
|
||||
*
|
||||
* \sa lower(), upperWithNullDiag()
|
||||
*/
|
||||
template<typename Derived>
|
||||
inline const Triangular<Lower|NullDiagBit, Derived> MatrixBase<Derived>::lowerWithNullDiag(void) const
|
||||
{
|
||||
return Triangular<Lower|NullDiagBit, Derived>(derived());
|
||||
return derived();
|
||||
}
|
||||
|
||||
/** \returns true if *this is approximately equal to an upper triangular matrix,
|
||||
* within the precision given by \a prec.
|
||||
*
|
||||
* \sa isLower(), upper()
|
||||
* \sa isLower(), extract(), part(), marked()
|
||||
*/
|
||||
template<typename Derived>
|
||||
bool MatrixBase<Derived>::isUpper(RealScalar prec) const
|
||||
@ -281,7 +127,7 @@ bool MatrixBase<Derived>::isUpper(RealScalar prec) const
|
||||
/** \returns true if *this is approximately equal to a lower triangular matrix,
|
||||
* within the precision given by \a prec.
|
||||
*
|
||||
* \sa isUpper(), upper()
|
||||
* \sa isUpper(), extract(), part(), marked()
|
||||
*/
|
||||
template<typename Derived>
|
||||
bool MatrixBase<Derived>::isLower(RealScalar prec) const
|
||||
@ -300,4 +146,4 @@ bool MatrixBase<Derived>::isLower(RealScalar prec) const
|
||||
return true;
|
||||
}
|
||||
|
||||
#endif // EIGEN_TRIANGULAR_H
|
||||
#endif // EIGEN_EXTRACT_H
|
||||
|
||||
@ -1,104 +0,0 @@
|
||||
// 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/>.
|
||||
|
||||
#ifndef EIGEN_TRIANGULAR_ASSIGN_H
|
||||
#define EIGEN_TRIANGULAR_ASSIGN_H
|
||||
|
||||
template<typename Derived1, typename Derived2, int UnrollCount, int Mode>
|
||||
struct ei_triangular_assign_unroller
|
||||
{
|
||||
enum {
|
||||
col = (UnrollCount-1) / Derived1::RowsAtCompileTime,
|
||||
row = (UnrollCount-1) % Derived1::RowsAtCompileTime
|
||||
};
|
||||
|
||||
inline static void run(Derived1 &dst, const Derived2 &src)
|
||||
{
|
||||
ei_triangular_assign_unroller<Derived1, Derived2,
|
||||
(Mode & Lower) ?
|
||||
((row==col) ? UnrollCount-1-row : UnrollCount-1)
|
||||
: ((row==0) ? UnrollCount-1-Derived1::ColsAtCompileTime+col : UnrollCount-1),
|
||||
Mode>::run(dst, src);
|
||||
dst.coeffRef(row, col) = src.coeff(row, col);
|
||||
}
|
||||
};
|
||||
|
||||
template<typename Derived1, typename Derived2, int Mode>
|
||||
struct ei_triangular_assign_unroller<Derived1, Derived2, 1, Mode>
|
||||
{
|
||||
inline static void run(Derived1 &dst, const Derived2 &src)
|
||||
{
|
||||
dst.coeffRef(0, 0) = src.coeff(0, 0);
|
||||
}
|
||||
};
|
||||
|
||||
// prevent buggy user code from causing an infinite recursion
|
||||
template<typename Derived1, typename Derived2, int Mode>
|
||||
struct ei_triangular_assign_unroller<Derived1, Derived2, 0, Mode>
|
||||
{
|
||||
inline static void run(Derived1 &, const Derived2 &) {}
|
||||
};
|
||||
|
||||
template<typename Derived1, typename Derived2, int Mode>
|
||||
struct ei_triangular_assign_unroller<Derived1, Derived2, Dynamic, Mode>
|
||||
{
|
||||
inline static void run(Derived1 &, const Derived2 &) {}
|
||||
};
|
||||
|
||||
|
||||
template <typename Derived, typename OtherDerived, bool DummyVectorize>
|
||||
struct ei_assignment_impl<Derived, OtherDerived, DummyVectorize, true>
|
||||
{
|
||||
static void execute(Derived & dst, const OtherDerived & src)
|
||||
{
|
||||
assert(src.rows()==src.cols());
|
||||
assert(dst.rows() == src.rows() && dst.cols() == src.cols());
|
||||
|
||||
const bool unroll = Derived::SizeAtCompileTime * OtherDerived::CoeffReadCost <= EIGEN_UNROLLING_LIMIT;
|
||||
|
||||
if(unroll)
|
||||
{
|
||||
ei_triangular_assign_unroller
|
||||
<Derived, OtherDerived, unroll ? int(Derived::SizeAtCompileTime) : Dynamic, Derived::Flags>::run
|
||||
(dst.derived(), src.derived());
|
||||
}
|
||||
else
|
||||
{
|
||||
if (Derived::Flags & Lower)
|
||||
{
|
||||
for(int j = 0; j < dst.cols(); j++)
|
||||
for(int i = j; i < dst.rows(); i++)
|
||||
dst.coeffRef(i, j) = src.coeff(i, j);
|
||||
}
|
||||
else
|
||||
{
|
||||
for(int j = 0; j < dst.cols(); j++)
|
||||
for(int i = 0; i <= j; i++)
|
||||
dst.coeffRef(i, j) = src.coeff(i, j);
|
||||
}
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
#endif // EIGEN_TRIANGULAR_ASSIGN_H
|
||||
@ -39,15 +39,13 @@ const unsigned int VectorizableBit = 0x10; ///< means the expression might be v
|
||||
const unsigned int VectorizableBit = 0x0;
|
||||
#endif
|
||||
const unsigned int Like1DArrayBit = 0x20; ///< means the expression can be seen as 1D vector (used for explicit vectorization)
|
||||
const unsigned int NullDiagBit = 0x40; ///< means all diagonal coefficients are equal to 0
|
||||
const unsigned int ZeroDiagBit = 0x40; ///< means all diagonal coefficients are equal to 0
|
||||
const unsigned int UnitDiagBit = 0x80; ///< means all diagonal coefficients are equal to 1
|
||||
const unsigned int NullLowerBit = 0x200; ///< means the strictly triangular lower part is 0
|
||||
const unsigned int NullUpperBit = 0x400; ///< means the strictly triangular upper part is 0
|
||||
const unsigned int SelfAdjointBit = 0x100; ///< means the matrix is selfadjoint (M=M*).
|
||||
const unsigned int UpperTriangularBit = 0x200; ///< means the strictly triangular lower part is 0
|
||||
const unsigned int LowerTriangularBit = 0x400; ///< means the strictly triangular upper part is 0
|
||||
const unsigned int DirectAccessBit = 0x800; ///< means the underlying matrix data can be direclty accessed
|
||||
const unsigned int TemporaryBit = 0x1000; ///< means the expression should be copied by value when nested
|
||||
|
||||
enum { Upper=NullLowerBit, Lower=NullUpperBit };
|
||||
enum { Aligned=0, UnAligned=1 };
|
||||
const unsigned int NestByValueBit = 0x1000; ///< means the expression should be copied by value when nested
|
||||
|
||||
// list of flags that are inherited by default
|
||||
const unsigned int HereditaryBits = RowMajorBit
|
||||
@ -55,6 +53,22 @@ const unsigned int HereditaryBits = RowMajorBit
|
||||
| EvalBeforeAssigningBit
|
||||
| LargeBit;
|
||||
|
||||
// Possible values for the PartType parameter of part() and the ExtractType parameter of extract()
|
||||
const unsigned int Upper = UpperTriangularBit;
|
||||
const unsigned int StrictlyUpper = UpperTriangularBit | ZeroDiagBit;
|
||||
const unsigned int Lower = LowerTriangularBit;
|
||||
const unsigned int StrictlyLower = LowerTriangularBit | ZeroDiagBit;
|
||||
|
||||
// additional possible values for the PartType parameter of part()
|
||||
const unsigned int SelfAdjoint = SelfAdjointBit;
|
||||
|
||||
// additional possible values for the ExtractType parameter of extract()
|
||||
const unsigned int UnitUpper = UpperTriangularBit | UnitDiagBit;
|
||||
const unsigned int UnitLower = LowerTriangularBit | UnitDiagBit;
|
||||
|
||||
|
||||
|
||||
enum { Aligned=0, UnAligned=1 };
|
||||
enum { ConditionalJumpCost = 5 };
|
||||
enum CornerType { TopLeft, TopRight, BottomLeft, BottomRight };
|
||||
enum DirectionType { Vertical, Horizontal };
|
||||
|
||||
@ -46,9 +46,10 @@ template<typename Lhs, typename Rhs, int EvalMode=ei_product_eval_mode<Lhs,Rhs>:
|
||||
template<typename CoeffsVectorType> class DiagonalMatrix;
|
||||
template<typename MatrixType> class DiagonalCoeffs;
|
||||
template<typename MatrixType> class Map;
|
||||
// template<typename Derived> class Eval;
|
||||
template<int Direction, typename UnaryOp, typename MatrixType> class PartialRedux;
|
||||
template<int Mode, typename MatrixType> class Triangular;
|
||||
template<typename MatrixType, unsigned int Mode> class Part;
|
||||
template<typename MatrixType, unsigned int Mode> class Extract;
|
||||
|
||||
|
||||
template<typename Scalar> struct ei_scalar_sum_op;
|
||||
template<typename Scalar> struct ei_scalar_difference_op;
|
||||
|
||||
@ -135,7 +135,7 @@ typedef typename Eigen::NumTraits<Scalar>::Real RealScalar; \
|
||||
typedef typename Base::PacketScalar PacketScalar; \
|
||||
typedef typename Eigen::ei_nested<Derived>::type Nested; \
|
||||
typedef typename Eigen::ei_eval<Derived>::type Eval; \
|
||||
typedef Eigen::Flagged<Derived, Eigen::TemporaryBit, 0> Temporary; \
|
||||
typedef Eigen::Flagged<Derived, NestByValueBit, 0> Temporary; \
|
||||
enum { RowsAtCompileTime = Base::RowsAtCompileTime, \
|
||||
ColsAtCompileTime = Base::ColsAtCompileTime, \
|
||||
MaxRowsAtCompileTime = Base::MaxRowsAtCompileTime, \
|
||||
|
||||
@ -192,15 +192,11 @@ template<typename T> struct ei_unref<T&> { typedef T type; };
|
||||
template<typename T> struct ei_unconst { typedef T type; };
|
||||
template<typename T> struct ei_unconst<const T> { typedef T type; };
|
||||
|
||||
template<typename T> struct ei_is_temporary
|
||||
{
|
||||
enum { ret = int(ei_traits<T>::Flags) & TemporaryBit };
|
||||
};
|
||||
|
||||
template<typename T, int n=1> struct ei_nested
|
||||
{
|
||||
typedef typename ei_meta_if<
|
||||
ei_is_temporary<T>::ret,
|
||||
ei_traits<T>::Flags & NestByValueBit,
|
||||
T,
|
||||
typename ei_meta_if<
|
||||
int(ei_traits<T>::Flags) & EvalBeforeNestingBit
|
||||
|
||||
@ -71,11 +71,11 @@ template<typename Derived>
|
||||
typename ei_traits<Derived>::Scalar MatrixBase<Derived>::determinant() const
|
||||
{
|
||||
assert(rows() == cols());
|
||||
if (Derived::Flags & (NullLowerBit | NullUpperBit))
|
||||
if (Derived::Flags & (UpperTriangularBit | LowerTriangularBit))
|
||||
{
|
||||
if (Derived::Flags & UnitDiagBit)
|
||||
return 1;
|
||||
else if (Derived::Flags & NullDiagBit)
|
||||
else if (Derived::Flags & ZeroDiagBit)
|
||||
return 0;
|
||||
else
|
||||
return derived().diagonal().redux(ei_scalar_product_op<Scalar>());
|
||||
|
||||
@ -111,7 +111,7 @@ template<typename MatrixType>
|
||||
typename QR<MatrixType>::RMatrixType QR<MatrixType>::matrixR(void) const
|
||||
{
|
||||
int cols = m_qr.cols();
|
||||
RMatrixType res = m_qr.block(0,0,cols,cols).upperWithNullDiag();
|
||||
RMatrixType res = m_qr.block(0,0,cols,cols).strictlyUpper();
|
||||
res.diagonal() = m_norms;
|
||||
return res;
|
||||
}
|
||||
|
||||
@ -47,8 +47,8 @@ template<typename MatrixType> void triangular(const MatrixType& m)
|
||||
v2 = VectorType::random(rows),
|
||||
vzero = VectorType::zero(rows);
|
||||
|
||||
MatrixType m1up = m1.upper();
|
||||
MatrixType m2up = m2.upper();
|
||||
MatrixType m1up = m1.template extract<Eigen::Upper>();
|
||||
MatrixType m2up = m2.template extract<Eigen::Upper>();
|
||||
|
||||
if (rows*cols>1)
|
||||
{
|
||||
@ -62,26 +62,26 @@ template<typename MatrixType> void triangular(const MatrixType& m)
|
||||
// test overloaded operator+=
|
||||
r1.setZero();
|
||||
r2.setZero();
|
||||
r1.upper() += m1;
|
||||
r1.template part<Eigen::Upper>() += m1;
|
||||
r2 += m1up;
|
||||
VERIFY_IS_APPROX(r1,r2);
|
||||
|
||||
// test overloaded operator=
|
||||
m1.setZero();
|
||||
m1.upper() = (m2.transpose() * m2).lazy();
|
||||
m1.template part<Eigen::Upper>() = (m2.transpose() * m2).lazy();
|
||||
m3 = m2.transpose() * m2;
|
||||
VERIFY_IS_APPROX(m3.lower().transpose(), m1);
|
||||
VERIFY_IS_APPROX(m3.template extract<Eigen::Lower>().transpose(), m1);
|
||||
|
||||
// test overloaded operator=
|
||||
m1.setZero();
|
||||
m1.lower() = (m2.transpose() * m2).lazy();
|
||||
VERIFY_IS_APPROX(m3.lower(), m1);
|
||||
m1.template part<Eigen::Lower>() = (m2.transpose() * m2).lazy();
|
||||
VERIFY_IS_APPROX(m3.template extract<Eigen::Lower>(), m1);
|
||||
|
||||
// test back and forward subsitution
|
||||
m1 = MatrixType::random(rows, cols);
|
||||
VERIFY_IS_APPROX(m1.upper() * (m1.upper().inverseProduct(m2)), m2);
|
||||
VERIFY_IS_APPROX(m1.lower() * (m1.lower().inverseProduct(m2)), m2);
|
||||
VERIFY((m1.upper() * m2.upper()).isUpper());
|
||||
VERIFY_IS_APPROX(m1.template extract<Eigen::Upper>() * (m1.template extract<Eigen::Upper>().inverseProduct(m2)), m2);
|
||||
VERIFY_IS_APPROX(m1.template extract<Eigen::Lower>() * (m1.template extract<Eigen::Lower>().inverseProduct(m2)), m2);
|
||||
VERIFY((m1.template extract<Eigen::Upper>() * m2.template extract<Eigen::Upper>()).isUpper());
|
||||
|
||||
}
|
||||
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user