Compare commits

...

12 Commits

Author SHA1 Message Date
Benoit Jacob
907fba9ea9 bump 2010-07-16 22:25:08 -04:00
Gael Guennebaud
7e0b7b1f25 uncomment tests (sorry) 2010-07-16 11:50:40 +02:00
Gael Guennebaud
5ac00624ed adapted pactch from Nick Lewycky fixing a couple of issue with CLang 2010-07-15 22:48:44 +02:00
Gael Guennebaud
c0a0d2d181 add unit tests for 0 matrix 2010-07-15 20:00:50 +02:00
Gael Guennebaud
ec39a39cb6 fix LLT for zero matrix 2010-07-15 20:00:34 +02:00
Gael Guennebaud
ccc6731f86 fix use of rank in QR 2010-07-15 19:59:21 +02:00
Gael Guennebaud
b09bb50aeb fix QR solving with m>n 2010-07-15 19:52:11 +02:00
Gael Guennebaud
c7b8de77c0 fix aligned_delete for null pointers 2010-07-15 09:28:29 +02:00
Gael Guennebaud
c44bbabdcc fix LU and QR solve when rank==0, fix LLT when the matrix is purely 0 2010-07-12 16:42:38 +02:00
Gael Guennebaud
468863f3a0 fix bad Map in row-vector by matrix products 2010-07-09 20:27:05 +02:00
Gael Guennebaud
81f36b0e21 disable MSVC optimization when the underlying compiler is ICC 2010-07-09 19:36:43 +02:00
Benoit Jacob
c196a49f67 Added tag 2.0.14 for changeset e18e51d891 2010-06-22 22:23:39 -04:00
16 changed files with 61 additions and 28 deletions

View File

@@ -7,7 +7,7 @@ set(INCLUDE_INSTALL_DIR
"The directory where we install the header files" "The directory where we install the header files"
FORCE) FORCE)
set(EIGEN_VERSION_NUMBER "2.0.14") set(EIGEN_VERSION_NUMBER "2.0.15")
set(EIGEN_VERSION "${EIGEN_VERSION_NUMBER}") set(EIGEN_VERSION "${EIGEN_VERSION_NUMBER}")
set(CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake) set(CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake)

View File

@@ -139,7 +139,7 @@ inline bool MatrixBase<Derived>::any() const
template<typename Derived> template<typename Derived>
inline int MatrixBase<Derived>::count() const inline int MatrixBase<Derived>::count() const
{ {
return this->cast<bool>().cast<int>().sum(); return this->cast<bool>().template cast<int>().sum();
} }
#endif // EIGEN_ALLANDANY_H #endif // EIGEN_ALLANDANY_H

View File

@@ -132,11 +132,12 @@ void LLT<MatrixType>::compute(const MatrixType& a)
m_isInitialized = true; m_isInitialized = true;
return; return;
} }
m_matrix.col(0).end(size-1) = a.row(0).end(size-1).adjoint() / ei_real(m_matrix.coeff(0,0)); if(ei_real(m_matrix.coeff(0,0))>0)
m_matrix.col(0).end(size-1) = a.row(0).end(size-1).adjoint() / ei_real(m_matrix.coeff(0,0));
for (int j = 1; j < size; ++j) for (int j = 1; j < size; ++j)
{ {
x = ei_real(a.coeff(j,j)) - m_matrix.row(j).start(j).squaredNorm(); x = ei_real(a.coeff(j,j)) - m_matrix.row(j).start(j).squaredNorm();
if (x < cutoff) if (x <= cutoff)
{ {
m_isPositiveDefinite = false; m_isPositiveDefinite = false;
continue; continue;

View File

@@ -32,7 +32,7 @@ namespace Eigen
{ {
#define EIGEN_INSTANTIATE_PRODUCT(TYPE) \ #define EIGEN_INSTANTIATE_PRODUCT(TYPE) \
template static void ei_cache_friendly_product<TYPE>( \ template void ei_cache_friendly_product<TYPE>( \
int _rows, int _cols, int depth, \ int _rows, int _cols, int depth, \
bool _lhsRowMajor, const TYPE* _lhs, int _lhsStride, \ bool _lhsRowMajor, const TYPE* _lhs, int _lhsStride, \
bool _rhsRowMajor, const TYPE* _rhs, int _rhsStride, \ bool _rhsRowMajor, const TYPE* _rhs, int _rhsStride, \

View File

@@ -35,16 +35,6 @@ template<typename T> inline T ei_random_amplitude()
else return static_cast<T>(10); else return static_cast<T>(10);
} }
template<typename T> inline T ei_hypot(T x, T y)
{
T _x = ei_abs(x);
T _y = ei_abs(y);
T p = std::max(_x, _y);
T q = std::min(_x, _y);
T qp = q/p;
return p * ei_sqrt(T(1) + qp*qp);
}
/************** /**************
*** int *** *** int ***
**************/ **************/
@@ -292,4 +282,14 @@ inline bool ei_isApproxOrLessThan(long double a, long double b, long double prec
return a <= b || ei_isApprox(a, b, prec); return a <= b || ei_isApprox(a, b, prec);
} }
template<typename T> inline T ei_hypot(T x, T y)
{
T _x = ei_abs(x);
T _y = ei_abs(y);
T p = std::max(_x, _y);
T q = std::min(_x, _y);
T qp = q/p;
return p * ei_sqrt(T(1) + qp*qp);
}
#endif // EIGEN_MATHFUNCTIONS_H #endif // EIGEN_MATHFUNCTIONS_H

View File

@@ -614,7 +614,7 @@ struct ei_cache_friendly_product_selector<ProductType,1,LhsOrder,LhsAccess,RhsCo
else else
{ {
_res = ei_aligned_stack_new(Scalar, res.size()); _res = ei_aligned_stack_new(Scalar, res.size());
Map<Matrix<Scalar,DestDerived::SizeAtCompileTime,1,ColMajor> >(_res, res.size()) = res; Map<Matrix<Scalar,1,DestDerived::SizeAtCompileTime,ColMajor> >(_res, res.size()) = res;
} }
ei_cache_friendly_product_colmajor_times_vector(res.size(), ei_cache_friendly_product_colmajor_times_vector(res.size(),
&product.rhs().const_cast_derived().coeffRef(0,0), product.rhs().stride(), &product.rhs().const_cast_derived().coeffRef(0,0), product.rhs().stride(),
@@ -622,7 +622,7 @@ struct ei_cache_friendly_product_selector<ProductType,1,LhsOrder,LhsAccess,RhsCo
if (!EvalToRes) if (!EvalToRes)
{ {
res = Map<Matrix<Scalar,DestDerived::SizeAtCompileTime,1,ColMajor> >(_res, res.size()); res = Map<Matrix<Scalar,1,DestDerived::SizeAtCompileTime,ColMajor> >(_res, res.size());
ei_aligned_stack_delete(Scalar, _res, res.size()); ei_aligned_stack_delete(Scalar, _res, res.size());
} }
} }
@@ -677,7 +677,7 @@ struct ei_cache_friendly_product_selector<ProductType,1,LhsOrder,LhsAccess,RhsCo
else else
{ {
_lhs = ei_aligned_stack_new(Scalar, product.lhs().size()); _lhs = ei_aligned_stack_new(Scalar, product.lhs().size());
Map<Matrix<Scalar,Lhs::SizeAtCompileTime,1,ColMajor> >(_lhs, product.lhs().size()) = product.lhs(); Map<Matrix<Scalar,1,Lhs::SizeAtCompileTime,ColMajor> >(_lhs, product.lhs().size()) = product.lhs();
} }
ei_cache_friendly_product_rowmajor_times_vector(&product.rhs().const_cast_derived().coeffRef(0,0), product.rhs().stride(), ei_cache_friendly_product_rowmajor_times_vector(&product.rhs().const_cast_derived().coeffRef(0,0), product.rhs().stride(),
_lhs, product.lhs().size(), res); _lhs, product.lhs().size(), res);

View File

@@ -35,7 +35,7 @@ template<typename Lhs, typename Rhs,
? UpperTriangular ? UpperTriangular
: -1, : -1,
int StorageOrder = ei_is_part<Lhs>::value ? -1 // this is to solve ambiguous specializations int StorageOrder = ei_is_part<Lhs>::value ? -1 // this is to solve ambiguous specializations
: int(Lhs::Flags) & (RowMajorBit|SparseBit) : int(Lhs::Flags) & int(RowMajorBit|SparseBit)
> >
struct ei_solve_triangular_selector; struct ei_solve_triangular_selector;

View File

@@ -114,13 +114,13 @@ template<> EIGEN_STRONG_INLINE void ei_pstoreu<float>(float* to, const __m128&
template<> EIGEN_STRONG_INLINE void ei_pstoreu<double>(double* to, const __m128d& from) { _mm_storeu_pd(to, from); } template<> EIGEN_STRONG_INLINE void ei_pstoreu<double>(double* to, const __m128d& from) { _mm_storeu_pd(to, from); }
template<> EIGEN_STRONG_INLINE void ei_pstoreu<int>(int* to, const __m128i& from) { _mm_storeu_si128(reinterpret_cast<__m128i*>(to), from); } template<> EIGEN_STRONG_INLINE void ei_pstoreu<int>(int* to, const __m128i& from) { _mm_storeu_si128(reinterpret_cast<__m128i*>(to), from); }
#if defined(_MSC_VER) && (_MSC_VER <= 1500) && defined(_WIN64) #if defined(_MSC_VER) && (_MSC_VER <= 1500) && defined(_WIN64) && !defined(__INTEL_COMPILER)
// The temporary variable fixes an internal compilation error. // The temporary variable fixes an internal compilation error.
// Direct of the struct members fixed bug #62. // Direct of the struct members fixed bug #62.
template<> EIGEN_STRONG_INLINE float ei_pfirst<__m128>(const __m128& a) { return a.m128_f32[0]; } template<> EIGEN_STRONG_INLINE float ei_pfirst<__m128>(const __m128& a) { return a.m128_f32[0]; }
template<> EIGEN_STRONG_INLINE double ei_pfirst<__m128d>(const __m128d& a) { return a.m128d_f64[0]; } template<> EIGEN_STRONG_INLINE double ei_pfirst<__m128d>(const __m128d& a) { return a.m128d_f64[0]; }
template<> EIGEN_STRONG_INLINE int ei_pfirst<__m128i>(const __m128i& a) { int x = _mm_cvtsi128_si32(a); return x; } template<> EIGEN_STRONG_INLINE int ei_pfirst<__m128i>(const __m128i& a) { int x = _mm_cvtsi128_si32(a); return x; }
#elif defined(_MSC_VER) && (_MSC_VER <= 1500) #elif defined(_MSC_VER) && (_MSC_VER <= 1500) && !defined(__INTEL_COMPILER)
// The temporary variable fixes an internal compilation error. // The temporary variable fixes an internal compilation error.
template<> EIGEN_STRONG_INLINE float ei_pfirst<__m128>(const __m128& a) { float x = _mm_cvtss_f32(a); return x; } template<> EIGEN_STRONG_INLINE float ei_pfirst<__m128>(const __m128& a) { float x = _mm_cvtss_f32(a); return x; }
template<> EIGEN_STRONG_INLINE double ei_pfirst<__m128d>(const __m128d& a) { double x = _mm_cvtsd_f64(a); return x; } template<> EIGEN_STRONG_INLINE double ei_pfirst<__m128d>(const __m128d& a) { double x = _mm_cvtsd_f64(a); return x; }

View File

@@ -30,7 +30,7 @@
#define EIGEN_WORLD_VERSION 2 #define EIGEN_WORLD_VERSION 2
#define EIGEN_MAJOR_VERSION 0 #define EIGEN_MAJOR_VERSION 0
#define EIGEN_MINOR_VERSION 14 #define EIGEN_MINOR_VERSION 15
#define EIGEN_VERSION_AT_LEAST(x,y,z) (EIGEN_WORLD_VERSION>x || (EIGEN_WORLD_VERSION>=x && \ #define EIGEN_VERSION_AT_LEAST(x,y,z) (EIGEN_WORLD_VERSION>x || (EIGEN_WORLD_VERSION>=x && \
(EIGEN_MAJOR_VERSION>y || (EIGEN_MAJOR_VERSION>=y && \ (EIGEN_MAJOR_VERSION>y || (EIGEN_MAJOR_VERSION>=y && \

View File

@@ -188,7 +188,8 @@ template<> inline void ei_conditional_aligned_free<false>(void *ptr)
template<typename T> inline void ei_destruct_elements_of_array(T *ptr, std::size_t size) template<typename T> inline void ei_destruct_elements_of_array(T *ptr, std::size_t size)
{ {
// always destruct an array starting from the end. // always destruct an array starting from the end.
while(size) ptr[--size].~T(); if(ptr)
while(size) ptr[--size].~T();
} }
/** \internal delete objects constructed with ei_aligned_new /** \internal delete objects constructed with ei_aligned_new

View File

@@ -515,9 +515,10 @@ bool LU<MatrixType>::solve(
if(!ei_isMuchSmallerThan(c.coeff(row,col), biggest_in_c, m_precision)) if(!ei_isMuchSmallerThan(c.coeff(row,col), biggest_in_c, m_precision))
return false; return false;
} }
m_lu.corner(TopLeft, m_rank, m_rank) if(m_rank>0)
.template marked<UpperTriangular>() m_lu.corner(TopLeft, m_rank, m_rank)
.solveTriangularInPlace(c.corner(TopLeft, m_rank, c.cols())); .template marked<UpperTriangular>()
.solveTriangularInPlace(c.corner(TopLeft, m_rank, c.cols()));
// Step 4 // Step 4
result->resize(m_lu.cols(), b.cols()); result->resize(m_lu.cols(), b.cols());

View File

@@ -270,12 +270,18 @@ bool QR<MatrixType>::solve(
ei_assert(m_isInitialized && "QR is not initialized."); ei_assert(m_isInitialized && "QR is not initialized.");
const int rows = m_qr.rows(); const int rows = m_qr.rows();
ei_assert(b.rows() == rows); ei_assert(b.rows() == rows);
result->resize(rows, b.cols()); // enforce the computation of the rank
rank();
result->resize(m_qr.cols(), b.cols());
// TODO(keir): There is almost certainly a faster way to multiply by // TODO(keir): There is almost certainly a faster way to multiply by
// Q^T without explicitly forming matrixQ(). Investigate. // Q^T without explicitly forming matrixQ(). Investigate.
*result = matrixQ().transpose()*b; *result = matrixQ().transpose()*b;
if(m_rank==0)
return result->isZero();
if(!isSurjective()) if(!isSurjective())
{ {
// is result is in the image of R ? // is result is in the image of R ?

View File

@@ -124,4 +124,8 @@ void test_cholesky()
CALL_SUBTEST( cholesky(MatrixXf(17,17)) ); CALL_SUBTEST( cholesky(MatrixXf(17,17)) );
CALL_SUBTEST( cholesky(MatrixXd(33,33)) ); CALL_SUBTEST( cholesky(MatrixXd(33,33)) );
} }
MatrixXf m = MatrixXf::Zero(10,10);
VectorXf b = VectorXf::Zero(10);
VERIFY(!m.llt().isPositiveDefinite());
} }

View File

@@ -132,4 +132,10 @@ void test_lu()
CALL_SUBTEST( lu_invertible<MatrixXcf>() ); CALL_SUBTEST( lu_invertible<MatrixXcf>() );
CALL_SUBTEST( lu_invertible<MatrixXcd>() ); CALL_SUBTEST( lu_invertible<MatrixXcd>() );
} }
MatrixXf m = MatrixXf::Zero(10,10);
VectorXf b = VectorXf::Zero(10);
VectorXf x = VectorXf::Random(10);
VERIFY(m.lu().solve(b,&x));
VERIFY(x.isZero());
} }

View File

@@ -48,4 +48,11 @@ void test_product_large()
MatrixXf a = MatrixXf::Random(10,4), b = MatrixXf::Random(4,10), c = a; MatrixXf a = MatrixXf::Random(10,4), b = MatrixXf::Random(4,10), c = a;
VERIFY_IS_APPROX((a = a * b), (c * b).eval()); VERIFY_IS_APPROX((a = a * b), (c * b).eval());
} }
{
MatrixXf mat1(10,10); mat1.setRandom();
MatrixXf mat2(32,10); mat2.setRandom();
MatrixXf result = mat1.row(2)*mat2.transpose();
VERIFY_IS_APPROX(result, (mat1.row(2)*mat2.transpose()).eval());
}
} }

View File

@@ -75,4 +75,11 @@ void test_qr()
mat << 1, 1, 1, 2, 2, 2, 1, 2, 3; mat << 1, 1, 1, 2, 2, 2, 1, 2, 3;
VERIFY(!mat.qr().isFullRank()); VERIFY(!mat.qr().isFullRank());
} }
{
MatrixXf m = MatrixXf::Zero(10,10);
VectorXf b = VectorXf::Zero(10);
VectorXf x = VectorXf::Random(10);
VERIFY(m.qr().solve(b,&x));
VERIFY(x.isZero());
}
} }