特殊 (besselj) 函数矩阵

Question

我是 julia 的新手，所以我欢迎提出一些建议来改进以下功能，

using SpecialFunctions

function rb(x, nu_max)

  bj = Array{Complex64}(length(x), nu_max)

  nu = 0.5 + (0:nu_max)
  # somehow dot broadcast isn't happy
  # bj .= [ besselj(_nu,_x)*sqrt(pi/2*_x) for _nu in nu, _x in x]

  bj   = [ besselj(_nu,_x)*sqrt(pi/2*_x) for _nu in nu, _x in x]

end

rb(1.0:0.1:2.0, 500)

基本上，我不太确定在这两个参数（x 和 nu）上获得矩阵的推荐方法是什么。 documentation doesn't offer much information, but I understand that the underlying fortran routine internally loops over nu，所以为了性能，我宁愿不再这样做。

编辑： 我被问及目标；它是为 $x$ 和 $\nu$ 的多个值计算 Riccati-Bessel 函数 $j_1(x,\nu),h_1(x,\nu)$。

我已经从原始版本中剥离了文体问题，以专注于这个核心问题。

Answer 1

详细说明我上面的评论。乍一看，一般来说，通过预先分配数组并就地填充它们（例如使用 dot broadcasting）来尝试避免临时分配。也可以使用 @inbounds.

给大家留下印象，之后

using SpecialFunctions
x = 1.0
nu_max = 3
nu = 0.5 + (0:nu_max)
f(nu,x) = besselj.(nu,x).*sqrt.(pi/2*x)

比较（使用BenchmarkTools）

的性能（和分配）

bj   = hcat([ besselj.(_nu,x).*sqrt.(pi/2*x) for _nu in nu]...)

和

f.(nu,x)

（技术上输出不相同，您必须使用上面的 vcat，但无论如何）

更新（在 OP 净化他的代码之后）：

好的，我想我（终于）现在明白你真正的问题了（抱歉）。我上面所说的是关于优化您的原始代码如何调用 besselj 并有效地处理它的输出（请参阅@Matt B. 的 post 以获得完整的广播解决方案）。

IIUC，你想在给定 nu 和 [=18 的 besselj 的计算中利用这个事实（我不知道也没有检查这是否真的是真的） =] 在内部对 nu 进行求和。换句话说，您想使用此内部求和的中间结果来避免冗余计算。

由于 SpecialFunctions 的 besselj 似乎只是调用 Fortran 例程（可能 here），我怀疑您是否可以访问任何此类信息。不幸的是我不能在这里帮助你（我可能会寻找 besselj 的纯 Julia 实现）。

Answer 2

这是一个很好的例子，您可以在其中充分利用广播。看起来您想要 x 和 nu 之间的笛卡尔积，其中行由 nu 的值填充，列为 x。这正是广播可以做的——你只需要重塑 x 使其成为跨多列的单行：

julia> using SpecialFunctions

julia> x = 1.0:0.1:2.0
1.0:0.1:2.0

julia> nu = 0.5 + (0:500)
0.5:1.0:500.5

 # this shows how broadcast works — these are the arguments and their location in the matrix
julia> tuple.(nu, reshape(x, 1, :))
501×11 Array{Tuple{Float64,Float64},2}:
 (0.5, 1.0)    (0.5, 1.1)    …  (0.5, 1.9)    (0.5, 2.0)
 (1.5, 1.0)    (1.5, 1.1)       (1.5, 1.9)    (1.5, 2.0)
 (2.5, 1.0)    (2.5, 1.1)       (2.5, 1.9)    (2.5, 2.0)
 (3.5, 1.0)    (3.5, 1.1)       (3.5, 1.9)    (3.5, 2.0)
 ⋮                           ⋱                ⋮
 (497.5, 1.0)  (497.5, 1.1)     (497.5, 1.9)  (497.5, 2.0)
 (498.5, 1.0)  (498.5, 1.1)     (498.5, 1.9)  (498.5, 2.0)
 (499.5, 1.0)  (499.5, 1.1)     (499.5, 1.9)  (499.5, 2.0)
 (500.5, 1.0)  (500.5, 1.1)  …  (500.5, 1.9)  (500.5, 2.0)

julia> bj = besselj.(nu,reshape(x, 1, :)).*sqrt.(pi/2*reshape(x, 1, :))
501×11 Array{Float64,2}:
 0.841471    0.891207     0.932039     …  0.9463       0.909297
 0.301169    0.356592     0.414341        0.821342     0.870796
 0.0620351   0.0813173    0.103815        0.350556     0.396896
 0.00900658  0.0130319    0.0182194       0.101174     0.121444
 ⋮                                     ⋱               ⋮
 0.0         0.0          0.0             0.0          0.0
 0.0         0.0          0.0             0.0          0.0
 0.0         0.0          0.0             0.0          0.0
 0.0         0.0          0.0          …  0.0          0.0

特殊 (besselj) 函数矩阵

matrix of special (besselj) functions

bessel-functions

julia