Changing size of an array within a loop. Good practice?

I am quite new to Fortran, so excuse me if you find my doubt a bit naïve.
I need to implement a loop (do) in which, at each instance, a vector array changes its length. I have implemented this by using an allocatable integer array wich is allocated once we enter into the loop to its corresponding length, and then, deallocated after finishing the work done within the loop.
My question is: Is this a good practice in Fortran?
I am interested in fast performance, and I can think in other ways of doing this, e.g. keeping the size of the array constant, initializing all the elements to -1 (my elements are always positive) and modifying just the first elements needed at each instance. Then, the code would be a bit more complex as it would have to check if the next element is -1 or not, but this should not be a major problem.
Does this higher complexity pay for the fact of avoiding the “allocate”/“deallocate” issue?

It may depend on how many times you do allocate / deallocate.
Ten would be ok, but a million then perhaps not.
Either way, the allocate overhead might not be too onerous.

If you require array syntax within the loop, you can still allocate the array to the largest size required for the loop, then perform each loop operation in a called routine. This can give you the ability to define the used array size for each loop.

  integer :: n,nm
  real, allocatable :: array(:)
  real :: res_n

    nm = 100
    allocate ( array(nm) )
    do n = 1,nm
      call work_on_array (array,n,res_n)
    end do
  end
  
  subroutine work_on_array (array,n,res)
  integer :: i,n
  real    :: array(n), res
  
  ! now have an array size = n to work on 
    array = 0
    do i = 1, size(array)
      array(i) = i
    end do
    res = dot_product ( array, array )
  end subroutine work_on_array

@jmcano,

Welcome to the forum.

If you are interested in performance, it may be better if you can share your code, or a small example that mimics it.

Readers will then be able to point out different ways to author the code, you will find out there are quite a few options.

You may then want evaluate the options via careful instrumentation to test CPU usage and then decide what works best for your needs.

1 Like

Also remember the maxime: Premature optimisation is the root of all evil.

Or in a positive sense: The biggest optimisation step is that you get from a non-working program to a working program. Anything else is mere peanuts compared to that :slight_smile:

But I agree with @FortranFan: show the code (or at least a representative example of what you want to achieve) and we can give better and more dedicated advice.

Would this approach imply that adoption of OpenMP for a working program is peanuts ?

Well, if you add OpenMP the peanuts are possibly gigantic, but I was referring to the usefulness (or performance) going from zilch to a positive value. And you can’t really beat a ratio of infinity, can you? :slight_smile:

Or you can use a block construct:

program changing_array
    implicit none

    integer :: n = 10
    integer :: sz
    integer :: i

    do i = 1, n
        sz = i*(i-1)
        block
            real :: changing(sz)
            integer :: k
            changing = [(k, k=1,sz)]
            print *,size(changing), sum(changing)
        end block
    enddo
end program

Whether or not it is a good idea to reallocate a work array within a loop, it depends on several factors. How much effort is done with the array each cycle of the loop? If there is a lot of effort, then the allocation/deallocation effort will be insignificant. If there is only a little work, say something like a dot product as shown in one of the above examples, then the dot product effort will be swamped by the allocation/deallocation effort. And there can be anything in between those two extremes.

My programming dates back to pre f77 days where allocate/deallocate were not available, so it was typical to dimension arrays to their maximum size and then use just a subset of the elements. Thus if I had this kind of problem, where I knew ahead of time the max array size, I would just allocate it to the max size and use the appropriate subset within each do loop iteration. That can be done several ways now with modern fortran: you can use (1:n) array slice notation, you can pass the array slice to a subroutine (perhaps even an internal subroutine) and then use whole array notation within the subroutine, or you can use an ASSOCIATE block inline in your code and use whole array notation within that block.

2 Likes

You may want to use a fixed size array instead of iteratively reallocating it if both these two conditions are met:

  • this part of your program represent a significant part of the total runtime
  • there is little calculation in the loop, so that the allocations/deallocations represent a significant part of this part (as mentioned by @RonShepard )
1 Like

Have you thought of using a linked list?

You’d be able to pass in the list at what ever point you want and the number of elements you want your code to work on. When the subroutine returns, just make sure that the first and last elements of your sub list are linked back into the larger list. The total number of links could even change, it does not matter!

type t_node, pointer :: p_root, p_cuta , p_cutb, p_cutc, p_cutd
… code that identifies node BEFORE cut as p_cuta and AFTER as p_cutd …
p_cutb => p_cuta&p_next
p_cutc => p_cutd&p_prev
p_cutb&p_prev => null()
p_cutc&p_next => null()
call datasubroutine(p_cutb, p_cutc)
p_cuta&p_next => p_cutb
p_cutb&p_prev => p_cuta
p_cutc&p_next => p_cutd
p_cutd&p_prev => p_cutc

Maybe give that a try?

Knarfnarf

Edit: need to ensure that the outer nodes are connected properly!

Linked lists can be powerful tools, but generally not pertinent for performances compared to arrays. And the OP stated that he was interested in “fast performances”.

The truth is that very often fixed size arrays can do, at least for rank 1 (this can be an allocatable, but allocated only once), even when the needed size varies along the execution. We tend to use the reallocation on assignement feature mainly because it’s more elegant.

A mid-term solution when one doesn’t want to allocate an unneeded oversized array is to implement a kind of smart reallocation where the physical size is always a power of 2, aiming at minimizing the number of reallocations while avoiding too much oversized arrays:

type smart_reallocatable_int
    integer, allocatable :: a(:)
    integer :: n = 0    ! requested number of elements
end type

subroutine smart_allocate(x,n,keep)
    type(smart_reallocatable_int) :: x
    integer, intent(in) :: n
    logical, intent(in) :: keep

    integer, allocatable :: tmp(:)
    integer :: n2

    if (n == 0) then        
        x%a = [integer::]
    else if (.not.allocated(x%a) .or. size(x%a) == 0) then
        n2 = 1
        do while (n2 < n)
            n2 = 2*n2
        end do
        allocate( x%a(n2) )
    else 
        n2 = size(x%a)
        do while (n2 < n)
            n2 = 2 * n2
        end do
        do while (n2 >= 2*n)
            n2 = n2 / 2
        end do
        if (n2 /= size(x%a)) then
            if (keep) then
                allocate( tmp(n2) )
                tmp(1:min(x%n,n2)) = x%a(1:min(x%n),n2))
                move_alloc(tmp,x%a)
            else
                deallocate( x%a )
                allocate( x%a(n2) )
            end if
            allocate( x%a(n2) )
        end if            
    end if
    x%n = n
end subroutine smart_allocate