Allow wrapping of vectors, extend Kronecker product (#136)

Author: dkarrasch

Project.toml

@@ -1,6 +1,6 @@
 name = "LinearMaps"
 uuid = "7a12625a-238d-50fd-b39a-03d52299707e"
-version = "3.2.4"
+version = "3.3.0"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

docs/src/history.md

@@ -1,5 +1,13 @@
 # Version history
 
+## What's new in v3.3
+
+* `AbstractVector`s can now be wrapped by a `LinearMap` just like `AbstractMatrix``
+  typed objects. Upon wrapping, there are not implicitly reshaped to matrices. This
+  feature might be helpful, for instance, in the lazy representation of rank-1
+  operators `kron(LinearMap(u), v') == ⊗(u, v') == u ⊗ v'` for vectors `u` and `v`.
+  The action on vectors,`(u⊗v')*x`, is implemented optimally via `u*(v'x)`.
+
 ## What's new in v3.2
 
 * In-place left-multiplication `mul!(Y, X, A::LinearMap)` is now allowed for

src/LinearMaps.jl

@@ -34,9 +34,9 @@ MulStyle(::FiveArg, ::ThreeArg) = ThreeArg()
 MulStyle(::ThreeArg, ::ThreeArg) = ThreeArg()
 MulStyle(::LinearMap) = ThreeArg() # default
 @static if VERSION ≥ v"1.3.0-alpha.115"
-    MulStyle(::AbstractMatrix) = FiveArg()
+    MulStyle(::AbstractVecOrMat) = FiveArg()
 else
-    MulStyle(::AbstractMatrix) = ThreeArg()
+    MulStyle(::AbstractVecOrMat) = ThreeArg()
 end
 MulStyle(A::LinearMap, As::LinearMap...) = MulStyle(MulStyle(A), MulStyle(As...))
 
@@ -68,8 +68,8 @@ function check_dim_mul(C, A, B)
     return nothing
 end
 
-# conversion of AbstractMatrix to LinearMap
-convert_to_lmaps_(A::AbstractMatrix) = LinearMap(A)
+# conversion of AbstractVecOrMat to LinearMap
+convert_to_lmaps_(A::AbstractVecOrMat) = LinearMap(A)
 convert_to_lmaps_(A::LinearMap) = A
 convert_to_lmaps() = ()
 convert_to_lmaps(A) = (convert_to_lmaps_(A),)
@@ -231,8 +231,8 @@ function _generic_mapmat_mul!(Y, A, X, α=true, β=false)
     return Y
 end
 
-_unsafe_mul!(y, A::MapOrMatrix, x) = mul!(y, A, x)
-_unsafe_mul!(y, A::AbstractMatrix, x, α, β) = mul!(y, A, x, α, β)
+_unsafe_mul!(y, A::MapOrVecOrMat, x) = mul!(y, A, x)
+_unsafe_mul!(y, A::AbstractVecOrMat, x, α, β) = mul!(y, A, x, α, β)
 function _unsafe_mul!(y::AbstractVecOrMat, A::LinearMap, x::AbstractVector, α, β)
     return _generic_mapvec_mul!(y, A, x, α, β)
 end
@@ -261,11 +261,11 @@ include("show.jl") # show methods for LinearMap objects
 
 """
     LinearMap(A::LinearMap; kwargs...)::WrappedMap
-    LinearMap(A::AbstractMatrix; kwargs...)::WrappedMap
+    LinearMap(A::AbstractVecOrMat; kwargs...)::WrappedMap
     LinearMap(J::UniformScaling, M::Int)::UniformScalingMap
     LinearMap{T=Float64}(f, [fc,], M::Int, N::Int = M; kwargs...)::FunctionMap
 
-Construct a linear map object, either from an existing `LinearMap` or `AbstractMatrix` `A`,
+Construct a linear map object, either from an existing `LinearMap` or `AbstractVecOrMat` `A`,
 with the purpose of redefining its properties via the keyword arguments `kwargs`;
 a `UniformScaling` object `J` with specified (square) dimension `M`; from a `Number`
 object to lazily represent filled matrices; or
@@ -293,7 +293,7 @@ For the function-based constructor, there is one more keyword argument:
     The default value is guessed by looking at the number of arguments of the first
     occurrence of `f` in the method table.
 """
-LinearMap(A::MapOrMatrix; kwargs...) = WrappedMap(A; kwargs...)
+LinearMap(A::MapOrVecOrMat; kwargs...) = WrappedMap(A; kwargs...)
 LinearMap(J::UniformScaling, M::Int) = UniformScalingMap(J.λ, M)
 LinearMap(f, M::Int; kwargs...) = LinearMap{Float64}(f, M; kwargs...)
 LinearMap(f, M::Int, N::Int; kwargs...) = LinearMap{Float64}(f, M, N; kwargs...)

src/blockmap.jl

@@ -38,15 +38,15 @@ function rowcolranges(maps, rows)
     rowstart = 1
     for (i, row) in enumerate(rows)
         mapind += 1
-        rowend = rowstart + size(maps[mapind], 1) - 1
+        rowend = rowstart + Int(size(maps[mapind], 1))::Int - 1
         rowranges = Base.setindex(rowranges, rowstart:rowend, i)
         colstart = 1
-        colend = size(maps[mapind], 2)
+        colend = Int(size(maps[mapind], 2))::Int
         colranges = Base.setindex(colranges, colstart:colend, mapind)
         for colind in 2:row
             mapind += 1
             colstart = colend + 1
-            colend += size(maps[mapind], 2)
+            colend += Int(size(maps[mapind], 2))::Int
             colranges = Base.setindex(colranges, colstart:colend, mapind)
         end
         rowstart = rowend + 1
@@ -226,9 +226,7 @@ function check_dim(A, dim, n)
     return nothing
 end
 
-promote_to_lmaps_(n::Int, dim, A::AbstractMatrix) = (check_dim(A, dim, n); LinearMap(A))
-promote_to_lmaps_(n::Int, dim, A::AbstractVector) =
-    (check_dim(A, dim, n); LinearMap(reshape(A, length(A), 1)))
+promote_to_lmaps_(n::Int, dim, A::AbstractVecOrMat) = (check_dim(A, dim, n); LinearMap(A))
 promote_to_lmaps_(n::Int, dim, J::UniformScaling) = UniformScalingMap(J.λ, n)
 promote_to_lmaps_(n::Int, dim, A::LinearMap) = (check_dim(A, dim, n); A)
 promote_to_lmaps(n, k, dim) = ()

src/conversion.jl

@@ -47,10 +47,10 @@ SparseArrays.SparseMatrixCSC(A::LinearMap) = sparse(A)
 # special cases
 
 # ScaledMap
-Base.Matrix{T}(A::ScaledMap{<:Any, <:Any, <:MatrixMap}) where {T} =
+Base.Matrix{T}(A::ScaledMap{<:Any, <:Any, <:VecOrMatMap}) where {T} =
     convert(Matrix{T}, A.λ * A.lmap.lmap)
-SparseArrays.sparse(A::ScaledMap{<:Any, <:Any, <:MatrixMap}) =
-    convert(SparseMatrixCSC, A.λ * A.lmap.lmap)
+SparseArrays.sparse(A::ScaledMap{<:Any, <:Any, <:VecOrMatMap}) =
+    A.λ * sparse(A.lmap.lmap)
 
 # UniformScalingMap
 Base.Matrix{T}(J::UniformScalingMap) where {T} = Matrix{T}(J.λ*I, size(J))
@@ -59,6 +59,8 @@ Base.convert(::Type{AbstractMatrix}, J::UniformScalingMap) = Diagonal(fill(J.λ,
 # WrappedMap
 Base.Matrix{T}(A::WrappedMap) where {T} = Matrix{T}(A.lmap)
 Base.convert(::Type{T}, A::WrappedMap) where {T<:Matrix} = convert(T, A.lmap)
+Base.Matrix{T}(A::VectorMap) where {T} = copyto!(Matrix{eltype(T)}(undef, size(A)), A.lmap)
+Base.convert(::Type{T}, A::VectorMap) where {T<:Matrix} = T(A)
 Base.convert(::Type{AbstractMatrix}, A::WrappedMap) = convert(AbstractMatrix, A.lmap)
 SparseArrays.sparse(A::WrappedMap) = sparse(A.lmap)
 Base.convert(::Type{SparseMatrixCSC}, A::WrappedMap) = convert(SparseMatrixCSC, A.lmap)
@@ -70,19 +72,19 @@ for (T, t) in ((AdjointMap, adjoint), (TransposeMap, transpose))
 end
 
 # LinearCombination
-function Base.Matrix{T}(ΣA::LinearCombination{<:Any, <:Tuple{Vararg{MatrixMap}}}) where {T}
+function Base.Matrix{T}(ΣA::LinearCombination{<:Any, <:Tuple{Vararg{VecOrMatMap}}}) where {T}
     maps = ΣA.maps
     mats = map(A->getfield(A, :lmap), maps)
     return Matrix{T}(sum(mats))
 end
-function SparseArrays.sparse(ΣA::LinearCombination{<:Any, <:Tuple{Vararg{MatrixMap}}})
+function SparseArrays.sparse(ΣA::LinearCombination{<:Any, <:Tuple{Vararg{VecOrMatMap}}})
     maps = ΣA.maps
     mats = map(A->getfield(A, :lmap), maps)
     return convert(SparseMatrixCSC, sum(mats))
 end
 
 # CompositeMap
-function Base.Matrix{T}(AB::CompositeMap{<:Any, <:Tuple{MatrixMap, LinearMap}}) where {T}
+function Base.Matrix{T}(AB::CompositeMap{<:Any, <:Tuple{VecOrMatMap, LinearMap}}) where {T}
     B, A = AB.maps
     require_one_based_indexing(B)
     Y = Matrix{eltype(AB)}(undef, size(AB))
@@ -91,36 +93,36 @@ function Base.Matrix{T}(AB::CompositeMap{<:Any, <:Tuple{MatrixMap, LinearMap}})
     end
     return Y
 end
-for ((TA, fieldA), (TB, fieldB)) in (((MatrixMap, :lmap), (MatrixMap, :lmap)),
-                                     ((MatrixMap, :lmap), (UniformScalingMap, :λ)),
-                                     ((UniformScalingMap, :λ), (MatrixMap, :lmap)))
+for ((TA, fieldA), (TB, fieldB)) in (((VecOrMatMap, :lmap), (VecOrMatMap, :lmap)),
+                                     ((VecOrMatMap, :lmap), (UniformScalingMap, :λ)),
+                                     ((UniformScalingMap, :λ), (VecOrMatMap, :lmap)))
     @eval function Base.convert(::Type{AbstractMatrix},
                                 AB::CompositeMap{<:Any,<:Tuple{$TB,$TA}})
         B, A = AB.maps
         return A.$fieldA*B.$fieldB
     end
 end
-function Base.Matrix{T}(AB::CompositeMap{<:Any, <:Tuple{MatrixMap, MatrixMap}}) where {T}
+function Base.Matrix{T}(AB::CompositeMap{<:Any, <:Tuple{VecOrMatMap, VecOrMatMap}}) where {T}
     B, A = AB.maps
     return convert(Matrix{T}, A.lmap*B.lmap)
 end
-function SparseArrays.sparse(AB::CompositeMap{<:Any, <:Tuple{MatrixMap, MatrixMap}})
+function SparseArrays.sparse(AB::CompositeMap{<:Any, <:Tuple{VecOrMatMap, VecOrMatMap}})
     B, A = AB.maps
     return convert(SparseMatrixCSC, A.lmap*B.lmap)
 end
-function Base.Matrix{T}(λA::CompositeMap{<:Any, <:Tuple{MatrixMap, UniformScalingMap}}) where {T}
+function Base.Matrix{T}(λA::CompositeMap{<:Any, <:Tuple{VecOrMatMap, UniformScalingMap}}) where {T}
     A, J = λA.maps
     return convert(Matrix{T}, J.λ*A.lmap)
 end
-function SparseArrays.sparse(λA::CompositeMap{<:Any, <:Tuple{MatrixMap, UniformScalingMap}})
+function SparseArrays.sparse(λA::CompositeMap{<:Any, <:Tuple{VecOrMatMap, UniformScalingMap}})
     A, J = λA.maps
     return convert(SparseMatrixCSC, J.λ*A.lmap)
 end
-function Base.Matrix{T}(Aλ::CompositeMap{<:Any, <:Tuple{UniformScalingMap, MatrixMap}}) where {T}
+function Base.Matrix{T}(Aλ::CompositeMap{<:Any, <:Tuple{UniformScalingMap, VecOrMatMap}}) where {T}
     J, A = Aλ.maps
     return convert(Matrix{T}, A.lmap*J.λ)
 end
-function SparseArrays.sparse(Aλ::CompositeMap{<:Any, <:Tuple{UniformScalingMap, MatrixMap}})
+function SparseArrays.sparse(Aλ::CompositeMap{<:Any, <:Tuple{UniformScalingMap, VecOrMatMap}})
     J, A = Aλ.maps
     return convert(SparseMatrixCSC, A.lmap*J.λ)
 end

src/kronecker.jl

@@ -64,18 +64,18 @@ Base.kron(A::KroneckerMap, B::ScaledMap) = kron(A, B.lmap) * B.λ
 # generic definitions
 Base.kron(A::LinearMap, B::LinearMap, C::LinearMap, Ds::LinearMap...) =
     kron(kron(A, B), C, Ds...)
-Base.kron(A::AbstractMatrix, B::LinearMap) = kron(LinearMap(A), B)
-Base.kron(A::LinearMap, B::AbstractMatrix) = kron(A, LinearMap(B))
+Base.kron(A::AbstractVecOrMat, B::LinearMap) = kron(LinearMap(A), B)
+Base.kron(A::LinearMap, B::AbstractVecOrMat) = kron(A, LinearMap(B))
 # promote AbstractMatrix arguments to LinearMaps, then take LinearMap-Kronecker product
 for k in 3:8 # is 8 sufficient?
-    Is = ntuple(n->:($(Symbol(:A, n))::AbstractMatrix), Val(k-1))
+    Is = ntuple(n->:($(Symbol(:A, n))::AbstractVecOrMat), Val(k-1))
     # yields (:A1, :A2, :A3, ..., :A(k-1))
     L = :($(Symbol(:A, k))::LinearMap)
     # yields :Ak::LinearMap
     mapargs = ntuple(n -> :(LinearMap($(Symbol(:A, n)))), Val(k-1))
     # yields (:LinearMap(A1), :LinearMap(A2), ..., :LinearMap(A(k-1)))
 
-    @eval Base.kron($(Is...), $L, As::MapOrMatrix...) =
+    @eval Base.kron($(Is...), $L, As::MapOrVecOrMat...) =
         kron($(mapargs...), $(Symbol(:A, k)), convert_to_lmaps(As...)...)
 end
 
@@ -139,7 +139,7 @@ end
     !isone(A.λ) && rmul!(y, A.λ)
     return y
 end
-@inline function _kronmul!(y, B, x, A::MatrixMap, _)
+@inline function _kronmul!(y, B, x, A::VecOrMatMap, _)
     ma, na = size(A)
     mb, nb = size(B)
     X = reshape(x, (nb, na))
@@ -153,10 +153,13 @@ end
     elseif nb*ma <= mb*na
         _unsafe_mul!(Y, B, X * At)
     else
-        _unsafe_mul!(Y, Matrix(B*X), At)
+        _unsafe_mul!(Y, Matrix(B * X), At)
     end
     return y
 end
+const VectorMap{T} = WrappedMap{T,<:AbstractVector}
+const AdjOrTransVectorMap{T} = WrappedMap{T,<:LinearAlgebra.AdjOrTransAbsVec}
+@inline _kronmul!(y, B::AdjOrTransVectorMap, x, a::VectorMap, _) = mul!(y, a.lmap, B.lmap * x)
 
 #################
 # multiplication with vectors

src/show.jl

@@ -8,7 +8,7 @@ end
 Base.show(io::IO, A::LinearMap) = print(io, map_show(io, A, 0))
 
 map_show(io::IO, A::LinearMap, i) = ' '^i * map_summary(A) * _show(io, A, i)
-map_show(io::IO, A::AbstractMatrix, i) = ' '^i * summary(A)
+map_show(io::IO, A::AbstractVecOrMat, i) = ' '^i * summary(A)
 _show(io::IO, ::LinearMap, _) = ""
 function _show(io::IO, A::FunctionMap{T,F,Nothing}, _) where {T,F}
     "($(A.f); ismutating=$(A._ismutating), issymmetric=$(A._issymmetric), ishermitian=$(A._ishermitian), isposdef=$(A._isposdef))"

src/wrappedmap.jl

@@ -1,43 +1,49 @@
-struct WrappedMap{T, A<:MapOrMatrix} <: LinearMap{T}
+struct WrappedMap{T, A<:MapOrVecOrMat} <: LinearMap{T}
     lmap::A
     _issymmetric::Bool
     _ishermitian::Bool
     _isposdef::Bool
 end
-function WrappedMap(lmap::MapOrMatrix{T};
-                    issymmetric::Bool = issymmetric(lmap),
-                    ishermitian::Bool = ishermitian(lmap),
-                    isposdef::Bool = isposdef(lmap)) where {T}
-    WrappedMap{T, typeof(lmap)}(lmap, issymmetric, ishermitian, isposdef)
-end
 function WrappedMap{T}(lmap::MapOrMatrix;
                         issymmetric::Bool = issymmetric(lmap),
                         ishermitian::Bool = ishermitian(lmap),
                         isposdef::Bool = isposdef(lmap)) where {T}
     WrappedMap{T, typeof(lmap)}(lmap, issymmetric, ishermitian, isposdef)
 end
+WrappedMap(lmap::MapOrMatrix{T}; kwargs...) where {T} = WrappedMap{T}(lmap; kwargs...)
+function WrappedMap{T}(lmap::AbstractVector;
+                        issymmetric::Bool = false,
+                        ishermitian::Bool = false,
+                        isposdef::Bool = false) where {T}
+    WrappedMap{T, typeof(lmap)}(lmap,
+                                length(lmap) == 1 && isreal(T),
+                                length(lmap) == 1 && isreal(T),
+                                length(lmap) == 1 && isposdef(first(lmap)))
+end
+WrappedMap(lmap::AbstractVector{T}; kwargs...) where {T} = WrappedMap{T}(lmap; kwargs...)
 
-const MatrixMap{T} = WrappedMap{T,<:AbstractMatrix}
+const VecOrMatMap{T} = WrappedMap{T,<:AbstractVecOrMat}
 
-MulStyle(A::WrappedMap) = MulStyle(A.lmap)
+MulStyle(A::VecOrMatMap) = MulStyle(A.lmap)
 
-LinearAlgebra.transpose(A::MatrixMap{T}) where {T} =
+LinearAlgebra.transpose(A::VecOrMatMap{T}) where {T} =
     WrappedMap{T}(transpose(A.lmap);
                     issymmetric = A._issymmetric,
                     ishermitian = A._ishermitian,
                     isposdef = A._isposdef)
-LinearAlgebra.adjoint(A::MatrixMap{T}) where {T} =
+LinearAlgebra.adjoint(A::VecOrMatMap{T}) where {T} =
     WrappedMap{T}(adjoint(A.lmap);
                     issymmetric = A._issymmetric,
                     ishermitian = A._ishermitian,
                     isposdef = A._isposdef)
 
-Base.:(==)(A::MatrixMap, B::MatrixMap) =
+Base.:(==)(A::VecOrMatMap, B::VecOrMatMap) =
     (eltype(A) == eltype(B) && A.lmap == B.lmap && A._issymmetric == B._issymmetric &&
      A._ishermitian == B._ishermitian && A._isposdef == B._isposdef)
 
 # properties
 Base.size(A::WrappedMap) = size(A.lmap)
+Base.size(A::WrappedMap{<:Any,<:AbstractVector}) = (Int(length(A.lmap))::Int, 1)
 LinearAlgebra.issymmetric(A::WrappedMap) = A._issymmetric
 LinearAlgebra.ishermitian(A::WrappedMap) = A._ishermitian
 LinearAlgebra.isposdef(A::WrappedMap) = A._isposdef
@@ -61,9 +67,9 @@ for (In, Out) in ((AbstractVector, AbstractVecOrMat), (AbstractMatrix, AbstractM
     end
 end
 
-mul!(Y::AbstractMatrix, X::AbstractMatrix, A::MatrixMap) = mul!(Y, X, A.lmap)
+mul!(Y::AbstractMatrix, X::AbstractMatrix, A::VecOrMatMap) = mul!(Y, X, A.lmap)
 # the following method is needed for disambiguation with left-multiplication
-mul!(Y::AbstractMatrix, X::TransposeAbsVecOrMat, A::MatrixMap) = mul!(Y, X, A.lmap)
+mul!(Y::AbstractMatrix, X::TransposeAbsVecOrMat, A::VecOrMatMap) = mul!(Y, X, A.lmap)
 
 if VERSION ≥ v"1.3.0-alpha.115"
     for (In, Out) in ((AbstractVector, AbstractVecOrMat), (AbstractMatrix, AbstractMatrix))
@@ -87,14 +93,14 @@ if VERSION ≥ v"1.3.0-alpha.115"
         end
     end
 
-    mul!(X::AbstractMatrix, Y::AbstractMatrix, A::MatrixMap, α::Number, β::Number) =
+    mul!(X::AbstractMatrix, Y::AbstractMatrix, A::VecOrMatMap, α::Number, β::Number) =
         mul!(X, Y, A.lmap, α, β)
     # the following method is needed for disambiguation with left-multiplication
-    function mul!(Y::AbstractMatrix{<:RealOrComplex}, X::AbstractMatrix{<:RealOrComplex}, A::MatrixMap{<:RealOrComplex},
+    function mul!(Y::AbstractMatrix{<:RealOrComplex}, X::AbstractMatrix{<:RealOrComplex}, A::VecOrMatMap{<:RealOrComplex},
                     α::RealOrComplex, β::RealOrComplex)
         return mul!(Y, X, A.lmap, α, β)
     end
-    function mul!(Y::AbstractMatrix{<:RealOrComplex}, X::TransposeAbsVecOrMat{<:RealOrComplex}, A::MatrixMap{<:RealOrComplex},
+    function mul!(Y::AbstractMatrix{<:RealOrComplex}, X::TransposeAbsVecOrMat{<:RealOrComplex}, A::VecOrMatMap{<:RealOrComplex},
                     α::RealOrComplex, β::RealOrComplex)
         return mul!(Y, X, A.lmap, α, β)
     end

test/blockmap.jl

@@ -26,8 +26,10 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays, BenchmarkTools, Interactive
             @test (@which [I I I A11 A11 A11]).module != LinearMaps
             @test (@which hcat(I, I, I)).module != LinearMaps
             @test (@which hcat(I, I, I, LinearMap(A11), A11, A11)).module == LinearMaps
+            maps = @inferred LinearMaps.promote_to_lmaps(ntuple(i->10, 7), 1, 1, I, I, I, LinearMap(A11), A11, A11, v)
+            @inferred LinearMaps.rowcolranges(maps, (7,))
             L = @inferred hcat(I, I, I, LinearMap(A11), A11, A11, v)
-            @test L == [I I I LinearMap(A11) LinearMap(A11) LinearMap(A11) LinearMap(reshape(v, :, 1))]
+            @test L == [I I I LinearMap(A11) LinearMap(A11) LinearMap(A11) LinearMap(v)]
             x = rand(elty, 61)
             @test L isa LinearMaps.BlockMap{elty}
             @test L * x ≈ A * x

test/conversion.jl

@@ -10,6 +10,7 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays, Quaternions
     β = rand()
     M = @inferred LinearMap(A)
     N = @inferred LinearMap(M)
+    U = @inferred LinearMap(v)
 
     @test Matrix(M) == A
     @test Array(M) == A
@@ -25,6 +26,8 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays, Quaternions
     @test convert(AbstractMatrix, M') isa Adjoint
     @test convert(Matrix, M*3I) == A*3
     @test convert(Matrix, M+M) == A + A
+    @test all(Matrix(U) .== v)
+    @test convert(Matrix{ComplexF32}, U) isa Matrix{ComplexF32}
 
     # UniformScalingMap
     J = LinearMap(α*I, 10)

test/kronecker.jl

@@ -22,6 +22,10 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays
             @test @inferred size(LK, i) == size(K, i)
         end
         @test LK isa LinearMaps.KroneckerMap{ComplexF64}
+        L = ones(3) ⊗ ones(ComplexF64, 4)'
+        v = rand(4)
+        @test Matrix(L) == ones(3,4)
+        @test L*v ≈ fill(sum(v), 3)
 
         for transform in (identity, transpose, adjoint)
             @test Matrix(transform(LK)) ≈ transform(Matrix(LK)) ≈ transform(K)