Re: [eigen] 3.0.1 released!

[ Thread Index | Date Index | More Archives ]

I've done few tests on Eigen + mpreal usage. 

I've checked only LU & SVD decompositions. 
Both crashes on memory allocation when matrix is not square. 

Debugger (MSVC 2010) traces problem to 
-> for (size_t i=0; i < size; ++i) ::new (ptr + i) T;

called through 

-> if(NumTraits<T>::RequireInitialization)
        Eigen::internal::construct_elements_of_array(m_ptr, size);

Simple example:

typedef Eigen::Matrix<mpfr::mpreal, Eigen::Dynamic, Eigen::Dynamic> MatrixXXmp;

MatrixXXmp A(2,3);

A << 8, 1, 6,
             3, 5, 7;
Eigen::JacobiSVD<MatrixXXmp> svd(A); // <- program crashes

Eigen 3.0.0 was working fine - 3.0.1 has new memory allocation.

I would like to express some critic regarding changes introduced by David H. Bailey in commit 4091 and others (These particular examples are taken from JacobiSVD.h):

RealScalar d = ....
rot1.s() = d > 0 ? 1 : -1; // was
rot1.s() = d > static_cast<RealScalar>(0) ? static_cast<RealScalar>(1) : static_cast<RealScalar>(-1) ; // became

rot1.c() = RealScalar(1) / sqrt(1 + abs2(u)); //was
rot1.c() = RealScalar(1) / sqrt(static_cast<RealScalar>(1) + abs2(u)); //became

rot1.c() = 0; //was 
rot1.c() = static_cast<RealScalar>(0); //became

Well designed custom scalar types include routines specifically optimized for the cases when one of the operands is of standard type.. 
In these cases math. operations can be done much faster. This is also true for conditional operations (<, >, etc.) and copy constructors.

For example, when we want to add integer to arbitrary precision number: 

mpreal a;

// compiler assumes 1 is of int type and calls optimized +=(int) operator, 
// it is very easy to add simple int to arbitrary precision number.
// no additional memory allocations, constructor call, etc.
a += 1; 

// unnecessary call of constructor/memory allocations, 
// but the worst - now we have to use heavy artillery to add numbers, 
// since they both of arbitrary precision now - optimization is disabled
a += mpreal(1); 

Actually huge part of arbitrary precision libraries are devoted to such optimized cases (just look to MPFR low level API) - otherwise we get significant drop in speed. 

There is only one case when we need such conversion, when operand is floating point literal, e.g.

a += 3.01;            // incorrect, since 3.01 is converted to double first.
a += mpreal("3.01");  // correct. string are being translated directly in multi-precision number. 

In upcoming C++0x string wrapping won't be necessary anymore - we will be able to create literals of custom type.

So, I'm back to 3.0.0 for a moment.

On Tue, May 31, 2011 at 12:02 PM, Pavel Holoborodko <pavel@xxxxxxxxxxxxxxx> wrote:
The later includes exceptions safety and the automatic uses of the
math functions declared in the scalar type's namespace.

Eigen 3.0.0 called math function(s) from std namespace (e.g. std::ceil ) with RealScalar as a parameter -  which wasn't right for custom scalars.
I was going to report this bug today, but 3.0.1 already fixed this - thank you so much!!

I'm testing Eigen step-by-step regarding arbitrary precision support - I'll let you know if something will come up. 

There are other question I want to discuss:

I needed (Moore-Penrose) pseudo inverse lately.. There are two versions available:

Both of them seems little outdated/flawed, (first doesn't support custom scalar types, second doesn't work with Eigen 3.x and both are not optimized) .  

To cure the situation I've added two functions to JacobiSVD class:

MatrixType pinverse(const RealScalar& tolerance) const; // use user-defined tolerance
MatrixType pinverse() const; // use default tolerance

Please see modified JacobiSVD.h in attach, whether this  pinverse() worth to be added in Eigen.

I see there is some work was done on porting LAPACK routines on Eigen natively (not just linking to Fortran binaries), e.g. getrf (\lapack\lu.cpp).

Although this is huge chunk of work it is the only decent way for porting LAPACK to C++ - to be implemented in Eigen natively, benefiting from its generic architecture, expressiveness, etc.

Besides I am sure this is the only and correct way to enable LAPACK with multi-precision. 

I'm very interested in this direction of development, to say the least.

Are there any plans/priorities in this direction?

Just quick question - maybe I'm missing something.
Is there any way to select matrix coefficients using custom binary operation as a criterion?

For example, to find sum of all elements greater than tolerance:

(m > tolerance).sum();

or choose elements based on arbitrary custom rule (proxy)


where rule can be an unary function which accepts one coefficient at a time, returning true/false (choose/not choose).
E.g. to select coefficients with sin > 0.5:

 if sin(x)> 0.5 return true; else return false;

Could you update URL of  MPFR C++ to  in "unsupported/Eigen/MPRealSupport"?
I've changed numerical URLs on my site to be more meaningful. 

I've been working with and studying Eigen for a few months. I must say this is absolutely the best C++ library I've ever saw. 
Both in design and masterful implementation. I really enjoy working with.
Thank you.

Pavel Holoborodko

On Mon, May 30, 2011 at 10:43 PM, Gael Guennebaud <gael.guennebaud@xxxxxxxxx> wrote:

Eigen 3.0.1 is released!

The source archive is at:

In addition to various minor bug fixes, this release brings official
support for gcc 4.6 and ARM NEON as well as an improved support for
custom scalar types.
The later includes exceptions safety and the automatic uses of the
math functions declared in the scalar type's namespace.
The support for ARM NEON has been possible thanks to the GCC Compile
Farm Project.

Complete changelog:


Mail converted by MHonArc 2.6.19+