Fix issues from API review

.git-blame-ignore-revs

@@ -0,0 +1,2 @@
+# runic formatting
+43e2848

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

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