Fortran: ошибка сегментации

Question

Я знаю, что когда-то я делал аналогичную тему, но это было другое. На этот раз добавление инструкции печати не изменится, получаю ли я segfault.

call omp_set_num_threads(omp_get_max_threads()) 
    !$omp parallel do & 
    !$omp default(firstprivate) & 
    !$omp private(present_surface) & 
    !$omp lastprivate(fermi) 
    do m = 1, fourth 
     do n = 1, third 
      do j = 1, second 
       do i = 1, first 
        !current angle is phi[i,j,ii,jj] 
        !we have to find the current fermi surface 
        present_surface = 0. 
        do a = 1, fourth  
         if (angle(i,j,n,m) == angle_array(a)) then 
          present_surface = surface(a) 
         end if 
        end do 
        if (radius(i,j,n,m) >= present_surface) then 
         fermi(i,j,n,m) = 0. 
        else 
         fermi(i,j,n,m) = 1. 
        end if   
       end do  
      end do  
     end do 
    end do 
    !$omp end parallel do 

Я не уверен в заявлении lastprivate (ферми), но в данный момент это не имеет значения. shared дает такое же поведение.

Итак, я запускаю этот скрипт с увеличением «первого, второго, третьего, четвертого». Типичный выход:

[10] Time elapsed is 0.001 [s]. 
[15] Time elapsed is 0.002 [s]. 
[20] Time elapsed is 0.005 [s]. 
./compile: line 2: 4310 Segmentation fault  python fortran_test.py 

Ну, как действовать. Я посмотрел на python gdb; запустить fortran_test.py, и обнаружили:

(gdb) run fortran_test.py 
Starting program: /usr/bin/python fortran_test.py 
[Thread debugging using libthread_db enabled] 
[New Thread 0x7fffed9b0700 (LWP 4251)] 
[New Thread 0x7fffe8ba5700 (LWP 4252)] 
[New Thread 0x7fffe83a4700 (LWP 4253)] 
[New Thread 0x7fffe7ba3700 (LWP 4254)] 
[10] Time elapsed is 0.008 [s]. 
[15] Time elapsed is 0.004 [s]. 
[20] Time elapsed is 0.005 [s]. 

Program received signal SIGSEGV, Segmentation fault. 
[Switching to Thread 0x7fffe7ba3700 (LWP 4254)] 
0x00007fffe8dc0bb7 in __populate_MOD_fermi_integrand._omp_fn.0() at  populate.f90:31 
31        do n = 1, third 

Если я что-то изменить во внутреннем цикле - например, удаление J, I петель и положить их на константу - тогда я просто получить Segfault на другой линии ,

I думаю, это имеет какое-то отношение к памяти, поскольку оно вызывает увеличение N. Тем не менее, то, что я пробовал (экспорт GOMP_STACKSIZE, OMP_STACKSIZE, ulimit), не исправил его, и я не вижу разницы в их использовании. (В настоящее время я удалил их).

Наконец, команда для компиляции (это в f2py):

f2py --fcompiler=gfortran --f90flags="-fopenmp -g" -lgomp -m -c populate populate.f90 

Как вы можете видеть, я совсем застрял. Я надеюсь, что некоторые из вас знают, как решить проблему.

Моя цель - быстро запустить N = 100 (отсюда функция openMP fortran), но это не должно иметь значения для кода fortran, не так ли?

Если вам интересно, 4 ГБ оперативной памяти и мой обмен 3.1G (Linux Swap/Solaris).

Спасибо!

Answer 1

default(firstprivate) lastprivate(fermi) означает, что каждый поток получает собственную копию angleradius, и fermi, с каждым существом размера first x second x third x fourth. По моему опыту, частные массивы всегда выделяются в стеке, даже когда предоставляется возможность автоматически превращать большие массивы в массивы кучи. Вероятно, ваш код сбой из-за недостаточного пространства стека.

Вместо того чтобы иметь все личное, вы должны действительно взглянуть на зависимость данных и выбрать правильные классы совместного использования данных. angle, angle_array, surface и radius никогда не записываются, поэтому все они должны быть shared. present_surface изменяется и должно быть private. fermi написано, но не в одном месте более чем с одним потоком, поэтому оно также должно быть shared. Также обратите внимание, что lastprivate не будет иметь ожидаемого результата, так как он дает доступное значение соответствующей переменной, имеющейся в потоке, который выполняет логически последнюю итерацию m -loop. first, second, third и fourth просто рассматриваются как константы и должны быть shared. Переменные цикла, конечно же, private, но это автоматизировано компилятором.

!$omp parallel do private(present_surface) 
do m = 1, fourth 
    do n = 1, third 
     do j = 1, second 
      do i = 1, first 
       !current angle is phi[i,j,ii,jj] 
       !we have to find the current fermi surface 
       present_surface = 0. 
       do a = 1, fourth  
        if (angle(i,j,n,m) == angle_array(a)) then 
         present_surface = surface(a) 
        end if 
       end do 
       if (radius(i,j,n,m) >= present_surface) then 
        fermi(i,j,n,m) = 0. 
       else 
        fermi(i,j,n,m) = 1. 
       end if   
      end do  
     end do  
    end do 
end do 
!$omp end parallel do