在Julia NFFT中，用單個陣列方法替換矢量和矩陣方法

在NFFT包中，矢量和矩陣有快速專用功能，通用陣列有慢速功能。我想盡快地完成一般功能和專門的功能，但我遇到了麻煩。專門的功能在Julia NFFT中，用單個陣列方法替換矢量和矩陣方法

一個家庭基本上是這個樣子（offset計算在實際的代碼，但它的值並不重要，這個問題）：

myfunc!(f::Vector, g::Vector) 
    offset = 5 
    n = length(g) 
    N = length(f) 

    for l = 1:N 
     g[ ((l+offset)%n)+1 ] = f[l] 
    end 
end

和

myfunc!(f::Matrix, g::Matrix) 
    offsetx = 5 
    offsety = 5 
    n1, n2 = size(g) 
    N1, N2 = size(f) 

    for ly = 1:N2 
     for lx = 1:N1 
      g[ ((lx+offsetx)%n1)+1, ((ly+offsety)%n2)+1 ] = f[lx, ly] 
     end 
    end 
end

一般功能可以這樣寫：

myfunc!{D}(f::Array{D}, g::Array{D}) 
    offset = ntuple(d -> 5, D) 
    n = size(g) 

    for l in CartesianRange(size(f)) 
     idx = CartesianIndex{D}(ntuple(d -> ((l[d]+offset[d])%n[d])+1, D)) 
     g[ idx ] = f[l] 
    end 
end

不幸的是，這是非常緩慢的。大部分時間都花在循環中的ntuple上。

爲idx其他可能性包括讓idx = Array{Int}(D)，使內環看起來像

for d = 1:D 
    idx[d] = ((l[d]+offset[d])%n[d])+1 
end 
g[idx...] = f[l]

這也慢。

我在想，既然尺寸D是一種說法，一個@generated功能可用於計算idx進行，但我無法弄清楚如何做到這一點（或者如果有更好的方法）。

我正在使用Julia v0.4.5。

來源

2016-08-11 Robert

如何使用生成的函數執行此操作的答案是以使用Base.Cartesian輔助宏。

using Base.Cartesian 

@generated function myfunc!{T,D}(f::Array{T,D}, g::Array{T,D}) 
    quote 
     @nexprs $D d->offset_d=5 #Declase offset_1=5, offset_2=5 etc 

     @nloops $D l f begin 
      (@nref $D g d->(l_d+offset_d)%size(g,d)+1) = @nref $D f l 
     end 
    end 
end

我已經確認這是正確的，至少對2D來說。我把它作爲一個消息給讀者來分析它。它或多或少不分配。

來源

2016-08-11 07:43:42

非常好的解決方案！用這種方法我簡化了OP。我也有第三個輸入'h :: Vector {Vector {Float64}}'，並且（在2D情況下）內循環從 g [（（lx + offsetx）％n1）+1，（（ly + offsty）％n2）+1] = f [lx，ly]' 至 g [（（lx + offsetx）％n1）+1，（（ly + offsety）％n2）+1] = f [1x ，ly] * h [1] [lx] * h [2] [ly]'。這些術語也可以使用'Base.Cartesian'包含在內。它看起來像嵌套'@ nref'是一個不行。 – Robert

當然：將最內部分（'@ nref'）改爲： 'v = @nref $ D f l'; '@nexprs $ D d-> v * = h [d] [l_d]'; （@nref $ D gd - >（l_d + offset_d）％size（g，d）+1）= v'; –

在Julia NFFT中，用單個陣列方法替換矢量和矩陣方法

回答

相關問題