CFD Online Discussion Forums - floating point exception error for eigenvalue calculation in applications/test/tensor

I am using OpenFOAM 4.x and found that applications/test/tensor make a floating point exception error.

in the below code, vector e = eigenValues(t6) makes the error

Code:

    tensor t6(1,0,-4,0,5,4,-4,4,3);

    //tensor t6(1,2,0,2,5,0,0,0,0);



    Info<< "tensor " << t6 << endl;

    vector e = eigenValues(t6);

    Info<< "eigenvalues " << e << endl;

problematic part is in the member function eigenValues(const tensor& t) of src/OpenFOAM/primitives/Tensor/tensor/tensor.C

in the code below, Q makes zero for the test matrix (1 0 -4 0 5 4 -4 4 3) which has eigen values (-3, 3, 9). So when checking mag(PPP/QQ -1), it leads to floating point exception error before going to the third conditional PPP > QQ.

The same test code works for OpenFOAM 2.3.0 but not work for OpenFOAM 4.x (also 3.0.x). I don't know much about the OpenFOAM code but it looks like a bug after version 3.0.x. Maybe Henry Weller can fix this behavior.

Code:

    // non-diagonal matrix

    else

    {

        // Coefficients of the characteristic polynmial

        // x^3 + a*x^2 + b*x + c = 0

        scalar a =

           - t.xx() - t.yy() - t.zz();



        scalar b =

            t.xx()*t.yy() + t.xx()*t.zz() + t.yy()*t.zz()

          - t.xy()*t.yx() - t.yz()*t.zy() - t.zx()*t.xz();



        scalar c =

          - t.xx()*t.yy()*t.zz()

          - t.xy()*t.yz()*t.zx() - t.xz()*t.zy()*t.yx()

          + t.xx()*t.yz()*t.zy() + t.yy()*t.zx()*t.xz() + t.zz()*t.xy()*t.yx();



        // Auxillary variables

        scalar aBy3 = a/3;



        scalar P = (a*a - 3*b)/9; // == -p_wikipedia/3

        scalar PPP = P*P*P;



        scalar Q = (2*a*a*a - 9*a*b + 27*c)/54; // == q_wikipedia/2

        scalar QQ = Q*Q;



        // Three identical roots

        if (mag(P) < SMALL && mag(Q) < SMALL)

        {

            return vector(- aBy3, - aBy3, - aBy3);

        }



        // Two identical roots and one distinct root

        else if (mag(PPP/QQ - 1) < SMALL)

        {

            scalar sqrtP = sqrt(P);

            scalar signQ = sign(Q);



            i = ii = signQ*sqrtP - aBy3;

            iii = - 2*signQ*sqrtP - aBy3;

        }