2015-02-07 3 views
7

я следующая функция:Сравнение производительности Fortran, Numpy, Cython и Numexpr

def get_denom(n_comp,qs,x,cp,cs): 
''' 
len(n_comp) = 1 # number of proteins 
len(cp) = n_comp # protein concentration 
len(qp) = n_comp # protein capacity 
len(x) = 3*n_comp + 1 # fit parameters 
len(cs) = 1 

''' 
    k = x[0:n_comp] 
    sigma = x[n_comp:2*n_comp] 
    z = x[2*n_comp:3*n_comp] 

    a = (sigma + z)*(k*(qs/cs)**(z-1))*cp 
    denom = np.sum(a) + cs 
    return denom 

Я сравниваю его против Fortran реализации (Моя первая Fortran функция когда-либо):

subroutine get_denom (qs,x,cp,cs,n_comp,denom) 

! Calculates the denominator in the SMA model (Brooks and Cramer 1992) 
! The function is called at a specific salt concentration and isotherm point 
! I loops over the number of components 

implicit none 

! declaration of input variables 
integer, intent(in) :: n_comp ! number of components 
double precision, intent(in) :: cs,qs ! salt concentration, free ligand concentration 
double precision, dimension(n_comp), INTENT(IN) ::cp ! protein concentration 
double precision, dimension(3*n_comp + 1), INTENT(IN) :: x ! parameters 

! declaration of local variables 
double precision, dimension(n_comp) :: k,sigma,z 
double precision :: a 
integer :: i 

! declaration of outpur variables 
double precision, intent(out) :: denom 

k = x(1:n_comp) ! equlibrium constant 
sigma = x(n_comp+1:2*n_comp) ! steric hindrance factor 
z = x(2*n_comp+1:3*n_comp) ! charge of protein 

a = 0. 
do i = 1,n_comp 
    a = a + (sigma(i) + z(i))*(k(i)*(qs/cs)**(z(i)-1.))*cp(i) 
end do 

denom = a + cs 

end subroutine get_denom 

Я скомпилировал. F95 файл с помощью:

1) f2py -c -m get_denom get_denom.f95 --fcompiler=gfortran

2) f2py -c -m get_denom_vec get_denom.f95 --fcompiler=gfortran --f90flags='-msse2' (последний вариант должен включить авто-векторизации)

Я проверить функции по:

import numpy as np 
import get_denom as fort_denom 
import get_denom_vec as fort_denom_vec 
from matplotlib import pyplot as plt 
%matplotlib inline 

def get_denom(n_comp,qs,x,cp,cs): 
    k = x[0:n_comp] 
    sigma = x[n_comp:2*n_comp] 
    z = x[2*n_comp:3*n_comp] 
    # calculates the denominator in Equ 14a - 14c (Brooks & Cramer 1992) 
    a = (sigma + z)*(k*(qs/cs)**(z-1))*cp 
    denom = np.sum(a) + cs 
    return denom 

n_comp = 100 
cp = np.tile(1.243,n_comp) 
cs = 100. 
qs = np.tile(1100.,n_comp) 
x= np.random.rand(3*n_comp+1) 
denom = np.empty(1) 
%timeit get_denom(n_comp,qs,x,cp,cs) 
%timeit fort_denom.get_denom(qs,x,cp,cs,n_comp) 
%timeit fort_denom_vec.get_denom(qs,x,cp,cs,n_comp) 

Я добавил следующий Cython код:

import cython 
# import both numpy and the Cython declarations for numpy 
import numpy as np 
cimport numpy as np 

@cython.boundscheck(False) 
@cython.wraparound(False) 
def get_denom(int n_comp,np.ndarray[double, ndim=1, mode='c'] qs, np.ndarray[double, ndim=1, mode='c'] x,np.ndarray[double, ndim=1, mode='c'] cp, double cs): 

    cdef int i 
    cdef double a 
    cdef double denom 
    cdef double[:] k = x[0:n_comp] 
    cdef double[:] sigma = x[n_comp:2*n_comp] 
    cdef double[:] z = x[2*n_comp:3*n_comp] 
    # calculates the denominator in Equ 14a - 14c (Brooks & Cramer 1992) 
    a = 0. 
    for i in range(n_comp): 
    #a += (sigma[i] + z[i])*(pow(k[i]*(qs[i]/cs), (z[i]-1)))*cp[i] 
     a += (sigma[i] + z[i])*(k[i]*(qs[i]/cs)**(z[i]-1))*cp[i] 

    denom = a + cs 

    return denom 

EDIT:

Добавлено Numexpr , используя одну резьбу:

def get_denom_numexp(n_comp,qs,x,cp,cs): 
    k = x[0:n_comp] 
    sigma = x[n_comp:2*n_comp] 
    z = x[2*n_comp:3*n_comp] 
    # calculates the denominator in Equ 14a - 14c (Brooks & Cramer 1992) 
    a = ne.evaluate('(sigma + z)*(k*(qs/cs)**(z-1))*cp') 
    return cs + np.sum(a) 

ne.set_num_threads(1) # using just 1 thread 
%timeit get_denom_numexp(n_comp,qs,x,cp,cs) 

Результат (меньше, тем лучше):

enter image description here

Почему это скорость Fortran ближе к Numpy с увеличением размера массивов? И как я мог ускорить Cython? Использование указателей?

+13

Я не смотрел ваш код в деталях. Но, как правило, подобное наблюдение связано с тем, что оба метода имеют * сравнимую скорость обработки, как только данные находятся в нужном месте в памяти, в правильном формате. Однако получить его там требуется другое время. Обычно это называется «накладными расходами». Таким образом, возможно, накладные расходы на изготовление больше в решении numpy, чем в решении fortran. Эти различия становятся все менее значительными с увеличением размера полезной нагрузки. –

+0

Опять же, обратите внимание на некоторые накладные расходы в коде. Кажется, что существует некоторое постоянное смещение, требующее значительного времени для Cython. –

+0

Возможно, также рассмотрим 'numexpr', особенно если у вас есть MKL, и ваши массивы большие. Это проще в использовании, чем Cython и f2py, и вы получаете многопоточность бесплатно. –

ответ

3

Sussed It.

ОК, наконец, нам разрешили установить Numpy и т. Д. В одном из наших ящиков, и это позволило полностью объяснить ваше первоначальное сообщение.

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

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

Часть A: Ответ на OP

Первых вещей первый, имеет дело с вопросом в исходном сообщении: Вы можете вспомнить, я мог бы прокомментировать только WRT на Fortran стороны, так как наша политика строга о том, что программное обеспечении могут быть установлены и где на наших машинах, и у нас не было Python и т. д. (до сих пор). Я также неоднократно заявлял, что характер вашего результата был интересен с точки зрения того, что мы можем назвать его изогнутым характером или, возможно, «вогнутостью».

Кроме того, работая исключительно с «относительными» результатами (поскольку вы не публиковали абсолютные тайминги, и в то время у меня не было Numpy), я несколько раз указывал, что может быть скрыта какая-то важная информация в нем.

Именно так.

Во-первых, я хотел быть уверенным, что могу воспроизвести ваши результаты, поскольку мы не используем Python/F2py нормально, неясно, какие параметры компилятора и т. Д. Подразумеваются в ваших результатах, поэтому я выполнил множество тестов для убедитесь, что мы говорим о яблоках-яблоках (мой апрельский ответ показал, что Debug vs Release/O2 имеет большое значение).

На рисунке 1 показаны результаты моего Python только для трех случаев: исходная внутренняя подпрограмма Python/Numpy (назовите это P, я просто вырезал/вставил ваш оригинал), ваш оригинал Do based Fortran s/r, который у вас был (назовите этот FDo, я как раз копировал/вставлял ваш оригинал, как и раньше), и один из вариантов, которые я предположил ранее, полагаясь на разделы массива и, таким образом, требуя Sum() (назовите это FSAS, созданный путем редактирования исходного FDo) , На рисунке 1 показаны абсолютные тайминги через timeit. Figure 1

На рисунке 2 показаны относительные результаты, основанные на вашей относительной стратегии деления на тайминги Python/Numpy (P). Отображаются только два (относительных) варианта Fortran. Figure 2

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

Теперь ваш первоначальный вопрос: «Почему скорость Fortran приближается к Numpy с увеличением размера массивов?».

Фактически, это не так. Это чисто артефакт/искажение, основанное исключительно на «относительных» таймингах, которые могут дать такое впечатление.

Посмотрите на рисунок 1, с абсолютными таймингами, ясно, что Numpy и Fortran расходятся. Итак, на самом деле, результаты Fortran отходят от Numpy или наоборот, если хотите.

Лучший вопрос, который неоднократно возникал в моих предыдущих комментариях, является причиной того, почему эти восходящие изгибы в первую очередь (например, линейные)? Мои предыдущие результаты, полученные только в Fortran, показали коэффициент относительной производительности «в основном» (желтые строки в диаграмме/ответе на апрель 14 апреля), и поэтому я предположил, что на стороне Python было что-то интересное, и предложил проверить это.

Один из способов показать это с еще одним видом относительной меры. Я разделил каждую (абсолютную) серию с ее собственным самым высоким значением (т. Е. На n_comp = 10k), чтобы увидеть, как разворачивается эта «внутренняя относительная» производительность (они называются значениями 10k, представляющими тайминги для n_comp = 10 000) , На рисунке 3 показаны эти результаты для P, FDo и FSAS как P/P10k, FDo/FDo10k и FSAS/FSAS10k. Для ясности ось y имеет логарифмический масштаб. Понятно, что варианты Fortran заготовки относительно намного лучше с уменьшением n_comp c.f. результаты Р (например, аннотированный красный круг). Figure 3

По-другому, Fortran очень (нелинейно) более эффективен для уменьшения размера массива. Непонятно, почему Python будет намного хуже с уменьшением n_comp ... но там он и может быть проблемой с внутренними накладными расходами/настройкой и т. Д., А также различиями между интерпретаторами и компиляторами и т. Д.

Так что это не то, что Fortran «догоняет» Python, совершенно наоборот, он продолжает дистанцироваться от Python (см. Рис. 1). Однако различия в нелинейностях между Python и Fortran с уменьшением n_comp производят «относительные» тайминги с явно противоречивым и нелинейным характером.

Таким образом, поскольку n_comp увеличивается, и каждый метод «стабилизируется» в более или менее линейном режиме, кривые выравниваются, чтобы показать, что их различия растут линейно, а относительные отношения согласуются с приблизительной константой (игнорируя конфликт памяти, smp проблемы и т. д.) ... это проще, если разрешить n_comp> 10k, но желтая строка в ответе 14 апреля уже показывает это для только для файлов Fortran.

Кроме этого: Мое предпочтение заключается в создании моих собственных подпрограмм/функций времени. timeit кажется удобным, но внутри этого «черного ящика» многое происходит. Настройка собственных циклов и структур, а также наличие свойств/разрешения ваших функций синхронизации важны для правильной оценки.

+0

Я не могу комментировать дискуссии между Владимиром и вами, но мне нравится ваш подход, и ответ легко понять. Поэтому я отмечаю это как принятый ответ. – Moritz

+2

Поскольку ситуация слишком личная, и она перешла в дискуссию, и это не дискуссионный сервер, я просто добавлю техническое замечание о том, как работает SO. Это не дискуссионный сервер. Предпочтительным здесь является тот факт, что я добавил свой ответ вместо того, чтобы начать другой. Не рекомендуется начинать новые ответы здесь, когда старый может быть отредактирован. Я даже рекомендовал бы @ user3024046 удалить старые и сделать один окончательный ответ и отполировать его до формы его предпочтения. –

+1

Я просто добавлю одну ссылку, она написана одним разработчиком gfortran, который, к сожалению, покинул это сообщество. Я знаю его имя, это было в его ник, прежде чем он ушел. Вы можете многому научиться. http://stackoverflow.com/a/29104905/721644 Я фактически сторонник предполагаемых массивов размеров, и я использую их везде, где могу, я просто использую атрибут 'прилежащий', когда это необходимо, чтобы преодолеть возможные накладные расходы. Я не использую их для повышения эффективности, потому что они не могут этого сделать.Только вы приняли разъяснительный комментарий как личное нападение. –

-2

Существует не достаточно информации в примечаниях, но некоторые из следующих может помочь:

1) Fortran оптимизировал внутренние функции, такие как «Сумма()» и «Dot_Product», который вы можете захотеть использовать вместо цикла Do для суммирования и т. д.

В некоторых случаях (не обязательно для обработки) здесь может быть «лучше» использовать ForAll или что-то другое для создания массивов «meta», которые должны быть суммированы, а затем применить суммирование на «мета» массивах.

Однако Fortran разрешает разделы массива, поэтому вам не нужно создавать автоматические/промежуточные массивы sigma, k и z и освобожденные служебные данные. Вместо этого может быть что-то вроде

n_compP1 = n_comp+1 
n_compT2 = n_comp*2 
a = Sum(x(1:n_comp)+2*x(n_compP1,n_compT2)) ! ... just for example 

2) Иногда (в зависимости от компилятора, машины и т.д.), может быть «раздоры памяти», если размеры массива не в определенных «двоичных интервалов» (например, 1024 против 1000) и т. д.

Возможно, вы захотите повторить свои эксперименты еще в нескольких точках вашей диаграммы (т. е. в разных других «n_comps») и особенно вблизи таких «границ».

3) Не удается определить, используете ли вы полную оптимизацию компилятора (флаги) для компиляции кода fortran. Вы можете посмотреть различные флажки «-o» и т. Д.

4) Возможно, вы захотите включить директиву OpemMP (или, по крайней мере, включить openmp в свои флаги и т. Д.). Иногда это может улучшить некоторые служебные проблемы, даже если они явно не полагаются на директивы OpenMP в ваших циклах и т. Д.

5) Общие сведения: Это, скорее всего, относится к каждому из ваших методов, где петли используются

а) «постоянные операции» в «суммировании» формула может быть выполнены вне цикла, например. создайте что-то вроде qsDcs = qs/cs и используйте qsDcs в цикле.

b) Точно так же иногда полезно создавать что-то вроде zM1 (:) = z (:) - 1 и вместо этого использовать zM1 (:) в цикле.

+0

1) сомнительно, компиляторы Fortran преуспевают в оптимизации циклов 'DO' и имеют реальные проблемы с' FORALL' и некоторыми назначениями массива. 3) это '-O' не' -o' 4), что вы имеете в виду? –

+1

На самом деле это не отвечает на вопрос, вы предложили несколько общих советов по оптимизации. Вопрос был: * Почему скорость Fortran приближается к Numpy с увеличением размера массивов? И как я мог ускорить Cython? С помощью указателей? * –

+0

Дорогой Владимир: По первому вопросу: Да, Фортран хорош в деле. Однако, поскольку пользователь не уверен, почему он наблюдает за поведением, которое он видит, кажется разумным протестировать другие варианты реализации Fortran, чтобы проверить, остается ли сравнение с другими методами как в оригинале. Повторите свой второй пункт: разумно попросить больше информации, чтобы увидеть, действительно ли изменение с размером массива подразумевается гигантским скачком от 1000 до 10000 ... мы не можем быть уверены, что поведение так же просто, как подразумевается. – DrOli

0

В дополнение к моему предыдущему ответу и слабым спекуляциям Владимира я установил два s/r: один как исходный, и один с использованием разделов массива и Sum(). Я также хотел продемонстрировать, что замечания Владимира о оптимизации цикла Do являются слабыми.

Кроме того, точка, которую я обычно делаю для бенчмаркинга, размер n_comp в приведенном выше примере слишком мал. В приведенных ниже результатах каждый из вариантов «оригинал» и «лучший» SumArraySection (SAS) повторяется в циклах, повторяющихся 1000 раз внутри вызовов синхронизации, поэтому результаты для 1000 вычислений каждого s/r. Если ваши тайминги составляют доли секунды, они, вероятно, ненадежны.

Существует ряд вариантов, которые стоит учитывать, но нет явных указателей. Один вариант используется для этих иллюстраций

subroutine get_denomSAS (qs,x,cp,cs,n_comp,denom) 

! Calculates the denominator in the SMA model (Brooks and Cramer 1992) 
! The function is called at a specific salt concentration and isotherm point 
! I loops over the number of components 

implicit none 

! declaration of input variables 
integer, intent(in) :: n_comp ! number of components 
double precision, intent(in) :: cs,qs ! salt concentration, free ligand concentration 
double precision, Intent(In)   :: cp(:) 
double precision, Intent(In)   :: x(:) 

! declaration of local variables 
integer :: i 

! declaration of outpur variables 
double precision, intent(out) :: denom 
! 
! 
double precision      :: qsDcs 
! 
! 
qsDcs = qs/cs 
! 
denom = Sum((x(n_comp+1:2*n_comp) + x(2*n_comp+1:3*n_comp))*(x(1:n_comp)*(qsDcs) & 
              **(x(2*n_comp+1:3*n_comp)-1))*cp(1:n_comp)) + cs 
! 
! 
end subroutine get_denomSAS 

Основные отличия:

а) прошли массивы (:) б) нет массива заданий в S/R, вместо того, чтобы использовать секции массива (что эквивалентно «эффективные» указатели). c) Используйте Sum() вместо Do

Затем попробуйте две различные оптимизации компилятора, чтобы продемонстрировать последствия.

Как показывают две диаграммы, код оригинала (синие бриллианты) намного медленнее c.f. SAS (красные квадраты) с низкой оптимизацией. SAS все еще лучше с высокой оптимизацией, но они приближаются. Частично это объясняется тем, что Sum() «оптимизирован», когда используется низкая оптимизация компилятора.

enter image description here

Желтые линии показывают соотношение между этими двумя тайминги S/R в. Игнорируйте значение желтой строки в «0» в верхнем изображении (слишком маленький n_comp вызвал одно из таймингов, которые нужно уйти)

Поскольку у меня нет исходных данных пользователя в отношении к Numpy, я могу сделать только что кривая SAS на его диаграмме должна лежать ниже его текущих результатов Фортрана и, возможно, быть более плоской или даже вниз по тренду.

Иными словами, на самом деле не может существовать расхождение, наблюдаемое в оригинальной публикации, или, по крайней мере, не в такой степени.

... хотя больше экспериментов может быть полезно, чтобы продемонстрировать и другие комментарии/ответы, уже предоставленные.

Уважаемый Мориц: Ой, я забыл упомянуть и относясь к вашему вопросу о указателях. Как и раньше, ключевой причиной улучшения с изменением SAS является то, что он лучше использует «эффективные указатели» в том смысле, что он устраняет необходимость переназначить массив x() на три новых локальных массива (т. Е. Поскольку x передается по ref , использование секций массива - это своего рода подход указателя, встроенный в Fortran, и, следовательно, нет необходимости в явных указателях), но затем требуется Sum() или Dot_Product() или что-то еще.

Вместо этого вы можете сохранить Do и достичь чего-то подобного, изменив x либо на 2D-массив n_compx3, либо напрямую передайте три явных массива 1D порядка n_comp. Это решение, скорее всего, будет определяться размером и сложностью вашего кода, поскольку для этого потребуется переписать операторы call/sr и т. Д., А в другом месте используется x(). Некоторые из наших проектов составляют> 300 000 строк кода, поэтому в этом случае гораздо дешевле менять код локально, например, на SAS и т. Д.

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

Конечно, комментарии о «правильном» бенчмаркинге и т. Д., А также вопрос о том, какие компиляторные переключатели подразумеваются по вашему использованию fpy, по-прежнему применяются, поскольку они могут значительно изменить характер ваших результатов.

Мне было бы интересно увидеть ваши результаты, если они были обновлены для этих перестановок.

+0

Если вы заинтересованы, я могу вставить код в pastebin или аналогичный. Спасибо всем за комментарии и ответы, я многому учусь. – Moritz

+0

Спасибо за использование моего имени в вашем ответе. Если вы хотите опровергнуть мои комментарии, вам нужно протестировать версию, которая использует do loops вместо разделов массивов. Я этого не вижу. Я не имею в виду оригинал. Я имею в виду, что выражение Sum было перезаписано для создания циклов. –

+0

BTW O2 - это не слишком высокая оптимизация. Направьте для -O5 или -Ofast, вы хотите увидеть векторизацию и другие вещи, которые нужно использовать. –

2

Будучи названным в другом ответе, я должен ответить.

Я знаю, что это действительно не отвечает на исходный вопрос, но оригинальный плакат поощрял преследовать это направление в своих комментариях.

Мои пункты таковы:

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

Я перевел предлагаемые массивы SAS в обычную петлю. Я называю это DOS. Я демонстрирую, что петли DO ничуть не медленнее, обе подпрограммы приводят к более или менее одному и тому же коду в этом случае.

qsDcs = qs/cs 

denom = 0 
do j = 1, n_comp 
    denom = denom + (x(n_comp+j) + x(2*n_comp+j)) * (x(j)*(qsDcs)**(x(2*n_comp+j)-1))*cp(j) 
end do 

denom = denom + cs 

Важно сказать, что я не верю, что эта версия менее читаемая только потому, что он имеет еще одну или две строки. На самом деле довольно просто понять, что здесь происходит.

Теперь тайминги для этих

f2py -c -m sas sas.f90 --opt='-Ofast' 
f2py -c -m dos dos.f90 --opt='-Ofast' 


In [24]: %timeit test_sas(10000) 
1000 loops, best of 3: 796 µs per loop 

In [25]: %timeit test_sas(10000) 
1000 loops, best of 3: 793 µs per loop 

In [26]: %timeit test_dos(10000) 
1000 loops, best of 3: 795 µs per loop 

In [27]: %timeit test_dos(10000) 
1000 loops, best of 3: 797 µs per loop 

Они так же. В массиве intrinsics и арифметике выражения массива нет скрытой магии производительности. В этом случае они просто переводятся на петли под капотом.

Если вы проверяете сгенерированный код GIMPLE, как SAS и DOS переводятся на тот же основной блок оптимизированного кода, не магическая версия SUM() не называется здесь:

<bb 8>: 
    # val.8_59 = PHI <val.8_49(9), 0.0(7)> 
    # ivtmp.18_123 = PHI <ivtmp.18_122(9), 0(7)> 
    # ivtmp.25_121 = PHI <ivtmp.25_120(9), ivtmp.25_117(7)> 
    # ivtmp.28_116 = PHI <ivtmp.28_115(9), ivtmp.28_112(7)> 
    _111 = (void *) ivtmp.25_121; 
    _32 = MEM[base: _111, index: _106, step: 8, offset: 0B]; 
    _36 = MEM[base: _111, index: _99, step: 8, offset: 0B]; 
    _37 = _36 + _32; 
    _40 = MEM[base: _111, offset: 0B]; 
    _41 = _36 - 1.0e+0; 
    _42 = __builtin_pow (qsdcs_18, _41); 
    _97 = (void *) ivtmp.28_116; 
    _47 = MEM[base: _97, offset: 0B]; 
    _43 = _40 * _47; 
    _44 = _43 * _42; 
    _48 = _44 * _37; 
    val.8_49 = val.8_59 + _48; 
    ivtmp.18_122 = ivtmp.18_123 + 1; 
    ivtmp.25_120 = ivtmp.25_121 + _118; 
    ivtmp.28_115 = ivtmp.28_116 + _113; 
    if (ivtmp.18_122 == _96) 
    goto <bb 10>; 
    else 
    goto <bb 9>; 

    <bb 9>: 
    goto <bb 8>; 

    <bb 10>: 
    # val.8_13 = PHI <val.8_49(8), 0.0(6)> 
    _51 = val.8_13 + _17; 
    *denom_52(D) = _51; 

код функционально идентичен версия do loop, просто имя переменных отличается.


2. Они предположили, что аргументы массива формы (:) имеют потенциал, чтобы снизить производительность.Принимая во внимание, что аргумент, полученный в аргументе предполагаемого размера (*) или явный размер (n), всегда просто смежный, предполагаемая форма теоретически не обязательно должна быть, и компилятор должен быть подготовлен к этому. Поэтому я всегда рекомендую использовать атрибут contiguous для ваших предполагаемых аргументов формы везде, где вы знаете, что вы назовете его непрерывными массивами.

В другом ответе изменение было совершенно бессмысленным, поскольку оно не использовало никаких преимуществ предполагаемых аргументов формы. А именно, что вам не нужно передавать аргументы с размерами массива, и вы можете использовать встроенные функции, такие как size() и shape().


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

import numpy as np 
import dos as dos 
import sas as sas 
from matplotlib import pyplot as plt 
import timeit 
import numexpr as ne 

#%matplotlib inline 



ne.set_num_threads(1) 

def test_n(n_comp): 

    cp = np.tile(1.243,n_comp) 
    cs = 100. 
    qs = np.tile(1100.,n_comp) 
    x= np.random.rand(3*n_comp+1) 

    def test_dos(): 
     denom = np.empty(1) 
     dos.get_denomsas(qs,x,cp,cs,n_comp) 


    def test_sas(): 
     denom = np.empty(1) 
     sas.get_denomsas(qs,x,cp,cs,n_comp) 

    def get_denom(): 
     k = x[0:n_comp] 
     sigma = x[n_comp:2*n_comp] 
     z = x[2*n_comp:3*n_comp] 
     # calculates the denominator in Equ 14a - 14c (Brooks & Cramer 1992) 
     a = (sigma + z)*(k*(qs/cs)**(z-1))*cp 
     denom = np.sum(a) + cs 
     return denom 

    def get_denom_numexp(): 
     k = x[0:n_comp] 
     sigma = x[n_comp:2*n_comp] 
     z = x[2*n_comp:3*n_comp] 
     loc_cp = cp 
     loc_cs = cs 
     loc_qs = qs 
     # calculates the denominator in Equ 14a - 14c (Brooks & Cramer 1992) 
     a = ne.evaluate('(sigma + z)*(k*(loc_qs/loc_cs)**(z-1))*loc_cp') 
     return cs + np.sum(a) 

    print 'py', timeit.Timer(get_denom).timeit(1000000/n_comp) 
    print 'dos', timeit.Timer(test_dos).timeit(1000000/n_comp) 
    print 'sas', timeit.Timer(test_sas).timeit(1000000/n_comp) 
    print 'ne', timeit.Timer(get_denom_numexp).timeit(1000000/n_comp) 


def test(): 
    for n in [10,100,1000,10000,100000,1000000]: 
     print "-----" 
     print n 
     test_n(n) 

Результаты:

  py    dos    sas    numexpr 
10   11.2188110352 1.8704519272 1.8659651279 28.6881871223 
100   1.6688809395 0.6675260067 0.667083025  3.4943861961 
1000  0.7014708519 0.5406000614 0.5441288948 0.9069931507 
10000  0.5825948715 0.5269498825 0.5309231281 0.6178650856 
100000  0.5736029148 0.526198864  0.5304090977 0.5886831284 
1000000  0.6355218887 0.5294830799 0.5366530418 0.5983200073 
10000000 0.7903120518 0.5301260948 0.5367569923 0.6030929089 

speed comparison

Вы можете видеть, что существует очень небольшая разница между этими двумя версиями Фортрана. Синтаксис массива немного медленнее, но на самом деле ничего не говорится.

Заключение 1: В этом сравнении накладных расходах для всех должно быть справедливым, и вы увидите, что для идеальной длиной вектора Numpy и Numexpr CAN почти достигает производительность FORtran, но когда вектор слишком мал, или, возможно, даже слишком велико преобладают решения Python.

Заключение 2: версия SAS с более высокой производительностью в другом сравнении обусловлена ​​сравнением с оригинальной версией OP, которая не является эквивалентной. Эквивалентная оптимизированная версия цикла DO включена в мой ответ.