Loop in fortran with openmp and allocatable arrays

Question

I like to do this:

program main
  implicit none
  integer l
  integer, allocatable, dimension(:) :: array

  allocate(array(10))

  array = 0

  !$omp parallel do private(array)
  do l = 1, 10
    array(l) = l
  enddo
  !$omp end parallel do

  print *, array

  deallocate(array)

end

But I am running into error messages:

* glibc detected * ./a.out: munmap_chunk(): invalid pointer: 0x00007fff25d05a40 *

This seems to be a bug in ifort according to some discussions at intel forums but should be resolved in the version I am using (11.1.073 - Linux). This is a MASSIVE downscaled version of my code! I unfortunately can not use static arrays to have a workaround.

If I put the print into the loop, I get other errors:

* glibc detected ./a.out: double free or corruption (out): 0x00002b22a0c016f0 **

Your code snippet is too incomplete and inconsistent with your description (only j and l are private according to the code snippet) to be useful. Chop your actual code down further to a small example that still exhibits the problem, and then show that, ideally in its entirety. While you're at it, make sure your threads have sufficient stack space and pop a DEFAULT(NONE) clause in on the OMP directive. Along the way chances are you'll figure out what's wrong... — IanH
@IanH I reformulated the question with better code. I also found a hint to where the error may be found. — DaPhil
What version of ifort are you using? Your code runs here with 13.1.0. — IanH
Does that version support OpenMP 3.0 (necessary for allocatables)? Regardless, might be time for an update. — IanH

mjswartz mjswartz · Accepted Answer · 2013-02-20T16:09:54

I just ran your code with ifort and openmp and it spewed 0d0's. I had to manually quit the execution. What is your expected output? I'm not a big fan of unnecessarily dynamically allocating arrays. You know what you're going to allocate your matrices as, so just make parameters and statically do it. I'll mess with some stuff and edit this response in a few.

Ok, so here's my edits:

  program main
  implicit none

  integer :: l, j
  integer, parameter :: lmax = 15e3
  integer, parameter :: jmax = 25
  integer, parameter :: nk = 300
  complex*16, dimension(9*nk) :: x0, xin, xout
  complex*16, dimension(lmax) :: e_pump, e_probe
  complex*16 :: e_pumphlp, e_probehlp
  character*25 :: problemtype
  real*8 :: m

  ! OpenMP variables
  integer :: myid, nthreads, omp_get_num_threads, omp_get_thread_num

  x0 = 0.0d0

  problemtype = 'type1'
  if (problemtype .ne. 'type1') then
     write(*,*) 'Problem type not specified. Quitting'
     stop
  else
     ! Spawn a parallel region explicitly scoping all variables
     !$omp parallel 
        myid = omp_get_thread_num()
        if (myid .eq. 0) then
           nthreads = omp_get_num_threads()
           write(*,*) 'Starting program with', nthreads, 'threads'
        endif

        !$omp do private(j,l,m,e_pumphlp,e_probehlp,e_pump,e_probe)
        do j = 1, jmax - 1
           do l = 1, lmax

              call electricfield(0.0d0, 0.0d0, e_pumphlp, &
                                 e_probehlp, 0.0d0)
              !   print *, e_pumphlp, e_probehlp

              e_pump(l) = e_pumphlp
              e_probe(l) = e_probehlp
              print *, e_pump(l), e_probe(l)

           end do
        end do
     !$omp end parallel
  end if

  end program main

Notice I removed your use of a module since it was unnecessary. You have an external module containing a subroutine, so just make it an external subroutine. Also, I changed your matrices to be statically allocated. Case statements are a fancy and expensive version of if statements. You were casing 15e3*25 times rather than once (expensive), so I moved those outside. I changed the OpenMP calls, but only semantically. I gave you some output so that you know what OpenMP is actually doing.

Here is the new subroutine:

  subroutine electricfield(t, tdelay, e_pump, e_probe, phase)
  implicit none

  real*8, intent(in) :: t, tdelay
  complex*16, intent(out) :: e_pump, e_probe
  real*8, optional, intent (in) :: phase

  e_pump = 0.0d0
  e_probe = 0.0d0

  return

  end subroutine electricfield

I just removed the module shell around it and changed some of your variable names. Fortran is not case sensitive, so don't torture yourself by doing caps and having to repeat it throughout.

I compiled this with

ifort -o diffeq diffeq.f90 electricfield.f90 -openmp

and ran with

./diffeq > output

to catch the program vomiting 0's and to see how many threads I was using:

(0.000000000000000E+000,0.000000000000000E+000) (0.000000000000000E+000,0.000000000000000E+000) (0.000000000000000E+000,0.000000000000000E+000) (0.000000000000000E+000,0.000000000000000E+000) (0.000000000000000E+000,0.000000000000000E+000) (0.000000000000000E+000,0.000000000000000E+000) (0.000000000000000E+000,0.000000000000000E+000) (0.000000000000000E+000,0.000000000000000E+000) (0.000000000000000E+000,0.000000000000000E+000) (0.000000000000000E+000,0.000000000000000E+000)
Starting program with 32 threads
(0.000000000000000E+000,0.000000000000000E+000) (0.000000000000000E+000,0.000000000000000E+000) (0.000000000000000E+000,0.000000000000000E+000) (0.000000000000000E+000,0.000000000000000E+000) (0.000000000000000E+000,0.000000000000000E+000) (0.000000000000000E+000,0.000000000000000E+000)

Hope this helps!

Loop in fortran with openmp and allocatable arrays

3 Answers