GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-08-13 20:26:03 +08:00
Code Issues Packages Projects Releases Wiki Activity

merge with default branch

Browse Source
This commit is contained in:
Gael Guennebaud 2014-09-14 17:34:54 +02:00
commit 749b56f6af
34 changed files with 426 additions and 111 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
Eigen/src/Core/ArrayWrapper.h
View File

@ -29,6 +29,11 @@ struct traits<ArrayWrapper<ExpressionType> >
: public traits<typename remove_all<typename ExpressionType::Nested>::type >
{
typedef ArrayXpr XprKind;
// Let's remove NestByRefBit
enum {
Flags0 = traits<typename remove_all<typename ExpressionType::Nested>::type >::Flags,
Flags = Flags0 & ~NestByRefBit
};
};
}
@ -167,6 +172,11 @@ struct traits<MatrixWrapper<ExpressionType> >
: public traits<typename remove_all<typename ExpressionType::Nested>::type >
{
typedef MatrixXpr XprKind;
// Let's remove NestByRefBit
enum {
Flags0 = traits<typename remove_all<typename ExpressionType::Nested>::type >::Flags,
Flags = Flags0 & ~NestByRefBit
};
};
}

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

@ -264,7 +264,7 @@ EIGEN_DONT_INLINE void outer_product_selector_run(const ProductType& prod, Dest&
// FIXME not very good if rhs is real and lhs complex while alpha is real too
const Index cols = dest.cols();
for (Index j=0; j<cols; ++j)
func(dest.col(j), prod.rhs().coeff(j) * prod.lhs());
func(dest.col(j), prod.rhs().coeff(0,j) * prod.lhs());
}
// Row major
@ -275,7 +275,7 @@ EIGEN_DONT_INLINE void outer_product_selector_run(const ProductType& prod, Dest&
// FIXME not very good if lhs is real and rhs complex while alpha is real too
const Index rows = dest.rows();
for (Index i=0; i<rows; ++i)
func(dest.row(i), prod.lhs().coeff(i) * prod.rhs());
func(dest.row(i), prod.lhs().coeff(i,0) * prod.rhs());
}
template<typename Lhs, typename Rhs>

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

@ -12,6 +12,15 @@
namespace Eigen {
// On WINCE, std::abs is defined for int only, so let's defined our own overloads:
// This issue has been confirmed with MSVC 2008 only, but the issue might exist for more recent versions too.
#if defined(_WIN32_WCE) && defined(_MSC_VER) && _MSC_VER<=1500
long abs(long x) { return (labs(x)); }
double abs(double x) { return (fabs(x)); }
float abs(float x) { return (fabsf(x)); }
long double abs(long double x) { return (fabsl(x)); }
#endif
namespace internal {
/** \internal \struct global_math_functions_filtering_base
@ -308,10 +317,17 @@ struct hypot_impl
using std::sqrt;
RealScalar _x = abs(x);
RealScalar _y = abs(y);
RealScalar p = (max)(_x, _y);
if(p==RealScalar(0)) return 0;
RealScalar q = (min)(_x, _y);
RealScalar qp = q/p;
Scalar p, qp;
if(_x>_y)
{
p = _x;
qp = _y / p;
}
else
{
p = _y;
qp = _x / p;
}
return p * sqrt(RealScalar(1) + qp*qp);
}
};

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

@ -83,6 +83,7 @@ template<typename Derived> class MatrixBase
using Base::operator-=;
using Base::operator*=;
using Base::operator/=;
using Base::operator*;
typedef typename Base::CoeffReturnType CoeffReturnType;
typedef typename Base::ConstTransposeReturnType ConstTransposeReturnType;

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

@ -230,7 +230,7 @@ struct redux_impl<Func, Derived, LinearVectorizedTraversal, NoUnrolling>
const Index packetSize = packet_traits<Scalar>::size;
const Index alignedStart = internal::first_aligned(mat);
enum {
alignment = bool(Derived::Flags & DirectAccessBit) || bool(Derived::Flags & AlignedBit)
alignment = (bool(Derived::Flags & DirectAccessBit) && bool(packet_traits<Scalar>::AlignedOnScalar)) || bool(Derived::Flags & AlignedBit)
? Aligned : Unaligned
};
const Index alignedSize2 = ((size-alignedStart)/(2*packetSize))*(2*packetSize);

10
Eigen/src/Core/SelfCwiseBinaryOp.h
View File

@ -256,15 +256,9 @@ inline Derived& ArrayBase<Derived>::operator-=(const Scalar& other)
template<typename Derived>
inline Derived& DenseBase<Derived>::operator/=(const Scalar& other)
{
typedef typename internal::conditional<NumTraits<Scalar>::IsInteger,
internal::scalar_quotient_op<Scalar>,
internal::scalar_product_op<Scalar> >::type BinOp;
typedef typename Derived::PlainObject PlainObject;
SelfCwiseBinaryOp<BinOp, Derived, typename PlainObject::ConstantReturnType> tmp(derived());
Scalar actual_other;
if(NumTraits<Scalar>::IsInteger) actual_other = other;
else actual_other = Scalar(1)/other;
tmp = PlainObject::Constant(rows(),cols(), actual_other);
SelfCwiseBinaryOp<internal::scalar_quotient_op<Scalar>, Derived, typename PlainObject::ConstantReturnType> tmp(derived());
tmp = PlainObject::Constant(rows(),cols(), other);
return derived();
}
#endif

22
Eigen/src/Core/StableNorm.h
View File

@ -20,7 +20,7 @@ inline void stable_norm_kernel(const ExpressionType& bl, Scalar& ssq, Scalar& sc
using std::max;
Scalar maxCoeff = bl.cwiseAbs().maxCoeff();
if (maxCoeff>scale)
if(maxCoeff>scale)
{
ssq = ssq * numext::abs2(scale/maxCoeff);
Scalar tmp = Scalar(1)/maxCoeff;
@ -29,12 +29,21 @@ inline void stable_norm_kernel(const ExpressionType& bl, Scalar& ssq, Scalar& sc
invScale = NumTraits<Scalar>::highest();
scale = Scalar(1)/invScale;
}
else if(maxCoeff>NumTraits<Scalar>::highest()) // we got a INF
{
invScale = Scalar(1);
scale = maxCoeff;
}
else
{
scale = maxCoeff;
invScale = tmp;
}
}
else if(maxCoeff!=maxCoeff) // we got a NaN
{
scale = maxCoeff;
}
// TODO if the maxCoeff is much much smaller than the current scale,
// then we can neglect this sub vector
@ -55,7 +64,7 @@ blueNorm_impl(const EigenBase<Derived>& _vec)
using std::abs;
const Derived& vec(_vec.derived());
static bool initialized = false;
static RealScalar b1, b2, s1m, s2m, overfl, rbig, relerr;
static RealScalar b1, b2, s1m, s2m, rbig, relerr;
if(!initialized)
{
int ibeta, it, iemin, iemax, iexp;
@ -84,7 +93,6 @@ blueNorm_impl(const EigenBase<Derived>& _vec)
iexp = - ((iemax+it)/2);
s2m = RealScalar(pow(RealScalar(ibeta),RealScalar(iexp))); // scaling factor for upper range
overfl = rbig*s2m; // overflow boundary for abig
eps = RealScalar(pow(double(ibeta), 1-it));
relerr = sqrt(eps); // tolerance for neglecting asml
initialized = true;
@ -101,13 +109,13 @@ blueNorm_impl(const EigenBase<Derived>& _vec)
else if(ax < b1) asml += numext::abs2(ax*s1m);
else amed += numext::abs2(ax);
}
if(amed!=amed)
return amed; // we got a NaN
if(abig > RealScalar(0))
{
abig = sqrt(abig);
if(abig > overfl)
{
return rbig;
}
if(abig > rbig) // overflow, or *this contains INF values
return abig; // return INF
if(amed > RealScalar(0))
{
abig = abig/s2m;

4
Eigen/src/Core/arch/AVX/PacketMath.h
View File

@ -141,7 +141,7 @@ template<> EIGEN_STRONG_INLINE Packet8f pmadd(const Packet8f& a, const Packet8f&
// so let's enforce it to generate a vfmadd231ps instruction since the most common use case is to accumulate
// the result of the product.
Packet8f res = c;
asm("vfmadd231ps %[a], %[b], %[c]" : [c] "+x" (res) : [a] "x" (a), [b] "x" (b));
__asm__("vfmadd231ps %[a], %[b], %[c]" : [c] "+x" (res) : [a] "x" (a), [b] "x" (b));
return res;
#else
return _mm256_fmadd_ps(a,b,c);
@ -151,7 +151,7 @@ template<> EIGEN_STRONG_INLINE Packet4d pmadd(const Packet4d& a, const Packet4d&
#if defined(__clang__) || defined(__GNUC__)
// see above
Packet4d res = c;
asm("vfmadd231pd %[a], %[b], %[c]" : [c] "+x" (res) : [a] "x" (a), [b] "x" (b));
__asm__("vfmadd231pd %[a], %[b], %[c]" : [c] "+x" (res) : [a] "x" (a), [b] "x" (b));
return res;
#else
return _mm256_fmadd_pd(a,b,c);

2
Eigen/src/Core/arch/NEON/PacketMath.h
View File

@ -57,7 +57,7 @@ typedef uint32x4_t Packet4ui;
#elif defined __pld
#define EIGEN_ARM_PREFETCH(ADDR) __pld(ADDR)
#elif !defined(__aarch64__)
#define EIGEN_ARM_PREFETCH(ADDR) asm volatile ( " pld [%[addr]]\n" :: [addr] "r" (ADDR) : "cc" );
#define EIGEN_ARM_PREFETCH(ADDR) __asm__ __volatile__ ( " pld [%[addr]]\n" :: [addr] "r" (ADDR) : "cc" );
#else
// by default no explicit prefetching
#define EIGEN_ARM_PREFETCH(ADDR)

14
Eigen/src/Core/functors/BinaryFunctors.h
View File

@ -167,9 +167,17 @@ template<typename Scalar> struct scalar_hypot_op {
EIGEN_USING_STD_MATH(max);
EIGEN_USING_STD_MATH(min);
using std::sqrt;
Scalar p = (max)(_x, _y);
Scalar q = (min)(_x, _y);
Scalar qp = q/p;
Scalar p, qp;
if(_x>_y)
{
p = _x;
qp = _y / p;
}
else
{
p = _y;
qp = _x / p;
}
return p * sqrt(Scalar(1) + qp*qp);
}
};

2
Eigen/src/Core/products/GeneralMatrixMatrix_MKL.h
View File

@ -53,6 +53,8 @@ template< \
int RhsStorageOrder, bool ConjugateRhs> \
struct general_matrix_matrix_product<Index,EIGTYPE,LhsStorageOrder,ConjugateLhs,EIGTYPE,RhsStorageOrder,ConjugateRhs,ColMajor> \
{ \
typedef gebp_traits<EIGTYPE,EIGTYPE> Traits; \
\
static void run(Index rows, Index cols, Index depth, \
const EIGTYPE* _lhs, Index lhsStride, \
const EIGTYPE* _rhs, Index rhsStride, \

4
Eigen/src/Core/products/TriangularMatrixMatrix_MKL.h
View File

@ -109,7 +109,7 @@ struct product_triangular_matrix_matrix_trmm<EIGTYPE,Index,Mode,true, \
/* Non-square case - doesn't fit to MKL ?TRMM. Fall to default triangular product or call MKL ?GEMM*/ \
if (rows != depth) { \
\
int nthr = mkl_domain_get_max_threads(MKL_BLAS); \
int nthr = mkl_domain_get_max_threads(EIGEN_MKL_DOMAIN_BLAS); \
\
if (((nthr==1) && (((std::max)(rows,depth)-diagSize)/(double)diagSize < 0.5))) { \
/* Most likely no benefit to call TRMM or GEMM from MKL*/ \
@ -223,7 +223,7 @@ struct product_triangular_matrix_matrix_trmm<EIGTYPE,Index,Mode,false, \
/* Non-square case - doesn't fit to MKL ?TRMM. Fall to default triangular product or call MKL ?GEMM*/ \
if (cols != depth) { \
\
int nthr = mkl_domain_get_max_threads(MKL_BLAS); \
int nthr = mkl_domain_get_max_threads(EIGEN_MKL_DOMAIN_BLAS); \
\
if ((nthr==1) && (((std::max)(cols,depth)-diagSize)/(double)diagSize < 0.5)) { \
/* Most likely no benefit to call TRMM or GEMM from MKL*/ \

32
Eigen/src/Core/util/MKL_support.h
View File

@ -76,6 +76,38 @@
#include <mkl_lapacke.h>
#define EIGEN_MKL_VML_THRESHOLD 128
/* MKL_DOMAIN_BLAS, etc are defined only in 10.3 update 7 */
/* MKL_BLAS, etc are not defined in 11.2 */
#ifdef MKL_DOMAIN_ALL
#define EIGEN_MKL_DOMAIN_ALL MKL_DOMAIN_ALL
#else
#define EIGEN_MKL_DOMAIN_ALL MKL_ALL
#endif
#ifdef MKL_DOMAIN_BLAS
#define EIGEN_MKL_DOMAIN_BLAS MKL_DOMAIN_BLAS
#else
#define EIGEN_MKL_DOMAIN_BLAS MKL_BLAS
#endif
#ifdef MKL_DOMAIN_FFT
#define EIGEN_MKL_DOMAIN_FFT MKL_DOMAIN_FFT
#else
#define EIGEN_MKL_DOMAIN_FFT MKL_FFT
#endif
#ifdef MKL_DOMAIN_VML
#define EIGEN_MKL_DOMAIN_VML MKL_DOMAIN_VML
#else
#define EIGEN_MKL_DOMAIN_VML MKL_VML
#endif
#ifdef MKL_DOMAIN_PARDISO
#define EIGEN_MKL_DOMAIN_PARDISO MKL_DOMAIN_PARDISO
#else
#define EIGEN_MKL_DOMAIN_PARDISO MKL_PARDISO
#endif
namespace Eigen {
typedef std::complex<double> dcomplex;

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

@ -279,7 +279,7 @@ namespace Eigen {
#if !defined(EIGEN_ASM_COMMENT)
#if (defined __GNUC__) && ( defined(__i386__) || defined(__x86_64__) )
#define EIGEN_ASM_COMMENT(X) asm("#" X)
#define EIGEN_ASM_COMMENT(X) __asm__("#" X)
#else
#define EIGEN_ASM_COMMENT(X)
#endif

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

@ -506,13 +506,21 @@ struct special_scalar_op_base<Derived,Scalar,OtherScalar,true> : public DenseCo
const CwiseUnaryOp<scalar_multiple2_op<Scalar,OtherScalar>, Derived>
operator*(const OtherScalar& scalar) const
{
#ifdef EIGEN_SPECIAL_SCALAR_MULTIPLE_PLUGIN
EIGEN_SPECIAL_SCALAR_MULTIPLE_PLUGIN
#endif
return CwiseUnaryOp<scalar_multiple2_op<Scalar,OtherScalar>, Derived>
(*static_cast<const Derived*>(this), scalar_multiple2_op<Scalar,OtherScalar>(scalar));
}
inline friend const CwiseUnaryOp<scalar_multiple2_op<Scalar,OtherScalar>, Derived>
operator*(const OtherScalar& scalar, const Derived& matrix)
{ return static_cast<const special_scalar_op_base&>(matrix).operator*(scalar); }
{
#ifdef EIGEN_SPECIAL_SCALAR_MULTIPLE_PLUGIN
EIGEN_SPECIAL_SCALAR_MULTIPLE_PLUGIN
#endif
return static_cast<const special_scalar_op_base&>(matrix).operator*(scalar);
}
};
template<typename XprType, typename CastType> struct cast_return_type

4
Eigen/src/Geometry/Homogeneous.h
View File

@ -158,7 +158,7 @@ template<typename MatrixType,int _Direction> class Homogeneous
* Example: \include MatrixBase_homogeneous.cpp
* Output: \verbinclude MatrixBase_homogeneous.out
*
* \sa class Homogeneous
* \sa VectorwiseOp::homogeneous(), class Homogeneous
*/
template<typename Derived>
inline typename MatrixBase<Derived>::HomogeneousReturnType
@ -175,7 +175,7 @@ MatrixBase<Derived>::homogeneous() const
* Example: \include VectorwiseOp_homogeneous.cpp
* Output: \verbinclude VectorwiseOp_homogeneous.out
*
* \sa MatrixBase::homogeneous() */
* \sa MatrixBase::homogeneous(), class Homogeneous */
template<typename ExpressionType, int Direction>
inline Homogeneous<ExpressionType,Direction>
VectorwiseOp<ExpressionType,Direction>::homogeneous() const

20
Eigen/src/SVD/JacobiSVD.h
View File

@ -425,18 +425,19 @@ void real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q,
JacobiRotation<RealScalar> rot1;
RealScalar t = m.coeff(0,0) + m.coeff(1,1);
RealScalar d = m.coeff(1,0) - m.coeff(0,1);
if(t == RealScalar(0))
if(d == RealScalar(0))
{
rot1.c() = RealScalar(0);
rot1.s() = d > RealScalar(0) ? RealScalar(1) : RealScalar(-1);
rot1.s() = RealScalar(0);
rot1.c() = RealScalar(1);
}
else
{
RealScalar t2d2 = numext::hypot(t,d);
rot1.c() = abs(t)/t2d2;
rot1.s() = d/t2d2;
if(t<RealScalar(0))
rot1.s() = -rot1.s();
// If d!=0, then t/d cannot overflow because the magnitude of the
// entries forming d are not too small compared to the ones forming t.
RealScalar u = t / d;
rot1.s() = RealScalar(1) / sqrt(RealScalar(1) + numext::abs2(u));
rot1.c() = rot1.s() * u;
}
m.applyOnTheLeft(0,1,rot1);
j_right->makeJacobi(m,0,1);
@ -726,7 +727,8 @@ JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsig
EIGEN_USING_STD_MATH(max);
RealScalar threshold = (max)(considerAsZero, precision * (max)(abs(m_workMatrix.coeff(p,p)),
abs(m_workMatrix.coeff(q,q))));
if((max)(abs(m_workMatrix.coeff(p,q)),abs(m_workMatrix.coeff(q,p))) > threshold)
// We compare both values to threshold instead of calling max to be robust to NaN (See bug 791)
if(abs(m_workMatrix.coeff(p,q))>threshold || abs(m_workMatrix.coeff(q,p)) > threshold)
{
finished = false;

29
Eigen/src/SVD/UpperBidiagonalization.h
View File

@ -154,14 +154,19 @@ void upperbidiagonalization_blocked_helper(MatrixType& A,
typename MatrixType::RealScalar *diagonal,
typename MatrixType::RealScalar *upper_diagonal,
typename MatrixType::Index bs,
Ref<Matrix<typename MatrixType::Scalar, Dynamic, Dynamic> > X,
Ref<Matrix<typename MatrixType::Scalar, Dynamic, Dynamic> > Y)
Ref<Matrix<typename MatrixType::Scalar, Dynamic, Dynamic,
traits<MatrixType>::Flags & RowMajorBit> > X,
Ref<Matrix<typename MatrixType::Scalar, Dynamic, Dynamic,
traits<MatrixType>::Flags & RowMajorBit> > Y)
{
typedef typename MatrixType::Index Index;
typedef typename MatrixType::Scalar Scalar;
typedef Ref<Matrix<Scalar, Dynamic, 1> > SubColumnType;
typedef Ref<Matrix<Scalar, 1, Dynamic>, 0, InnerStride<> > SubRowType;
typedef Ref<Matrix<Scalar, Dynamic, Dynamic> > SubMatType;
enum { StorageOrder = traits<MatrixType>::Flags & RowMajorBit };
typedef InnerStride<int(StorageOrder) == int(ColMajor) ? 1 : Dynamic> ColInnerStride;
typedef InnerStride<int(StorageOrder) == int(ColMajor) ? Dynamic : 1> RowInnerStride;
typedef Ref<Matrix<Scalar, Dynamic, 1>, 0, ColInnerStride> SubColumnType;
typedef Ref<Matrix<Scalar, 1, Dynamic>, 0, RowInnerStride> SubRowType;
typedef Ref<Matrix<Scalar, Dynamic, Dynamic, StorageOrder > > SubMatType;
Index brows = A.rows();
Index bcols = A.cols();
@ -288,8 +293,18 @@ void upperbidiagonalization_inplace_blocked(MatrixType& A, BidiagType& bidiagona
Index cols = A.cols();
Index size = (std::min)(rows, cols);
Matrix<Scalar,MatrixType::RowsAtCompileTime,Dynamic,ColMajor,MatrixType::MaxRowsAtCompileTime> X(rows,maxBlockSize);
Matrix<Scalar,MatrixType::ColsAtCompileTime,Dynamic,ColMajor,MatrixType::MaxColsAtCompileTime> Y(cols,maxBlockSize);
// X and Y are work space
enum { StorageOrder = traits<MatrixType>::Flags & RowMajorBit };
Matrix<Scalar,
MatrixType::RowsAtCompileTime,
Dynamic,
StorageOrder,
MatrixType::MaxRowsAtCompileTime> X(rows,maxBlockSize);
Matrix<Scalar,
MatrixType::ColsAtCompileTime,
Dynamic,
StorageOrder,
MatrixType::MaxColsAtCompileTime> Y(cols,maxBlockSize);
Index blockSize = (std::min)(maxBlockSize,size);
Index k = 0;

18
Eigen/src/plugins/ArrayCwiseUnaryOps.h
View File

@ -29,6 +29,9 @@ abs2() const
}
/** \returns an expression of the coefficient-wise exponential of *this.
*
* This function computes the coefficient-wise exponential. The function MatrixBase::exp() in the
* unsupported module MatrixFunctions computes the matrix exponential.
*
* Example: \include Cwise_exp.cpp
* Output: \verbinclude Cwise_exp.out
@ -43,6 +46,9 @@ exp() const
}
/** \returns an expression of the coefficient-wise logarithm of *this.
*
* This function computes the coefficient-wise logarithm. The function MatrixBase::log() in the
* unsupported module MatrixFunctions computes the matrix logarithm.
*
* Example: \include Cwise_log.cpp
* Output: \verbinclude Cwise_log.out
@ -57,6 +63,9 @@ log() const
}
/** \returns an expression of the coefficient-wise square root of *this.
*
* This function computes the coefficient-wise square root. The function MatrixBase::sqrt() in the
* unsupported module MatrixFunctions computes the matrix square root.
*
* Example: \include Cwise_sqrt.cpp
* Output: \verbinclude Cwise_sqrt.out
@ -71,6 +80,9 @@ sqrt() const
}
/** \returns an expression of the coefficient-wise cosine of *this.
*
* This function computes the coefficient-wise cosine. The function MatrixBase::cos() in the
* unsupported module MatrixFunctions computes the matrix cosine.
*
* Example: \include Cwise_cos.cpp
* Output: \verbinclude Cwise_cos.out
@ -86,6 +98,9 @@ cos() const
/** \returns an expression of the coefficient-wise sine of *this.
*
* This function computes the coefficient-wise sine. The function MatrixBase::sin() in the
* unsupported module MatrixFunctions computes the matrix sine.
*
* Example: \include Cwise_sin.cpp
* Output: \verbinclude Cwise_sin.out
@ -155,6 +170,9 @@ atan() const
}
/** \returns an expression of the coefficient-wise power of *this to the given exponent.
*
* This function computes the coefficient-wise power. The function MatrixBase::pow() in the
* unsupported module MatrixFunctions computes the matrix power.
*
* Example: \include Cwise_pow.cpp
* Output: \verbinclude Cwise_pow.out

33
bench/bench_norm.cpp
View File

@ -6,19 +6,25 @@ using namespace Eigen;
using namespace std;
template<typename T>
EIGEN_DONT_INLINE typename T::Scalar sqsumNorm(const T& v)
EIGEN_DONT_INLINE typename T::Scalar sqsumNorm(T& v)
{
return v.norm();
}
template<typename T>
EIGEN_DONT_INLINE typename T::Scalar hypotNorm(const T& v)
EIGEN_DONT_INLINE typename T::Scalar stableNorm(T& v)
{
return v.stableNorm();
}
template<typename T>
EIGEN_DONT_INLINE typename T::Scalar hypotNorm(T& v)
{
return v.hypotNorm();
}
template<typename T>
EIGEN_DONT_INLINE typename T::Scalar blueNorm(const T& v)
EIGEN_DONT_INLINE typename T::Scalar blueNorm(T& v)
{
return v.blueNorm();
}
@ -217,20 +223,21 @@ EIGEN_DONT_INLINE typename T::Scalar pblueNorm(const T& v)
}
#define BENCH_PERF(NRM) { \
float af = 0; double ad = 0; std::complex<float> ac = 0; \
Eigen::BenchTimer tf, td, tcf; tf.reset(); td.reset(); tcf.reset();\
for (int k=0; k<tries; ++k) { \
tf.start(); \
for (int i=0; i<iters; ++i) NRM(vf); \
for (int i=0; i<iters; ++i) { af += NRM(vf); } \
tf.stop(); \
} \
for (int k=0; k<tries; ++k) { \
td.start(); \
for (int i=0; i<iters; ++i) NRM(vd); \
for (int i=0; i<iters; ++i) { ad += NRM(vd); } \
td.stop(); \
} \
/*for (int k=0; k<std::max(1,tries/3); ++k) { \
tcf.start(); \
for (int i=0; i<iters; ++i) NRM(vcf); \
for (int i=0; i<iters; ++i) { ac += NRM(vcf); } \
tcf.stop(); \
} */\
std::cout << #NRM << "\t" << tf.value() << " " << td.value() << " " << tcf.value() << "\n"; \
@ -316,14 +323,17 @@ int main(int argc, char** argv)
std::cout << "\n";
}
y = 1;
std::cout.precision(4);
std::cerr << "Performance (out of cache):\n";
int s1 = 1024*1024*32;
std::cerr << "Performance (out of cache, " << s1 << "):\n";
{
int iters = 1;
VectorXf vf = VectorXf::Random(1024*1024*32) * y;
VectorXd vd = VectorXd::Random(1024*1024*32) * y;
VectorXcf vcf = VectorXcf::Random(1024*1024*32) * y;
VectorXf vf = VectorXf::Random(s1) * y;
VectorXd vd = VectorXd::Random(s1) * y;
VectorXcf vcf = VectorXcf::Random(s1) * y;
BENCH_PERF(sqsumNorm);
BENCH_PERF(stableNorm);
BENCH_PERF(blueNorm);
BENCH_PERF(pblueNorm);
BENCH_PERF(lapackNorm);
@ -332,13 +342,14 @@ int main(int argc, char** argv)
BENCH_PERF(bl2passNorm);
}
std::cerr << "\nPerformance (in cache):\n";
std::cerr << "\nPerformance (in cache, " << 512 << "):\n";
{
int iters = 100000;
VectorXf vf = VectorXf::Random(512) * y;
VectorXd vd = VectorXd::Random(512) * y;
VectorXcf vcf = VectorXcf::Random(512) * y;
BENCH_PERF(sqsumNorm);
BENCH_PERF(stableNorm);
BENCH_PERF(blueNorm);
BENCH_PERF(pblueNorm);
BENCH_PERF(lapackNorm);

1
doc/CMakeLists.txt
View File

@ -97,6 +97,7 @@ add_dependencies(doc-unsupported-prerequisites unsupported_snippets unsupported_
add_custom_target(doc ALL
COMMAND doxygen
COMMAND doxygen Doxyfile-unsupported
COMMAND ${CMAKE_COMMAND} -E copy ${Eigen_BINARY_DIR}/doc/html/group__TopicUnalignedArrayAssert.html ${Eigen_BINARY_DIR}/doc/html/TopicUnalignedArrayAssert.html
COMMAND ${CMAKE_COMMAND} -E rename html eigen-doc
COMMAND ${CMAKE_COMMAND} -E remove eigen-doc/eigen-doc.tgz
COMMAND ${CMAKE_COMMAND} -E tar cfz eigen-doc/eigen-doc.tgz eigen-doc

7
doc/snippets/DirectionWise_hnormalized.cpp Normal file
View File

@ -0,0 +1,7 @@
typedef Matrix<double,4,Dynamic> Matrix4Xd;
Matrix4Xd M = Matrix4Xd::Random(4,5);
Projective3d P(Matrix4d::Random());
cout << "The matrix M is:" << endl << M << endl << endl;
cout << "M.colwise().hnormalized():" << endl << M.colwise().hnormalized() << endl << endl;
cout << "P*M:" << endl << P*M << endl << endl;
cout << "(P*M).colwise().hnormalized():" << endl << (P*M).colwise().hnormalized() << endl << endl;

6
doc/snippets/MatrixBase_hnormalized.cpp Normal file
View File

@ -0,0 +1,6 @@
Vector4d v = Vector4d::Random();
Projective3d P(Matrix4d::Random());
cout << "v = " << v.transpose() << "]^T" << endl;
cout << "v.hnormalized() = " << v.hnormalized().transpose() << "]^T" << endl;
cout << "P*v = " << (P*v).transpose() << "]^T" << endl;
cout << "(P*v).hnormalized() = " << (P*v).hnormalized().transpose() << "]^T" << endl;

6
doc/snippets/MatrixBase_homogeneous.cpp Normal file
View File

@ -0,0 +1,6 @@
Vector3d v = Vector3d::Random(), w;
Projective3d P(Matrix4d::Random());
cout << "v = [" << v.transpose() << "]^T" << endl;
cout << "h.homogeneous() = [" << v.homogeneous().transpose() << "]^T" << endl;
cout << "(P * v.homogeneous()) = [" << (P * v.homogeneous()).transpose() << "]^T" << endl;
cout << "(P * v.homogeneous()).hnormalized() = [" << (P * v.homogeneous()).eval().hnormalized().transpose() << "]^T" << endl;

7
doc/snippets/VectorwiseOp_homogeneous.cpp Normal file
View File

@ -0,0 +1,7 @@
typedef Matrix<double,3,Dynamic> Matrix3Xd;
Matrix3Xd M = Matrix3Xd::Random(3,5);
Projective3d P(Matrix4d::Random());
cout << "The matrix M is:" << endl << M << endl << endl;
cout << "M.colwise().homogeneous():" << endl << M.colwise().homogeneous() << endl << endl;
cout << "P * M.colwise().homogeneous():" << endl << P * M.colwise().homogeneous() << endl << endl;
cout << "P * M.colwise().homogeneous().hnormalized(): " << endl << (P * M.colwise().homogeneous()).colwise().hnormalized() << endl << endl;

52
test/jacobisvd.cpp
View File

@ -315,16 +315,30 @@ void jacobisvd_inf_nan()
VERIFY(sub(some_inf, some_inf) != sub(some_inf, some_inf));
svd.compute(MatrixType::Constant(10,10,some_inf), ComputeFullU | ComputeFullV);
Scalar some_nan = zero<Scalar>() / zero<Scalar>();
VERIFY(some_nan != some_nan);
svd.compute(MatrixType::Constant(10,10,some_nan), ComputeFullU | ComputeFullV);
Scalar nan = std::numeric_limits<Scalar>::quiet_NaN();
VERIFY(nan != nan);
svd.compute(MatrixType::Constant(10,10,nan), ComputeFullU | ComputeFullV);
MatrixType m = MatrixType::Zero(10,10);
m(internal::random<int>(0,9), internal::random<int>(0,9)) = some_inf;
svd.compute(m, ComputeFullU | ComputeFullV);
m = MatrixType::Zero(10,10);
m(internal::random<int>(0,9), internal::random<int>(0,9)) = some_nan;
m(internal::random<int>(0,9), internal::random<int>(0,9)) = nan;
svd.compute(m, ComputeFullU | ComputeFullV);
// regression test for bug 791
m.resize(3,3);
m << 0, 2*NumTraits<Scalar>::epsilon(), 0.5,
0, -0.5, 0,
nan, 0, 0;
svd.compute(m, ComputeFullU | ComputeFullV);
m.resize(4,4);
m << 1, 0, 0, 0,
0, 3, 1, 2e-308,
1, 0, 1, nan,
0, nan, nan, 0;
svd.compute(m, ComputeFullU | ComputeFullV);
}
@ -340,11 +354,33 @@ void jacobisvd_underoverflow()
Matrix2d M;
M << -7.90884e-313, -4.94e-324,
0, 5.60844e-313;
JacobiSVD<Matrix2d> svd;
svd.compute(M,ComputeFullU|ComputeFullV);
jacobisvd_check_full(M,svd);
VectorXd value_set(9);
value_set << 0, 1, -1, 5.60844e-313, -5.60844e-313, 4.94e-324, -4.94e-324, -4.94e-223, 4.94e-223;
Array4i id(0,0,0,0);
int k = 0;
do
{
M << value_set(id(0)), value_set(id(1)), value_set(id(2)), value_set(id(3));
svd.compute(M,ComputeFullU|ComputeFullV);
jacobisvd_check_full(M,svd);
id(k)++;
if(id(k)>=value_set.size())
{
while(k<3 && id(k)>=value_set.size()) id(++k)++;
id.head(k).setZero();
k=0;
}
} while((id<int(value_set.size())).all());
#if defined __INTEL_COMPILER
#pragma warning pop
#endif
JacobiSVD<Matrix2d> svd;
svd.compute(M); // just check we don't loop indefinitely
// Check for overflow:
Matrix3d M3;
@ -353,7 +389,8 @@ void jacobisvd_underoverflow()
-8.7190887618028355e+307, -7.3453213709232193e+307, -2.4367363684472105e+307;
JacobiSVD<Matrix3d> svd3;
svd3.compute(M3); // just check we don't loop indefinitely
svd3.compute(M3,ComputeFullU|ComputeFullV); // just check we don't loop indefinitely
jacobisvd_check_full(M3,svd3);
}
void jacobisvd_preallocate()
@ -437,6 +474,7 @@ void test_jacobisvd()
// Test on inf/nan matrix
CALL_SUBTEST_7( jacobisvd_inf_nan<MatrixXf>() );
CALL_SUBTEST_10( jacobisvd_inf_nan<MatrixXd>() );
}
CALL_SUBTEST_7(( jacobisvd<MatrixXf>(MatrixXf(internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/2), internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/2))) ));

26
test/linearstructure.cpp
View File

@ -2,11 +2,16 @@
// for linear algebra.
//
// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com>
// Copyright (C) 2014 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
static bool g_called;
#define EIGEN_SPECIAL_SCALAR_MULTIPLE_PLUGIN { g_called = true; }
#include "main.h"
template<typename MatrixType> void linearStructure(const MatrixType& m)
@ -68,6 +73,24 @@ template<typename MatrixType> void linearStructure(const MatrixType& m)
VERIFY_IS_APPROX(m1.block(0,0,rows,cols) * s1, m1 * s1);
}
// Make sure that complex * real and real * complex are properly optimized
template<typename MatrixType> void real_complex(DenseIndex rows = MatrixType::RowsAtCompileTime, DenseIndex cols = MatrixType::ColsAtCompileTime)
{
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
RealScalar s = internal::random<RealScalar>();
MatrixType m1 = MatrixType::Random(rows, cols);
g_called = false;
VERIFY_IS_APPROX(s*m1, Scalar(s)*m1);
VERIFY(g_called && "real * matrix<complex> not properly optimized");
g_called = false;
VERIFY_IS_APPROX(m1*s, m1*Scalar(s));
VERIFY(g_called && "matrix<complex> * real not properly optimized");
}
void test_linearstructure()
{
for(int i = 0; i < g_repeat; i++) {
@ -80,5 +103,8 @@ void test_linearstructure()
CALL_SUBTEST_7( linearStructure(MatrixXi (internal::random<int>(1,EIGEN_TEST_MAX_SIZE), internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) );
CALL_SUBTEST_8( linearStructure(MatrixXcd(internal::random<int>(1,EIGEN_TEST_MAX_SIZE/2), internal::random<int>(1,EIGEN_TEST_MAX_SIZE/2))) );
CALL_SUBTEST_9( linearStructure(ArrayXXf (internal::random<int>(1,EIGEN_TEST_MAX_SIZE), internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) );
CALL_SUBTEST_10( real_complex<Matrix4cd>() );
CALL_SUBTEST_10( real_complex<MatrixXcf>(10,10) );
}
}

8
test/product.h
View File

@ -139,4 +139,12 @@ template<typename MatrixType> void product(const MatrixType& m)
// inner product
Scalar x = square2.row(c) * square2.col(c2);
VERIFY_IS_APPROX(x, square2.row(c).transpose().cwiseProduct(square2.col(c2)).sum());
// outer product
VERIFY_IS_APPROX(m1.col(c) * m1.row(r), m1.block(0,c,rows,1) * m1.block(r,0,1,cols));
VERIFY_IS_APPROX(m1.row(r).transpose() * m1.col(c).transpose(), m1.block(r,0,1,cols).transpose() * m1.block(0,c,rows,1).transpose());
VERIFY_IS_APPROX(m1.block(0,c,rows,1) * m1.row(r), m1.block(0,c,rows,1) * m1.block(r,0,1,cols));
VERIFY_IS_APPROX(m1.col(c) * m1.block(r,0,1,cols), m1.block(0,c,rows,1) * m1.block(r,0,1,cols));
VERIFY_IS_APPROX(m1.leftCols(1) * m1.row(r), m1.block(0,0,rows,1) * m1.block(r,0,1,cols));
VERIFY_IS_APPROX(m1.col(c) * m1.topRows(1), m1.block(0,c,rows,1) * m1.block(0,0,1,cols));
}

69
test/stable_norm.cpp
View File

@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
// Copyright (C) 2009-2014 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
@ -14,6 +14,21 @@ template<typename T> bool isNotNaN(const T& x)
return x==x;
}
template<typename T> bool isNaN(const T& x)
{
return x!=x;
}
template<typename T> bool isInf(const T& x)
{
return x > NumTraits<T>::highest();
}
template<typename T> bool isMinusInf(const T& x)
{
return x < NumTraits<T>::lowest();
}
// workaround aggressive optimization in ICC
template<typename T> EIGEN_DONT_INLINE T sub(T a, T b) { return a - b; }
@ -106,6 +121,58 @@ template<typename MatrixType> void stable_norm(const MatrixType& m)
VERIFY_IS_APPROX(vrand.rowwise().stableNorm(), vrand.rowwise().norm());
VERIFY_IS_APPROX(vrand.rowwise().blueNorm(), vrand.rowwise().norm());
VERIFY_IS_APPROX(vrand.rowwise().hypotNorm(), vrand.rowwise().norm());
// test NaN, +inf, -inf
MatrixType v;
Index i = internal::random<Index>(0,rows-1);
Index j = internal::random<Index>(0,cols-1);
// NaN
{
v = vrand;
v(i,j) = RealScalar(0)/RealScalar(0);
VERIFY(!isFinite(v.squaredNorm())); VERIFY(isNaN(v.squaredNorm()));
VERIFY(!isFinite(v.norm())); VERIFY(isNaN(v.norm()));
VERIFY(!isFinite(v.stableNorm())); VERIFY(isNaN(v.stableNorm()));
VERIFY(!isFinite(v.blueNorm())); VERIFY(isNaN(v.blueNorm()));
VERIFY(!isFinite(v.hypotNorm())); VERIFY(isNaN(v.hypotNorm()));
}
// +inf
{
v = vrand;
v(i,j) = RealScalar(1)/RealScalar(0);
VERIFY(!isFinite(v.squaredNorm())); VERIFY(isInf(v.squaredNorm()));
VERIFY(!isFinite(v.norm())); VERIFY(isInf(v.norm()));
VERIFY(!isFinite(v.stableNorm())); VERIFY(isInf(v.stableNorm()));
VERIFY(!isFinite(v.blueNorm())); VERIFY(isInf(v.blueNorm()));
VERIFY(!isFinite(v.hypotNorm())); VERIFY(isInf(v.hypotNorm()));
}
// -inf
{
v = vrand;
v(i,j) = RealScalar(-1)/RealScalar(0);
VERIFY(!isFinite(v.squaredNorm())); VERIFY(isInf(v.squaredNorm()));
VERIFY(!isFinite(v.norm())); VERIFY(isInf(v.norm()));
VERIFY(!isFinite(v.stableNorm())); VERIFY(isInf(v.stableNorm()));
VERIFY(!isFinite(v.blueNorm())); VERIFY(isInf(v.blueNorm()));
VERIFY(!isFinite(v.hypotNorm())); VERIFY(isInf(v.hypotNorm()));
}
// mix
{
Index i2 = internal::random<Index>(0,rows-1);
Index j2 = internal::random<Index>(0,cols-1);
v = vrand;
v(i,j) = RealScalar(-1)/RealScalar(0);
v(i2,j2) = RealScalar(0)/RealScalar(0);
VERIFY(!isFinite(v.squaredNorm())); VERIFY(isNaN(v.squaredNorm()));
VERIFY(!isFinite(v.norm())); VERIFY(isNaN(v.norm()));
VERIFY(!isFinite(v.stableNorm())); VERIFY(isNaN(v.stableNorm()));
VERIFY(!isFinite(v.blueNorm())); VERIFY(isNaN(v.blueNorm()));
VERIFY(!isFinite(v.hypotNorm())); VERIFY(isNaN(v.hypotNorm()));
}
}
void test_stable_norm()

2
test/upperbidiagonalization.cpp
View File

@ -35,7 +35,7 @@ void test_upperbidiagonalization()
CALL_SUBTEST_1( upperbidiag(MatrixXf(3,3)) );
CALL_SUBTEST_2( upperbidiag(MatrixXd(17,12)) );
CALL_SUBTEST_3( upperbidiag(MatrixXcf(20,20)) );
CALL_SUBTEST_4( upperbidiag(MatrixXcd(16,15)) );
CALL_SUBTEST_4( upperbidiag(Matrix<std::complex<double>,Dynamic,Dynamic,RowMajor>(16,15)) );
CALL_SUBTEST_5( upperbidiag(Matrix<float,6,4>()) );
CALL_SUBTEST_6( upperbidiag(Matrix<float,5,5>()) );
CALL_SUBTEST_7( upperbidiag(Matrix<double,4,3>()) );

21
unsupported/Eigen/MatrixFunctions
View File

@ -82,7 +82,9 @@ const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cos() const
\param[in] M a square matrix.
\returns expression representing \f$ \cos(M) \f$.
This function calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::cos().
This function computes the matrix cosine. Use ArrayBase::cos() for computing the entry-wise cosine.
The implementation calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::cos().
\sa \ref matrixbase_sin "sin()" for an example.
@ -123,6 +125,9 @@ differential equations: the solution of \f$ y' = My \f$ with the
initial condition \f$ y(0) = y_0 \f$ is given by
\f$ y(t) = \exp(M) y_0 \f$.
The matrix exponential is different from applying the exp function to all the entries in the matrix.
Use ArrayBase::exp() if you want to do the latter.
The cost of the computation is approximately \f$ 20 n^3 \f$ for
matrices of size \f$ n \f$. The number 20 depends weakly on the
norm of the matrix.
@ -177,6 +182,9 @@ the scalar logarithm, the equation \f$ \exp(X) = M \f$ may have
multiple solutions; this function returns a matrix whose eigenvalues
have imaginary part in the interval \f$ (-\pi,\pi] \f$.
The matrix logarithm is different from applying the log function to all the entries in the matrix.
Use ArrayBase::log() if you want to do the latter.
In the real case, the matrix \f$ M \f$ should be invertible and
it should have no eigenvalues which are real and negative (pairs of
complex conjugate eigenvalues are allowed). In the complex case, it
@ -232,7 +240,8 @@ const MatrixPowerReturnValue<Derived> MatrixBase<Derived>::pow(RealScalar p) con
The matrix power \f$ M^p \f$ is defined as \f$ \exp(p \log(M)) \f$,
where exp denotes the matrix exponential, and log denotes the matrix
logarithm.
logarithm. This is different from raising all the entries in the matrix
to the p-th power. Use ArrayBase::pow() if you want to do the latter.
If \p p is complex, the scalar type of \p M should be the type of \p
p . \f$ M^p \f$ simply evaluates into \f$ \exp(p \log(M)) \f$.
@ -391,7 +400,9 @@ const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sin() const
\param[in] M a square matrix.
\returns expression representing \f$ \sin(M) \f$.
This function calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::sin().
This function computes the matrix sine. Use ArrayBase::sin() for computing the entry-wise sine.
The implementation calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::sin().
Example: \include MatrixSine.cpp
Output: \verbinclude MatrixSine.out
@ -428,7 +439,9 @@ const MatrixSquareRootReturnValue<Derived> MatrixBase<Derived>::sqrt() const
The matrix square root of \f$ M \f$ is the matrix \f$ M^{1/2} \f$
whose square is the original matrix; so if \f$ S = M^{1/2} \f$ then
\f$ S^2 = M \f$.
\f$ S^2 = M \f$. This is different from taking the square root of all
the entries in the matrix; use ArrayBase::sqrt() if you want to do the
latter.
In the <b>real case</b>, the matrix \f$ M \f$ should be invertible and
it should have no eigenvalues which are real and negative (pairs of

19
unsupported/Eigen/src/NonLinearOptimization/LevenbergMarquardt.h
View File

@ -45,18 +45,24 @@ namespace LevenbergMarquardtSpace {
template<typename FunctorType, typename Scalar=double>
class LevenbergMarquardt
{
static Scalar sqrt_epsilon()
{
using std::sqrt;
return sqrt(NumTraits<Scalar>::epsilon());
}
public:
LevenbergMarquardt(FunctorType &_functor)
: functor(_functor) { nfev = njev = iter = 0; fnorm = gnorm = 0.; useExternalScaling=false; }
typedef DenseIndex Index;
struct Parameters {
Parameters()
: factor(Scalar(100.))
, maxfev(400)
, ftol(sqrt_(NumTraits<Scalar>::epsilon()))
, xtol(sqrt_(NumTraits<Scalar>::epsilon()))
, ftol(sqrt_epsilon())
, xtol(sqrt_epsilon())
, gtol(Scalar(0.))
, epsfcn(Scalar(0.)) {}
Scalar factor;
@ -72,7 +78,7 @@ public:
LevenbergMarquardtSpace::Status lmder1(
FVectorType &x,
const Scalar tol = sqrt_(NumTraits<Scalar>::epsilon())
const Scalar tol = sqrt_epsilon()
);
LevenbergMarquardtSpace::Status minimize(FVectorType &x);
@ -83,12 +89,12 @@ public:
FunctorType &functor,
FVectorType &x,
Index *nfev,
const Scalar tol = sqrt_(NumTraits<Scalar>::epsilon())
const Scalar tol = sqrt_epsilon()
);
LevenbergMarquardtSpace::Status lmstr1(
FVectorType &x,
const Scalar tol = sqrt_(NumTraits<Scalar>::epsilon())
const Scalar tol = sqrt_epsilon()
);
LevenbergMarquardtSpace::Status minimizeOptimumStorage(FVectorType &x);
@ -109,7 +115,6 @@ public:
Scalar lm_param(void) { return par; }
private:
static Scalar sqrt_(const Scalar& x) { using std::sqrt; return sqrt(x); }
FunctorType &functor;
Index n;

22
unsupported/test/NonLinearOptimization.cpp
View File

@ -1022,7 +1022,8 @@ void testNistLanczos1(void)
VERIFY_IS_EQUAL(lm.nfev, 79);
VERIFY_IS_EQUAL(lm.njev, 72);
// check norm^2
VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.430899764097e-25); // should be 1.4307867721E-25, but nist results are on 128-bit floats
std::cout.precision(30);
VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.4290986055242372e-25); // should be 1.4307867721E-25, but nist results are on 128-bit floats
// check x
VERIFY_IS_APPROX(x[0], 9.5100000027E-02);
VERIFY_IS_APPROX(x[1], 1.0000000001E+00);
@ -1043,7 +1044,7 @@ void testNistLanczos1(void)
VERIFY_IS_EQUAL(lm.nfev, 9);
VERIFY_IS_EQUAL(lm.njev, 8);
// check norm^2
VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.428595533845e-25); // should be 1.4307867721E-25, but nist results are on 128-bit floats
VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.430571737783119393e-25); // should be 1.4307867721E-25, but nist results are on 128-bit floats
// check x
VERIFY_IS_APPROX(x[0], 9.5100000027E-02);
VERIFY_IS_APPROX(x[1], 1.0000000001E+00);
@ -1262,8 +1263,8 @@ void testNistBoxBOD(void)
// check return value
VERIFY_IS_EQUAL(info, 1);
VERIFY_IS_EQUAL(lm.nfev, 31);
VERIFY_IS_EQUAL(lm.njev, 25);
VERIFY(lm.nfev < 31); // 31
VERIFY(lm.njev < 25); // 25
// check norm^2
VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.1680088766E+03);
// check x
@ -1342,10 +1343,6 @@ void testNistMGH17(void)
lm.parameters.maxfev = 1000;
info = lm.minimize(x);
// check return value
VERIFY_IS_EQUAL(info, 2);
VERIFY_IS_EQUAL(lm.nfev, 602 );
VERIFY_IS_EQUAL(lm.njev, 545 );
// check norm^2
VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.4648946975E-05);
// check x
@ -1354,6 +1351,11 @@ void testNistMGH17(void)
VERIFY_IS_APPROX(x[2], -1.4646871366E+00);
VERIFY_IS_APPROX(x[3], 1.2867534640E-02);
VERIFY_IS_APPROX(x[4], 2.2122699662E-02);
// check return value
VERIFY_IS_EQUAL(info, 2);
VERIFY(lm.nfev < 650); // 602
VERIFY(lm.njev < 600); // 545
/*
* Second try
@ -1832,8 +1834,8 @@ void test_NonLinearOptimization()
// NIST tests, level of difficulty = "Average"
CALL_SUBTEST/*_5*/(testNistHahn1());
CALL_SUBTEST/*_6*/(testNistMisra1d());
// CALL_SUBTEST/*_7*/(testNistMGH17());
// CALL_SUBTEST/*_8*/(testNistLanczos1());
CALL_SUBTEST/*_7*/(testNistMGH17());
CALL_SUBTEST/*_8*/(testNistLanczos1());
// // NIST tests, level of difficulty = "Higher"
CALL_SUBTEST/*_9*/(testNistRat42());

44
unsupported/test/levenberg_marquardt.cpp
View File

@ -787,16 +787,17 @@ void testNistMGH10(void)
LevenbergMarquardt<MGH10_functor> lm(functor);
info = lm.minimize(x);
// check return value
VERIFY_IS_EQUAL(info, 1);
VERIFY_IS_EQUAL(lm.nfev(), 284 );
VERIFY_IS_EQUAL(lm.njev(), 249 );
// check norm^2
VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 8.7945855171E+01);
// check x
VERIFY_IS_APPROX(x[0], 5.6096364710E-03);
VERIFY_IS_APPROX(x[1], 6.1813463463E+03);
VERIFY_IS_APPROX(x[2], 3.4522363462E+02);
// check return value
//VERIFY_IS_EQUAL(info, 1);
VERIFY_IS_EQUAL(lm.nfev(), 284 );
VERIFY_IS_EQUAL(lm.njev(), 249 );
/*
* Second try
@ -805,16 +806,17 @@ void testNistMGH10(void)
// do the computation
info = lm.minimize(x);
// check return value
VERIFY_IS_EQUAL(info, 1);
VERIFY_IS_EQUAL(lm.nfev(), 126);
VERIFY_IS_EQUAL(lm.njev(), 116);
// check norm^2
VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 8.7945855171E+01);
// check x
VERIFY_IS_APPROX(x[0], 5.6096364710E-03);
VERIFY_IS_APPROX(x[1], 6.1813463463E+03);
VERIFY_IS_APPROX(x[2], 3.4522363462E+02);
// check return value
//VERIFY_IS_EQUAL(info, 1);
VERIFY_IS_EQUAL(lm.nfev(), 126);
VERIFY_IS_EQUAL(lm.njev(), 116);
}
@ -866,15 +868,16 @@ void testNistBoxBOD(void)
lm.setFactor(10);
info = lm.minimize(x);
// check return value
VERIFY_IS_EQUAL(info, 1);
VERIFY_IS_EQUAL(lm.nfev(), 31);
VERIFY_IS_EQUAL(lm.njev(), 25);
// check norm^2
VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.1680088766E+03);
// check x
VERIFY_IS_APPROX(x[0], 2.1380940889E+02);
VERIFY_IS_APPROX(x[1], 5.4723748542E-01);
// check return value
VERIFY_IS_EQUAL(info, 1);
VERIFY(lm.nfev() < 31); // 31
VERIFY(lm.njev() < 25); // 25
/*
* Second try
@ -948,10 +951,6 @@ void testNistMGH17(void)
lm.setMaxfev(1000);
info = lm.minimize(x);
// check return value
// VERIFY_IS_EQUAL(info, 2); //FIXME Use (lm.info() == Success)
// VERIFY_IS_EQUAL(lm.nfev(), 602 );
VERIFY_IS_EQUAL(lm.njev(), 545 );
// check norm^2
VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.4648946975E-05);
// check x
@ -960,6 +959,11 @@ void testNistMGH17(void)
VERIFY_IS_APPROX(x[2], -1.4646871366E+00);
VERIFY_IS_APPROX(x[3], 1.2867534640E-02);
VERIFY_IS_APPROX(x[4], 2.2122699662E-02);
// check return value
// VERIFY_IS_EQUAL(info, 2); //FIXME Use (lm.info() == Success)
VERIFY(lm.nfev() < 700 ); // 602
VERIFY(lm.njev() < 600 ); // 545
/*
* Second try
@ -1035,10 +1039,6 @@ void testNistMGH09(void)
lm.setMaxfev(1000);
info = lm.minimize(x);
// check return value
VERIFY_IS_EQUAL(info, 1);
VERIFY_IS_EQUAL(lm.nfev(), 490 );
VERIFY_IS_EQUAL(lm.njev(), 376 );
// check norm^2
VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 3.0750560385E-04);
// check x
@ -1046,6 +1046,10 @@ void testNistMGH09(void)
VERIFY_IS_APPROX(x[1], 0.19126423573); // should be 1.9128232873E-01
VERIFY_IS_APPROX(x[2], 0.12305309914); // should be 1.2305650693E-01
VERIFY_IS_APPROX(x[3], 0.13605395375); // should be 1.3606233068E-01
// check return value
VERIFY_IS_EQUAL(info, 1);
VERIFY(lm.nfev() < 510 ); // 490
VERIFY(lm.njev() < 400 ); // 376
/*
* Second try