Allocate on assignment to grow an array.

I have been experimenting with allocate on assigment and the following code uses this to to increment the dimensions of an array:

program p
implicit none
real, allocatable :: a(:)
integer i

do i = 1, 4
  call inc_alloc(a,real(i))
  print*, size(a), a
end do

contains

subroutine inc_alloc (a,next)
real, allocatable, intent(inout) :: a(:)
real, intent(in) :: next
if (allocated(a)) then
  a = [a,next]
else
  a = [next]
end if
end subroutine inc_alloc

end program p

With the 64 bit compiler (release and debug) the terminal output is correct:

                     1     1.00000
                     2     1.00000         2.00000
                     3     1.00000         2.00000         3.00000
                     4     1.00000         2.00000         3.00000         4.00000

With the 32 bit compiler (release and debug) the terminal output is incorrect:

            1     1.00000
            2     0.00000         0.00000
            3     0.00000         0.00000         0.00000
            4     0.00000         0.00000         0.00000         0.00000

With both the 64 and 32 bit compilers (using checkmate) a run time error occurs (nonconformant arrays) at the line:

a = [a,next]

while the arrays are nonconformant, I don't think this should be an error in this case.

Ken

Thank you for the bug report.

I am guessing that the checkmate failure may have to be tollerated in this context. You can use /inhibit_check 19 20 for this purpose. This inhibits both 19 and 20 which provide bounds checking.

Ken

Here is a temporary work-around for Win32:

program p
implicit none
real, allocatable :: a(:)
integer i

do i = 1, 4
  call inc_alloc(a,real(i))
  print*, size(a), a
end do

contains

subroutine inc_alloc (a,next)
real, allocatable, intent(inout) :: a(:)
real tmp(size(a)+1)
real, intent(in) :: next
if (allocated(a)) then
!>>>> a = [a,next]
  tmp = [a,next]
  a = tmp 
else
  a = [next]
end if
end subroutine inc_alloc

end program p

Thanks Paul,

In addition using a temporary does not cause the checkmate failure, which is a big plus 😃

As I don't want to inhibit bounds checking (with checkmate) I think I just need to avoid the a = [a,next] construction.

This allocate on assignement construct turns out to be very difficult to handle for Win32. This is because of the way in which the memory is allocated for the internal temporary array. With the existing method, this memory has already been released at the critical point where the assignment takes place.

For the time being FTN95 will report a compilation failure in this context.

Paul,

Thanks for the example. I have been interested in using auto allocate to extend an allocatable array. This approach requires that there is sufficient free memory for the temporary array, ie memory for twice the new array size.

I tried to learn from your example and made the vector extension size 'chunk(n)' larger to test n = 11 initially worked well in my example.

When n = 1000001 ( 1+10^6 ), it did not work with FTN95 (in Plato) it could be stack size or 'chunk = [ ( i, i=1,n ) ]'

I then changed chunk, tmp to be allocatable, so not overflow the stack, but this still doesn't work with FTN95. it works in Gfortran, but not FTN95 /64

My next attempt would be to have a(:,:) and chunk(6,n) to see what syntax is needed for 2d arrays to be extended.

Any thoughts ?

program p
  implicit none
  integer, parameter :: n = 1000001
  real, allocatable :: a(:)
  real, allocatable :: chunk(:)
  real, allocatable :: larger(:)
  integer :: i,k

  write (*,10) 'step = 1'
 10 format ( /a,i0 )
  do i = 1, 4
    call inc_alloc ( a, [real(i)] )
    print*, size(a), a
  end do
  deallocate ( a )

  chunk = [ ( i, i=1,n ) ]
  write (*,10) 'step = ',size(chunk)
  
  do i = 1, 4
    call inc_alloc ( a, chunk )
    print*, size(a), (a(k),k=1,size(a),size(chunk)/2)
  end do
  deallocate ( a )

  larger = [ ( i, i=1,size(chunk)+50 ) ]
  write (*,10) 'step = ',size(larger)
  
  do i = 1, 4
    call inc_alloc ( a, larger )
    print*, size(a), (a(k),k=1,size(a),size(larger)/2)
  end do
  deallocate ( a )

contains

  subroutine inc_alloc (a, next )
    real, allocatable, intent(inout) :: a(:)
    real, intent(in) :: next(:)
!z    real :: tmp( size(a)+size(next) )
    real, allocatable :: tmp(:)       ! this is safer for large vectors

    if ( allocated(a) ) then
    !>>>> a = [a,next]
      tmp = [ a, next ]
      a = tmp
      write (*,*) 'extend a',size(a)
    else
      a = next
      write (*,*) 'new a',size(a)
    end if

  end subroutine inc_alloc

end program p

Update: if I change line 3 to : integer, parameter :: n = 1022-50 ! 1000001 line 26 to : larger = [ ( i, i=1,n+50 ) ] then the program works,

but n+50 >= 1024 it does not work in Plato ( FTN95 x64 Debug )

This is a lot smaller than what is expected for stack overflow

Also if n = 1024, program fails at 17 : chunk = [ ( i, i=1,n ) ]

( still using FTN95 v8.97.2 )

John

There are lots of different issues here and much depends on whether it's 32 bits or 64 bits.

If you want to make the best use of available memory then avoid 'allocate on assignment' and don't use a construction that forces the compiler to provide a temporary array.

Basically do everything using explicit calls to ALLOCATE. This results in calls to HeapAlloc.

In theory, calling HeapReAlloc might avoid the copying process and we could consider if there is a simple way to extend ALLOCATE so that it calls HeapReAlloc when the user wants to grow the memory allocation.

With this in mind, there might be a way to use (or extend the use of) the SOURCE keyword for ALLOCATE.

Paul,

I've opened another can of Ken's worms!

For the original problem, I wondered about the possibility of using the 2003 intrinsic MOVE_ALLOC. Unfortunately, I cannot get this to work, having tried variations on the following which runs OK with both gFortran and Intel.

program p
implicit none
real, allocatable :: a(:)
integer i

do i = 1, 4
  call inc_alloc(a,real(i))
  print*, size(a), a
end do

contains

subroutine inc_alloc (a,next)
real, allocatable, intent(inout) :: a(:)
real, allocatable :: tmp(:)
real, intent(in) :: next
integer n_new

if (allocated(a)) then
  allocate(tmp(1:size(a)+1))
  tmp(1:size(a)) = a
  tmp(size(a)+1:) = next
  
  call move_alloc(tmp,a)  ! moves the allocation from from TMP to A. 
                          ! TMP will become deallocated in the process.
                          ! FTN complier says: error 1267
                          ! The arguments of MOVE_ALLOC must be
                          ! ALLOCATABLE and have the same type.
                          ! But this code is OK with gFortran and INTEL
else
  a = [next]
end if
end subroutine inc_alloc

end program p

Ken,

The following uses allocate and deallocate exp-licitly and as Paul suggests, hopefully avoids 'allocate on assignment' and doesn't use a construction that forces the compiler to provide a temporary array.

program p2
  implicit none
  real, allocatable :: a(:)
  integer i

  do i = 1, 4
    call inc_alloc ( a, [ real(i),real(i+1) ] )
    print*, size(a), a
  end do

contains

  subroutine inc_alloc ( a, next )
    real, allocatable, intent(inout) :: a(:)
    real, intent(in) :: next(:)
    real, allocatable :: tmp(:)
    integer na, nn

    nn = size(next)
    if ( allocated(a) ) then
       na = size(a)
       write (*,fmt='(a,i0,a,i0)') 'Extend A from ',na,' to ',na+nn
       allocate ( tmp(na+nn) )
       tmp(1:na) = a
       tmp(na+1:) = next
 
       deallocate ( a )
       allocate ( a(na+nn) )
       a = tmp
       deallocate ( tmp )
    else
       write (*,fmt='(a,i0)') 'Create A as ',nn
       allocate ( a(nn) )
       a = next
    end if

  end subroutine inc_alloc

end program p2

It may be a bit more economical (in terms of both memory used and number of bytes copied) to do, instead,

...
    if ( allocated(a) ) then
       na = size(a)
       write (*,fmt='(a,i0,a,i0)') 'Extend A from ',na,' to ',na+nn
       allocate ( tmp(na) )
       tmp(1:na) = a

       deallocate ( a )
       allocate ( a(na+nn) )
       a(1:na) = tmp
       a(na+1:) = next

       deallocate ( tmp )
    else
...

You can see that the contents of 'next' are copied only once, and to the final destination.

The initial issues raised by Ken have now been fixed for the next release of FTN95. His sample program will work for both 64 and 32 bits and runtime bounds checking (e.g. via checkmate) will not cause the program to fail at runtime.

An initial test of MOVE_ALLOC suggests that the FTN95 compile time type checking for the arguments is too strong but /ignore 1267 provides a work-around.

Pail,

My simple auto-allocate code 'a = [ (i,i=1,n) ]' is now not working at all in Plato with /64 /debug (but ok with Gfortran)

 use iso_fortran_env
      real, allocatable :: a(:)
      integer, allocatable :: sizes(:)
      integer n,i,k

      write (*,*) compiler_version ()
      sizes = [ 10,100,1000,1023,1024,1025 ]
      do k = 1,size(sizes)
        n = sizes(k)
        write (*,11) ' '
        write (*,11) 'attempting array size = ', n
!        allocate ( a(n) )
        a = [ (i,i=1,n) ]
        write (*,11) 'resulting array size  = ', size(a)
!        deallocate ( a )
      end do
   11 format (a,i0)
      end 

mecej4,

Thanks, I do think your example is the best approach, as it removes auto-allocate. My final approach was to transfer the size of next and zero the extension. alternatives could be:

  subroutine extend_alloc ( a, nn )     ! with no auto-allocate

!  extend vector A by nn without auto-allocate
  
    real, allocatable, intent(inout) :: a(:)
    integer, intent(in) :: nn
    real, allocatable :: copy(:)
    integer na

    if ( allocated (a) ) then
       na = size(a)
       write (*,fmt='(a,i0,a,i0)') 'Extend A from ',na,' to ',na+nn
       allocate ( copy(na) )
       copy = a
       deallocate ( a )
       allocate ( a(na+nn) )
       a(1:na)  = copy
       a(na+1:) = 0
       deallocate ( copy )
    else
       write (*,fmt='(a,i0)') 'Create A as ',nn
       allocate ( a(nn) )
       a = 0
    end if

  end subroutine extend_alloc

  subroutine ext_alloc ( a, nn )     ! using auto-allocate

!  extend vector A by nn using auto-allocate

    real, allocatable, intent(inout) :: a(:)
    integer, intent(in) :: nn
    real, allocatable :: extra(:)
    integer na

    allocate ( extra(nn) )
    extra = 0     !  or = [ (0,i=1,nn) ]
    if ( allocated (a) ) then
       na = size(a)
       write (*,fmt='(a,i0,a,i0)') 'Extend A from ',na,' to ',na+nn
       a = [ a, extra ]
    else
       write (*,fmt='(a,i0)') 'Create A as ',nn
       a = extra
    end if

  end subroutine ext_alloc

Saturday devilry !! I put the full test back together and it works without a problem !!

program extend
  use iso_fortran_env
  implicit none
  real, allocatable :: a(:)
  integer i

  write (*,*) compiler_version ()

  do i = 1, 4
    select case ( mod(i,2) )
      case (0)
        call auto_alloc ( a, 5000 )
      case (1)
        call extend_alloc ( a, 500 )
    end select
    print*, size(a), a(size(a))
  end do

contains

  subroutine extend_alloc ( a, nn )     ! with no auto-allocate

!  extend vector A by nn without auto-allocate
  
    real, allocatable, intent(inout) :: a(:)
    integer, intent(in) :: nn
    real, allocatable :: copy(:)
    integer na

    if ( allocated (a) ) then
       na = size(a)
       write (*,fmt='(a,i0,a,i0)') 'Extend A from ',na,' to ',na+nn
       allocate ( copy(na) )
       copy = a
       deallocate ( a )
       allocate ( a(na+nn) )
       a(1:na)  = copy
       a(na+1:) = 0
       deallocate ( copy )
    else
       write (*,fmt='(a,i0)') 'Create A as ',nn
       allocate ( a(nn) )
       a = 0
    end if

  end subroutine extend_alloc

  subroutine auto_alloc ( a, nn )     ! using auto-allocate

!  extend vector A by nn using auto-allocate

    real, allocatable, intent(inout) :: a(:)
    integer, intent(in) :: nn
    real, allocatable :: extra(:)
    integer na

!    allocate ( extra(nn) )
    extra = [ (0.0, i=1,nn) ]
    if ( allocated (a) ) then
       na = size(a)
       write (*,fmt='(a,i0,a,i0)') 'Auto-Extend A from ',na,' to ',na+nn
       a = [ a, extra ]
    else
       write (*,fmt='(a,i0)') 'Auto-Create A as ',nn
       a = extra
    end if

  end subroutine auto_alloc

end program extend

However, the 'simple auto-allocate' code in the previous post still fails ?

There is a problem with :

• auto allocate, or
• implied do in array fill or
• both where the array has not been allocated memory ?

perhaps ' a = [ (i,i=1,n) ] ' fails for some values of n or more probably a problem with the order of the statements.

John, Mecej4, Paul, thanks for the suggestions, improvements and feedback on my queries.

John, the following simple example only runs as expected with checkmate for both win32 and x64, so I suspect there is an issue with allocate on assignment for a = [ (i,i=1,n) ]. All other combinations of 32 vs 64 and debug/release generate an exception report at this point.

The error I am seeing may already have been fixed. Paul do you see a problem with this example?

real, allocatable :: a(:)
integer n, i
n = 10 
a = [ (i,i=1,n) ]
print*, allocated(a)
end

Ken

It fails for me so I have added it to the list.

This runs correctly:

real, allocatable :: a(:)
integer n, i
n = 10
a = [ (real(i),i=1,n) ]
print*, allocated(a)
print*, a
end

FTN95 has been fixed for the next release so that this will work...

real, allocatable :: a(:)
integer n, i
n = 10
a = [ (i,i=1,n) ]
print*, allocated(a)
end

Thanks Paul, I think that will fix John's issue.

Thanks Paul for finding the source of the error.