Fortran-like arrays such as FArray(Float64, -1:1,-7:7,-128:512)

[Alexander Samoilov]

Dear developers,

I'm aware about several discussions in this group and julia-users for zero-based Julia arrays, proposals to have negative indices for Julia arrays, like Fortran has, and so on.

Generally having 1-based array for Julia is a good decision, but sometimes it is desirable to have Fortran-like array with indices that span some subranges of ℤ:

    julia> x = FArray(Float64, -1:1,-7:7,-128:512)

    FArray((-1:1,-7:7,-128:512),3x15x641 Array{Float64,3}:

    [:, :, 1] = ...

where it would be useful:

in the codes accompanying the book Numerical Solution of Hyperbolic Partial Differential Equations by prof. John A. Trangenstein these negative indices are used intensively for ghost cells for boundary conditions.

The same is true for Clawpack (stands for “Conservation Laws Package”) by prof. Randall J. LeVeque http://depts.washington.edu/clawpack/ and there are many other codes where such indices would be natural.

So such auxiliary class would be useful for speedy translation of such codes.

I just started to implement such an auxiliary type but as I'm quite new to Julia your help would be greatly appreciated.

I started with:

type FArray

    ranges

    array::Array

    function FArray(T, r::Range1{Int}...)

        dims = map((x) -> length(x), r)

        array = Array(T, dims)

        new(r, array)

    end

end

and the Julia REPL output is:

> julia> include ("FortranArray.jl")

> 

> julia> x = FArray(Float64, -1:1,-7:7,-128:512)

> FArray((-1:1,-7:7,-128:512),3x15x641 Array{Float64,3}:

> [:, :, 1] =

>  6.90321e-310  2.6821e-316   1.96042e-316  0.0  0.0  0.0  9.84474e-317  …  1.83233e-316  2.63285e-316  0.0           9.61618e-317  0.0        

>  6.90321e-310  6.32404e-322  2.63285e-316  0.0  0.0  0.0  2.63292e-316     2.67975e-316  0.0           0.0           0.0           0.0        

>  2.6821e-316   6.32404e-322  0.0           0.0  0.0  0.0  2.3428e-316      1.87441e-316  0.0           2.12342e-316  0.0           2.1135e-316

Questions:

1. I'd like to add type signature to ranges, say ranges::Vector{Range1{Int}} instead of typeless ranges (is it more efficient than tuple for index computations?) but can't convert from tuple to vector (in general it is impossible as tuple can store different types).
2. I'd like to specialize array more, don't know how
3. getindex, setindex are not written yet, any recommendations would be appreciated.
4. Any other comments, proposals would be greatly appreciated.

Thanks!

Alexander

[Tim Holy]

I don't know if it makes your FArray unnecessary, but you can already
construct a SubArray with such unconventional indices. In fact, all of the
questions you asked may be answered by a careful study of base/subarray.jl;
it's a pretty good example of the type system and of how you write access
functions for a new array type.

Best,
--Tim

> 1. I'd like to add type signature to ranges, say
> *ranges::Vector{Range1{Int}} *instead of typeless *ranges *(is it more

> efficient than tuple for index computations?) but can't convert from tuple
> to vector (in general it is impossible as tuple can store different types).

> 2. I'd like to specialize array more, don't know how
> 3. getindex, setindex are not written yet, any recommendations would be
> appreciated.
> 4. Any other comments, proposals would be greatly appreciated.
>
> Thanks!
> Alexander

[Kevin Squire]

Just as a concrete example (mostly because I hadn't explore this and find it cool):

julia> a = [-10:10];

julia> b = sub(a, 12:21)

10-element SubArray{Int64,1,Array{Int64,1},(Range1{Int64},)}:

  1

  2

  3

  4

  5

  6

  7

  8

  9

 10

julia> b[1]

1

julia> b[0]

0

julia> b[-10]

-10

julia> b[10]

10

julia> b[-11]

ERROR: BoundsError()

 in getindex at subarray.jl:176

julia> b[11]

ERROR: BoundsError()

 in getindex at subarray.jl:176

Cheers!  Kevin

[Alexander Samoilov]

Thanks a lot for the hint!

Got to the following convenience function draft:

julia> function farray(T, r::Range1{Int}...)

       dims = map((x) -> length(x), r)

       array = Array(T, dims)

       sub_indices = map((x) -> -minimum(x) + 2 : maximum(x), r)

       sub(array, sub_indices)

       end

farray (generic function with 1 method)

julia> x = farray(Float64, -1:1,-7:7,-128:512)

0x0x383 SubArray{Float64,3,Array{Float64,3},(Range1{Int64},Range1{Int64},Range1{Int64})}

julia> x[-1,-7,-128]

6.9202491077408e-310

julia> x[-2,-7,-128]

ERROR: BoundsError()

 in getindex at subarray.jl:182

julia> x[-1,-8,-128]

ERROR: BoundsError()

 in getindex at subarray.jl:182

Guess a serious performance penalty would be incurred due to indirect indexing.

Cheers,

Alexander

[Alexander Samoilov]

Unfortunately indeed there is a serious performance penalty for such an attempt of using sub as an replacement for negative/zero based arrays:

Original Fortran code:

time ./main_f 

real 0m0.624s

user 0m0.616s

sys 0m0.004s

Translated to Julia with lower bounds of the array shifted:

julia main.jl 

elapsed time: 0.440328351 seconds

Ha, not bad :-)

julia main_farr.jl 

elapsed time: 44.954791721 seconds

Performance degraded by 100 times - 2 orders of magnitude :-(

The original Fortran code uses 3 arrays

      double precision

     &  u(-2:ncells+1),

     &  x(0:ncells),

     &  flux(0:ncells)

in main.jl the arrays are just translated to:

      u    = Array(Float64, ncells+4)

      x    = Array(Float64, ncells+1)

      flux = Array(Float64, ncells+1)

and in the all index expressions an appropriate shift was added.

and in the new code where sub is used:

function farray(T, r::Range1{Int}...)

    dims = map((x) -> length(x), r)

    array = Array(T, dims)

    sub_indices = map((x) -> -minimum(x) + 2 : maximum(x), r)

    sub(array, sub_indices)

end

function main()

[...]

      u    = farray(Float64, -2:ncells+1)

      x    = farray(Float64,  0:ncells)

      flux = farray(Float64,  0:ncells)

all index expression were reverted to original Fortran ones.

Thanks,

Alexander

[Tim Holy]

There are well-known performance limitations with SubArrays, but they don't
occur everywhere and they may not be due to the things you're expecting. It's
very unlikely that simply adding a constant shift to the indexing takes so
long. Most _optimized_ operations on arrays tend to be cache-limited, and an
extra add or multiply here and there tends to make almost no difference.
Fundamentally it comes down to the speed of light; operations within the CPU
are fast, and operations that need to go fetch something from another chip on
your motherboard are slow.

Unfortunately, a poorly-optimized routine can be a different story. For
SubArrays, linear indexing is currently really slow (search the issues for
more detail). Since many core algorithms in Julia, like sum(), use linear
indexing, this can be a problem. But using carestian indexing, e.g., `S[3, 5]`
for a 2-d subarray `S`, should be pretty good (max 2x speed reduction).

I didn't see a full definition of main_farr in your email (I may have
overlooked it or not understood what I was seeing), so it's hard to say where
your 100x slowdown occurs. You might try the profiler to learn more about
what's happening.

--Tim

> > julia> x[*-2*,-7,-128]

> > ERROR: BoundsError()
> >
> > in getindex at subarray.jl:182
> >

> > julia> x[-1,*-8*,-128]

[Alexander Samoilov]

Hello Tim,

Fully agree with you here - surprised to see such huge performance degradation.

The complete definition is about 140 lines - too lengthy for the letter.

I will put all the 3 codes (fortran, Array based and SubArray based) somewhere and prepare a shorter testcase.

I profiled the codes:

For the fast Array based code I got:

 cat main.jl.profile 

elapsed time: 0.290131068 seconds

elapsed time: 0.305984925 seconds

288 boot.jl; include; line: 238

    288 profile.jl; anonymous; line: 14

     4   ...law/PROGRAM0/main.jl; main; line: 101

     46  ...law/PROGRAM0/main.jl; main; line: 108

     55  ...law/PROGRAM0/main.jl; main; line: 109

     183 ...law/PROGRAM0/main.jl; main; line: 115

For slow SubArray based code:

cat main_farr.jl.profile 

elapsed time: 31.571042152 seconds

elapsed time: 32.855328249 seconds

30239 boot.jl; include; line: 238

    30239 profile.jl; anonymous; line: 14

     1     ...OGRAM0/main_farr.jl; main; line: 99

     46    ...OGRAM0/main_farr.jl; main; line: 101

     1     ...OGRAM0/main_farr.jl; main; line: 107

     5     ...OGRAM0/main_farr.jl; main; line: 110

     4     ...OGRAM0/main_farr.jl; main; line: 114

     3     ...OGRAM0/main_farr.jl; main; line: 119

     6497  ...OGRAM0/main_farr.jl; main; line: 121

       26   float.jl; *; line: 136

       24   subarray.jl; getindex; line: 176

       32   subarray.jl; setindex!; line: 449

     23678 ...OGRAM0/main_farr.jl; main; line: 126

       63    float.jl; -; line: 134

       21    float.jl; /; line: 138

       141   subarray.jl; getindex; line: 176

       29    subarray.jl; setindex!; line: 449

     2     ...OGRAM0/main_farr.jl; main; line: 130

     2     ...OGRAM0/main_farr.jl; main; line: 137

27    float.jl; *; line: 136

54    float.jl; -; line: 134

17    float.jl; /; line: 138

651   subarray.jl; getindex; line: 176

567   subarray.jl; setindex!; line: 449

Timing is different - another machine: 0.29s for the fast version Array-based version and 31s for slow SubArray-based.

I noticed that for the fast version there are no getindex, setindex calls in the profile log, looks like they were inlined and effectively eliminated as Julia has an efficient intrinsic Array support.

The difference between fast and slow codes (not full, just a few excerpts):

diff -u main.jl main_farr.jl |less

+

+function farray(T, r::Range1{Int}...)

+    dims = map((x) -> length(x), r)

+    array = Array(T, dims)

+    sub_indices = map((x) -> -minimum(x) + 2 : maximum(x), r)

+    sub(array, sub_indices)

+end

+

 function main()

@@ -33,9 +41,13 @@

    #  &  x(0:ncells),

    #  &  flux(0:ncells)

 

-      u    = Array(Float64, ncells+4)

-      x    = Array(Float64, ncells+1)

-      flux = Array(Float64, ncells+1)

+      u    = farray(Float64, -2:ncells+1)

+      x    = farray(Float64,  0:ncells)

+      flux = farray(Float64,  0:ncells)

[...]

#       conservative difference

         for ic=ifirst:ilast

-          u[ic+3] -= (flux[ic+2]-flux[ic+1]) / (x[ic+2]-x[ic+1])

+          u[ic] -= (flux[ic+1]-flux[ic]) / (x[ic+1]-x[ic])

         end

[...]

Thanks,

Alexander

[Tim Holy]

My best guess is that you have a type problem somewhere. In the profile, notice
that there are two lines showing thousands of samples (line 121 and 126), but
that the sum of samples in the "children" of these lines comes nowhere close
to the parent. That's almost surely a sign that the garbage-collector is
running a lot (you can check with ProfileView or C=true to be sure), and one
possible source of this is that some computation has ambiguous type and forces
allocation of "generic" variables. You can read about this in the manual,
particularly in the FAQ.

--Tim

[Alexander Samoilov]

Thanks for your comments.

lines 121 and 126 look innocuous:

#       upwind fluxes times dt (ie, flux time integral over cell side)

#       assumes velocity > 0

        vdt=velocity*dt

        for ie=ifirst:ilast+1

          flux[ie]=vdt*u[ie-1]  # line 121

        end

#       conservative difference

        for ic=ifirst:ilast

          u[ic] -= (flux[ic+1]-flux[ic]) / (x[ic+1]-x[ic]) # line 126

        end

In original Julia code the expressions are similar, but index expressions for u,x,flux arrays are shifted to adopt 1-based arrays.

I extremely simplified farray function and converted it into a macro:

diff -u main_farr.jl main_farrm.jl 

--- main_farr.jl        2013-12-08 13:41:48.955932728 +0400

+++ main_farrm.jl       2013-12-08 17:28:27.331785730 +0400

@@ -21,11 +21,8 @@

 #  Use the software at your own risk.

 #***********************************************************************

 

-function farray(T, r::Range1{Int}...)

-    dims = map((x) -> length(x), r)

-    array = Array(T, dims)

-    sub_indices = map((x) -> -minimum(x) + 2 : maximum(x), r)

-    sub(array, sub_indices)

+macro farray(T, r)

+    :(sub(Array($T, length($r)), -minimum($r) + 2 : maximum($r)))

 end

 

 function main()

@@ -45,9 +42,9 @@

 #      x    = Array(Float64, ncells+1)

 #      flux = Array(Float64, ncells+1)

 

-      u    = farray(Float64, -2:ncells+1)

-      x    = farray(Float64,  0:ncells)

-      flux = farray(Float64,  0:ncells)

+      u    = @farray(Float64, -2:ncells+1)

+      x    = @farray(Float64,  0:ncells)

+      flux = @farray(Float64,  0:ncells)

There is no performance improvement after this change (actually I didn't expect it but t worth trying to exclude farray from consideration).

Any ideas how to achieve a significant performance improvement (100x performance degradation is too severe),

or SubArray is not suitable at all for the tricks with the shifted  array base?

Thanks,

Alexander

[Stefan Karpinski]

SubArray is notoriously not sufficiently fast, but I'm working on a replacement that will hopefully be much faster: https://github.com/JuliaLang/julia/pull/5003.

[Tim Holy]

On Sunday, December 08, 2013 06:41:49 AM Alexander Samoilov wrote:
> Any ideas how to achieve a significant performance improvement (100x
> performance degradation is too severe),
> or SubArray is not suitable at all for the tricks with the shifted array
> base?

I really doubt that's the issue here. Copy the code below into a file and
include it:

function mysum(A)
s = zero(eltype(A))
for i = 1:length(A)-1
s += A[i]
end
s
end

A = rand(10^8)
B = sub(A, 1:length(A)-1)
Bs = sub(A, 2:length(A))

# JIT-compile
mysum(A)
mysum(B)
mysum(Bs)
@time mysum(A) # once to compile functions called by @time

# Now the real test
@time mysum(A)
@time mysum(B)
@time mysum(Bs)



Results:
julia> include("/tmp/shifted.jl")
elapsed time: 0.179086382 seconds (6740 bytes allocated)
elapsed time: 0.180591451 seconds (64 bytes allocated)
elapsed time: 0.312678742 seconds (80 bytes allocated)
elapsed time: 0.315888863 seconds (80 bytes allocated)
5.000054323638211e7

It's within a factor of 2 or so from Arrays, with no difference between the two
SubArrays. So shifting indexes is not the problem.

The short answer is that you probably have a type problem, and Leah Hanson
wrote a terrific blog post
http://blog.leahhanson.us/julia-introspects.html
describing how to use the tools that might be helpful to debug what's going
wrong.

For example, compare the output of code_typed on mysum with what happens if
you change the initialization of s to `s=0`. You'll see that this code
suddenly gets more than 100x slower, and using code_typed you'll see
Union(Int64, Float64) in places where it doesn't belong (most notably, inside
the loop).

Best,
--Tim

[Alexander Samoilov]

Thanks a lot for the information on the branch arrayviews,

I'd like to play with the code, but unfortunately don't know Julia yet, could you please help?

julia> include ("base/array/view.jl")

ArrayView{T,n,m} (constructor with 3 methods)

julia> r = -7:7

-7:7

julia> a = Array(Float64, length(r));

So good so far.

julia> av = ArrayView(a, r)

ERROR: no method ArrayView{T,n,m}(Array{Float64,1}, Range1{Int64})

Found constructor

function ArrayView{n}(a::Array, R::NTuple{n,Ranges})

But how to coerce to NTuple (failed to find it in docs)?

Best regards,

Alex

[Alexander Samoilov]

Great, thanks a lot!

I'll try to figure out what's going wrong follow your hints.

[Alexander Samoilov]

Got it, though there is now bounds checking (yet:-)

julia> r = -7:7

-7:7

julia> a = Array(Float64, length(r));

julia> av = ArrayView(a, (r,))

ArrayView{Float64,1,1}([0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0],(15,),(1,),-7)

julia> av[1]

8.8874505e-317

julia> av[0]

8.8874505e-317

julia> av[-7]

8.8874505e-317

julia> av[-8]

8.8874505e-317

julia> av[-9]

8.8874505e-317

[Stefan Karpinski]

That branch is nowhere near using, so please don't interpret it as usable. As the title says, it's an exploratory branch.

t; > > >>> >    4. Any other comments, proposals would be greatly appreciated.
> > > >>> >
> > > >>> > Thanks!
> > > >>> > Alexander

>

[Alexander Samoilov]

Sure, I got, that it is not usable at the moment.

But the concept is great.

[Alexander Samoilov]

Hello Tim,

I succeeded in reproducing erratic behaviour:

The reduced testcase is below derived from your mysum kernel you have sent:

function mykern(r, s, a, x, y)

    const nsteps = 25 

    istep = 0

    ifirst = minimum(r)

    ilast  = maximum(r)

    while istep < nsteps 

        #for i = r

        for i = ifirst : ilast # here is the problem, here 100x perf degradation compared to normal array if shift is not 1

                               # if you change to line above all will be good, only 2x perf penalty as expected for shift other than 1

            s[i] = a * (x[i+1] - x[i]) / (y[i+1] - y[i])

        end

        istep=istep + 1

    end 

    s

end 

const size = 10^6

const a = e^pi - pi^e

s = rand(size)

x = rand(size)

y = rand(size)

s1 = sub(s, 1:(size - 1))

x1 = sub(x, 1:(size - 1))

y1 = sub(y, 1:(size - 1))

s2 = sub(s, 2:size)

x2 = sub(x, 2:size)

y2 = sub(y, 2:size)

# JIT-compile 

mykern( 2:(size - 3), s, a, x,  y  ) 

mykern( 2:(size - 3), s, a, x1, y1 ) 

mykern( 3:(size - 2), s, a, x2, y2 ) 

@time mykern( 2:(size - 3), s, a, x,  y )  # once to compile functions called by @time 

# Now the real test 

@time mykern( 2:(size - 3), s, a, x,  y  ) 

@time mykern( 2:(size - 3), s, a, x1, y1 ) 

@time mykern( 3:(size - 2), s, a, x2, y2 ) 

The root cause is the loop construct written like

for i = ifirst : ilast

with it

julia shifted2.jl 

elapsed time: 8.219997503 seconds (2799181116 bytes allocated)

elapsed time: 8.190918047 seconds (2799174504 bytes allocated)

elapsed time: 8.445548648 seconds (2799174504 bytes allocated)

elapsed time: 8.511262972 seconds (2799176104 bytes allocated)

the litmus paper - huge amount of bytes allocated reported.

The same 100x perf degradation compared to the case when

the loop written as

for i = r

the output here is:

julia shifted2.jl 

elapsed time: 0.186420902 seconds (6700 bytes allocated)

elapsed time: 0.187166614 seconds (88 bytes allocated)

elapsed time: 0.249333829 seconds (88 bytes allocated)

elapsed time: 0.252043335 seconds (88 bytes allocated)

Any ideas what is going wrong?

Thanks,

Alexander

On Sunday, December 8, 2013 9:19:57 PM UTC+4, Tim wrote:

[Stefan Karpinski]

One slightly questionable thing here is using the name "size" as a variable. This is fine, but if you do this, you're not going to be able to use the Base.size function.

[Tim Holy]

The problem is that the inferred return type of minimum(r::Range1{Int}) is
Union(Int64, Bool). This is a julia bug; can you please file an issue? Or, if
you want to fix it yourself, the place to look is line 100 of base/range.jl.

Good detective work!

Best,
--Tim

On Monday, December 09, 2013 07:46:40 AM Alexander Samoilov wrote:
> Hello Tim,
>
> I succeeded in reproducing erratic behaviour:
>
> The reduced testcase is below derived from your mysum kernel you have sent:
>
> function mykern(r, s, a, x, y)
> const nsteps = 25
> istep = 0
> ifirst = minimum(r)
> ilast = maximum(r)
> while istep < nsteps
> #for i = r

> *for i = ifirst : ilast # here is the problem, here 100x perf
> degradation compared to normal array if shift is not 1*
> * # if you change to line above all will be
> good, only 2x perf penalty as expected for shift other than 1*

[Stefan Karpinski]

Possibly the same issue: https://github.com/JuliaLang/julia/issues/5079.

[Alexander Samoilov]

Indeed maybe the same issue.

I'll file an issue nevertheless referring to the issue 5079 and changing size to say siz as Stefan proposed.

:%s/size/siz/gc

Best regards,

Alexander

[Alexander Samoilov]

I submitted an issue:

https://github.com/JuliaLang/julia/issues/5084

Thanks!

Alexander

On Monday, December 9, 2013 11:45:16 PM UTC+4, Stefan Karpinski wrote:

[Alexander Samoilov]

The issue #5084 was rapidly and precisely fixed by Stefan.

Created another issue as still have performance problems with sub:

https://github.com/JuliaLang/julia/issues/5100

Thanks!

[Kevin Squire]

Actually, to be clear: if you have a variable size and want to call Base.size(...), you can directly call Base.size(...).  You just can't call size(...) by itself.

Kevin

[Alexander Samoilov]

Thanks to all for their feedback.

Filed a new issue for SubArray performance improvement:

https://github.com/JuliaLang/julia/issues/5117

[Alexander Samoilov]

I came to the following solution for Fortran-like arrays that faster than existing SubArray-based and has about 2x performance penalty (checked for 1d vectors only).

type FArray{T<:Number, N, A<:AbstractArray} <: AbstractArray

    ranges  #::I

    offsets::NTuple{N,Int}

    array::A

    function FArray(r::Range1{Int}...)

        dims = map((x) -> length(x), r)

        array = Array(T, dims)

        offs = map((x) -> 1 - minimum(x), r)

        new(r, offs, array)

    end

end

FArray(T, r::Range1{Int}...) = FArray{T, length(r,), Array{T, length(r,)}}(r...)

getindex{T<:Number}(FA::FArray{T}, i1::Int) = FA.array[i1+FA.offsets[1]]

getindex{T<:Number}(FA::FArray{T}, i1::Int, i2::Int) = FA.array[i1+FA.offsets[1], i2+FA.offsets[2]]

getindex{T<:Number}(FA::FArray{T}, i1::Int, i2::Int, i3::Int) = FA.array[i1+FA.offsets[1], i2+FA.offsets[2], i3+FA.offsets[3]]

setindex!{T}(FA::FArray{T}, x, i1::Int) = arrayset(FA.array, convert(T,x), i1+FA.offsets[1])

setindex!{T}(FA::FArray{T}, x, i1::Int, i2::Int) = arrayset(FA.array, convert(T,x), i1+FA.offsets[1], i2+FA.offsets[2])

setindex!{T}(FA::FArray{T}, x, i1::Int, i2::Int, i3::Int) = arrayset(FA.array, convert(T,x), i1+FA.offsets[1], i2+FA.offsets[2], i3+FA.offsets[3])

And use case:

julia> y = FArray(Float64, -1:1,-7:7,-128:512)

julia> y[-1,-7,-128] = 777

julia> y[-1,-7,-128] + 88

865.0

julia> y[-2,-7,-128]

ERROR: BoundsError()

 in getindex at /home/asamoilov/work/projects/SD/julia.sandbox/FortranArray/src/FortranArray.jl:27

julia> y[2,-7,-128]

ERROR: BoundsError()

 in getindex at /home/asamoilov/work/projects/SD/julia.sandbox/FortranArray/src/FortranArray.jl:27

I understand that the constructor is not efficient as uses functional programming features such as maps, lambdas - Julia is inefficient here.

But the constructor is called only once during array construction and cannot influence significantly on the performance in my understanding. Am I right?

Could you please propose ways to improve efficiency of the above code?

Thanks!

Alexander

[Alexander Samoilov]

I created a gist for Julia implementation of Fortran-like arrays with arbitrary starting indices (negative or zero):

https://gist.github.com/alsam/8283205

Your comments would be greatly appreciated.

Thanks,

Alexander