From f94dbac514db2bf77083efc1b5920c3583fa03de Mon Sep 17 00:00:00 2001 From: Anton Pozharskiy Date: Thu, 11 Jun 2026 15:36:58 +0200 Subject: [PATCH 1/4] Constructors for dmat, smat, dvec, and svec --- Project.toml | 5 +++- src/BLASFEO.jl | 9 +++++- src/mats.jl | 78 ++++++++++++++++++++++++++++++++++++++++++++++++++ src/show.jl | 21 ++++++++++++++ src/vecs.jl | 70 ++++++++++++++++++++++++++++++++++++++++++++ 5 files changed, 181 insertions(+), 2 deletions(-) create mode 100644 src/mats.jl create mode 100644 src/show.jl create mode 100644 src/vecs.jl diff --git a/Project.toml b/Project.toml index b901066..8711a32 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,10 @@ name = "BLASFEO" uuid = "2d913c56-f67c-11e9-3afb-9119491b9f4c" -authors = ["Ian McInerney", "Imperial College London"] version = "0.1.0" +authors = ["Ian McInerney", "Imperial College London"] + +[deps] +blasfeo_jll = "6b574d4a-bb57-5a4b-b7e6-1c794c903646" [compat] julia = "1" diff --git a/src/BLASFEO.jl b/src/BLASFEO.jl index e7d76f7..3fe7bae 100644 --- a/src/BLASFEO.jl +++ b/src/BLASFEO.jl @@ -1,5 +1,12 @@ module BLASFEO +using blasfeo_jll +include("vecs.jl") +include("mats.jl") +include("show.jl") -greet() = print("Hello World!") +# Vectors +export BlasfeoDvec, BlasfeoSvec +# Matrices +export BlasfeoDmat, BlasfeoSmat end # module diff --git a/src/mats.jl b/src/mats.jl new file mode 100644 index 0000000..4a73171 --- /dev/null +++ b/src/mats.jl @@ -0,0 +1,78 @@ +# NOTE(@anton): Both of these structs are treated as mutable and Julia manages their +# lifetimes itself. Blasfeo also offers a way to preallocate memory but +# we still need the bitstype struct and to manage that memory so I think +# may as well use `blasfeo_allocate_*` and `blasfeo_free_*` + +# bits clone of panel major `blasfeo_dmat` +# blasfeo_jll compiles with PANELMAJ so we only need this one and not the column major version +mutable struct BlasfeoDmat + const mem::Ptr{Cdouble} # pointer to passed chunk of memory + const pA::Ptr{Cdouble} # pointer to a pm*pn array of doubles, the first is aligned to cache line size + const dA::Ptr{Cdouble} # pointer to a min(m,n) (or max???) array of doubles + const m::Cint # rows + const n::Cint # cols + const pm::Cint # packed number or rows + const cn::Cint # packed number or cols + const use_dA::Cint # flag to tell if dA can be used + const memsize::Cint # size of needed memory + + function BlasfeoDmat(m::Int,n::Int) + mat = new(C_NULL,C_NULL,C_NULL,0,0,0,0,0,0) + @ccall blasfeo.blasfeo_allocate_dmat(m::Cint, n::Cint, pointer_from_objref(mat)::Ptr{BlasfeoDmat})::Cvoid + function destructor(this) + @ccall blasfeo.blasfeo_free_dmat(pointer_from_objref(this)::Ptr{BlasfeoDmat})::Cvoid + end + return finalizer(destructor, mat) + end + + function BlasfeoDmat(other::Matrix{Cdouble}) + mat = new(C_NULL,C_NULL,C_NULL,0,0,0,0,0,0) + m,n = size(other) + @ccall blasfeo.blasfeo_allocate_dmat(m::Cint, n::Cint, pointer_from_objref(mat)::Ptr{BlasfeoDmat})::Cvoid + function destructor(this) + @ccall blasfeo.blasfeo_free_dmat(pointer_from_objref(this)::Ptr{BlasfeoDmat})::Cvoid + end + @ccall blasfeo.blasfeo_pack_dmat(m::Cint, n::Cint, + other::Ptr{Cdouble}, m::Cint, + pointer_from_objref(mat)::Ptr{BlasfeoDmat}, + 0::Cint, 0::Cint)::Cvoid + return finalizer(destructor, mat) + end +end + +# bits clone of panel major `blasfeo_smat` +# blasfeo_jll compiles with PANELMAJ so we only need this one and not the column major version +mutable struct BlasfeoSmat + const mem::Ptr{Cfloat} # pointer to passed chunk of memory + const pA::Ptr{Cfloat} # pointer to a pm*pn array of doubles, the first is aligned to cache line size + const dA::Ptr{Cfloat} # pointer to a min(m,n) (or max???) array of doubles + const m::Cint # rows + const n::Cint # cols + const pm::Cint # packed number or rows + const cn::Cint # packed number or cols + const use_dA::Cint # flag to tell if dA can be used + const memsize::Cint # size of needed memory + + function BlasfeoSmat(m::Int,n::Int) + mat = new(C_NULL,C_NULL,C_NULL,0,0,0,0,0,0) + @ccall blasfeo.blasfeo_allocate_smat(m::Cint, n::Cint, pointer_from_objref(mat)::Ptr{BlasfeoSmat})::Cvoid + function destructor(this) + @ccall blasfeo.blasfeo_free_smat(pointer_from_objref(this)::Ptr{BlasfeoSmat})::Cvoid + end + return finalizer(destructor, mat) + end + + function BlasfeoSmat(other::Matrix{Cfloat}) + mat = new(C_NULL,C_NULL,C_NULL,0,0,0,0,0,0) + m,n = size(other) + @ccall blasfeo.blasfeo_allocate_smat(m::Cint, n::Cint, pointer_from_objref(mat)::Ptr{BlasfeoSmat})::Cvoid + function destructor(this) + @ccall blasfeo.blasfeo_free_smat(pointer_from_objref(this)::Ptr{BlasfeoSmat})::Cvoid + end + @ccall blasfeo.blasfeo_pack_smat(m::Cint, n::Cint, + other::Ptr{Cdouble}, m::Cint, + pointer_from_objref(mat)::Ptr{BlasfeoSmat}, + 0::Cint, 0::Cint)::Cvoid + return finalizer(destructor, mat) + end +end diff --git a/src/show.jl b/src/show.jl new file mode 100644 index 0000000..91eb527 --- /dev/null +++ b/src/show.jl @@ -0,0 +1,21 @@ +# TODO(@anton) mxn or nxm + +function Base.show(io::IO, mat::BlasfeoDmat) + println("$(mat.m)x$(mat.n) BlasfeoDmat:") + @ccall blasfeo.blasfeo_print_dmat(mat.m::Cint, mat.n::Cint, pointer_from_objref(mat)::Ptr{BlasfeoDmat}, 0::Cint, 0::Cint)::Cvoid +end + +function Base.show(io::IO, mat::BlasfeoSmat) + println("$(mat.m)x$(mat.n) BlasfeoSmat:") + @ccall blasfeo.blasfeo_print_smat(mat.m::Cint, mat.n::Cint, pointer_from_objref(mat)::Ptr{BlasfeoSmat}, 0::Cint, 0::Cint)::Cvoid +end + +function Base.show(io::IO, vec::BlasfeoDvec) + println("$(vec.m)-element BlasfeoDvec:") + @ccall blasfeo.blasfeo_print_dvec(vec.m::Cint, pointer_from_objref(vec)::Ptr{BlasfeoDvec}, 0::Cint)::Cvoid +end + +function Base.show(io::IO, vec::BlasfeoSvec) + println("$(vec.m)-element BlasfeoSvec:") + @ccall blasfeo.blasfeo_print_svec(vec.m::Cint, pointer_from_objref(vec)::Ptr{BlasfeoSvec}, 0::Cint)::Cvoid +end diff --git a/src/vecs.jl b/src/vecs.jl new file mode 100644 index 0000000..771809b --- /dev/null +++ b/src/vecs.jl @@ -0,0 +1,70 @@ +# NOTE(@anton): Both of these structs are treated as mutable and Julia manages their +# lifetimes itself. Blasfeo also offers a way to preallocate memory but +# we still need the bitstype struct and to manage that memory so I think +# may as well use `blasfeo_allocate_*` and `blasfeo_free_*` + +# bits clone of panel major `blasfeo_dvec` +mutable struct BlasfeoDvec + mem::Ptr{Cdouble} # pointer to passed chunk of memory + pa::Ptr{Cdouble} # pointer to a pm array of doubles, the first is aligned to cache line size + m::Cint # size + pm::Cint # packed size + memsize::Cint # size of needed memory + + function BlasfeoDvec(m::Int) + vec = new(C_NULL,C_NULL,0,0,0) + @ccall blasfeo.blasfeo_allocate_dvec(m::Cint, pointer_from_objref(vec)::Ptr{BlasfeoDvec})::Cvoid + function destructor(this) + @ccall blasfeo.blasfeo_free_dvec(pointer_from_objref(this)::Ptr{BlasfeoDvec})::Cvoid + end + return finalizer(destructor, vec) + end + + function BlasfeoDvec(other::Vector{Cdouble}) + vec = new(C_NULL,C_NULL,0,0,0) + @ccall blasfeo.blasfeo_allocate_dvec(m::Cint, pointer_from_objref(vec)::Ptr{BlasfeoDvec})::Cvoid + function destructor(this) + @ccall blasfeo.blasfeo_free_dvec(pointer_from_objref(this)::Ptr{BlasfeoDvec})::Cvoid + end + + @ccall blasfeo.blasfeo_pack_dvec(m::Cint, + other::Ptr{Cdouble}, 1::Cint, + pointer_from_objref(vec)::Ptr{BlasfeoDvec}, + 0::Cint)::Cvoid + + return finalizer(destructor, vec) + end +end + +# bits clone of panel major `blasfeo_svec` +mutable struct BlasfeoSvec + mem::Ptr{Cfloat} # pointer to passed chunk of memory + pa::Ptr{Cfloat} # pointer to a pm array of floats, the first is aligned to cache line size + m::Cint # size + pm::Cint # packed size + memsize::Cint # size of needed memory + + function BlasfeoSvec(m::Int,n::Int) + vec = new(C_NULL,C_NULL,0,0,0) + @ccall blasfeo.blasfeo_allocate_svec(m::Cint, pointer_from_objref(vec)::Ptr{BlasfeoSvec})::Cvoid + function destructor(this) + @ccall blasfeo.blasfeo_free_svec(pointer_from_objref(this)::Ptr{BlasfeoSvec})::Cvoid + end + return finalizer(destructor, vec) + end + + function BlasfeoSvec(other::Vector{Cfloat}) + vec = new(C_NULL,C_NULL,0,0,0) + @ccall blasfeo.blasfeo_allocate_svec(m::Cint, pointer_from_objref(vec)::Ptr{BlasfeoSvec})::Cvoid + function destructor(this) + @ccall blasfeo.blasfeo_free_svec(pointer_from_objref(this)::Ptr{BlasfeoSvec})::Cvoid + end + + @ccall blasfeo.blasfeo_pack_svec(m::Cint, + other::Ptr{Cfloat}, 1::Cint, + pointer_from_objref(vec)::Ptr{BlasfeoSvec}, + 0::Cint)::Cvoid + + return finalizer(destructor, vec) + end +end From c9c6438d91c892b3a88118fac8c74490d9ce2f39 Mon Sep 17 00:00:00 2001 From: Anton Pozharskiy Date: Thu, 11 Jun 2026 18:00:44 +0200 Subject: [PATCH 2/4] working show and beginnings of basic AbstractMatrix/Vector work --- src/mats.jl | 22 ++++++++++++++++++---- src/show.jl | 45 +++++++++++++++++++++++++++++++++++++++++---- src/vecs.jl | 12 ++++++------ 3 files changed, 65 insertions(+), 14 deletions(-) diff --git a/src/mats.jl b/src/mats.jl index 4a73171..89a084c 100644 --- a/src/mats.jl +++ b/src/mats.jl @@ -5,7 +5,7 @@ # bits clone of panel major `blasfeo_dmat` # blasfeo_jll compiles with PANELMAJ so we only need this one and not the column major version -mutable struct BlasfeoDmat +mutable struct BlasfeoDmat <: AbstractMatrix{Cdouble} const mem::Ptr{Cdouble} # pointer to passed chunk of memory const pA::Ptr{Cdouble} # pointer to a pm*pn array of doubles, the first is aligned to cache line size const dA::Ptr{Cdouble} # pointer to a min(m,n) (or max???) array of doubles @@ -16,7 +16,7 @@ mutable struct BlasfeoDmat const use_dA::Cint # flag to tell if dA can be used const memsize::Cint # size of needed memory - function BlasfeoDmat(m::Int,n::Int) + function BlasfeoDmat(m::Integer,n::Integer) mat = new(C_NULL,C_NULL,C_NULL,0,0,0,0,0,0) @ccall blasfeo.blasfeo_allocate_dmat(m::Cint, n::Cint, pointer_from_objref(mat)::Ptr{BlasfeoDmat})::Cvoid function destructor(this) @@ -42,7 +42,7 @@ end # bits clone of panel major `blasfeo_smat` # blasfeo_jll compiles with PANELMAJ so we only need this one and not the column major version -mutable struct BlasfeoSmat +mutable struct BlasfeoSmat <: AbstractMatrix{Cfloat} const mem::Ptr{Cfloat} # pointer to passed chunk of memory const pA::Ptr{Cfloat} # pointer to a pm*pn array of doubles, the first is aligned to cache line size const dA::Ptr{Cfloat} # pointer to a min(m,n) (or max???) array of doubles @@ -53,7 +53,7 @@ mutable struct BlasfeoSmat const use_dA::Cint # flag to tell if dA can be used const memsize::Cint # size of needed memory - function BlasfeoSmat(m::Int,n::Int) + function BlasfeoSmat(m::Integer,n::Integer) mat = new(C_NULL,C_NULL,C_NULL,0,0,0,0,0,0) @ccall blasfeo.blasfeo_allocate_smat(m::Cint, n::Cint, pointer_from_objref(mat)::Ptr{BlasfeoSmat})::Cvoid function destructor(this) @@ -76,3 +76,17 @@ mutable struct BlasfeoSmat return finalizer(destructor, mat) end end + +# basic matrix operations +# TODO(@anton) Do we _need_ such tight typing on setindex? +for (type,flag) in [ + (:BlasfeoDmat, :d), + (:BlasfeoSmat, :s), + ] + @eval Base.size(A::$type) = (A.m, A.n) + blasfeo_geex1 = Symbol(:blasfeo_, flag, :geex1) + @eval Base.getindex(A::$type, I::Vararg{Int, 2}) = @ccall blasfeo.$blasfeo_geex1(pointer_from_objref(A)::Ptr{$type}, (I[1]-1)::Cint, (I[2]-1)::Cint)::eltype($type) + @eval Base.similar(A::$type) = $type(A.m, A.n) + @eval Base.similar(A::$type,dims::Dims{2}) = $type(dims...) + +end diff --git a/src/show.jl b/src/show.jl index 91eb527..add7b1f 100644 --- a/src/show.jl +++ b/src/show.jl @@ -1,21 +1,58 @@ # TODO(@anton) mxn or nxm -function Base.show(io::IO, mat::BlasfeoDmat) +function Base.show(io, mat::BlasfeoDmat) println("$(mat.m)x$(mat.n) BlasfeoDmat:") @ccall blasfeo.blasfeo_print_dmat(mat.m::Cint, mat.n::Cint, pointer_from_objref(mat)::Ptr{BlasfeoDmat}, 0::Cint, 0::Cint)::Cvoid end -function Base.show(io::IO, mat::BlasfeoSmat) +function Base.show(io, mat::BlasfeoSmat) println("$(mat.m)x$(mat.n) BlasfeoSmat:") @ccall blasfeo.blasfeo_print_smat(mat.m::Cint, mat.n::Cint, pointer_from_objref(mat)::Ptr{BlasfeoSmat}, 0::Cint, 0::Cint)::Cvoid end -function Base.show(io::IO, vec::BlasfeoDvec) +function Base.show(io, vec::BlasfeoDvec) println("$(vec.m)-element BlasfeoDvec:") @ccall blasfeo.blasfeo_print_dvec(vec.m::Cint, pointer_from_objref(vec)::Ptr{BlasfeoDvec}, 0::Cint)::Cvoid end -function Base.show(io::IO, vec::BlasfeoSvec) +function Base.show(io, vec::BlasfeoSvec) println("$(vec.m)-element BlasfeoSvec:") @ccall blasfeo.blasfeo_print_svec(vec.m::Cint, pointer_from_objref(vec)::Ptr{BlasfeoSvec}, 0::Cint)::Cvoid end + +# TODO(@anton) This is probably not needed, we just need `print_matrix`? this is horribly documented. +for (type,shortname) in [ + (:BlasfeoDvec, :dvec), + (:BlasfeoSvec, :svec), + ] + printer = Symbol(:blasfeo_print_,shortname) + @eval begin + function Base.show(io::IO, ::MIME"text/plain", vec::$type) + println("$(vec.m)-element $($type):") + @ccall blasfeo.$printer(vec.m::Cint, pointer_from_objref(vec)::Ptr{$type}, 0::Cint)::Cvoid + end + + function Base.show(io::IO, vec::$type) + println("$(vec.m)-element $($type):") + @ccall blasfeo.$printer(vec.m::Cint, pointer_from_objref(vec)::Ptr{$type}, 0::Cint)::Cvoid + end + end +end + +for (type,shortname) in [ + (:BlasfeoDmat, "dmat"), + (:BlasfeoSmat, "smat"), + ] + printer = Symbol(:blasfeo_print_,shortname) + @eval begin + function Base.show(io::IO, ::MIME"text/plain", mat::$type) + println("$(mat.m)x$(mat.n) $($type):") + @ccall blasfeo.$printer(mat.m::Cint, mat.n::Cint, pointer_from_objref(mat)::Ptr{$type}, 0::Cint, 0::Cint)::Cvoid + end + + function Base.show(io::IO, mat::$type) + println("$(mat.m)x$(mat.n) $($type):") + @ccall blasfeo.$printer(mat.m::Cint, mat.n::Cint, pointer_from_objref(mat)::Ptr{$type}, 0::Cint, 0::Cint)::Cvoid + end + end +end diff --git a/src/vecs.jl b/src/vecs.jl index 771809b..4fa1356 100644 --- a/src/vecs.jl +++ b/src/vecs.jl @@ -4,7 +4,7 @@ # may as well use `blasfeo_allocate_*` and `blasfeo_free_*` # bits clone of panel major `blasfeo_dvec` -mutable struct BlasfeoDvec +mutable struct BlasfeoDvec <: AbstractVector{Cdouble} mem::Ptr{Cdouble} # pointer to passed chunk of memory pa::Ptr{Cdouble} # pointer to a pm array of doubles, the first is aligned to cache line size m::Cint # size @@ -22,12 +22,12 @@ mutable struct BlasfeoDvec function BlasfeoDvec(other::Vector{Cdouble}) vec = new(C_NULL,C_NULL,0,0,0) - @ccall blasfeo.blasfeo_allocate_dvec(m::Cint, pointer_from_objref(vec)::Ptr{BlasfeoDvec})::Cvoid + @ccall blasfeo.blasfeo_allocate_dvec(length(other)::Cint, pointer_from_objref(vec)::Ptr{BlasfeoDvec})::Cvoid function destructor(this) @ccall blasfeo.blasfeo_free_dvec(pointer_from_objref(this)::Ptr{BlasfeoDvec})::Cvoid end - @ccall blasfeo.blasfeo_pack_dvec(m::Cint, + @ccall blasfeo.blasfeo_pack_dvec(vec.m::Cint, other::Ptr{Cdouble}, 1::Cint, pointer_from_objref(vec)::Ptr{BlasfeoDvec}, 0::Cint)::Cvoid @@ -37,7 +37,7 @@ mutable struct BlasfeoDvec end # bits clone of panel major `blasfeo_svec` -mutable struct BlasfeoSvec +mutable struct BlasfeoSvec <: AbstractVector{Cfloat} mem::Ptr{Cfloat} # pointer to passed chunk of memory pa::Ptr{Cfloat} # pointer to a pm array of floats, the first is aligned to cache line size m::Cint # size @@ -55,12 +55,12 @@ mutable struct BlasfeoSvec function BlasfeoSvec(other::Vector{Cfloat}) vec = new(C_NULL,C_NULL,0,0,0) - @ccall blasfeo.blasfeo_allocate_svec(m::Cint, pointer_from_objref(vec)::Ptr{BlasfeoSvec})::Cvoid + @ccall blasfeo.blasfeo_allocate_svec(length(other)::Cint, pointer_from_objref(vec)::Ptr{BlasfeoSvec})::Cvoid function destructor(this) @ccall blasfeo.blasfeo_free_svec(pointer_from_objref(this)::Ptr{BlasfeoSvec})::Cvoid end - @ccall blasfeo.blasfeo_pack_svec(m::Cint, + @ccall blasfeo.blasfeo_pack_svec(vec.m::Cint, other::Ptr{Cfloat}, 1::Cint, pointer_from_objref(vec)::Ptr{BlasfeoSvec}, 0::Cint)::Cvoid From 182d587a35ca233a579e0e3e48ff101c1721a5b4 Mon Sep 17 00:00:00 2001 From: Anton Pozharskiy Date: Fri, 12 Jun 2026 08:11:08 +0200 Subject: [PATCH 3/4] basic operations and testing basic operations --- src/mats.jl | 34 ++++++++++++++++++++++++++++++- src/vecs.jl | 50 ++++++++++++++++++++++++++++++++++++++++++++-- test/mat/basics.jl | 27 +++++++++++++++++++++++++ test/runtests.jl | 5 ++--- test/vec/basics.jl | 27 +++++++++++++++++++++++++ 5 files changed, 137 insertions(+), 6 deletions(-) create mode 100644 test/mat/basics.jl create mode 100644 test/vec/basics.jl diff --git a/src/mats.jl b/src/mats.jl index 89a084c..36c6fb5 100644 --- a/src/mats.jl +++ b/src/mats.jl @@ -78,15 +78,47 @@ mutable struct BlasfeoSmat <: AbstractMatrix{Cfloat} end # basic matrix operations -# TODO(@anton) Do we _need_ such tight typing on setindex? for (type,flag) in [ (:BlasfeoDmat, :d), (:BlasfeoSmat, :s), ] + # size @eval Base.size(A::$type) = (A.m, A.n) + + # getindex blasfeo_geex1 = Symbol(:blasfeo_, flag, :geex1) @eval Base.getindex(A::$type, I::Vararg{Int, 2}) = @ccall blasfeo.$blasfeo_geex1(pointer_from_objref(A)::Ptr{$type}, (I[1]-1)::Cint, (I[2]-1)::Cint)::eltype($type) + + # setindex! + blasfeo_gein1 = Symbol(:blasfeo_, flag, :gein1) + @eval Base.setindex!(A::$type, v::T, I::Vararg{Int, 2}) where {T <: Real} = @ccall blasfeo.$blasfeo_gein1(v::eltype($type),pointer_from_objref(A)::Ptr{$type}, (I[1]-1)::Cint, (I[2]-1)::Cint)::eltype($type) + + # similar @eval Base.similar(A::$type) = $type(A.m, A.n) @eval Base.similar(A::$type,dims::Dims{2}) = $type(dims...) + # copy + blasfeo_gecp = Symbol(:blasfeo_, flag, :gecp) + @eval function copy(A::$type) + B = similar(A) + A_ptr = pointer_from_objref(A) + B_ptr = pointer_from_objref(B) + @ccall blasfeo.blasfeo_gecp( + A.m::Cint, A.n::Cint, + A_ptr::Ptr{$type}, 0::Cint, 0::Cint, + B_ptr::Ptr{$type}, 0::Cint, 0::Cint + )::Cvoid + return B + end + + # fill + blasfeo_gese = Symbol(:blasfeo_, flag, :gese) + @eval function fill(A::$type, b::T) where {T <: Real} + A_ptr = pointer_from_objref(A) + @ccall blasfeo.blasfeo_gese( + A.m::Cint, A.n::Cint, + b::eltype($type), + A_ptr::Ptr{$type}, 0::Cint, 0::Cint, + )::Cvoid + end end diff --git a/src/vecs.jl b/src/vecs.jl index 4fa1356..7384649 100644 --- a/src/vecs.jl +++ b/src/vecs.jl @@ -11,7 +11,7 @@ mutable struct BlasfeoDvec <: AbstractVector{Cdouble} pm::Cint # packed size memsize::Cint # size of needed memory - function BlasfeoDvec(m::Int) + function BlasfeoDvec(m::T) where {T <: Integer} vec = new(C_NULL,C_NULL,0,0,0) @ccall blasfeo.blasfeo_allocate_dvec(m::Cint, pointer_from_objref(vec)::Ptr{BlasfeoDvec})::Cvoid function destructor(this) @@ -44,7 +44,7 @@ mutable struct BlasfeoSvec <: AbstractVector{Cfloat} pm::Cint # packed size memsize::Cint # size of needed memory - function BlasfeoSvec(m::Int,n::Int) + function BlasfeoSvec(m::T) where {T <: Integer} vec = new(C_NULL,C_NULL,0,0,0) @ccall blasfeo.blasfeo_allocate_svec(m::Cint, pointer_from_objref(vec)::Ptr{BlasfeoSvec})::Cvoid function destructor(this) @@ -68,3 +68,49 @@ mutable struct BlasfeoSvec <: AbstractVector{Cfloat} return finalizer(destructor, vec) end end + +# basic vector operations +for (type,flag) in [ + (:BlasfeoDvec, :d), + (:BlasfeoSvec, :s), + ] + # size + @eval Base.size(A::$type) = (A.m,) + + # getindex + blasfeo_vecex1 = Symbol(:blasfeo_, flag, :vecex1) + @eval Base.getindex(A::$type, I::Vararg{Int, 1}) = @ccall blasfeo.$blasfeo_vecex1(pointer_from_objref(A)::Ptr{$type}, (I[1]-1)::Cint)::eltype($type) + + # setindex! + blasfeo_vecin1 = Symbol(:blasfeo_, flag, :vecin1) + @eval Base.setindex!(A::$type, v::T, I::Vararg{Int, 1}) where {T <: Real} = @ccall blasfeo.$blasfeo_vecin1(v::eltype($type),pointer_from_objref(A)::Ptr{$type}, (I[1]-1)::Cint)::eltype($type) + + # similar + @eval Base.similar(A::$type) = $type(A.m) + @eval Base.similar(A::$type,dims::Dims{1}) = $type(dims...) + + # copy + blasfeo_veccp = Symbol(:blasfeo_, flag, :veccp) + @eval function copy(A::$type) + B = similar(A) + A_ptr = pointer_from_objref(A) + B_ptr = pointer_from_objref(B) + @ccall blasfeo.blasfeo_veccp( + A.m::Cint, + A_ptr::Ptr{$type}, 0::Cint, + B_ptr::Ptr{$type}, 0::Cint, + )::Cvoid + return B + end + + # fill + blasfeo_vecse = Symbol(:blasfeo_, flag, :vecse) + @eval function fill(A::$type, b::T) where {T <: Real} + A_ptr = pointer_from_objref(A) + @ccall blasfeo.blasfeo_vecse( + A.m::Cint, + b::eltype($type), + A_ptr::Ptr{$type}, 0::Cint, + )::Cvoid + end +end diff --git a/test/mat/basics.jl b/test/mat/basics.jl new file mode 100644 index 0000000..d8b1e02 --- /dev/null +++ b/test/mat/basics.jl @@ -0,0 +1,27 @@ +@testset "Matrix Basic Operations" begin + @testset for MAT in (BlasfeoDmat, BlasfeoSmat) + A = rand(eltype(MAT), 100, 100) + A_blasfeo = MAT(A) + @test A == A_blasfeo + @test A[10,20] == A_blasfeo[10,20] + + A[10,20] = 30.0; A_blasfeo[10,20] = 30.0 + A[20,30] = 50; A_blasfeo[20,30] = 50 + @test A == A_blasfeo + + B_blasfeo = similar(A_blasfeo) + @test size(A_blasfeo) == size(B_blasfeo) + # Note(@anton) Blasfeo _always_ clears memory. + @test A_blasfeo != B_blasfeo + + C_blasfeo = copy(A_blasfeo) + @test A_blasfeo == C_blasfeo + + C_blasfeo[15,15] = 30.0 + @test C_blasfeo[15,15] == 30.0 + @test A_blasfeo != C_blasfeo + + fill!(C_blasfeo, 100.0) + @test all(C_blasfeo .== 100.0) + end +end diff --git a/test/runtests.jl b/test/runtests.jl index f3ad109..63df0b2 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,6 +1,5 @@ using BLASFEO using Test -@testset "BLASFEO.jl" begin - # Write your own tests here. -end +include("mat/basics.jl") +include("vec/basics.jl") diff --git a/test/vec/basics.jl b/test/vec/basics.jl new file mode 100644 index 0000000..90084ee --- /dev/null +++ b/test/vec/basics.jl @@ -0,0 +1,27 @@ +@testset "Vector Basic Operations" begin + @testset for VEC in (BlasfeoDvec, BlasfeoSvec) + A = rand(eltype(VEC), 100) + A_blasfeo = VEC(A) + @test A == A_blasfeo + @test A[10] == A_blasfeo[10] + + A[10] = 10.0; A_blasfeo[10] = 10.0 + A[20] = 20; A_blasfeo[20] = 20 + @test A == A_blasfeo + + B_blasfeo = similar(A_blasfeo) + @test size(A_blasfeo) == size(B_blasfeo) + # Note(@anton) Blasfeo _always_ clears memory. + @test A_blasfeo != B_blasfeo + + C_blasfeo = copy(A_blasfeo) + @test A_blasfeo == C_blasfeo + + C_blasfeo[15] = 15.0 + @test C_blasfeo[15] == 15.0 + @test A_blasfeo != C_blasfeo + + fill!(C_blasfeo, 100.0) + @test all(C_blasfeo .== 100.0) + end +end From 49d11ce213d2256e3db4a2f893b46a3f3151331a Mon Sep 17 00:00:00 2001 From: Anton Pozharskiy Date: Fri, 12 Jun 2026 18:18:05 +0200 Subject: [PATCH 4/4] start with lvl1 blas operations --- Project.toml | 2 ++ src/BLASFEO.jl | 5 ++++ src/level1.jl | 75 ++++++++++++++++++++++++++++++++++++++++++++++++++ src/mats.jl | 40 +++++++++++++++++++++++---- src/vecs.jl | 11 ++++---- 5 files changed, 123 insertions(+), 10 deletions(-) create mode 100644 src/level1.jl diff --git a/Project.toml b/Project.toml index 8711a32..01e1db6 100644 --- a/Project.toml +++ b/Project.toml @@ -4,9 +4,11 @@ version = "0.1.0" authors = ["Ian McInerney", "Imperial College London"] [deps] +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" blasfeo_jll = "6b574d4a-bb57-5a4b-b7e6-1c794c903646" [compat] +LinearAlgebra = "1.12.0" julia = "1" [extras] diff --git a/src/BLASFEO.jl b/src/BLASFEO.jl index 3fe7bae..d32c6ed 100644 --- a/src/BLASFEO.jl +++ b/src/BLASFEO.jl @@ -1,12 +1,17 @@ module BLASFEO +using LinearAlgebra using blasfeo_jll include("vecs.jl") include("mats.jl") include("show.jl") +include("level1.jl") # Vectors export BlasfeoDvec, BlasfeoSvec # Matrices export BlasfeoDmat, BlasfeoSmat +# Level 1 +export axpy, axpby + end # module diff --git a/src/level1.jl b/src/level1.jl new file mode 100644 index 0000000..a9fc235 --- /dev/null +++ b/src/level1.jl @@ -0,0 +1,75 @@ +for (type,flag) in [ + (:BlasfeoDvec, :d), + (:BlasfeoSvec, :s), + ] + # dot + blasfeo_dot = Symbol(:blasfeo_, flag, :dot) + @eval function LinearAlgebra.dot(a::$type, b::$type) + a_ptr=pointer_from_objref(a) + b_ptr=pointer_from_objref(b) + @ccall blasfeo.$blasfeo_dot( + a.m::Cint, a_ptr::Ptr{$type}, 0::Cint, + b_ptr::Ptr{$type}, 0::Cint, + )::eltype($type) + end + + # axpy! + blasfeo_axpy = Symbol(:blasfeo_, flag, :axpy) + @eval function LinearAlgebra.axpy!(α::T, x::$type, y::$type) where {T <: Real} + x_ptr=pointer_from_objref(x) + y_ptr=pointer_from_objref(y) + @ccall blasfeo.$blasfeo_axpy( + x.m::Cint, + α::Cdouble, x_ptr::Ptr{$type}, 0::Cint, + y_ptr::Ptr{$type}, 0::Cint, + y_ptr::Ptr{$type}, 0::Cint, + )::Cvoid + return y + end + + # 3 arg axpby + @eval function axpy(α::T, x::$type, y::$type, z::$type) where {T <: Real} + x_ptr=pointer_from_objref(x) + y_ptr=pointer_from_objref(y) + z_ptr=pointer_from_objref(z) + @ccall blasfeo.$blasfeo_axpy( + x.m::Cint, + α::Cdouble, x_ptr::Ptr{$type}, 0::Cint, + y_ptr::Ptr{$type}, 0::Cint, + z_ptr::Ptr{$type}, 0::Cint, + )::Cvoid + return z + end + + # axpby + blasfeo_axpby = Symbol(:blasfeo_, flag, :axpby) + @eval function LinearAlgebra.axpby!(α::T1, x::$type, β::T2, y::$type) where {T1 <: Real, T2 <: Real} + x_ptr=pointer_from_objref(x) + y_ptr=pointer_from_objref(y) + @ccall blasfeo.$blasfeo_axpby( + a.m::Cint, + α::Cdouble, x_ptr::Ptr{$type}, 0::Cint, + β::Cdouble, y_ptr::Ptr{$type}, 0::Cint, + y_ptr::Ptr{$type}, 0::Cint, + )::Cvoid + return y + end + + # 3 arg axpby + blasfeo_axpby = Symbol(:blasfeo_, flag, :axpby) + @eval function axpby(α::T1, x::$type, β::T2, y::$type, z::$type) where {T1 <: Real, T2 <: Real} + x_ptr=pointer_from_objref(x) + y_ptr=pointer_from_objref(y) + z_ptr=pointer_from_objref(z) + @ccall blasfeo.$blasfeo_axpby( + x.m::Cint, + α::Cdouble, x_ptr::Ptr{$type}, 0::Cint, + β::Cdouble, y_ptr::Ptr{$type}, 0::Cint, + z_ptr::Ptr{$type}, 0::Cint, + )::Cvoid + return z + end + + # :* + #@eval function Base.:* +end diff --git a/src/mats.jl b/src/mats.jl index 36c6fb5..126db22 100644 --- a/src/mats.jl +++ b/src/mats.jl @@ -99,11 +99,11 @@ for (type,flag) in [ # copy blasfeo_gecp = Symbol(:blasfeo_, flag, :gecp) - @eval function copy(A::$type) + @eval function Base.copy(A::$type) B = similar(A) A_ptr = pointer_from_objref(A) B_ptr = pointer_from_objref(B) - @ccall blasfeo.blasfeo_gecp( + @ccall blasfeo.$blasfeo_gecp( A.m::Cint, A.n::Cint, A_ptr::Ptr{$type}, 0::Cint, 0::Cint, B_ptr::Ptr{$type}, 0::Cint, 0::Cint @@ -111,14 +111,44 @@ for (type,flag) in [ return B end - # fill + # fill! blasfeo_gese = Symbol(:blasfeo_, flag, :gese) - @eval function fill(A::$type, b::T) where {T <: Real} + @eval function Base.fill!(A::$type, b::T) where {T <: Real} A_ptr = pointer_from_objref(A) - @ccall blasfeo.blasfeo_gese( + @ccall blasfeo.$blasfeo_gese( A.m::Cint, A.n::Cint, b::eltype($type), A_ptr::Ptr{$type}, 0::Cint, 0::Cint, )::Cvoid + return A + end + + # convert to Matrix and back + # TODO(@anton) maybe the performance overhead of this means maybe we want only explicit constructors? + blasfeo_unpack = Symbol(:blasfeo_unpack_, flag, :mat) + blasfeo_pack = Symbol(:blasfeo_pack_, flag, :mat) + @eval function Base.convert(::Type{Matrix{eltype($type)}}, A::$type) + A_ptr = pointer_from_objref(A) + B = Matrix{eltype($type)}(undef, size(A)) + @ccall blasfeo.$blasfeo_unpack( + A.m::Cint, A.n::Cint, + A_ptr::Ptr{$type}, + 0::Cint, 0::Cint, + B::Ptr{eltype($type)}, + A.m::Cint + )::Cvoid + return B + end + @eval function Base.convert(::Type{$type}, A::Matrix{eltype($type)}) + B = $type(size(A)...) + B_ptr = pointer_from_objref(B) + m,n = size(A) + @ccall blasfeo.$blasfeo_pack( + m::Cint, n::Cint, + A::Ptr{eltype($type)}, m::Cint, + B_ptr::Ptr{$type}, + 0::Cint, 0::Cint + )::Cvoid + return B end end diff --git a/src/vecs.jl b/src/vecs.jl index 7384649..509809e 100644 --- a/src/vecs.jl +++ b/src/vecs.jl @@ -91,11 +91,11 @@ for (type,flag) in [ # copy blasfeo_veccp = Symbol(:blasfeo_, flag, :veccp) - @eval function copy(A::$type) + @eval function Base.copy(A::$type) B = similar(A) A_ptr = pointer_from_objref(A) B_ptr = pointer_from_objref(B) - @ccall blasfeo.blasfeo_veccp( + @ccall blasfeo.$blasfeo_veccp( A.m::Cint, A_ptr::Ptr{$type}, 0::Cint, B_ptr::Ptr{$type}, 0::Cint, @@ -103,14 +103,15 @@ for (type,flag) in [ return B end - # fill + # fill! blasfeo_vecse = Symbol(:blasfeo_, flag, :vecse) - @eval function fill(A::$type, b::T) where {T <: Real} + @eval function Base.fill!(A::$type, b::T) where {T <: Real} A_ptr = pointer_from_objref(A) - @ccall blasfeo.blasfeo_vecse( + @ccall blasfeo.$blasfeo_vecse( A.m::Cint, b::eltype($type), A_ptr::Ptr{$type}, 0::Cint, )::Cvoid + return A end end