Julia LoopVectorization ArgumentError: недопустимый индекс: VectorizationBase.MM{2, 1, Int64}<2, 3> типа VectorizationBase.MM{2, 1, Int64}

avatar
crevell
1 июля 2021 в 17:35
78
1
1

Я пытаюсь решить 2D PDE на большой конечно-разностной сетке. У меня есть собственная небольшая функция для выполнения операции Лапласа над этой сеткой. В настоящее время модель довольно медленная, и я пытаюсь ускорить матричные операции. Я добавил макрос @turbo из LoopVectorization.jl, чтобы посмотреть, какая разница. Код Лапласа теперь выглядит следующим образом: 

module Laplacian

using LinearAlgebra
using LoopVectorization

@inline @views function ∇²!(∇²u, u, N, h, a)

    @turbo for j = 2+a:N+5-a
        for i = 2+a:N+5-a
            ∇²u[i,j] = (u[i-1,j] + u[i+1,j] + u[i,j-1] + u[i,j+1] - 4*u[i,j])/h^2
        end
    end

    return nothing

end

export ∇²!

end

∇²u — предварительно выделенный массив для хранения результата, например ∇²u = zeros(N+6,N+6)

u — это поле, к которому применяется лапласиан, например, u = rand(N+6,N+6)

N — количество точек сетки в u (минус 6 призрачных точек в каждом измерении, например, N=400

h — расстояние между точками сетки, например h=0.1

a — количество ложных точек, которые нужно пропустить на краях матрицы, например a=1

К сожалению, это приводит к следующей ошибке: ERROR: ArgumentError: invalid index: VectorizationBase.MM{2, 1, Int64}<4, 5> of type VectorizationBase.MM{2, 1, Int64}

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

Источник
crevell
1 июля 2021 в 17:45
0

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

mcabbott
1 июля 2021 в 18:18
2

Я думаю, проблема в макросе @views. Если вы запустите @macroexpand1 на этом, вы увидите, что то, что получает @turbo, выглядит совсем не так, как обычно ожидается.

Ответы (1)

avatar
cbk
1 июля 2021 в 18:33
2

Следуя комментарию @mcabbot, я могу подтвердить, что проблема связана с макросом @views. К счастью, этот макрос также совершенно не нужен в этом случае, поскольку вы не берете какие-либо срезы массива, которые в любом случае выделялись бы в любой точке этой функции. Так, например, следующее работает без ошибок в моей системе: 

using LinearAlgebra
using LoopVectorization

@inline function ∇²!(∇²u, u, N, h, a)

    @turbo for j = 2+a:N+5-a
        for i = 2+a:N+5-a
            ∇²u[i,j] = (u[i-1,j] + u[i+1,j] + u[i,j-1] + u[i,j+1] - 4*u[i,j])/h^2
        end
    end

    return ∇²u # There is no performance cost to returning the result, so might as well!

end

N=400
∇²u = zeros(N+6,N+6)
u = rand(N+6,N+6)
h=0.1
a=1

∇²!(∇²u, u, N, h, a)

Это также должно работать нормально, если вам нужно передать view в эту функцию в качестве аргумента — просто не помещайте макрос @views в определение функции:

julia> some_big_array = zeros(1000,1000);

julia> ∇²u = view(some_big_array, 1:N+6, 1:N+6); # or equivalently ∇²u = @views some_big_array[1:N+6, 1:N+6]

julia> ∇²!(∇²u, u, N, h, a)
406×406 view(::Matrix{Float64}, 1:406, 1:406) with eltype Float64:
  ...
crevell
1 июля 2021 в 19:40
1

Спасибо, это очень полезно. У меня есть дурная привычка использовать макрос @views для функций и предполагать, что это не повредит.