Problems with memory access/allocation when using LAPACK DSYEVR. Code used is FORTRAN

Question

I have been writing a code in FORTRAN but I am having problems using the lapack dsyevr:

http://netlib.sandia.gov/lapack/double/dsyevr.f

The problems I am getting seem to be linked to memory allocation issues, specifically I believe to do with the output arrays the dsyevr produces (including A which is an input and an output).

I have tried to write a simplified code to demonstrate the issues I am seeing. Please let me know if any of it needs clarification. The code is names prof1.f90 and calls the dsyevr function:

PROGRAM prog1

implicit none

real(kind=8), allocatable :: W(:) 
real(kind=8), allocatable :: Z(:,:) 
real(kind=8), allocatable :: A(:,:)
integer(kind=8) :: n, info, il, iu, m, lwork, liwork
integer(kind=8) :: i, k, p, q, nu
real(kind=8) :: abstol, vl, vu
real(kind=8), allocatable :: work(:)
integer, allocatable :: isuppz(:), iwork(:)

n = 3

allocate(W(3),Z(3,3),A(n,n),stat=info)
if (info .ne. 0) stop "error allocating arrays"

A(1,1)=3.78136524999999994E-003
A(1,2)=0.0000000000000000
A(1,3)=-7.92918150000000038E-004
A(2,1)=0.0000000000000000
A(2,2)=5.20293929999999984E-003
A(2,3)=0.0000000000000000
A(3,1)=-7.92918150000000038E-004
A(3,2)=0.0000000000000000
A(3,3)=3.78136524999999994E-003
vl = 1.06451084056294826E-313
vu = 0.0
il = 4294967297
iu = 8839891
m = 140733655445712
W(1) = 2.98844710000000001E-003  
W(2) = 4.57428340000000030E-003  
W(3) = 5.20293929999999984E-003
Z(1,1) = 8.65587596665713699E-317  
Z(1,2) = 8.65587596665713699E-317
Z(1,3) = 1.58101006669198894E-322  
Z(2,1) = 1.58101006669198894E-322   
Z(2,2) = 0.0000000000000000
Z(2,3) = 8.65569415049946741E-317  
Z(3,1) = 4.24400777097956191E-314  
Z(3,2) = 4.79243676466009148E-322
Z(3,3) = 3.51391740150311405E-316
lwork = -1
liwork = -1
abstol = 1d-5

allocate(work(1),iwork(1),isuppz(6))

call dsyevr('V','A','U',n,A,n,vl,vu,il,iu,abstol,m,W,Z,n,isuppz,work,lwork,iwork,liwork,info)
if (info .ne. 0) stop "error obtaining work array dimensions"

lwork = work(1)
liwork = iwork(1)
deallocate(work,iwork)
allocate(work(lwork),iwork(liwork),stat=info)

if (info .ne. 0) stop "error allocating work arrays"

call dsyevr('V','A','U',n,A,n,vl,vu,il,iu,abstol,m,W,Z,n,isuppz,work,lwork,iwork,liwork,info)
if (info .ne. 0) stop "error diagonalizing the hamiltonian"

deallocate(A,work,iwork,isuppz)

END PROGRAM prog1

In the code above the dsyevr function is called twice, the first time which is just to get the dimensions of the work etc... matrices runs correctly however the second time when it is called it returns the following error

*** glibc detected *** ./PROGRAM: munmap_chunk(): invalid pointer: 0x000000000134cc20 ***

There is also a Backtrace and MemoryMap which I can provide.

If it is useful the makefile I have been using is given below. The program is created using the line:

make PROGRAM

Makefile:

FC = gfortran
FCFLAGS = -g -fbounds-check
FCFLAGS = -O2
FCFLAGS += -I/usr/include

%: %.o
    $(FC) $(FCFLAGS) -o $@ $^ $(LDFLAGS)

%.o: %.f90
    $(FC) $(FCFLAGS) -c $< -fno-range-check

%.o: %.F90
    $(FC) $(FCFLAGS) -c $<

.PHONY: clean veryclean

PROGRAM:        prog1.f90 prog1.o
    $(FC) $(FCFLAGS) -o $@ prog1.o $(LIBS)  -Wl,--start-group -L$(MKLROOT)/lib/intel64 -lmkl_gf_ilp64 -lmkl_core -lmkl_sequential -Wl,--end-group -lpthread


clean:
    rm -f *.o *.mod *.MOD

veryclean: clean
    rm -f *~ $(PROGRAM)

where $MKLROOT is /opt/intel/composer_xe_2011_sp1.8.273/mkl

If it is useful I have used valgrind:

 valgrind --tool=memcheck --db-attach=yes ./PROGRAM

And I have found the following error:

==31069== 
==31069== Invalid write of size 8
==31069==    at 0x57B9F92: mkl_lapack_dsyevr (in /opt/intel/composer_xe_2011_sp1.8.273    /mkl/lib/intel64/libmkl_core.so)
==31069==    by 0x4D9A580: DSYEVR (in /opt/intel/composer_xe_2011_sp1.8.273/mkl/lib    /intel64/libmkl_gf_ilp64.so)
==31069==    by 0x401163: MAIN__ (in /home/j/workbook/Test4/PROGRAM)
==31069==    by 0x401F09: main (in /home/j/workbook/Test4/PROGRAM)
==31069==  Address 0x6afe0b0 is 0 bytes inside a block of size 4 alloc'd
==31069==    at 0x4A069EE: malloc (vg_replace_malloc.c:270)
==31069==    by 0x400FE4: MAIN__ (in /home/j/workbook/Test4/PROGRAM)
==31069==    by 0x401F09: main (in /home/j/workbook/Test4/PROGRAM)
==31069==

I am not sure if the line "Invalid write of size 8" is referring to the size of the some of the integers or reals (as mentioned by Jonathan Dursi) or referring to the size of an array being passed into dsyevr.

Invalid write of size 8 means write of a single element of 8 bytes. It could be both integer or real. — Vladimir F Героям слава
– Vladimir F Героям слава, Commented Aug 7, 2014 at 15:47
I have tried setting all the parameters to kind=4 however when I do this I get an error messager: "MKL ERROR: Parameter 20 was incorrect on entry to DSYEVR" — James
– James, Commented Aug 7, 2014 at 16:43
SO the MKL requires all (?) integers to be 8 byte. However the gfortran compiler requires the integers to be long, i.e. 4 byte (as mentioned by Jonathan Dursi). Is this a correct. Is there a way around this problem? — James
– James, Commented Aug 7, 2014 at 16:53

Jonathan Dursi · Accepted Answer · 2014-08-07 01:09:00Z

4

It's not obvious from the documentation, but isuppz has to be allocated, even for that initial call where lwork = -1. If you move the isuppz allocation to before the first call your code completes successfully.

After that, since you're using the ILP (8-byte integer) version of MKL, e.g. LAPACK with 8-byte integer indices, all of your integer parameters to the LAPACK routines will have to be the same long kind, so that you'll need

integer(kind=8), allocatable :: isuppz(:), iwork(:)

I'll note that for neither and integer, kind=8 is actually part of the standard, and you should really use either selected_int/real_kind or the iso_fortran_env module and int64, etc.

edited Aug 7, 2014 at 1:09

answered Aug 6, 2014 at 16:26

Jonathan Dursi

51.1k10 gold badges131 silver badges160 bronze badges

Sign up to request clarification or add additional context in comments.

5 Comments

James Over a year ago

Thank you very much for your help. I have made this correction (and updated my question) but unfortunately I am still seeing the same error.

Jonathan Dursi Over a year ago

I can't reproduce the error you're getting now - works w/ gfortran (or ifort) + MKL for me.

James Over a year ago

That is odd,is that using the exact same code and makefile as my question now has? Would you mind telling me if your $MKLROOT is the same as mine? (I have seen this function work on ifort myself however it am trying to compile it with gfortran so that it can be compiled with a bunch of other codes that I don't think will compile quickly on ifort)

Jonathan Dursi Over a year ago

One remaining issue (which weirdly, didn't manifest with ifort?) -- all of your integers have to be of the long kind.

James Over a year ago

I added some output from valgrind to my question which may relate to this

Collectives™ on Stack Overflow

Problems with memory access/allocation when using LAPACK DSYEVR. Code used is FORTRAN

1 Answer 1

5 Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

5 Comments

Your Answer

Sign up or log in

Post as a guest

Related