Matlab sub2ind in Fortran

1

Is there a Fortran implementation of the Matlab function sub2ind out there? Simple google search did not help. The function is very simple, it converts subscripts of a multidimensional array into linear indexes. For example, if A is a 2 dim array with n1 rows and n2 columns, the subscripts (i1,i2) are mapped to

Lin = (i2-1)*nx+i1

I wrote sub2ind in Fortran myself for the case of 2 dim arrays but I would like to generalize up to, say, 4 or 5 dimensions. What I plan to do is two implement sub2ind for 2,3,4,5 dimensions separatley and then use


Module mymod

Implicit none

interface sub2ind
module procedure sub2ind2
module procedure sub2ind3
Etc
end interface

contains

function sub2ind2(nvec,i1,i2) result(lin)

integer, intent(in) :: nvec(2)
integer, intent(in) :: i1,i2
Integer :: lin

Lin = (i2-1)*nvec(1) + i1


end function sub2ind2

End module mymod

However, I would like to have a well-tested code with error messages etc, so I figured that before I reinvent the wheel I should ask on this forum

2 Likes
2

I’ve been thinking about writing a recursive template for this using Fypp: Introduction — Fypp 3.2 documentation (see comment near the bottom of def about recursive invocation). You can also think of this as a fold operation.

As a side-note, not directly answering your question, but potentially offering an alternative solution, Fortran has the notion pointer array remapping:

real, target :: a(100) ! Pretend it's a 5-by-20 array
Real, pointer :: b(:,:) => null()

b(1:5,1:20) => a

Now you can reference elements of b without needing to do the linear index computation.

2 Likes
3

maybe useful?

1 Like
4

Here is my version but its not much different than the OPs version.

  PURE Subroutine sub2ind(isize, col, row, ind)

    Implicit NONE

    Integer,  Intent(IN)    :: isize(2)
    Integer,  Intent(IN)    :: col(:)
    Integer,  Intent(IN)    :: row(:)
    Integer,  Intent(INOUT) :: ind(:)

    Integer :: i, j, ie, ni

    ie = SIZE(col)
    ni = isize(1)

    Do i=1,ie
      ind(i) = col(i) + (row(i)-1)*ni
    End Do

  End Subroutine sub2ind
1 Like
5

Thanks! but this works only for two dimensions, right?

6

You may want to also check out the pm_matrixIndex :: getMatIndex() of the ParaMonte library. It has a comprehensive set of index-conversion routines for many of the major LAPACK packing formats, including the one that you are looking for (called lfpack in the library). If you would like to simply extract the module and use it in your codebase, try the following in a Git-aware terminal:

git clone https://github.com/cdslaborg/paramonte.git
cd paramonte
./install.sh --mod pm_matrixIndex

The last command will copy the pm_matrixIndex and all its dependencies for inclusion in your project in the installation directory ./bin/<library_build_specs>/fpp in the root directory of the Git project.

If you simply prefer to reimplement the algorithm in your codebase, here is the actual implementation. Even if you intend to reimplement, I highly recommend to work with preprocessed source files using the above commands instead of the actual implementation source code linked above.

2 Likes
7

Here is rather odd-looking way to do this in Fortran 2018 using reduce (and probably not very efficient based on what I saw in Compiler Explorer). A newer version of the Intel Fortran compiler is needed.

module fold
    implicit none
    public
contains
    function flatidx(extent,index) result(idx)
        integer, intent(in) :: extent(:), index(:)
        integer :: idx
        type :: pair
            integer :: e, i
        end type
        type(pair) :: tmp(size(extent))
        type(pair) :: res
        tmp%e = extent
        tmp%i = index
        res = reduce(tmp,op)
        idx = res%i
    contains
       pure function op(a,b) result(c)
        type(pair), intent(in) :: a, b
        type(pair) :: c
        c = pair(a%e*b%e, a%i + a%e*(b%i - 1))
       end function
    end function
end module

program test_fold
use fold, only: flatidx
implicit none
integer :: i, j

! 1-d
do i = 1, 8
    print *, i, flatidx([8],[i])
end do

! 2-d
do j = 1, 8
    do i = 1, 4
        print *, i, j, flatidx([4,8],[i,j])
    end do
end do
end program

The solution is inspired by the blog posts:

8

I guess this also works for any rank, but I did not test carefully.

Edit: Just checked against Octave Online and I do get the same results.


module array
   implicit none
   private
   public :: sub2ind

contains

   integer function sub2ind(a, indices) result(res)
      integer, intent(in) :: a(..)
      integer, intent(in) :: indices(:)

      integer, allocatable :: a_shape(:)
      integer :: a_rank, stride, i

      a_shape = shape(a)
      a_rank = rank(a)

      res = 1
      stride = 1
      do i = 1, a_rank
         res = res + (indices(i) - 1)*stride
         stride = stride*a_shape(i)
      end do
   end function

end module array

program test
   use array, only: sub2ind
   implicit none

   integer :: a(2, 3), b(3, 2, 4), c(3, 4, 5, 6)

   print *, sub2ind(a, [1, 3])
   print *, sub2ind(a, [2, 2])

   print *, sub2ind(b, [3, 1, 1])
   print *, sub2ind(b, [3, 2, 1])
   print *, sub2ind(b, [1, 2, 1])
   print *, sub2ind(b, [3, 2, 4])

   print *, sub2ind(c, [1, 2, 4, 6])

end program
4 Likes
9

I created my own library to shift between array and index representation:

A minimal example,

program example
  use MAC, only: container_specifier
  implicit none
  type(container_specifier) :: a

  !Create a 2 x 3 x 4 3-dimensional array handle.
  call a%specify(dimension_specifier=[2, 3, 4])
  print*, a%ind([2, 1, 3]) !14.
  print*, a%ind([3, 1, 3]) !Out of bounds (returns 0).

  !Create a 2 x 5 2-dimensional array handle.
  call a%specify(dimension_specifier=[2, 5])
  print*, a%ind([1,3]) !5.
end program example

I took the idea from this website.

2 Likes
10

You can also specify that the array is in row-major order

 call a%specify(dimension_specifier=[2, 5], layout="row")
 print*, a%ind([1,3]) !3.
1 Like
11

Nice code! I would add a check on the input arguments though. You want to make sure that if the array a has n dimensions, then the vector of subscripts has n elements

if (rank(a) /= size(indices)) then
error stop "some message"
endif