Tree structure using derived type - Slow allocation and recovery

[Arash_Z_]

Hi,

I am trying to read and store data in a tree structure using derived types to save memory. My code looks like something like shown below but here I just store a constant value and recover it for simplicity. The array dimensions in my actual problem are so large that I had to use jagged arrays using derived types to mange storage. The issue is that the code becomes very slow. I understand the slowness in the storing part due to multiple allocation since for derived types I had to allocate for each layer. However, the recovery part is also very slow. I was wondering if there is another way around this type of problem to increase the speed. I cannot define a big static array from the beginning since the size will be so large that I get memory errors. Thanks for helps and suggestions in advance. Here is my sample code:

      program arash

 
      integer i,j,k,l,m,n
      real a

          type :: Row5
  real, allocatable :: Row5(:)
          end type Row5
              
              
          type :: Row4
  type(Row5), allocatable  :: Row4(:)
          end type Row4
          
          type :: Row3
  type(Row4), allocatable  :: Row3(:)
          end type Row3
          
          type :: Row2
  type(Row3), allocatable  :: Row2(:)
          end type Row2
          
          type :: Row1
  type(Row2), allocatable  :: Row1(:)
          
          end type Row1
          
          type(Row1), allocatable, dimension(:) :: MainMatrix
 

              
          allocate (MainMatrix(100))
          
C---------------Storing data-----------------------------------   
          do i=1,100
                          if (.not. allocated(MainMatrix(i)%Row1)) allocate (MainMatrix(i)%Row1(8))
           do j=1,8
                           if (.not. allocated(MainMatrix(i)%Row1(j)%Row2)) allocate (MainMatrix(i)%Row1(j)%Row2(54))
            do k=1,54
                            if (.not. allocated(MainMatrix(i)%Row1(j)%Row2(k)%Row3)) 
     1      allocate (MainMatrix(i)%Row1(j)%Row2(k)%Row3(5))
             do l=1,5
                             if (.not. allocated(MainMatrix(i)%Row1(j)%Row2(k)%Row3(l)%Row4))
     1      allocate (MainMatrix(i)%Row1(j)%Row2(k)%Row3(l)%Row4(10))
              do m=1,10

            if (.not. allocated(MainMatrix(i)%Row1(j)%Row2(k)%Row3(l)%Row4(m)%Row5)) 
     1      allocate (MainMatrix(i)%Row1(j)%Row2(k)%Row3(l)%Row4(m)%Row5(21))
            do n=1,21
            MainMatrix(i)%Row1(j)%Row2(k)%Row3(l)%Row4(m)%Row5(n)=3.66
            enddo
              enddo
             enddo
            enddo
           enddo
          enddo
             
          
C---------------Recovering data-----------------------------------      
         do i=1,100
           do j=1,8
            do k=1,54
             do l=1,5
              do m=1,10
            do n=1,21
            a=MainMatrix(i)%Row1(j)%Row2(k)%Row3(l)%Row4(m)%Row5(n)
              enddo
              enddo
             enddo
            enddo
           enddo
         enddo
                    
                    
      deallocate (MainMatrix)
      
      end program arash

[jimdempseyatthecove]

Are your array dimensions static once known?

If so, then make a flat array (allocate single dimension array of size equal to product of individual indexes), then use a function to produce the linear index given the individual indices.

module foo
  real, allocatable :: MainMatrix(:)
  integer :: iCount, jCount, kCount, lCount, mCount
  contains
  subroutine AllocateMainMatrix(i,j,k,l,m)
    integer :: i,j,k,l,m
    allocate(MainMatrix(i*j*k*l*m)
    iCount = i; jCount=j; kCount=k; lCount=l; mCount=m
  end subroutine AllocateMainMatrix
  integer function IndexMainMatrix(i,j,k,l,m)
    integer :: i,j,k,l,m
    IndexMainMatrix = (((((m-1)*lCount + l-1)*kCount) + k-1)*jCount + j-1)*iCount + i
  end function IndexMainMatrix
end module foo

(untested code)

Then: a=MainMatrix(IndexMainMatrix(I,j,k,l,m))

Jim Dempsey

[Arash_Z_]

Thanks for the reply Jim. The array dimensions are not known. It reads the dimension for each branch of each node from the input files. The main issue is that if I try to allocate with the larges dimension for each node, it saturates the memory.

Thanks,

Arash

[FortranFan]

It appears you are looking for a memory-efficient data structure for a sparse matrix.  If that's correct, research this topic and you'll find a lot of good ideas you can implement in Fortran.  Also, look into the first book mentioned in Dr Fortran blog by Hanson and Hopkins, they offer good examples in it that you can review for your needs.

[andrew_4619]

I am not clear what you are trying to do from the example. You seem to have a 5D array of reals. Is it that you have a numbers of nodes and for each node you have a 4D array of reals, where the 4D dimensions vary from node to node? If that is that case it would probably make more sense to have that data in a large 1D array and then have have some integer index array that gives the start position in the real array and the four  dimensions of each node array so you can extract the data as and when needed.

The would eliminate all the allocating and deallocating delays. If you need to dynamically modify the data structure  a more complex system would be needed.

[jimdempseyatthecove]

>The array dimensions are not known. It reads the dimension for each branch of each node from the input files.

I understand that. I'd asked: Are your array dimensions static once known?

At some point of reading the input file(s) you know the(a) dimension(s), else you would not know how much to allocate. Your code snip illustrates that at some point you know the dimensions. If your sketch code is not right, and your dimensions are embedded throughout the data set, then this may require: a) two passes on the input files, b) estimating upper limit of total number of variables based on file size, c) when creating the data, write a header that contains the pertinent information.

Note, the 1-D array concept for the entire dataset can be extended (reduced) to a sub-set of your 5-D array (e.g. 1-D array of 1-D arrays representing a flat space of 4-D array).

Please note, using the current structure (array of arrays of arrays), that is the last (inner most) dimension is small, then the size of the array descriptor could exceed that of the allocation.

Jim Dempsey

[jimdempseyatthecove]

Additional information.

If you compile as 64-bit application, as long as the heap manager does not wipe/prefill the allocation, you can allocate an arbitrarily huge block of address space without consuming RAM nor swap file space. It is only when you first touch (typically write to) the virtual memory page that the RAM and/or page file gets mapped. Debug build might wipe, so run without the debug heap manager.

Jim Dempsey

[Arash_Z_]

Thanks for your responses app4619 and jimdempseyatthecove.

Yes I have a 5 node array and each node have a 4D array (let's call them sub nodes). I can not store to whole structure in a 1D array since the dimension for that 1D array i.e. I*j*k*l*m*n will exceeds the maximum allowable integer size on a 64-bit system.

Jim, yes at some point I know the dimension of arrays for each nodes. I will try to extend your 1D array concept sub-sets as you suggested and post the outcome.

Thanks,

Arash

[jimdempseyatthecove]

>>i.e. i*j*k*l*m*n will exceeds the maximum allowable integer size on a 64-bit system

Use INTEGER(8) variables and(or) at least assure that the expression i*j*k*l*m*n produces an INTEGER(8) result.

integer(8) :: bigSize
integer :: i,j,k,lo,m
...
bigSize = INT8(i) * INT8(j) * INT8(k) * INT8(l) * INT8(m) * INT8(n)
allocate(bigArray(bigSize))
...
 

Also use a similar technique to assure you produce an INTEGER(8) index in the function you write to convert the INTEGER(4) indicies into a flat space index

Note, you could probably get by with bisSize = (((((INT8(i))*j)*k)*l)*m)*n

But this requires that the parenthesis be honored (don't assume they are and will always be).

Jim Dempsey

[andrew_4619]

at #8, If the array index exceeds a 64 bit integer you would have a much bigger problem - the need for more memory than any version of windows currently supports..... > 32 exabytes (?)

Just for completeness -  the default integer is 32 bit even if you are making a 64 bit application. If you want 64 bit integers you explicitly need integer(8) or need to change default compiler options.