read_dense Subroutine

    subroutine read_dense(this, filename)
        class(dense_lMat_t), intent(inout) :: this
        character(*), intent(in) :: filename
        integer :: iunit, ierr
        integer(int64) :: indices(6)
        HElement_t(dp) :: matel
        character(*), parameter :: t_r = "readLMat"
        integer(int64) :: counter, index
        character(1024) :: char1024

        call this%alloc(this%lMat_size())

        if (tHDF5LMat) then
#ifdef USE_HDF_
            call this%read_hdf5_dense(filename)
#else
            call stop_all(t_r, "HDF5 integral files disabled at compile time")
#endif
        else
            if (iProcIndex_intra == 0) then
                iunit = get_free_unit()
                open(iunit, file=filename, status='old')
                ! If first line cannot be parsed as matrix element we skip it
                ! (it would contain the number of orbitals then).
                read(iunit, '(a)', iostat=ierr) char1024
                if (ierr/=0) call stop_all(t_r, "Error reading TCDUMP file")
                ! Only assume first two indices are on same line.
                read(char1024, *, iostat=ierr) matel, indices(1:2)
                if (ierr==0) rewind(iunit)
                counter = 0
                do
                    read(iunit,*, iostat = ierr) matel, indices
                    ! end of file reached?
                    if (ierr < 0) then
                        exit
                    else if (ierr > 0) then
                        ! error while reading?
                        call stop_all(t_r, "Error reading TCDUMP file")
                    else
                        ! permutational factor of 3 (not accounted for in integral calculation)
                        matel = 3.0_dp * matel
                        call freeze_lmat(matel,indices)
                        ! else assign the matrix element

                        if(abs(matel) > LMatEps) then
                            counter = counter + 1
                            index = this%indexFunc(&
                                indices(1),indices(2),indices(3),indices(4),indices(5), indices(6))
                            if(index > this%lMat_size()) then
                                counter = index
                                write(stdout, *) "Warning, exceeding size"
                            endif
                            this%lMat_vals%ptr(index) = matel
                        endif
                    endif

                end do

                write(stdout, *) "Sparsity of LMat", real(counter) / real(this%lMat_size())
                write(stdout, *) "Nonzero elements in LMat", counter
                write(stdout, *) "Allocated size of LMat", this%lMat_size()
                close(iunit)
            end if
            call MPIBcast(counter)
        end if

        ! Synchronize the shared resource
        call this%lMat_vals%sync()
    end subroutine read_dense