devtools/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2026-04-10 11:34:33 +08:00
Code Issues Packages Projects Releases Wiki Activity

synch with main devel branch

Browse Source
This commit is contained in:
Gael Guennebaud
2009-07-14 23:06:25 +02:00
parent 279cedc1ce f5d2317b12
commit 8120a5cecd
46 changed files with 1241 additions and 966 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
Eigen/Core
View File

@@ -88,6 +88,8 @@
#include <cstring>
#include <string>
#include <limits>
// for min/max:
#include <algorithm>
#if (defined(_CPPUNWIND) || defined(__EXCEPTIONS)) && !defined(EIGEN_NO_EXCEPTIONS)
#define EIGEN_EXCEPTIONS
@@ -162,6 +164,7 @@ namespace Eigen {
#include "src/Core/MapBase.h"
#include "src/Core/Map.h"
#include "src/Core/Block.h"
#include "src/Core/VectorBlock.h"
#include "src/Core/Minor.h"
#include "src/Core/Transpose.h"
#include "src/Core/DiagonalMatrix.h"
@@ -182,7 +185,6 @@ namespace Eigen {
#include "src/Core/SolveTriangular.h"
#include "src/Core/products/SelfadjointRank2Update.h"
#include "src/Core/products/TriangularMatrixVector.h"
#include "src/Core/BandMatrix.h"
} // namespace Eigen

1
Eigen/Geometry
View File

@@ -6,6 +6,7 @@
#include "src/Core/util/DisableMSVCWarnings.h"
#include "Array"
#include "SVD"
#include <limits>
#ifndef M_PI

1
Eigen/Sparse
View File

@@ -84,6 +84,7 @@ namespace Eigen {
#include "src/Sparse/SparseUtil.h"
#include "src/Sparse/SparseMatrixBase.h"
#include "src/Sparse/SparseNestByValue.h"
#include "src/Sparse/CompressedStorage.h"
#include "src/Sparse/AmbiVector.h"
#include "src/Sparse/RandomSetter.h"

17
Eigen/src/Array/Replicate.h
View File

@@ -140,21 +140,4 @@ VectorwiseOp<ExpressionType,Direction>::replicate(int factor) const
(_expression(),Direction==Vertical?factor:1,Direction==Horizontal?factor:1);
}
/** \nonstableyet
* \return an expression of the replication of each column (or row) of \c *this
*
* Example: \include DirectionWise_replicate.cpp
* Output: \verbinclude DirectionWise_replicate.out
*
* \sa VectorwiseOp::replicate(int), MatrixBase::replicate(), class Replicate
*/
template<typename ExpressionType, int Direction>
template<int Factor>
const Replicate<ExpressionType,(Direction==Vertical?Factor:1),(Direction==Horizontal?Factor:1)>
VectorwiseOp<ExpressionType,Direction>::replicate(int factor) const
{
return Replicate<ExpressionType,Direction==Vertical?Factor:1,Direction==Horizontal?Factor:1>
(_expression(),Direction==Vertical?factor:1,Direction==Horizontal?factor:1);
}
#endif // EIGEN_REPLICATE_H

52
Eigen/src/Array/VectorwiseOp.h
View File

@@ -179,6 +179,11 @@ template<typename ExpressionType, int Direction> class VectorwiseOp
> Type;
};
enum {
IsVertical = (Direction==Vertical) ? 1 : 0,
IsHorizontal = (Direction==Horizontal) ? 1 : 0
};
protected:
/** \internal
@@ -222,9 +227,17 @@ template<typename ExpressionType, int Direction> class VectorwiseOp
/** \internal */
inline const ExpressionType& _expression() const { return m_matrix; }
/** \returns a row or column vector expression of \c *this reduxed by \a func
*
* The template parameter \a BinaryOp is the type of the functor
* of the custom redux operator. Note that func must be an associative operator.
*
* \sa class VectorwiseOp, MatrixBase::colwise(), MatrixBase::rowwise()
*/
template<typename BinaryOp>
const typename ReduxReturnType<BinaryOp>::Type
redux(const BinaryOp& func = BinaryOp()) const;
redux(const BinaryOp& func = BinaryOp()) const
{ return typename ReduxReturnType<BinaryOp>::Type(_expression(), func); }
/** \returns a row (or column) vector expression of the smallest coefficient
* of each column (or row) of the referenced expression.
@@ -319,16 +332,26 @@ template<typename ExpressionType, int Direction> class VectorwiseOp
*
* \sa MatrixBase::reverse() */
const Reverse<ExpressionType, Direction> reverse() const
{
return Reverse<ExpressionType, Direction>( _expression() );
}
{ return Reverse<ExpressionType, Direction>( _expression() ); }
const Replicate<ExpressionType,Direction==Vertical?Dynamic:1,Direction==Horizontal?Dynamic:1>
replicate(int factor) const;
template<int Factor>
const Replicate<ExpressionType,(Direction==Vertical?Factor:1),(Direction==Horizontal?Factor:1)>
replicate(int factor = Factor) const;
/** \nonstableyet
* \return an expression of the replication of each column (or row) of \c *this
*
* Example: \include DirectionWise_replicate.cpp
* Output: \verbinclude DirectionWise_replicate.out
*
* \sa VectorwiseOp::replicate(int), MatrixBase::replicate(), class Replicate
*/
// NOTE implemented here because of sunstudio's compilation errors
template<int Factor> const Replicate<ExpressionType,(IsVertical?Factor:1),(IsHorizontal?Factor:1)>
replicate(int factor = Factor) const
{
return Replicate<ExpressionType,Direction==Vertical?Factor:1,Direction==Horizontal?Factor:1>
(_expression(),Direction==Vertical?factor:1,Direction==Horizontal?factor:1);
}
/////////// Artithmetic operators ///////////
@@ -466,19 +489,4 @@ MatrixBase<Derived>::rowwise()
return derived();
}
/** \returns a row or column vector expression of \c *this reduxed by \a func
*
* The template parameter \a BinaryOp is the type of the functor
* of the custom redux operator. Note that func must be an associative operator.
*
* \sa class VectorwiseOp, MatrixBase::colwise(), MatrixBase::rowwise()
*/
template<typename ExpressionType, int Direction>
template<typename BinaryOp>
const typename VectorwiseOp<ExpressionType,Direction>::template ReduxReturnType<BinaryOp>::Type
VectorwiseOp<ExpressionType,Direction>::redux(const BinaryOp& func) const
{
return typename ReduxReturnType<BinaryOp>::Type(_expression(), func);
}
#endif // EIGEN_PARTIAL_REDUX_H

249
Eigen/src/Core/Block.h
View File

@@ -271,13 +271,19 @@ class Block<MatrixType,BlockRows,BlockCols,PacketAccess,HasDirectAccess>
inline int stride(void) const { return m_matrix.stride(); }
#ifndef __SUNPRO_CC
// FIXME sunstudio is not friendly with the above friend...
protected:
#endif
#ifndef EIGEN_PARSED_BY_DOXYGEN
/** \internal used by allowAligned() */
inline Block(const MatrixType& matrix, const Scalar* data, int blockRows, int blockCols)
: Base(data, blockRows, blockCols), m_matrix(matrix)
{}
#endif
protected:
const typename MatrixType::Nested m_matrix;
};
@@ -314,249 +320,6 @@ inline const typename BlockReturnType<Derived>::Type MatrixBase<Derived>
return typename BlockReturnType<Derived>::Type(derived(), startRow, startCol, blockRows, blockCols);
}
/** \returns a dynamic-size expression of a segment (i.e. a vector block) in *this.
*
* \only_for_vectors
*
* \addexample SegmentIntInt \label How to reference a sub-vector (dynamic size)
*
* \param start the first coefficient in the segment
* \param size the number of coefficients in the segment
*
* Example: \include MatrixBase_segment_int_int.cpp
* Output: \verbinclude MatrixBase_segment_int_int.out
*
* \note Even though the returned expression has dynamic size, in the case
* when it is applied to a fixed-size vector, it inherits a fixed maximal size,
* which means that evaluating it does not cause a dynamic memory allocation.
*
* \sa class Block, segment(int)
*/
template<typename Derived>
inline typename BlockReturnType<Derived>::SubVectorType MatrixBase<Derived>
::segment(int start, int size)
{
EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
return typename BlockReturnType<Derived>::SubVectorType(derived(), RowsAtCompileTime == 1 ? 0 : start,
ColsAtCompileTime == 1 ? 0 : start,
RowsAtCompileTime == 1 ? 1 : size,
ColsAtCompileTime == 1 ? 1 : size);
}
/** This is the const version of segment(int,int).*/
template<typename Derived>
inline const typename BlockReturnType<Derived>::SubVectorType
MatrixBase<Derived>::segment(int start, int size) const
{
EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
return typename BlockReturnType<Derived>::SubVectorType(derived(), RowsAtCompileTime == 1 ? 0 : start,
ColsAtCompileTime == 1 ? 0 : start,
RowsAtCompileTime == 1 ? 1 : size,
ColsAtCompileTime == 1 ? 1 : size);
}
/** \returns a dynamic-size expression of the first coefficients of *this.
*
* \only_for_vectors
*
* \param size the number of coefficients in the block
*
* \addexample BlockInt \label How to reference a sub-vector (fixed-size)
*
* Example: \include MatrixBase_start_int.cpp
* Output: \verbinclude MatrixBase_start_int.out
*
* \note Even though the returned expression has dynamic size, in the case
* when it is applied to a fixed-size vector, it inherits a fixed maximal size,
* which means that evaluating it does not cause a dynamic memory allocation.
*
* \sa class Block, block(int,int)
*/
template<typename Derived>
inline typename BlockReturnType<Derived,Dynamic>::SubVectorType
MatrixBase<Derived>::start(int size)
{
EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
return Block<Derived,
RowsAtCompileTime == 1 ? 1 : Dynamic,
ColsAtCompileTime == 1 ? 1 : Dynamic>
(derived(), 0, 0,
RowsAtCompileTime == 1 ? 1 : size,
ColsAtCompileTime == 1 ? 1 : size);
}
/** This is the const version of start(int).*/
template<typename Derived>
inline const typename BlockReturnType<Derived,Dynamic>::SubVectorType
MatrixBase<Derived>::start(int size) const
{
EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
return Block<Derived,
RowsAtCompileTime == 1 ? 1 : Dynamic,
ColsAtCompileTime == 1 ? 1 : Dynamic>
(derived(), 0, 0,
RowsAtCompileTime == 1 ? 1 : size,
ColsAtCompileTime == 1 ? 1 : size);
}
/** \returns a dynamic-size expression of the last coefficients of *this.
*
* \only_for_vectors
*
* \param size the number of coefficients in the block
*
* \addexample BlockEnd \label How to reference the end of a vector (fixed-size)
*
* Example: \include MatrixBase_end_int.cpp
* Output: \verbinclude MatrixBase_end_int.out
*
* \note Even though the returned expression has dynamic size, in the case
* when it is applied to a fixed-size vector, it inherits a fixed maximal size,
* which means that evaluating it does not cause a dynamic memory allocation.
*
* \sa class Block, block(int,int)
*/
template<typename Derived>
inline typename BlockReturnType<Derived,Dynamic>::SubVectorType
MatrixBase<Derived>::end(int size)
{
EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
return Block<Derived,
RowsAtCompileTime == 1 ? 1 : Dynamic,
ColsAtCompileTime == 1 ? 1 : Dynamic>
(derived(),
RowsAtCompileTime == 1 ? 0 : rows() - size,
ColsAtCompileTime == 1 ? 0 : cols() - size,
RowsAtCompileTime == 1 ? 1 : size,
ColsAtCompileTime == 1 ? 1 : size);
}
/** This is the const version of end(int).*/
template<typename Derived>
inline const typename BlockReturnType<Derived,Dynamic>::SubVectorType
MatrixBase<Derived>::end(int size) const
{
EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
return Block<Derived,
RowsAtCompileTime == 1 ? 1 : Dynamic,
ColsAtCompileTime == 1 ? 1 : Dynamic>
(derived(),
RowsAtCompileTime == 1 ? 0 : rows() - size,
ColsAtCompileTime == 1 ? 0 : cols() - size,
RowsAtCompileTime == 1 ? 1 : size,
ColsAtCompileTime == 1 ? 1 : size);
}
/** \returns a fixed-size expression of a segment (i.e. a vector block) in \c *this
*
* \only_for_vectors
*
* The template parameter \a Size is the number of coefficients in the block
*
* \param start the index of the first element of the sub-vector
*
* Example: \include MatrixBase_template_int_segment.cpp
* Output: \verbinclude MatrixBase_template_int_segment.out
*
* \sa class Block
*/
template<typename Derived>
template<int Size>
inline typename BlockReturnType<Derived,Size>::SubVectorType
MatrixBase<Derived>::segment(int start)
{
EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
return Block<Derived, (RowsAtCompileTime == 1 ? 1 : Size),
(ColsAtCompileTime == 1 ? 1 : Size)>
(derived(), RowsAtCompileTime == 1 ? 0 : start,
ColsAtCompileTime == 1 ? 0 : start);
}
/** This is the const version of segment<int>(int).*/
template<typename Derived>
template<int Size>
inline const typename BlockReturnType<Derived,Size>::SubVectorType
MatrixBase<Derived>::segment(int start) const
{
EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
return Block<Derived, (RowsAtCompileTime == 1 ? 1 : Size),
(ColsAtCompileTime == 1 ? 1 : Size)>
(derived(), RowsAtCompileTime == 1 ? 0 : start,
ColsAtCompileTime == 1 ? 0 : start);
}
/** \returns a fixed-size expression of the first coefficients of *this.
*
* \only_for_vectors
*
* The template parameter \a Size is the number of coefficients in the block
*
* \addexample BlockStart \label How to reference the start of a vector (fixed-size)
*
* Example: \include MatrixBase_template_int_start.cpp
* Output: \verbinclude MatrixBase_template_int_start.out
*
* \sa class Block
*/
template<typename Derived>
template<int Size>
inline typename BlockReturnType<Derived,Size>::SubVectorType
MatrixBase<Derived>::start()
{
EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
return Block<Derived, (RowsAtCompileTime == 1 ? 1 : Size),
(ColsAtCompileTime == 1 ? 1 : Size)>(derived(), 0, 0);
}
/** This is the const version of start<int>().*/
template<typename Derived>
template<int Size>
inline const typename BlockReturnType<Derived,Size>::SubVectorType
MatrixBase<Derived>::start() const
{
EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
return Block<Derived, (RowsAtCompileTime == 1 ? 1 : Size),
(ColsAtCompileTime == 1 ? 1 : Size)>(derived(), 0, 0);
}
/** \returns a fixed-size expression of the last coefficients of *this.
*
* \only_for_vectors
*
* The template parameter \a Size is the number of coefficients in the block
*
* Example: \include MatrixBase_template_int_end.cpp
* Output: \verbinclude MatrixBase_template_int_end.out
*
* \sa class Block
*/
template<typename Derived>
template<int Size>
inline typename BlockReturnType<Derived,Size>::SubVectorType
MatrixBase<Derived>::end()
{
EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
return Block<Derived, RowsAtCompileTime == 1 ? 1 : Size,
ColsAtCompileTime == 1 ? 1 : Size>
(derived(),
RowsAtCompileTime == 1 ? 0 : rows() - Size,
ColsAtCompileTime == 1 ? 0 : cols() - Size);
}
/** This is the const version of end<int>.*/
template<typename Derived>
template<int Size>
inline const typename BlockReturnType<Derived,Size>::SubVectorType
MatrixBase<Derived>::end() const
{
EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
return Block<Derived, RowsAtCompileTime == 1 ? 1 : Size,
ColsAtCompileTime == 1 ? 1 : Size>
(derived(),
RowsAtCompileTime == 1 ? 0 : rows() - Size,
ColsAtCompileTime == 1 ? 0 : cols() - Size);
}
/** \returns a dynamic-size expression of a corner of *this.
*
* \param type the type of corner. Can be \a Eigen::TopLeft, \a Eigen::TopRight,

15
Eigen/src/Core/CwiseUnaryView.h
View File

@@ -42,13 +42,13 @@ struct ei_traits<CwiseUnaryView<ViewOp, MatrixType> >
: ei_traits<MatrixType>
{
typedef typename ei_result_of<
ViewOp(typename MatrixType::Scalar)
ViewOp(typename ei_traits<MatrixType>::Scalar)
>::type Scalar;
typedef typename MatrixType::Nested MatrixTypeNested;
typedef typename ei_unref<MatrixTypeNested>::type _MatrixTypeNested;
typedef typename ei_cleantype<MatrixTypeNested>::type _MatrixTypeNested;
enum {
Flags = (_MatrixTypeNested::Flags & (HereditaryBits | LinearAccessBit | AlignedBit)),
CoeffReadCost = _MatrixTypeNested::CoeffReadCost + ei_functor_traits<ViewOp>::Cost
Flags = (ei_traits<_MatrixTypeNested>::Flags & (HereditaryBits | LinearAccessBit | AlignedBit)),
CoeffReadCost = ei_traits<_MatrixTypeNested>::CoeffReadCost + ei_functor_traits<ViewOp>::Cost
};
};
@@ -62,7 +62,7 @@ class CwiseUnaryView : ei_no_assignment_operator,
inline CwiseUnaryView(const MatrixType& mat, const ViewOp& func = ViewOp())
: m_matrix(mat), m_functor(func) {}
EIGEN_INHERIT_ASSIGNMENT_OPERATORS(CwiseUnaryView)
EIGEN_STRONG_INLINE int rows() const { return m_matrix.rows(); }
@@ -77,7 +77,7 @@ class CwiseUnaryView : ei_no_assignment_operator,
{
return m_functor(m_matrix.coeff(index));
}
EIGEN_STRONG_INLINE Scalar& coeffRef(int row, int col)
{
return m_functor(m_matrix.const_cast_derived().coeffRef(row, col));
@@ -89,7 +89,8 @@ class CwiseUnaryView : ei_no_assignment_operator,
}
protected:
const typename MatrixType::Nested m_matrix;
// FIXME changed from MatrixType::Nested because of a weird compilation error with sun CC
const typename ei_nested<MatrixType>::type m_matrix;
const ViewOp m_functor;
};

138
Eigen/src/Core/Dot.h
View File

@@ -295,7 +295,7 @@ inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real MatrixBase<
/** \returns the \em l2 norm of \c *this using a numerically more stable
* algorithm.
*
* \sa norm(), dot(), squaredNorm()
* \sa norm(), dot(), squaredNorm(), blueNorm()
*/
template<typename Derived>
inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real
@@ -304,6 +304,142 @@ MatrixBase<Derived>::stableNorm() const
return this->cwise().abs().redux(ei_scalar_hypot_op<RealScalar>());
}
/** \internal Computes ibeta^iexp by binary expansion of iexp,
* exact if ibeta is the machine base */
template<typename T> inline T bexp(int ibeta, int iexp)
{
T tbeta = T(ibeta);
T res = 1.0;
int n = iexp;
if (n<0)
{
n = - n;
tbeta = 1.0/tbeta;
}
for(;;)
{
if ((n % 2)==0)
res = res * tbeta;
n = n/2;
if (n==0) return res;
tbeta = tbeta*tbeta;
}
return res;
}
/** \returns the \em l2 norm of \c *this using the Blue's algorithm.
* A Portable Fortran Program to Find the Euclidean Norm of a Vector,
* ACM TOMS, Vol 4, Issue 1, 1978.
*
* \sa norm(), dot(), squaredNorm(), stableNorm()
*/
template<typename Derived>
inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real
MatrixBase<Derived>::blueNorm() const
{
static int nmax;
static Scalar b1, b2, s1m, s2m, overfl, rbig, relerr;
int n;
Scalar ax, abig, amed, asml;
if(nmax <= 0)
{
int nbig, ibeta, it, iemin, iemax, iexp;
Scalar abig, eps;
// This program calculates the machine-dependent constants
// bl, b2, slm, s2m, relerr overfl, nmax
// from the "basic" machine-dependent numbers
// nbig, ibeta, it, iemin, iemax, rbig.
// The following define the basic machine-dependent constants.
// For portability, the PORT subprograms "ilmaeh" and "rlmach"
// are used. For any specific computer, each of the assignment
// statements can be replaced
nbig = std::numeric_limits<int>::max(); // largest integer
ibeta = NumTraits<Scalar>::Base; // base for floating-point numbers
it = NumTraits<Scalar>::Mantissa; // number of base-beta digits in mantissa
iemin = std::numeric_limits<Scalar>::min_exponent; // minimum exponent
iemax = std::numeric_limits<Scalar>::max_exponent; // maximum exponent
rbig = std::numeric_limits<Scalar>::max(); // largest floating-point number
// Check the basic machine-dependent constants.
if(iemin > 1 - 2*it || 1+it>iemax || (it==2 && ibeta<5)
|| (it<=4 && ibeta <= 3 ) || it<2)
{
ei_assert(false && "the algorithm cannot be guaranteed on this computer");
}
iexp = -((1-iemin)/2);
b1 = bexp<Scalar>(ibeta, iexp); // lower boundary of midrange
iexp = (iemax + 1 - it)/2;
b2 = bexp<Scalar>(ibeta,iexp); // upper boundary of midrange
iexp = (2-iemin)/2;
s1m = bexp<Scalar>(ibeta,iexp); // scaling factor for lower range
iexp = - ((iemax+it)/2);
s2m = bexp<Scalar>(ibeta,iexp); // scaling factor for upper range
overfl = rbig*s2m; // overfow boundary for abig
eps = bexp<Scalar>(ibeta, 1-it);
relerr = ei_sqrt(eps); // tolerance for neglecting asml
abig = 1.0/eps - 1.0;
if (Scalar(nbig)>abig) nmax = abig; // largest safe n
else nmax = nbig;
}
n = size();
if(n==0)
return 0;
asml = Scalar(0);
amed = Scalar(0);
abig = Scalar(0);
for(int j=0; j<n; ++j)
{
ax = ei_abs(coeff(j));
if(ax > b2) abig += ei_abs2(ax*s2m);
else if(ax < b2) asml += ei_abs2(ax*s1m);
else amed += ei_abs2(ax);
}
if(abig > Scalar(0))
{
abig = ei_sqrt(abig);
if(abig > overfl)
{
ei_assert(false && "overflow");
return rbig;
}
if(amed > Scalar(0))
{
abig = abig/s2m;
amed = ei_sqrt(amed);
}
else
{
return abig/s2m;
}
}
else if(asml > Scalar(0))
{
if (amed > Scalar(0))
{
abig = ei_sqrt(amed);
amed = ei_sqrt(asml) / s1m;
}
else
{
return ei_sqrt(asml)/s1m;
}
}
else
{
return ei_sqrt(amed);
}
asml = std::min(abig, amed);
abig = std::max(abig, amed);
if(asml <= abig*relerr)
return abig;
else
return abig * ei_sqrt(Scalar(1) + ei_abs2(asml/abig));
}
/** \returns an expression of the quotient of *this by its own norm.
*
* \only_for_vectors

2
Eigen/src/Core/MathFunctions.h
View File

@@ -151,7 +151,7 @@ inline bool ei_isApproxOrLessThan(float a, float b, float prec = precision<float
*** double ***
**************/
template<> inline double precision<double>() { return 1e-11; }
template<> inline double precision<double>() { return 1e-12; }
template<> inline double machine_epsilon<double>() { return 2.220e-16; }
inline double ei_real(double x) { return x; }

4
Eigen/src/Core/Matrix.h
View File

@@ -124,7 +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;
@@ -218,7 +218,7 @@ class Matrix
*
* This method is intended for dynamic-size matrices, although it is legal to call it on any
* matrix as long as fixed dimensions are left unchanged. If you only want to change the number
* of rows and/or of columns, you can use resize(NoChange_t, int), resize(int, NoChange_t).
* of rows and/or of columns, you can use resize(NoChange_t, int), resize(int, NoChange_t).
*
* If the current number of coefficients of \c *this exactly matches the
* product \a rows * \a cols, then no memory allocation is performed and

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

@@ -163,10 +163,14 @@ template<typename Derived> class MatrixBase
* constructed from this one. See the \ref flags "list of flags".
*/
CoeffReadCost = ei_traits<Derived>::CoeffReadCost
CoeffReadCost = ei_traits<Derived>::CoeffReadCost,
/**< This is a rough measure of how expensive it is to read one coefficient from
* this expression.
*/
#ifndef EIGEN_PARSED_BY_DOXYGEN
_HasDirectAccess = (int(Flags)&DirectAccessBit) ? 1 : 0 // workaround sunCC
#endif
};
/** Default constructor. Just checks at compile-time for self-consistency of the flags. */
@@ -230,7 +234,7 @@ template<typename Derived> class MatrixBase
/** \internal the return type of coeff()
*/
typedef typename ei_meta_if<bool(int(Flags)&DirectAccessBit), const Scalar&, Scalar>::ret CoeffReturnType;
typedef typename ei_meta_if<_HasDirectAccess, const Scalar&, Scalar>::ret CoeffReturnType;
/** \internal Represents a matrix with all coefficients equal to one another*/
typedef CwiseNullaryOp<ei_scalar_constant_op<Scalar>,Derived> ConstantReturnType;
@@ -426,8 +430,9 @@ template<typename Derived> class MatrixBase
template<typename OtherDerived>
Scalar dot(const MatrixBase<OtherDerived>& other) const;
RealScalar squaredNorm() const;
RealScalar norm() const;
RealScalar stableNorm() const;
RealScalar norm() const;
RealScalar stableNorm() const;
RealScalar blueNorm() const;
const PlainMatrixType normalized() const;
void normalize();
@@ -450,14 +455,14 @@ template<typename Derived> class MatrixBase
const typename BlockReturnType<Derived>::Type
block(int startRow, int startCol, int blockRows, int blockCols) const;
typename BlockReturnType<Derived>::SubVectorType segment(int start, int size);
const typename BlockReturnType<Derived>::SubVectorType segment(int start, int size) const;
VectorBlock<Derived> segment(int start, int size);
const VectorBlock<Derived> segment(int start, int size) const;
typename BlockReturnType<Derived,Dynamic>::SubVectorType start(int size);
const typename BlockReturnType<Derived,Dynamic>::SubVectorType start(int size) const;
VectorBlock<Derived> start(int size);
const VectorBlock<Derived> start(int size) const;
typename BlockReturnType<Derived,Dynamic>::SubVectorType end(int size);
const typename BlockReturnType<Derived,Dynamic>::SubVectorType end(int size) const;
VectorBlock<Derived> end(int size);
const VectorBlock<Derived> end(int size) const;
typename BlockReturnType<Derived>::Type corner(CornerType type, int cRows, int cCols);
const typename BlockReturnType<Derived>::Type corner(CornerType type, int cRows, int cCols) const;
@@ -472,14 +477,14 @@ template<typename Derived> class MatrixBase
template<int CRows, int CCols>
const typename BlockReturnType<Derived, CRows, CCols>::Type corner(CornerType type) const;
template<int Size> typename BlockReturnType<Derived,Size>::SubVectorType start(void);
template<int Size> const typename BlockReturnType<Derived,Size>::SubVectorType start() const;
template<int Size> VectorBlock<Derived,Size> start(void);
template<int Size> const VectorBlock<Derived,Size> start() const;
template<int Size> typename BlockReturnType<Derived,Size>::SubVectorType end();
template<int Size> const typename BlockReturnType<Derived,Size>::SubVectorType end() const;
template<int Size> VectorBlock<Derived,Size> end();
template<int Size> const VectorBlock<Derived,Size> end() const;
template<int Size> typename BlockReturnType<Derived,Size>::SubVectorType segment(int start);
template<int Size> const typename BlockReturnType<Derived,Size>::SubVectorType segment(int start) const;
template<int Size> VectorBlock<Derived,Size> segment(int start);
template<int Size> const VectorBlock<Derived,Size> segment(int start) const;
Diagonal<Derived,0> diagonal();
const Diagonal<Derived,0> diagonal() const;
@@ -696,6 +701,7 @@ template<typename Derived> class MatrixBase
const PartialLU<PlainMatrixType> partialLu() const;
const PlainMatrixType inverse() const;
void computeInverse(PlainMatrixType *result) const;
bool computeInverseWithCheck( PlainMatrixType *result ) const;
Scalar determinant() const;
/////////// Cholesky module ///////////
@@ -705,7 +711,7 @@ template<typename Derived> class MatrixBase
/////////// QR module ///////////
const QR<PlainMatrixType> qr() const;
const HouseholderQR<PlainMatrixType> householderQr() const;
EigenvaluesReturnType eigenvalues() const;
RealScalar operatorNorm() const;

8
Eigen/src/Core/NumTraits.h
View File

@@ -70,7 +70,9 @@ template<> struct NumTraits<float>
HasFloatingPoint = 1,
ReadCost = 1,
AddCost = 1,
MulCost = 1
MulCost = 1,
Base = 2,
Mantissa = 23
};
};
@@ -83,7 +85,9 @@ template<> struct NumTraits<double>
HasFloatingPoint = 1,
ReadCost = 1,
AddCost = 1,
MulCost = 1
MulCost = 1,
Base = 2,
Mantissa = 52
};
};

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

@@ -155,13 +155,14 @@ template<typename Lhs, typename Rhs> struct ei_product_mode
typedef typename ei_blas_traits<Lhs>::ActualXprType ActualLhs;
typedef typename ei_blas_traits<Rhs>::ActualXprType ActualRhs;
enum{
value = Lhs::MaxColsAtCompileTime == Dynamic
&& ( Lhs::MaxRowsAtCompileTime == Dynamic
|| Rhs::MaxColsAtCompileTime == Dynamic )
&& (!(Rhs::IsVectorAtCompileTime && (Lhs::Flags&RowMajorBit) && (!(ActualLhs::Flags&DirectAccessBit))))
&& (!(Lhs::IsVectorAtCompileTime && (!(Rhs::Flags&RowMajorBit)) && (!(ActualRhs::Flags&DirectAccessBit))))
&& (ei_is_same_type<typename Lhs::Scalar, typename Rhs::Scalar>::ret)
// workaround sun studio:
LhsIsVectorAtCompileTime = ei_traits<Lhs>::ColsAtCompileTime==1 || ei_traits<Rhs>::ColsAtCompileTime==1,
value = ei_traits<Lhs>::MaxColsAtCompileTime == Dynamic
&& ( ei_traits<Lhs>::MaxRowsAtCompileTime == Dynamic
|| ei_traits<Rhs>::MaxColsAtCompileTime == Dynamic )
&& (!(Rhs::IsVectorAtCompileTime && (ei_traits<Lhs>::Flags&RowMajorBit) && (!(ei_traits<Lhs>::Flags&DirectAccessBit))))
&& (!(LhsIsVectorAtCompileTime && (!(ei_traits<Rhs>::Flags&RowMajorBit)) && (!(ei_traits<Rhs>::Flags&DirectAccessBit))))
&& (ei_is_same_type<typename ei_traits<Lhs>::Scalar, typename ei_traits<Rhs>::Scalar>::ret)
? CacheFriendlyProduct
: NormalProduct };
};
@@ -577,7 +578,7 @@ struct ei_product_packet_impl<ColMajor, Dynamic, Lhs, Rhs, PacketScalar, LoadMod
// Forward declarations
template<typename Scalar, bool ConjugateLhs, bool ConjugateRhs>
void ei_cache_friendly_product(
static void ei_cache_friendly_product(
int _rows, int _cols, int depth,
bool _lhsRowMajor, const Scalar* _lhs, int _lhsStride,
bool _rhsRowMajor, const Scalar* _rhs, int _rhsStride,

2
Eigen/src/Core/SolveTriangular.h
View File

@@ -180,7 +180,7 @@ struct ei_triangular_solver_unroller<Lhs,Rhs,Mode,Index,Size,false> {
template<typename Lhs, typename Rhs, int Mode, int Index, int Size>
struct ei_triangular_solver_unroller<Lhs,Rhs,Mode,Index,Size,true> {
static void run(const Lhs& lhs, Rhs& rhs) {}
static void run(const Lhs&, Rhs&) {}
};
template<typename Lhs, typename Rhs, int Mode, int StorageOrder>

311
Eigen/src/Core/VectorBlock.h Normal file
View File

@@ -0,0 +1,311 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008-2009 Gael Guennebaud <g.gael@free.fr>
// Copyright (C) 2006-2008 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/>.
#ifndef EIGEN_VECTORBLOCK_H
#define EIGEN_VECTORBLOCK_H
/** \class VectorBlock
*
* \brief Expression of a fixed-size or dynamic-size sub-vector
*
* \param VectorType the type of the object in which we are taking a sub-vector
* \param Size size of the sub-vector we are taking at compile time (optional)
* \param _PacketAccess allows to enforce aligned loads and stores if set to ForceAligned.
* The default is AsRequested. This parameter is internaly used by Eigen
* in expressions such as \code mat.segment() += other; \endcode and most of
* the time this is the only way it is used.
*
* This class represents an expression of either a fixed-size or dynamic-size sub-vector.
* It is the return type of MatrixBase::segment(int,int) and MatrixBase::segment<int>(int) and
* most of the time this is the only way it is used.
*
* However, if you want to directly maniputate sub-vector expressions,
* for instance if you want to write a function returning such an expression, you
* will need to use this class.
*
* Here is an example illustrating the dynamic case:
* \include class_VectorBlock.cpp
* Output: \verbinclude class_VectorBlock.out
*
* \note Even though this expression has dynamic size, in the case where \a VectorType
* has fixed size, this expression inherits a fixed maximal size which means that evaluating
* it does not cause a dynamic memory allocation.
*
* Here is an example illustrating the fixed-size case:
* \include class_FixedVectorBlock.cpp
* Output: \verbinclude class_FixedVectorBlock.out
*
* \sa class Block, MatrixBase::segment(int,int,int,int), MatrixBase::segment(int,int)
*/
template<typename VectorType, int Size, int _PacketAccess>
struct ei_traits<VectorBlock<VectorType, Size, _PacketAccess> >
: public ei_traits<Block<VectorType,
ei_traits<VectorType>::RowsAtCompileTime==1 ? 1 : Size,
ei_traits<VectorType>::ColsAtCompileTime==1 ? 1 : Size,
_PacketAccess> >
{
};
template<typename VectorType, int Size, int PacketAccess> class VectorBlock
: public Block<VectorType,
ei_traits<VectorType>::RowsAtCompileTime==1 ? 1 : Size,
ei_traits<VectorType>::ColsAtCompileTime==1 ? 1 : Size,
PacketAccess>
{
typedef Block<VectorType,
ei_traits<VectorType>::RowsAtCompileTime==1 ? 1 : Size,
ei_traits<VectorType>::ColsAtCompileTime==1 ? 1 : Size,
PacketAccess> Base;
enum {
IsColVector = ei_traits<VectorType>::ColsAtCompileTime==1
};
public:
using Base::operator=;
using Base::operator+=;
using Base::operator-=;
using Base::operator*=;
using Base::operator/=;
/** Dynamic-size constructor
*/
inline VectorBlock(const VectorType& vector, int start, int size)
: Base(vector,
IsColVector ? start : 0, IsColVector ? 0 : start,
IsColVector ? size : 1, IsColVector ? 1 : size)
{
EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorBlock);
}
/** Fixed-size constructor
*/
inline VectorBlock(const VectorType& vector, int start)
: Base(vector, IsColVector ? start : 0, IsColVector ? 0 : start)
{
EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorBlock);
}
};
/** \returns a dynamic-size expression of a segment (i.e. a vector block) in *this.
*
* \only_for_vectors
*
* \addexample VectorBlockIntInt \label How to reference a sub-vector (dynamic size)
*
* \param start the first coefficient in the segment
* \param size the number of coefficients in the segment
*
* Example: \include MatrixBase_segment_int_int.cpp
* Output: \verbinclude MatrixBase_segment_int_int.out
*
* \note Even though the returned expression has dynamic size, in the case
* when it is applied to a fixed-size vector, it inherits a fixed maximal size,
* which means that evaluating it does not cause a dynamic memory allocation.
*
* \sa class Block, segment(int)
*/
template<typename Derived>
inline VectorBlock<Derived> MatrixBase<Derived>
::segment(int start, int size)
{
EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
return VectorBlock<Derived>(derived(), start, size);
}
/** This is the const version of segment(int,int).*/
template<typename Derived>
inline const VectorBlock<Derived>
MatrixBase<Derived>::segment(int start, int size) const
{
EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
return VectorBlock<Derived>(derived(), start, size);
}
/** \returns a dynamic-size expression of the first coefficients of *this.
*
* \only_for_vectors
*
* \param size the number of coefficients in the block
*
* \addexample BlockInt \label How to reference a sub-vector (fixed-size)
*
* Example: \include MatrixBase_start_int.cpp
* Output: \verbinclude MatrixBase_start_int.out
*
* \note Even though the returned expression has dynamic size, in the case
* when it is applied to a fixed-size vector, it inherits a fixed maximal size,
* which means that evaluating it does not cause a dynamic memory allocation.
*
* \sa class Block, block(int,int)
*/
template<typename Derived>
inline VectorBlock<Derived>
MatrixBase<Derived>::start(int size)
{
EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
return VectorBlock<Derived>(derived(), 0, size);
}
/** This is the const version of start(int).*/
template<typename Derived>
inline const VectorBlock<Derived>
MatrixBase<Derived>::start(int size) const
{
EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
return VectorBlock<Derived>(derived(), 0, size);
}
/** \returns a dynamic-size expression of the last coefficients of *this.
*
* \only_for_vectors
*
* \param size the number of coefficients in the block
*
* \addexample BlockEnd \label How to reference the end of a vector (fixed-size)
*
* Example: \include MatrixBase_end_int.cpp
* Output: \verbinclude MatrixBase_end_int.out
*
* \note Even though the returned expression has dynamic size, in the case
* when it is applied to a fixed-size vector, it inherits a fixed maximal size,
* which means that evaluating it does not cause a dynamic memory allocation.
*
* \sa class Block, block(int,int)
*/
template<typename Derived>
inline VectorBlock<Derived>
MatrixBase<Derived>::end(int size)
{
EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
return VectorBlock<Derived>(derived(), this->size() - size, size);
}
/** This is the const version of end(int).*/
template<typename Derived>
inline const VectorBlock<Derived>
MatrixBase<Derived>::end(int size) const
{
EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
return VectorBlock<Derived>(derived(), this->size() - size, size);
}
/** \returns a fixed-size expression of a segment (i.e. a vector block) in \c *this
*
* \only_for_vectors
*
* The template parameter \a Size is the number of coefficients in the block
*
* \param start the index of the first element of the sub-vector
*
* Example: \include MatrixBase_template_int_segment.cpp
* Output: \verbinclude MatrixBase_template_int_segment.out
*
* \sa class Block
*/
template<typename Derived>
template<int Size>
inline VectorBlock<Derived,Size>
MatrixBase<Derived>::segment(int start)
{
EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
return VectorBlock<Derived,Size>(derived(), start);
}
/** This is the const version of segment<int>(int).*/
template<typename Derived>
template<int Size>
inline const VectorBlock<Derived,Size>
MatrixBase<Derived>::segment(int start) const
{
EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
return VectorBlock<Derived,Size>(derived(), start);
}
/** \returns a fixed-size expression of the first coefficients of *this.
*
* \only_for_vectors
*
* The template parameter \a Size is the number of coefficients in the block
*
* \addexample BlockStart \label How to reference the start of a vector (fixed-size)
*
* Example: \include MatrixBase_template_int_start.cpp
* Output: \verbinclude MatrixBase_template_int_start.out
*
* \sa class Block
*/
template<typename Derived>
template<int Size>
inline VectorBlock<Derived,Size>
MatrixBase<Derived>::start()
{
EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
return VectorBlock<Derived,Size>(derived(), 0);
}
/** This is the const version of start<int>().*/
template<typename Derived>
template<int Size>
inline const VectorBlock<Derived,Size>
MatrixBase<Derived>::start() const
{
EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
return VectorBlock<Derived,Size>(derived(), 0);
}
/** \returns a fixed-size expression of the last coefficients of *this.
*
* \only_for_vectors
*
* The template parameter \a Size is the number of coefficients in the block
*
* Example: \include MatrixBase_template_int_end.cpp
* Output: \verbinclude MatrixBase_template_int_end.out
*
* \sa class Block
*/
template<typename Derived>
template<int Size>
inline VectorBlock<Derived,Size>
MatrixBase<Derived>::end()
{
EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
return VectorBlock<Derived, Size>(derived(), size() - Size);
}
/** This is the const version of end<int>.*/
template<typename Derived>
template<int Size>
inline const VectorBlock<Derived,Size>
MatrixBase<Derived>::end() const
{
EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
return VectorBlock<Derived, Size>(derived(), size() - Size);
}
#endif // EIGEN_VECTORBLOCK_H

2
Eigen/src/Core/util/Constants.h
View File

@@ -185,7 +185,7 @@ const unsigned int HereditaryBits = RowMajorBit
| EvalBeforeNestingBit
| EvalBeforeAssigningBit
| SparseBit;
// Possible values for the Mode parameter of part()
const unsigned int UpperTriangular = UpperTriangularBit;
const unsigned int StrictlyUpperTriangular = UpperTriangularBit | ZeroDiagBit;

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

@@ -42,6 +42,7 @@ template<typename MatrixType> class Minor;
template<typename MatrixType, int BlockRows=Dynamic, int BlockCols=Dynamic, int PacketAccess=AsRequested,
int _DirectAccessStatus = ei_traits<MatrixType>::Flags&DirectAccessBit ? DirectAccessBit
: ei_traits<MatrixType>::Flags&SparseBit> class Block;
template<typename MatrixType, int Size=Dynamic, int PacketAccess=AsRequested> class VectorBlock;
template<typename MatrixType> class Transpose;
template<typename MatrixType> class Conjugate;
template<typename NullaryOp, typename MatrixType> class CwiseNullaryOp;
@@ -111,7 +112,7 @@ template<typename MatrixType, int Direction = BothDirections> class Reverse;
template<typename MatrixType> class LU;
template<typename MatrixType> class PartialLU;
template<typename MatrixType> class QR;
template<typename MatrixType> class HouseholderQR;
template<typename MatrixType> class SVD;
template<typename MatrixType, int UpLo = LowerTriangular> class LLT;
template<typename MatrixType> class LDLT;

16
Eigen/src/Core/util/Macros.h
View File

@@ -51,7 +51,8 @@
#define EIGEN_GCC3_OR_OLDER 0
#endif
#if !EIGEN_GCC_AND_ARCH_DOESNT_WANT_ALIGNMENT && !EIGEN_GCC3_OR_OLDER
// FIXME vectorization + alignment is completely disabled with sun studio
#if !EIGEN_GCC_AND_ARCH_DOESNT_WANT_ALIGNMENT && !EIGEN_GCC3_OR_OLDER && !defined(__SUNPRO_CC)
#define EIGEN_ARCH_WANTS_ALIGNMENT 1
#else
#define EIGEN_ARCH_WANTS_ALIGNMENT 0
@@ -104,7 +105,7 @@
/** Allows to disable some optimizations which might affect the accuracy of the result.
* Such optimization are enabled by default, and set EIGEN_FAST_MATH to 0 to disable them.
* They currently include:
* - single precision Cwise::sin() and Cwise::cos() when SSE vectorization is enabled.
* - single precision Cwise::sin() and Cwise::cos() when SSE vectorization is enabled.
*/
#ifndef EIGEN_FAST_MATH
#define EIGEN_FAST_MATH 1
@@ -206,13 +207,16 @@ using Eigen::ei_cos;
* vectorized and non-vectorized code.
*/
#if !EIGEN_ALIGN
#define EIGEN_ALIGN_128
#define EIGEN_ALIGN_128
#elif (defined __GNUC__)
#define EIGEN_ALIGN_128 __attribute__((aligned(16)))
#define EIGEN_ALIGN_128 __attribute__((aligned(16)))
#elif (defined _MSC_VER)
#define EIGEN_ALIGN_128 __declspec(align(16))
#define EIGEN_ALIGN_128 __declspec(align(16))
#elif (defined __SUNPRO_CC)
// FIXME not sure about this one:
#define EIGEN_ALIGN_128 __attribute__((aligned(16)))
#else
#error Please tell me what is the equivalent of __attribute__((aligned(16))) for your compiler
#error Please tell me what is the equivalent of __attribute__((aligned(16))) for your compiler
#endif
#define EIGEN_RESTRICT __restrict

12
Eigen/src/Core/util/Memory.h
View File

@@ -27,7 +27,17 @@
#ifndef EIGEN_MEMORY_H
#define EIGEN_MEMORY_H
#if defined(__APPLE__) || defined(_WIN64) || defined (__FreeBSD__)
// FreeBSD 6 seems to have 16-byte aligned malloc
// See http://svn.freebsd.org/viewvc/base/stable/6/lib/libc/stdlib/malloc.c?view=markup
// FreeBSD 7 seems to have 16-byte aligned malloc except on ARM and MIPS architectures
// See http://svn.freebsd.org/viewvc/base/stable/7/lib/libc/stdlib/malloc.c?view=markup
#if defined(__FreeBSD__) && !defined(__arm__) && !defined(__mips__)
#define EIGEN_FREEBSD_MALLOC_ALREADY_ALIGNED 1
#else
#define EIGEN_FREEBSD_MALLOC_ALREADY_ALIGNED 0
#endif
#if defined(__APPLE__) || defined(_WIN64) || EIGEN_FREEBSD_MALLOC_ALREADY_ALIGNED
#define EIGEN_MALLOC_ALREADY_ALIGNED 1
#else
#define EIGEN_MALLOC_ALREADY_ALIGNED 0

6
Eigen/src/Core/util/XprHelper.h
View File

@@ -217,7 +217,7 @@ template<typename Derived,typename Scalar,typename OtherScalar,
bool EnableIt = !ei_is_same_type<Scalar,OtherScalar>::ret >
struct ei_special_scalar_op_base
{
// dummy operator* so that the
// dummy operator* so that the
// "using ei_special_scalar_op_base::operator*" compiles
void operator*() const;
};
@@ -237,8 +237,6 @@ struct ei_special_scalar_op_base<Derived,Scalar,OtherScalar,true>
* TODO: could be a good idea to define a big ReturnType struct ??
*/
template<typename ExpressionType, int RowsOrSize=Dynamic, int Cols=Dynamic> struct BlockReturnType {
typedef Block<ExpressionType, (ei_traits<ExpressionType>::RowsAtCompileTime == 1 ? 1 : RowsOrSize),
(ei_traits<ExpressionType>::ColsAtCompileTime == 1 ? 1 : RowsOrSize)> SubVectorType;
typedef Block<ExpressionType, RowsOrSize, Cols> Type;
};
@@ -251,7 +249,7 @@ template<typename ExpressionType> struct HNormalizedReturnType {
typedef Block<ExpressionType,
ei_traits<ExpressionType>::ColsAtCompileTime==1 ? SizeMinusOne : 1,
ei_traits<ExpressionType>::ColsAtCompileTime==1 ? 1 : SizeMinusOne> StartMinusOne;
typedef CwiseUnaryOp<ei_scalar_quotient1_op<typename ei_traits<ExpressionType>::Scalar>,
typedef CwiseUnaryOp<ei_scalar_quotient1_op<typename ei_traits<ExpressionType>::Scalar>,
NestByValue<StartMinusOne> > Type;
};

45
Eigen/src/Geometry/Quaternion.h
View File

@@ -303,7 +303,9 @@ inline Quaternion<Scalar>& Quaternion<Scalar>::operator=(const MatrixBase<Derive
return *this;
}
/** Convert the quaternion to a 3x3 rotation matrix */
/** Convert the quaternion to a 3x3 rotation matrix. The quaternion is required to
* be normalized, otherwise the result is undefined.
*/
template<typename Scalar>
inline typename Quaternion<Scalar>::Matrix3
Quaternion<Scalar>::toRotationMatrix(void) const
@@ -340,11 +342,15 @@ Quaternion<Scalar>::toRotationMatrix(void) const
return res;
}
/** Sets *this to be a quaternion representing a rotation sending the vector \a a to the vector \a b.
/** Sets \c *this to be a quaternion representing a rotation between
* the two arbitrary vectors \a a and \a b. In other words, the built
* rotation represent a rotation sending the line of direction \a a
* to the line of direction \a b, both lines passing through the origin.
*
* \returns a reference to *this.
* \returns a reference to \c *this.
*
* Note that the two input vectors do \b not have to be normalized.
* Note that the two input vectors do \b not have to be normalized, and
* do not need to have the same norm.
*/
template<typename Scalar>
template<typename Derived1, typename Derived2>
@@ -354,21 +360,26 @@ inline Quaternion<Scalar>& Quaternion<Scalar>::setFromTwoVectors(const MatrixBas
Vector3 v1 = b.normalized();
Scalar c = v0.dot(v1);
// if dot == 1, vectors are the same
if (ei_isApprox(c,Scalar(1)))
// if dot == -1, vectors are nearly opposites
// => accuraletly compute the rotation axis by computing the
// intersection of the two planes. This is done by solving:
// x^T v0 = 0
// x^T v1 = 0
// under the constraint:
// ||x|| = 1
// which yields a singular value problem
if (c < Scalar(-1)+precision<Scalar>())
{
// set to identity
this->w() = 1; this->vec().setZero();
return *this;
}
// if dot == -1, vectors are opposites
if (ei_isApprox(c,Scalar(-1)))
{
this->vec() = v0.unitOrthogonal();
this->w() = 0;
return *this;
}
c = std::max<Scalar>(c,-1);
Matrix<Scalar,2,3> m; m << v0.transpose(), v1.transpose();
SVD<Matrix<Scalar,2,3> > svd(m);
Vector3 axis = svd.matrixV().col(2);
Scalar w2 = (Scalar(1)+c)*Scalar(0.5);
this->w() = ei_sqrt(w2);
this->vec() = axis * ei_sqrt(Scalar(1) - w2);
return *this;
}
Vector3 axis = v0.cross(v1);
Scalar s = ei_sqrt((Scalar(1)+c)*Scalar(2));
Scalar invs = Scalar(1)/s;

183
Eigen/src/LU/Inverse.h
View File

@@ -29,41 +29,45 @@
*** Part 1 : optimized implementations for fixed-size 2,3,4 cases ***
********************************************************************/
template<typename MatrixType>
void ei_compute_inverse_in_size2_case(const MatrixType& matrix, MatrixType* result)
template<typename XprType, typename MatrixType>
inline void ei_compute_inverse_size2_helper(
const XprType& matrix, const typename MatrixType::Scalar& invdet,
MatrixType* result)
{
typedef typename MatrixType::Scalar Scalar;
const Scalar invdet = Scalar(1) / matrix.determinant();
result->coeffRef(0,0) = matrix.coeff(1,1) * invdet;
result->coeffRef(1,0) = -matrix.coeff(1,0) * invdet;
result->coeffRef(0,1) = -matrix.coeff(0,1) * invdet;
result->coeffRef(1,1) = matrix.coeff(0,0) * invdet;
}
template<typename MatrixType>
inline void ei_compute_inverse_size2(const MatrixType& matrix, MatrixType* result)
{
typedef typename MatrixType::Scalar Scalar;
const Scalar invdet = Scalar(1) / matrix.determinant();
ei_compute_inverse_size2_helper( matrix, invdet, result );
}
template<typename XprType, typename MatrixType>
bool ei_compute_inverse_in_size2_case_with_check(const XprType& matrix, MatrixType* result)
bool ei_compute_inverse_size2_with_check(const XprType& matrix, MatrixType* result)
{
typedef typename MatrixType::Scalar Scalar;
const Scalar det = matrix.determinant();
if(ei_isMuchSmallerThan(det, matrix.cwise().abs().maxCoeff())) return false;
const Scalar invdet = Scalar(1) / det;
result->coeffRef(0,0) = matrix.coeff(1,1) * invdet;
result->coeffRef(1,0) = -matrix.coeff(1,0) * invdet;
result->coeffRef(0,1) = -matrix.coeff(0,1) * invdet;
result->coeffRef(1,1) = matrix.coeff(0,0) * invdet;
ei_compute_inverse_size2_helper( matrix, invdet, result );
return true;
}
template<typename MatrixType>
void ei_compute_inverse_in_size3_case(const MatrixType& matrix, MatrixType* result)
template<typename XprType, typename MatrixType>
void ei_compute_inverse_size3_helper(
const XprType& matrix,
const typename MatrixType::Scalar& invdet,
const typename MatrixType::Scalar& det_minor00,
const typename MatrixType::Scalar& det_minor10,
const typename MatrixType::Scalar& det_minor20,
MatrixType* result)
{
typedef typename MatrixType::Scalar Scalar;
const Scalar det_minor00 = matrix.minor(0,0).determinant();
const Scalar det_minor10 = matrix.minor(1,0).determinant();
const Scalar det_minor20 = matrix.minor(2,0).determinant();
const Scalar invdet = Scalar(1) / ( det_minor00 * matrix.coeff(0,0)
- det_minor10 * matrix.coeff(1,0)
+ det_minor20 * matrix.coeff(2,0) );
result->coeffRef(0, 0) = det_minor00 * invdet;
result->coeffRef(0, 1) = -det_minor10 * invdet;
result->coeffRef(0, 2) = det_minor20 * invdet;
@@ -75,8 +79,24 @@ void ei_compute_inverse_in_size3_case(const MatrixType& matrix, MatrixType* resu
result->coeffRef(2, 2) = matrix.minor(2,2).determinant() * invdet;
}
template<bool Check, typename XprType, typename MatrixType>
bool ei_compute_inverse_size3(const XprType& matrix, MatrixType* result)
{
typedef typename MatrixType::Scalar Scalar;
const Scalar det_minor00 = matrix.minor(0,0).determinant();
const Scalar det_minor10 = matrix.minor(1,0).determinant();
const Scalar det_minor20 = matrix.minor(2,0).determinant();
const Scalar det = ( det_minor00 * matrix.coeff(0,0)
- det_minor10 * matrix.coeff(1,0)
+ det_minor20 * matrix.coeff(2,0) );
if(Check) if(ei_isMuchSmallerThan(det, matrix.cwise().abs().maxCoeff())) return false;
const Scalar invdet = Scalar(1) / det;
ei_compute_inverse_size3_helper( matrix, invdet, det_minor00, det_minor10, det_minor20, result );
return true;
}
template<typename MatrixType>
bool ei_compute_inverse_in_size4_case_helper(const MatrixType& matrix, MatrixType* result)
bool ei_compute_inverse_size4_helper(const MatrixType& matrix, MatrixType* result)
{
/* Let's split M into four 2x2 blocks:
* (P Q)
@@ -94,7 +114,7 @@ bool ei_compute_inverse_in_size4_case_helper(const MatrixType& matrix, MatrixTyp
typedef Block<MatrixType,2,2> XprBlock22;
typedef typename MatrixBase<XprBlock22>::PlainMatrixType Block22;
Block22 P_inverse;
if(ei_compute_inverse_in_size2_case_with_check(matrix.template block<2,2>(0,0), &P_inverse))
if(ei_compute_inverse_size2_with_check(matrix.template block<2,2>(0,0), &P_inverse))
{
const Block22 Q = matrix.template block<2,2>(0,2);
const Block22 P_inverse_times_Q = P_inverse * Q;
@@ -104,7 +124,7 @@ bool ei_compute_inverse_in_size4_case_helper(const MatrixType& matrix, MatrixTyp
const XprBlock22 S = matrix.template block<2,2>(2,2);
const Block22 X = S - R_times_P_inverse_times_Q;
Block22 Y;
ei_compute_inverse_in_size2_case(X, &Y);
ei_compute_inverse_size2(X, &Y);
result->template block<2,2>(2,2) = Y;
result->template block<2,2>(2,0) = - Y * R_times_P_inverse;
const Block22 Z = P_inverse_times_Q * Y;
@@ -118,13 +138,13 @@ bool ei_compute_inverse_in_size4_case_helper(const MatrixType& matrix, MatrixTyp
}
}
template<typename MatrixType>
void ei_compute_inverse_in_size4_case(const MatrixType& matrix, MatrixType* result)
template<typename XprType, typename MatrixType>
bool ei_compute_inverse_size4_with_check(const XprType& matrix, MatrixType* result)
{
if(ei_compute_inverse_in_size4_case_helper(matrix, result))
if(ei_compute_inverse_size4_helper(matrix, result))
{
// good ! The topleft 2x2 block was invertible, so the 2x2 blocks approach is successful.
return;
return true;
}
else
{
@@ -134,17 +154,17 @@ void ei_compute_inverse_in_size4_case(const MatrixType& matrix, MatrixType* resu
MatrixType m(matrix);
m.row(0).swap(m.row(2));
m.row(1).swap(m.row(3));
if(ei_compute_inverse_in_size4_case_helper(m, result))
if(ei_compute_inverse_size4_helper(m, result))
{
// good, the topleft 2x2 block of m is invertible. Since m is different from matrix in that some
// rows were permuted, the actual inverse of matrix is derived from the inverse of m by permuting
// the corresponding columns.
result->col(0).swap(result->col(2));
result->col(1).swap(result->col(3));
return true;
}
else
{
// last possible case. Since matrix is assumed to be invertible, this last case has to work.
// first, undo the swaps previously made
m.row(0).swap(m.row(2));
m.row(1).swap(m.row(3));
@@ -154,13 +174,23 @@ void ei_compute_inverse_in_size4_case(const MatrixType& matrix, MatrixType* resu
// swap row 1 with the the row among 2 and 3 that has the biggest 2 first coeffs
int swap1with = ei_abs(m.coeff(2,0))+ei_abs(m.coeff(2,1))>ei_abs(m.coeff(3,0))+ei_abs(m.coeff(3,1)) ? 2 : 3;
m.row(1).swap(m.row(swap1with));
ei_compute_inverse_in_size4_case_helper(m, result);
result->col(1).swap(result->col(swap1with));
result->col(0).swap(result->col(swap0with));
if( ei_compute_inverse_size4_helper(m, result) )
{
result->col(1).swap(result->col(swap1with));
result->col(0).swap(result->col(swap0with));
return true;
}
else
{
// non-invertible matrix
return false;
}
}
}
}
/***********************************************
*** Part 2 : selector and MatrixBase methods ***
***********************************************/
@@ -189,7 +219,7 @@ struct ei_compute_inverse<MatrixType, 2>
{
static inline void run(const MatrixType& matrix, MatrixType* result)
{
ei_compute_inverse_in_size2_case(matrix, result);
ei_compute_inverse_size2(matrix, result);
}
};
@@ -198,7 +228,7 @@ struct ei_compute_inverse<MatrixType, 3>
{
static inline void run(const MatrixType& matrix, MatrixType* result)
{
ei_compute_inverse_in_size3_case(matrix, result);
ei_compute_inverse_size3<false, MatrixType, MatrixType>(matrix, result);
}
};
@@ -207,7 +237,7 @@ struct ei_compute_inverse<MatrixType, 4>
{
static inline void run(const MatrixType& matrix, MatrixType* result)
{
ei_compute_inverse_in_size4_case(matrix, result);
ei_compute_inverse_size4_with_check(matrix, result);
}
};
@@ -215,14 +245,15 @@ struct ei_compute_inverse<MatrixType, 4>
*
* Computes the matrix inverse of this matrix.
*
* \note This matrix must be invertible, otherwise the result is undefined.
* \note This matrix must be invertible, otherwise the result is undefined. If you need an invertibility check, use
* computeInverseWithCheck().
*
* \param result Pointer to the matrix in which to store the result.
*
* Example: \include MatrixBase_computeInverse.cpp
* Output: \verbinclude MatrixBase_computeInverse.out
*
* \sa inverse()
* \sa inverse(), computeInverseWithCheck()
*/
template<typename Derived>
inline void MatrixBase<Derived>::computeInverse(PlainMatrixType *result) const
@@ -236,7 +267,8 @@ inline void MatrixBase<Derived>::computeInverse(PlainMatrixType *result) const
*
* \returns the matrix inverse of this matrix.
*
* \note This matrix must be invertible, otherwise the result is undefined.
* \note This matrix must be invertible, otherwise the result is undefined. If you need an invertibility check, use
* computeInverseWithCheck().
*
* \note This method returns a matrix by value, which can be inefficient. To avoid that overhead,
* use computeInverse() instead.
@@ -244,7 +276,7 @@ inline void MatrixBase<Derived>::computeInverse(PlainMatrixType *result) const
* Example: \include MatrixBase_inverse.cpp
* Output: \verbinclude MatrixBase_inverse.out
*
* \sa computeInverse()
* \sa computeInverse(), computeInverseWithCheck()
*/
template<typename Derived>
inline const typename MatrixBase<Derived>::PlainMatrixType MatrixBase<Derived>::inverse() const
@@ -254,4 +286,81 @@ inline const typename MatrixBase<Derived>::PlainMatrixType MatrixBase<Derived>::
return result;
}
/********************************************
* Compute inverse with invertibility check *
*******************************************/
template<typename MatrixType, int Size = MatrixType::RowsAtCompileTime>
struct ei_compute_inverse_with_check
{
static inline bool run(const MatrixType& matrix, MatrixType* result)
{
typedef typename MatrixType::Scalar Scalar;
LU<MatrixType> lu( matrix );
if( !lu.isInvertible() ) return false;
lu.computeInverse(result);
return true;
}
};
template<typename MatrixType>
struct ei_compute_inverse_with_check<MatrixType, 1>
{
static inline bool run(const MatrixType& matrix, MatrixType* result)
{
if( 0 == result->coeffRef(0,0) ) return false;
typedef typename MatrixType::Scalar Scalar;
result->coeffRef(0,0) = Scalar(1) / matrix.coeff(0,0);
return true;
}
};
template<typename MatrixType>
struct ei_compute_inverse_with_check<MatrixType, 2>
{
static inline bool run(const MatrixType& matrix, MatrixType* result)
{
return ei_compute_inverse_size2_with_check(matrix, result);
}
};
template<typename MatrixType>
struct ei_compute_inverse_with_check<MatrixType, 3>
{
static inline bool run(const MatrixType& matrix, MatrixType* result)
{
return ei_compute_inverse_size3<true, MatrixType, MatrixType>(matrix, result);
}
};
template<typename MatrixType>
struct ei_compute_inverse_with_check<MatrixType, 4>
{
static inline bool run(const MatrixType& matrix, MatrixType* result)
{
return ei_compute_inverse_size4_with_check(matrix, result);
}
};
/** \lu_module
*
* Computation of matrix inverse, with invertibility check.
*
* \returns true if the matrix is invertible, false otherwise.
*
* \param result Pointer to the matrix in which to store the result.
*
* \sa inverse(), computeInverse()
*/
template<typename Derived>
inline bool MatrixBase<Derived>::computeInverseWithCheck(PlainMatrixType *result) const
{
ei_assert(rows() == cols());
EIGEN_STATIC_ASSERT(NumTraits<Scalar>::HasFloatingPoint,NUMERIC_TYPE_MUST_BE_FLOATING_POINT)
return ei_compute_inverse_with_check<PlainMatrixType>::run(eval(), result);
}
#endif // EIGEN_INVERSE_H

158
Eigen/src/QR/QR.h
View File

@@ -28,18 +28,20 @@
/** \ingroup QR_Module
* \nonstableyet
*
* \class QR
* \class HouseholderQR
*
* \brief QR decomposition of a matrix
* \brief Householder QR decomposition of a matrix
*
* \param MatrixType the type of the matrix of which we are computing the QR decomposition
*
* This class performs a QR decomposition using Householder transformations. The result is
* stored in a compact way compatible with LAPACK.
*
* Note that no pivoting is performed. This is \b not a rank-revealing decomposition.
*
* \sa MatrixBase::qr()
*/
template<typename MatrixType> class QR
template<typename MatrixType> class HouseholderQR
{
public:
@@ -53,88 +55,23 @@ template<typename MatrixType> class QR
* \brief Default Constructor.
*
* The default constructor is useful in cases in which the user intends to
* perform decompositions via QR::compute(const MatrixType&).
* perform decompositions via HouseholderQR::compute(const MatrixType&).
*/
QR() : m_qr(), m_hCoeffs(), m_isInitialized(false) {}
HouseholderQR() : m_qr(), m_hCoeffs(), m_isInitialized(false) {}
QR(const MatrixType& matrix)
HouseholderQR(const MatrixType& matrix)
: m_qr(matrix.rows(), matrix.cols()),
m_hCoeffs(matrix.cols()),
m_isInitialized(false)
{
compute(matrix);
}
/** \deprecated use isInjective()
* \returns whether or not the matrix is of full rank
*
* \note Since the rank is computed only once, i.e. the first time it is needed, this
* method almost does not perform any further computation.
*/
EIGEN_DEPRECATED bool isFullRank() const
{
ei_assert(m_isInitialized && "QR is not initialized.");
return rank() == m_qr.cols();
}
/** \returns the rank of the matrix of which *this is the QR decomposition.
*
* \note Since the rank is computed only once, i.e. the first time it is needed, this
* method almost does not perform any further computation.
*/
int rank() const;
/** \returns the dimension of the kernel of the matrix of which *this is the QR decomposition.
*
* \note Since the rank is computed only once, i.e. the first time it is needed, this
* method almost does not perform any further computation.
*/
inline int dimensionOfKernel() const
{
ei_assert(m_isInitialized && "QR is not initialized.");
return m_qr.cols() - rank();
}
/** \returns true if the matrix of which *this is the QR decomposition represents an injective
* linear map, i.e. has trivial kernel; false otherwise.
*
* \note Since the rank is computed only once, i.e. the first time it is needed, this
* method almost does not perform any further computation.
*/
inline bool isInjective() const
{
ei_assert(m_isInitialized && "QR is not initialized.");
return rank() == m_qr.cols();
}
/** \returns true if the matrix of which *this is the QR decomposition represents a surjective
* linear map; false otherwise.
*
* \note Since the rank is computed only once, i.e. the first time it is needed, this
* method almost does not perform any further computation.
*/
inline bool isSurjective() const
{
ei_assert(m_isInitialized && "QR is not initialized.");
return rank() == m_qr.rows();
}
/** \returns true if the matrix of which *this is the QR decomposition is invertible.
*
* \note Since the rank is computed only once, i.e. the first time it is needed, this
* method almost does not perform any further computation.
*/
inline bool isInvertible() const
{
ei_assert(m_isInitialized && "QR is not initialized.");
return isInjective() && isSurjective();
}
/** \returns a read-only expression of the matrix R of the actual the QR decomposition */
const TriangularView<NestByValue<MatrixRBlockType>, UpperTriangular>
matrixR(void) const
{
ei_assert(m_isInitialized && "QR is not initialized.");
ei_assert(m_isInitialized && "HouseholderQR is not initialized.");
int cols = m_qr.cols();
return MatrixRBlockType(m_qr, 0, 0, cols, cols).nestByValue().template part<UpperTriangular>();
}
@@ -148,58 +85,35 @@ template<typename MatrixType> class QR
* Resized if necessary, so that result->rows()==A.cols() and result->cols()==b.cols().
* If no solution exists, *result is left with undefined coefficients.
*
* \returns true if any solution exists, false if no solution exists.
*
* \note If there exist more than one solution, this method will arbitrarily choose one.
* If you need a complete analysis of the space of solutions, take the one solution obtained
* by this method and add to it elements of the kernel, as determined by kernel().
*
* \note The case where b is a matrix is not yet implemented. Also, this
* code is space inefficient.
*
* Example: \include QR_solve.cpp
* Output: \verbinclude QR_solve.out
*
* \sa MatrixBase::solveTriangular(), kernel(), computeKernel(), inverse(), computeInverse()
* Example: \include HouseholderQR_solve.cpp
* Output: \verbinclude HouseholderQR_solve.out
*/
template<typename OtherDerived, typename ResultType>
bool solve(const MatrixBase<OtherDerived>& b, ResultType *result) const;
void solve(const MatrixBase<OtherDerived>& b, ResultType *result) const;
MatrixType matrixQ(void) const;
/** \returns a reference to the matrix where the Householder QR decomposition is stored
* in a LAPACK-compatible way.
*/
const MatrixType& matrixQR() const { return m_qr; }
void compute(const MatrixType& matrix);
protected:
MatrixType m_qr;
VectorType m_hCoeffs;
mutable int m_rank;
mutable bool m_rankIsUptodate;
bool m_isInitialized;
};
/** \returns the rank of the matrix of which *this is the QR decomposition. */
template<typename MatrixType>
int QR<MatrixType>::rank() const
{
ei_assert(m_isInitialized && "QR is not initialized.");
if (!m_rankIsUptodate)
{
RealScalar maxCoeff = m_qr.diagonal().cwise().abs().maxCoeff();
int n = m_qr.cols();
m_rank = 0;
while(m_rank<n && !ei_isMuchSmallerThan(m_qr.diagonal().coeff(m_rank), maxCoeff))
++m_rank;
m_rankIsUptodate = true;
}
return m_rank;
}
#ifndef EIGEN_HIDE_HEAVY_CODE
template<typename MatrixType>
void QR<MatrixType>::compute(const MatrixType& matrix)
void HouseholderQR<MatrixType>::compute(const MatrixType& matrix)
{
m_rankIsUptodate = false;
m_qr = matrix;
m_hCoeffs.resize(matrix.cols());
@@ -262,12 +176,12 @@ void QR<MatrixType>::compute(const MatrixType& matrix)
template<typename MatrixType>
template<typename OtherDerived, typename ResultType>
bool QR<MatrixType>::solve(
void HouseholderQR<MatrixType>::solve(
const MatrixBase<OtherDerived>& b,
ResultType *result
) const
{
ei_assert(m_isInitialized && "QR is not initialized.");
ei_assert(m_isInitialized && "HouseholderQR is not initialized.");
const int rows = m_qr.rows();
ei_assert(b.rows() == rows);
result->resize(rows, b.cols());
@@ -276,27 +190,17 @@ bool QR<MatrixType>::solve(
// Q^T without explicitly forming matrixQ(). Investigate.
*result = matrixQ().transpose()*b;
if(!isSurjective())
{
// is result is in the image of R ?
RealScalar biggest_in_res = result->corner(TopLeft, m_rank, result->cols()).cwise().abs().maxCoeff();
for(int col = 0; col < result->cols(); ++col)
for(int row = m_rank; row < result->rows(); ++row)
if(!ei_isMuchSmallerThan(result->coeff(row,col), biggest_in_res))
return false;
}
m_qr.corner(TopLeft, m_rank, m_rank)
const int rank = std::min(result->rows(), result->cols());
m_qr.corner(TopLeft, rank, rank)
.template marked<UpperTriangular>()
.solveTriangularInPlace(result->corner(TopLeft, m_rank, result->cols()));
return true;
.solveTriangularInPlace(result->corner(TopLeft, rank, result->cols()));
}
/** \returns the matrix Q */
template<typename MatrixType>
MatrixType QR<MatrixType>::matrixQ() const
MatrixType HouseholderQR<MatrixType>::matrixQ() const
{
ei_assert(m_isInitialized && "QR is not initialized.");
ei_assert(m_isInitialized && "HouseholderQR is not initialized.");
// compute the product Q_0 Q_1 ... Q_n-1,
// where Q_k is the k-th Householder transformation I - h_k v_k v_k'
// and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...]
@@ -319,15 +223,15 @@ MatrixType QR<MatrixType>::matrixQ() const
#endif // EIGEN_HIDE_HEAVY_CODE
/** \return the QR decomposition of \c *this.
/** \return the Householder QR decomposition of \c *this.
*
* \sa class QR
* \sa class HouseholderQR
*/
template<typename Derived>
const QR<typename MatrixBase<Derived>::PlainMatrixType>
MatrixBase<Derived>::qr() const
const HouseholderQR<typename MatrixBase<Derived>::PlainMatrixType>
MatrixBase<Derived>::householderQr() const
{
return QR<PlainMatrixType>(eval());
return HouseholderQR<PlainMatrixType>(eval());
}

606
Eigen/src/SVD/SVD.h
View File

@@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
// Copyright (C) 2008-2009 Gael Guennebaud <g.gael@free.fr>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
@@ -34,9 +34,7 @@
*
* \param MatrixType the type of the matrix of which we are computing the SVD decomposition
*
* This class performs a standard SVD decomposition of a real matrix A of size \c M x \c N
* with \c M \>= \c N.
*
* This class performs a standard SVD decomposition of a real matrix A of size \c M x \c N.
*
* \sa MatrixBase::SVD()
*/
@@ -55,13 +53,13 @@ template<typename MatrixType> class SVD
typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> ColVector;
typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> RowVector;
typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MinSize> MatrixUType;
typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> MatrixUType;
typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime> MatrixVType;
typedef Matrix<Scalar, MinSize, 1> SingularValuesType;
typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> SingularValuesType;
public:
/**
/**
* \brief Default Constructor.
*
* The default constructor is useful in cases in which the user intends to
@@ -70,9 +68,9 @@ template<typename MatrixType> class SVD
SVD() : m_matU(), m_matV(), m_sigma(), m_isInitialized(false) {}
SVD(const MatrixType& matrix)
: m_matU(matrix.rows(), std::min(matrix.rows(), matrix.cols())),
: m_matU(matrix.rows(), matrix.rows()),
m_matV(matrix.cols(),matrix.cols()),
m_sigma(std::min(matrix.rows(),matrix.cols())),
m_sigma(matrix.cols()),
m_isInitialized(false)
{
compute(matrix);
@@ -81,22 +79,22 @@ template<typename MatrixType> class SVD
template<typename OtherDerived, typename ResultType>
bool solve(const MatrixBase<OtherDerived> &b, ResultType* result) const;
const MatrixUType& matrixU() const
{
ei_assert(m_isInitialized && "SVD is not initialized.");
return m_matU;
}
const SingularValuesType& singularValues() const
const MatrixUType& matrixU() const
{
ei_assert(m_isInitialized && "SVD is not initialized.");
return m_sigma;
return m_matU;
}
const MatrixVType& matrixV() const
const SingularValuesType& singularValues() const
{
ei_assert(m_isInitialized && "SVD is not initialized.");
return m_matV;
return m_sigma;
}
const MatrixVType& matrixV() const
{
ei_assert(m_isInitialized && "SVD is not initialized.");
return m_matV;
}
void compute(const MatrixType& matrix);
@@ -111,6 +109,23 @@ template<typename MatrixType> class SVD
template<typename ScalingType, typename RotationType>
void computeScalingRotation(ScalingType *positive, RotationType *unitary) const;
protected:
// Computes (a^2 + b^2)^(1/2) without destructive underflow or overflow.
inline static Scalar pythag(Scalar a, Scalar b)
{
Scalar abs_a = ei_abs(a);
Scalar abs_b = ei_abs(b);
if (abs_a > abs_b)
return abs_a*ei_sqrt(Scalar(1.0)+ei_abs2(abs_b/abs_a));
else
return (abs_b == Scalar(0.0) ? Scalar(0.0) : abs_b*ei_sqrt(Scalar(1.0)+ei_abs2(abs_a/abs_b)));
}
inline static Scalar sign(Scalar a, Scalar b)
{
return (b >= Scalar(0.0) ? ei_abs(a) : -ei_abs(a));
}
protected:
/** \internal */
MatrixUType m_matU;
@@ -123,380 +138,271 @@ template<typename MatrixType> class SVD
/** Computes / recomputes the SVD decomposition A = U S V^* of \a matrix
*
* \note this code has been adapted from JAMA (public domain)
* \note this code has been adapted from Numerical Recipes, third edition.
*/
template<typename MatrixType>
void SVD<MatrixType>::compute(const MatrixType& matrix)
{
const int m = matrix.rows();
const int n = matrix.cols();
const int nu = std::min(m,n);
m_matU.resize(m, nu);
m_matU.resize(m, m);
m_matU.setZero();
m_sigma.resize(std::min(m,n));
m_sigma.resize(n);
m_matV.resize(n,n);
RowVector e(n);
ColVector work(m);
MatrixType matA(matrix);
const bool wantu = true;
const bool wantv = true;
int i=0, j=0, k=0;
int max_iters = 30;
// Reduce A to bidiagonal form, storing the diagonal elements
// in s and the super-diagonal elements in e.
int nct = std::min(m-1,n);
int nrt = std::max(0,std::min(n-2,m));
for (k = 0; k < std::max(nct,nrt); ++k)
MatrixVType& V = m_matV;
MatrixType A = matrix;
SingularValuesType& W = m_sigma;
bool flag;
int i,its,j,jj,k,l,nm;
Scalar anorm, c, f, g, h, s, scale, x, y, z;
bool convergence = true;
Scalar eps = precision<Scalar>();
Matrix<Scalar,Dynamic,1> rv1(n);
g = scale = anorm = 0;
// Householder reduction to bidiagonal form.
for (i=0; i<n; i++)
{
if (k < nct)
l = i+2;
rv1[i] = scale*g;
g = s = scale = 0.0;
if (i < m)
{
// Compute the transformation for the k-th column and
// place the k-th diagonal in m_sigma[k].
m_sigma[k] = matA.col(k).end(m-k).norm();
if (m_sigma[k] != 0.0) // FIXME
scale = A.col(i).end(m-i).cwise().abs().sum();
if (scale != Scalar(0))
{
if (matA(k,k) < 0.0)
m_sigma[k] = -m_sigma[k];
matA.col(k).end(m-k) /= m_sigma[k];
matA(k,k) += 1.0;
}
m_sigma[k] = -m_sigma[k];
}
for (j = k+1; j < n; ++j)
{
if ((k < nct) && (m_sigma[k] != 0.0))
{
// Apply the transformation.
Scalar t = matA.col(k).end(m-k).dot(matA.col(j).end(m-k)); // FIXME dot product or cwise prod + .sum() ??
t = -t/matA(k,k);
matA.col(j).end(m-k) += t * matA.col(k).end(m-k);
}
// Place the k-th row of A into e for the
// subsequent calculation of the row transformation.
e[j] = matA(k,j);
}
// Place the transformation in U for subsequent back multiplication.
if (wantu & (k < nct))
m_matU.col(k).end(m-k) = matA.col(k).end(m-k);
if (k < nrt)
{
// Compute the k-th row transformation and place the
// k-th super-diagonal in e[k].
e[k] = e.end(n-k-1).norm();
if (e[k] != 0.0)
{
if (e[k+1] < 0.0)
e[k] = -e[k];
e.end(n-k-1) /= e[k];
e[k+1] += 1.0;
}
e[k] = -e[k];
if ((k+1 < m) & (e[k] != 0.0))
{
// Apply the transformation.
work.end(m-k-1) = matA.corner(BottomRight,m-k-1,n-k-1) * e.end(n-k-1);
for (j = k+1; j < n; ++j)
matA.col(j).end(m-k-1) += (-e[j]/e[k+1]) * work.end(m-k-1);
}
// Place the transformation in V for subsequent back multiplication.
if (wantv)
m_matV.col(k).end(n-k-1) = e.end(n-k-1);
}
}
// Set up the final bidiagonal matrix or order p.
int p = std::min(n,m+1);
if (nct < n)
m_sigma[nct] = matA(nct,nct);
if (m < p)
m_sigma[p-1] = 0.0;
if (nrt+1 < p)
e[nrt] = matA(nrt,p-1);
e[p-1] = 0.0;
// If required, generate U.
if (wantu)
{
for (j = nct; j < nu; ++j)
{
m_matU.col(j).setZero();
m_matU(j,j) = 1.0;
}
for (k = nct-1; k >= 0; k--)
{
if (m_sigma[k] != 0.0)
{
for (j = k+1; j < nu; ++j)
for (k=i; k<m; k++)
{
Scalar t = m_matU.col(k).end(m-k).dot(m_matU.col(j).end(m-k)); // FIXME is it really a dot product we want ?
t = -t/m_matU(k,k);
m_matU.col(j).end(m-k) += t * m_matU.col(k).end(m-k);
A(k, i) /= scale;
s += A(k, i)*A(k, i);
}
m_matU.col(k).end(m-k) = - m_matU.col(k).end(m-k);
m_matU(k,k) = Scalar(1) + m_matU(k,k);
if (k-1>0)
m_matU.col(k).start(k-1).setZero();
}
else
{
m_matU.col(k).setZero();
m_matU(k,k) = 1.0;
f = A(i, i);
g = -sign( ei_sqrt(s), f );
h = f*g - s;
A(i, i)=f-g;
for (j=l-1; j<n; j++)
{
s = A.col(i).end(m-i).dot(A.col(j).end(m-i));
f = s/h;
A.col(j).end(m-i) += f*A.col(i).end(m-i);
}
A.col(i).end(m-i) *= scale;
}
}
}
// If required, generate V.
if (wantv)
{
for (k = n-1; k >= 0; k--)
W[i] = scale * g;
g = s = scale = 0.0;
if (i+1 <= m && i+1 != n)
{
if ((k < nrt) & (e[k] != 0.0))
scale = A.row(i).end(n-l+1).cwise().abs().sum();
if (scale != Scalar(0))
{
for (j = k+1; j < nu; ++j)
for (k=l-1; k<n; k++)
{
Scalar t = m_matV.col(k).end(n-k-1).dot(m_matV.col(j).end(n-k-1)); // FIXME is it really a dot product we want ?
t = -t/m_matV(k+1,k);
m_matV.col(j).end(n-k-1) += t * m_matV.col(k).end(n-k-1);
A(i, k) /= scale;
s += A(i, k)*A(i, k);
}
f = A(i,l-1);
g = -sign(ei_sqrt(s),f);
h = f*g - s;
A(i,l-1) = f-g;
rv1.end(n-l+1) = A.row(i).end(n-l+1)/h;
for (j=l-1; j<m; j++)
{
s = A.row(j).end(n-l+1).dot(A.row(i).end(n-l+1));
A.row(j).end(n-l+1) += s*rv1.end(n-l+1).transpose();
}
A.row(i).end(n-l+1) *= scale;
}
}
anorm = std::max( anorm, (ei_abs(W[i])+ei_abs(rv1[i])) );
}
// Accumulation of right-hand transformations.
for (i=n-1; i>=0; i--)
{
//Accumulation of right-hand transformations.
if (i < n-1)
{
if (g != Scalar(0.0))
{
for (j=l; j<n; j++) //Double division to avoid possible underflow.
V(j, i) = (A(i, j)/A(i, l))/g;
for (j=l; j<n; j++)
{
s = A.row(i).end(n-l).dot(V.col(j).end(n-l));
V.col(j).end(n-l) += s * V.col(i).end(n-l);
}
}
m_matV.col(k).setZero();
m_matV(k,k) = 1.0;
V.row(i).end(n-l).setZero();
V.col(i).end(n-l).setZero();
}
V(i, i) = 1.0;
g = rv1[i];
l = i;
}
// Main iteration loop for the singular values.
int pp = p-1;
int iter = 0;
Scalar eps = ei_pow(Scalar(2),ei_is_same_type<Scalar,float>::ret ? Scalar(-23) : Scalar(-52));
while (p > 0)
// Accumulation of left-hand transformations.
for (i=std::min(m,n)-1; i>=0; i--)
{
int k=0;
int kase=0;
// Here is where a test for too many iterations would go.
// This section of the program inspects for
// negligible elements in the s and e arrays. On
// completion the variables kase and k are set as follows.
// kase = 1 if s(p) and e[k-1] are negligible and k<p
// kase = 2 if s(k) is negligible and k<p
// kase = 3 if e[k-1] is negligible, k<p, and
// s(k), ..., s(p) are not negligible (qr step).
// kase = 4 if e(p-1) is negligible (convergence).
for (k = p-2; k >= -1; --k)
l = i+1;
g = W[i];
if (n-l>0)
A.row(i).end(n-l).setZero();
if (g != Scalar(0.0))
{
if (k == -1)
break;
if (ei_abs(e[k]) <= eps*(ei_abs(m_sigma[k]) + ei_abs(m_sigma[k+1])))
g = Scalar(1.0)/g;
for (j=l; j<n; j++)
{
e[k] = 0.0;
break;
s = A.col(i).end(m-l).dot(A.col(j).end(m-l));
f = (s/A(i,i))*g;
A.col(j).end(m-i) += f * A.col(i).end(m-i);
}
}
if (k == p-2)
{
kase = 4;
A.col(i).end(m-i) *= g;
}
else
A.col(i).end(m-i).setZero();
++A(i,i);
}
// Diagonalization of the bidiagonal form: Loop over
// singular values, and over allowed iterations.
for (k=n-1; k>=0; k--)
{
for (its=0; its<max_iters; its++)
{
int ks;
for (ks = p-1; ks >= k; --ks)
flag = true;
for (l=k; l>=0; l--)
{
if (ks == k)
break;
Scalar t = (ks != p ? ei_abs(e[ks]) : Scalar(0)) + (ks != k+1 ? ei_abs(e[ks-1]) : Scalar(0));
if (ei_abs(m_sigma[ks]) <= eps*t)
// Test for splitting.
nm = l-1;
// Note that rv1[1] is always zero.
//if ((double)(ei_abs(rv1[l])+anorm) == anorm)
if (l==0 || ei_abs(rv1[l]) <= eps*anorm)
{
m_sigma[ks] = 0.0;
flag = false;
break;
}
//if ((double)(ei_abs(W[nm])+anorm) == anorm)
if (ei_abs(W[nm]) <= eps*anorm)
break;
}
if (ks == k)
if (flag)
{
kase = 3;
c = 0.0; //Cancellation of rv1[l], if l > 0.
s = 1.0;
for (i=l ;i<k+1; i++)
{
f = s*rv1[i];
rv1[i] = c*rv1[i];
//if ((double)(ei_abs(f)+anorm) == anorm)
if (ei_abs(f) <= eps*anorm)
break;
g = W[i];
h = pythag(f,g);
W[i] = h;
h = Scalar(1.0)/h;
c = g*h;
s = -f*h;
for (j=0; j<m; j++)
{
y = A(j,nm);
z = A(j,i);
A(j,nm) = y*c + z*s;
A(j,i) = z*c - y*s;
}
}
}
else if (ks == p-1)
z = W[k];
if (l == k) //Convergence.
{
kase = 1;
if (z < 0.0) { // Singular value is made nonnegative.
W[k] = -z;
V.col(k) = -V.col(k);
}
break;
}
else
if (its+1 == max_iters)
{
kase = 2;
k = ks;
convergence = false;
}
x = W[l]; // Shift from bottom 2-by-2 minor.
nm = k-1;
y = W[nm];
g = rv1[nm];
h = rv1[k];
f = ((y-z)*(y+z) + (g-h)*(g+h))/(Scalar(2.0)*h*y);
g = pythag(f,1.0);
f = ((x-z)*(x+z) + h*((y/(f+sign(g,f)))-h))/x;
c = s = 1.0;
//Next QR transformation:
for (j=l; j<=nm; j++)
{
i = j+1;
g = rv1[i];
y = W[i];
h = s*g;
g = c*g;
z = pythag(f,h);
rv1[j] = z;
c = f/z;
s = h/z;
f = x*c + g*s;
g = g*c - x*s;
h = y*s;
y *= c;
for (jj=0; jj<n; jj++)
{
x = V(jj,j);
z = V(jj,i);
V(jj,j) = x*c + z*s;
V(jj,i) = z*c - x*s;
}
z = pythag(f,h);
W[j] = z;
// Rotation can be arbitrary if z = 0.
if (z!=Scalar(0))
{
z = Scalar(1.0)/z;
c = f*z;
s = h*z;
}
f = c*g + s*y;
x = c*y - s*g;
for (jj=0; jj<m; jj++)
{
y = A(jj,j);
z = A(jj,i);
A(jj,j) = y*c + z*s;
A(jj,i) = z*c - y*s;
}
}
rv1[l] = 0.0;
rv1[k] = f;
W[k] = x;
}
}
// sort the singular values:
{
for (int i=0; i<n; i++)
{
int k;
W.end(n-i).minCoeff(&k);
if (k != i)
{
std::swap(W[k],W[i]);
A.col(i).swap(A.col(k));
V.col(i).swap(V.col(k));
}
}
++k;
// Perform the task indicated by kase.
switch (kase)
{
// Deflate negligible s(p).
case 1:
{
Scalar f(e[p-2]);
e[p-2] = 0.0;
for (j = p-2; j >= k; --j)
{
Scalar t(ei_hypot(m_sigma[j],f));
Scalar cs(m_sigma[j]/t);
Scalar sn(f/t);
m_sigma[j] = t;
if (j != k)
{
f = -sn*e[j-1];
e[j-1] = cs*e[j-1];
}
if (wantv)
{
for (i = 0; i < n; ++i)
{
t = cs*m_matV(i,j) + sn*m_matV(i,p-1);
m_matV(i,p-1) = -sn*m_matV(i,j) + cs*m_matV(i,p-1);
m_matV(i,j) = t;
}
}
}
}
break;
// Split at negligible s(k).
case 2:
{
Scalar f(e[k-1]);
e[k-1] = 0.0;
for (j = k; j < p; ++j)
{
Scalar t(ei_hypot(m_sigma[j],f));
Scalar cs( m_sigma[j]/t);
Scalar sn(f/t);
m_sigma[j] = t;
f = -sn*e[j];
e[j] = cs*e[j];
if (wantu)
{
for (i = 0; i < m; ++i)
{
t = cs*m_matU(i,j) + sn*m_matU(i,k-1);
m_matU(i,k-1) = -sn*m_matU(i,j) + cs*m_matU(i,k-1);
m_matU(i,j) = t;
}
}
}
}
break;
// Perform one qr step.
case 3:
{
// Calculate the shift.
Scalar scale = std::max(std::max(std::max(std::max(
ei_abs(m_sigma[p-1]),ei_abs(m_sigma[p-2])),ei_abs(e[p-2])),
ei_abs(m_sigma[k])),ei_abs(e[k]));
Scalar sp = m_sigma[p-1]/scale;
Scalar spm1 = m_sigma[p-2]/scale;
Scalar epm1 = e[p-2]/scale;
Scalar sk = m_sigma[k]/scale;
Scalar ek = e[k]/scale;
Scalar b = ((spm1 + sp)*(spm1 - sp) + epm1*epm1)/Scalar(2);
Scalar c = (sp*epm1)*(sp*epm1);
Scalar shift = 0.0;
if ((b != 0.0) || (c != 0.0))
{
shift = ei_sqrt(b*b + c);
if (b < 0.0)
shift = -shift;
shift = c/(b + shift);
}
Scalar f = (sk + sp)*(sk - sp) + shift;
Scalar g = sk*ek;
// Chase zeros.
for (j = k; j < p-1; ++j)
{
Scalar t = ei_hypot(f,g);
Scalar cs = f/t;
Scalar sn = g/t;
if (j != k)
e[j-1] = t;
f = cs*m_sigma[j] + sn*e[j];
e[j] = cs*e[j] - sn*m_sigma[j];
g = sn*m_sigma[j+1];
m_sigma[j+1] = cs*m_sigma[j+1];
if (wantv)
{
for (i = 0; i < n; ++i)
{
t = cs*m_matV(i,j) + sn*m_matV(i,j+1);
m_matV(i,j+1) = -sn*m_matV(i,j) + cs*m_matV(i,j+1);
m_matV(i,j) = t;
}
}
t = ei_hypot(f,g);
cs = f/t;
sn = g/t;
m_sigma[j] = t;
f = cs*e[j] + sn*m_sigma[j+1];
m_sigma[j+1] = -sn*e[j] + cs*m_sigma[j+1];
g = sn*e[j+1];
e[j+1] = cs*e[j+1];
if (wantu && (j < m-1))
{
for (i = 0; i < m; ++i)
{
t = cs*m_matU(i,j) + sn*m_matU(i,j+1);
m_matU(i,j+1) = -sn*m_matU(i,j) + cs*m_matU(i,j+1);
m_matU(i,j) = t;
}
}
}
e[p-2] = f;
iter = iter + 1;
}
break;
// Convergence.
case 4:
{
// Make the singular values positive.
if (m_sigma[k] <= 0.0)
{
m_sigma[k] = m_sigma[k] < Scalar(0) ? -m_sigma[k] : Scalar(0);
if (wantv)
m_matV.col(k).start(pp+1) = -m_matV.col(k).start(pp+1);
}
// Order the singular values.
while (k < pp)
{
if (m_sigma[k] >= m_sigma[k+1])
break;
Scalar t = m_sigma[k];
m_sigma[k] = m_sigma[k+1];
m_sigma[k+1] = t;
if (wantv && (k < n-1))
m_matV.col(k).swap(m_matV.col(k+1));
if (wantu && (k < m-1))
m_matU.col(k).swap(m_matU.col(k+1));
++k;
}
iter = 0;
p--;
}
break;
} // end big switch
} // end iterations
}
m_matU.setZero();
if (m>=n)
m_matU.block(0,0,m,n) = A;
else
m_matU = A.block(0,0,m,m);
m_isInitialized = true;
}
@@ -554,6 +460,8 @@ bool SVD<MatrixType>::solve(const MatrixBase<OtherDerived> &b, ResultType* resul
const int rows = m_matU.rows();
ei_assert(b.rows() == rows);
result->resize(m_matV.rows(), b.cols());
Scalar maxVal = m_sigma.cwise().abs().maxCoeff();
for (int j=0; j<b.cols(); ++j)
{

5
Eigen/src/Sparse/SparseMatrix.h
View File

@@ -443,9 +443,8 @@ class SparseMatrix
// two passes algorithm:
// 1 - compute the number of coeffs per dest inner vector
// 2 - do the actual copy/eval
// Since each coeff of the rhs has to be evaluated twice, let's evauluate it if needed
//typedef typename ei_nested<OtherDerived,2>::type OtherCopy;
typedef typename ei_eval<OtherDerived>::type OtherCopy;
// Since each coeff of the rhs has to be evaluated twice, let's evaluate it if needed
typedef typename ei_nested<OtherDerived,2>::type OtherCopy;
typedef typename ei_cleantype<OtherCopy>::type _OtherCopy;
OtherCopy otherCopy(other.derived());

10
Eigen/src/Sparse/SparseMatrixBase.h
View File

@@ -102,8 +102,10 @@ template<typename Derived> class SparseMatrixBase
/** \internal the return type of MatrixBase::imag() */
typedef SparseCwiseUnaryOp<ei_scalar_imag_op<Scalar>, Derived> ImagReturnType;
/** \internal the return type of MatrixBase::adjoint() */
typedef SparseTranspose</*NestByValue<*/typename ei_cleantype<ConjugateReturnType>::type> /*>*/
AdjointReturnType;
typedef typename ei_meta_if<NumTraits<Scalar>::IsComplex,
SparseCwiseUnaryOp<ei_scalar_conjugate_op<Scalar>, SparseNestByValue<Eigen::SparseTranspose<Derived> > >,
SparseTranspose<Derived>
>::ret AdjointReturnType;
#ifndef EIGEN_PARSED_BY_DOXYGEN
/** This is the "real scalar" type; if the \a Scalar type is already real numbers
@@ -357,7 +359,7 @@ template<typename Derived> class SparseMatrixBase
SparseTranspose<Derived> transpose() { return derived(); }
const SparseTranspose<Derived> transpose() const { return derived(); }
// void transposeInPlace();
const AdjointReturnType adjoint() const { return conjugate()/*.nestByValue()*/; }
const AdjointReturnType adjoint() const { return transpose().nestByValue(); }
// sub-vector
SparseInnerVectorSet<Derived,1> row(int i);
@@ -529,7 +531,7 @@ template<typename Derived> class SparseMatrixBase
*/
// inline int stride(void) const { return derived().stride(); }
// inline const NestByValue<Derived> nestByValue() const;
inline const SparseNestByValue<Derived> nestByValue() const;
ConjugateReturnType conjugate() const;

84
Eigen/src/Sparse/SparseNestByValue.h Normal file
View File

@@ -0,0 +1,84 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008-2009 Gael Guennebaud <g.gael@free.fr>
// Copyright (C) 2006-2008 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/>.
#ifndef EIGEN_SPARSENESTBYVALUE_H
#define EIGEN_SPARSENESTBYVALUE_H
/** \class SparseNestByValue
*
* \brief Expression which must be nested by value
*
* \param ExpressionType the type of the object of which we are requiring nesting-by-value
*
* This class is the return type of MatrixBase::nestByValue()
* and most of the time this is the only way it is used.
*
* \sa SparseMatrixBase::nestByValue(), class NestByValue
*/
template<typename ExpressionType>
struct ei_traits<SparseNestByValue<ExpressionType> > : public ei_traits<ExpressionType>
{};
template<typename ExpressionType> class SparseNestByValue
: public SparseMatrixBase<SparseNestByValue<ExpressionType> >
{
public:
typedef typename ExpressionType::InnerIterator InnerIterator;
EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseNestByValue)
inline SparseNestByValue(const ExpressionType& matrix) : m_expression(matrix) {}
EIGEN_STRONG_INLINE int rows() const { return m_expression.rows(); }
EIGEN_STRONG_INLINE int cols() const { return m_expression.cols(); }
operator const ExpressionType&() const { return m_expression; }
protected:
const ExpressionType m_expression;
};
/** \returns an expression of the temporary version of *this.
*/
template<typename Derived>
inline const SparseNestByValue<Derived>
SparseMatrixBase<Derived>::nestByValue() const
{
return SparseNestByValue<Derived>(derived());
}
// template<typename MatrixType>
// class SparseNestByValue<MatrixType>::InnerIterator : public MatrixType::InnerIterator
// {
// typedef typename MatrixType::InnerIterator Base;
// public:
//
// EIGEN_STRONG_INLINE InnerIterator(const SparseNestByValue& expr, int outer)
// : Base(expr.m_expression, outer)
// {}
// };
#endif // EIGEN_SPARSENESTBYVALUE_H

3
Eigen/src/Sparse/SparseUtil.h
View File

@@ -106,6 +106,7 @@ template<typename _Scalar, int _Flags = 0> class DynamicSparseMatrix;
template<typename _Scalar, int _Flags = 0> class SparseVector;
template<typename _Scalar, int _Flags = 0> class MappedSparseMatrix;
template<typename MatrixType> class SparseNestByValue;
template<typename MatrixType> class SparseTranspose;
template<typename MatrixType, int Size> class SparseInnerVectorSet;
template<typename Derived> class SparseCwise;
@@ -146,4 +147,6 @@ template<typename T> class ei_eval<T,IsSparse>
typedef SparseMatrix<_Scalar, _Flags> type;
};
template<typename T> struct ei_must_nest_by_value<SparseNestByValue<T> > { enum { ret = true }; };
#endif // EIGEN_SPARSEUTIL_H