<|start_filename|>test/misc.jl<|end_filename|>

@noinline function dostuff(A; B = I)
    StrideArrays.@avx for i in eachindex(A)
        A[i] += 1
    end
    C = A * B
    s = zero(eltype(C))
    StrideArrays.@avx for i in eachindex(C)
        s += C[i]
    end
    s
end
function gc_preserve_test()
    A = @StrideArray rand(8, $(2<<2));
    B = @StrideArray rand($(1<<3), 8);
    @gc_preserve dostuff(A, B = B)
end

@testset "Miscellaneous" begin
    x = @StrideArray rand(127);
    @test maximum(abs, x) == maximum(abs, Array(x))
    y = @StrideArray rand(3);
    @test maximum(abs, y) == maximum(abs, Array(y))
    @test iszero(@allocated gc_preserve_test())

    A = @StrideArray rand($(5<<1),$(1<<3));
    @test StrideArrays.size(A) === (StaticInt(10), StaticInt(8))
    A_u = view(A, StaticInt(1):StaticInt(6), :)
    A_l = view(A, StaticInt(7):StaticInt(10), :)
    @test A == @inferred(vcat(A_u, A_l))
    @test StrideArrays.size(A) === StrideArrays.size(vcat(A_u, A_l))

    # On 1.5 tests fail if you don't do this first.
    @test pointer(view(A, 1:StaticInt(6), :)) == pointer(A)
    A_u = view(A, 1:StaticInt(6), :);
    A_l = view(A, StaticInt(7):StaticInt(10), :);
    @test A == @inferred(vcat(A_u, A_l))
    @test StrideArrays.size(A) == StrideArrays.size(vcat(A_u, A_l))
    @test StrideArrays.size(A) !== StrideArrays.size(vcat(A_u, A_l))

    Aa = Array(A)
    @test sum(A) ≈ sum(Aa)
    @test maximum(A) == maximum(Aa) == maximum(abs, @. -A)
    @test mapreduce(abs2, +, A) ≈ mapreduce(abs2, +, Aa)

    @test A[1,1] == A[1] == vec(A)[1]
    A[1,1] = 3
    @test A[1,1] == A[1] == vec(A)[1] == 3
    A[2] = 4
    @test A[2,1] == A[2] == vec(A)[2] == 4
    vec(A)[3] = 5
    @test A[3,1] == A[3] == vec(A)[3] == 5

    @test ArrayInterface.stride_rank(StrideArrays.similar_layout(A')) === ArrayInterface.stride_rank(StrideArrays.similar_layout(Aa')) === (StaticInt(2),StaticInt(1))
    @test ArrayInterface.contiguous_axis(StrideArrays.similar_layout(A')) === ArrayInterface.contiguous_axis(StrideArrays.similar_layout(Aa')) === StrideArrays.StaticInt{2}()
    
    let B = @StrideArray rand(10,10)
        @test exp(B) ≈ exp(Array(B))
    end
end



<|start_filename|>test/matmul_tests.jl<|end_filename|>

# const TRAVIS_SKIP = VERSION.minor != 4 && !isnothing(get(ENV, "TRAVIS_BRANCH", nothing))

function test_fixed_size(A, At, B, Bt, Aa, Aat, Ba, Bat)::NTuple{4,Float64}
    t0 = @elapsed(@test Aa * Ba ≈ A * B)
    t1 = @elapsed(@test Aa * Bat ≈ A * Bt)
    t2 = @elapsed(@test Aat * Ba ≈ At * B)
    t3 = @elapsed(@test Aat * Bat ≈ At * Bt)
    t0, t1, t2, t3
end
    

function test_fixed_size(M, K, N)
    A = @StrideArray rand(M,K);
    B = @StrideArray rand(K,N);
    At = (@StrideArray rand(K,M))';
    Bt = (@StrideArray rand(N,K))';
    Aa = Array(A); Ba = Array(B);
    Aat = Array(At); Bat = Array(Bt);
    t0 = test_fixed_size(A, At, B, Bt, Aa, Aat, Ba, Bat)
    t1 = test_fixed_size(StrideArrays.make_dynamic(A), StrideArrays.make_dynamic(At), B, Bt, Aa, Aat, Ba, Bat)
    t2 = test_fixed_size(A, At, StrideArrays.make_dynamic(B), StrideArrays.make_dynamic(Bt), Aa, Aat, Ba, Bat)
    t3 = test_fixed_size(StrideArrays.make_dynamic(A), StrideArrays.make_dynamic(At), StrideArrays.make_dynamic(B), StrideArrays.make_dynamic(Bt), Aa, Aat, Ba, Bat)
    gflops = let gflop = 2e-9M*K*N
        map((t0, t1, t2, t3)) do t
            gflop ./ t
        end
    end
    @show (M,K,N), gflops
    nothing
end

test_fixed_size(M) = test_fixed_size(M, M, M)

# gmul(A, B) = LinearAlgebra.generic_matmatmul('N','N', A, B)

@testset "MatrixMultiply" begin
    r = 2:12
    for M ∈ r, K ∈ r, N ∈ r
        test_fixed_size(M, K, N)
    end
    r = 13:33
    for M ∈ r
        test_fixed_size(M)
    end
    M = K = N = 80
    A = @StrideArray rand(M,K);
    B = @StrideArray rand(K,N);
    C = StrideArray{Float64}(undef, (StaticInt(M),StaticInt(N)));
    M, K, N = 23, 37, 19
    @views begin
        Av = A[1:M, 1:K]; 
        Bv = B[1:K, 1:N]; 
        Cv = C[1:M, 1:N]; 
        Avsl = A[StaticInt(1):M, StaticInt(1):K]; 
        Bvsl = B[StaticInt(1):K, StaticInt(1):N]; 
        Cvsl = C[StaticInt(1):M, StaticInt(1):N]; 
        Avsr = A[StaticInt(1):StaticInt(M), StaticInt(1):StaticInt(K)]; 
        Bvsr = B[StaticInt(1):StaticInt(K), StaticInt(1):StaticInt(N)]; 
        Cvsr = C[StaticInt(1):StaticInt(M), StaticInt(1):StaticInt(N)]; 
    end
    Creference = Array(Av) * Array(Bvsl);
    time = @elapsed mul!(Cv, Av, Bv)
    @test Creference ≈ Cv
    @show M, K, N, time
    time = @elapsed mul!(Cvsl, Avsl, Bvsl)
    @test Creference ≈ Cv
    @show M, K, N, time
    time = @elapsed mul!(Cvsr, Avsr, Bvsr)
    @test Creference ≈ Cv
    @show M, K, N, time

    Av = A[1:M, 1:K]; 
    Bv = B[1:K, 1:N]; 
    Cv = C[1:M, 1:N]; 
    Avsl = A[StaticInt(1):M, StaticInt(1):K]; 
    Bvsl = B[StaticInt(1):K, StaticInt(1):N]; 
    Cvsl = C[StaticInt(1):M, StaticInt(1):N]; 
    Avsr = A[StaticInt(1):StaticInt(M), StaticInt(1):StaticInt(K)]; 
    Bvsr = B[StaticInt(1):StaticInt(K), StaticInt(1):StaticInt(N)]; 
    Cvsr = C[StaticInt(1):StaticInt(M), StaticInt(1):StaticInt(N)]; 
    Creference = Array(Av) * Array(Bvsl);
    time = @elapsed mul!(Cv, Av, Bv)
    @test Creference ≈ Cv
    @show M, K, N, time
    time = @elapsed mul!(Cvsl, Avsl, Bvsl)
    @test Creference ≈ Cv
    @show M, K, N, time
    time = @elapsed mul!(Cvsr, Avsr, Bvsr)
    @test Creference ≈ Cv
    @show M, K, N, time
end




<|start_filename|>docs/make.jl<|end_filename|>
using StrideArrays
using Documenter

makedocs(;
    modules=[StrideArrays],
    authors="<NAME>",
    repo="https://github.com/chriselrod/StrideArrays.jl/blob/{commit}{path}#L{line}",
    sitename="StrideArrays.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://chriselrod.github.io/StrideArrays.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Getting Started" => "getting_started.md",
        "Architecture Benchmarks" => [
            "arches/cascadelake.md",
            "arches/tigerlake.md",
            "arches/haswell.md"
        ],
        "Random Number Generation" => "rng.md",
        "Broadcasting" => "broadcasting.md",
        "Stack Allocattion" => "stack_allocation.md"        
    ],
    strict=false,
)

deploydocs(;
    repo="github.com/chriselrod/StrideArrays.jl",
)


<|start_filename|>src/blas.jl<|end_filename|>

@inline LinearAlgebra.mul!(C::AbstractStrideMatrix, A::StridedMatrix, B::StridedMatrix) = (@gc_preserve(matmul!(C, A, B)); return C)
@inline LinearAlgebra.mul!(C::StridedMatrix, A::AbstractStrideMatrix, B::StridedMatrix) = (@gc_preserve(matmul!(C, A, B)); return C)
@inline LinearAlgebra.mul!(C::StridedMatrix, A::StridedMatrix, B::AbstractStrideMatrix) = (@gc_preserve(matmul!(C, A, B)); return C)
@inline LinearAlgebra.mul!(C::AbstractStrideMatrix, A::AbstractStrideMatrix, B::StridedMatrix) = (@gc_preserve(matmul!(C, A, B)); return C)
@inline LinearAlgebra.mul!(C::AbstractStrideMatrix, A::StridedMatrix, B::AbstractStrideMatrix) = (@gc_preserve(matmul!(C, A, B)); return C)
@inline LinearAlgebra.mul!(C::StridedMatrix, A::AbstractStrideMatrix, B::AbstractStrideMatrix) = (@gc_preserve(matmul!(C, A, B)); return C)
@inline LinearAlgebra.mul!(C::AbstractStrideMatrix, A::AbstractStrideMatrix, B::AbstractStrideMatrix) = (@gc_preserve(matmul!(C, A, B)); return C)

@inline function stridematmul(
    A::AbstractStrideMatrix{<:Any,<:Any,TA},
    B::AbstractStrideMatrix{<:Any,<:Any,TB}
) where {TA<:Base.HWReal,TB<:Base.HWReal}
    M, KA = size(A)
    KB, N = size(B)
    @assert KA == KB "Size mismatch."
    K = Octavian._select(KA, KB)
    TC = promote_type(TA,TB)
    C = StrideArray{TC}(undef, (M, N))
    @gc_preserve matmul!(C, A, B, One(), Zero(), nothing, (M, K, N), ArrayInterface.contiguous_axis(C))
    return C
end
@inline Base.:*(A::AbstractStrideMatrix{<:Any,<:Any,TA}, B::StridedMatrix{TB}) where {TA<:Base.HWReal,TB<:Base.HWReal} = @gc_preserve(stridematmul(A, B))
@inline Base.:*(A::AbstractStrideMatrix{<:Any,<:Any,TA}, B::StridedMatrix{TB}) where {TA<:Union{Float32,Float64},TB<:Union{Float32,Float64}} = @gc_preserve(stridematmul(A, B))
@inline Base.:*(A::StridedMatrix{TA}, B::AbstractStrideMatrix{<:Any,<:Any,TB}) where {TA<:Base.HWReal,TB<:Base.HWReal} = @gc_preserve(stridematmul(A, B))
@inline Base.:*(A::StridedMatrix{TA}, B::AbstractStrideMatrix{<:Any,<:Any,TB}) where {TA<:Union{Float32,Float64},TB<:Union{Float32,Float64}} = @gc_preserve(stridematmul(A, B))
@inline Base.:*(A::AbstractStrideMatrix{<:Any,<:Any,TA}, B::AbstractStrideMatrix{<:Any,<:Any,TB}) where {TA<:Base.HWReal,TB<:Base.HWReal} = @gc_preserve(stridematmul(A, B))
@inline Base.:*(A::AbstractStrideMatrix{<:Any,<:Any,TA}, B::AbstractStrideMatrix{<:Any,<:Any,TB}) where {TA<:Union{Float32,Float64},TB<:Union{Float32,Float64}} = @gc_preserve(stridematmul(A, B))

@inline extract_λ(a) = a
@inline extract_λ(a::UniformScaling) = a.λ
@inline function Base.:*(A::AbstractStrideArray{S,D,T}, bλ::Union{Tb,UniformScaling{Tb}}) where {S,D,T<:VectorizationBase.NativeTypes,Tb <: Real}
    mv = similar(A)
    b = T(extract_λ(bλ))
    @avx for i ∈ eachindex(A)
        mv[i] = A[i] * b
    end
    mv
end

function LinearAlgebra.mul!(
    C::AbstractStrideMatrix{<:Any,<:Any,T},
    A::LinearAlgebra.Diagonal{T,<:AbstractStrideVector{<:Any,T}},
    B::StridedMatrix{T}
) where {T}
    M, K, N = matmul_axes(C, A, B)
    MandK = ArrayInterface._pick_range(M, K)
    vA = parent(A)
    @avx for n ∈ N, m ∈ 1:MandK
        C[m,n] = vA[m] * B[m,n]
    end
    C
end
@inline function Base.:*(
    A::LinearAlgebra.Diagonal{T,<:AbstractVector{T}},
    B::AbstractStrideMatrix{<:Any,<:Any,T}
) where {T}
    mul!(similar(B), A, B)
end


<|start_filename|>benchmark/staticarraybench.jl<|end_filename|>

using StrideArrays, StaticArrays, LinearAlgebra, BenchmarkTools

# BLAS.set_num_threads(1)

# For laptops that thermally throttle, you can set the `JULIA_SLEEP_BENCH` environment variable for #seconds to sleep before each `@belapsed`
const SLEEPTIME = parse(Float64, get(ENV, "JULIA_SLEEP_BENCH", "0"))
maybe_sleep() = iszero(SLEEPTIME) || sleep(SLEEPTIME)
# BenchmarkTools.DEFAULT_PARAMETERS.samples = 1_000_000
# BenchmarkTools.DEFAULT_PARAMETERS.seconds = 10

matrix_sizes(x::Integer) = (x,x,x)
matrix_sizes(x::NTuple{3}) = x

function runbenches(sr, ::Type{T} = Float64) where {T}
    bench_results = Matrix{T}(undef, length(sr), 5)
    for (i,s) ∈ enumerate(sr)
        M, N, K = matrix_sizes(s)
        if true#s ≤ 20
            Astatic = @SMatrix rand(T, M, K);
            Bstatic = @SMatrix rand(T, K, N);
            maybe_sleep()
            bench_results[i,1] = @belapsed $(Ref(Astatic))[] * $(Ref(Bstatic))[]
            Amutable = MArray(Astatic);
            Bmutable = MArray(Bstatic);
            Cmutable = similar(Amutable);
        else
            bench_results[i,1] = Inf
            Amutable = MMatrix{M,K,T}(undef)
            Bmutable = MMatrix{K,N,T}(undef)
            Cmutable = MMatrix{M,N,T}(undef)
            @inbounds for i ∈ 1:s^2
                Amutable[i] = rand()
                Bmutable[i] = rand()
            end
        end
        maybe_sleep()
        bench_results[i,2] = @belapsed mul!($Cmutable, $Amutable, $Bmutable)
        Afixed = StrideArray(Amutable)
        Bfixed = StrideArray(Bmutable)
        Cfixed = StrideArray{T}(undef, (StaticInt(M),StaticInt(N)))
        maybe_sleep()
        bench_results[i,3] = @belapsed mul!($Cfixed, $Afixed, $Bfixed)
        Cfixed2 = similar(Cfixed);
        Aptr = PtrArray(Afixed); Bptr = PtrArray(Bfixed); Cptr = PtrArray(Cfixed2);
        maybe_sleep()
        GC.@preserve Afixed Bfixed Cfixed2 begin
            bench_results[i,4] = @belapsed mul!($Cptr, $Aptr, $Bptr)
        end
        A = Array(Afixed); B = Array(Bfixed); C = Matrix{T}(undef, M, N);
        maybe_sleep()
        bench_results[i,5] = @belapsed jmul!($C, $A, $B)
        # @show Array(Cmutable) Cfixed Cfixed2 C
        @assert Array(Cmutable) ≈ Cfixed ≈ Cfixed2 ≈ C
        v = @view(bench_results[i,:])'
        @show s, v
    end
    bench_results
end

sizerange = 2:48
br = runbenches(sizerange);
using DataFrames, VegaLite

gflops = @. 2e-9 * (sizerange) ^ 3 / br;

df = DataFrame(gflops);
matmulmethodnames = [:SMatrix, :MMatrix, :StrideArray, :PtrArray, :matmul!];
rename!(df, matmulmethodnames);
df.Size = sizerange

function pick_suffix(desc = "")
    suffix = if StrideArrays.VectorizationBase.AVX512F
        "AVX512"
    elseif StrideArrays.VectorizationBase.AVX2
        "AVX2"
    elseif StrideArrays.VectorizationBase.REGISTER_SIZE == 32
        "AVX"
    else
        "REGSUZE$(StrideArrays.VectorizationBase.REGISTER_SIZE)"
    end
    if desc != ""
        suffix *= '_' * desc
    end
    "$(Sys.CPU_NAME)_$suffix"
end

dfs = stack(df, matmulmethodnames, variable_name = :MatMulType, value_name = :GFLOPS);
p = dfs |> @vlplot(:line, x = :Size, y = :GFLOPS, width = 900, height = 600, color = {:MatMulType});
save(joinpath(pkgdir(StrideArrays), "docs/src/assets/sizedarraybenchmarks_$(pick_suffix()).svg"), p)



