Re: [eigen] Scalar vs. int / double in SelfAdjointEigenSolver.h

[Pavel Holoborodko <pavel@xxxxxxxxxxxxxxx>]

Of course, we still need explicit conversions for constants:

const Scalar s_inv3 = Scalar(1.0)/Scalar(3.0);
const Scalar s_sqrt3 = sqrt(Scalar(3.0));

But not in expressions involving custom-type variables (e.g. just 1/x is fine) nor in conditionals.

On Thu, Aug 11, 2016 at 3:40 PM, Pavel Holoborodko <pavel@xxxxxxxxxxxxxxx> wrote:

...which may not the case for user defined type.

From my experience of creating custom scalar type (mpreal) - 

I can say that keeping code without explicit conversions to Scalar is very beneficial.

Just because operations involving standard types can be greatly optimized.

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

mpreal a;

// compiler 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 GMP/MPFR low level API).

Otherwise we get significant drop in speed. 

This is also true for conditional operations (<, >, etc.).

Please let the custom scalar types handle such special cases on its own, so that it can produce fast code!

 

On Thu, Aug 11, 2016 at 1:22 AM, Peter <list@xxxxxxxxxxxxxxxxx> wrote:

Dear All,

SelfAdjointEigenSolver.h  assumes that integer will implicitly converted to Scalar,
which may not the case for user defined type.

E.g.
 if(td==0)
 if(e2==0)
 while (start>0 && m_subdiag[start-1]!=0)
....

Shouldn't those lines read

 if(td==Scalar(0))
 if(e2==0)
 while (start>Scalar(0) && m_subdiag[start-1]!=Scalar(0))
....

In addition  void computeRoots(const MatrixType& m, VectorType& roots)
assume that double can explicitly converted to Scalar.


const Scalar s_inv3 = Scalar(1.0)/Scalar(3.0);
const Scalar s_sqrt3 = sqrt(Scalar(3.0));


Is there a reason not using

const Scalar s_inv3  = Scalar(1)/Scalar(3);
const Scalar s_sqrt3 = sqrt(Scalar(3));

However, there's also RealScalar(0.5).
I'm not sure whether one should use 1/RealScalar(2) here.

The attached version contains the corresponding changes based on eigen-eigen-dc6cfdf9bcec .

Best regards,
Peter