Рекомендации по производительности

В следующих разделах мы кратко рассмотрим несколько техник, которые могут помочь сделать ваш код Julia максимально быстрым.

Избегайте глобальных переменных

Глобальная переменная может изменить свое значение, а следовательно, и свой тип в любой момент. Это затрудняет компилятору оптимизацию кода, использующего глобальные переменные. Переменные должны быть локальными или передаваться в качестве аргументов функциям, когда это возможно.

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

Мы обнаружили, что глобальные имена часто являются константами, и объявление их как констант значительно улучшает производительность:

const DEFAULT_VAL = 0

Использование неконстантных глобальных переменных можно оптимизировать, добавив аннотации типов в месте использования:

global x
y = f(x::Int + 1)

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

ПРИМЕЧАНИЕ: Весь код в REPL вычисляется в глобальной области видимости, поэтому переменная, определенная и присвоенная на верхнем уровне, будет глобальной переменной.

В следующей сессии REPL:

julia> x = 1.0

эквивалентно:

julia> global x = 1.0

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

Измеряйте производительность с помощью @time и обратите внимание на выделение памяти

Самый полезный инструмент для измерения производительности — макрос @time. Следующий пример демонстрирует хороший стиль работы:

julia> function f(n)
           s = 0
           for i = 1:n
               s += i/2
           end
           s
       end
f (generic function with 1 method)

julia> @time f(1)
elapsed time: 0.004710563 seconds (93504 bytes allocated)
0.5

julia> @time f(10^6)
elapsed time: 0.04123202 seconds (32002136 bytes allocated)
2.5000025e11

При первом вызове (@time f(1)), f компилируется. (Если вы еще не использовали @time в этой сессии, он также скомпилирует функции, необходимые для измерения времени.) Результаты этого выполнения не следует рассматривать всерьёз. При втором выполнении, помимо отчёта о времени, также будет указано, что было выделено большое количество памяти. Это единственное преимущество @time по сравнению с функциями, такими как tic() и toc(), которые сообщают только о времени.

Неожиданное выделение памяти почти всегда является признаком какой-либо проблемы с вашим кодом, обычно проблемы со стабильностью типов. Соответственно, помимо самого выделения, очень вероятно, что сгенерированный для вашей функции код далёк от оптимального. Возьмите такие указания всерьёз и следуйте приведенным ниже рекомендациям.

В качестве затравки обратите внимание, что улучшенная версия этой функции не выделяет памяти (кроме передачи результата обратно в REPL) и имеет на порядок более высокую скорость выполнения после первого вызова:

julia> @time f_improved(1)   # first call
elapsed time: 0.003702172 seconds (78944 bytes allocated)
0.5

julia> @time f_improved(10^6)
elapsed time: 0.004313644 seconds (112 bytes allocated)
2.5000025e11

Ниже вы узнаете, как определить проблему с f и как ее исправить.

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

Инструменты

Julia и ее экосистема пакетов включают инструменты, которые могут помочь вам диагностировать проблемы и улучшить производительность вашего кода:

• Профилирование позволяет измерить производительность вашего работающего кода и определить строки, являющиеся узкими местами. Для сложных проектов пакет ProfileView может помочь визуализировать результаты профилирования.
• Неожиданно большие выделения памяти, как сообщалось в @time, @allocated или в профилировщике (через вызовы процедур сборки мусора), указывают на то, что могут быть проблемы с вашим кодом. Если вы не видите другой причины для выделения памяти, подозревайте проблему с типом. Вы также можете запустить Julia с опцией --track-allocation=user и изучить полученные *.mem файлы, чтобы получить информацию о том, где произошли эти выделения. См. Анализ выделения памяти.
• @code_warntype генерирует представление вашего кода, которое может быть полезным для поиска выражений, которые приводят к неопределённости типа. См. @code_warntype ниже.
• Пакеты Lint и TypeCheck также могут предупредить вас об определенных типах программистских ошибок.

Избегайте контейнеров с параметрами абстрактного типа

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

Рассмотрим следующее:

a = Real[]    # typeof(a) = Array{Real,1}
if (f = rand()) < .8
    push!(a, f)
end

Поскольку a является массивом абстрактного типа Real, он должен уметь хранить любое вещественное значение. Поскольку объекты Real могут иметь произвольный размер и структуру, a должен быть представлен массивом указателей на индивидуально выделенные объекты Real. Поскольку f всегда будет Float64, мы должны вместо этого использовать:

a = Float64[] # typeof(a) = Array{Float64,1}

который создаст непрерывный блок значений с плавающей запятой 64-битной точности, которые можно эффективно обрабатывать.

См. также обсуждение в Параметризованные типы.

Объявления типов

Во многих языках с необязательными объявлениями типов добавление объявлений — основной способ сделать код быстрее. Это не так в Julia. В Julia компилятор, как правило, знает типы всех аргументов функций, локальных переменных и выражений. Однако есть несколько конкретных случаев, когда объявления полезны.

Избегайте полей с абстрактным типом

Типы могут быть объявлены без указания типов их полей:

julia> type MyAmbiguousType
           a
       end

Это позволяет a быть любого типа. Это часто бывает полезным, но у этого есть недостаток: для объектов типа MyAmbiguousType, компилятор не сможет сгенерировать высокопроизводительный код. Причина в том, что компилятор использует типы объектов, а не их значения, для определения того, как строить код. К сожалению, очень мало можно вывести о объекте типа MyAmbiguousType:

julia> b = MyAmbiguousType("Hello")
MyAmbiguousType("Hello")

julia> c = MyAmbiguousType(17)
MyAmbiguousType(17)

julia> typeof(b)
MyAmbiguousType

julia> typeof(c)
MyAmbiguousType

b и c имеют один и тот же тип, но их внутреннее представление данных в памяти очень отличается. Даже если вы храните только числовые значения в поле a, тот факт, что представление в памяти UInt8 отличается от Float64 также означает, что процессору нужно обрабатывать их с помощью двух различных типов инструкций. Поскольку необходимая информация недоступна в типе, такие решения должны приниматься во время выполнения. Это снижает производительность.

Вы можете добиться лучшего результата, объявив тип a. Здесь мы сосредоточены на случае, когда a может быть любым из нескольких типов, в этом случае естественным решением является использование параметров. Например:

julia> type MyType{T<:AbstractFloat}
         a::T
       end

Это лучший выбор, чем

julia> type MyStillAmbiguousType
         a::AbstractFloat
       end

потому что в первом варианте тип a определяется типом объекта-обертки. Например:

julia> m = MyType(3.2)
MyType{Float64}(3.2)

julia> t = MyStillAmbiguousType(3.2)
MyStillAmbiguousType(3.2)

julia> typeof(m)
MyType{Float64}

julia> typeof(t)
MyStillAmbiguousType

Тип поля a легко определяется из типа m, но не из типа t. Действительно, в t можно изменить тип поля a:

julia> typeof(t.a)
Float64

julia> t.a = 4.5f0
4.5f0

julia> typeof(t.a)
Float32

В отличие от этого, после того, как m создан, тип m.a не может измениться:

julia> m.a = 4.5f0
4.5f0

julia> typeof(m.a)
Float64

Тот факт, что тип m.a известен из типа m (в сочетании с тем, что его тип не может измениться во время работы функции), позволяет компилятору генерировать высокооптимизированный код для объектов, таких как m , но не для объектов, таких как t.

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

julia> m = MyType{AbstractFloat}(3.2)
MyType{AbstractFloat}(3.2)

julia> typeof(m.a)
Float64

julia> m.a = 4.5f0
4.5f0

julia> typeof(m.a)
Float32

Практически такие объекты ведут себя идентично объектам MyStillAmbiguousType.

Вполне поучительно сравнить объём кода, генерируемого для простой функции

func(m::MyType) = m.a+1

с помощью

code_llvm(func,(MyType{Float64},))
code_llvm(func,(MyType{AbstractFloat},))
code_llvm(func,(MyType,))

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

Избегайте полей с абстрактными контейнерами

Те же лучшие практики также работают для типов контейнеров:

julia> type MySimpleContainer{A<:AbstractVector}
         a::A
       end

julia> type MyAmbiguousContainer{T}
         a::AbstractVector{T}
       end

Например:

julia> c = MySimpleContainer(1:3);

julia> typeof(c)
MySimpleContainer{UnitRange{Int64}}

julia> c = MySimpleContainer([1:3;]);

julia> typeof(c)
MySimpleContainer{Array{Int64,1}}

julia> b = MyAmbiguousContainer(1:3);

julia> typeof(b)
MyAmbiguousContainer{Int64}

julia> b = MyAmbiguousContainer([1:3;]);

julia> typeof(b)
MyAmbiguousContainer{Int64}

Для MySimpleContainer, объект полностью определён своим типом и параметрами, поэтому компилятор может генерировать оптимизированные функции. В большинстве случаев этого будет достаточно.

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

function sumfoo(c::MySimpleContainer)
    s = 0
for x in c.a
    s += foo(x)
end
s
end

foo(x::Integer) = x
foo(x::AbstractFloat) = round(x)

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

Однако есть случаи, когда вам может потребоваться объявить разные версии внешней функции для разных типов элементов a. Вы могли бы сделать это так:

function myfun{T<:AbstractFloat}(c::MySimpleContainer{Vector{T}})
    ...
end
function myfun{T<:Integer}(c::MySimpleContainer{Vector{T}})
    ...
end

Это работает хорошо для Vector{T}, но нам также пришлось бы написать явные версии для UnitRange{T} или других абстрактных типов. Чтобы избежать такой рутины, вы можете использовать два параметра в объявлении MyContainer:

type MyContainer{T, A<:AbstractVector}
    a::A
end
MyContainer(v::AbstractVector) = MyContainer{eltype(v), typeof(v)}(v)

julia> b = MyContainer(1.3:5);

julia> typeof(b)
MyContainer{Float64,UnitRange{Float64}}

Обратите внимание на несколько удивительный факт, что T не отображается в объявлении поля a, а этот момент мы рассмотрим чуть позже. С помощью этого подхода можно написать функции, такие как:

function myfunc{T<:Integer, A<:AbstractArray}(c::MyContainer{T,A})
    return c.a[1]+1
end
# Note: because we can only define MyContainer for
# A<:AbstractArray, and any unspecified parameters are arbitrary,
# the previous could have been written more succinctly as
#     function myfunc{T<:Integer}(c::MyContainer{T})

function myfunc{T<:AbstractFloat}(c::MyContainer{T})
    return c.a[1]+2
end

function myfunc{T<:Integer}(c::MyContainer{T,Vector{T}})
    return c.a[1]+3
end

julia> myfunc(MyContainer(1:3))
2

julia> myfunc(MyContainer(1.0:3))
3.0

julia> myfunc(MyContainer([1:3]))
4

Как вы можете видеть, с помощью этого подхода можно специализироваться как на типе элемента T, так и на типе массива A.

Однако, остается одна проблема: мы не гарантируем, что A имеет тип элемента T, поэтому вполне возможно создание объекта такого вида:

julia> b = MyContainer{Int64, UnitRange{Float64}}(1.3:5);

julia> typeof(b)
MyContainer{Int64,UnitRange{Float64}}

Чтобы этого избежать, можно добавить внутренний конструктор:

type MyBetterContainer{T<:Real, A<:AbstractVector}
    a::A

    MyBetterContainer(v::AbstractVector{T}) = new(v)
end
MyBetterContainer(v::AbstractVector) = MyBetterContainer{eltype(v),typeof(v)}(v)


julia> b = MyBetterContainer(1.3:5);

julia> typeof(b)
MyBetterContainer{Float64,UnitRange{Float64}}

julia> b = MyBetterContainer{Int64, UnitRange{Float64}}(1.3:5);
ERROR: no method MyBetterContainer(UnitRange{Float64},)

Внутренний конструктор требует, чтобы тип элемента A был T.

Аннотирование значений, взятых из нетипизированных расположений

Часто удобно работать со структурами данных, которые могут содержать значения любого типа (массивы типа Array{Any}). Но если вы используете одну из этих структур и знаете тип элемента, то полезно поделиться этой информацией с компилятором:

function foo(a::Array{Any,1})
    x = a[1]::Int32
    b = x+1
    ...
end

В данном случае мы знаем, что первый элемент a будет Int32. Добавление такой аннотации имеет дополнительное преимущество: она вызовет ошибку во время выполнения, если значение не имеет ожидаемого типа, потенциально обнаруживая некоторые ошибки на более ранних этапах.

Объявление типов ключевых аргументов

Ключевые аргументы могут иметь объявленные типы:

function with_keyword(x; name::Int = 1)
    ...
end

Функции специализируются на типах ключевых аргументов, поэтому эти объявления не повлияют на производительность кода внутри функции. Однако они уменьшат накладные расходы при вызовах функции, которые включают ключевые аргументы.

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

Передача динамических списков ключевых аргументов, как в f(x; keywords...), может быть медленной и следует избегать в коде, чувствительном к производительности.

Разделение функций на несколько определений

Разбиение функции на несколько небольших определений позволяет компилятору напрямую вызывать наиболее подходящий код или даже встраивать его.

Вот пример «составной функции», которую следует писать как несколько определений:

function norm(A)
    if isa(A, Vector)
        return sqrt(real(dot(A,A)))
    elseif isa(A, Matrix)
        return max(svd(A)[2])
    else
        error("norm: invalid argument")
    end
end

Это можно написать более лаконично и эффективно как:

norm(x::Vector) = sqrt(real(dot(x,x)))
norm(A::Matrix) = max(svd(A)[2])

Написание «стабильных по типу» функций

В случае возможности, полезно гарантировать, что функция всегда возвращает значение одного и того же типа. Рассмотрим следующее определение:

pos(x) = x < 0 ? 0 : x

Хотя это кажется безобидным, проблема в том, что 0 — это целое число (типа Int), а x может быть любого типа. Таким образом, в зависимости от значения x, эта функция может возвращать значение одного из двух типов. Это поведение разрешено и может быть желательным в некоторых случаях. Но это легко исправить следующим образом:

pos(x) = x < 0 ? zero(x) : x

Также существует функция one() и более общая функция oftype(x,y), которая возвращает y, преобразованное к типу x.

Избегайте изменения типа переменной

Аналогичная проблема «стабильности типа» существует для переменных, многократно используемых в функции:

function foo()
    x = 1
    for i = 1:10
        x = x/bar()
    end
    return x
end

Локальная переменная x начинается как целое число, а после одной итерации цикла становится числом с плавающей точкой (результат оператора /). Это затрудняет компилятору оптимизацию тела цикла. Существует несколько возможных исправлений:

• Инициализировать x значением x = 1.0
• Объявить тип x: x::Float64 = 1
• Использовать явное преобразование: x = one(T)

Отделение функций ядра (также известные как барьеры функций)

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

function strange_twos(n)
    a = Array(rand(Bool) ? Int64 : Float64, n)
    for i = 1:n
        a[i] = 2
    end
    return a
end

Это следует написать как:

function fill_twos!(a)
    for i=1:length(a)
        a[i] = 2
    end
end

function strange_twos(n)
    a = Array(rand(Bool) ? Int64 : Float64, n)
    fill_twos!(a)
    return a
end

Компилятор Julia специализирует код для типов аргументов на границах функций, поэтому в исходной реализации он не знает тип a во время цикла (так как он выбирается случайно). Поэтому второй вариант обычно быстрее, так как внутренний цикл может быть перекомпилирован как часть fill_twos! для разных типов a.

Второй вариант также часто имеет лучший стиль и может привести к большей повторной переиспользованию кода.

Этот шаблон используется в нескольких местах стандартной библиотеки. Например, см. hvcat_fill в abstractarray.jl или функцию fill!, которую можно было бы использовать вместо написания собственной fill_twos!.

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

Типы со значениями-параметрами

Предположим, что вы хотите создать N-мерный массив, размер которого равен 3 по каждой оси. Такие массивы можно создать следующим образом:

A = fill(5.0, (3, 3))

Этот подход работает очень хорошо: компилятор может определить, что A является Array{Float64,2}, потому что он знает тип значения заполнения (5.0::Float64) и размерность ((3, 3)::NTuple{2,Int}). Это означает, что компилятор может генерировать очень эффективный код для любого последующего использования A в той же функции.

Но теперь предположим, что вы хотите написать функцию, которая создаёт 3×3×... массив в произвольных измерениях; вы, возможно, захотите написать функцию

function array3(fillval, N)
    fill(fillval, ntuple(d->3, N))
end

Это работает, но (как вы можете проверить самостоятельно с помощью @code_warntype array3(5.0, 2)) проблема в том, что тип результата не может быть выведен: аргумент N — это значение типа Int, и вычисление типов не может (и не может) предсказать его значение заранее. Это означает, что код, использующий результат этой функции, должен быть консервативным, проверяя тип при каждом обращении к A; такой код будет очень медленным.

Теперь, очень хороший способ решения таких проблем — использовать технику «барьера функций». Однако в некоторых случаях вы можете захотеть полностью устранить нестабильность типов. В таких случаях один из подходов заключается в передаче размерности в качестве параметра, например, через Val{T} (см. «Типы со значениями»):

function array3{N}(fillval, ::Type{Val{N}})
    fill(fillval, ntuple(d->3, Val{N}))
end

Julia имеет специализированную версию ntuple , которая принимает Val{::Int} в качестве второго параметра; передавая N в качестве параметра типа, вы делаете его «значение» известным компилятору. В результате эта версия array3 позволяет компилятору предсказать тип возвращаемого значения.

Однако использование таких техник может быть удивительно тонким. Например, это не поможет, если вы вызовете array3 из функции такого вида:

function call_array3(fillval, n)
    A = array3(fillval, Val{n})
end

Здесь вы создали ту же проблему снова: компилятор не может угадать тип n, поэтому он не знает тип Val{n} . Попытка использовать Val, но сделать это неправильно, может легко привести к ухудшению производительности во многих ситуациях. (Только в ситуациях, когда вы фактически объединяете Val с техникой барьера функций, чтобы сделать функцию ядра более эффективной, код, подобный вышеуказанному, следует использовать.)

Пример правильного использования Val был бы:

function filter3{T,N}(A::AbstractArray{T,N})
    kernel = array3(1, Val{N})
    filter(A, kernel)
end

В этом примере N передаётся в качестве параметра, поэтому его «значение» известно компилятору. По существу, Val{T} работает только тогда, когда T либо жёстко закодирован (Val{3} ), либо уже задан в области типов.

Опасности злоупотребления множественным диспетчерированием (а также больше о типах со значениями-параметрами)

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

immutable Car{Make,Model}
    year::Int
    ...more fields...
end

а затем вызываете диспетчеризацию на объектах, таких как Car{:Honda,:Accord}(year, args...).

Это может быть полезным, если выполняются следующие условия:

• Вам необходимы CPU-ёмкие вычисления на каждом Car, и это становится гораздо более эффективным, если вы знаете Make и Model во время компиляции.
• У вас есть однородные списки объектов одного и того же типа Car для обработки, так что вы можете сохранить их все в Array{Car{:Honda,:Accord},N}.

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

Когда это не так, то, скорее всего, вы не получите никакой выгоды; что хуже, возникающий «комбинаторный взрыв типов» будет контрпродуктивным. Если items[i+1] имеет другой тип, чем item[i], Julia должна найти тип во время выполнения, найти соответствующий метод в таблицах методов, решить (через пересечение типов), какой из них подходит, определить, был ли он уже JIT-скомпилирован (и сделать это, если нет), и затем сделать вызов. По сути, вы просите всю систему типов и механизм JIT-компиляции выполнить примерно эквивалент инструкции switch или поиска в словаре в вашем собственном коде.

Некоторые тесты производительности во время выполнения, сравнивающие (1) диспетчерирование по типам, (2) поиск в словаре и (3) инструкцию «switch», можно найти на почтовой рассылке.

Возможно, еще хуже, чем влияние на время выполнения, — это влияние на время компиляции: Julia будет компилировать специализированные функции для каждого разного Car{Make, Model}; если у вас сотни или тысячи таких типов, то каждая функция, которая принимает такой объект в качестве параметра (от пользовательской get_year функции, которую вы можете написать сами, до универсальной push! функции в стандартной библиотеке), будет иметь сотни или тысячи вариантов, скомпилированных для нее. Каждый из них увеличивает размер кэша скомпилированного кода, длину внутренних списков методов и т. д. Чрезмерный энтузиазм по поводу значений в качестве параметров может легко привести к огромным потерям ресурсов.

Доступ к массивам в порядке памяти по столбцам

Многомерные массивы в Julia хранятся в порядке следования столбцов. Это означает, что массивы выстраиваются по одному столбцу за раз. Это можно проверить, используя функцию vec или синтаксис [:], как показано ниже (обратите внимание, что массив упорядочен [1 3 2 4], а не [1 2 3 4]):

julia> x = [1 2; 3 4]
2×2 Array{Int64,2}:
 1  2
 3  4

julia> x[:]
4-element Array{Int64,1}:
 1
 3
 2
 4

Эта конвенция для упорядочения массивов широко используется во многих языках, таких как Fortran, Matlab и R (и не только). Альтернативой порядку следования столбцов является порядок следования строк, который является конвенцией, принятой C и Python (numpy) среди других языков. Учёт порядка массивов может существенно повлиять на производительность при циклическом проходе по массивам. Правило большого пальца, которое следует помнить, заключается в том, что с массивами в порядке следования столбцов первый индекс изменяется быстрее всего. По существу, это означает, что циклический проход будет быстрее, если индекс внутреннего цикла — первый, который появляется в выражении среза.

Рассмотрим следующий искусственный пример. Предположим, что мы хотим написать функцию, которая принимает Vector и возвращает квадратный Matrix с заполненными строками или столбцами копиями входного вектора. Предположим, что порядок заполнения — строки или столбцы — не имеет значения (возможно, остальная часть кода может быть легко адаптирована соответственно). Мы могли бы сделать это по меньшей мере четырьмя способами (кроме рекомендуемого вызова встроенной функции repmat()):

function copy_cols{T}(x::Vector{T})
    n = size(x, 1)
    out = Array{T}(n, n)
    for i=1:n
        out[:, i] = x
    end
    out
end

function copy_rows{T}(x::Vector{T})
    n = size(x, 1)
    out = Array{T}(n, n)
    for i=1:n
        out[i, :] = x
    end
    out
end

function copy_col_row{T}(x::Vector{T})
    n = size(x, 1)
    out = Array{T}(n, n)
    for col=1:n, row=1:n
        out[row, col] = x[row]
    end
    out
end

function copy_row_col{T}(x::Vector{T})
    n = size(x, 1)
    out = Array{T}(n, n)
    for row=1:n, col=1:n
        out[row, col] = x[col]
    end
    out
end

Теперь мы измерим время работы каждой из этих функций, используя тот же случайный 10000 на 1 входной вектор:

julia> x = randn(10000);

julia> fmt(f) = println(rpad(string(f)*": ", 14, ' '), @elapsed f(x))

julia> map(fmt, Any[copy_cols, copy_rows, copy_col_row, copy_row_col]);
copy_cols:    0.331706323
copy_rows:    1.799009911
copy_col_row: 0.415630047
copy_row_col: 1.721531501

Обратите внимание, что copy_cols намного быстрее, чем copy_rows. Это ожидаемо, так как copy_cols учитывает структуру хранения массива в порядке следования столбцов и заполняет его по одному столбцу за раз. Кроме того, copy_col_row намного быстрее, чем copy_row_col, поскольку он следует нашему правилу большого пальца, что первый элемент, появляющийся в выражении среза, должен быть связан с внутренним циклом.

Предварительное выделение памяти

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

Иногда можно избежать необходимости выделения памяти на каждом вызове функции, предварительно выделив выходные данные. В качестве тривиального примера сравните

function xinc(x)
    return [x, x+1, x+2]
end

function loopinc()
    y = 0
    for i = 1:10^7
        ret = xinc(i)
        y += ret[2]
    end
    y
end

с

function xinc!{T}(ret::AbstractVector{T}, x::T)
    ret[1] = x
    ret[2] = x+1
    ret[3] = x+2
    nothing
end

function loopinc_prealloc()
    ret = Array{Int}(3)
    y = 0
    for i = 1:10^7
        xinc!(ret, i)
        y += ret[2]
    end
    y
end

Результаты измерения времени:

julia> @time loopinc()
elapsed time: 1.955026528 seconds (1279975584 bytes allocated)
50000015000000

julia> @time loopinc_prealloc()
elapsed time: 0.078639163 seconds (144 bytes allocated)
50000015000000

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

Доведенное до крайности, предварительное выделение может сделать ваш код менее удобочитаемым, поэтому могут потребоваться измерения производительности и некоторое суждение. Однако для «векторизованных» (элементных) функций удобный синтаксис x .= f.(y) может использоваться для операций на месте с объединенными циклами и без временных массивов (Точечная нотация для векторизации функций).

Избегайте интерполяции строк для ввода-вывода

При записи данных в файл (или другое устройство ввода-вывода) формирование дополнительных промежуточных строк является источником накладных расходов. Вместо:

println(file, "$a $b")

используйте:

println(file, a, " ", b)

В первом варианте кода формируется строка, затем она записывается в файл, а во втором значения записываются непосредственно в файл. Обратите также внимание, что в некоторых случаях интерполяция строк может быть сложнее для чтения. Рассмотрите:

println(file, "$(f(a))$(f(b))")

по сравнению с:

println(file, f(a), f(b))

Оптимизация сетевого ввода-вывода во время параллельного выполнения

При выполнении удаленной функции в параллельном режиме:

responses = Vector{Any}(nworkers())
@sync begin
    for (idx, pid) in enumerate(workers())
        @async responses[idx] = remotecall_fetch(pid, foo, args...)
    end
end

быстрее, чем:

refs = Vector{Any}(nworkers())
for (idx, pid) in enumerate(workers())
    refs[idx] = @spawnat pid foo(args...)
end
responses = [fetch(r) for r in refs]

Первый вариант приводит к одному сетевому обмену с каждым рабочим процессом, а второй — к двум сетевым вызовам: сначала от @spawnat , а второй — из-за fetch (или даже из-за wait). fetch/wait также выполняется последовательно, что приводит к снижению общей производительности.

Исправление предупреждений о устаревании

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

Корректировки

Вот некоторые незначительные моменты, которые могут помочь в узких циклах.

• Избегайте ненужных массивов. Например, вместо sum([x,y,z]) используйте x+y+z.
• Используйте abs2(z) вместо abs(z)^2 для комплексных z. В общем, старайтесь переписать код, чтобы использовать abs2() вместо abs() для комплексных аргументов.
• Используйте div(x,y) для усечения целочисленного деления вместо trunc(x/y), fld(x,y) вместо floor(x/y) и cld(x,y) вместо ceil(x/y).

Аннотации производительности

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

• Используйте @inbounds для устранения проверки границ массива в выражениях. Убедитесь в этом перед тем, как сделать это. Если индексы выходят за пределы границ, могут возникнуть сбои или неявное повреждение данных.
• Используйте @fastmath для разрешения оптимизации операций с плавающей точкой, которые верны для вещественных чисел, но приводят к различиям для чисел IEEE. Будьте осторожны при использовании этого, так как это может изменить численные результаты. Это соответствует опции -ffast-math clang.
• Напишите @simd перед for циклами, которые допускают векторизацию. Эта функция находится в стадии тестирования и может измениться или исчезнуть в будущих версиях Julia.

Примечание: хотя @simd необходимо разместить непосредственно перед циклом, @inbounds и @fastmath можно применять к нескольким операторам одновременно, например, используя begin ... end, или даже к всей функции.

Вот пример с обоих видов разметки @inbounds и @simd.

function inner( x, y )
    s = zero(eltype(x))
    for i=1:length(x)
        @inbounds s += x[i]*y[i]
    end
    s
end

function innersimd( x, y )
    s = zero(eltype(x))
    @simd for i=1:length(x)
        @inbounds s += x[i]*y[i]
    end
    s
end

function timeit( n, reps )
    x = rand(Float32,n)
    y = rand(Float32,n)
    s = zero(Float64)
    time = @elapsed for j in 1:reps
        s+=inner(x,y)
    end
    println("GFlop/sec        = ",2.0*n*reps/time*1E-9)
    time = @elapsed for j in 1:reps
        s+=innersimd(x,y)
    end
    println("GFlop/sec (SIMD) = ",2.0*n*reps/time*1E-9)
end

timeit(1000,1000)

На компьютере с процессором Intel Core i5 с частотой 2,4 ГГц это даёт:

GFlop/sec        = 1.9467069505224963
GFlop/sec (SIMD) = 17.578554163920018

(GFlop/sec измеряет производительность, и более высокие значения лучше.) Диапазон для цикла @simd for должен быть одномерным диапазоном. Переменная, используемая для накопления, такая как s в примере, называется переменной накопления. Используя @simd, вы утверждаете несколько свойств цикла:

• Итерации можно безопасно выполнять в произвольном или перекрывающемся порядке, с учётом специальных переменных накопления.
• Операции с плавающей точкой на переменных накопления могут быть переупорядочены, что, возможно, приведёт к другим результатам, чем без @simd.
• Ни одна итерация не ожидает завершения другой итерации, чтобы сделать дальнейший прогресс.

Цикл, содержащий break, continue, или @goto, вызовет ошибку во время компиляции.

Использование @simd просто даёт компилятору разрешение на векторизацию. Реально ли это будет сделано, зависит от компилятора. Чтобы получить реальную выгоду от текущей реализации, ваш цикл должен обладать следующими дополнительными свойствами:

• Цикл должен быть внутренним циклом.
• Тело цикла должно быть последовательным кодом. Вот почему @inbounds в настоящее время требуется для всех обращений к массивам. Компилятор иногда может преобразовать короткие &&, ||, и ?: выражения в последовательный код, если безопасно вычислить все операнды безусловно. Рассмотрите использование ifelse() вместо ?: в цикле, если это безопасно.
• Доступы должны иметь структуру шага и не могут быть «сборками» (считываниями с произвольными индексами) или «разбросами» (записями с произвольными индексами).
• Шаг должен быть единичным.
• В некоторых простых случаях, например, с 2-3 массивами, обращёнными в цикле, LLVM может автоматически выполнить автовекторизацию, что не приведёт к дальнейшему ускорению с @simd.

Вот пример со всеми тремя видами разметки. Эта программа сначала вычисляет конечную разность одномерного массива, а затем вычисляет L2-норму результата:

function init!(u)
    n = length(u)
    dx = 1.0 / (n-1)
    @fastmath @inbounds @simd for i in 1:n
        u[i] = sin(2pi*dx*i)
    end
end

function deriv!(u, du)
    n = length(u)
    dx = 1.0 / (n-1)
    @fastmath @inbounds du[1] = (u[2] - u[1]) / dx
    @fastmath @inbounds @simd for i in 2:n-1
        du[i] = (u[i+1] - u[i-1]) / (2*dx)
    end
    @fastmath @inbounds du[n] = (u[n] - u[n-1]) / dx
end

function norm(u)
    n = length(u)
    T = eltype(u)
    s = zero(T)
    @fastmath @inbounds @simd for i in 1:n
        s += u[i]^2
    end
    @fastmath @inbounds return sqrt(s/n)
end

function main()
    n = 2000
    u = Array{Float64}(n)
    init!(u)
    du = similar(u)

    deriv!(u, du)
    nu = norm(du)

    @time for i in 1:10^6
        deriv!(u, du)
        nu = norm(du)
    end

    println(nu)
end

main()

На компьютере с процессором Intel Core i7 с частотой 2,7 ГГц это даёт:

$ julia wave.jl;
elapsed time: 1.207814709 seconds (0 bytes allocated)

$ julia --math-mode=ieee wave.jl;
elapsed time: 4.487083643 seconds (0 bytes allocated)

Здесь опция --math-mode=ieee отключает макрос @fastmath, чтобы мы могли сравнить результаты.

В этом случае ускорение, обусловленное @fastmath, составляет примерно 3,7. Это необычно велико – обычно ускорение будет меньше. (В этом конкретном примере рабочий набор бенчмарка достаточно мал, чтобы поместиться в кэш L1 процессора, так что задержка доступа к памяти не играет роли, а время вычислений определяется использованием ЦП. Во многих реальных программах это не так.) Также в этом случае данная оптимизация не меняет результат – обычно результат будет немного отличаться. В некоторых случаях, особенно для числово неустойчивых алгоритмов, результат может быть сильно другим.

Аннотация @fastmath переупорядочивает выражения с плавающей точкой, например, изменяя порядок вычислений или предполагая, что некоторые особые случаи (inf, nan) не могут возникнуть. В этом случае (и на этом конкретном компьютере) основное различие заключается в том, что выражение 1 / (2*dx) в функции deriv вынесено за цикл (т. е. вычислено вне цикла), как если бы было написано idx = 1 / (2*dx). В цикле выражение ... / (2*dx) затем становится ... * idx, что гораздо быстрее для вычисления. Разумеется, как сама применяемая компилятором оптимизация, так и полученное ускорение в значительной степени зависят от аппаратного обеспечения. Вы можете просмотреть изменения в сгенерированном коде, используя функцию Julia code_native().

Обработка субнормальных чисел как нулей

Субнормальные числа, ранее называвшиеся денормализованными числами, полезны во многих контекстах, но на некоторых аппаратных платформах влекут за собой потери производительности. Вызов set_zero_subnormals(true) дает разрешение операциям с плавающей точкой обрабатывать субнормальные входные или выходные данные как нули, что может повысить производительность на некоторых аппаратных платформах. Вызов set_zero_subnormals(false) обеспечивает строгое соответствие стандарту IEEE для субнормальных чисел.

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

function timestep{T}( b::Vector{T}, a::Vector{T}, Δt::T )
    @assert length(a)==length(b)
    n = length(b)
    b[1] = 1                            # Boundary condition
    for i=2:n-1
        b[i] = a[i] + (a[i-1] - T(2)*a[i] + a[i+1]) * Δt
    end
    b[n] = 0                            # Boundary condition
end

function heatflow{T}( a::Vector{T}, nstep::Integer )
    b = similar(a)
    for t=1:div(nstep,2)                # Assume nstep is even
        timestep(b,a,T(0.1))
        timestep(a,b,T(0.1))
    end
end

heatflow(zeros(Float32,10),2)           # Force compilation
for trial=1:6
    a = zeros(Float32,1000)
    set_zero_subnormals(iseven(trial))  # Odd trials use strict IEEE arithmetic
    @time heatflow(a,1000)
end

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

Обработку субнормалей как нулей следует использовать с осторожностью, поскольку при этом нарушаются некоторые тождества, например, x-y==0 подразумевает x==y:

julia> x=3f-38; y=2f-38;

julia> set_zero_subnormals(false); (x-y,x==y)
(1.0000001f-38,false)

julia> set_zero_subnormals(true); (x-y,x==y)
(0.0f0,false)

В некоторых приложениях альтернативой нулевым субнормальным числам является введение небольшого шума. Например, вместо инициализации a нулями, инициализируйте её:

a = rand(Float32,1000) * 1.f-9

@code_warntype

Макрос @code_warntype (или его функциональный вариант code_warntype()) иногда может быть полезен для диагностики проблем, связанных с типами. Вот пример:

pos(x) = x < 0 ? 0 : x

function f(x)
    y = pos(x)
    sin(y*x+1)
end

julia> @code_warntype f(3.2)
Variables:
  x::Float64
  y::UNION(INT64,FLOAT64)
  _var0::Float64
  _var3::Tuple{Int64}
  _var4::UNION(INT64,FLOAT64)
  _var1::Float64
  _var2::Float64

Body:
  begin  # none, line 2:
      _var0 = (top(box))(Float64,(top(sitofp))(Float64,0))
      unless (top(box))(Bool,(top(or_int))((top(lt_float))(x::Float64,_var0::Float64)::Bool,(top(box))(Bool,(top(and_int))((top(box))(Bool,(top(and_int))((top(eq_float))(x::Float64,_var0::Float64)::Bool,(top(lt_float))(_var0::Float64,9.223372036854776e18)::Bool)),(top(slt_int))((top(box))(Int64,(top(fptosi))(Int64,_var0::Float64)),0)::Bool)))) goto 1
      _var4 = 0
      goto 2
      1:
      _var4 = x::Float64
      2:
      y = _var4::UNION(INT64,FLOAT64) # line 3:
      _var1 = y::UNION(INT64,FLOAT64) * x::Float64::Float64
      _var2 = (top(box))(Float64,(top(add_float))(_var1::Float64,(top(box))(Float64,(top(sitofp))(Float64,1))))
      return (GlobalRef(Base.Math,:nan_dom_err))((top(ccall))($(Expr(:call1, :(top(tuple)), "sin", GlobalRef(Base.Math,:libm))),Float64,$(Expr(:call1, :(top(tuple)), :Float64)),_var2::Float64,0)::Float64,_var2::Float64)::Float64
  end::Float64

Интерпретация вывода @code_warntype, как и его аналогов @code_lowered, @code_typed, @code_llvm и @code_native, требует некоторой практики. Ваш код представлен в форме, которая была частично обработана на пути к генерации скомпилированного машинного кода. Большинство выражений снабжены аннотацией типа, обозначенной ::T (где T может быть, например, Float64). Наиболее важной характеристикой @code_warntype является то, что неконкретные типы отображаются красным цветом; в приведенном выше примере такой вывод показан заглавными буквами.

Верхняя часть вывода обобщает информацию о типе для различных внутренних переменных функции. Вы можете увидеть, что y, одна из созданных вами переменных, является Union{Int64,Float64}, из-за неустойчивости типа pos. Существует еще одна переменная _var4, которая также имеет тот же тип.

Следующие строки представляют тело f . Строки, начинающиеся с числа, за которым следует двоеточие (1:, 2: ), являются метками и представляют собой целевые адреса для переходов (через goto) в вашем коде. Рассматривая тело, можно увидеть, что pos было вложено в f—всё до 2: происходит из кода, определенного в pos.

Начиная с 2:, переменная y определена и снова аннотирована как тип Union . Далее мы видим, что компилятор создал временную переменную _var1 для хранения результата y*x . Поскольку Float64 умножается либо на Int64 , либо на Float64, результат представляет собой Float64, вся неустойчивость типа заканчивается здесь. В итоге f(x::Float64) не будет иметь неустойчивость типа в своём выводе, даже если некоторые промежуточные вычисления являются неустойчивыми по типу.

Как вы будете использовать эту информацию, зависит от вас. Очевидно, лучше всего было бы исправить pos , чтобы он был устойчив к типу: если вы это сделаете, все переменные в f будут конкретными, и его производительность будет оптимальной. Однако существуют обстоятельства, когда такая временная неустойчивость типа может не иметь большого значения: например, если pos никогда не используется изолированно, тот факт, что вывод f устойчив к типу (для входных значений Float64), защитит последующий код от распространения эффектов неустойчивости типа. Это особенно актуально в тех случаях, когда исправление неустойчивости типа затруднено или невозможно: например, в настоящее время невозможно определить тип возвращаемого значения анонимной функции. В таких случаях указанные выше советы (например, добавление аннотаций типов и/или разделение функций) являются лучшими инструментами для предотвращения «повреждений» от неустойчивости типа.

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

• Тело функции заканчивается на end::Union{T1,T2})
  • Интерпретация: функция с неустойчивым типом возврата
  • Предложение: сделайте тип возвращаемого значения устойчивым к типу, даже если придётся его аннотировать
• f(x::T)::Union{T1,T2}
  • Интерпретация: вызов функции с неустойчивым типом
  • Предложение: исправить функцию или, при необходимости, аннотировать тип возврата
• (top(arrayref))(A::Array{Any,1},1)::Any
  • Интерпретация: доступ к элементам массивов с неустойчивым типом
  • Предложение: использовать массивы с более определёнными типами или, при необходимости, аннотировать тип отдельных элементов массива
• (top(getfield))(A::ArrayContainer{Float64},:data)::Array{Float64,N}
  • Интерпретация: получение поля, являющегося типом с неустойчивым типом. В этом случае ArrayContainer имело поле data::Array{T} . Но Array также нуждается в измерении N для того, чтобы иметь конкретный тип.
  • Предложение: использовать конкретные типы, такие как Array{T,3} или Array{T,N}, где N теперь является параметром ArrayContainer