From 8b4c1bcad1cd4eddec4564a503da1cd3ca11b7a6 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Tue, 28 Apr 2026 06:42:23 -0500 Subject: [PATCH 1/3] Add .git-blame-ignore-revs for runic formatting commit --- .git-blame-ignore-revs | 2 ++ 1 file changed, 2 insertions(+) create mode 100644 .git-blame-ignore-revs diff --git a/.git-blame-ignore-revs b/.git-blame-ignore-revs new file mode 100644 index 0000000..4e35f23 --- /dev/null +++ b/.git-blame-ignore-revs @@ -0,0 +1,2 @@ +# runic formatting +43e2848 From f8a0c7bc017500c417e7deb557a7be7a6eb26e68 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Tue, 28 Apr 2026 06:58:11 -0500 Subject: [PATCH 2/3] Correctness: `any` and `issymmetric` `any` had been implemented for two functions, `isnan` and `isinf`, because those has an efficient implementation in terms of the factors. More broadly, though, it's inconsistent to support `any` for some functions but not others. This adds a fallback that first computes the dense matrix. `issymmetric` had been implemented as a simple check of the factors, but that checks only *structural* symmetry, not *numerical* symmetry. This deletes the method. --- src/FactoredMatrices.jl | 9 ++++----- test/runtests.jl | 3 +++ 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/src/FactoredMatrices.jl b/src/FactoredMatrices.jl index cd264d7..dd5d43e 100644 --- a/src/FactoredMatrices.jl +++ b/src/FactoredMatrices.jl @@ -89,6 +89,9 @@ function Base.show(io::IO, A::FactoredMatrix) return print(io, "FactoredMatrix of size ", size(A), " with rank ", size(A.U, 2)) end +Base.any(f::Union{typeof(isinf), typeof(isnan)}, A::FactoredMatrix) = any(f, A.U) || any(f, A.V) +Base.any(f, A::FactoredMatrix) = any(f, Matrix(A)) + 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)) @@ -222,9 +225,7 @@ end *(A::FactoredMatrix, B) = A.U * (A.V * B) *(A, B::FactoredMatrix) = (A * B.U) * B.V -# --- Queries --- - -LinearAlgebra.issymmetric(A::FactoredMatrix) = A.U == A.V' +# --- misc --- function dot(A::FactoredMatrix, B::FactoredMatrix) M1 = B.U' * A.U @@ -239,6 +240,4 @@ Compute the sum of squared differences between `A` and `B` without forming the f """ 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) - end # module diff --git a/test/runtests.jl b/test/runtests.jl index b45e868..f13eb10 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -64,6 +64,9 @@ end @test E * M' ≈ E * Mf' @test !any(isnan, Mf) @test !any(isinf, Mf) + V1, U1 = [1; 1;;], [1 -1] + M1 = FactoredMatrix(U1, V1) + @test any(iszero, M1) # Workspace (cache=) path: correctness + zero allocations ws = FactoredMatrices.Workspace(Mf, 5) From a32b54bfbb7aa5591613fb6a8bec8e7bc3ac192f Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Tue, 28 Apr 2026 07:18:33 -0500 Subject: [PATCH 3/3] Add 5-arg mul methods The downside is that it adds a lot of ambiguities against banded and `AbstractTriangular` methods. --- src/FactoredMatrices.jl | 48 ++++++++++++++++++++++++++++++++++++---- test/runtests.jl | 49 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 93 insertions(+), 4 deletions(-) diff --git a/src/FactoredMatrices.jl b/src/FactoredMatrices.jl index dd5d43e..c6f6ab8 100644 --- a/src/FactoredMatrices.jl +++ b/src/FactoredMatrices.jl @@ -110,14 +110,14 @@ const _multypes = (FactoredMatrix, Adjoint{T, <:FactoredMatrix{T}} where {T}, Tr 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::AbstractMatrix, 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::AbstractMatrix, A::$T, B::AbstractMatrix; cache::Union{Nothing, Workspace} = nothing) = _mul!(C, rewrap(A), B, cache) + mul!(C::AbstractVecOrMat, A::$T, b::AbstractVector; cache::Union{Nothing, Workspace} = nothing) = _mul!(C, rewrap(A), b, cache) + mul!(C::AbstractMatrix, A::AbstractMatrix, B::$T; cache::Union{Nothing, Workspace} = nothing) = _mul!(C, A, rewrap(B), cache) end end @@ -139,6 +139,24 @@ for CT in _multypes end end +# 5-arg mul!: C = α*A*B + β*C +for AT in _multypes, BT in _multypes + @eval begin + mul!(C::AbstractMatrix, A::$AT, B::$BT, α::Number, β::Number; cache::Union{Nothing, Workspace} = nothing) = + _mul!(C, rewrap(A), rewrap(B), α, β, cache) + end +end +for T in _multypes + @eval begin + mul!(C::AbstractMatrix, A::$T, B::AbstractMatrix, α::Number, β::Number; cache::Union{Nothing, Workspace} = nothing) = + _mul!(C, rewrap(A), B, α, β, cache) + mul!(C::AbstractVector, A::$T, b::AbstractVector, α::Number, β::Number; cache::Union{Nothing, Workspace} = nothing) = + _mul!(C, rewrap(A), b, α, β, cache) + mul!(C::AbstractMatrix, A::AbstractMatrix, B::$T, α::Number, β::Number; cache::Union{Nothing, Workspace} = nothing) = + _mul!(C, A, rewrap(B), α, β, cache) + end +end + # Internals # When the output is an AbstractMatrix: @@ -179,6 +197,28 @@ function _mul!(C, A::AbstractMatrix, B::FactoredMatrix, cache::Union{Nothing, Wo return C end +# 5-arg _mul!: C = α*A*B + β*C +function _mul!(C, A::FactoredMatrix, B::FactoredMatrix, α::Number, β::Number, cache::Union{Nothing, Workspace}) + tmp = cache === nothing ? A.V * B.U : mul!(cache.templeft, A.V, B.U) + k, p = size(tmp) + 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, α::Number, β::Number, cache::Union{Nothing, Workspace}) + tmp = cache === nothing ? A.V * B : mul!(cache.templeft, A.V, B) + return mul!(C, A.U, tmp, α, β) +end + +function _mul!(C, A::FactoredMatrix, b::AbstractVector, α::Number, β::Number, cache::Union{Nothing, Workspace}) + tmp = cache === nothing ? A.V * b : mul!(view(cache.templeft, :, 1), A.V, b) + return mul!(C, A.U, tmp, α, β) +end + +function _mul!(C, A::AbstractMatrix, B::FactoredMatrix, α::Number, β::Number, cache::Union{Nothing, Workspace}) + tmp = cache === nothing ? A * B.U : mul!(cache.tempright, A, B.U) + return mul!(C, tmp, B.V, α, β) +end + # When the output is a FactoredMatrix: function _mul!(C::FactoredMatrix, A::FactoredMatrix, B::FactoredMatrix, cache::Union{Nothing, Workspace}) tmp = if cache === nothing diff --git a/test/runtests.jl b/test/runtests.jl index f13eb10..dc3b88a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -198,5 +198,54 @@ end mul!(c, Mf, b; cache = ws_vec); @test c ≈ M * b check_vec_alloc(Mf, b, ws_vec, c) end + + # 5-arg mul!: C = α*A*B + β*C + # Mf is 15×10, M is its materialization + # C is 10×5 (right for Mf), D is 15×5 (right for Mf'), E is 5×10 (left for Mf') + let + α = T(2) + β = T(3) + R15x5 = myrand(T, 15, 5) # result shape for FM×AM, AM×adj(FM) is 5×15 + R10x5 = myrand(T, 10, 5) # result shape for adj(FM)×AM + R5x10 = myrand(T, 5, 10) # result shape for AM×FM + R5x15 = myrand(T, 5, 15) # result shape for E×adj(Mf) + E_pre = myrand(T, 5, 15) # 5×15 matrix: E_pre * Mf (15×10) → 5×10 + b5 = vec(myrand(T, 10, 1)) + c5 = vec(myrand(T, 15, 1)) + + # FM(15×10) × AM(10×5) → 15×5 + R = copy(R15x5); mul!(R, Mf, C, α, β); @test R ≈ α * M * C + β * R15x5 + # adj(FM)(10×15) × AM(15×5) → 10×5 + R = copy(R10x5); mul!(R, Adjoint(Mf), D, α, β); @test R ≈ α * M' * D + β * R10x5 + # transpose(FM)(10×15) × AM(15×5) → 10×5 + R = copy(R10x5); mul!(R, Transpose(Mf), D, α, β); @test R ≈ α * transpose(M) * D + β * R10x5 + # AM(5×15) × FM(15×10) → 5×10 + R = copy(R5x10); mul!(R, E_pre, Mf, α, β); @test R ≈ α * E_pre * M + β * R5x10 + # AM(5×10) × adj(FM)(10×15) → 5×15 + R = copy(R5x15); mul!(R, E, Adjoint(Mf), α, β); @test R ≈ α * E * M' + β * R5x15 + # AM(5×10) × transpose(FM)(10×15) → 5×15 + R = copy(R5x15); mul!(R, E, Transpose(Mf), α, β); @test R ≈ α * E * transpose(M) + β * R5x15 + # FM(15×10) × vector(10) → 15 + R = copy(c5); mul!(R, Mf, b5, α, β); @test R ≈ α * M * b5 + β * c5 + # adj(FM)(10×15) × vector(15) → 10 + b15 = vec(myrand(T, 15, 1)) + c10 = vec(myrand(T, 10, 1)) + R = copy(c10); mul!(R, Adjoint(Mf), b15, α, β); @test R ≈ α * M' * b15 + β * c10 + + # FM × FM + Mf_A = FactoredMatrix(myrand(T, 10, 3), myrand(T, 3, 8)); MA = Array(Mf_A) + Mf_B = FactoredMatrix(myrand(T, 15, 3), myrand(T, 3, 8)); MB = Array(Mf_B) + R15x8 = myrand(T, 15, 8) + R10x8 = myrand(T, 10, 8) + R = copy(R15x8); mul!(R, Mf, Mf_A, α, β); @test R ≈ α * M * MA + β * R15x8 + R = copy(R10x8); mul!(R, Adjoint(Mf), Mf_B, α, β); @test R ≈ α * M' * MB + β * R10x8 + R = copy(R10x8); mul!(R, Transpose(Mf), Mf_B, α, β); @test R ≈ α * transpose(M) * MB + β * R10x8 + + # with cache + ws5 = FactoredMatrices.Workspace(Mf, 5) + R = copy(R15x5); mul!(R, Mf, C, α, β; cache = ws5); @test R ≈ α * M * C + β * R15x5 + R = copy(R10x5); mul!(R, Adjoint(Mf), D, α, β; cache = ws5); @test R ≈ α * M' * D + β * R10x5 + R = copy(c5); mul!(R, Mf, b5, α, β; cache = ws5); @test R ≈ α * M * b5 + β * c5 + end end end