Почему мой код Fortran завернут в f2py, используя так много памяти?

Я пытаюсь рассчитать все расстояния между примерно сотней тысяч точек. У меня есть следующий код, написанный на Фортране и скомпилированный с использованием f2py:

C         1         2         3         4         5         6         7
C123456789012345678901234567890123456789012345678901234567890123456789012
        subroutine distances(coor,dist,n)
        double precision coor(n,3),dist(n,n)
        integer n
        double precision x1,y1,z1,x2,y2,z2,diff2

cf2py intent(in) :: coor,dist
cf2py intent(in,out):: dist
cf2py intent(hide)::n
cf2py intent(hide)::x1,y1,z1,x2,y2,z2,diff2

        do 200,i=1,n-1
            x1=coor(i,1)
            y1=coor(i,2)
            z1=coor(i,3)
            do 100,j=i+1,n
                x2=coor(j,1)
                y2=coor(j,2)
                z2=coor(j,3)
                diff2=(x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+(z1-z2)*(z1-z2)
                dist(i,j)=sqrt(diff2)
  100       continue
  200   continue

        end

Я компилирую код fortran, используя следующий код Python setup_collision.py:

# System imports
from distutils.core import *
from distutils      import sysconfig

# Third-party modules
import numpy
from numpy.distutils.core import Extension, setup

# Obtain the numpy include directory.  This logic works across numpy versions.
try:
    numpy_include = numpy.get_include()
except AttributeError:
    numpy_include = numpy.get_numpy_include()

# simple extension module
collision = Extension(name="collision",sources=['./collision.f'],
               include_dirs = [numpy_include],
               )

# NumyTypemapTests setup
setup(  name        = "COLLISION",
        description = "Module calculates collision energies",
        author      = "Stvn66",
        version     = "0.1",
        ext_modules = [collision]
        )

Затем запустите его следующим образом:

import numpy as np
import collision

coor = np.loadtxt('coordinates.txt')
n_atoms = len(coor)
dist = np.zeros((n_atoms, n_atoms), dtype=np.float16) # float16 reduces memory
n_dist = n_atoms*(n_atoms-1)/2
n_GB = n_dist * 2 / float(2**30) # 1 kB = 1024 B
n_Gb = n_dist * 2 / 1E9          # 1 kB = 1000 B
print 'calculating %d distances between %d atoms' % (n_dist, n_atoms)
print 'should use between %f and %f GB of memory' % (n_GB, n_Gb)
dist = collision.distances(coor, dist)

Используя этот код с 30 000 атомов, который должен использовать около 1 ГБ памяти для хранения расстояний, вместо этого он использует 10 ГБ. С этой разницей для выполнения этого расчета со 100 000 атомов потребуется 100 ГБ вместо 10 ГБ. У меня всего 20 Гб на компе.

Я упустил что-то, связанное с передачей данных между Python и Fortran? Огромная разница указывает на серьезный недостаток в реализации.


person Steven C. Howell    schedule 12.05.2015    source источник
comment
Вероятно, проблема с перетасовкой данных через границу между numpy и fortran? Однако я недостаточно знаю ни о том, ни о другом, чтобы строить догадки.   -  person Patrick Collins    schedule 12.05.2015
comment
Вам действительно нужно полное взаимодействие NxN? Вероятно, вы могли бы использовать схему разбиения (сетка, BSP, октодерево и т. д.), чтобы ускорить процесс...   -  person Alexander Vogt    schedule 12.05.2015
comment
Мне нужно знать ближайшие расстояния. Мне нужно определить, находятся ли какие-либо атомы на определенном расстоянии друг от друга, то есть сталкиваются (у меня есть аналогичная процедура, которая немедленно прерывается, если этот порог пересекается).   -  person Steven C. Howell    schedule 12.05.2015
comment
Взгляните на книгу Эриксона Обнаружение столкновений в реальном времени. Вы найдете несколько алгоритмов, которые значительно уменьшают сложность (у вас есть (On^2)!)   -  person Alexander Vogt    schedule 12.05.2015


Ответы (1)


Вы передаете подпрограмме Fortran массивы двойной точности. Каждый элемент двойной точности требует 8 байт памяти. Для N=30,000 это делает

coor(n,3) => 30,000*3*8 ~ 0.7 MB
dist(n,n) => 30,000^2*8 ~ 6.7 GB

Поскольку для Python дополнительно требуются числа с половинной точностью, это составляет еще 1-2 ГБ. Таким образом, общее требование составляет 9-10 ГБ.

То же самое верно и для N=100,000, для которого потребуется около 75 ГБ только для части Fortran.

Вместо double precision с плавающей запятой вы должны использовать одинарную точность reals - если этого достаточно для ваших расчетов. Это приведет к половине требований к памяти. [У меня нет опыта в этом, но я предполагаю, что если обе части используют одинаковую точность, Python может работать с данными напрямую...]

Как отметил @VladimirF в своем комментарии, "обычные компиляторы не поддерживают 2-байтовые вещественные числа". Я проверил с gfortran и ifort, и у них обоих нет. Поэтому вам нужно использовать как минимум одинарную точность.

person Alexander Vogt    schedule 12.05.2015
comment
Я думаю, что обычные компиляторы не поддерживают 2-байтовые числа. selected_real_kind вернет значение по умолчанию. - person Vladimir F; 12.05.2015
comment
Код, который я разместил, использовал половинную точность в питоне, но не в фортране (пропустил это). Я попытаюсь изменить это на половину или, по крайней мере, на один. - person Steven C. Howell; 12.05.2015
comment
Есть ли способ заставить Python и Fortran использовать одни и те же блоки памяти, чтобы избежать двойного требования? - person Steven C. Howell; 13.05.2015
comment
@stvn66: Да, вы можете использовать ctypes для этой задачи. См. мой ответ здесь. - person Alexander Vogt; 13.05.2015
comment
@ stvn66 Если вы создаете свои массивы numpy как order='F' (поэтому они хранятся так, как ожидает Fortran), f2py не должен их копировать. Я не уверен, как это взаимодействует с intent(in,out), поэтому вам, возможно, придется поэкспериментировать с этим, чтобы увидеть. - person DavidW; 15.05.2015