From ed08bb30351e390ab2d0a339055b697e7332388d Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Tue, 28 Apr 2026 06:37:34 -0500 Subject: [PATCH] Apply runic formatting --- src/FactoredMatrices.jl | 71 +++++++++++----------- test/runtests.jl | 130 ++++++++++++++++++++-------------------- 2 files changed, 101 insertions(+), 100 deletions(-) diff --git a/src/FactoredMatrices.jl b/src/FactoredMatrices.jl index fd2b49d..cd264d7 100644 --- a/src/FactoredMatrices.jl +++ b/src/FactoredMatrices.jl @@ -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 """ @@ -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) @@ -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! --- @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -228,7 +229,7 @@ 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 """ @@ -236,8 +237,8 @@ end 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 diff --git a/test/runtests.jl b/test/runtests.jl index fe7a579..b45e868 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,29 +1,29 @@ using Test, LinearAlgebra using FactoredMatrices -myrand(::Type{T}, m, n) where {T<:Real} = rand(T, m, n) -myrand(::Type{Complex{T}}, m, n) where {T<:Real} = rand(T, m, n) + im*rand(T, m, n) +myrand(::Type{T}, m, n) where {T <: Real} = rand(T, m, n) +myrand(::Type{Complex{T}}, m, n) where {T <: Real} = rand(T, m, n) + im * rand(T, m, n) # Wrapper so the compiler sees ws as a concrete positional argument and can elide # the keyword NamedTuple that would otherwise allocate in Julia 1.10's kwarg machinery. function check_fmfm_fm_alloc(Mf, Mf_A, ws, Cf15x8) - @test @allocated(mul!(Cf15x8, Mf, Mf_A; cache=ws)) == 0 + return @test @allocated(mul!(Cf15x8, Mf, Mf_A; cache = ws)) == 0 end function check_vec_alloc(Mf, b, ws, c) - @test @allocated(mul!(c, Mf, b; cache=ws)) == 0 + return @test @allocated(mul!(c, Mf, b; cache = ws)) == 0 end function check_allocs(Mf, ws, C, D, E, r15x5, r10x5, r5x10, r5x15) - @test @allocated(mul!(r15x5, Mf, C; cache=ws)) == 0 - @test @allocated(mul!(r10x5, Adjoint(Mf), D; cache=ws)) == 0 - @test @allocated(mul!(r10x5, Transpose(Mf), D; cache=ws)) == 0 - @test @allocated(mul!(r5x10, D', Mf; cache=ws)) == 0 - @test @allocated(mul!(r5x10, transpose(D), Mf; cache=ws)) == 0 - @test @allocated(mul!(r15x5, Mf, E'; cache=ws)) == 0 - @test @allocated(mul!(r15x5, Mf, transpose(E); cache=ws)) == 0 - @test @allocated(mul!(r5x15, E, Adjoint(Mf); cache=ws)) == 0 - @test @allocated(mul!(r5x15, E, Transpose(Mf); cache=ws)) == 0 + @test @allocated(mul!(r15x5, Mf, C; cache = ws)) == 0 + @test @allocated(mul!(r10x5, Adjoint(Mf), D; cache = ws)) == 0 + @test @allocated(mul!(r10x5, Transpose(Mf), D; cache = ws)) == 0 + @test @allocated(mul!(r5x10, D', Mf; cache = ws)) == 0 + @test @allocated(mul!(r5x10, transpose(D), Mf; cache = ws)) == 0 + @test @allocated(mul!(r15x5, Mf, E'; cache = ws)) == 0 + @test @allocated(mul!(r15x5, Mf, transpose(E); cache = ws)) == 0 + @test @allocated(mul!(r5x15, E, Adjoint(Mf); cache = ws)) == 0 + return @test @allocated(mul!(r5x15, E, Transpose(Mf); cache = ws)) == 0 end @testset "FactoredMatrices" begin @@ -33,7 +33,7 @@ end V = myrand(T, 3, 10) Mf = FactoredMatrix(U, V) M = Matrix(Mf) - @test M == U*V + @test M == U * V @test M == Array(Mf) Mf2 = FactoredMatrix(copy(U), copy(V)) @test Mf == Mf2 @@ -49,46 +49,46 @@ end C = myrand(T, 10, 5) # for Mf*C (15×10)*(10×5) → 15×5 D = myrand(T, 15, 5) # for Mf'*D (10×15)*(15×5) → 10×5 E = myrand(T, 5, 10) # for E*Mf' (5×10)*(10×15) → 5×15 - @test M*C ≈ Mf*C - @test M'*D ≈ Mf'*D ≈ Adjoint(Mf)*D - @test transpose(M)*D ≈ transpose(Mf)*D ≈ Transpose(Mf)*D - @test D'*M ≈ D'*Mf - @test transpose(D)*M ≈ transpose(D)*Mf - @test M*adjoint(E) ≈ Mf*E' - @test M*transpose(E) ≈ Mf*transpose(E) - @test E*adjoint(M) ≈ E*Mf' - @test E*transpose(M) ≈ E*transpose(Mf) ≈ E*Transpose(Mf) - @test M'*D ≈ Mf'*D - @test D'*M ≈ D'*Mf - @test M*E' ≈ Mf*E' - @test E*M' ≈ E*Mf' + @test M * C ≈ Mf * C + @test M' * D ≈ Mf' * D ≈ Adjoint(Mf) * D + @test transpose(M) * D ≈ transpose(Mf) * D ≈ Transpose(Mf) * D + @test D' * M ≈ D' * Mf + @test transpose(D) * M ≈ transpose(D) * Mf + @test M * adjoint(E) ≈ Mf * E' + @test M * transpose(E) ≈ Mf * transpose(E) + @test E * adjoint(M) ≈ E * Mf' + @test E * transpose(M) ≈ E * transpose(Mf) ≈ E * Transpose(Mf) + @test M' * D ≈ Mf' * D + @test D' * M ≈ D' * Mf + @test M * E' ≈ Mf * E' + @test E * M' ≈ E * Mf' @test !any(isnan, Mf) @test !any(isinf, Mf) # Workspace (cache=) path: correctness + zero allocations ws = FactoredMatrices.Workspace(Mf, 5) - r15x5 = Matrix{T}(undef, 15, 5) # Mf*C, Mf*E', Mf*transpose(E) - r10x5 = Matrix{T}(undef, 10, 5) # Mf'*D, transpose(Mf)*D - r5x10 = Matrix{T}(undef, 5, 10) # D'*Mf, transpose(D)*Mf - r5x15 = Matrix{T}(undef, 5, 15) # E*Mf', E*transpose(Mf) - - mul!(r15x5, Mf, C; cache=ws); @test r15x5 ≈ M*C - mul!(r10x5, Mf', D; cache=ws); @test r10x5 ≈ M'*D - mul!(r10x5, transpose(Mf), D; cache=ws); @test r10x5 ≈ transpose(M)*D - mul!(r5x10, D', Mf; cache=ws); @test r5x10 ≈ D'*M - mul!(r5x10, transpose(D), Mf; cache=ws); @test r5x10 ≈ transpose(D)*M - mul!(r15x5, Mf, E'; cache=ws); @test r15x5 ≈ M*E' - mul!(r15x5, Mf, transpose(E); cache=ws); @test r15x5 ≈ M*transpose(E) - mul!(r5x15, E, Mf'; cache=ws); @test r5x15 ≈ E*M' - mul!(r5x15, E, transpose(Mf); cache=ws); @test r5x15 ≈ E*transpose(M) + r15x5 = Matrix{T}(undef, 15, 5) # Mf*C, Mf*E', Mf*transpose(E) + r10x5 = Matrix{T}(undef, 10, 5) # Mf'*D, transpose(Mf)*D + r5x10 = Matrix{T}(undef, 5, 10) # D'*Mf, transpose(D)*Mf + r5x15 = Matrix{T}(undef, 5, 15) # E*Mf', E*transpose(Mf) + + mul!(r15x5, Mf, C; cache = ws); @test r15x5 ≈ M * C + mul!(r10x5, Mf', D; cache = ws); @test r10x5 ≈ M' * D + mul!(r10x5, transpose(Mf), D; cache = ws); @test r10x5 ≈ transpose(M) * D + mul!(r5x10, D', Mf; cache = ws); @test r5x10 ≈ D' * M + mul!(r5x10, transpose(D), Mf; cache = ws); @test r5x10 ≈ transpose(D) * M + mul!(r15x5, Mf, E'; cache = ws); @test r15x5 ≈ M * E' + mul!(r15x5, Mf, transpose(E); cache = ws); @test r15x5 ≈ M * transpose(E) + mul!(r5x15, E, Mf'; cache = ws); @test r5x15 ≈ E * M' + mul!(r5x15, E, transpose(Mf); cache = ws); @test r5x15 ≈ E * transpose(M) check_allocs(Mf, ws, C, D, E, r15x5, r10x5, r5x10, r5x15) U1 = myrand(T, 15, 4) V1 = myrand(T, 4, 10) - M1 = U1*V1 + M1 = U1 * V1 M1f = FactoredMatrix(U1, V1) - @test sum(abs2, M-M1) ≈ FactoredMatrices.ssd(Mf, M1f) + @test sum(abs2, M - M1) ≈ FactoredMatrices.ssd(Mf, M1f) # FM × FM (all 9 op combinations); Mf is 15×10 # Compatible right-hand FMs: @@ -99,18 +99,18 @@ end let 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) - Mf_C = FactoredMatrix(myrand(T, 8, 3), myrand(T, 3, 10)); MC = Array(Mf_C) - Mf_D = FactoredMatrix(myrand(T, 8, 3), myrand(T, 3, 15)); MD = Array(Mf_D) + Mf_C = FactoredMatrix(myrand(T, 8, 3), myrand(T, 3, 10)); MC = Array(Mf_C) + Mf_D = FactoredMatrix(myrand(T, 8, 3), myrand(T, 3, 15)); MD = Array(Mf_D) r15x8 = Matrix{T}(undef, 15, 8) r10x8 = Matrix{T}(undef, 10, 8) - @test Mf * Mf_A ≈ M * MA - @test Adjoint(Mf) * Mf_B ≈ M' * MB - @test transpose(Mf) * Mf_B ≈ transpose(M) * MB - @test Mf * Adjoint(Mf_C) ≈ M * MC' - @test Mf * Transpose(Mf_C) ≈ M * transpose(MC) - @test Adjoint(Mf) * Mf_D' ≈ M' * MD' - @test Adjoint(Mf) * Transpose(Mf_D) ≈ M' * transpose(MD) - @test Transpose(Mf) * Mf_D' ≈ transpose(M) * MD' + @test Mf * Mf_A ≈ M * MA + @test Adjoint(Mf) * Mf_B ≈ M' * MB + @test transpose(Mf) * Mf_B ≈ transpose(M) * MB + @test Mf * Adjoint(Mf_C) ≈ M * MC' + @test Mf * Transpose(Mf_C) ≈ M * transpose(MC) + @test Adjoint(Mf) * Mf_D' ≈ M' * MD' + @test Adjoint(Mf) * Transpose(Mf_D) ≈ M' * transpose(MD) + @test Transpose(Mf) * Mf_D' ≈ transpose(M) * MD' @test Transpose(Mf) * Transpose(Mf_D) ≈ transpose(M) * transpose(MD) mul!(r15x8, Mf, Mf_A); @test r15x8 ≈ M * MA mul!(r10x8, Adjoint(Mf), Mf_B); @test r10x8 ≈ M' * MB @@ -130,9 +130,9 @@ end r15x8 = Matrix{T}(undef, 15, 8) r10x8 = Matrix{T}(undef, 10, 8) ws3x3 = FactoredMatrices.Workspace{T}(3, 3) - mul!(r15x8, Mf, Mf_A; cache=ws3x3); @test r15x8 ≈ M * MA - mul!(r10x8, Adjoint(Mf), Mf_B; cache=ws3x3); @test r10x8 ≈ M' * MB - mul!(r10x8, Transpose(Mf), Mf_B; cache=ws3x3); @test r10x8 ≈ transpose(M) * MB + mul!(r15x8, Mf, Mf_A; cache = ws3x3); @test r15x8 ≈ M * MA + mul!(r10x8, Adjoint(Mf), Mf_B; cache = ws3x3); @test r10x8 ≈ M' * MB + mul!(r10x8, Transpose(Mf), Mf_B; cache = ws3x3); @test r10x8 ≈ transpose(M) * MB end # mul! with FactoredMatrix output: FM×AM, AM×FM, AM×AM, adj/transpose variants @@ -168,21 +168,21 @@ end let 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) - Mf_C = FactoredMatrix(myrand(T, 8, 3), myrand(T, 3, 10)); MC = Array(Mf_C) + Mf_C = FactoredMatrix(myrand(T, 8, 3), myrand(T, 3, 10)); MC = Array(Mf_C) ws3x3 = FactoredMatrices.Workspace{T}(3, 3) Cf15x8 = FactoredMatrix(Matrix{T}(undef, 15, 3), Matrix{T}(undef, 3, 8)) Cf10x8 = FactoredMatrix(Matrix{T}(undef, 10, 3), Matrix{T}(undef, 3, 8)) - mul!(Cf15x8, Mf, Mf_A); @test Array(Cf15x8) ≈ M * MA - mul!(Cf10x8, Adjoint(Mf), Mf_B); @test Array(Cf10x8) ≈ M' * MB + mul!(Cf15x8, Mf, Mf_A); @test Array(Cf15x8) ≈ M * MA + mul!(Cf10x8, Adjoint(Mf), Mf_B); @test Array(Cf10x8) ≈ M' * MB mul!(Cf10x8, Transpose(Mf), Mf_B); @test Array(Cf10x8) ≈ transpose(M) * MB - mul!(Cf15x8, Mf, Adjoint(Mf_C)); @test Array(Cf15x8) ≈ M * MC' + mul!(Cf15x8, Mf, Adjoint(Mf_C)); @test Array(Cf15x8) ≈ M * MC' mul!(Cf15x8, Mf, Transpose(Mf_C)); @test Array(Cf15x8) ≈ M * transpose(MC) - mul!(Cf15x8, Mf, Mf_A; cache=ws3x3); @test Array(Cf15x8) ≈ M * MA - mul!(Cf10x8, Adjoint(Mf), Mf_B; cache=ws3x3); @test Array(Cf10x8) ≈ M' * MB - mul!(Cf10x8, Transpose(Mf), Mf_B; cache=ws3x3); @test Array(Cf10x8) ≈ transpose(M) * MB - mul!(Cf15x8, Mf, Adjoint(Mf_C); cache=ws3x3); @test Array(Cf15x8) ≈ M * MC' + mul!(Cf15x8, Mf, Mf_A; cache = ws3x3); @test Array(Cf15x8) ≈ M * MA + mul!(Cf10x8, Adjoint(Mf), Mf_B; cache = ws3x3); @test Array(Cf10x8) ≈ M' * MB + mul!(Cf10x8, Transpose(Mf), Mf_B; cache = ws3x3); @test Array(Cf10x8) ≈ transpose(M) * MB + mul!(Cf15x8, Mf, Adjoint(Mf_C); cache = ws3x3); @test Array(Cf15x8) ≈ M * MC' check_fmfm_fm_alloc(Mf, Mf_A, ws3x3, Cf15x8) end @@ -192,7 +192,7 @@ end c = Vector{T}(undef, 15) mul!(c, Mf, b); @test c ≈ M * b ws_vec = FactoredMatrices.Workspace(Mf, 1) - mul!(c, Mf, b; cache=ws_vec); @test c ≈ M * b + mul!(c, Mf, b; cache = ws_vec); @test c ≈ M * b check_vec_alloc(Mf, b, ws_vec, c) end end