特殊 (besselj) 函数矩阵
matrix of special (besselj) functions
我是 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)$。
我已经从原始版本中剥离了文体问题,以专注于这个核心问题。
详细说明我上面的评论。乍一看,一般来说,通过预先分配数组并就地填充它们(例如使用 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 实现)。
这是一个很好的例子,您可以在其中充分利用广播。看起来您想要 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
我是 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)$。
我已经从原始版本中剥离了文体问题,以专注于这个核心问题。
详细说明我上面的评论。乍一看,一般来说,通过预先分配数组并就地填充它们(例如使用 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 实现)。
这是一个很好的例子,您可以在其中充分利用广播。看起来您想要 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