Open Jutho opened 9 years ago
It's a tough call; on one hand it seems like "internal implementation detail" but on the other hand we export StridedArray
itself. I lean towards filling this interface out, but I guess the other approach would be to stop exporting StridedArray
and parent
. (Actually, given that many AbstractArrays wrap another array, parent
seems more fundamental, but OTOH it's not very meaningful without more details about the wrapper.)
Indeed, I could also see the other proposal work. Stop exporting StridedArray
, maybe even rename it to BlasArray
in the context of BLAS, just like BlasFloat
etc. Or probably just BlasVector
and BlasMatrix
. I guess it would even be possible to restrict BlasMatrix
such that it is automatically stride 1 in the first index, (e.g. union of DenseArray
and subarrays with UnitRange
indexing in the first index?). I don't know anything about the requirements of FFTW. But then, I think there is also no point of exporting strides
, or maybe it does, but it should be moved in the C interface section of the docs, and it should not be defined for general AbstractArray
.
@Jutho The value of strided array is not limited to BLAS. There are many C libraries out there that require strided or contiguous memory layout. Unexporting such types would discourage new packages that work with such C libraries.
I think we should clearly define the semantics of a strided array, and standardize the interface.
From a syntax standpoint, all strided array types should provide these functions (in addition to those that provided by all arrays): strides
and pointer
. From a semantics standard point, the underlying memory layout should conform to what these methods say.
I agree that the concept of strided arrays is very useful. Both for C libraries, but maybe also for pure julia methods. The main point about this issue is that the interface for accessing a strided array is not complete.
Either strided array is only for C libraries, then pointer
and strides
is sufficient as interface but it should belong in the C interop part of the standard library. If pure julia methods should also be allowed to use it, then the interface is not complete and too conflated with linear indexing.
Ref: https://github.com/JuliaLang/julia/pull/25321
This added documentation for a strided array interface consisting of overriding strides(A)
and Base.unsafe_convert(::Type{Ptr{T}}, A)
. I didn't need to use a trait to do this.
I may put together another pull request to deprecate StridedArray
by adding support for a trait that specifies how an arrays memory is layed out, allowing support for banded matrices and triangular matrices in addition to strided arrays. What I'm thinking is something like the following:
abstract type ArrayMemoryLayout{T} end
struct UnknownLayout{T} <: ArrayMemoryLayout{T} end
struct StridedLayout{T,trans} <: ArrayMemoryLayout{T} end
struct SymmetricStridedLayout{T} <: ArrayMemoryLayout{T} end
struct UpperTriangularLayout{T} <: ArrayMemoryLayout{T} end
memorylayout(::Type{A}) where A <: AbstractArray{T} where T = UnknownLayout{T}()
memorylayout(::Type{A}) where A <: Array{T} where T = StridedLayout{T,:N}()
memorylayout(A::AbstractArray) = memorylayout(typeof(A))
...
Then linear algebra routines will dispatch to BLAS and LAPACK when memorylayout
returns StridedLayout{<:BlasFloat}
, etc. Also, BandedMatrices.jl can add a BandedLayout
and SymBandedLayout
, and enable automatic dispatching to the appropriate LAPACK routines.
Note that this construction is not quite the same as a trait/interface as memory layout is always unique, whereas a single array can have multiple traits/interfaces.
I still need to think how to handle adjoints/transposes, but maybe something like this will work:
adjoint(::ArrayMemoryLayout{T}) where T = UnknownLayout{T}()
adjoint(::StridedLayout{T,:N}) where T = StridedLayout{T,:C}()
memorylayout(::Type{A}) where A <: Adjoint{T, AA} where {T, AA<:AbstractArray} = adjoint(memorylayout(AA))
Yes, the big question is what the kinds of properties are important enough to make it into the trait hierarchy and how that gets put together. In general we'd been thinking that this trait can layer onto 1.x as a new feature, but I'd most value a second set of eyes in making sure that things are indeed structured in such a way to support that.
My current thought is that this will require multiple traits that work together for different orthogonal dimensions of variance in how arrays get stored. Here's a start of a document on what I've been thinking: https://gist.github.com/mbauman/3fa8a1ac00f574480968a4f47be2c361
Let me reiterate my main point: "trait" isn't the right word here since memory storage is unique, hence we can ask a type what it is. So instead of querying isstrided(A)
(which doesn't scale well when you have isbanded(A)
, isuppertriangular(A)
, etc.) we can query memorylayout(A)
which tells us the memory layout.
I'm imagining that instead of the current code
mul!(y::StridedVector{T}, A::StridedVecOrMat{T}, x::StridedVector{T}) where {T<:BlasFloat} = gemv!(y, 'N', A, x)
mul!(y::StridedVector{T}, transA::Transpose{<:Any,<:StridedVecOrMat{T}}, x::StridedVector{T}) where {T<:BlasFloat} = (A = transA.parent; gemv!(y, 'T', A, x))
we would have
mul!(y, A, x) = _mul!(y, A, x, memoryformat(y), memoryformat(A), memoryformat(x))
_mul!(y::AbstractVector{T}, A::AbstractMatrix{T}, x::AbstractVector{T}, ::StridedFormat{T,:N}, ::StridedFormat{T,:N}, ::StridedFormat{T,:N}) where T<:BlasFloat = gemv!(y, 'N', A, x)
_mul!(y::AbstractVector{T}, A::AbstractMatrix{T}, x::AbstractVector{T}, ::StridedFormat{T,:N}, ::StridedFormat{T,:T}, ::StridedFormat{T,:N}) where T<:BlasFloat = gemv!(y, 'T', A, x)
_mul!(y::AbstractVector, A::AbstractMatrix, x::AbstractVector, ::UnknownFormat, ::StridedFormat, ::StridedFormat)= # native julia fallback
(We don't need the transA.parent
call in the proposed version since we would add pointer(A::Transpose) = pointer(A.parent)
and strides(A::Transpose) = strides(A.parent)
.)
Yup, you have roughly the same vision I do. I would just call memoryformat
— exactly as you've written it — a trait. I don't see the distinction you're trying to make. The fact that you can ask a type (and not just values) about a property like this is precisely what makes it an effective technique; it'd be type-unstable and slow otherwise. That's the core of a trait.
It’s only a point about semantics. “Traits” implies possibly multiple traits at the same time (a matrix can be banded and upper triangular at the same time) whereas the actual storage in memory is unique (the entries of a banded and upper triangular matrix could be stored either in eitherBandedLayout
Or UpperTriangularLayout
, but only one of these is valid for a given implementation)
Also mentioned in #10889: There might be some advantage to making it really a hierarchy instead of just a bunch of mutually exclusive types, e.g.
abstract type ArrayMemoryLayout{T} end
abstract type UnknownMemoryLayout{T} <: AbstractMemoryLayout{T} end
abstract type StridedFormat{T,...} <: AbstractMemoryLayout{T} end
abstract type UnitStridedFormat{T,...} <: StridedFormat{T} end
struct Contiguous{T,...} <: UnitStridedFormat{T} end
since some routines just require strided, other require strided with unit first stride, and e.g. linear indexing requires contiguous layout. The downside is that you need to use the type itself instead of the singleton instance as trait (unless you also make a concrete struct
for every level in the hierarchy), such that the method definitions would need to use arguments of the form ::Type{<:UnitStridedFormat}
.
Note added: I think that's also the contents of gist of @mbaumann, and so there should be an additional level between StridedFormat
and UnitStrided
which is MonotonousStrided
or something alike.
@dlfivefifty , thanks for all your nice work for completing the strided interface. I do have a question/request, after reading https://docs.julialang.org/en/latest/manual/interfaces/#man-interface-strided-arrays-1
It would be good if information about the strided memory layout is not only something to link to external C or Fortran libraries, but can also be used in Julia algorithms, i.e. to write native Julia algorithms that access memory in an optimized way. For plain Array
s, one would just use linear indexing using indices that you computed using the stride information. But what for other types that have only strides
and unsafe_convert(Ptr,..)
? For SubArray
s, there is parent
and first_index
to basically use linear indexing into the parent array. parent
actually also works on Array
and returns itself, first_index
does not.
Maybe parent
and first_index
can become part of the Strided interface, so that in a Julia algorithm, where parent
should always return a LinearIndexing
array.
Is your suggestion that parent(A::Array) = A
?
That is already the case in Base. But first_index(A::Array) = 1
does not exist.
It would be good if information about the strided memory layout is not only something to link to external C or Fortran libraries, but can also be used in Julia algorithms
It is in PR #25558: the order memory is accessed by generic_matvecmul!
in Julia-native code depends on whether A
is ColumnMajor
or RowMajor
.
I'm not sure parent
is actually a useful generic function. The only way you know what to do with the parent is if you also know what the wrapper was. It's used quite a bit across the package ecosystem, but in every single case I've examined, its argument is constrained to a specific implementation with the enclosing function.
first_index
is not exported, and I think should remain so.
@Jutho How do you use parent
and first_index
? Is it accomplishing something beyond what pointer
and stride[s]
can do?
I've always been a bit suspicious of the array parent
function.
Here’s an example: both A
and view(A’,:,:)’
have the same entries stored in the same order, but their parents are two different matrices.
Amusingly, I was looking for a solution to a broken test for one of my PRs and found an attempt to use parent
generically — by myself. The key was to recurse until you found a fixed point.
https://github.com/JuliaLang/julia/issues/15209#issuecomment-193971647
Even then, it's not sufficient for that task as some arrays have multiple "parents" — like Bidiagonal's two vectors or (in the context of aliasing detection) SubArray's indices.
My main use case is the algorithms TensorOperations.jl
, and a new package Strided.jl
that I am working on (both on https://github.com/Jutho). If you want to permutedims!
one multidimensional array into another (or any joint operation between two (or even three or more) multidimensional arrays with strides that are completely different and, importantly, non-monotously increasing) you want to write some kind of blocked algorithm to exploit memory locality (cache and all that), similar to transpose and matrix multiplication.
The easiest way to do this is to work with strides directly and compute some kind of linear index (memory offset) for each of the arrays involved. But how should you then use this index, in a Julia algorithm. Doing unsafe_load(pointer(A), i)
does not seem very julian, and, last time I checked,
@simd
)The fact that pointer(A)
exists, means that there always exists some underlying dense array, not all entries of which might be valid, but at least the once that you correctly compute using the strides should be. So maybe I should just use my own 'parent' function, let's call it, say, denseparent(A) = unsafe_wrap(Array, pointer(A), dims; own = false)
.
C.f. ArrayInterface.jl (and things that build on it like StrideArrays.jl)
Many methods in Julia base, especially those that interface with external C libraries like BLAS, LAPACK and FFTW are only defined for
StridedArray
, which currently is aUnion
type. It has been argued before by @timholy that properties likestrided
andcontiguous
should actually be a trait (a THT).However, my main question is about the interface for a
StridedArray
, which is currently incomplete and somewhat inconsistent. The manual describesStridedArray
as being an arrayA
such that, increasing indexk
by 1 is equivalent to linear indexing where the linear index is increased bystride(A,k)
. But that is of course not true, since linear indexing is expected to be such that the linear index increases byprod(size(A)[1:k-1])
. A simple example isA=rand(5,5);B=sub(A,1:3,1:3)
. Clearly,stride(B,2)=5
, but linear indexing withi+5*(j-1)
would not give me element(i,j)
, you need to index withi+3*(j-1)
for that. In that respect,stride(A,k)
is also defined for a generalAbstractArray
and just returnsprod(size(A)[1:k-1])
, which does not correspond to its description in the manual as "distance in memory".So what is the role of
StridedArray
in Julia. Is it only to interface with C libraries and to be used in combination withpointer
? It is certainly useful to also exploit the stridedness in pure Julia methods, see e.g.permutedims!
. Currently, this is dealt with by checking explicitly whether theStridedArray
is aSubArray
, and in that case, to extract the underlying parent array and the offset, and then apply linear indexing to that. I noticed thatparent
is now part of the exported method list, but there is no corresponding method to get theoffset
orfirst_index
. Maybe it would be good to have this method as well and then define these methods (stride
andstrides
,parent
,offset
?) as the interface for a strided array (as type or as trait).