changing array dimensions in fortran

Question

There are basically two ways to pass arrays to a subroutine in Fortran 90/95:

PROGRAM ARRAY
INTEGER, ALLOCATABLE :: A(:,:)
INTEGER :: N
ALLOCATE(A(N,N))
CALL ARRAY_EXPLICIT(A,N)
! or
CALL ARRAY_ASSUMED(A)
END PROGRAM ARRAY

SUBROUTINE ARRAY_EXPLICIT(A,N)
INTEGER :: N
INTEGER :: A(N,N)
! bla bla
END SUBROUTINE ARRAY_EXPLICIT

SUBROUTINE ARRAY_ASSUMED(A)
INTEGER, ALLOCATABLE :: A(:,:)
N=SIZE(A,1)
! bla bla
END SUBROUTINE ARRAY_ASSUMED

where you need an explicit interface for the second, usually through the use of a module.

From FORTRAN77, I'm used to the first alternative, and I read this is also the most efficient if you pass the whole array.

The nice thing with the explicit shape is that I can also call a subroutine and treat the array as a vector instead of a matrix:

SUBROUTINE ARRAY_EXPLICIT(A,N)
INTEGER :: N
INTEGER :: A(N**2)
! bla bla
END SUBROUTINE ARRAY_EXPLICIT

I wondered if there is a nice way to do that kind of thing using the second, assumed shape interface, without copying it.

Emirate

"without copying it." Without copying what? – Drabble

I think that @Emirate means changing the shape of the array in place, rather than copying it to a 1D array. – Cowbird

thnx @M. S. B. for the clarification, that's what I meant. So not the reshape intrinsic solutions mentioned below – Emirate

@ M. S. B. I suspected it. Than reshape intrinsic is not a solution. Well, I think that things depends heavily on "! bla, bla" part of the story. =) You always can just iterate over your 2D array as it is 1D using loops. – Drabble

@kemiisto true, you can always control how you index the thing, but I wondered if this 'manual' treatment is the only way. Say, I have a subroutine accepting a NxN matrix and the matrix I want to pass can be most easily filled by using 4 index-loops on an MxMxMxM shape, then I can elegantly use an explicit shape dummy. – Emirate

Answer 1

See the RESHAPE intrinsic, e.g.

http://gcc.gnu.org/onlinedocs/gfortran/RESHAPE.html

Alternatively, if you want to avoid the copy (in some cases an optimizing compiler might be able to do a reshape without copying, e.g. if the RHS array is not used afterwards, but I wouldn't count on it), as of Fortran 2003 you can assign pointers to targets of different rank, using bounds remapping. E.g. something like

program ptrtest
  real, pointer :: a(:)
  real, pointer :: b(:,:)
  integer :: n = 10
  allocate(a(n**2))
  a = 42
  b (1:n, 1:n) => a
end program ptrtest

Flamsteed

+1, however, your example doesn't work in my case, I compiled with ifort 11. The reported size of b is 2*n**2, double the size of a. – Emirate

Only the very latest compilers support this feature of Fortran 2003, called pointer bounds remapping. I think gfortran 4.6 but not 4.5, and Intel 12 but not 11.1. – Cowbird

@M. S. B. thnx! BTW, now that you say, it was indeed ifort 12, not 11. I still didn't get it to work though, very strange. Could someone check this with ifort 12? – Emirate

I just tested the ptrtest program above using Intel's fortran compiler on Mac OSX 10.11.6 (ifort 17.0.7 20180403). I compiled once with and once without the following flags: debug all -g -traceback -gen-interfaces -warn interfaces -check bounds -check arg_temp_created -check uninit -check all,noarg_temp_created. In both cases, it compiles and runs without warning or error. – Baler

Answer 2

I was looking to do the same thing and came across this discussion. None of the solutions suited my purposes, but I found that there is a way to reshape an array without copying the data using iso_c_binding if you are using the fortran 2003 standard which current fortran 90/95 compilers tend to support. I know the discussion is old, but I figured I would add what I came up with for the benefit of others with this question.

The key is to use the function C_LOC to convert an array to an array pointer, and then use C_F_POINTER to convert this back into a fortran array pointer with the desired shape. One challenge with using C_LOC is that C_LOC only works for array that have a directly specified shape. This is because arrays in fortran with an incomplete size specification (i.e., that use a : for some dimension) include an array descriptor along with the array data. C_LOC does not give you the memory location of the array data, but the location of the descriptor. So an allocatable array or a pointer array don't work with C_LOC (unless you want the location of the compiler specific array descriptor data structure). The solution is to create a subroutine or function that receives the array as an array of fixed size (the size really doesn't matter). This causes the array variable in the function (or subroutine) to point to the location of the array data rather than the location of the array descriptor. You then use C_LOC to get a pointer to the array data location and C_F_POINTER to convert this pointer back into an array with the desired shape. The desired shape must be passed into this function to be used with C_F_POINTER. Below is an example:

program arrayresize
  implicit none
  integer, allocatable :: array1(:)
  integer, pointer :: array2(:,:)

  ! allocate and initialize array1
  allocate(array1(6))
  array1 = (/1,2,3,4,5,6/)

  ! This starts out initialized to 2
  print *, 'array1(2) = ', array1(2)

  ! Point array2 to same data as array1. The shape of array2
  ! is passed in as an array of intergers because C_F_POINTER
  ! uses and array of intergers as a SIZE parameter.
  array2 => getArray(array1, (/2,3/))

  ! Change the value at array2(2,1) (same as array1(2))
  array2(2,1) = 5

  ! Show that data in array1(2) was modified by changing
  ! array2(2,1)
  print *, 'array(2,1) = array1(2) = ', array1(2)

contains

  function getArray(array, shape_) result(aptr)
    use iso_c_binding, only: C_LOC, C_F_POINTER
    ! Pass in the array as an array of fixed size so that there
    ! is no array descriptor associated with it. This means we
    ! can get a pointer to the location of the data using C_LOC
    integer, target :: array(1)
    integer :: shape_(:)
    integer, pointer :: aptr(:,:)

    ! Use C_LOC to get the start location of the array data, and
    ! use C_F_POINTER to turn this into a fortran pointer (aptr).
    ! Note that we need to specify the shape of the pointer using an
    ! integer array.
    call C_F_POINTER(C_LOC(array), aptr, shape_)
  end function
end program

Rectory

this is similar to janneb's answer: "as of Fortran 2003 you can assign pointers to targets of different rank", and it also relies on F2003 features, but only more complicated... – Emirate

Answer 3

@janneb has already answered re RESHAPE. RESHAPE is a function -- usually used in an assignment statement so there will be a copy operation. Perhaps it can be done without copying using pointers. Unless the array is huge, it is probably better to use RESHAPE.

I'm skeptical that the explicit shape array is more efficient than the assumed shape, in terms of runtime. My inclination is to use the features of the Fortran >=90 language and use assumed shape declarations ... that way you don't have to bother passing the dimensions.

EDIT: I tested the sample program of @janneb with ifort 11, gfortran 4.5 and gfortran 4.6. Of these three, it only works in gfortran 4.6. Interestingly, to go the other direction and connect a 1-D array to an existing 2-D array requires another new feature of Fortran 2008, the "contiguous" attribute -- at least according to gfortran 4.6.0 20110318. Without this attribute in the declaration, there is a compile time error.

    program test_ptrs

   implicit none

   integer :: i, j

   real, dimension (:,:), pointer, contiguous :: array_twod
   real, dimension (:), pointer :: array_oned

   allocate ( array_twod (2,2) )

   do i=1,2
      do j=1,2
         array_twod (i,j) = i*j
      end do
   end do

   array_oned (1:4) => array_twod

   write (*, *) array_oned

   stop

end program test_ptrs

Cowbird

Yes, I know the reshape function, but this is, like you say, problematic for huge arrays, e.g. if 2 of them don't fit in memory. The fact that explicit shape is more efficient, I read that here software.intel.com/file/6397 on page40-41. You are right about using the new features of Fortran>=90, that's why I want to use assumed shape for dummys and use an explicit interface. That's why I wanted to know about the reshaping, if it's possible in-place. – Emirate

Answer 4

You can use assumed-size arrays, but it can mean multiple layers of wrapper routines:

program test

  implicit none

  integer :: test_array(10,2)

  test_array(:,1) = (/1,   2,  3,  4,  5,  6,  7,  8,  9, 10/)
  test_array(:,2) = (/11, 12, 13, 14, 15, 16, 17, 18, 19, 20/)

  write(*,*) "Original array:"
  call print_a(test_array)

  write(*,*) "Reshaped array:"
  call print_reshaped(test_array, size(test_array))

contains

  subroutine print_reshaped(a, n)
  integer, intent(in) :: a(*)
  integer, intent(in) :: n
  call print_two_dim(a, 2, n/2)
  end subroutine

  subroutine print_two_dim(a, n1, n2)
  integer, intent(in) :: a(1:n1,1:*)
  integer, intent(in) :: n1, n2
  call print_a(a(1:n1,1:n2))
  end subroutine

  subroutine print_a(a)
  integer, intent(in) :: a(:,:)
  integer :: i
  write(*,*) "shape:", shape(a)
  do i = 1, size(a(1,:))
      write(*,*) a(:,i)
  end do
  end subroutine

end program test

Tomekatomes

Answer 5

I am using ifort 14.0.3 and 2D to 1D conversion, I could use an allocatable array for 2D array and a pointer array for 1D:

integer,allocatable,target :: A(:,:)
integer,pointer :: AP(:)

allocate(A(3,N))
AP(1:3*N) => A

As @M.S.B mentioned, in case both A and AP have the pointer attribute, I had to use contiguous attribute for A to guarantee the consistency of the conversion.

Fouquet

Answer 6

Gfortran is a bit paranoid with interfaces. It not only wants to know the type, kind, rank and number of arguments, but also the shape, the target attribute and the intent (although I agree with the intent part). I encountered a similar problem.

With gfortran, there are three different dimension definition:
1. Fixed
2. Variable
3. Assumed-size

With ifort, categories 1 and 2 are considered the same, so you can do just define any dimension size as 0 in the interface and it works.

program test

  implicit none

  integer, dimension(:), allocatable :: ownlist

  interface
    subroutine blueprint(sz,arr)
      integer, intent(in) :: sz
      integer, dimension(0), intent(in) :: arr
      ! This zero means that the size does not matter,
      ! as long as it is a one-dimensional integer array.
    end subroutine blueprint
  end interface

  procedure(blueprint), pointer :: ptr

  allocate(ownlist(3))
  ownlist = (/3,4,5/)
  ptr => rout1
  call ptr(3,ownlist)
  deallocate(ownlist)

  allocate(ownlist(0:10))
  ownlist = (/3,4,5,6,7,8,9,0,1,2,3/)
  ptr => rout2
  call ptr(3,ownlist)
  deallocate(ownlist)

contains

  ! This one has a dimension size as input.
  subroutine rout1(sz,arr)
    implicit none
    integer, intent(in) :: sz
    integer, dimension(sz), intent(in) :: arr
    write(*,*) arr
    write(*,*) arr(1)
  end subroutine rout1

  ! This one has a fixed dimension size.
  subroutine rout2(sz,arr)
    implicit none
    integer, intent(in) :: sz
    integer, dimension(0:10), intent(in) :: arr
    write(*,*) "Ignored integer: ",sz
    write(*,*) arr
    write(*,*) arr(1)
  end subroutine rout2

end program test

Gfortran complains about the interface. Changing the 0 into 'sz' solves the problem four 'rout1', but not for 'rout2'.

However, you can fool gfortran around and say dimension(0:10+0*sz) instead of dimension(0:10) and gfortran compiles and gives the same result as ifort.

This is a stupid trick and it relies on the existence of the integer 'sz' that may not be there. Another program:

program difficult_test

  implicit none

  integer, dimension(:), allocatable :: ownlist

  interface
    subroutine blueprint(arr)
      integer, dimension(0), intent(in) :: arr
    end subroutine blueprint
  end interface

  procedure(blueprint), pointer :: ptr

  allocate(ownlist(3))
  ownlist = (/3,4,5/)
  ptr => rout1
  call ptr(ownlist)
  deallocate(ownlist)

  allocate(ownlist(0:10))
  ownlist = (/3,4,5,6,7,8,9,0,1,2,3/)
  ptr => rout2
  call ptr(ownlist)
  deallocate(ownlist)

contains

  subroutine rout1(arr)
    implicit none
    integer, dimension(3), intent(in) :: arr
    write(*,*) arr
    write(*,*) arr(1)
  end subroutine rout1

  subroutine rout2(arr)
    implicit none
    integer, dimension(0:10), intent(in) :: arr
    write(*,*) arr
    write(*,*) arr(1)
  end subroutine rout2

end program difficult_test

This works under ifort for the same reasons as the previous example, but gfortran complains about the interface. I do not know how I can fix it.

The only thing I want to tell gfortran is 'I do not know the dimension size yet, but we will fix it.'. But this needs a spare integer arguemnt (or something else that we can turn into an integer) to fool gfortran around.

Cramped

It is quite long, but I am not sure if you are answering to the question actually asked by the OP. – Grau