Use LazyArrays to simplify linear algebra

.travis.yml

@@ -4,7 +4,6 @@ os:
   - linux
   - osx
 julia:
-  - 0.6
   - 0.7
   - 1.0
   - nightly

REQUIRE

@@ -1,3 +1,3 @@
-julia 0.6
-Compat 0.66
-FillArrays 0.1
+julia 0.7
+FillArrays 0.2
+LazyArrays 0.1.1

docs/src/index.md

@@ -104,6 +104,11 @@ interface consists of the following:
 | `isbanded(A)`    | Override to return `true` |
 | `inbands_getindex(A, k, j)` | Unsafe: return `A[k,j]`, without the need to check if we are inside the bands |
 | `inbands_setindex!(A, v, k, j)` | Unsafe: set `A[k,j] = v`, without the need to check if we are inside the bands |
+| --------------- | --------------- |
+| Optional methods | Brief description |
+| --------------- | --------------- |
+| `BandedMatrices.MemoryLayout(A)` | Override to get banded lazy linear algebra, e.g. `y .= Mul(A,x)` |
+| `BandedMatrices.bandeddata(A)` | Override to return a matrix of the entries in BLAS format. Required if `MemoryLayout(A)` returns `BandedColumnMajor` |
 
 Note that certain `SubArray`s of `BandedMatrix` are also banded matrices.
 The banded matrix interface is implemented for such `SubArray`s to take advantage of this.

examples/clarrays.jl

@@ -20,7 +20,7 @@ function expliciteuler(L, u₀, n)
     u = copy(u₀)
     v = copy(u)
     for _ = 1:n
-        A_mul_B!(v, L, u)
+        mul!(v, L, u)
         u,v = v,u
     end
     u

src/BandedMatrices.jl

@@ -1,78 +1,42 @@
 __precompile__()
 
 module BandedMatrices
-using Base, Compat, FillArrays
-if VERSION ≥ v"0.7-"
-    using LinearAlgebra, SparseArrays, Compat.Random
-    using LinearAlgebra.LAPACK
-    import Base: axes1, getproperty, iterate
-    import LinearAlgebra: BlasInt,
-                        BlasReal,
-                        BlasFloat,
-                        BlasComplex,
-                        axpy!,
-                        copy_oftype,
-                        checksquare,
-                        adjoint,
-                        transpose
-   import LinearAlgebra.BLAS: libblas
-   import LinearAlgebra.LAPACK: liblapack
-   import LinearAlgebra: cholesky, cholesky!, norm, diag, eigvals!, eigvals,
-                qr, axpy!, ldiv!, mul!
-   import SparseArrays: sparse
-
-   const lufact = LinearAlgebra.lu # TODO: Remove once 0.6 is dropped
-else
-    import Base.LinAlg: BlasInt,
-                        BlasReal,
-                        BlasFloat,
-                        BlasComplex,
-                        axpy!,
-                        A_mul_B!,
-                        Ac_mul_B,
-                        Ac_mul_B!,
-                        A_mul_Bc,
-                        A_mul_Bc!,
-                        Ac_mul_Bc,
-                        Ac_mul_Bc!,
-                        At_mul_B,
-                        At_mul_B!,
-                        A_mul_Bt,
-                        A_mul_Bt!,
-                        At_mul_Bt,
-                        At_mul_Bt!,
-                        A_ldiv_B!,
-                        At_ldiv_B!,
-                        Ac_ldiv_B!,
-                        copy_oftype,
-                        checksquare
-   using Base.LAPACK
-   import Base.BLAS: libblas
-   import Base.LAPACK: liblapack
-   import Base: lufact, cholfact, cholfact!, norm, diag, eigvals!, eigvals,
-                At_mul_B, Ac_mul_B, A_mul_B!, qr, qrfact, axpy!
-   import Base: sparse, indices1, indices
-
-   rmul!(A::AbstractArray, b::Number) = scale!(A, b)
-   lmul!(a::Number, B::AbstractArray) = scale!(a, B)
-   parentindices(A) = parentindexes(A)
-   const axes1 = indices1
-   const adjoint = Base.ctranspose
-   const cholesky = Base.cholfact
-   const cholesky! = Base.cholfact!
-end
+using Base, FillArrays, LazyArrays
+
+using LinearAlgebra, SparseArrays, Random
+using LinearAlgebra.LAPACK
+import Base: axes1, getproperty, iterate
+import LinearAlgebra: BlasInt,
+                    BlasReal,
+                    BlasFloat,
+                    BlasComplex,
+                    axpy!,
+                    copy_oftype,
+                    checksquare,
+                    adjoint,
+                    transpose,
+                    AdjOrTrans
+import LinearAlgebra.BLAS: libblas
+import LinearAlgebra.LAPACK: liblapack
+import LinearAlgebra: cholesky, cholesky!, norm, diag, eigvals!, eigvals,
+            qr, axpy!, ldiv!, mul!, lu, lu!
+import SparseArrays: sparse
 
 import Base: getindex, setindex!, *, +, -, ==, <, <=, >,
                 >=, /, ^, \, transpose, showerror, reindex, checkbounds, @propagate_inbounds
 
 import Base: convert, size, view, unsafe_indices,
                 first, last, size, length, unsafe_length, step,
-                to_indices, to_index, show, fill!, copy!, promote_op,
+                to_indices, to_index, show, fill!, promote_op,
                 MultiplicativeInverses, OneTo, ReshapedArray,
                                similar, copy, convert, promote_rule, rand,
-                            IndexStyle, real, imag, Slice, pointer
+                            IndexStyle, real, imag, Slice, pointer, unsafe_convert, copyto!
 
+import Base.Broadcast: BroadcastStyle, AbstractArrayStyle, DefaultArrayStyle, Broadcasted, broadcasted
 
+import LazyArrays: MemoryLayout, blasmul!, @blasmatvec, @blasmatmat, @lazymul,
+                    AbstractStridedLayout, AbstractColumnMajor,
+                    _copyto!, BMatVec, BMixedMatVec, BMixedMatMat, transposelayout, ConjLayout, conjlayout
 import FillArrays: AbstractFill
 
 export BandedMatrix,
@@ -94,16 +58,18 @@ export BandedMatrix,
        Eye
 
 
-
-include("memorylayout.jl")
 include("blas.jl")
 include("lapack.jl")
 
+
 include("generic/AbstractBandedMatrix.jl")
+include("generic/broadcast.jl")
+include("generic/matmul.jl")
 include("generic/Band.jl")
 include("generic/utils.jl")
 include("generic/interface.jl")
 
+
 include("banded/BandedMatrix.jl")
 include("banded/BandedLU.jl")
 include("banded/BandedQR.jl")
@@ -118,8 +84,4 @@ include("interfaceimpl.jl")
 
 include("deprecate.jl")
 
-
-# include("precompile.jl")
-# _precompile_()
-
 end #module

src/banded/BandedLU.jl

@@ -23,7 +23,7 @@ size(A::BandedLU, k::Integer) = k <= 0 ? error("dimension out of range") :
                                 k == 2 ? size(A.data, 2) : 1
 
 # LU factorisation with pivoting. This makes a copy!
-function lufact(A::BandedMatrix{T}) where {T<:Number}
+function lu(A::BandedMatrix{T}) where {T<:Number}
     # copy to a blas type that allows calculation of the factorisation in place.
     S = _promote_to_blas_type(T, T)
     # copy into larger array of size (2l+u*1)×n, i.e. l additional rows
@@ -33,38 +33,40 @@ function lufact(A::BandedMatrix{T}) where {T<:Number}
     data, ipiv = gbtrf!(A.l, A.u, m, data)
     BandedLU{S}(data, ipiv, A.l, A.u, m)
 end
-lufact(F::BandedLU) = F # no op
+lu(F::BandedLU) = F # no op
 
-if VERSION ≥ v"0.7-"
-    adjoint(F::BandedLU) = Adjoint(F)
-    transpose(F::BandedLU) = Transpose(F)
+lu(A::AbstractBandedMatrix) = lu(convert(BandedMatrix, A))
+lu(A::Transpose{<:Any,<:AbstractBandedMatrix}) = lu(convert(BandedMatrix, A))
+lu(A::Adjoint{<:Any,<:AbstractBandedMatrix}) = lu(convert(BandedMatrix, A))
 
-    # TODO: Finish. There was something weird about the pivoting
-    # function Base.getproperty(F::BandedLU{T}, d::Symbol) where T
-    #     if d == :L
-    #         m, n = size(F)
-    #         l,u = F.l, F.u
-    #         return _BandedMatrix([Ones{T}(1,n); view(F.data,l+u+2:2l+u+1, :)], m, l,0)
-    #     elseif d == :U
-    #         m, n = size(F)
-    #         l,u = F.l, F.u
-    #         return _BandedMatrix(F.data[l+1:l+u+1, :], n, 0, u)
-    #     elseif d == :p
-    #         return LinearAlgebra.ipiv2perm(getfield(F, :ipiv), m)
-    #     elseif d == :P
-    #         m, n = size(F)
-    #         return Matrix{T}(I, m, m)[:,invperm(F.p)]
-    #     else
-    #         getfield(F, d)
-    #     end
-    # end
-    #
-    # # iteration for destructuring into components
-    # Base.iterate(S::BandedLU) = (S.L, Val(:U))
-    # Base.iterate(S::BandedLU, ::Val{:U}) = (S.U, Val(:p))
-    # Base.iterate(S::BandedLU, ::Val{:p}) = (S.p, Val(:done))
-    # Base.iterate(S::BandedLU, ::Val{:done}) = nothing
-end
+adjoint(F::BandedLU) = Adjoint(F)
+transpose(F::BandedLU) = Transpose(F)
+
+# TODO: Finish. There was something weird about the pivoting
+# function Base.getproperty(F::BandedLU{T}, d::Symbol) where T
+#     if d == :L
+#         m, n = size(F)
+#         l,u = F.l, F.u
+#         return _BandedMatrix([Ones{T}(1,n); view(F.data,l+u+2:2l+u+1, :)], m, l,0)
+#     elseif d == :U
+#         m, n = size(F)
+#         l,u = F.l, F.u
+#         return _BandedMatrix(F.data[l+1:l+u+1, :], n, 0, u)
+#     elseif d == :p
+#         return LinearAlgebra.ipiv2perm(getfield(F, :ipiv), m)
+#     elseif d == :P
+#         m, n = size(F)
+#         return Matrix{T}(I, m, m)[:,invperm(F.p)]
+#     else
+#         getfield(F, d)
+#     end
+# end
+#
+# # iteration for destructuring into components
+# Base.iterate(S::BandedLU) = (S.L, Val(:U))
+# Base.iterate(S::BandedLU, ::Val{:U}) = (S.U, Val(:p))
+# Base.iterate(S::BandedLU, ::Val{:p}) = (S.p, Val(:done))
+# Base.iterate(S::BandedLU, ::Val{:done}) = nothing
 
 ## Utilities
 
@@ -84,32 +86,15 @@ function _promote_to_blas_type(::Type{T}, ::Type{S}) where {T<:Number, S<:Number
 end
 
 # Converts A and b to the narrowest blas type
-if VERSION < v"0.7-"
-    for typ in (BandedMatrix, BandedLU)
-        @eval function _convert_to_blas_type(A::$typ, B::AbstractVecOrMat{T}) where {T<:Number}
-            TS = _promote_to_blas_type(eltype(A), eltype(B))
-            AA = convert($typ{TS}, A)
-            BB = convert(Array{TS, ndims(B)}, B)
-            AA, BB # one of these two might make a copy
-        end
-    end
-else
-    for typ in (BandedMatrix, BandedLU, (Transpose{V,BandedLU{V}} where V), (Adjoint{V,BandedLU{V}} where V))
-        @eval function _convert_to_blas_type(A::$typ, B::AbstractVecOrMat{T}) where {T<:Number}
-            TS = _promote_to_blas_type(eltype(A), eltype(B))
-            AA = convert($typ{TS}, A)
-            BB = convert(Array{TS, ndims(B)}, B)
-            AA, BB # one of these two might make a copy
-        end
+for typ in (BandedMatrix, BandedLU, (Transpose{V,BandedLU{V}} where V), (Adjoint{V,BandedLU{V}} where V))
+    @eval function _convert_to_blas_type(A::$typ, B::AbstractVecOrMat{T}) where {T<:Number}
+        TS = _promote_to_blas_type(eltype(A), eltype(B))
+        AA = convert($typ{TS}, A)
+        BB = convert(Array{TS, ndims(B)}, B)
+        AA, BB # one of these two might make a copy
     end
 end
 
-
 # basic interface
-if VERSION < v"0.7-"
-    (\)(A::Union{BandedLU{T}, BandedMatrix{T}}, B::StridedVecOrMat{T}) where {T<:BlasFloat} =
-        A_ldiv_B!(A, copy(B)) # makes a copy
-else
-    (\)(A::Union{BandedLU{T}, BandedMatrix{T}}, B::StridedVecOrMat{T}) where {T<:BlasFloat} =
-        ldiv!(A, copy(B)) # makes a copy
-end
+(\)(A::Union{BandedLU{T}, BandedMatrix{T}}, B::AbstractVecOrMat{T}) where {T<:BlasFloat} =
+    ldiv!(A, copy(B)) # makes a copy