Hi! Suppose I want to code a function that returns a `SVector` but each component

General Usage
1

Hi! Suppose I want to code a function that returns a SVector but each component is computed inside the function without knowing the size. Should I create a MVector and then cast that to a SVector? This is what I am thinking:

function test(::Val{D}) where {D}
    mv = @MVector zeros(D)

    for i in 1:D
        mv[i] = i + 1
    end

    sv = convert(SVector, mv)
    return sv
end

Then:

test(Val{4}())

does not allocate and is type stable. But I don’t know if this is the correct approach. Thanks!

Note that the original poster on Slack cannot see your response here on Discourse. Consider transcribing the appropriate answer back to Slack, or pinging the poster here on Discourse so they can follow this thread.
(Original message :slack:) (More Info)

2

It seems that the function was not translated correctly from Slack:

function test(::Val{D}) where {D}
    mv = @MVector zeros(D)
    for i in 1:D
        mv[i] = i + 1
    end
    sv = convert(SVector, mv)
    return sv
end

Then:

test(Val{4}())
1 Like
3

Fixed :slight_smile:

1 Like
4

This is shorter, simpler, and probably faster:

test2(n::Val) = SVector(ntuple(i -> float(i + 1), n))

The key is that an SVector can be constructed from a tuple, and a tuple can be constructed efficiently using ntuple with a Val argument so that the size is known at compile time.

4 Likes
5

Hi @stevengj, thanks for your reply!

I have a more complicated case, where each value in the vector depends on more complicated expressions. This is the actual function that I am coding:

@inline function drift(u, p::LMMP{false,D,D,true,Terminal}, t) where {D}
    @unpack Tenors, τ, ρ = p
    σ = p.σ(t)

    du = @MVector zeros(eltype(u), D)

    @inbounds begin
        # 'i' ranges from 2 to D-1 because L₁ and LN have μ = 0
        for i in 2:D-1
            if t ≤ Tenors[i]
                for j in i+1:D
                    du[i] += (ρ[i,j] * τ[j] * σ[j] * u[j]) / (1 + τ[j] * u[j])
                end
                du[i] *= (-σ[i] * u[i])
            end
        end
    end

    return convert(SVector, m)
end

Should I also implement the tuple strategy?

1 Like
6

You could do

return SVector(ntuple(Val{D}()) do i
    du = zero(eltype(u))
    if t ≤ Tenors[i]
        ... same as above...
    end
    du
end)
2 Likes
7

This is great, thank you so much! I will try it!

1 Like
8

I have tested both and using MVector is faster… This is a MWE:

using BenchmarkTools
using StaticArrays
using UnPack

N = 10
τ = @SVector [0.25 for i in 1:N]
Tenors = vcat(zero(eltype(τ)), cumsum(τ))
σ(t) = @SVector ones(10)
ρ = @SMatrix rand(N,N)
u = @SVector ones(N)
p = (τ = τ, Tenors = Tenors, σ = σ, ρ = ρ)
t = 0.1
val = Val{N}()

@inline function driftlmm1(u, p, t, ::Val{D}) where {D}
    @unpack Tenors, τ, ρ = p
    σ = p.σ(t)

    du = @MVector zeros(eltype(u), D)

    @inbounds begin
        # 'i' ranges from 2 to D-1 because L₁ and LD have μ = 0
        for i in 2:D-1
            if t ≤ Tenors[i]
                for j in i+1:D
                    du[i] += (ρ[i,j] * τ[j] * σ[j] * u[j]) / (1 + τ[j] * u[j])
                end
                du[i] *= (-σ[i] * u[i])
            end
        end
    end

    return convert(SVector, du)
end

@btime driftlmm1($u, $p, $t, $val)
# 152.777 ns (0 allocations: 0 bytes)

@inline function driftlmm2(u, p, t, ::Val{D}) where {D}
    @unpack Tenors, τ, ρ = p
    σ = p.σ(t)

    @inbounds begin
        return SVector(ntuple(Val{D}()) do i
            du = zero(eltype(u))
            if t ≤ Tenors[i]
                for j in i+1:D
                    du += (ρ[i,j] * τ[j] * σ[j] * u[j]) / (1 + τ[j] * u[j])
                end
                du *= (-σ[i] * u[i])
            end
            du
        end)
    end
end

@btime driftlmm2($u, $p, $t, $val)
# 178.298 ns (0 allocations: 0 bytes)

Any ideas?

9

I’m not sure if the @inbounds propagates to the inside of the do block.

1 Like
10

I have removed @inbounds and @inline in both cases and the one using MVector is still faster.

Any way, the method that you show me is great.

I was wondering if there is a faster way!