特殊(贝塞尔)函数矩阵

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)上获取矩阵的推荐方法是什么。该文档没有提供太多信息，但我知道底层 fortran 例程在 nu 上内部循环，因此为了性能，我宁愿不再这样做。

---

编辑：我被问及目标;它是计算 Riccati-Bessel 函数 $j_1(x，u)，h_1(x，u)$ 的多个值 $x$ 和 $u$。

我从原始版本中剥离了风格问题，以专注于这个核心问题。

Answer 1

这是一个很好的例子，您可以在其中充分利用广播。 看起来你想要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

Answer 2

详细说明我上面的评论。乍一看，一般来说，尽量通过预分配数组并就地填充它们(例如使用点广播)来避免临时分配。也可以使用@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)

比较(使用基准工具)性能(和分配)

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

和

f.(nu,x)

(从技术上讲，输出并不相同，您必须使用上面的vcat，但无论如何)

更新(在OP净化了他的代码之后)：

好的，我想我现在(终于)看到了你的真正问题(对此感到抱歉)。我上面所说的是关于优化你的原始代码，关于它如何调用besselj并有效地处理它的输出(请参阅@Matt B.的文章，了解这里关于完整的广播解决方案)。

IIUC，您想利用这一事实(我不知道也没有检查这是否属实)来计算给定nu的besselj，并且x内部有一个对nu的总和.换句话说，您希望使用此内部求和的中间结果来避免冗余计算。

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