maybe trying to simplify the program is wrong, here is the whole thing. ReconstructBlock which doesn't appear here doesn't seem to be the problem since commenting it out doesn't change anything and leaving only it solves the issue. The issue is with the call to EigenVecs. it computes a whole lot of eigenvectors of an operator using arpack. All the libraries used seem to be thread safe as far as i can tell. changing compiler doesn't help either. If I use more then one thread usually dsaupd reports memory corruption and ocationally either it or one of the deallocates segfaults. I also get the occasional abort signal. I tried making a minimal example but i can't reproduce the problem.

integer function PeronaMalikIteration(in, out, pprms, eprms)

real*8, dimension(r, c), intent(IN) :: in

real*8, dimension(r, c), intent(OUT) :: out

type(PeronaParams), intent(IN) :: pprms

type(EigParams), intent(IN) :: eprms

! Storage for the eigenvectors/eigenvalues

real*8, allocatable :: v(:,

, d(:)

!DEC$ ATTRIBUTES ALIGN: 16 :: v, d

! Parameters for window calculation

integer :: winLeft, winBot, colOffset, rowOffset, row, col

integer :: threadId, OMP_GET_THREAD_NUM

integer :: ret

!$OMP PARALLEL SHARED(in, out, eprms, pprms) PRIVATE(ret, winLeft, winBot, colOffset, rowOffset, col, row, v, d, threadId) NUM_THREADS(2)

allocate(v(eprms%ldv, eprms%ncv), d(eprms%neigs))

!$OMP DO SCHEDULE(DYNAMIC, 1)

! TODO: handle the last sub-block

do col = eprms%blockR + 1, eprms%c - eprms%blockR, eprms%blockSz

do row = eprms%blockR + 1, eprms%r - eprms%blockR, eprms%blockSz

threadId = OMP_GET_THREAD_NUM()

print *, "col: ", col, "row: ", row, "thread: ", threadId

! The window has to be centered around the pixel in question IF

! possible but still inside the boundaries of the image

winLeft = max(col - eprms%windowR, 1)

winLeft = min(winLeft, eprms%c - eprms%windowSz + 1)

colOffset = col - eprms%blockR - winLeft

winBot = max(row - eprms%windowR, 1)

winBot = min(winBot, eprms%r - eprms%windowSz + 1)

rowOffset = row - eprms%blockR - winBot

!PeronaMalikIteration =

ret = EigenVecs( &

in(winBot:winBot+eprms%windowSz-1, winLeft:winLeft+eprms%windowSz-1), &

v, d, pprms, eprms)

call ReconstructBlock( &

in(winBot:winBot+eprms%windowSz-1, winLeft:winLeft+eprms%windowSz-1), &

out(winBot+rowOffset:winBot+rowOffset+eprms%blockSz-1, &

winLeft+colOffset:winLeft+colOffset+eprms%blockSz-1), &

colOffset, rowOffset, v, d, pprms, eprms)

end do

end do

!$OMP END DO NOWAIT

deallocate(v, d)

!$OMP END PARALLEL

PeronaMalikIteration = 0

end function PeronaMalikIteration

! ------------------------------------------------------

! calculate eigenvectors

! ------------------------------------------------------

integer function EigenVecs(img, v, d, pprms, eprms)

type(PeronaParams), intent(IN) :: pprms

type(EigParams), intent(IN) :: eprms

real*8, intent(IN) :: img(eprms%windowSz, eprms%windowSz), &

v(eprms%ldv, eprms%ncv), d(eprms%neigs)

! Algorithm parameters

! Relative Tolerance, set to default - machine precision

real*8 :: tol = 0.0d0

real*8 :: sigma = 0.0d0

! Algorithm variables

integer :: ido, lworkl, info

integer, dimension(11) :: iparam, ipntr

real*8, allocatable, dimension(:) :: resid, workd, workl

!DEC$ ATTRIBUTES ALIGN: 16 :: resid, workd, workl

logical :: rvec = .true.

logical, allocatable, dimension(:) :: select

character(len=1) :: B = 'I'

character(len=2) :: which = 'LA'

! Set parameter values

! Default return, no error

EigenVecs = 0

! Initial parameters to the dsaupd routine

ido = 0

info = 0

! (/ ISHIFT, -, MXITER, NB, NCONV, -, MODE, NP, NUMOP, NUMOPB, NUMREO /)

iparam = (/ 1 , 0, 1000*eprms%neigs, 1 , 0 , 0, 1 , 0 , 0 , 0 , 0 /)

! Workspace

lworkl = eprms%ncv*(eprms%ncv + 8)

allocate(resid(eprms%n), workd(3*eprms%n), workl(lworkl))

do

call dsaupd(ido, B, eprms%n, which, eprms%neigs, tol, resid, &

eprms%ncv, v, eprms%ldv, iparam, ipntr, workd, workl, lworkl, info)

select case (ido)

case (-1,1)

! multiply y <- OP*x

! where x = workd(ipntr(1))

! y = workd(ipntr(2))

call PeronaMalikResponse(workd(ipntr(1)), workd(ipntr(2)), img, &

eprms, pprms)

case default

exit

end select

end do

! Error !!!

if (info /= 0) then

EigenVecs = info

goto 99999

end if

allocate(select(eprms%ncv))

call dseupd(rvec, 'A', select, d, v, eprms%ldv, sigma, B , eprms%n, &

which, eprms%neigs, tol, resid, eprms%ncv, v, eprms%ldv, iparam, &

ipntr, workd, workl, lworkl, info)

deallocate(select)

99999 deallocate(resid, workd, workl)

end function EigenVecs