Сложность получения OpenMP для работы с f2py

Question

Я работаю над некоторыми симуляционными работами для своих исследований и столкнулся с тем, что импортировал fortran в мои сценарии python. В качестве фона я работал с Python уже несколько лет, и, когда возникла необходимость, они только играли вокруг Fortran.

Я проделал некоторую работу в прошлом, когда Fortran реализовал некоторые простые функции OpenMP. Я не являюсь экспертом в этом, но у меня уже есть основы.

Теперь я использую f2py для создания библиотеки, которую я могу вызвать из моего сценария python. Когда я пытаюсь скомпилировать openmp, он компилируется правильно и заканчивается, но есть нулевое улучшение скорости и, глядя вверх, я вижу, что использование ЦП указывает, что работает только один поток.

Я просмотрел документацию для f2py (которая не очень хорошо документирована), а также провела обычную проверку сети для ответов. Я включил код Fortran, который я компилирую, а также простой скрипт python, который вызывает его. Я также бросаю команду компиляции, которую я использую.

В настоящее время я срубил моделирование до 10^4 в качестве хорошего теста. В моей системе требуется 3 секунды. В конечном счете мне нужно запустить несколько 10^6 моделирования частиц, хотя, поэтому мне нужно немного сократить время.

Если кто-нибудь может указать мне в направлении того, как заставить мой код работать, это было бы оценено по достоинству. Я также могу попытаться включить любую подробную информацию о системе по мере необходимости.

Приветствия, Rylkan

---

1) обобщать

f2py -c --f90flags='-fopenmp' -lgomp -m calc_accel_jerk calc_accel_jerk.f90 

---

2) Python скрипт для вызова

import numpy as N 
import calc_accel_jerk 

# a is a (1e5,7) array with M,r,v information 
a  = N.load('../test.npy') 
a  = a[:1e4] 

out = calc_accel_jerk.calc(a,a.shape[0]) 
print out[:10] 

---

3) Fortran код

subroutine calc (input_array, nrow, output_array) 
implicit none 
!f2py threadsafe 
include "omp_lib.h" 

integer, intent(in) :: nrow 
double precision, dimension(nrow,7), intent(in) :: input_array 
double precision, dimension(nrow,2), intent(out) :: output_array 

! Calculation parameters with set values 
double precision,parameter :: psr_M=1.55*1.3267297e20 
double precision,parameter :: G_Msun=1.3267297e20 
double precision,parameter :: pc_to_m=3.08e16 

! Vector declarations 
integer :: irow 
double precision :: vfac 
double precision, dimension(nrow) :: drx,dry,drz,dvx,dvy,dvz,rmag,jfac,az,jz 

! Break up the input array for faster access 
double precision,dimension(nrow) :: input_M 
double precision,dimension(nrow) :: input_rx 
double precision,dimension(nrow) :: input_ry 
double precision,dimension(nrow) :: input_rz 
double precision,dimension(nrow) :: input_vx 
double precision,dimension(nrow) :: input_vy 
double precision,dimension(nrow) :: input_vz 

input_M(:) = input_array(:,1)*G_Msun 
input_rx(:) = input_array(:,2)*pc_to_m 
input_ry(:) = input_array(:,3)*pc_to_m 
input_rz(:) = input_array(:,4)*pc_to_m 
input_vx(:) = input_array(:,5)*1000 
input_vy(:) = input_array(:,6)*1000 
input_vz(:) = input_array(:,7)*1000 

!$OMP PARALLEL DO private(vfac,drx,dry,drz,dvx,dvy,dvz,rmag,jfac,az,jz) shared(output_array) NUM_THREADS(2) 
DO irow = 1,nrow 
    ! Get the i-th iteration 
    vfac = sqrt(input_M(irow)/psr_M) 
    drx = (input_rx-input_rx(irow)) 
    dry = (input_ry-input_ry(irow)) 
    drz = (input_rz-input_rz(irow)) 
    dvx = (input_vx-input_vx(irow)*vfac) 
    dvy = (input_vy-input_vy(irow)*vfac) 
    dvz = (input_vz-input_vz(irow)*vfac) 
    rmag = sqrt(drx**2+dry**2+drz**2) 
    jfac = -3*drz/(drx**2+dry**2+drz**2) 

    ! Calculate the acceleration and jerk 
    az = input_M*(drz/rmag**3) 
    jz = (input_M/rmag**3)*((dvx*drx*jfac)+(dvy*dry*jfac)+(dvz+dvz*drz*jfac)) 

    ! Remove bad index 
    az(irow) = 0 
    jz(irow) = 0 

    output_array(irow,1) = sum(az) 
    output_array(irow,2) = sum(jz) 
END DO 
!$OMP END PARALLEL DO 

END subroutine calc 

Answer 1

Вот простая проверка, чтобы увидеть, кастрированный баран нить OpenMP действительно видна в Fortran коды:

module OTmod 
    !$ use omp_lib 
    implicit none 

    public :: get_threads 

contains 

    function get_threads() result(nt) 
    integer :: nt 

    nt = 0 
    !$ nt = omp_get_max_threads() 

    end function get_threads 

end module OTmod 

Компиляция:

> f2py -m OTfor --fcompiler=gfortran --f90flags='-fopenmp' -lgomp -c OTmod.f90 

Execution :

> python 
>>> from OTfor import otmod 
>>> otmod.get_threads() 
12