Degenerate Conic

Algorithms • Modern Fortran Programming • Orbital Mechanics

Aug 07, 2016

Dynamically Sizing Arrays

Often the need arises to add (or subtract) elements from an array on the fly. Fortran 2008 allows for this to be easily done using standard allocatable arrays. An example for integer arrays is shown here:

integer,dimension(:),allocatable :: x

x = [1,2,3]
x = [x,[4,5,6]] ! x is now [1,2,3,4,5,6]
x = x(1:4) ! x is now [1,2,3,4]

Note that, if using the Intel compiler, this behavior is not enabled by default for computational efficiency reasons. To enable it you have to use the -assume realloc_lhs compiler flag.

Resizing an array like this carries a performance penalty. When adding a new element, the compiler will likely have to make a temporary copy of the array, deallocate the original and resize it, and then copy over the original elements and the new one. A simple test case is shown here (compiled with gfortran 6.1.0 with -O3 optimization enabled):

program test

implicit none

integer,dimension(:),allocatable :: x
integer :: i

x = [0]
do i=1,100000
    x = [x,i]
end do

end program test

This requires 2.828986 seconds on my laptop (or 35,348 assignments per second). Now, that may be good enough for some applications. However, performance can be improved significantly by allocating the array in chunks, as shown in the following example, where we allocate in chunks of 100 elements, and then resize it to the correct size at the end:

program test

    implicit none

    integer, dimension(:), allocatable :: x
    integer :: i, n

    integer, parameter :: chunk_size = 100

    n = 0
    do i = 0, 100000
        call add_to(x, i, n, chunk_size, finished=i == 100000)
    end do

contains

    pure subroutine add_to(vec, val, n, chunk_size, finished)

        implicit none

        integer, dimension(:), allocatable, intent(inout) :: vec
            !! the vector to add to
        integer, intent(in) :: val
            !! the value to add
        integer, intent(inout) :: n
            !! counter for last element added to vec.
            !! must be initialized to size(vec)
            !! (or 0 if not allocated) before first call
        integer, intent(in) :: chunk_size
            !! allocate vec in blocks of this size (>0)
        logical, intent(in) :: finished
            !! set to true to return vec
            !! as its correct size (n)

        integer, dimension(:), allocatable :: tmp

        if (allocated(vec)) then
            if (n == size(vec)) then
                ! have to add another chunk:
                allocate (tmp(size(vec) + chunk_size))
                tmp(1:size(vec)) = vec
                call move_alloc(tmp, vec)
            end if
            n = n + 1
        else
            ! the first element:
            allocate (vec(chunk_size))
            n = 1
        end if

        vec(n) = val

        if (finished) then
            ! set vec to actual size (n):
            if (allocated(tmp)) deallocate (tmp)
            allocate (tmp(n))
            tmp = vec(1:n)
            call move_alloc(tmp, vec)
        end if

    end subroutine add_to

end program test

This requires only 0.022938 seconds (or 4,359,577 assignments per second) which is nearly 123 times faster. Note that we are using the Fortran 2003 move_alloc intrinsic function, which saves us an extra copy operation when the array is resized.

Increasing the chunk size can improve performance even more:

results_cases_per_sec

Depending on the specific application, a linked list is another option for dynamically-sized objects.