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 diff --git a/src/FactoredMatrices.jl b/src/FactoredMatrices.jl index cd264d7..c6f6ab8 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)) @@ -107,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 @@ -136,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: @@ -176,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 @@ -222,9 +265,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 +280,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..dc3b88a 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) @@ -195,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