| Re: [eigen] 3.0.1 released! |
[ Thread Index |
Date Index
| More lists.tuxfamily.org/eigen 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
construct_elements_of_array()
..
-> for (size_t i=0; i < size; ++i) ::new (ptr + i) T;
..
called through
aligned_stack_memory_handler()
...
-> 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.
II.
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):
(a)
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
(b)
rot1.c() = RealScalar(1) / sqrt(1 + abs2(u)); //was
rot1.c() = RealScalar(1) / sqrt(static_cast<RealScalar>(1) + abs2(u)); //became
(c)
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
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.
Hi,
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.
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.
II.
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?
III.
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)
(m.choose(rule))
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:
rule(x)
{
if sin(x)> 0.5 return true; else return false;
}
IV.
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