Skip to content

Apply runic formatting #12

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We鈥檒l occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
timholy merged 1 commit into from
Apr 28, 2026
Merged

Apply runic formatting #12

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
71 changes: 36 additions & 35 deletions src/FactoredMatrices.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,26 +16,27 @@ Multiplications exploit the factored structure without forming the full dense ma
For allocation-free repeated multiplications, create a [`Workspace`](@ref) and pass it
as `mul!(C, A, B; cache=ws)`.
"""
struct FactoredMatrix{T, UT<:AbstractMatrix{T}, VT<:AbstractMatrix{T}} <: Factorization{T}
struct FactoredMatrix{T, UT <: AbstractMatrix{T}, VT <: AbstractMatrix{T}} <: Factorization{T}
U::UT # m 脳 k
V::VT # k 脳 n

function FactoredMatrix{T, UT, VT}(U::UT, V::VT) where {T, UT, VT}
# This deliberately requires an exact match to ensure these buffers will be used
size(U, 2) == size(V, 1) ||
throw(DimensionMismatch("U is $(size(U,1))脳$(size(U,2)) but V is $(size(V,1))脳$(size(V,2))"))
new{T, UT, VT}(U, V)
throw(DimensionMismatch("U is $(size(U, 1))脳$(size(U, 2)) but V is $(size(V, 1))脳$(size(V, 2))"))
return new{T, UT, VT}(U, V)
end
end
FactoredMatrix{T}(U::AbstractMatrix{T}, V::AbstractMatrix{T}) where T = FactoredMatrix{T, typeof(U), typeof(V)}(U, V)
FactoredMatrix(U::AbstractMatrix{T}, V::AbstractMatrix{T}) where T = FactoredMatrix{T}(U, V)
FactoredMatrix{T}(U::AbstractMatrix{T}, V::AbstractMatrix{T}) where {T} = FactoredMatrix{T, typeof(U), typeof(V)}(U, V)
FactoredMatrix(U::AbstractMatrix{T}, V::AbstractMatrix{T}) where {T} = FactoredMatrix{T}(U, V)

function FactoredMatrix(U::Matrix{T}, V::Matrix{T}, ::Integer) where T
function FactoredMatrix(U::Matrix{T}, V::Matrix{T}, ::Integer) where {T}
Base.depwarn(
"FactoredMatrix(U, V, j) is deprecated; use FactoredMatrix(U, V) and " *
"create a Workspace(fm, j) for allocation-free mul!",
:FactoredMatrix)
FactoredMatrix{T}(U, V)
"create a Workspace(fm, j) for allocation-free mul!",
:FactoredMatrix
)
return FactoredMatrix{T}(U, V)
end

"""
Expand Down Expand Up @@ -67,9 +68,9 @@ struct Workspace{T}
tempright::Matrix{T} # p 脳 k: intermediate for A 脳 (right FM or its adjoint/transpose)
end

Workspace{T}(k::Int, p::Int) where T = Workspace{T}(Matrix{T}(undef, k, p), Matrix{T}(undef, p, k))
Workspace{T}(k::Int, p::Int) where {T} = Workspace{T}(Matrix{T}(undef, k, p), Matrix{T}(undef, p, k))

Workspace(A::FactoredMatrix{T}, p::Int) where T = Workspace{T}(size(A.U, 2), p)
Workspace(A::FactoredMatrix{T}, p::Int) where {T} = Workspace{T}(size(A.U, 2), p)

size(A::FactoredMatrix) = size(A.U, 1), size(A.V, 2)
size(A::FactoredMatrix, d::Integer) = d == 1 ? size(A.U, 1) : (d == 2 ? size(A.V, 2) : 1)
Expand All @@ -85,11 +86,11 @@ const FMhash = Int === Int64 ? 0x99ac4c2c56e5bd6e : 0x70c3ac8c
hash(A::FactoredMatrix, h::UInt) = hash(A.U, hash(A.V, hash(FMhash, h)))

function Base.show(io::IO, A::FactoredMatrix)
print(io, "FactoredMatrix of size ", size(A), " with rank ", size(A.U, 2))
return print(io, "FactoredMatrix of size ", size(A), " with rank ", size(A.U, 2))
end

adjoint(A::FactoredMatrix{T}) where T = FactoredMatrix{T}(adjoint(A.V), adjoint(A.U))
transpose(A::FactoredMatrix{T}) where T = FactoredMatrix{T}(transpose(A.V), transpose(A.U))
adjoint(A::FactoredMatrix{T}) where {T} = FactoredMatrix{T}(adjoint(A.V), adjoint(A.U))
transpose(A::FactoredMatrix{T}) where {T} = FactoredMatrix{T}(transpose(A.V), transpose(A.U))

# --- mul! ---

Expand All @@ -102,53 +103,53 @@ rewrap(A::FactoredMatrix) = A
rewrap(A::Transpose{T, <:FactoredMatrix{T}}) where {T} = FactoredMatrix{T}(transpose(parent(A).V), transpose(parent(A).U))
rewrap(A::Adjoint{T, <:FactoredMatrix{T}}) where {T} = FactoredMatrix{T}(adjoint(parent(A).V), adjoint(parent(A).U))

const _multypes = (FactoredMatrix, Adjoint{T, <:FactoredMatrix{T}} where T, Transpose{T, <:FactoredMatrix{T}} where T)
const _multypes = (FactoredMatrix, Adjoint{T, <:FactoredMatrix{T}} where {T}, Transpose{T, <:FactoredMatrix{T}} where {T})

for AT in _multypes, BT in _multypes
@eval begin
mul!(C, A::$AT, B::$BT; cache::Union{Nothing,Workspace}=nothing) = _mul!(C, rewrap(A), rewrap(B), cache)
mul!(C, A::$AT, B::$BT; cache::Union{Nothing, Workspace} = nothing) = _mul!(C, rewrap(A), rewrap(B), cache)
end
end
for T in _multypes
@eval begin
mul!(C, A::$T, B::AbstractMatrix; cache::Union{Nothing,Workspace}=nothing) = _mul!(C, rewrap(A), B, cache)
mul!(C, A::$T, b::AbstractVector; cache::Union{Nothing,Workspace}=nothing) = _mul!(C, rewrap(A), b, cache)
mul!(C, A::AbstractMatrix, B::$T; cache::Union{Nothing,Workspace}=nothing) = _mul!(C, A, rewrap(B), cache)
mul!(C, A::$T, B::AbstractMatrix; cache::Union{Nothing, Workspace} = nothing) = _mul!(C, rewrap(A), B, cache)
mul!(C, A::$T, b::AbstractVector; cache::Union{Nothing, Workspace} = nothing) = _mul!(C, rewrap(A), b, cache)
mul!(C, A::AbstractMatrix, B::$T; cache::Union{Nothing, Workspace} = nothing) = _mul!(C, A, rewrap(B), cache)
end
end

# Now with a FactoredMatrix output
for CT in _multypes, AT in _multypes, BT in _multypes
@eval begin
mul!(C::$CT, A::$AT, B::$BT; cache::Union{Nothing,Workspace}=nothing) = _mul!(C, rewrap(A), rewrap(B), cache)
mul!(C::$CT, A::$AT, B::$BT; cache::Union{Nothing, Workspace} = nothing) = _mul!(C, rewrap(A), rewrap(B), cache)
end
end
for CT in _multypes, T in _multypes
@eval begin
mul!(C::$CT, A::$T, B::AbstractMatrix; cache::Union{Nothing,Workspace}=nothing) = _mul!(C, rewrap(A), B, nothing)
mul!(C::$CT, A::AbstractMatrix, B::$T; cache::Union{Nothing,Workspace}=nothing) = _mul!(C, A, rewrap(B), nothing)
mul!(C::$CT, A::$T, B::AbstractMatrix; cache::Union{Nothing, Workspace} = nothing) = _mul!(C, rewrap(A), B, nothing)
mul!(C::$CT, A::AbstractMatrix, B::$T; cache::Union{Nothing, Workspace} = nothing) = _mul!(C, A, rewrap(B), nothing)
end
end
for CT in _multypes
@eval begin
mul!(C::$CT, A::AbstractMatrix, B::AbstractMatrix; cache::Union{Nothing,Workspace}=nothing) = _mul!(C, A, B, nothing)
mul!(C::$CT, A::AbstractMatrix, B::AbstractMatrix; cache::Union{Nothing, Workspace} = nothing) = _mul!(C, A, B, nothing)
end
end

# Internals

# When the output is an AbstractMatrix:
function _mul!(C, A::FactoredMatrix, B::FactoredMatrix, cache::Union{Nothing,Workspace})
function _mul!(C, A::FactoredMatrix, B::FactoredMatrix, cache::Union{Nothing, Workspace})
tmp = if cache === nothing
A.V * B.U
else
mul!(cache.templeft, A.V, B.U)
end
k, p = size(tmp)
return k <= p ? mul!(C, A.U * tmp, B.V) : mul!(C, A.U, tmp * B.V)
return k <= p ? mul!(C, A.U * tmp, B.V) : mul!(C, A.U, tmp * B.V)
end

function _mul!(C, A::FactoredMatrix, B::AbstractMatrix, cache::Union{Nothing,Workspace})
function _mul!(C, A::FactoredMatrix, B::AbstractMatrix, cache::Union{Nothing, Workspace})
if cache === nothing
mul!(C, A.U, A.V * B)
else
Expand All @@ -157,7 +158,7 @@ function _mul!(C, A::FactoredMatrix, B::AbstractMatrix, cache::Union{Nothing,Wor
return C
end

function _mul!(C, A::FactoredMatrix, b::AbstractVector, cache::Union{Nothing,Workspace})
function _mul!(C, A::FactoredMatrix, b::AbstractVector, cache::Union{Nothing, Workspace})
if cache === nothing
mul!(C, A.U, A.V * b)
else
Expand All @@ -166,7 +167,7 @@ function _mul!(C, A::FactoredMatrix, b::AbstractVector, cache::Union{Nothing,Wor
return C
end

function _mul!(C, A::AbstractMatrix, B::FactoredMatrix, cache::Union{Nothing,Workspace})
function _mul!(C, A::AbstractMatrix, B::FactoredMatrix, cache::Union{Nothing, Workspace})
if cache === nothing
mul!(C, A * B.U, B.V)
else
Expand All @@ -176,7 +177,7 @@ function _mul!(C, A::AbstractMatrix, B::FactoredMatrix, cache::Union{Nothing,Wor
end

# When the output is a FactoredMatrix:
function _mul!(C::FactoredMatrix, A::FactoredMatrix, B::FactoredMatrix, cache::Union{Nothing,Workspace})
function _mul!(C::FactoredMatrix, A::FactoredMatrix, B::FactoredMatrix, cache::Union{Nothing, Workspace})
tmp = if cache === nothing
A.V * B.U
else
Expand All @@ -195,21 +196,21 @@ function _mul!(C::FactoredMatrix, A::FactoredMatrix, B::FactoredMatrix, cache::U
return C
end

function _mul!(C::FactoredMatrix, A::FactoredMatrix, B::AbstractMatrix, ::Union{Nothing,Workspace})
function _mul!(C::FactoredMatrix, A::FactoredMatrix, B::AbstractMatrix, ::Union{Nothing, Workspace})
# cache is unused
copyto!(C.U, A.U)
mul!(C.V, A.V, B)
return C
end

function _mul!(C::FactoredMatrix, A::AbstractMatrix, B::FactoredMatrix, ::Union{Nothing,Workspace})
function _mul!(C::FactoredMatrix, A::AbstractMatrix, B::FactoredMatrix, ::Union{Nothing, Workspace})
# cache is unused
mul!(C.U, A, B.U)
copyto!(C.V, B.V)
return C
end

function _mul!(C::FactoredMatrix, A::AbstractMatrix, B::AbstractMatrix, ::Union{Nothing,Workspace})
function _mul!(C::FactoredMatrix, A::AbstractMatrix, B::AbstractMatrix, ::Union{Nothing, Workspace})
copyto!(C.U, A)
copyto!(C.V, B)
return C
Expand All @@ -228,16 +229,16 @@ LinearAlgebra.issymmetric(A::FactoredMatrix) = A.U == A.V'
function dot(A::FactoredMatrix, B::FactoredMatrix)
M1 = B.U' * A.U
M2 = B.V * A.V'
sum(M1 .* conj(M2))
return sum(M1 .* conj(M2))
end

"""
ssd(A::FactoredMatrix, B::FactoredMatrix)

Compute the sum of squared differences between `A` and `B` without forming the full matrices.
"""
ssd(A::FactoredMatrix, B::FactoredMatrix) = dot(A, A) - 2*real(dot(A, B)) + dot(B, B)
ssd(A::FactoredMatrix, B::FactoredMatrix) = dot(A, A) - 2 * real(dot(A, B)) + dot(B, B)

Base.any(f::Union{typeof(isinf),typeof(isnan)}, A::FactoredMatrix) = any(f, A.U) || any(f, A.V)
Base.any(f::Union{typeof(isinf), typeof(isnan)}, A::FactoredMatrix) = any(f, A.U) || any(f, A.V)

end # module
Loading
Loading