Segmentation fault in gesv when matrices are too large

[scriabin3]

The following Fortran program uses gesv to solve the equation A x = b for random complex matrices. I have set A to be a 1600 x 1600 matrix and b to be 1600 x 160 (i.e there are 160 right hand sides in the equation). As is, the program runs correctly, but when I increase the 160 to 170I get a segmentation fault. Is this a bug or have I hit some limit I'm not aware of?

! Example program that uses the Intel MKL subroutine library to solve the matrix

! equation

! A x = b

!

! for random complex matrices: A is n x n, while x and b are n x nrhs.

!

! To call routines from the LAPACK one needs to "use lapack95" module

! For 64 bit computing:

! export FPATH=/opt/intel/Compiler/11.1/080/Frameworks/mkl/include/em64t/lp64

! ifort -r8 test.f90 -lmkl_core -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_lapack95_lp64

program main

use lapack95, only: gesv

implicit none

integer :: n, nrhs

complex, allocatable :: a(:,:), aa(:,:), b(:,:), bb(:,:)

real, allocatable :: ar(:,:),ai(:,:),br(:,:),bi(:,:)

n=1600; nrhs=170

allocate( ar(n,n),ai(n,n),a(n,n), aa(n,n), br(n,nrhs), bi(n,nrhs), b(n,nrhs), bb(n,nrhs) )

call random_number(ar) ! n x n real matrix of random numbers

call random_number(ai) ! n x n real matrix of random numbers

call random_number(br) ! n x nrhs real matrix of random numbers

call random_number(bi) ! n x nrhs real matrix of random numbers

aa=ar+(0.,1.)*ai ! n x n complex matrix of random numbers

bb=br+(0.,1.)*bi ! n x n complex matrix of random numbers

a=aa; b=bb ! make copies

call gesv(a,b) ! solve a x = b

write(*,*) matmul(aa,b)-bb ! check the result: output should be a zero matrix

end program main

[TimP]

The first suspicion would be that you didn't give attention to the stack limit setting in your shell, or to the KMP_STACKSIZE (thread stack size limit) environment variable, which defaults to 4MB in the intel64 library.

[Gennady_F_Intel]

- the simialr symptoms may happened when you linked ilp64 libs without /4I8 compiler options

--Gennady

[scriabin3]

Thank you gentlemen. Unfortunately I have never been educated about stacks. Looking on the web, I found that for OSX one can increase stack size up to the "hard limit" with

$ ulimit -s hard

I tried this, and indeed it worked, allowing me to increase nrhs from 160 to 1300. But with nrhs=1400, the program bombs with an "illegal instruction" message. I guess it's reached the hard limit? To go beyond the hard limit I read that one needs to modify and recompile the Mac OSX kernel. Now that's going beyond the call of duty for me. Is there are an easier way? I know I can get around the problem simply by solving Ax = b multiple times with different b's in order to keep nrhs small, but that seems like an ugly solution. Is there are better way? Is there even a recommended way from the MKL team?

[Vladimir_Koldakov__I · Accepted solution]

Hi,

Try to avoid temporary matrices that are placed on the stack.

Replace

write(*,*) matmul(aa,b)-bb ! check the result: output should be a zero matrix

With a few explicit assignments:

mm = matmul(aa,b)
mm2 = mm - bb
write(*,*) mm2 ! check the result: output should be a zero matrix

Of course, declare:

complex, allocatable :: mm(:,:), mm2(:,:)

and allocate:

allocate( mm(n,nrhs), mm2(n,nrhs) )

The following code passes without changing thestack size (not needed: $ ulimit -s hard):

program main
use lapack95, only: gesv
implicit none
integer :: n, nrhs
complex, allocatable :: a(:,:), aa(:,:), b(:,:), bb(:,:)
complex, allocatable :: mm(:,:), mm2(:,:)
real, allocatable :: ar(:,:),ai(:,:),br(:,:),bi(:,:)
n=1600; nrhs=1700
allocate( ar(n,n),ai(n,n),a(n,n), aa(n,n), br(n,nrhs), bi(n,nrhs), b(n,nrhs), bb(n,nrhs) )
allocate( mm(n,nrhs), mm2(n,nrhs) )
call random_number(ar) ! n x n real matrix of random numbers
call random_number(ai) ! n x n real matrix of random numbers
call random_number(br) ! n x nrhs real matrix of random numbers
call random_number(bi) ! n x nrhs real matrix of random numbers
aa=ar+(0.,1.)*ai ! n x n complex matrix of random numbers
bb=br+(0.,1.)*bi ! n x n complex matrix of random numbers
a=aa; b=bb ! make copies
call gesv(a,b) ! solve a x = b
mm = matmul(aa,b)
mm2 = mm - bb
write(*,*) mm2 ! check the result: output should be a zero matrix
end program main

Thanks,

Vladimir

[scriabin3]

Thank you Valodya, that was very helpful and educational too.