dbcsr_binary_read Subroutine

   SUBROUTINE dbcsr_binary_read(filepath, distribution, matrix_new)
      !! Reads a DBCSR matrix from a file

      IMPLICIT NONE

      CHARACTER(len=*), INTENT(IN)                :: filepath
         !! path to the file
      TYPE(dbcsr_distribution_obj), INTENT(IN)     :: distribution
         !! row and column distribution
      TYPE(dbcsr_type), INTENT(INOUT)               :: matrix_new
         !! DBCSR matrix

      CHARACTER(LEN=*), PARAMETER :: routineN = 'dbcsr_binary_read'

      INTEGER                               :: nblkrows_total, nblkcols_total, &
                                               nblks, darea_size, data_type, type_size, &
                                               globalinfo_size, &
                                               size_of_pgrid, &
                                               i, j, &
                                               nblocks, &
                                               share_size, order, cur_blks, &
                                               job_count, start_index, end_index, &
                                               localinfo_length, blockdata_length, &
                                               worker_id, group_list_size, handle, linfo_length
      CHARACTER                             :: matrix_type
      CHARACTER(LEN=default_string_length)  :: matrix_name
      INTEGER, PARAMETER                    :: version_len = 10
      CHARACTER(LEN=version_len)            :: version
      CHARACTER(LEN=80)                     :: matrix_name_v_1_0
      CHARACTER(LEN=version_len), PARAMETER :: version_v_1_0 = "DBCSRv_1.0"
      TYPE(mp_comm_type)                    :: group_id
      TYPE(mp_file_type)                    :: thefile

      INTEGER, DIMENSION(:), POINTER        :: row_p, col_i, blk_p, &
                                               proc_nblks, proc_darea_sizes
      INTEGER, DIMENSION(4)                 :: values
      INTEGER, ALLOCATABLE, DIMENSION(:)    :: linfo_lens, bdata_lens
      INTEGER, ALLOCATABLE, DIMENSION(:), TARGET     :: ginfo_vec, linfo_vec, &
                                                        rowp, coli, blkp
      INTEGER, ALLOCATABLE, DIMENSION(:, :), TARGET   :: val_data
      INTEGER, DIMENSION(:), POINTER, CONTIGUOUS     :: row_blk_size, col_blk_size
      TYPE(dbcsr_mp_obj)                             :: mp_env
      TYPE(dbcsr_data_obj)                           :: dblk
      REAL(sp)                                       :: rsp_dummy(1)
      REAL(dp)                                       :: rdp_dummy(1)
      COMPLEX(sp)                                    :: csp_dummy(1)
      COMPLEX(dp)                                    :: cdp_dummy(1)
      REAL(sp), ALLOCATABLE, DIMENSION(:), TARGET     :: rsp
      REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET     :: rdp
      COMPLEX(sp), ALLOCATABLE, DIMENSION(:), TARGET  :: csp
      COMPLEX(dp), ALLOCATABLE, DIMENSION(:), TARGET  :: cdp
      INTEGER(kind=file_offset), ALLOCATABLE, DIMENSION(:)   :: subh2_offsets, &
                                                                subh3_disps, subh3_offsets, &
                                                                bdata_disps, bdata_offsets
      INTEGER(kind=file_offset), PARAMETER    :: BOF = 0
      INTEGER(kind=file_offset)               :: offset, subh2_start, subh3_start, bdata_start, file_size, &
                                                 localinfo_offset, blockdata_offset, sum_nblks, subh3_length, data_area_size
      INTEGER, PARAMETER                      :: char_count = 1 + version_len + default_string_length

      CALL timeset(routineN, handle)

      mp_env = dbcsr_distribution_mp(distribution)
      group_id = dbcsr_mp_group(mp_env)

      CALL mp_environ(group_list_size, worker_id, group_id)
      CALL mp_file_open(group_id, thefile, filepath, file_amode_rdonly)

! read version, matrix name and matrix type
      CALL mp_file_read_at_all(thefile, BOF, version)

      IF (version /= version_v_1_0) &
         DBCSR_WARN("Trying to read an unknown version of the matrix data file. Good luck!")

      CALL mp_file_read_at_all(thefile, BOF + version_len*mpi_character_size, matrix_name_v_1_0)
      matrix_name = matrix_name_v_1_0

      CALL mp_file_read_at_all(thefile, BOF + (version_len + default_string_length)*mpi_character_size, matrix_type)
! read 4 integer values form sub-header1
      CALL mp_file_read_at_all(thefile, BOF + char_count*mpi_character_size, values)
      size_of_pgrid = values(1)
      data_type = values(2)
      nblkrows_total = values(3)
      nblkcols_total = values(4)
! read 2 vectors, row_blk_size and col_blk_size, from sub-header1
      globalinfo_size = nblkrows_total + nblkcols_total
      ALLOCATE (ginfo_vec(globalinfo_size))
      CALL mp_file_read_at_all(thefile, BOF + char_count*mpi_character_size + 4*mpi_integer_size, ginfo_vec)
      row_blk_size => ginfo_vec(1:nblkrows_total)
      col_blk_size => ginfo_vec(nblkrows_total + 1:globalinfo_size)

! compute the offsets where sub-header2 and sub-header3 start
      subh2_start = (4 + globalinfo_size)*mpi_integer_size + char_count*mpi_character_size
      subh3_start = subh2_start + 2*size_of_pgrid*mpi_integer_size

! compute the offsets in sub-header2 and read 2 integers nblocks, data_area_size
      ! number of data chunks from sub-header 2 and 3 to be read by every node rounded up
      ! to the next integer to make it even for all the nodes in the specified mpi group
      share_size = CEILING(REAL(size_of_pgrid, KIND=dp)/group_list_size)

      ALLOCATE (subh2_offsets(share_size))
      subh2_offsets = BOF
      DO i = 1, share_size
         offset = subh2_start + mpi_integer_size*2*(worker_id + (i - 1)*group_list_size)
         IF (offset .GE. subh3_start) EXIT
         subh2_offsets(i) = offset
      END DO

      ALLOCATE (val_data(3, share_size))
      val_data(:, :) = 0
      DO i = 1, share_size
         CALL mp_file_read_at_all(thefile, subh2_offsets(i), values, msglen=2)
         nblocks = values(1)
         data_area_size = values(2)
         IF (subh2_offsets(i) .EQ. 0) EXIT
         val_data(1, i) = nblocks
         IF (data_area_size >= HUGE(val_data(2, i))) &
            DBCSR_ABORT("Data area too large, fix code.")
         val_data(2, i) = INT(data_area_size)
         val_data(3, i) = worker_id + (i - 1)*group_list_size + 1 ! order
         ! order = indices of an array of length size_of_pgrid to be accessed by the current node
      END DO
      nblks = SUM(val_data(1, :))
      darea_size = SUM(val_data(2, :))
      proc_nblks => val_data(1, :) ! to be passed to dbcsr_datablock_redistribute
      proc_darea_sizes => val_data(2, :) ! to be passed to dbcsr_datablock_redistribute

! compute the offsets in sub-header3 and read 3 vectors row_p, col_i, blk_p
      ! actual number of chunks to be read by the current node
      job_count = COUNT(val_data(3, :) .NE. 0)
      CALL mp_file_get_size(thefile, file_size)

      ALLOCATE (linfo_lens(size_of_pgrid))
      ALLOCATE (subh3_disps(size_of_pgrid))
      ALLOCATE (subh3_offsets(size_of_pgrid))
      linfo_lens = 0; subh3_disps = 0
      DO i = 1, size_of_pgrid
         DO j = 1, share_size
            order = val_data(3, j)
            IF (i .EQ. order) linfo_lens(order) = &
               1 + nblkrows_total + 2*val_data(1, j)
         END DO
      END DO
      CALL mp_sum(linfo_lens, group_id)
      CALL cumsum_l(INT(linfo_lens, kind=file_offset), subh3_disps)
      subh3_disps(:) = CSHIFT(subh3_disps, shift=-1)
      subh3_disps(1) = BOF
      subh3_offsets(:) = subh3_start + subh3_disps*mpi_integer_size

      sum_nblks = INT(nblks, kind=file_offset)
      CALL mp_sum(sum_nblks, group_id)
      subh3_length = size_of_pgrid*INT(1 + nblkrows_total, KIND=file_offset) + 2*sum_nblks

      linfo_length = nblkrows_total + 1 + 2*MAXVAL(val_data(1, :))

      ALLOCATE (linfo_vec(linfo_length))
      ALLOCATE (rowp((nblkrows_total + 1)*job_count))
      ALLOCATE (coli(nblks))
      ALLOCATE (blkp(nblks))
      DO i = 1, share_size
         order = val_data(3, i)
         cur_blks = val_data(1, i)
         IF (order .EQ. 0) THEN
            localinfo_offset = file_size
            localinfo_length = 0
         ELSE
            localinfo_offset = subh3_offsets(order)
            localinfo_length = linfo_lens(order)
         END IF
         CALL mp_file_read_at_all(thefile, localinfo_offset, linfo_vec, msglen=localinfo_length)
         IF (localinfo_length .EQ. 0) EXIT

         rowp((i - 1)*(nblkrows_total + 1) + 1:i*(nblkrows_total + 1)) = linfo_vec(1:nblkrows_total + 1)
         start_index = SUM(val_data(1, 1:i - 1)) + 1
         end_index = SUM(val_data(1, 1:i))
         coli(start_index:end_index) = &
            linfo_vec(nblkrows_total + 2:cur_blks + nblkrows_total + 1)
         blkp(start_index:end_index) = &
            linfo_vec(cur_blks + nblkrows_total + 2:2*cur_blks + nblkrows_total + 1)
      END DO
      row_p => rowp
      col_i => coli
      blk_p => blkp

! compute the offsets and read block data
      ALLOCATE (bdata_lens(size_of_pgrid))
      ALLOCATE (bdata_disps(size_of_pgrid))
      ALLOCATE (bdata_offsets(size_of_pgrid))
      bdata_lens = 0
      DO i = 1, size_of_pgrid
         DO j = 1, share_size
            order = val_data(3, j)
            IF (i .EQ. order) bdata_lens(order) = val_data(2, j)
         END DO
      END DO
      CALL mp_sum(bdata_lens, group_id)
      CALL cumsum_l(INT(bdata_lens, kind=file_offset), bdata_disps)
      bdata_disps(:) = CSHIFT(bdata_disps, shift=-1)
      bdata_disps(1) = BOF

      bdata_start = subh3_start + subh3_length*mpi_integer_size
      SELECT CASE (data_type)
      CASE (dbcsr_type_real_4)
         type_size = real_4_size
      CASE (dbcsr_type_real_8)
         type_size = real_8_size
      CASE (dbcsr_type_complex_4)
         type_size = 2*real_4_size
      CASE (dbcsr_type_complex_8)
         type_size = 2*real_8_size
      END SELECT
      bdata_offsets(:) = bdata_start + bdata_disps*type_size

      SELECT CASE (data_type)
      CASE (dbcsr_type_real_4)
         ALLOCATE (rsp(darea_size))
         DO i = 1, share_size
            order = val_data(3, i)
            ! use dummy one-sized data array as buffer in place of empty array
            ! when nothing is supposed to be read (order = 0)
            IF (order .EQ. 0) THEN
               blockdata_offset = file_size
               CALL mp_file_read_at_all(thefile, blockdata_offset, rsp_dummy)
            ELSE
               start_index = SUM(val_data(2, 1:i - 1)) + 1
               end_index = SUM(val_data(2, 1:i))
               blockdata_length = bdata_lens(order)
               blockdata_offset = bdata_offsets(order)
               CALL mp_file_read_at_all(thefile, blockdata_offset, rsp(start_index:end_index), &
                                        msglen=blockdata_length)
            END IF
         END DO
      CASE (dbcsr_type_real_8)
         ALLOCATE (rdp(darea_size))
         DO i = 1, share_size
            order = val_data(3, i)
            IF (order .EQ. 0) THEN
               blockdata_offset = file_size
               CALL mp_file_read_at_all(thefile, blockdata_offset, rdp_dummy)
            ELSE
               start_index = SUM(val_data(2, 1:i - 1)) + 1
               end_index = SUM(val_data(2, 1:i))
               blockdata_length = bdata_lens(order)
               blockdata_offset = bdata_offsets(order)
               CALL mp_file_read_at_all(thefile, blockdata_offset, rdp(start_index:end_index), &
                                        msglen=blockdata_length)
            END IF
         END DO
      CASE (dbcsr_type_complex_4)
         ALLOCATE (csp(darea_size))
         DO i = 1, share_size
            order = val_data(3, i)
            IF (order .EQ. 0) THEN
               blockdata_offset = file_size
               CALL mp_file_read_at_all(thefile, blockdata_offset, csp_dummy)
            ELSE
               start_index = SUM(val_data(2, 1:i - 1)) + 1
               end_index = SUM(val_data(2, 1:i))
               blockdata_length = bdata_lens(order)
               blockdata_offset = bdata_offsets(order)
               CALL mp_file_read_at_all(thefile, blockdata_offset, csp(start_index:end_index), &
                                        msglen=blockdata_length)
            END IF
         END DO
      CASE (dbcsr_type_complex_8)
         ALLOCATE (cdp(darea_size))
         DO i = 1, share_size
            order = val_data(3, i)
            IF (order .EQ. 0) THEN
               blockdata_offset = file_size
               CALL mp_file_read_at_all(thefile, blockdata_offset, cdp_dummy)
            ELSE
               start_index = SUM(val_data(2, 1:i - 1)) + 1
               end_index = SUM(val_data(2, 1:i))
               blockdata_length = bdata_lens(order)
               blockdata_offset = bdata_offsets(order)
               CALL mp_file_read_at_all(thefile, blockdata_offset, cdp(start_index:end_index), &
                                        msglen=blockdata_length)
            END IF
         END DO
      END SELECT
      CALL dbcsr_data_init(dblk)
      CALL dbcsr_data_new(dblk, data_type)
      IF (ALLOCATED(rdp)) dblk%d%r_dp => rdp
      IF (ALLOCATED(rsp)) dblk%d%r_sp => rsp
      IF (ALLOCATED(cdp)) dblk%d%c_dp => cdp
      IF (ALLOCATED(csp)) dblk%d%c_sp => csp

      CALL mp_file_close(thefile)

      CALL dbcsr_create(matrix_new, matrix_name, distribution, matrix_type, &
                        row_blk_size, col_blk_size, nze=darea_size, &
                        data_type=data_type)
      CALL dbcsr_datablock_redistribute(dblk, row_p, col_i, blk_p, proc_nblks, proc_darea_sizes, matrix_new)

      DEALLOCATE (subh2_offsets, subh3_offsets, bdata_offsets)
      DEALLOCATE (subh3_disps, bdata_disps)
      DEALLOCATE (linfo_lens, bdata_lens)
      DEALLOCATE (val_data, ginfo_vec, linfo_vec)
      DEALLOCATE (rowp, coli, blkp)
      IF (ALLOCATED(rdp)) DEALLOCATE (rdp)
      IF (ALLOCATED(rsp)) DEALLOCATE (rsp)
      IF (ALLOCATED(cdp)) DEALLOCATE (cdp)
      IF (ALLOCATED(csp)) DEALLOCATE (csp)
      CALL dbcsr_data_clear_pointer(dblk)
      DEALLOCATE (dblk%d)

      CALL timestop(handle)
   CONTAINS
      SUBROUTINE cumsum_l(arr, cumsum)
         INTEGER(kind=file_offset), DIMENSION(:), &
            INTENT(IN)                                      :: arr
         INTEGER(kind=file_offset), DIMENSION(SIZE(arr)), &
            INTENT(OUT)                                     :: cumsum

         INTEGER                                            :: i

         cumsum(1) = arr(1)
         DO i = 2, SIZE(arr)
            cumsum(i) = cumsum(i - 1) + arr(i)
         END DO
      END SUBROUTINE cumsum_l

   END SUBROUTINE dbcsr_binary_read