Porting Fortran code, negative indices

General Usage
,
1

I’m porting a low-level program from Fortran, and ideally I would only make minimal changes in the conversion (do loop becoming for kind of thing) to help track down any changes between the two versions. (Later on, once the program produces correct results, I might look into making it more julia-idiomatic, though it’s mostly matrix operations and for loops over many indices, so I don’t think there’s much to change).
At any rate, I’m finding out that Fortran has a nice feature which is arbitrary array indexing, such as -n:n, which translates quite nicely from the mathematical description in this particular use-case.

A moment on google suggests long discussions about the topic and a reference to OffsetArrays, which implements this kind of arbitrary indexing. Now, maybe I’m overthinking this, but I’m reluctant to convert my arrays into OffsetArray type: I fear it will make it harder to work with special types, such as sparse matrices, StaticArrays, BlockArrays, etc., (I’m not sure how the interaction would work between all these types), and, just as importantly, for the most part I’m very happy to reason with the standard 1-based indexing. It’s just a few places where I’d prefer to stick to the original code’s negative-index conventions.

I’m wondering if a safe and straight-forward compromise would be to stick with standard 1-based indexing, but have a helper function convert the original Fortran indices to their 1-based equivalent. I believe that’s the reverse of OffsetArrays.IdOffsetRange,

function fortran_index(range = 1:5, start = 1)
    return range .- (start -1)
end

# Fortran loop 
# do i=-5,5
# A(i) = ...
# end do

becomes

ro = fortran_index(range = -5:5, start = -5)
# use it in julia loop with 1-based arrays
for i in ro[-5:5]
  A[i] = ...
end

Does this make sense at all, or am I shooting myself in the foot/overthinking it all?

2

this is not right:

julia> ro = fortran_index(-5:5,-5)
1:11

julia> ro[-5:5]
ERROR: BoundsError: attempt to access 11-element UnitRange{Int64} at index [-5:5]
Stacktrace:

fear not, multi-dispatch is born for this, if Fortran didn’t struggle, I would bet it takes minimal effort for OffsetArray to work. OffsetArray wraps any Array, it only changes the getindex() behavior, so you don’t need to worry about underlying Array as long as you don’t overly constraint your function argument type.

3

Oops, good point – goes to show my brain cannot be trusted with indices. I guess I meant,

function fortran_index(range = -5:5)
 OffsetArray(1:length(range), range)           
end

ro = fortran_index(range = -5:5)
for i in ro[-5:5]
  A[i] = ...
end

That’s fair, but in this instance part of me still worries that some package might have its own (incompatible) getindex() behavior, e.g. for sparse matrices.

4

when you index a sparse matrix, S[idx], it gets lowered to getindex(S, idx), when you wrap it inside an OffsetArray, and do OA[idx2], all it is doing is to translate idx2 → idx3 according to your offset and then get index S[idx3]. So there’s not much that can go wrong, and there’s no “incompatible” because everything is just getindex

4 Likes
5

We never know how familiar with Julia the OP is, but just to point out that there are many alternatives to write the loops in a index-style-agnostic way:

for i in eachindex(v)
end

for (i,val) in pairs(v)
end

for val in v
end

such that if loops run on the complete arrays, these options can be used without relying on how the original array was indexed.

5 Likes
6

Thanks but I think this wouldn’t change the issue of offset indexing of arrays – I’m happy with the for loops to be either for i in -5:5 or those more elegant alternatives, but eventually when I get/set array slices with index A[i] that’s where the offset needs to be taken into account one way or another.

The Fortran code is doing a lot of computations on indices with nested for loops (Wigner 3j symbols kind of thing) and assigning values in arrays that are best described (mathematically) with a symmetric index -n:n (spherical harmonic order),

1 Like
7

Forgot this other option:

julia> x = OffsetArray(1:5,-2:2)
1:5 with indices -2:2

julia> for i in firstindex(x):lastindex(x)
           @show i
       end
i = -2
i = -1
i = 0
i = 1
i = 2
1 Like
8

Offsetarrays do indeed have certain limitations, eg. Linear algebra, that often requires 1-based indices. However, indexing is not one of them. Offsetarrays always pass the indices to the parent array after applying a shift, so any special indexing behaviour of the parent (eg. Sparse, static) is respected.

OffsetArrays seem ideally suited for this application, and given that it’s trivial to switch between the OffsetArray, its parent, and a non-offset view, there’s nothing to lose by using an OffsetArray here.

3 Likes
9

Thanks – there is some linear algebra in various parts of the code, can you please elaborate a little on the kinds of gotchas/incompatibilities in this area?

10

The easy version is to specialize it for two OffsetArray inputs with parents that have 1-based indexing, in which case you can strip the wrappers, call the ordinary routines, and then re-wrap.

Perhaps requiring the index to start at 1 is part of the LAPACK, BLAS API.

https://github.com/JuliaLang/julia/search?q=require_one_based_indexing

1 Like
11

Thanks!

For my own peace of mind I might stick with the little conversion function between offset and natural indices, because at some point this will come to bite me. The code definitely uses quite a bit of BLAS/Lapack and I can’t rely on my programming brain to spot where I should switch back to regular arrays.