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

dlfivefifty · 2017-12-29T16:28:49Z

This pull request resolves Issues #17812 and #25247:

Deprecate stride(::AbstractArray, i::Int) and strides(::AbstractArray).
Re-implement above for the special cases of StridedArray.
Removes restrictions to StridedArray in BLAS routines, so that user types that are memory-equivalent to strided arrays can access these routines.
Removes the restrictions to StridedArray in sparse/sparsevector.jl. I'm not sure why these were there to begin.

I still need to weaken the types for ther LAPack routines. This is straightforward but will take some time to get the tests implemented.

stevengj · 2017-12-29T18:06:54Z

So you need a subtype of StridedArray to use BLAS? Wouldn't it be better to pull away from subtypes here and push towards traits?

dlfivefifty · 2017-12-29T18:09:05Z

That’s exactly what this is doing: the strides matrix interface consists of overriding conversion to Ptr and implementing stride.

Down the line a trait could be added as well but this is the first step.

stevengj · 2017-12-29T18:15:44Z

base/abstractarray.jl

-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


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

stevengj · 2017-12-29T18:19:11Z

doc/src/manual/interfaces.md

+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`


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

stevengj · 2017-12-29T20:32:51Z

base/abstractarray.jl

@@ -251,8 +251,8 @@ julia> strides(A)
 """
 function strides end

-# the definition of strides for Array is the cumsum of the product of sizes
-function _cumsumprodsizes(a::AbstractArray, i)
+# the definition of strides for Array is cumprod([1;collect(size(A))[1:end-1]])


well, technically cumprod would return an array, whereas strides returns a tuple

Hmm, that's not quite right since:

julia> A = Array{Float64,0}() 0-dimensional Array{Float64,0}: 0.0 julia> strides(A) ()

How about saying:

the defition of strides for Array{T,N} is tuple() if N = 0, otherwise it is a tuple containing 1 and a cumulative product of the first N-1 sizes

timholy

This is looking good, thanks for tackling this.

timholy · 2017-12-30T11:49:44Z

base/deprecated.jl

@@ -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)


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.

timholy · 2017-12-30T11:55:39Z

base/linalg/linalg.jl

@@ -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
+"""


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.

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).

timholy · 2017-12-30T11:58:03Z

base/subarray.jl

@@ -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)


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).

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.

timholy · 2017-12-30T12:02:29Z

doc/src/manual/interfaces.md

+A typical example of a user-defined strided array is one that wraps a standard `Array`
+with additional structure.
+
+


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

…tion message and docs

dlfivefifty · 2017-12-30T19:59:22Z

I've changed this a bit so that the default override is strides(A) and so one does not need to override stride(A,i). This is consistent with the abstract array interface, where one overrides size(A), not size(A,i).

I've also removed the restriction to StridedArray in the LAPACK routines, in the process I found and fixed some missing chkstride1 calls, and a bug in tzrzf! which used size instead of stride.

dlfivefifty · 2018-01-02T08:11:40Z

There is an issue with type-inferrence with the new substrides implementation:

julia> versioninfo()
Julia Version 0.7.0-DEV.3204
Commit d1086e0091* (2017-12-30 21:09 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin17.3.0)
  CPU: Intel(R) Core(TM) i7-7700HQ CPU @ 2.80GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell)
  LAPACK: libopenblas64_
  LIBM: libopenlibm
  LLVM: libLLVM-3.9.1 (ORCJIT, skylake)
Environment:
  JULIA_VERSION = 0.6

julia> A = rand(5,5,5,5);

julia> V = view(A, 1:1 ,:, 1:3, :);

julia> using Test

julia> @inferred strides(A)
(1, 5, 25, 125)

julia> @inferred strides(V)
ERROR: return type NTuple{4,Int64} does not match inferred return type Tuple{Int64,Int64,Int64,Any}
Stacktrace:
 [1] error(::String) at ./error.jl:33
 [2] top-level scope

This seems to be a limitation on type inference and I'm not sure how to resolve it.

The relevant function is

substrides(parent, I::Tuple) = 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]))"))

dlfivefifty · 2018-01-02T11:38:38Z

I've worked around the issue with substrides by adding an ::Int assertion.

I think this PR is ready to be merged, though I'm happy to make any other changes.

There is one test failure but it appears to be a download failure, which shouldn't be caused by this PR. It doesn't look like I have the ability to restart the test.

stevengj · 2018-01-03T15:10:31Z

An ::Int assertion may just be hiding the underlying problem; I'm worried that an Any variable is still heap-allocated somewhere and then converted to Int?

dlfivefifty · 2018-01-03T15:53:29Z

Good point. I'm not sure how to proceed as it appears to be a breakdown of type inference.

Perhaps putting the type assertion in the arguments is better? That is, for example:

substrides(parent, strds::NTuple{N,Int}, I::Tuple{Slice, Vararg{Any}}) where N = (first(strds), substrides(parent, tail(strds), tail(I))...)

Edit: "assertion" isn't the right word here.

timholy · 2018-01-04T15:38:56Z

Yeah, this does seem to be an inference bug. I get different answers depending on sequence of operations. For example:

@code_warntype substrides(A, (1, 5, 25, 125), (1:1, Base.Slice(Base.OneTo(5)), 1:3, Base.Slice(Base.OneTo(5))))  # "fails"
@code_warntype substrides(A, (5, 25, 125), (Base.Slice(Base.OneTo(5)), 1:3, Base.Slice(Base.OneTo(5))))  # "succeeds"
@code_warntype substrides(A, (1, 5, 25, 125), (1:1, Base.Slice(Base.OneTo(5)), 1:3, Base.Slice(Base.OneTo(5)))) # "succeeds"

timholy

Presuming tests pass, one tiny tweak and I think this is done. Many, many thanks for seeing this through!

timholy · 2018-01-04T19:22:38Z

test/linalg/blas.jl

+Base.setindex!(A::WrappedArray, v, i::Int) = setindex!(A.A, v, i)
+Base.setindex!(A::WrappedArray{T, N}, v, I::Vararg{Int, N}) where {T, N} = setindex!(A.A, v, I...)
+Base.unsafe_convert(::Type{Ptr{T}}, A::WrappedArray{T}) where T = Base.unsafe_convert(Ptr{T}, A.A)
+Base.stride(A::WrappedArray, i::Int) = stride(A.A, i)


Perhaps before this line you can test that strides(A::WrappedArray) and stride(A::WrappedArray, i) both throw an error, and that afterward it does not? In a sense, the main contribution of this PR is removing a fallback method, so we should test that the fallback doesn't accidentally get reintroduced.

Or you could define WrappedArray2 and test for the error.

I assume by "throws error" you mean @test_deprecated since the fallback is still there?

Right indeed. 👍

dlfivefifty · 2018-01-05T07:10:43Z

I think it’s ready to be merged now (again). The only failed test was a download error, which I don’t think this PR could cause.

mbauman · 2018-01-05T16:39:20Z

This is really nice! 👍

timholy · 2018-01-08T02:05:32Z

I belatedly realized that I should have asked you to add the deprecation to NEWS.md. Any chance you could do that now?

dlfivefifty · 2018-01-08T07:53:30Z

Yes I’ll try to do this soon.

Just like the deprecation of strides (#25321), the default implementation for `elsize` is generally not true. This deprecates the fallback, restructures it to be defined on the type (and not the value), and implements it where possible for the builtin array types and wrappers.

Nosferican · 2018-09-13T22:15:06Z

The deprecation method for this ain't very helpful on 0.7, given this message what should the appropriate method to update code be?

 Warning: The default `strides(a::AbstractArray)` implementation is deprecated for general arrays.
│ Specialize `strides(::Adjoint)` if `Adjoint` indeed uses a strided representation in memory.
│ Warning: inappropriately implementing this method for an array type that does not use strided
│ storage may lead to incorrect results or segfaults.
│   caller = stride at abstractarray.jl:350 [inlined]
└ @ Core ./abstractarray.jl:350

timholy · 2018-09-14T00:01:44Z

The alternative is to use cartesian indexing, i.e., A[i,j,k] rather than A[i] for a 3d array.

dlfivefifty · 2018-09-14T14:10:38Z

If your custom array is strided in memory, then it should implement the interface and override strides(::MyArray).

If it's not strided in memory, you probably actually mean to use size(A,1), not stride(A,2).

mbauman · 2018-09-14T15:21:34Z

If you're not using strides as a computation for pointer-loading (or passing to a C-API), then ideally you'd try to use eachindex or CartesianIndices to generate the loops you want. That requires a refactor, though. In the meantime, the operation strides is performing for Array is cumprod([1, size(A)...])[1:end-1]. It might be worth implementing cumprod(::Tuple) to make this easier and faster.

timholy · 2018-09-14T16:07:30Z

The xref above (added after my response) is from PyPlot, which probably means it is a pointer operation. Tread carefully in terms of following the advice above; I see there's an Adjoint in the message which means that the fallback implementation of cumprod([1, size(A)...])[1:end-1] may give the wrong answer.

The best advice might be to read up on strided arrays; once you understand the concept you will know if your array is strided and how to best implement the fix.

mbauman · 2018-09-14T16:34:28Z

It may be that #29135 will solve the issue. A second thumbs up there would be good and then it can be merged.

Sheehan Olver and others added 6 commits December 27, 2017 13:38

deprecate stride, support AbstractArray in blas routines

5c5d369

move size_to_strides back to abstractarray.jl

13eb90f

use parents stride in substrides

acaf1c2

Add strides for StridedReshapedArray

93998ad

Move stride, stride1 to dense.jl

0858576

move strides tests to test/dense.jl

903c42e

stevengj reviewed Dec 29, 2017

View reviewed changes

_cumsumprodsizes -> _defaultstride, addresses stevengl comments

4c2885d

stevengj reviewed Dec 29, 2017

View reviewed changes

Remove _defaultstride since it's only used one place

d353ab0

timholy reviewed Dec 30, 2017

View reviewed changes

dlfivefifty added 2 commits December 30, 2017 16:27

Remove restriction to StridedArray in LAPACK

31d4ce6

Make strides(A) the default override, make sterner warning in depreca…

9c8fb44

…tion message and docs

dlfivefifty changed the title ~~WIP: Deprecate generic stride implementation (#17812), allow AbstractArrays to use BLAS routines (#25247)~~ Deprecate generic stride implementation (#17812), allow AbstractArrays to use BLAS routines (#25247) Dec 30, 2017

dlfivefifty changed the title ~~Deprecate generic stride implementation (#17812), allow AbstractArrays to use BLAS routines (#25247)~~ Deprecate generic stride implementation (#17812), allow AbstractArrays to use BLAS/LAPACK routines (#25247) Dec 30, 2017

O(N) substrides implementation

d1086e0

ararslan added deprecation This change introduces or involves a deprecation linear algebra Linear algebra labels Dec 30, 2017

Deprecate old substrides

579ceec

work around type inference failure in substrides

337785d

state argument type in substrides so inference is correct

a558d76

timholy mentioned this pull request Jan 4, 2018

Allow transfer of array views backed by contiguous memory (fix #76) JuliaGPU/CUDAdrv.jl#77

Merged

timholy approved these changes Jan 4, 2018

View reviewed changes

dlfivefifty added 2 commits January 4, 2018 21:07

merge master

7ce3036

test_deprecated for strides before implementing

346605b

timholy merged commit 0ddf9ed into JuliaLang:master Jan 5, 2018

dlfivefifty mentioned this pull request Jan 5, 2018

RFC: Completing the StridedArray interface JuliaLang/LinearAlgebra.jl#186

Open

This was referenced Jan 11, 2018

Add strides and stride deprecation to NEWS #25510

Merged

Rewrite mul! to dispatch based on memory layout, not matrix type #25558

Closed

mbauman mentioned this pull request Mar 8, 2018

Deprecate fallback elsize for AbstractArray #26379

Merged

cormullion mentioned this pull request Mar 18, 2018

NEWS.md is getting a bit untidy #26508

Closed

Nosferican mentioned this pull request Sep 14, 2018

pcolor(): Warning: The default strides(a::AbstractArray) implementation is deprecated ....... JuliaPy/PyPlot.jl#380

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

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

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

dlfivefifty commented Dec 29, 2017

stevengj commented Dec 29, 2017

dlfivefifty commented Dec 29, 2017

stevengj Dec 29, 2017

stevengj Dec 29, 2017

stevengj Dec 29, 2017

dlfivefifty Dec 29, 2017

timholy left a comment

timholy Dec 30, 2017

timholy Dec 30, 2017

dlfivefifty Dec 30, 2017

timholy Dec 30, 2017

dlfivefifty Dec 30, 2017 •

edited

Loading

timholy Dec 30, 2017

dlfivefifty commented Dec 30, 2017

dlfivefifty commented Jan 2, 2018

dlfivefifty commented Jan 2, 2018

stevengj commented Jan 3, 2018

dlfivefifty commented Jan 3, 2018 •

edited

Loading

timholy commented Jan 4, 2018

timholy left a comment

timholy Jan 4, 2018

dlfivefifty Jan 4, 2018

timholy Jan 4, 2018

dlfivefifty commented Jan 5, 2018

mbauman commented Jan 5, 2018

timholy commented Jan 8, 2018

dlfivefifty commented Jan 8, 2018

Nosferican commented Sep 13, 2018

timholy commented Sep 14, 2018

dlfivefifty commented Sep 14, 2018

mbauman commented Sep 14, 2018

timholy commented Sep 14, 2018

mbauman commented Sep 14, 2018

		A typical example of a user-defined strided array is one that wraps a standard `Array`
		with additional structure.

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

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

Conversation

dlfivefifty commented Dec 29, 2017

stevengj commented Dec 29, 2017

dlfivefifty commented Dec 29, 2017

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

timholy left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dlfivefifty Dec 30, 2017 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dlfivefifty commented Dec 30, 2017

dlfivefifty commented Jan 2, 2018

dlfivefifty commented Jan 2, 2018

stevengj commented Jan 3, 2018

dlfivefifty commented Jan 3, 2018 • edited Loading

timholy commented Jan 4, 2018

timholy left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dlfivefifty commented Jan 5, 2018

mbauman commented Jan 5, 2018

timholy commented Jan 8, 2018

dlfivefifty commented Jan 8, 2018

Nosferican commented Sep 13, 2018

timholy commented Sep 14, 2018

dlfivefifty commented Sep 14, 2018

mbauman commented Sep 14, 2018

timholy commented Sep 14, 2018

mbauman commented Sep 14, 2018

dlfivefifty Dec 30, 2017 •

edited

Loading

dlfivefifty commented Jan 3, 2018 •

edited

Loading