# Porting Fortran code, negative indices

**URL:** <https://discourse.julialang.org/t/porting-fortran-code-negative-indices/72386>\
**Category:** General Usage\
**Tags:** fortran, indexing\
**Created:** [December 1, 2021, 4:16pm UTC](https://discourse.julialang.org/t/porting-fortran-code-negative-indices/72386 "2021-12-01T16:16:57Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![baptnz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baptnz/32/24519_2.png) [@baptnz](https://discourse.julialang.org/u/baptnz)\
**Post date:** [December 1, 2021, 4:16pm UTC](https://discourse.julialang.org/t/porting-fortran-code-negative-indices/72386/1 "2021-12-01T16:16:57Z")

</div>

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

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

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

```

becomes

```julia
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?

---

<div class="post-metadata">

**Author:** ![jling](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jling/32/212909_2.png) [@jling](https://discourse.julialang.org/u/jling)\
**Post date:** [December 1, 2021, 4:20pm UTC](https://discourse.julialang.org/t/porting-fortran-code-negative-indices/72386/2 "2021-12-01T16:20:20Z")

</div>

> [@baptnz](#):
>
> `ro[-5:5]`

this is not right:

```julia
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:

```

> [@baptnz](#):
>
> into `OffsetArray` type: I fear it will make it harder to work with special types

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.

---

<div class="post-metadata">

**Author:** ![baptnz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baptnz/32/24519_2.png) [@baptnz](https://discourse.julialang.org/u/baptnz)\
**Post date:** [December 1, 2021, 4:49pm UTC](https://discourse.julialang.org/t/porting-fortran-code-negative-indices/72386/3 "2021-12-01T16:49:05Z")

</div>

> [@jling](#):
>
> this is not right:

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

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

```

> [@jling](#):
>
> fear not, multi-dispatch is born for this

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.

---

<div class="post-metadata">

**Author:** ![jling](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jling/32/212909_2.png) [@jling](https://discourse.julialang.org/u/jling)\
**Post date:** [December 1, 2021, 4:56pm UTC](https://discourse.julialang.org/t/porting-fortran-code-negative-indices/72386/4 "2021-12-01T16:56:39Z")

</div>

> [@baptnz](#):
>
> `getindex()` behavior, e.g. for sparse matrices.

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`

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [December 1, 2021, 8:28pm UTC](https://discourse.julialang.org/t/porting-fortran-code-negative-indices/72386/5 "2021-12-01T20:28:21Z")

</div>

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:

```julia
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.

---

<div class="post-metadata">

**Author:** ![baptnz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baptnz/32/24519_2.png) [@baptnz](https://discourse.julialang.org/u/baptnz)\
**Post date:** [December 2, 2021, 12:25pm UTC](https://discourse.julialang.org/t/porting-fortran-code-negative-indices/72386/6 "2021-12-02T12:25:39Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [December 2, 2021, 12:30pm UTC](https://discourse.julialang.org/t/porting-fortran-code-negative-indices/72386/7 "2021-12-02T12:30:55Z")

</div>

Forgot this other option:

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

```

---

<div class="post-metadata">

**Author:** ![jishnub](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jishnub/32/33620_2.png) [@jishnub](https://discourse.julialang.org/u/jishnub)\
**Post date:** [December 2, 2021, 12:40pm UTC](https://discourse.julialang.org/t/porting-fortran-code-negative-indices/72386/8 "2021-12-02T12:40:53Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![baptnz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baptnz/32/24519_2.png) [@baptnz](https://discourse.julialang.org/u/baptnz)\
**Post date:** [December 2, 2021, 1:13pm UTC](https://discourse.julialang.org/t/porting-fortran-code-negative-indices/72386/9 "2021-12-02T13:13:38Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![woclass](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/woclass/32/212699_2.png) [@woclass](https://discourse.julialang.org/u/woclass)\
**Post date:** [December 2, 2021, 2:06pm UTC](https://discourse.julialang.org/t/porting-fortran-code-negative-indices/72386/10 "2021-12-02T14:06:00Z")

</div>

> <https://github.com/JuliaArrays/OffsetArrays.jl/issues/91>
>
> \`\`\`
> julia\> A = OffsetArray(\[1 2; 3 4\], 0:1, 0:1)
> 2×2 OffsetArray(::Array{Int64…,2}, 0:1, 0:1) with eltype Int64 with indices 0:1×0:1:
> 1 2
> 3 4
> 
> julia\> b = OffsetArray(\[5; 6\], 0:1)
> 2-element OffsetArray(::Array{Int64,1}, 0:1) with eltype Int64 with indices 0:1:
> 5
> 6
> 
> julia\> A\*b
> ERROR: ArgumentError: offset arrays are not supported but got an array with index other than 1
> Stacktrace:
> \[1\] require\_one\_based\_indexing at .\\abstractarray.jl:89 \[inlined\]
> \[2\] generic\_matvecmul!(::OffsetArray{Int64,1,Array{Int64,1}}, ::Char, ::OffsetArray{Int64,2,Array{Int64,2}}, ::OffsetArray{Int64,1,Array{Int64,1}}) at C:\\cygwin\\home\\Administrator\\buildbot\\worker\\package\_win64\\build\\usr\\share\\julia\\stdlib\\v1.2\\LinearAlgebra\\src\\matmul.jl:501
> \[3\] mul! at C:\\cygwin\\home\\Administrator\\buildbot\\worker\\package\_win64\\build\\usr\\share\\julia\\stdlib\\v1.2\\LinearAlgebra\\src\\matmul.jl:77 \[inlined\]
> \[4\] \*(::OffsetArray{Int64,2,Array{Int64,2}}, ::OffsetArray{Int64,1,Array{Int64,1}}) at C:\\cygwin\\home\\Administrator\\buildbot\\worker\\package\_win64\\build\\usr\\share\\julia\\stdlib\\v1.2\\LinearAlgebra\\src\\matmul.jl:51
> \[5\] top-level scope at none:0
> \`\`\`

> 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](https://github.com/JuliaLang/julia/search?q=require_one_based_indexing)

---

<div class="post-metadata">

**Author:** ![baptnz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baptnz/32/24519_2.png) [@baptnz](https://discourse.julialang.org/u/baptnz)\
**Post date:** [December 2, 2021, 2:32pm UTC](https://discourse.julialang.org/t/porting-fortran-code-negative-indices/72386/11 "2021-12-02T14:32:04Z")

</div>

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.