I have trouble when I use different versions of Lapack to calculate the same issue. I compared the results of the programs, compiling with Lapack 3.0, 3.1.1, and 3.2.1. Furthermore, I performed the same test with Intel MKL.
After compared these different versions, Lapack 3.0 seems more stable than the other versions, and I think it is possible that I encountered some bug in Lapack 3.1.1 and 3.2.1.
You tried this calculation on 4 different linear algebra libraries, and 3 out of 4 times you
got a result that, if I understand you correctly, is out of bounds from a theoretical point of
view. (And all the results apart from the one you report from Lapack 3.0 agree to within a
margin of machine noise.) That strongly suggests that there's a bug in your code.
To debug: put in print statements to see where the results diverge. If you can narrow it
down to a particular function, check the documentation for that function carefully:
Are you feeding it valid input? Are you interpreting the return value correctly?