unknown size matrix in Fortran

Question

I want to add elements to a 1d matrix mat, subject to a condition as in the test program below. In Fortran 2003 you can add an element

mat=[mat,i]

as mentioned in the related question Fortran array automatically growing when adding a value. Unfortunately, this is very slow for large matrices. So I tried to overcome this, by writing the matrix elements in an unformatted file and reading them afterwards. This turned out to be way faster than using mat=[mat,i]. For example for n=2000000_ilong the run time is 5.1078133666666661 minutes, whereas if you store the matrix elements in the file the run time drops to 3.5234166666666665E-003 minutes.
The problem is that for large matrix sizes the file storage.dat can be hundreds of GB... Any ideas?

program test


implicit none

integer, parameter :: ndig=8
integer, parameter :: ilong=selected_int_kind(ndig)
integer (ilong), allocatable :: mat(:)
integer (ilong), parameter :: n=2000000_ilong
integer (ilong) :: i, cn 
logical, parameter :: store=.false.
real(8) :: z, START_CLOCK, STOP_CLOCK


open(1, file='storage.dat',form='unformatted')

call cpu_time(START_CLOCK)

if(store) then 

 cn=0
 do i=1,n
   call random_number(z)
   if (z<0.5d0) then 
     write(1) i
     cn=cn+1
   end if 
 end do


 rewind(1); allocate(mat(cn)); mat=0

 do i=1,cn
  read(1) mat(i)
 end do


else 

 allocate(mat(1)); mat=0
 do i=1,n
   call random_number(z)
   if (z<0.5d0) then 
    mat=[mat,i]
   end if 
 end do



end if




call cpu_time(STOP_CLOCK)
print *, 'run took:', (STOP_CLOCK - START_CLOCK)/60.0d0, 'minutes.'


end program test

Answer 1

Maybe you need a dynamic container such as C++'s std::vector, with a push_back() function.

The following is a simplified version. You probably ought to check the allocation to make sure that you don't run out of addressable memory.

Note the need for random_seed.

module container
   use iso_fortran_env
   implicit none

   type array
      integer(int64), allocatable :: A(:)
      integer(int64) num
   contains
      procedure push_back
      procedure print
   end type array

   interface array                               ! additional constructors
      procedure array_constructor      
   end interface array

contains

   !----------------------------------------------

   function array_constructor() result( this )   ! performs initial allocation
      type(array) this
      allocate( this%A(1) )
      this%num = 0
   end function array_constructor

   !----------------------------------------------

   subroutine push_back( this, i )
      class(array), intent(inout) :: this
      integer(int64) i
      integer(int64), allocatable :: temp(:)

      if ( size(this%A) == this%num ) then       ! Need to resize
         allocate( temp( 2 * this%num ) )        ! <==== for example
         temp(1:this%num ) = this%A
         call move_alloc( temp, this%A )
!        print *, "Resized to ", size( this%A )  ! debugging only!!!
      end if

      this%num = this%num + 1
      this%A(this%num) = i

   end subroutine push_back

   !----------------------------------------------

   subroutine print( this )
      class(array), intent(in) :: this
      write( *, "( *( i0, 1x ) )" ) ( this%A(1:this%num) )
   end subroutine print

end module container

!=======================================================================

program test
   use iso_fortran_env
   use container
   implicit none
   type(array) mat
   integer(int64) :: n = 2000000_int64
   integer(int64) i
   real(real64) z, START_CLOCK, STOP_CLOCK

   mat = array()                       ! initial trivial allocation

   call random_seed                    ! you probably need this
   call cpu_time(START_CLOCK)
   do i = 1, n
      call random_number( z )
      if ( z < 0.5_real64 ) call mat%push_back( i )
   end do

   call cpu_time(STOP_CLOCK)
   print *, 'Run took ', ( STOP_CLOCK - START_CLOCK ) / 60.0_real64, ' minutes.'

!  call mat%print                      ! debugging only!!!

end program test