Hello,

I have an embedded system running VxWorks on PowerPC. I've been using Eigen in our algorithms but now need to use an extended data type for finding the roots of 16th order polynomial for better precision. Our PowerPC implements long doubles as the IBM 128bit floats causing the Polynomial Solver to always fail (returns 0 roots). The same algorithm works without problems on x86_64/x86 Linux, where long doubles are implemented as the gcc 80 bit extended precision. We've verified the problem by giving the Polynomial Solver on the PowerPC target a set of simple coefficients with a known set of roots, but it always returns 0 roots found.

Is there a way to get the Polynomial Solver working on PowerPC using extended precision? I should note that simply using doubles in the Polynomial Solver works as expected on PowerPC but lacks the precision necessary for this task.

Thanks,
Ira

I cannot help much here. Maybe you could start by running a modified version of the eigensolver_generic unit test on this platform. See file test/eigensolver_generic.cpp, and line 97 add:

CALL_SUBTEST_2( eigensolver(Matrix<long double,Dynamic,Dynamic>(s,s)) );

If it works, then do the same with the polynomial solver unit test (in unsupported/test) and see...

Also if you have a selfcontained test, I could try it with emulated 128bits doubles.

Thank you for the reply! I had some difficulty cross-compiling the unit test so I instead just create a similar self-test.

Given the following random matrix A:
0.814723686393179, 0.913375856139019, 0.27849821886704
0.905791937075619, 0.63235924622541, 0.546881519204984
0.126986816293506, 0.09754040499941, 0.957506835434298

EigenSolver< Matrix<long double, 4, 4> > ei0(A) returns the following result on x86 Linux.
(-0.187943833218033,0)
(1.75267034254289,0)
(0.839863258728032,0)

On PowerPC though, ei0.info() = 2 and no eigenvalues are found if I use long double. Any suggestions?

Thanks,
Ira

For reference, I'm using a cross-compile gcc compiler for the PowerPC;
powerpc-wrs-vxworks-g++ -v
Using built-in specs.
Target: powerpc-wrs-vxworks
Configured with: ../gcc-4.7.0/configure --prefix=/usr/local/share/powerpc-vxworks --with-headers=/home/darkknight/Code/crossCompile/WindRiver/vxworks-6.3/target/h --with-cpu-64=e300c2 --target=powerpc-wrs-vxworks --with-gnu-ld --with-gnu-as --enable-languages=c,c++ --disable-multilib --disable-shared
Thread model: vxworks
gcc version 4.7.0 (GCC)

I assume that "EigenSolver< Matrix<long double, 4, 4> > ei0(A)" is a typo and the actual code is:

EigenSolver< Matrix<long double, 3, 3> > ei0(A)

since you have a 3x3 matrix?

Anyway, the next step is thus to test the HessenbergDecomposition:

HessenbergDecomposition<Matrix<long double, 3, 3> > hess(A);

and looked at:

Matrix<long double, 3, 3> Q = hess.matrixQ();
Matrix<long double, 3, 3> H = hess.matrixH();

in particular we should have:

(A - Q * H * Q.transpose()).norm()

close to 0.

Yes, sorry! I copied the output but retyped that source line incorrectly.

I didn't have access to the hardware today, but I will test this out tomorrow. Also as a sanity check, I'm going to recheck some basic math operations using the 128 bit floats on PowerPC.

Thanks again.

Here are the results on PowerPC:
HessenbergDecomposition<Matrix<long double, 3, 3> > hess(A);
Matrix<long double, 3, 3> Q = hess.matrixQ();
Matrix<long double, 3, 3> H = hess.matrixH();
Q:
1 0 0
0 -0.99031531553884413199 -0.13883650747983855567
0 -0.13883650747983855567 0.99031531553884430963

H:
0.81472368639317913619 -0.94319581922472988111 0.14899113761156998592
-0.91465003404775693241 0.72722946069452536477 -0.57916508178232328774
0 -0.12982396757674949939 0.86263662096518256561

(A - Q * H * Q.transpose()).norm():
8.5461011298038620509e-16

Thanks.

OK, so the last test to isolate the issue is to try the RealSchur decomposition:

RealSchur<Matrix<long double, 3, 3> > shur(A);
Matrix<long double, 3, 3> U = shur.matrixU();
Matrix<long double, 3, 3> T = shur.matrixT();

look at U and T and check that U * T * U.transpose() is close to zero!