getindex for LayoutVector (#19)

* getindex for LayoutVector

* v0.3.2

* Better print

* move print_matrix_row to FillArrays.jl

* Diagonal overload improvements

* Update muladd.jl

* Tridiagonal tests

* Update Project.toml

Author: dlfivefifty

Project.toml

@@ -1,19 +1,20 @@
 name = "ArrayLayouts"
 uuid = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.3.1"
+version = "0.3.2"
 
 [deps]
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 
 [compat]
-FillArrays = "0.8"
+FillArrays = "0.8.10"
 julia = "1"
 
 [extras]
+Base64 = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "Random"]
+test = ["Test", "Random", "Base64"]

src/ArrayLayouts.jl

@@ -2,7 +2,7 @@ module ArrayLayouts
 using Base, Base.Broadcast, LinearAlgebra, FillArrays
 import LinearAlgebra.BLAS
 
-import Base: AbstractArray, AbstractMatrix, AbstractVector, 
+import Base: AbstractArray, AbstractMatrix, AbstractVector,
         ReinterpretArray, ReshapedArray, AbstractCartesianIndex, Slice,
              RangeIndex, BroadcastStyle, copyto!, length, broadcastable, axes,
              getindex, eltype, tail, IndexStyle, IndexLinear, getproperty,
@@ -33,33 +33,33 @@ import Base.Broadcast: BroadcastStyle, AbstractArrayStyle, Broadcasted, broadcas
                         combine_eltypes, DefaultArrayStyle, instantiate, materialize,
                         materialize!, eltypes
 
-import LinearAlgebra: AbstractTriangular, AbstractQ, checksquare, pinv, fill!, tilebufsize, Abuf, Bbuf, Cbuf, dot, factorize, qr, lu, cholesky, 
-                        norm2, norm1, normInf, normMinusInf, qr, lu, qr!, lu!
+import LinearAlgebra: AbstractTriangular, AbstractQ, checksquare, pinv, fill!, tilebufsize, Abuf, Bbuf, Cbuf, dot, factorize, qr, lu, cholesky,
+                        norm2, norm1, normInf, normMinusInf, qr, lu, qr!, lu!, AdjOrTrans, HermOrSym
 
 import LinearAlgebra.BLAS: BlasFloat, BlasReal, BlasComplex
 
-import FillArrays: AbstractFill, getindex_value
+import FillArrays: AbstractFill, getindex_value, axes_print_matrix_row
 
 if VERSION < v"1.2-"
     import Base: has_offset_axes
     require_one_based_indexing(A...) = !has_offset_axes(A...) || throw(ArgumentError("offset arrays are not supported but got an array with index other than 1"))
 else
-    import Base: require_one_based_indexing    
-end     
+    import Base: require_one_based_indexing
+end
 
 export materialize, materialize!, MulAdd, muladd!, Ldiv, Rdiv, Lmul, Rmul, lmul, rmul, ldiv, rdiv, mul, MemoryLayout, AbstractStridedLayout,
         DenseColumnMajor, ColumnMajor, ZerosLayout, FillLayout, AbstractColumnMajor, RowMajor, AbstractRowMajor, UnitStride,
-        DiagonalLayout, ScalarLayout, SymTridiagonalLayout, HermitianLayout, SymmetricLayout, TriangularLayout, 
+        DiagonalLayout, ScalarLayout, SymTridiagonalLayout, HermitianLayout, SymmetricLayout, TriangularLayout,
         UnknownLayout, AbstractBandedLayout, ApplyBroadcastStyle, ConjLayout, AbstractFillLayout,
         colsupport, rowsupport, layout_getindex, QLayout, LayoutArray, LayoutMatrix, LayoutVector
 
 struct ApplyBroadcastStyle <: BroadcastStyle end
-@inline function copyto!(dest::AbstractArray, bc::Broadcasted{ApplyBroadcastStyle}) 
+@inline function copyto!(dest::AbstractArray, bc::Broadcasted{ApplyBroadcastStyle})
     @assert length(bc.args) == 1
     copyto!(dest, first(bc.args))
 end
 
-# Subtypes of LayoutArray default to 
+# Subtypes of LayoutArray default to
 # ArrayLayouts routines
 abstract type LayoutArray{T,N} <: AbstractArray{T,N} end
 const LayoutMatrix{T} = LayoutArray{T,2}
@@ -76,15 +76,15 @@ strides(A::Transpose) = _transpose_strides(strides(parent(A))...)
 
 represents that the entry is the complex-conjugate of the pointed to entry.
 """
-struct ConjPtr{T} 
+struct ConjPtr{T}
     ptr::Ptr{T}
 end
 
 unsafe_convert(::Type{Ptr{T}}, A::Adjoint{<:Real}) where T<:Real = unsafe_convert(Ptr{T}, parent(A))
 unsafe_convert(::Type{Ptr{T}}, A::Transpose) where T = unsafe_convert(Ptr{T}, parent(A))
 # work-around issue with complex conjugation of pointer
 unsafe_convert(::Type{Ptr{T}}, Ac::Adjoint{<:Complex}) where T<:Complex = unsafe_convert(ConjPtr{T}, parent(Ac))
-unsafe_convert(::Type{ConjPtr{T}}, Ac::Adjoint{<:Complex}) where T<:Complex = unsafe_convert(Ptr{T}, parent(Ac))    
+unsafe_convert(::Type{ConjPtr{T}}, Ac::Adjoint{<:Complex}) where T<:Complex = unsafe_convert(Ptr{T}, parent(Ac))
 function unsafe_convert(::Type{ConjPtr{T}}, V::SubArray{T,2}) where {T,N,P}
     kr, jr = parentindices(V)
     unsafe_convert(Ptr{T}, view(parent(V)', jr, kr))
@@ -125,32 +125,34 @@ end
 
 @layoutmatrix LayoutMatrix
 
-_copyto!(_, _, dest::AbstractArray{T,N}, src::AbstractArray{V,N}) where {T,V,N} = 
+getindex(A::LayoutVector, kr::AbstractVector) = layout_getindex(A, kr)
+
+_copyto!(_, _, dest::AbstractArray{T,N}, src::AbstractArray{V,N}) where {T,V,N} =
     Base.invoke(copyto!, Tuple{AbstractArray{T,N},AbstractArray{V,N}}, dest, src)
 
-    
-copyto!(dest::LayoutArray{<:Any,N}, src::LayoutArray{<:Any,N}) where N = 
+
+copyto!(dest::LayoutArray{<:Any,N}, src::LayoutArray{<:Any,N}) where N =
     _copyto!(MemoryLayout(dest), MemoryLayout(src), dest, src)
-copyto!(dest::AbstractArray{<:Any,N}, src::LayoutArray{<:Any,N}) where N = 
+copyto!(dest::AbstractArray{<:Any,N}, src::LayoutArray{<:Any,N}) where N =
     _copyto!(MemoryLayout(dest), MemoryLayout(src), dest, src)
-copyto!(dest::LayoutArray{<:Any,N}, src::AbstractArray{<:Any,N}) where N = 
+copyto!(dest::LayoutArray{<:Any,N}, src::AbstractArray{<:Any,N}) where N =
     _copyto!(MemoryLayout(dest), MemoryLayout(src), dest, src)
 
-copyto!(dest::SubArray{<:Any,N,<:LayoutArray}, src::SubArray{<:Any,N,<:LayoutArray}) where N = 
+copyto!(dest::SubArray{<:Any,N,<:LayoutArray}, src::SubArray{<:Any,N,<:LayoutArray}) where N =
     _copyto!(MemoryLayout(dest), MemoryLayout(src), dest, src)
-copyto!(dest::SubArray{<:Any,N,<:LayoutArray}, src::LayoutArray{<:Any,N}) where N = 
+copyto!(dest::SubArray{<:Any,N,<:LayoutArray}, src::LayoutArray{<:Any,N}) where N =
     _copyto!(MemoryLayout(dest), MemoryLayout(src), dest, src)
-copyto!(dest::LayoutArray{<:Any,N}, src::SubArray{<:Any,N,<:LayoutArray}) where N = 
+copyto!(dest::LayoutArray{<:Any,N}, src::SubArray{<:Any,N,<:LayoutArray}) where N =
     _copyto!(MemoryLayout(dest), MemoryLayout(src), dest, src)
-copyto!(dest::SubArray{<:Any,N,<:LayoutArray}, src::AbstractArray{<:Any,N}) where N = 
+copyto!(dest::SubArray{<:Any,N,<:LayoutArray}, src::AbstractArray{<:Any,N}) where N =
+    _copyto!(MemoryLayout(dest), MemoryLayout(src), dest, src)
+copyto!(dest::AbstractArray{<:Any,N}, src::SubArray{<:Any,N,<:LayoutArray}) where N =
     _copyto!(MemoryLayout(dest), MemoryLayout(src), dest, src)
-copyto!(dest::AbstractArray{<:Any,N}, src::SubArray{<:Any,N,<:LayoutArray}) where N = 
-    _copyto!(MemoryLayout(dest), MemoryLayout(src), dest, src)    
 
 
 
 zero!(A::AbstractArray{T}) where T = fill!(A,zero(T))
-function zero!(A::AbstractArray{<:AbstractArray}) 
+function zero!(A::AbstractArray{<:AbstractArray})
     for a in A
         zero!(a)
     end
@@ -213,4 +215,35 @@ end
     return A
 end
 
+###
+# printing
+###
+
+layout_replace_in_print_matrix(_, A, i, j, s) =
+    i in colsupport(A,j) ? s : Base.replace_with_centered_mark(s)
+
+Base.replace_in_print_matrix(A::Union{LayoutMatrix,
+                                      UpperTriangular{<:Any,<:LayoutMatrix},
+                                      UnitUpperTriangular{<:Any,<:LayoutMatrix},
+                                      LowerTriangular{<:Any,<:LayoutMatrix},
+                                      UnitLowerTriangular{<:Any,<:LayoutMatrix},
+                                      AdjOrTrans{<:Any,<:LayoutMatrix},
+                                      HermOrSym{<:Any,<:LayoutMatrix},
+                                      SubArray{<:Any,2,<:LayoutMatrix}}, i::Integer, j::Integer, s::AbstractString) =
+    layout_replace_in_print_matrix(MemoryLayout(A), A, i, j, s)                
+
+Base.print_matrix_row(io::IO,
+        X::Union{LayoutMatrix,
+        LayoutVector,
+        AbstractTriangular{<:Any,<:LayoutMatrix},
+        AdjOrTrans{<:Any,<:LayoutMatrix},
+        HermOrSym{<:Any,<:LayoutMatrix},
+        SubArray{<:Any,2,<:LayoutMatrix},
+        Diagonal{<:Any,<:LayoutVector}}, A::Vector,
+        i::Integer, cols::AbstractVector, sep::AbstractString) =
+        axes_print_matrix_row(axes(X), io, X, A, i, cols, sep)
+
+
+
+
 end

src/diagonal.jl

@@ -16,8 +16,8 @@ end
 
 copy(M::Lmul{<:DiagonalLayout,<:DiagonalLayout}) = Diagonal(M.A.diag .* M.B.diag)
 copy(M::Lmul{<:DiagonalLayout}) = M.A.diag .* M.B
-copy(M::Lmul{<:DiagonalLayout{<:AbstractFillLayout}}) = getindex_value(M.A.diag) .* M.B
-copy(M::Lmul{<:DiagonalLayout{<:AbstractFillLayout},<:DiagonalLayout}) = Diagonal(getindex_value(M.A.diag) .* M.B.diag)
+copy(M::Lmul{<:DiagonalLayout{<:AbstractFillLayout}}) = getindex_value(M.A.diag) * M.B
+copy(M::Lmul{<:DiagonalLayout{<:AbstractFillLayout},<:DiagonalLayout}) = Diagonal(getindex_value(M.A.diag) * M.B.diag)
 
 # Diagonal multiplication never changes structure
 similar(M::Rmul{<:Any,<:DiagonalLayout}, ::Type{T}, axes) where T = similar(M.A, T, axes)

src/memorylayout.jl

@@ -497,9 +497,11 @@ triangulardata(A::SubArray{<:Any,2,<:Any,<:Tuple{<:Union{Slice,Base.OneTo},<:Uni
 
 
 abstract type AbstractBandedLayout <: MemoryLayout end
+abstract type AbstractTridiagonalLayout <: MemoryLayout end
 
 struct DiagonalLayout{ML} <: AbstractBandedLayout end
-struct SymTridiagonalLayout{ML} <: AbstractBandedLayout end
+struct SymTridiagonalLayout{ML} <: AbstractTridiagonalLayout end
+struct TridiagonalLayout{ML} <: AbstractTridiagonalLayout end
 
 diagonallayout(_) = DiagonalLayout{UnknownLayout}()
 diagonallayout(::ML) where ML<:AbstractStridedLayout = DiagonalLayout{ML}()
@@ -508,13 +510,30 @@ diagonaldata(D::Diagonal) = parent(D)
 
 MemoryLayout(::Type{SymTridiagonal{T,P}}) where {T,P} = SymTridiagonalLayout{typeof(MemoryLayout(P))}()
 diagonaldata(D::SymTridiagonal) = D.dv
-offdiagonaldata(D::SymTridiagonal) = D.ev
+supdiagonaldata(D::SymTridiagonal) = D.ev
+subdiagonaldata(D::SymTridiagonal) = D.ev
+
+MemoryLayout(::Type{Tridiagonal{T,P}}) where {T,P} = TridiagonalLayout{typeof(MemoryLayout(P))}()
+diagonaldata(D::Tridiagonal) = D.d
+subdiagonaldata(D::Tridiagonal) = D.dl
+supdiagonaldata(D::Tridiagonal) = D.du
+
 
 transposelayout(ml::DiagonalLayout) = ml
 transposelayout(ml::SymTridiagonalLayout) = ml
+transposelayout(ml::TridiagonalLayout) = ml
 transposelayout(ml::ConjLayout{DiagonalLayout}) = ml
 
 adjointlayout(::Type{<:Real}, ml::SymTridiagonalLayout) = ml
+adjointlayout(::Type{<:Real}, ml::TridiagonalLayout) = ml
+
+subdiagonaldata(D::Transpose) = supdiagonaldata(parent(D))
+supdiagonaldata(D::Transpose) = subdiagonaldata(parent(D))
+diagonaldata(D::Transpose) = diagonaldata(parent(D))
+
+subdiagonaldata(D::Adjoint{<:Real}) = supdiagonaldata(parent(D))
+supdiagonaldata(D::Adjoint{<:Real}) = subdiagonaldata(parent(D))
+diagonaldata(D::Adjoint{<:Real}) = diagonaldata(parent(D))
 
 ###
 # Fill

src/muladd.jl

@@ -216,6 +216,14 @@ function default_blasmul!(α, A::AbstractMatrix, B::AbstractVector, β, C::Abstr
 end
 
 function materialize!(M::MatMulMatAdd)
+    α, A, B, β, C = M.α, M.A, M.B, M.β, M.C
+    if C ≡ B
+        B = copy(B)
+    end
+    default_blasmul!(α, A, B, iszero(β) ? false : β, C)
+end
+
+function materialize!(M::MatMulMatAdd{<:AbstractStridedLayout,<:AbstractStridedLayout,<:AbstractStridedLayout})
     α, A, B, β, C = M.α, M.A, M.B, M.β, M.C
     if C ≡ B
         B = copy(B)
@@ -379,11 +387,12 @@ function materialize!(M::MulAdd{<:Any,<:DiagonalLayout{<:AbstractFillLayout}})
     M.C
 end
 
-copy(M::MulAdd{<:DiagonalLayout{<:AbstractFillLayout}}) = (M.α * getindex_value(M.A.diag)) .* M.B .+ M.β .* M.C
-copy(M::MulAdd{<:DiagonalLayout{<:AbstractFillLayout},<:Any,ZerosLayout}) = (M.α * getindex_value(M.A.diag)) .* M.B
+copy(M::MulAdd{<:DiagonalLayout{<:AbstractFillLayout}}) = (M.α * getindex_value(M.A.diag)) * M.B .+ M.β * M.C
+copy(M::MulAdd{<:DiagonalLayout{<:AbstractFillLayout},<:Any,ZerosLayout}) = (M.α * getindex_value(M.A.diag)) * M.B
 copy(M::MulAdd{<:AbstractFillLayout,<:AbstractFillLayout,<:AbstractFillLayout}) = M.α*M.A*M.B + M.β*M.C
-copy(M::MulAdd{<:Any,<:DiagonalLayout{<:AbstractFillLayout}}) = (M.α * getindex_value(M.B.diag)) .* M.A .+ M.β .* M.C
-copy(M::MulAdd{<:Any,<:DiagonalLayout{<:AbstractFillLayout},ZerosLayout}) = (M.α * getindex_value(M.B.diag)) .* M.A
+copy(M::MulAdd{<:Any,<:DiagonalLayout{<:AbstractFillLayout}}) = (M.α * getindex_value(M.B.diag)) * M.A .+ M.β * M.C
+copy(M::MulAdd{<:Any,<:DiagonalLayout{<:AbstractFillLayout},ZerosLayout}) = (M.α * getindex_value(M.B.diag)) * M.A
+
 
 BroadcastStyle(::Type{<:MulAdd}) = ApplyBroadcastStyle()
 

test/runtests.jl

@@ -1,4 +1,4 @@
-using ArrayLayouts, Random, FillArrays, Test
+using ArrayLayouts, Random, FillArrays, Test, Base64
 import ArrayLayouts: MemoryLayout, @_layoutlmul, triangulardata
 
 Random.seed!(0)
@@ -18,51 +18,73 @@ Base.strides(A::MyMatrix) = strides(A.A)
 Base.unsafe_convert(::Type{Ptr{T}}, A::MyMatrix) where T = Base.unsafe_convert(Ptr{T}, A.A)
 MemoryLayout(::Type{MyMatrix}) = DenseColumnMajor()
 
-@testset "LayoutMatrix" begin
-    A = MyMatrix(randn(5,5))
-    for (kr,jr) in ((1:2,2:3), (:,:), (:,1:2), (2:3,:), ([1,2],3:4), (:,[1,2]), ([2,3],:))
-        @test A[kr,jr] == A.A[kr,jr]
-    end
-    b = randn(5)
-    for Tri in (UpperTriangular, UnitUpperTriangular, LowerTriangular, UnitLowerTriangular)
-        @test ldiv!(Tri(A), copy(b)) ≈ ldiv!(Tri(A.A), copy(b))
-        @test lmul!(Tri(A), copy(b)) ≈ lmul!(Tri(A.A), copy(b))
+struct MyVector <: LayoutVector{Float64}
+    A::Vector{Float64}
+end
+
+Base.getindex(A::MyVector, k::Int) = A.A[k]
+Base.setindex!(A::MyVector, v, k::Int) = setindex!(A.A, v, k)
+Base.size(A::MyVector) = size(A.A)
+Base.strides(A::MyVector) = strides(A.A)
+Base.unsafe_convert(::Type{Ptr{T}}, A::MyVector) where T = Base.unsafe_convert(Ptr{T}, A.A)
+MemoryLayout(::Type{MyVector}) = DenseColumnMajor()
+
+@testset "LayoutArray" begin
+    @testset "LayoutVector" begin
+        A = MyVector(randn(5))
+        @test A == A.A == Vector(A)
+        @test A[1:3] == A.A[1:3]
     end
 
-    @test copyto!(MyMatrix(Array{Float64}(undef,5,5)), A) == A
-    @test copyto!(Array{Float64}(undef,5,5), A) == A
-    @test copyto!(MyMatrix(Array{Float64}(undef,5,5)), A.A) == A
-    @test copyto!(view(MyMatrix(Array{Float64}(undef,5,5)),1:3,1:3), view(A,1:3,1:3)) == A[1:3,1:3]
-    @test copyto!(view(MyMatrix(Array{Float64}(undef,5,5)),:,:), A) == A
-    @test copyto!(MyMatrix(Array{Float64}(undef,3,3)), view(A,1:3,1:3)) == A[1:3,1:3]
-    @test copyto!(view(MyMatrix(Array{Float64}(undef,5,5)),:,:), A.A) == A
-    @test copyto!(Array{Float64}(undef,3,3), view(A,1:3,1:3)) == A[1:3,1:3]
-
-    @test qr(A).factors ≈ qr(A.A).factors
-    @test qr(A,Val(true)).factors ≈ qr(A.A,Val(true)).factors
-    @test lu(A).factors ≈ lu(A.A).factors
-    @test lu(A,Val(true)).factors ≈ lu(A.A,Val(true)).factors
-    @test_throws ErrorException qr!(A)
-    @test_throws ErrorException lu!(A)
-
-    @test qr(A) isa LinearAlgebra.QRCompactWY
-    @test inv(A) ≈ inv(A.A)
-
-    Bin = randn(5,5)
-    B = MyMatrix(copy(Bin))
-    muladd!(1.0, A, A, 2.0, B)
-    @test all(B .=== A.A^2 + 2Bin)
-
-    # tiled_blasmul!
-    B = MyMatrix(copy(Bin))
-    muladd!(1.0, Ones(5,5), A, 2.0, B)
-
-    #generic_blasmul!
-    A = BigFloat.(randn(5,5))
-    Bin = BigFloat.(randn(5,5))
-    B = copy(Bin)
-    muladd!(1.0, Ones(5,5), A, 2.0, B)
-    @test B == Ones(5,5)*A + 2.0Bin
+    @testset "LayoutMatrix" begin
+        A = MyMatrix(randn(5,5))
+        for (kr,jr) in ((1:2,2:3), (:,:), (:,1:2), (2:3,:), ([1,2],3:4), (:,[1,2]), ([2,3],:))
+            @test A[kr,jr] == A.A[kr,jr]
+        end
+        b = randn(5)
+        for Tri in (UpperTriangular, UnitUpperTriangular, LowerTriangular, UnitLowerTriangular)
+            @test ldiv!(Tri(A), copy(b)) ≈ ldiv!(Tri(A.A), copy(b))
+            @test lmul!(Tri(A), copy(b)) ≈ lmul!(Tri(A.A), copy(b))
+        end
+
+        @test copyto!(MyMatrix(Array{Float64}(undef,5,5)), A) == A
+        @test copyto!(Array{Float64}(undef,5,5), A) == A
+        @test copyto!(MyMatrix(Array{Float64}(undef,5,5)), A.A) == A
+        @test copyto!(view(MyMatrix(Array{Float64}(undef,5,5)),1:3,1:3), view(A,1:3,1:3)) == A[1:3,1:3]
+        @test copyto!(view(MyMatrix(Array{Float64}(undef,5,5)),:,:), A) == A
+        @test copyto!(MyMatrix(Array{Float64}(undef,3,3)), view(A,1:3,1:3)) == A[1:3,1:3]
+        @test copyto!(view(MyMatrix(Array{Float64}(undef,5,5)),:,:), A.A) == A
+        @test copyto!(Array{Float64}(undef,3,3), view(A,1:3,1:3)) == A[1:3,1:3]
+
+        @test qr(A).factors ≈ qr(A.A).factors
+        @test qr(A,Val(true)).factors ≈ qr(A.A,Val(true)).factors
+        @test lu(A).factors ≈ lu(A.A).factors
+        @test lu(A,Val(true)).factors ≈ lu(A.A,Val(true)).factors
+        @test_throws ErrorException qr!(A)
+        @test_throws ErrorException lu!(A)
+
+        @test qr(A) isa LinearAlgebra.QRCompactWY
+        @test inv(A) ≈ inv(A.A)
+
+        Bin = randn(5,5)
+        B = MyMatrix(copy(Bin))
+        muladd!(1.0, A, A, 2.0, B)
+        @test all(B .=== A.A^2 + 2Bin)
+
+        # tiled_blasmul!
+        B = MyMatrix(copy(Bin))
+        muladd!(1.0, Ones(5,5), A, 2.0, B)
+
+        #generic_blasmul!
+        A = BigFloat.(randn(5,5))
+        Bin = BigFloat.(randn(5,5))
+        B = copy(Bin)
+        muladd!(1.0, Ones(5,5), A, 2.0, B)
+        @test B == Ones(5,5)*A + 2.0Bin
+
+        C = MyMatrix([1 2; 3 4])
+        @test stringmime("text/plain", C) == "2×2 MyMatrix:\n 1.0  2.0\n 3.0  4.0"
+    end
 end
 
 struct MyUpperTriangular{T} <: AbstractMatrix{T}