I am trying to solve a random linear system with a large square system matrix using Octave and Julia. Because the syntax of Octave and Julia are quite similar I run the following code in both a Octave shell and a Julia shell:
N = 5000;
A = rand(N, N);
b = rand(N, 1);
x = A\b;
r = norm(A*x - b)/norm(b)
Octave returns r in the neighborhood of 1e-12. Julia on the other hand returns an error:
ERROR: stack overflow
in getrf! at linalg/lapack.jl:342
in LU at linalg/factorization.jl:134
in \ at linalg/dense.jl:518
The backslash operator does work in Julia for smaller systems (e.g. 10 x 10), however a 50 x 50 system already gives an error. As far as I know both Octave and Julia use BLAS and LAPACK, so I am rather confused why Julia is unable to perform this task. Can someone please tell me how I can fix this?
System information
EDIT
The problem has been solved now that OpenBLAS 0.2.7 is out. When re-compiling Julia make sure that Julia either uses a system version of OpenBLAS >=0.2.7 or that Julia internally compiles its own version of OpenBLAS >=0.2.7.
As I mentioned in the issue (https://github.com/JuliaLang/julia/issues/3630), this is most likely the same openblas threading bug as discussed in https://github.com/xianyi/OpenBLAS/issues/221.
There is a tentative fix on the openblas develop branch, which sets a larger stack size.
For now, do blas_set_num_threads(1)
.
Now that the new version of OpenBLAS: 0.2.7 is out I have compiled Julia anew. Unfortunately this did not amount to anything as Julia still uses OpenBLAS 0.2.6. However it is possible to use a system version of OpenBLAS while compiling Julia, instead of letting Julia download a version and compile it on its own. This way I made Julia use 0.2.7 instead of 0.2.6 and now the problem I was having is solved. No more stack overflows.
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With