Update for LazyArrays 0.6 (#89)

* update for LazyArrays 0.6

* update Mul usage

* avoid upper triangular ambiguity

* Fix NaN bug in Julia 1.1

* Bump LazyArrays

REQUIRE

@@ -1,3 +1,3 @@
 julia 0.7
 FillArrays 0.3
-LazyArrays 0.5
+LazyArrays 0.6

src/BandedMatrices.jl

@@ -29,11 +29,11 @@ import LazyArrays: MemoryLayout, @lazymul, @lazylmul, @lazyldiv,
                     _copyto!, MatMulVec, MatMulMat, transposelayout, triangulardata,
                     ConjLayout, conjlayout, SymmetricLayout, symmetriclayout, symmetricdata,
                     triangularlayout, MatMulVec, MatLdivVec, TriangularLayout,
-                    AbstractBandedLayout, DiagonalLayout,
+                    AbstractBandedLayout, DiagonalLayout, LayoutApplyStyle,
                     ArrayMulArrayStyle, HermitianLayout, hermitianlayout, hermitiandata,
                     MulAdd, materialize!, BlasMatMulMat, BlasMatMulVec, VcatLayout, ZerosLayout,
                     AbstractColumnMajor, MulLayout, colsupport, rowsupport,
-                    DenseColumnMajor, DenseRowMajor
+                    DenseColumnMajor, DenseRowMajor, ArrayLdivArray
 import FillArrays: AbstractFill
 
 export BandedMatrix,

src/banded/linalg.jl

@@ -2,22 +2,21 @@
 
 # Direct and transposed algorithms
 
-function _copyto!(_, dest::AbstractVecOrMat, L::Ldiv{<:BandedColumnMajor})
-    Ai, B = L.factors
-    A = parent(Ai)
+function _copyto!(_, dest::AbstractVecOrMat, L::ArrayLdivArray{<:BandedColumnMajor})
+    A, B = L.args
     checksquare(A)
     dest ≡ B || copyto!(dest, B)
     ldiv!(factorize(A), dest)
 end
 
-function _copyto!(_, dest::AbstractVecOrMat, L::Ldiv{<:BandedRowMajor})
-    Ai, B = L.factors
-    copyto!(dest, Mul(transpose(factorize(transpose(parent(Ai)))), B))
+function _copyto!(_, dest::AbstractVecOrMat, L::ArrayLdivArray{<:BandedRowMajor})
+    A, B = L.args
+    copyto!(dest, Mul(transpose(factorize(transpose(A))), B))
 end
 
-function _copyto!(_, dest::AbstractVecOrMat, L::Ldiv{<:ConjLayout{<:BandedRowMajor}})
-    Ai, B = L.factors
-    copyto!(dest, Mul(factorize(parent(Ai)')', B))
+function _copyto!(_, dest::AbstractVecOrMat, L::ArrayLdivArray{<:ConjLayout{<:BandedRowMajor}})
+    A, B = L.args
+    copyto!(dest, Mul(factorize(A')', B))
 end
 
 

src/generic/matmul.jl

@@ -3,7 +3,7 @@ BroadcastStyle(::BandedStyle, M::ArrayMulArrayStyle) = M
 
 @lazymul AbstractBandedMatrix
 
-bandwidths(M::Mul) = prodbandwidths(M.factors...)
+bandwidths(M::Mul) = prodbandwidths(M.args...)
 
 similar(M::MatMulMat{<:AbstractBandedLayout,<:AbstractBandedLayout}, ::Type{T}) where T =
     similar(M, T, axes(M))
@@ -65,7 +65,7 @@ function materialize!(M::BlasMatMulVec{<:BandedColumnMajor,<:AbstractStridedLayo
     (length(y) ≠ m || length(x) ≠ n) && throw(DimensionMismatch("*"))
     l, u = bandwidths(A)
     if -l > u # no bands
-        lmul!(β, y)
+        _fill_lmul!(β, y)
     elseif l < 0
         materialize!(MulAdd(α, view(A, :, 1-l:n), view(x, 1-l:n), β, y))
     elseif u < 0
@@ -123,23 +123,23 @@ end
 
 
 @inline function _copyto!(_, C::AbstractVector{T}, M::MatMulVec{<:BandedColumnMajor,<:AbstractStridedLayout,T,T}) where T<:BlasFloat
-    A,B = M.factors
+    A,B = M.args
     materialize!(MulAdd(one(T), A, B, zero(T), C))
 end
 
 @inline function _copyto!(_, C::AbstractVector{T}, M::MatMulVec{<:BandedRowMajor,<:AbstractStridedLayout,T,T}) where T<:BlasFloat
-    A,B = M.factors
+    A,B = M.args
     materialize!(MulAdd(one(T), A, B, zero(T), C))
 end
 
 @inline function _copyto!(_, C::AbstractVector{T}, M::MatMulVec{<:ConjLayout{<:BandedRowMajor},<:AbstractStridedLayout,T,T}) where T<:BlasFloat
-    A,B = M.factors
+    A,B = M.args
     materialize!(MulAdd(one(T), A, B, zero(T), C))
 end
 
 
 @inline function _copyto!(_, c::AbstractVector, M::MatMulVec{<:AbstractBandedLayout})
-    A,b = M.factors
+    A,b = M.args
     for k = 1:length(c)
         c[k] = zero(eltype(A)) * b[1]
     end
@@ -150,7 +150,7 @@ end
 end
 
 @inline function _copyto!(_, c::AbstractVector, M::MatMulVec{<:BandedRowMajor})
-    At,b = M.factors
+    At,b = M.args
     A = transpose(At)
     c .= zero.(c)
     @inbounds for j = 1:size(A,2), k = colrange(A,j)
@@ -160,7 +160,7 @@ end
 end
 
 @inline function _copyto!(_, c::AbstractVector, M::MatMulVec{<:ConjLayout{<:BandedRowMajor}})
-    Ac,b = M.factors
+    Ac,b = M.args
     A = Ac'
     c .= zero.(c)
     @inbounds for j = 1:size(A,2), k = colrange(A,j)
@@ -208,22 +208,27 @@ const ConjOrBandedLayout = Union{AbstractBandedLayout,ConjLayout{<:AbstractBande
 const ConjOrBandedColumnMajor = Union{<:BandedColumnMajor,ConjLayout{<:BandedColumnMajor}}
 
 @inline function _copyto!(_, C::AbstractMatrix{T}, M::MatMulMat{<:ConjOrBandedColumnMajor,<:ConjOrBandedColumnMajor,T,T}) where T<:BlasFloat
-    A,B = M.factors
+    A,B = M.args
     materialize!(MulAdd(one(T), A, B, zero(T), C))
 end
 
 @inline function _copyto!(_, C::AbstractMatrix, M::MatMulMat{<:ConjOrBandedLayout,<:ConjOrBandedLayout})
-     A, B = M.factors
+     A, B = M.args
      banded_mul!(C, A, B)
 end
 
 @inline function _copyto!(_, C::AbstractMatrix, M::MatMulMat{<:ConjOrBandedLayout})
-     A, B = M.factors
+     A, B = M.args
      banded_mul!(C, A, B)
 end
 
 @inline function _copyto!(_, C::AbstractMatrix, M::MatMulMat{<:Any,<:ConjOrBandedLayout})
-     A, B = M.factors
+     A, B = M.args
+     banded_mul!(C, A, B)
+end
+
+@inline function _copyto!(_, C::AbstractMatrix, M::MatMulMat{<:TriangularLayout,<:ConjOrBandedLayout})
+     A, B = M.args
      banded_mul!(C, A, B)
 end
 

src/interfaceimpl.jl

@@ -71,7 +71,7 @@ end
 
 function _copyto!(::VcatLayout{<:Tuple{<:Any,ZerosLayout}}, y::AbstractVector,
                  M::MatMulVec{<:AbstractBandedLayout,<:VcatLayout{<:Tuple{<:Any,ZerosLayout}}})
-    A,x = M.factors
+    A,x = M.args
     length(y) == size(A,1) || throw(DimensionMismatch())
     length(x) == size(A,2) || throw(DimensionMismatch())
 
@@ -86,7 +86,7 @@ function _copyto!(::VcatLayout{<:Tuple{<:Any,ZerosLayout}}, y::AbstractVector,
 end
 
 function similar(M::MatMulVec{<:AbstractBandedLayout,<:VcatLayout{<:Tuple{<:Any,ZerosLayout}}}, ::Type{T}) where T
-    A,x = M.factors
+    A,x = M.args
     xf,_ = x.arrays
     n = max(0,min(length(xf) + bandwidth(A,1),length(M)))
     Vcat(Vector{T}(undef, n), Zeros{T}(size(A,1)-n))
@@ -97,10 +97,10 @@ end
 # MulMatrix
 ###
 
-bandwidths(M::MulMatrix) = bandwidths(M.mul)
-isbanded(M::MulMatrix) = all(isbanded, M.mul.factors)
+bandwidths(M::MulMatrix) = bandwidths(M.applied)
+isbanded(M::MulMatrix) = all(isbanded, M.applied.args)
 
-const MulBandedMatrix{T} = MulMatrix{T, <:Mul{<:Tuple{Vararg{<:AbstractBandedLayout}}}}
+const MulBandedMatrix{T} = MulMatrix{T, <:Mul{<:LayoutApplyStyle{<:Tuple{Vararg{<:AbstractBandedLayout}}}}}
 
 BroadcastStyle(::Type{<:MulBandedMatrix}) = BandedStyle()
 
@@ -120,10 +120,10 @@ end
 
 
 getindex(M::MatMulMat{<:AbstractBandedLayout,<:AbstractBandedLayout}, k::Integer, j::Integer) =
-    _banded_mul_getindex(eltype(M), M.factors, k, j)
+    _banded_mul_getindex(eltype(M), M.args, k, j)
 
-getindex(M::Mul{<:Tuple{Vararg{<:AbstractBandedLayout}}}, k::Integer, j::Integer) =
-    _banded_mul_getindex(eltype(M), (first(M.factors), Mul(tail(M.factors)...)), k, j)
+getindex(M::Mul{<:LayoutApplyStyle{<:Tuple{Vararg{<:AbstractBandedLayout}}}}, k::Integer, j::Integer) =
+    _banded_mul_getindex(eltype(M), (first(M.args), Mul(tail(M.args)...)), k, j)
 
 
 @inline _sub_materialize(::MulLayout{<:Tuple{Vararg{<:AbstractBandedLayout}}}, V) = BandedMatrix(V)

src/tribanded.jl

@@ -46,23 +46,23 @@ Base.replace_in_print_matrix(A::Union{LowerTriangular{<:Any,<:AbstractBandedMatr
 @inline function _copyto!(::AbstractStridedLayout, dest::AbstractVector{T},
          M::MatMulVec{<:TriangularLayout{'U',UNIT,<:BandedColumnMajor},
                                    <:AbstractStridedLayout,T,T}) where {UNIT,T <: BlasFloat}
-    A,x = M.factors
+    A,x = M.args
     x ≡ dest || copyto!(dest, x)
     tbmv!('U', 'N', UNIT, size(A,1), bandwidth(A,2), tribandeddata(A), dest)
 end
 
 @inline function _copyto!(::AbstractStridedLayout, dest::AbstractVector{T},
          M::MatMulVec{<:TriangularLayout{'L',UNIT,<:BandedColumnMajor},
                                    <:AbstractStridedLayout,T,T}) where {UNIT,T <: BlasFloat}
-    A,x = M.factors
+    A,x = M.args
     x ≡ dest || copyto!(dest, x)
     tbmv!('L', 'N', UNIT, size(A,1), bandwidth(A,1), tribandeddata(A), dest)
 end
 
 @inline function _copyto!(::AbstractStridedLayout, dest::AbstractVector{T},
          M::MatMulVec{<:TriangularLayout{UPLO,UNIT,BandedRowMajor},
                                    <:AbstractStridedLayout,T,T}) where {UPLO,UNIT,T <: BlasFloat}
-    A,x = M.factors
+    A,x = M.args
     x ≡ dest || copyto!(dest, x)
     tbmv!(UPLO, 'T', UNIT, transpose(tribandeddata(A)), dest)
 end
@@ -71,7 +71,7 @@ end
 @inline function _copyto!(::AbstractStridedLayout, dest::AbstractVector{T},
          M::MatMulVec{<:TriangularLayout{UPLO,UNIT,ConjLayout{BandedRowMajor}},
                                    <:AbstractStridedLayout, T, T}) where {UPLO,UNIT,T <: BlasFloat}
-    A,x = M.factors
+    A,x = M.args
     x ≡ dest || copyto!(dest, x)
     tbmv!(UPLO, 'C', UNIT, tribandeddata(A)', dest)
 end
@@ -85,26 +85,23 @@ end
 @inline function _copyto!(::AbstractStridedLayout, dest::AbstractVector{T},
          M::MatLdivVec{<:TriangularLayout{'U',UNIT,<:BandedColumnMajor},
                                    <:AbstractStridedLayout,T,T}) where {UNIT,T <: BlasFloat}
-    Ai,x = M.factors
-    A = inv(Ai)
+    A,x = M.args
     x ≡ dest || copyto!(dest, x)
     tbsv!('U', 'N', UNIT, size(A,1), bandwidth(A,2), tribandeddata(A), dest)
 end
 
 @inline function _copyto!(::AbstractStridedLayout, dest::AbstractVector,
          M::MatLdivVec{<:TriangularLayout{'L',UNIT,<:BandedColumnMajor},
                                    <:AbstractStridedLayout,T,T}) where {UNIT,T <: BlasFloat}
-    Ai,x = M.factors
-    A = inv(Ai)
+    A,x = M.args
     x ≡ dest || copyto!(dest, x)
     tbsv!('L', 'N', UNIT, size(A,1), bandwidth(A,1), tribandeddata(A), dest)
 end
 
 @inline function _copyto!(::AbstractStridedLayout, dest::AbstractVector,
          M::MatLdivVec{<:TriangularLayout{UPLO,UNIT,BandedRowMajor},
                                    <:AbstractStridedLayout,T,T}) where {UPLO,UNIT,T <: BlasFloat}
-    Ai,x = M.factors
-    A = inv(Ai)
+    A,x = M.args
     x ≡ dest || copyto!(dest, x)
     tbsv!(UPLO, 'T', UNIT, transpose(tribandeddata(A)), dest)
 end
@@ -113,8 +110,7 @@ end
 @inline function _copyto!(::AbstractStridedLayout, dest::AbstractVector,
          M::MatLdivVec{<:TriangularLayout{UPLO,UNIT,ConjLayout{BandedRowMajor}},
                                    <:AbstractStridedLayout,T,T}) where {UPLO,UNIT,T <: BlasFloat}
-    Ai,x = M.factors
-    A = inv(Ai)
+    A,x = M.args
     x ≡ dest || copyto!(dest, x)
     tbsv!(UPLO, 'C', UNIT, tribandeddata(A)', dest)
 end

test/test_bandedlu.jl

@@ -9,60 +9,60 @@ tldiv!(A, b) = ldiv!(transpose(A), b)
 cldiv!(A, b) = ldiv!(adjoint(A), b)
 
 
-@testset "Conversion to blas type" begin
-    @testset "_promote_to_blas_type" begin
-        typ = Float64
-        @test BandedMatrices._promote_to_blas_type(typ, ComplexF64) == ComplexF64
-        @test BandedMatrices._promote_to_blas_type(typ, ComplexF32)  == ComplexF64
-        @test BandedMatrices._promote_to_blas_type(typ, Float64)    == typ
-        @test BandedMatrices._promote_to_blas_type(typ, Float32)    == typ
-        @test BandedMatrices._promote_to_blas_type(typ, Int64)      == typ
-        @test BandedMatrices._promote_to_blas_type(typ, Int32)      == typ
-        @test BandedMatrices._promote_to_blas_type(typ, Int16)      == typ
-        @test BandedMatrices._promote_to_blas_type(typ, Int8)       == typ
-
-        typ = Float32
-        @test BandedMatrices._promote_to_blas_type(typ, ComplexF64) == ComplexF64
-        @test BandedMatrices._promote_to_blas_type(typ, ComplexF32)  == ComplexF32
-        @test BandedMatrices._promote_to_blas_type(typ, Float64)    == Float64
-        @test BandedMatrices._promote_to_blas_type(typ, Float32)    == typ
-        @test BandedMatrices._promote_to_blas_type(typ, Int64)      == typ
-        @test BandedMatrices._promote_to_blas_type(typ, Int32)      == typ
-        @test BandedMatrices._promote_to_blas_type(typ, Int16)      == typ
-        @test BandedMatrices._promote_to_blas_type(typ, Int8)       == typ
-
-        @test_throws ErrorException BandedMatrices._promote_to_blas_type(_foo, Float64)
-    end
+@testset "Banded A\\b" begin
+    @testset "Conversion to blas type" begin
+        @testset "_promote_to_blas_type" begin
+            typ = Float64
+            @test BandedMatrices._promote_to_blas_type(typ, ComplexF64) == ComplexF64
+            @test BandedMatrices._promote_to_blas_type(typ, ComplexF32)  == ComplexF64
+            @test BandedMatrices._promote_to_blas_type(typ, Float64)    == typ
+            @test BandedMatrices._promote_to_blas_type(typ, Float32)    == typ
+            @test BandedMatrices._promote_to_blas_type(typ, Int64)      == typ
+            @test BandedMatrices._promote_to_blas_type(typ, Int32)      == typ
+            @test BandedMatrices._promote_to_blas_type(typ, Int16)      == typ
+            @test BandedMatrices._promote_to_blas_type(typ, Int8)       == typ
+
+            typ = Float32
+            @test BandedMatrices._promote_to_blas_type(typ, ComplexF64) == ComplexF64
+            @test BandedMatrices._promote_to_blas_type(typ, ComplexF32)  == ComplexF32
+            @test BandedMatrices._promote_to_blas_type(typ, Float64)    == Float64
+            @test BandedMatrices._promote_to_blas_type(typ, Float32)    == typ
+            @test BandedMatrices._promote_to_blas_type(typ, Int64)      == typ
+            @test BandedMatrices._promote_to_blas_type(typ, Int32)      == typ
+            @test BandedMatrices._promote_to_blas_type(typ, Int16)      == typ
+            @test BandedMatrices._promote_to_blas_type(typ, Int8)       == typ
+
+            @test_throws ErrorException BandedMatrices._promote_to_blas_type(_foo, Float64)
+        end
 
-    @testset "ldiv!" begin
-        As   = Any[_BandedMatrix(rand(1:10, 3, 5), 5, 1, 1),
-                   _BandedMatrix(rand(3, 5)*im,    5, 1, 1),
-                   _BandedMatrix(rand(3, 5),       5, 1, 1)
-                  ]
-        bs   = Any[rand(1:10,   5),
-                   rand(1:10,   5),
-                   rand(1:10.0, 5)*im,
-                  ]
-        typs = Any[Float64,
-                   ComplexF64,
-                   ComplexF64]
-
-        for (A, b, typ) in zip(As, bs, typs)
-            AA,   bb   = BandedMatrices._convert_to_blas_type(A,         b)
-            AAlu, bblu = BandedMatrices._convert_to_blas_type(lu(A), b)
-            @test eltype(AA) == eltype(bb) == eltype(AAlu) == eltype(bblu) == typ
-            @test Matrix(A)\copy(b)             ≈ A\copy(b)
-            @test Matrix(A)\copy(b)             ≈ lu(A)\copy(b)
-            @test Matrix(A)\copy(b)             ≈ ldiv!(A, copy(b))
-            @test Matrix(A)\copy(b)             ≈ ldiv!(lu(A), copy(b))
-            @test transpose(Matrix(A))\copy(b)  ≈ tldiv!(lu(A), copy(b))
-            @test transpose(Matrix(A))\copy(b)  ≈ ldiv!(lu(transpose(A)), copy(b))
-            @test Matrix(A)'\copy(b)            ≈ cldiv!(lu(A), copy(b))
+        @testset "ldiv!" begin
+            As   = Any[_BandedMatrix(rand(1:10, 3, 5), 5, 1, 1),
+                       _BandedMatrix(rand(3, 5)*im,    5, 1, 1),
+                       _BandedMatrix(rand(3, 5),       5, 1, 1)
+                      ]
+            bs   = Any[rand(1:10,   5),
+                       rand(1:10,   5),
+                       rand(1:10.0, 5)*im,
+                      ]
+            typs = Any[Float64,
+                       ComplexF64,
+                       ComplexF64]
+
+            for (A, b, typ) in zip(As, bs, typs)
+                AA,   bb   = BandedMatrices._convert_to_blas_type(A,         b)
+                AAlu, bblu = BandedMatrices._convert_to_blas_type(lu(A), b)
+                @test eltype(AA) == eltype(bb) == eltype(AAlu) == eltype(bblu) == typ
+                @test Matrix(A)\copy(b)             ≈ A\copy(b)
+                @test Matrix(A)\copy(b)             ≈ lu(A)\copy(b)
+                @test Matrix(A)\copy(b)             ≈ ldiv!(A, copy(b))
+                @test Matrix(A)\copy(b)             ≈ ldiv!(lu(A), copy(b))
+                @test transpose(Matrix(A))\copy(b)  ≈ tldiv!(lu(A), copy(b))
+                @test transpose(Matrix(A))\copy(b)  ≈ ldiv!(lu(transpose(A)), copy(b))
+                @test Matrix(A)'\copy(b)            ≈ cldiv!(lu(A), copy(b))
+            end
         end
     end
-end
 
-@testset "Banded A\b" begin
     @testset "banded" begin
         A = brand(5, 1, 1)
         b = Float64[1, 2, 3, 4, 5]