add matrix powers

Author: dlfivefifty

src/InfiniteLinearAlgebra.jl

@@ -1,17 +1,17 @@
 module InfiniteLinearAlgebra
 using BlockArrays, BlockBandedMatrices, BandedMatrices, LazyArrays, FillArrays, InfiniteArrays, MatrixFactorizations, LinearAlgebra
 
-import Base: +, -, *, /, \, OneTo, getindex, promote_op, _unsafe_getindex, print_matrix_row, size,
+import Base: +, -, *, /, \, ^, OneTo, getindex, promote_op, _unsafe_getindex, print_matrix_row, size,
             AbstractMatrix, AbstractArray, Matrix, Array, Vector, AbstractVector,
             show, getproperty
+import Base.Broadcast: BroadcastStyle
+
 import InfiniteArrays: OneToInf, InfUnitRange, Infinity, InfStepRange
 import FillArrays: AbstractFill
 import BandedMatrices: BandedMatrix, _BandedMatrix, bandeddata
-import LinearAlgebra: lmul!,  ldiv!, matprod, qr, QRPackedQ, AbstractTriangular, AbstractQ, adjoint, transpose,
-                        QR
-import LazyArrays: CachedArray, DenseColumnMajor, FillLayout, ApplyMatrix, check_mul_axes
-import MatrixFactorizations: ql, ql!, QLPackedQ, getL, reflector!, reflectorApply!,
-                            QL
+import LinearAlgebra: lmul!,  ldiv!, matprod, qr, QRPackedQ, AbstractTriangular, AbstractQ, adjoint, transpose, QR
+import LazyArrays: CachedArray, DenseColumnMajor, FillLayout, ApplyMatrix, check_mul_axes, ApplyStyle, LazyArrayApplyStyle, LazyArrayStyle
+import MatrixFactorizations: ql, ql!, QLPackedQ, getL, reflector!, reflectorApply!, QL
 
 import BlockArrays: BlockSizes, cumulsizes, _find_block, AbstractBlockVecOrMat, sizes_from_blocks
 
@@ -21,6 +21,24 @@ import BlockBandedMatrices: _BlockSkylineMatrix, _BandedMatrix, AbstractBlockSiz
         BlockSkylineSizes, BlockSkylineMatrix, BlockBandedMatrix, _BlockBandedMatrix, BlockTridiagonal
 
 
+# Fix ∞ BandedMatrix
+ApplyStyle(::typeof(*), ::Type{<:BandedMatrix{<:Any,<:Any,<:OneToInf}}, _::Type{<:AbstractArray}...) =
+    LazyArrayApplyStyle()
+
+# BroadcastStyle(::Type{<:BandedMatrix{<:Any,<:Any,<:OneToInf}}) = LazyArrayStyle{2}()
+
+^(A::BandedMatrix{T,<:Any,<:OneToInf}, p::Integer) where T =
+    if p < 0 
+        inv(A)^(-p)
+    elseif p == 0
+        Eye{T}(∞)
+    elseif p == 1
+        copy(A)
+    else
+        A*A^(p-1)
+    end
+    
+
 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"))

test/runtests.jl

@@ -40,7 +40,13 @@ import BandedMatrices: bandeddata, _BandedMatrix
     At = BandedMatrix(transpose(A))
     @test Ac[1:10,1:10] ≈ (A')[1:10,1:10] ≈ A[1:10,1:10]'
     @test At[1:10,1:10] ≈ transpose(A)[1:10,1:10] ≈ transpose(A[1:10,1:10])
+
+    A = _BandedMatrix(Fill(1,4,∞),∞,1,2)
+    @test A*A isa ApplyArray
+    @test (A^2)[1:10,1:10] == (A*A)[1:10,1:10] == (A[1:100,1:100]^2)[1:10,1:10]
+    @test (A^3)[1:10,1:10] == (A*A*A)[1:10,1:10] == (A[1:100,1:100]^3)[1:10,1:10]
 end
+
 include("test_hessenbergq.jl")
 include("test_infql.jl")
 

test/test_infqr.jl

@@ -4,8 +4,10 @@ import BandedMatrices: _BandedMatrix, colsupport
 
 A = _BandedMatrix(Vcat(Ones(1,∞), (1:∞)', Ones(1,∞)), ∞, 1, 1)
 C = cache(A)
-V = view(C.data,1:10,1:11)
-V isa BandedMatrices.BandedSubBandedMatrix
+V = view(C.data,:,1:11)
+@V isa BandedMatrices.BandedSubBandedMatrix
+
+
 qr!(V)
 qr!(C.data)