Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Deprecate generic stride implementation (#17812), allow AbstractArrays to use BLAS/LAPACK routines (#25247) #25321

Merged
merged 16 commits into from
Jan 5, 2018
Merged

Deprecate generic stride implementation (#17812), allow AbstractArrays to use BLAS/LAPACK routines (#25247) #25321

Show file tree
Hide file tree
Changes from 6 commits
Commits
Show all changes
16 commits
Select commit Hold shift + click to select a range
5c5d369
deprecate stride, support AbstractArray in blas routines
Dec 27, 2017
13eb90f
move size_to_strides back to abstractarray.jl
dlfivefifty Dec 28, 2017
acaf1c2
use parents stride in substrides
dlfivefifty Dec 28, 2017
93998ad
Add strides for StridedReshapedArray
dlfivefifty Dec 28, 2017
0858576
Move stride, stride1 to dense.jl
dlfivefifty Dec 29, 2017
903c42e
move strides tests to test/dense.jl
dlfivefifty Dec 29, 2017
4c2885d
_cumsumprodsizes -> _defaultstride, addresses stevengl comments
dlfivefifty Dec 29, 2017
d353ab0
Remove _defaultstride since it's only used one place
dlfivefifty Dec 30, 2017
31d4ce6
Remove restriction to StridedArray in LAPACK
dlfivefifty Dec 30, 2017
9c8fb44
Make strides(A) the default override, make sterner warning in depreca…
dlfivefifty Dec 30, 2017
d1086e0
O(N) substrides implementation
dlfivefifty Dec 30, 2017
579ceec
Deprecate old substrides
dlfivefifty Dec 31, 2017
337785d
work around type inference failure in substrides
dlfivefifty Jan 2, 2018
a558d76
state argument type in substrides so inference is correct
dlfivefifty Jan 4, 2018
7ce3036
merge master
dlfivefifty Jan 4, 2018
346605b
test_deprecated for strides before implementing
dlfivefifty Jan 4, 2018
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
28 changes: 16 additions & 12 deletions base/abstractarray.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -234,16 +234,7 @@ julia> stride(A,3)
12
```
"""
function stride(a::AbstractArray, i::Integer)
if i > ndims(a)
return length(a)
end
s = 1
for n = 1:(i-1)
s *= size(a, n)
end
return s
end
function stride end

"""
strides(A)
Expand All @@ -258,7 +249,20 @@ julia> strides(A)
(1, 3, 12)
```
"""
strides(A::AbstractArray) = size_to_strides(1, size(A)...)
function strides end

# the definition of strides for Array is the cumsum of the product of sizes
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

cumprod of the size? But it's not quite that, it's cumprod([1;collect(size(A))[1:end-1]])

function _cumsumprodsizes(a::AbstractArray, i)
if i > ndims(a)
return length(a)
end
s = 1
for n = 1:(i-1)
s *= size(a, n)
end
return s
end

@inline size_to_strides(s, d, sz...) = (s, size_to_strides(s * d, sz...)...)
size_to_strides(s, d) = (s,)
size_to_strides(s) = ()
Expand Down Expand Up @@ -2043,4 +2047,4 @@ function hash(a::AbstractArray{T}, h::UInt) where T
end
!isequal(x2, x1) && (h = hash(x2, h))
return h
end
end
14 changes: 14 additions & 0 deletions base/deprecated.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3489,6 +3489,20 @@ end
@deprecate getq(F::Factorization) F.Q
end

# Issues #17812 Remove default stride implementation
function stride(a::AbstractArray, i::Integer)
depwarn("`stride(a::AbstractArray, i)` is deprecated for general arrays. Override `stride` for custom array types that implement the strided array interface.", :stride)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I wonder if this is insufficiently cautionary. Maybe

Specialize `stride` for custom array types that have the appropriate representation in memory. Warning: inappropriately implementing this method for an array type that does not use strided storage can lead to incorrect results or segfaults.

if i > ndims(a)
return length(a)
end
s = 1
for n = 1:(i-1)
s *= size(a, n)
end
return s
_stride(a, i)
end

# END 0.7 deprecations

# BEGIN 1.0 deprecations
Expand Down
181 changes: 105 additions & 76 deletions base/linalg/blas.jl
View file
Open in desktop

Large diffs are not rendered by default.

31 changes: 2 additions & 29 deletions base/linalg/dense.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -112,35 +112,8 @@ true
isposdef(A::AbstractMatrix) = ishermitian(A) && isposdef(cholfact(Hermitian(A)))
isposdef(x::Number) = imag(x)==0 && real(x) > 0

"""
stride1(A) -> Int

Return the distance between successive array elements
in dimension 1 in units of element size.

# Examples
```jldoctest
julia> A = [1,2,3,4]
4-element Array{Int64,1}:
1
2
3
4

julia> Base.LinAlg.stride1(A)
1

julia> B = view(A, 2:2:4)
2-element view(::Array{Int64,1}, 2:2:4) with eltype Int64:
2
4

julia> Base.LinAlg.stride1(B)
2
```
"""
stride1(x::Array) = 1
stride1(x::StridedVector) = stride(x, 1)::Int
stride(A::Union{DenseArray,StridedReshapedArray,StridedReinterpretArray}, i::Int) = Base._cumsumprodsizes(A, i)
strides(A::Union{DenseArray,StridedReshapedArray,StridedReinterpretArray}) = size_to_strides(1, size(A)...)

function norm(x::StridedVector{T}, rx::Union{UnitRange{TI},AbstractRange{TI}}) where {T<:BlasFloat,TI<:Integer}
if minimum(rx) < 1 || maximum(rx) > length(x)
Expand Down
37 changes: 34 additions & 3 deletions base/linalg/linalg.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,8 +36,8 @@ import Base: USE_BLAS64, abs, acos, acosh, acot, acoth, acsc, acsch, adjoint, as
cosh, cot, coth, csc, csch, eltype, exp, findmax, findmin, fill!, floor, getindex, hcat,
getproperty, imag, inv, isapprox, isone, IndexStyle, kron, length, log, map, ndims,
oneunit, parent, power_by_squaring, print_matrix, promote_rule, real, round, sec, sech,
setindex!, show, similar, sin, sincos, sinh, size, sqrt, tan, tanh, transpose, trunc,
typed_hcat, vec
setindex!, show, similar, sin, sincos, sinh, size, size_to_strides, sqrt, StridedReinterpretArray,
StridedReshapedArray, strides, stride, tan, tanh, transpose, trunc, typed_hcat, vec
using Base: hvcat_fill, iszero, IndexLinear, _length, promote_op, promote_typeof,
@propagate_inbounds, @pure, reduce, typed_vcat
# We use `_length` because of non-1 indices; releases after julia 0.5
Expand Down Expand Up @@ -216,9 +216,40 @@ end
# Check that stride of matrix/vector is 1
# Writing like this to avoid splatting penalty when called with multiple arguments,
# see PR 16416
"""
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

With the new IPO constant propagation do we need stride1? It seems that stride(x, dim) could check for dim == 1, and the branch may be elided now.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This happens automatically at compile time? Cool!

I'm not sure how to test for performance regression or that it's doing the right behaviour, so I might hold off on that change (since it's independent of this PR).

stride1(A) -> Int

Return the distance between successive array elements
in dimension 1 in units of element size.

# Examples
```jldoctest
julia> A = [1,2,3,4]
4-element Array{Int64,1}:
1
2
3
4

julia> Base.LinAlg.stride1(A)
1

julia> B = view(A, 2:2:4)
2-element view(::Array{Int64,1}, 2:2:4) with eltype Int64:
2
4

julia> Base.LinAlg.stride1(B)
2
```
"""
stride1(x) = stride(x,1)
stride1(x::Array) = 1
stride1(x::DenseArray) = stride(x, 1)::Int

@inline chkstride1(A...) = _chkstride1(true, A...)
@noinline _chkstride1(ok::Bool) = ok || error("matrix does not have contiguous columns")
@inline _chkstride1(ok::Bool, A, B...) = _chkstride1(ok & (stride(A, 1) == 1), B...)
@inline _chkstride1(ok::Bool, A, B...) = _chkstride1(ok & (stride1(A) == 1), B...)

"""
LinAlg.checksquare(A)
Expand Down
52 changes: 26 additions & 26 deletions base/sparse/sparsevector.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1354,8 +1354,8 @@ for (vop, fun, mode) in [(:_vadd, :+, 1),
(:_vmul, :*, 0)]
@eval begin
$(vop)(x::AbstractSparseVector, y::AbstractSparseVector) = _binarymap($(fun), x, y, $mode)
$(vop)(x::StridedVector, y::AbstractSparseVector) = _binarymap($(fun), x, y, $mode)
$(vop)(x::AbstractSparseVector, y::StridedVector) = _binarymap($(fun), x, y, $mode)
$(vop)(x::AbstractVector, y::AbstractSparseVector) = _binarymap($(fun), x, y, $mode)
$(vop)(x::AbstractSparseVector, y::AbstractVector) = _binarymap($(fun), x, y, $mode)
end
end

Expand All @@ -1368,31 +1368,31 @@ broadcast(::typeof(*), x::AbstractSparseVector{Bool}, y::BitVector) = _vmul(x, y
for (op, vop) in [(:+, :_vadd), (:-, :_vsub), (:*, :_vmul)]
op != :* && @eval begin
$(op)(x::AbstractSparseVector, y::AbstractSparseVector) = $(vop)(x, y)
$(op)(x::StridedVector, y::AbstractSparseVector) = $(vop)(x, y)
$(op)(x::AbstractSparseVector, y::StridedVector) = $(vop)(x, y)
$(op)(x::AbstractVector, y::AbstractSparseVector) = $(vop)(x, y)
$(op)(x::AbstractSparseVector, y::AbstractVector) = $(vop)(x, y)
end
@eval begin
broadcast(::typeof($op), x::AbstractSparseVector, y::AbstractSparseVector) = $(vop)(x, y)
broadcast(::typeof($op), x::StridedVector, y::AbstractSparseVector) = $(vop)(x, y)
broadcast(::typeof($op), x::AbstractSparseVector, y::StridedVector) = $(vop)(x, y)
broadcast(::typeof($op), x::AbstractVector, y::AbstractSparseVector) = $(vop)(x, y)
broadcast(::typeof($op), x::AbstractSparseVector, y::AbstractVector) = $(vop)(x, y)
end
end

# definition of other binary functions

broadcast(::typeof(min), x::SparseVector{<:Real}, y::SparseVector{<:Real}) = _binarymap(min, x, y, 2)
broadcast(::typeof(min), x::AbstractSparseVector{<:Real}, y::AbstractSparseVector{<:Real}) = _binarymap(min, x, y, 2)
broadcast(::typeof(min), x::StridedVector{<:Real}, y::AbstractSparseVector{<:Real}) = _binarymap(min, x, y, 2)
broadcast(::typeof(min), x::AbstractSparseVector{<:Real}, y::StridedVector{<:Real}) = _binarymap(min, x, y, 2)
broadcast(::typeof(min), x::AbstractVector{<:Real}, y::AbstractSparseVector{<:Real}) = _binarymap(min, x, y, 2)
broadcast(::typeof(min), x::AbstractSparseVector{<:Real}, y::AbstractVector{<:Real}) = _binarymap(min, x, y, 2)

broadcast(::typeof(max), x::SparseVector{<:Real}, y::SparseVector{<:Real}) = _binarymap(max, x, y, 2)
broadcast(::typeof(max), x::AbstractSparseVector{<:Real}, y::AbstractSparseVector{<:Real}) = _binarymap(max, x, y, 2)
broadcast(::typeof(max), x::StridedVector{<:Real}, y::AbstractSparseVector{<:Real}) = _binarymap(max, x, y, 2)
broadcast(::typeof(max), x::AbstractSparseVector{<:Real}, y::StridedVector{<:Real}) = _binarymap(max, x, y, 2)
broadcast(::typeof(max), x::AbstractVector{<:Real}, y::AbstractSparseVector{<:Real}) = _binarymap(max, x, y, 2)
broadcast(::typeof(max), x::AbstractSparseVector{<:Real}, y::AbstractVector{<:Real}) = _binarymap(max, x, y, 2)

complex(x::AbstractSparseVector{<:Real}, y::AbstractSparseVector{<:Real}) = _binarymap(complex, x, y, 1)
complex(x::StridedVector{<:Real}, y::AbstractSparseVector{<:Real}) = _binarymap(complex, x, y, 1)
complex(x::AbstractSparseVector{<:Real}, y::StridedVector{<:Real}) = _binarymap(complex, x, y, 1)
complex(x::AbstractVector{<:Real}, y::AbstractSparseVector{<:Real}) = _binarymap(complex, x, y, 1)
complex(x::AbstractSparseVector{<:Real}, y::AbstractVector{<:Real}) = _binarymap(complex, x, y, 1)

### Reduction

Expand Down Expand Up @@ -1438,7 +1438,7 @@ adjoint(sv::SparseVector) = Adjoint(sv)

# axpy

function LinAlg.axpy!(a::Number, x::SparseVectorUnion, y::StridedVector)
function LinAlg.axpy!(a::Number, x::SparseVectorUnion, y::AbstractVector)
length(x) == length(y) || throw(DimensionMismatch())
nzind = nonzeroinds(x)
nzval = nonzeros(x)
Expand Down Expand Up @@ -1479,7 +1479,7 @@ scale!(a::Complex, x::SparseVectorUnion) = (scale!(nonzeros(x), a); x)
(/)(x::SparseVectorUnion, a::Number) = SparseVector(length(x), copy(nonzeroinds(x)), nonzeros(x) / a)

# dot
function dot(x::StridedVector{Tx}, y::SparseVectorUnion{Ty}) where {Tx<:Number,Ty<:Number}
function dot(x::AbstractVector{Tx}, y::SparseVectorUnion{Ty}) where {Tx<:Number,Ty<:Number}
n = length(x)
length(y) == n || throw(DimensionMismatch())
nzind = nonzeroinds(y)
Expand Down Expand Up @@ -1543,7 +1543,7 @@ end
### BLAS-2 / dense A * sparse x -> dense y

# lowrankupdate (BLAS.ger! like)
function LinAlg.lowrankupdate!(A::StridedMatrix, x::StridedVector, y::SparseVectorUnion, α::Number = 1)
function LinAlg.lowrankupdate!(A::StridedMatrix, x::AbstractVector, y::SparseVectorUnion, α::Number = 1)
nzi = nonzeroinds(y)
nzv = nonzeros(y)
@inbounds for (j,v) in zip(nzi,nzv)
Expand All @@ -1565,10 +1565,10 @@ function (*)(A::StridedMatrix{Ta}, x::AbstractSparseVector{Tx}) where {Ta,Tx}
mul!(y, A, x)
end

mul!(y::StridedVector{Ty}, A::StridedMatrix, x::AbstractSparseVector{Tx}) where {Tx,Ty} =
mul!(y::AbstractVector{Ty}, A::StridedMatrix, x::AbstractSparseVector{Tx}) where {Tx,Ty} =
mul!(one(Tx), A, x, zero(Ty), y)

function mul!(α::Number, A::StridedMatrix, x::AbstractSparseVector, β::Number, y::StridedVector)
function mul!(α::Number, A::StridedMatrix, x::AbstractSparseVector, β::Number, y::AbstractVector)
m, n = size(A)
length(x) == n && length(y) == m || throw(DimensionMismatch())
m == 0 && return y
Expand Down Expand Up @@ -1603,10 +1603,10 @@ function *(transA::Transpose{<:Any,<:StridedMatrix{Ta}}, x::AbstractSparseVector
mul!(y, Transpose(A), x)
end

mul!(y::StridedVector{Ty}, transA::Transpose{<:Any,<:StridedMatrix}, x::AbstractSparseVector{Tx}) where {Tx,Ty} =
mul!(y::AbstractVector{Ty}, transA::Transpose{<:Any,<:StridedMatrix}, x::AbstractSparseVector{Tx}) where {Tx,Ty} =
(A = transA.parent; mul!(one(Tx), Transpose(A), x, zero(Ty), y))

function mul!(α::Number, transA::Transpose{<:Any,<:StridedMatrix}, x::AbstractSparseVector, β::Number, y::StridedVector)
function mul!(α::Number, transA::Transpose{<:Any,<:StridedMatrix}, x::AbstractSparseVector, β::Number, y::AbstractVector)
A = transA.parent
m, n = size(A)
length(x) == m && length(y) == n || throw(DimensionMismatch())
Expand Down Expand Up @@ -1662,10 +1662,10 @@ end

# * and mul!

mul!(y::StridedVector{Ty}, A::SparseMatrixCSC, x::AbstractSparseVector{Tx}) where {Tx,Ty} =
mul!(y::AbstractVector{Ty}, A::SparseMatrixCSC, x::AbstractSparseVector{Tx}) where {Tx,Ty} =
mul!(one(Tx), A, x, zero(Ty), y)

function mul!(α::Number, A::SparseMatrixCSC, x::AbstractSparseVector, β::Number, y::StridedVector)
function mul!(α::Number, A::SparseMatrixCSC, x::AbstractSparseVector, β::Number, y::AbstractVector)
m, n = size(A)
length(x) == n && length(y) == m || throw(DimensionMismatch())
m == 0 && return y
Expand Down Expand Up @@ -1695,21 +1695,21 @@ end

# * and *(Tranpose(A), B)

mul!(y::StridedVector{Ty}, transA::Transpose{<:Any,<:SparseMatrixCSC}, x::AbstractSparseVector{Tx}) where {Tx,Ty} =
mul!(y::AbstractVector{Ty}, transA::Transpose{<:Any,<:SparseMatrixCSC}, x::AbstractSparseVector{Tx}) where {Tx,Ty} =
(A = transA.parent; mul!(one(Tx), Transpose(A), x, zero(Ty), y))

mul!(α::Number, transA::Transpose{<:Any,<:SparseMatrixCSC}, x::AbstractSparseVector, β::Number, y::StridedVector) =
mul!(α::Number, transA::Transpose{<:Any,<:SparseMatrixCSC}, x::AbstractSparseVector, β::Number, y::AbstractVector) =
(A = transA.parent; _At_or_Ac_mul_B!(*, α, A, x, β, y))

mul!(y::StridedVector{Ty}, adjA::Adjoint{<:Any,<:SparseMatrixCSC}, x::AbstractSparseVector{Tx}) where {Tx,Ty} =
mul!(y::AbstractVector{Ty}, adjA::Adjoint{<:Any,<:SparseMatrixCSC}, x::AbstractSparseVector{Tx}) where {Tx,Ty} =
(A = adjA.parent; mul!(one(Tx), Adjoint(A), x, zero(Ty), y))

mul!(α::Number, adjA::Adjoint{<:Any,<:SparseMatrixCSC}, x::AbstractSparseVector, β::Number, y::StridedVector) =
mul!(α::Number, adjA::Adjoint{<:Any,<:SparseMatrixCSC}, x::AbstractSparseVector, β::Number, y::AbstractVector) =
(A = adjA.parent; _At_or_Ac_mul_B!(dot, α, A, x, β, y))

function _At_or_Ac_mul_B!(tfun::Function,
α::Number, A::SparseMatrixCSC, x::AbstractSparseVector,
β::Number, y::StridedVector)
β::Number, y::AbstractVector)
m, n = size(A)
length(x) == m && length(y) == n || throw(DimensionMismatch())
n == 0 && return y
Expand Down
8 changes: 4 additions & 4 deletions base/subarray.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -255,11 +255,11 @@ IndexStyle(::Type{<:SubArray}) = IndexCartesian()
# so they are well-defined even for non-linear memory layouts
strides(V::SubArray) = substrides(V.parent, V.indices)

substrides(parent, I::Tuple) = substrides(1, parent, 1, I)
substrides(parent, I::Tuple) = substrides(stride(parent, 1), parent, 1, I)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I imagine that many implementations of stride(x, dim) will involve a loop over d = 1:dim; since you call stride once per dimension, this is now O(N^2) where N = ndims(parent). It might be more efficient to write this in terms of strides(parent) which can be implemented with a single loop and thus be O(N).

Copy link
Contributor Author

@dlfivefifty dlfivefifty Dec 30, 2017 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good point. Maybe something like:

substrides(parent, I) = substrides(parent, strides(parent), I)
substrides(parent, strds, ::Tuple{}) = ()
substrides(parent, strds, I::Tuple{ScalarIndex, Vararg{Any}}) = (substrides(parent, tail(strds), tail(I))...,)
substrides(parent, strds, I::Tuple{Slice, Vararg{Any}}) = (first(strds), substrides(parent, tail(strds), tail(I))...)
substrides(parent, strds, I::Tuple{AbstractRange, Vararg{Any}}) = (first(strds)*step(I[1]), substrides(parent, tail(strds), tail(I))...)
substrides(parent, strds, I::Tuple{Any, Vararg{Any}}) = throw(ArgumentError("strides is invalid for SubArrays with indices of type $(typeof(I[1]))"))

Does substrides have to go through a deprecation process? It's not exported.

substrides(s, parent, dim, ::Tuple{}) = ()
substrides(s, parent, dim, I::Tuple{ScalarIndex, Vararg{Any}}) = (substrides(s*size(parent, dim), parent, dim+1, tail(I))...,)
substrides(s, parent, dim, I::Tuple{Slice, Vararg{Any}}) = (s, substrides(s*size(parent, dim), parent, dim+1, tail(I))...)
substrides(s, parent, dim, I::Tuple{AbstractRange, Vararg{Any}}) = (s*step(I[1]), substrides(s*size(parent, dim), parent, dim+1, tail(I))...)
substrides(s, parent, dim, I::Tuple{ScalarIndex, Vararg{Any}}) = (substrides(stride(parent, dim+1), parent, dim+1, tail(I))...,)
substrides(s, parent, dim, I::Tuple{Slice, Vararg{Any}}) = (s, substrides(stride(parent, dim+1), parent, dim+1, tail(I))...)
substrides(s, parent, dim, I::Tuple{AbstractRange, Vararg{Any}}) = (s*step(I[1]), substrides(stride(parent, dim+1), parent, dim+1, tail(I))...)
substrides(s, parent, dim, I::Tuple{Any, Vararg{Any}}) = throw(ArgumentError("strides is invalid for SubArrays with indices of type $(typeof(I[1]))"))

stride(V::SubArray, d::Integer) = d <= ndims(V) ? strides(V)[d] : strides(V)[end] * size(V)[end]
Expand Down
18 changes: 18 additions & 0 deletions doc/src/manual/interfaces.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -391,6 +391,24 @@ something other than 1), you should specialize `indices`. You should also specia
so that the `dims` argument (ordinarily a `Dims` size-tuple) can accept `AbstractUnitRange` objects,
perhaps range-types `Ind` of your own design. For more information, see [Arrays with custom indices](@ref).

## [Strided Arrays]

| Methods to implement |   | Brief description |
|:----------------------------------------------- |:-------------------------------------- |:------------------------------------------------------------------------------------- |
| `stride(A, i::Int)` |   | Return the distance in memory (in number of elements) between adjacent elements in dimension k. |
| `Base.unsafe_convert(::Type{Ptr{T}}, A)` |   | Return the native address of an array. |


A strided array is a sub-type of `AbstractArray` whose entries are stored in memory with fixed strides.
Provided the element type of the array is compatible with BLAS, a strided array can utilize BLAS and LAPACK routines
for more efficient linear algebra routines.

A typical example of a user defined strided array is one that wraps a standard `Array`
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

hyphenate "user-defined" (hyphenate compound adjectives). Don't hyphenate "subtype".

with additional structure.


Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'd add a warning here too about how dangerous it is to implement these methods if the underlying storage is not actually strided. Otherwise the following might happen:

  • user calls a function. Gets a MethodError due to missing implementation of strides
  • "helpfully" implements strides to fix the method error

It might be useful to give some concrete examples of arrays and categorize them as to whether they have strided memory representation. Examples:

1:5   # not strided (there is no storage associated with this array, could segfault if implemented `strides`)
collect(1:5)  # is strided
A = [1 5; 2 6; 3 7; 4 8]  # is strided
V = view(A, 1:2, :)   # is strided
V = view(A, 1:2:4, :)   # is strided
V = view(A, [1,2,4], :)   # is not strided



## [Broadcasting](@id man-interfaces-broadcasting)

| Methods to implement | Brief description |
Expand Down
7 changes: 0 additions & 7 deletions test/abstractarray.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -431,13 +431,6 @@ function test_primitives(::Type{T}, shape, ::Type{TestAbstractArray}) where T
# last(a)
@test last(B) == B[length(B)]

# strides(a::AbstractArray)
@inferred strides(B)
strides_B = strides(B)
for (i, _stride) in enumerate(collect(strides_B))
@test _stride == stride(B, i)
end

# isassigned(a::AbstractArray, i::Int...)
j = rand(1:length(B))
@test isassigned(B, j) == true
Expand Down
Loading