Stack overflow when solving large system using Julia

Question

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.

Only for small systems. Up until 33 x 33 works fine, however a 34 x 34 system gives a similar error. BTW for Julia [L,U] must be replaced by (L,U). — Aeronaelius, Jul 05 '13 at 02:26

score 6 · Answer 1 · answered Jul 05 '13 at 06:30

6

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

answered Jul 05 '13 at 06:30

ViralBShah

317
2
6

Hi, thanks for the fix. However in the discussion with Xianyi I saw that you were having this problem for also a random matrix multiplication. I don't have a problem with that. My problem is strictly with solving a large system: x=A\b. Setting blas_set_num_threads(1) however circumvented this issue, but I am guessing this kills performance? – Aeronaelius Jul 05 '13 at 22:33
I tried this in Julia 1.0.5, and works just fine. No need to change the threads or anything. – ViralBShah Sep 30 '19 at 19:55

score 1 · Answer 2 · answered Jul 22 '13 at 09:33

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.

answered Jul 22 '13 at 09:33

Aeronaelius

1,291
3
12
31

Given that this works out of the box, is it possible to accept the solution? I do not frequent stackoverflow, so I am not sure what is the equivalent of "closing or accepting" here. – ViralBShah Aug 23 '21 at 19:03

Stack overflow when solving large system using Julia

2 Answers2