Question

How can one combine bunch of large allocatable vectors in an array in Fortran? I'd like to avoid copying and reshaping using reshape as the arrays are large. The effect I want to achieve is just like Fortran's equivalence, to illustrate:

program test_equiv
    integer x(10), y(10), z(10), xyz(10,3)
    equivalence (x, xyz(1,1))
    equivalence (y, xyz(1,2))
    equivalence (z, xyz(1,3))

    x = 1
    y = 2
    z = 3

! and we can use just normal array syntax
    print *, xyz(3,:)
end program

This will, however, not work with allocatable arrays. If it's about accessing vectors of matrix it can be easily implemented via pointers. But how can one combine vectors in a two dimensional array? Till now I came only to array of pointers which has problems:

program test_arofpntrs
implicit none

integer :: i
integer, allocatable, target :: xs(:), ys(:), zs(:)

type t_p_xs
  integer, pointer :: p_xs(:)
end type t_p_xs

type(t_p_xs), allocatable :: coords(:)

allocate(coords(3), xs(10), ys(10), zs(10))

xs = 1
ys = 2
zs = 3


coords(1) % p_xs => xs
coords(2) % p_xs => ys
coords(3) % p_xs => zs

print *, coords(1) % p_xs(:)

! this fails:
!print *, coords(:) % p_xs(1)


end program

This is ugly and one cannot access xs(i), ys(i), zs(i). Is it possible to do what I want without copies?

Was it helpful?

Solution

This will not be possible, if you start with separate 1D arrays. The allocatable arrays may be anywhere in the memory. Although Fortran arrays don't have to be contiguous, there has to be some system of strides.

! this fails:
!print *, coords(:) % p_xs(1)

is prohibited by the Fortran standard, because one can not calculate the address of the next element simply. The 1D arrays are even not guaranteed to have the same length.

Also reshape as such does not have to be inefficient, it may just syntactically help with indexing, but don't touch the data at all.

Pointers are a great tool and may help here. You would have to use your 1D arrays bit differently. For example allocate a long 1D array and have 1D pointers for parts of it and a 2D pointer as a whole, or better the other way round:

real,allocatable,target :: xyz(:,:)
real,pointer :: x(:),y(:),z(:)


allocate(xyz(1:10,1:3))
x => xyz(:,1)
y => xyz(:,2)
z => xyz(:,3)

Even the other order of indexes is possible, i.e., xyz(3,10); x => xyz(1,:).

Also you can do

long1D(1:size(xyz)) => xyz

but note it is a Fortran 2008 feature in this direction (2003 otherwise).

Licensed under: CC-BY-SA with attribution
Not affiliated with StackOverflow
scroll top