Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Stack overflow when solving large system using Julia

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

  • Linux Mint 13 KDE, 64bit
  • Installed LLVM 3.2 and Clang 3.2 from a PPA: ppa:kxstudio-team/builds
  • Compiled Julia 0.2.0-2429.rb0a9ea79 from source

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.

like image 902
Aeronaelius Avatar asked Jul 04 '13 22:07

Aeronaelius


2 Answers

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).

like image 52
ViralBShah Avatar answered Sep 18 '22 13:09

ViralBShah


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.

like image 23
Aeronaelius Avatar answered Sep 19 '22 13:09

Aeronaelius