Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 7 additions & 5 deletions .github/workflows/ci.yml
Original file line number Diff line number Diff line change
@@ -1,7 +1,9 @@
name: CI
on:
- push
- pull_request
push:
branches:
- master
pull_request:
jobs:
test:
name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
Expand All @@ -10,9 +12,9 @@ jobs:
fail-fast: false
matrix:
version:
- '1.5'
- '1.6'
- '1'
- '^1.7.0-0'
- '^1.8.0-0'
os:
- ubuntu-latest
- macOS-latest
Expand Down Expand Up @@ -62,4 +64,4 @@ jobs:
- run: julia --project=docs docs/make.jl
env:
GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }}
DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }}
9 changes: 5 additions & 4 deletions Project.toml
Original file line number Diff line number Diff line change
@@ -1,22 +1,23 @@
name = "ArrayLayouts"
uuid = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
authors = ["Sheehan Olver <[email protected]>"]
version = "0.7.10"
version = "0.8"

[deps]
FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[compat]
FillArrays = "0.12.6, 0.13"
julia = "1.5"
FillArrays = "0.13.1"
julia = "1.6"

[extras]
Base64 = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[targets]
test = ["Test", "Random", "Base64", "Compat"]
test = ["Base64", "Compat", "Random", "StableRNGs", "Test"]
30 changes: 24 additions & 6 deletions src/ArrayLayouts.jl
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@ import Base: AbstractArray, AbstractMatrix, AbstractVector,
AbstractMatrix, AbstractArray, checkindex, unsafe_length, OneTo, one, zero,
to_shape, _sub2ind, print_matrix, print_matrix_row, print_matrix_vdots,
checkindex, Slice, @propagate_inbounds, @_propagate_inbounds_meta,
_in_range, _range, _rangestyle, Ordered,
_in_range, _range, Ordered,
ArithmeticWraps, floatrange, reverse, unitrange_last,
AbstractArray, AbstractVector, axes, (:), _sub2ind_recurse, broadcast, promote_eltypeof,
similar, @_gc_preserve_end, @_gc_preserve_begin,
Expand All @@ -35,12 +35,12 @@ import Base.Broadcast: BroadcastStyle, AbstractArrayStyle, Broadcasted, broadcas
materialize!, eltypes

import LinearAlgebra: AbstractTriangular, AbstractQ, checksquare, pinv, fill!, tilebufsize, factorize, qr, lu, cholesky,
norm2, norm1, normInf, normMinusInf, qr, lu, qr!, lu!, AdjOrTrans, HermOrSym, copy_oftype,
AdjointAbsVec, TransposeAbsVec, cholcopy, checknonsingular, _apply_ipiv_rows!, ipiv2perm, RealHermSymComplexHerm, chkfullrank
norm2, norm1, normInf, normMinusInf, qr, lu, qr!, lu!, AdjOrTrans, HermOrSym, AdjointAbsVec,
TransposeAbsVec, cholcopy, checknonsingular, _apply_ipiv_rows!, ipiv2perm, RealHermSymComplexHerm, chkfullrank

import LinearAlgebra.BLAS: BlasFloat, BlasReal, BlasComplex

import FillArrays: AbstractFill, getindex_value, axes_print_matrix_row
import FillArrays: AbstractFill, getindex_value, axes_print_matrix_row, _copy_oftype

import Base: require_one_based_indexing

Expand All @@ -58,6 +58,14 @@ if VERSION < v"1.7-"
const RowMaximum = Val{true}
const NoPivot = Val{false}
end

if VERSION < v"1.8-"
const CRowMaximum = Val{true}
const CNoPivot = Val{false}
else
const CRowMaximum = RowMaximum
const CNoPivot = NoPivot
end


struct ApplyBroadcastStyle <: BroadcastStyle end
Expand All @@ -81,7 +89,9 @@ strides(A::Transpose) = _transpose_strides(strides(parent(A))...)
"""
ConjPtr{T}

represents that the entry is the complex-conjugate of the pointed to entry.
A memory address referring to complex conjugated data of type T. However, there is no guarantee
that the memory is actually valid, or that it actually represents the complex conjugate of data of
the specified type.
"""
struct ConjPtr{T}
ptr::Ptr{T}
Expand All @@ -96,6 +106,14 @@ function unsafe_convert(::Type{ConjPtr{T}}, V::SubArray{T,2}) where {T,N,P}
unsafe_convert(Ptr{T}, view(parent(V)', jr, kr))
end

Base.elsize(::Type{<:Adjoint{<:Complex,P}}) where P<:AbstractVecOrMat = conjelsize(P)
conjelsize(::Type{<:Adjoint{<:Complex,P}}) where P<:AbstractVecOrMat = Base.elsize(P)
conjelsize(::Type{<:Transpose{<:Any, P}}) where {P<:AbstractVecOrMat} = conjelsize(P)
conjelsize(::Type{<:PermutedDimsArray{<:Any, <:Any, <:Any, <:Any, P}}) where {P} = conjelsize(P)
conjelsize(::Type{<:ReshapedArray{<:Any,<:Any,P}}) where {P} = conjelsize(P)
conjelsize(::Type{<:SubArray{<:Any,<:Any,P}}) where {P} = conjelsize(P)
conjelsize(A::AbstractArray) = conjelsize(typeof(A))

include("memorylayout.jl")
include("mul.jl")
include("muladd.jl")
Expand Down Expand Up @@ -367,4 +385,4 @@ Base.typed_vcat(::Type{T}, A::LayoutVecOrMats, B::LayoutVecOrMats, C::AbstractVe
Base.typed_hcat(::Type{T}, A::LayoutVecOrMats, B::LayoutVecOrMats, C::AbstractVecOrMat...) where T = typed_hcat(T, A, B, C...)
Base.typed_vcat(::Type{T}, A::AbstractVecOrMat, B::LayoutVecOrMats, C::AbstractVecOrMat...) where T = typed_vcat(T, A, B, C...)
Base.typed_hcat(::Type{T}, A::AbstractVecOrMat, B::LayoutVecOrMats, C::AbstractVecOrMat...) where T = typed_hcat(T, A, B, C...)
end # module
end # module
22 changes: 11 additions & 11 deletions src/diagonal.jl
Original file line number Diff line number Diff line change
Expand Up @@ -56,11 +56,11 @@ copy(M::Lmul{<:DiagonalLayout,<:TridiagonalLayout}) = M.A * convert(Tridiagonal,
copy(M::Rmul{<:SymTridiagonalLayout,<:DiagonalLayout}) = convert(SymTridiagonal, M.A) * M.B
copy(M::Lmul{<:DiagonalLayout,<:SymTridiagonalLayout}) = M.A * convert(SymTridiagonal, M.B)

copy(M::Lmul{DiagonalLayout{OnesLayout}}) = copy_oftype(M.B, eltype(M))
copy(M::Lmul{DiagonalLayout{OnesLayout},<:DiagonalLayout}) = Diagonal(copy_oftype(diagonaldata(M.B), eltype(M)))
copy(M::Lmul{<:DiagonalLayout,DiagonalLayout{OnesLayout}}) = Diagonal(copy_oftype(diagonaldata(M.A), eltype(M)))
copy(M::Lmul{DiagonalLayout{OnesLayout},DiagonalLayout{OnesLayout}}) = copy_oftype(M.B, eltype(M))
copy(M::Rmul{<:Any,DiagonalLayout{OnesLayout}}) = copy_oftype(M.A, eltype(M))
copy(M::Lmul{DiagonalLayout{OnesLayout}}) = _copy_oftype(M.B, eltype(M))
copy(M::Lmul{DiagonalLayout{OnesLayout},<:DiagonalLayout}) = Diagonal(_copy_oftype(diagonaldata(M.B), eltype(M)))
copy(M::Lmul{<:DiagonalLayout,DiagonalLayout{OnesLayout}}) = Diagonal(_copy_oftype(diagonaldata(M.A), eltype(M)))
copy(M::Lmul{DiagonalLayout{OnesLayout},DiagonalLayout{OnesLayout}}) = _copy_oftype(M.B, eltype(M))
copy(M::Rmul{<:Any,DiagonalLayout{OnesLayout}}) = _copy_oftype(M.A, eltype(M))

copy(M::Lmul{<:DiagonalLayout{<:AbstractFillLayout}}) = getindex_value(diagonaldata(M.A)) * M.B
copy(M::Lmul{<:DiagonalLayout{<:AbstractFillLayout},<:DiagonalLayout}) = getindex_value(diagonaldata(M.A)) * M.B
Expand All @@ -74,9 +74,9 @@ copy(M::Rmul{<:SymTridiagonalLayout,<:DiagonalLayout{<:AbstractFillLayout}}) = M
copy(M::Lmul{<:DiagonalLayout{<:AbstractFillLayout},<:SymTridiagonalLayout}) = getindex_value(diagonaldata(M.A)) * M.B


copy(M::Rmul{<:BidiagonalLayout,DiagonalLayout{OnesLayout}}) = copy_oftype(M.A, eltype(M))
copy(M::Lmul{DiagonalLayout{OnesLayout},<:BidiagonalLayout}) = copy_oftype(M.B, eltype(M))
copy(M::Rmul{<:TridiagonalLayout,DiagonalLayout{OnesLayout}}) = copy_oftype(M.A, eltype(M))
copy(M::Lmul{DiagonalLayout{OnesLayout},<:TridiagonalLayout}) = copy_oftype(M.B, eltype(M))
copy(M::Rmul{<:SymTridiagonalLayout,DiagonalLayout{OnesLayout}}) = copy_oftype(M.A, eltype(M))
copy(M::Lmul{DiagonalLayout{OnesLayout},<:SymTridiagonalLayout}) = copy_oftype(M.B, eltype(M))
copy(M::Rmul{<:BidiagonalLayout,DiagonalLayout{OnesLayout}}) = _copy_oftype(M.A, eltype(M))
copy(M::Lmul{DiagonalLayout{OnesLayout},<:BidiagonalLayout}) = _copy_oftype(M.B, eltype(M))
copy(M::Rmul{<:TridiagonalLayout,DiagonalLayout{OnesLayout}}) = _copy_oftype(M.A, eltype(M))
copy(M::Lmul{DiagonalLayout{OnesLayout},<:TridiagonalLayout}) = _copy_oftype(M.B, eltype(M))
copy(M::Rmul{<:SymTridiagonalLayout,DiagonalLayout{OnesLayout}}) = _copy_oftype(M.A, eltype(M))
copy(M::Lmul{DiagonalLayout{OnesLayout},<:SymTridiagonalLayout}) = _copy_oftype(M.B, eltype(M))
44 changes: 27 additions & 17 deletions src/factorizations.jl
Original file line number Diff line number Diff line change
Expand Up @@ -305,20 +305,32 @@ _qr!(layout, axes, A, args...; kwds...) = error("Overload _qr!(::$(typeof(layout
_lu(layout, axes, A; kwds...) = Base.invoke(lu, Tuple{AbstractMatrix{eltype(A)}}, A; kwds...)
_lu(layout, axes, A, pivot::P; kwds...) where P = Base.invoke(lu, Tuple{AbstractMatrix{eltype(A)},P}, A, pivot; kwds...)
_lu!(layout, axes, A, args...; kwds...) = error("Overload _lu!(::$(typeof(layout)), axes, A)")
_cholesky(layout, axes, A, ::Val{false}=Val(false); check::Bool = true) = cholesky!(cholcopy(A); check = check)
_cholesky(layout, axes, A, ::Val{true}; tol = 0.0, check::Bool = true) = cholesky!(cholcopy(A), Val(true); tol = tol, check = check)
_cholesky(layout, axes, A, ::CNoPivot=CNoPivot(); check::Bool = true) = cholesky!(cholcopy(A); check = check)
_cholesky(layout, axes, A, ::CRowMaximum; tol = 0.0, check::Bool = true) = cholesky!(cholcopy(A), CRowMaximum(); tol = tol, check = check)
_factorize(layout, axes, A) = qr(A) # Default to QR


_factorize(::AbstractStridedLayout, axes, A) = lu(A)
function _lu!(::AbstractColumnMajor, axes, A::AbstractMatrix{T}, pivot::Union{NoPivot, RowMaximum} = RowMaximum();
check::Bool = true) where T<:BlasFloat
if pivot === NoPivot()
return generic_lufact!(A, pivot; check = check)
if VERSION < v"1.8-"
function _lu!(::AbstractColumnMajor, axes, A::AbstractMatrix{T}, pivot::Union{NoPivot, RowMaximum} = RowMaximum();
check::Bool = true) where T<:BlasFloat
if pivot === NoPivot()
return generic_lufact!(A, pivot; check = check)
end
lpt = LAPACK.getrf!(A)
check && checknonsingular(lpt[3])
return LU{T,typeof(A)}(lpt[1], lpt[2], lpt[3])
end
else
function _lu!(::AbstractColumnMajor, axes, A::AbstractMatrix{T}, pivot::Union{NoPivot, RowMaximum} = RowMaximum();
check::Bool = true) where T<:BlasFloat
if pivot === NoPivot()
return generic_lufact!(A, pivot; check = check)
end
lpt = LAPACK.getrf!(A)
check && checknonsingular(lpt[3])
return LU{T,typeof(A),typeof(lpt[2])}(lpt[1], lpt[2], lpt[3])
end
lpt = LAPACK.getrf!(A)
check && checknonsingular(lpt[3])
return LU{T,typeof(A)}(lpt[1], lpt[2], lpt[3])
end

# for some reason only defined for StridedMatrix in LinearAlgebra
Expand Down Expand Up @@ -355,14 +367,14 @@ end

_chol!(_, A, UL) = LinearAlgebra._chol!(A, UL)

function _cholesky!(layout, axes, A::RealHermSymComplexHerm, ::Val{false}; check::Bool = true)
function _cholesky!(layout, axes, A::RealHermSymComplexHerm, ::CNoPivot; check::Bool = true)
C, info = _chol!(layout, A.data, A.uplo == 'U' ? UpperTriangular : LowerTriangular)
check && LinearAlgebra.checkpositivedefinite(info)
return Cholesky(C.data, A.uplo, info)
end

function _cholesky!(::SymmetricLayout{<:AbstractColumnMajor}, axes, A::AbstractMatrix{<:BlasReal},
::Val{true}; tol = 0.0, check::Bool = true)
::CRowMaximum; tol = 0.0, check::Bool = true)
AA, piv, rank, info = LAPACK.pstrf!(A.uplo, A.data, tol)
C = CholeskyPivoted{eltype(AA),typeof(AA)}(AA, A.uplo, piv, rank, tol, info)
check && chkfullrank(C)
Expand Down Expand Up @@ -390,9 +402,10 @@ end

macro _layoutfactorizations(Typ)
esc(quote
LinearAlgebra.cholesky(A::$Typ, args...; kwds...) = ArrayLayouts._cholesky(ArrayLayouts.MemoryLayout(A), axes(A), A, args...; kwds...)
LinearAlgebra.cholesky!(A::LinearAlgebra.RealHermSymComplexHerm{<:Real,<:$Typ}, v::Val{false}=Val(false); check::Bool = true) = ArrayLayouts._cholesky!(ArrayLayouts.MemoryLayout(A), axes(A), A, v; check=check)
LinearAlgebra.cholesky!(A::LinearAlgebra.RealHermSymComplexHerm{<:Real,<:$Typ}, v::Val{true}; check::Bool = true, tol = 0.0) = ArrayLayouts._cholesky!(ArrayLayouts.MemoryLayout(A), axes(A), A, v; check=check, tol=tol)
LinearAlgebra.cholesky(A::$Typ, v::CNoPivot = CNoPivot(); kwds...) = ArrayLayouts._cholesky(ArrayLayouts.MemoryLayout(A), axes(A), A, v; kwds...)
LinearAlgebra.cholesky(A::$Typ, v::CRowMaximum; kwds...) = ArrayLayouts._cholesky(ArrayLayouts.MemoryLayout(A), axes(A), A, v; kwds...)
LinearAlgebra.cholesky!(A::LinearAlgebra.RealHermSymComplexHerm{<:Real,<:$Typ}, v::CNoPivot = CNoPivot(); check::Bool = true) = ArrayLayouts._cholesky!(ArrayLayouts.MemoryLayout(A), axes(A), A, v; check=check)
LinearAlgebra.cholesky!(A::LinearAlgebra.RealHermSymComplexHerm{<:Real,<:$Typ}, v::CRowMaximum; check::Bool = true, tol = 0.0) = ArrayLayouts._cholesky!(ArrayLayouts.MemoryLayout(A), axes(A), A, v; check=check, tol=tol)
LinearAlgebra.qr(A::$Typ, args...; kwds...) = ArrayLayouts._qr(ArrayLayouts.MemoryLayout(A), axes(A), A, args...; kwds...)
LinearAlgebra.qr!(A::$Typ, args...; kwds...) = ArrayLayouts._qr!(ArrayLayouts.MemoryLayout(A), axes(A), A, args...; kwds...)
LinearAlgebra.lu(A::$Typ, pivot::Union{ArrayLayouts.NoPivot,ArrayLayouts.RowMaximum}; kwds...) = ArrayLayouts._lu(ArrayLayouts.MemoryLayout(A), axes(A), A, pivot; kwds...)
Expand Down Expand Up @@ -424,6 +437,3 @@ end

LinearAlgebra.ldiv!(L::LU{<:Any,<:LayoutMatrix}, B::LayoutVector) = ArrayLayouts.ldiv!(L, B)

if VERSION ≥ v"1.7-"
@deprecate qr(A::LayoutMatrix, ::Val{true}) qr(A, ColumnNorm())
end
2 changes: 1 addition & 1 deletion src/mul.jl
Original file line number Diff line number Diff line change
Expand Up @@ -309,4 +309,4 @@ for Typ in (:LayoutMatrix, :(Symmetric{<:Any,<:LayoutMatrix}), :(Hermitian{<:Any
@inline -(A::$Typ, Λ::UniformScaling) = _apply(MemoryLayout(A), size(A), -, A, Λ)
@inline -(Λ::UniformScaling, A::$Typ) = _apply(MemoryLayout(A), size(A), -, Λ, A)
end
end
end
4 changes: 2 additions & 2 deletions src/muladd.jl
Original file line number Diff line number Diff line change
Expand Up @@ -74,7 +74,7 @@ copy(M::MulAdd) = copyto!(similar(M), M)
_fill_copyto!(dest, C) = copyto!(dest, C)
_fill_copyto!(dest, C::Zeros) = zero!(dest) # exploit special fill! overload

@inline copyto!(dest::AbstractArray{T}, M::MulAdd) where T =
@inline copyto!(dest::AbstractArray{T}, M::MulAdd) where T =
muladd!(M.α, unalias(dest,M.A), unalias(dest,M.B), M.β, _fill_copyto!(dest, M.C))

# Modified from LinearAlgebra._generic_matmatmul!
Expand Down Expand Up @@ -420,4 +420,4 @@ end
const ZerosLayouts = Union{ZerosLayout,DualLayout{ZerosLayout}}
copy(M::MulAdd{<:ZerosLayouts, <:ZerosLayouts, <:ZerosLayouts}) = M.C
copy(M::MulAdd{<:ZerosLayouts, <:Any, <:ZerosLayouts}) = M.C
copy(M::MulAdd{<:Any, <:ZerosLayouts, <:ZerosLayouts}) = M.C
copy(M::MulAdd{<:Any, <:ZerosLayouts, <:ZerosLayouts}) = M.C
11 changes: 5 additions & 6 deletions test/test_layoutarray.jl
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
using ArrayLayouts, LinearAlgebra, FillArrays, Base64, Test
import ArrayLayouts: sub_materialize, MemoryLayout, ColumnNorm, RowMaximum
import ArrayLayouts: sub_materialize, MemoryLayout, ColumnNorm, RowMaximum, CRowMaximum


struct MyMatrix <: LayoutMatrix{Float64}
Expand All @@ -10,6 +10,7 @@ Base.getindex(A::MyMatrix, k::Int, j::Int) = A.A[k,j]
Base.setindex!(A::MyMatrix, v, k::Int, j::Int) = setindex!(A.A, v, k, j)
Base.size(A::MyMatrix) = size(A.A)
Base.strides(A::MyMatrix) = strides(A.A)
Base.elsize(::Type{MyMatrix}) = sizeof(Float64)
Base.unsafe_convert(::Type{Ptr{T}}, A::MyMatrix) where T = Base.unsafe_convert(Ptr{T}, A.A)
MemoryLayout(::Type{MyMatrix}) = DenseColumnMajor()
Base.copy(A::MyMatrix) = MyMatrix(copy(A.A))
Expand All @@ -22,6 +23,7 @@ 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.elsize(::Type{MyVector}) = sizeof(Float64)
Base.unsafe_convert(::Type{Ptr{T}}, A::MyVector) where T = Base.unsafe_convert(Ptr{T}, A.A)
MemoryLayout(::Type{MyVector}) = DenseColumnMajor()

Expand Down Expand Up @@ -100,17 +102,14 @@ MemoryLayout(::Type{MyVector}) = DenseColumnMajor()

@test qr(A) isa LinearAlgebra.QRCompactWY
@test inv(A) ≈ inv(A.A)
if VERSION ≥ v"1.7-"
@test qr(A, Val(true)) == qr(A, ColumnNorm())
end

S = Symmetric(MyMatrix(reshape(inv.(1:25),5,5) + 10I))
@test cholesky(S).U ≈ @inferred(cholesky!(deepcopy(S))).U
@test cholesky(S,Val(true)).U ≈ cholesky(Matrix(S),Val(true)).U
@test cholesky(S,CRowMaximum()).U ≈ cholesky(Matrix(S),CRowMaximum()).U

S = Symmetric(MyMatrix(reshape(inv.(1:25),5,5) + 10I),:L)
@test cholesky(S).U ≈ @inferred(cholesky!(deepcopy(S))).U
@test cholesky(S,Val(true)).U ≈ cholesky(Matrix(S),Val(true)).U
@test cholesky(S,CRowMaximum()).U ≈ cholesky(Matrix(S),CRowMaximum()).U
end
Bin = randn(5,5)
B = MyMatrix(copy(Bin))
Expand Down
Loading