timing array unformatted writing

timed.f90

module kinds
  save
  integer, parameter :: dp = selected_real_kind(15, 307)
end module kinds


program am_timed

    use kinds, only : dp

    implicit none
    integer, parameter       :: asz = 1e6
    real(dp), dimension(asz) :: rarr1
    real(dp), allocatable    :: rarr2(:)
    integer                  :: ierr,i,j
    real(dp)                 :: t1, t2, onek
    character(256)           :: fnam

    onek = 1.0e3_dp
    allocate( rarr2(asz) )

    do i = 1, asz ! fill the arrays with garbage
        rarr1(i) = real(i, kind=dp)
        rarr2(i) = 0.1_dp * real(i, kind=dp)
    end do

    write(*,'(2x,"Time to write arrays of size:", es8.1e2)') real(asz)
    write(*,'(2x,"Loop Type,  Array Type,  Time [ms]")')

    fnam = 'stackarr_do.dat'
    call cpu_time(t1)
    call do_write(rarr1, asz, fnam)
    call cpu_time(t2)
    write(*,1001) 'do', 'stack', (t2-t1)*onek
    call do_verify(rarr1, asz, fnam)

    fnam = 'heaparr_do.dat'
    call cpu_time(t1)
    call do_write(rarr2, asz, fnam)
    call cpu_time(t2)
    write(*,1001) 'do', 'heap', (t2-t1)*onek
    call do_verify(rarr2, asz, fnam)

    fnam = 'stackarr_seq.dat'
    call cpu_time(t1)
    call seq_write(rarr1, asz, fnam)
    call cpu_time(t2)
    write(*,1001) 'seq', 'stack', (t2-t1)*onek
    call seq_verify(rarr1, asz, fnam)

    fnam = 'heaparr_seq.dat'
    call cpu_time(t1)
    call seq_write(rarr2, asz, fnam)
    call cpu_time(t2)
    write(*,1001) 'seq', 'heap', (t2-t1)*onek
    call seq_verify(rarr2, asz, fnam)

    fnam = 'stackarr_duff.dat'
    call cpu_time(t1)
    call bloc_write(rarr1, asz, fnam)
    call cpu_time(t2)
    write(*,1001) 'block', 'stack', (t2-t1)*onek
    call do_verify(rarr1, asz, fnam)

    fnam = 'heaparr_duff.dat'
    call cpu_time(t1)
    call bloc_write(rarr2, asz, fnam)
    call cpu_time(t2)
    write(*,1001) 'block', 'heap', (t2-t1)*onek
    call do_verify(rarr2, asz, fnam)

1001 format(2x,a9,2x,a11,2x,f10.4)
end program am_timed


subroutine myopen(fnam, iun, ierr)
    ! simple open and error check
    implicit none
    character(256), intent(in) :: fnam
    integer, intent(in)        :: iun
    integer, intent(out)       :: ierr

    open(unit=iun, file=trim(fnam), form='unformatted', iostat=ierr)
    if ( ierr /= 0 ) then
        write(*,*) 'Unable to open '//trim(fnam)//', exiting'
        call exit(1)
    end if
end subroutine myopen


subroutine do_write(array, asz, fnam)
    ! write with a do-loop
    use kinds, only : dp
    implicit none
    ! Params
    integer, intent(in)  :: asz
    real(dp), intent(in) :: array(asz)
    character(256)       :: fnam
    ! Locals
    integer              :: i, iun, ierr

    iun = 17

    call myopen(fnam, iun, ierr)
    do i = 1, asz
        write(iun) array(i)
    end do
    close(iun)
end subroutine do_write


subroutine seq_write(array, asz, fnam)
    ! write with sequentially
    use kinds, only : dp
    implicit none
    ! Params
    integer, intent(in)  :: asz
    real(dp), intent(in) :: array(asz)
    character(256)       :: fnam
    ! Locals
    integer              :: i, iun, ierr

    iun = 17

    call myopen(fnam, iun, ierr)
    write(iun) array
    close(iun)
end subroutine seq_write


subroutine bloc_write(array, asz, fnam)
    ! write with unrolled do-loop & duff's device
    use kinds, only : dp
    implicit none
    ! Params
    integer, intent(in)  :: asz
    real(dp), intent(in) :: array(asz)
    character(256)       :: fnam
    ! Locals
    integer              :: i, j, iun, ierr
    integer              :: dufflim, dolim
    integer              :: bsz ! block size

    iun = 17
    call myopen(fnam, iun, ierr)

    bsz = 16
    if ( asz .lt. bsz ) then ! smaller than block case
        do i = 1, asz
            write(iun) array(i)
        end do
        close(iun)
        return
    end if

    ! Block expand do loop in bsz units
    do i = 1, asz-bsz+1, bsz
        write(iun) array(i)
        write(iun) array(i+1)
        write(iun) array(i+2)
        write(iun) array(i+3)
        write(iun) array(i+4)
        write(iun) array(i+5)
        write(iun) array(i+6)
        write(iun) array(i+7)
        write(iun) array(i+8)
        write(iun) array(i+9)
        write(iun) array(i+10)
        write(iun) array(i+11)
        write(iun) array(i+12)
        write(iun) array(i+13)
        write(iun) array(i+14)
        write(iun) array(i+15)
    end do

    ! finish the remainder with a duff-device
    dufflim = mod(asz, bsz) + 1 ! +1 because goto 0 maps to goto 1
    dolim = asz - bsz + 2
1399 go to ( 1300, 1301, 1302, 1303, 1304, 1305, 1306, 1307, &
             1308, 1309, 1310, 1311, 1312, 1313, 1314, 1315 ), dufflim
1315 write(iun) array(dolim)
1314 write(iun) array(dolim+1)
1313 write(iun) array(dolim+2)
1312 write(iun) array(dolim+3)
1311 write(iun) array(dolim+4)
1310 write(iun) array(dolim+5)
1309 write(iun) array(dolim+6)
1308 write(iun) array(dolim+7)
1307 write(iun) array(dolim+8)
1306 write(iun) array(dolim+9)
1305 write(iun) array(dolim+10)
1304 write(iun) array(dolim+11)
1303 write(iun) array(dolim+12)
1302 write(iun) array(dolim+13)
1301 write(iun) array(dolim+14)
1300 continue

    close(iun)
end subroutine bloc_write


subroutine do_verify(array, asz, fnam)
    ! verify array is EXACTLY the same as that stored in fnam
    use kinds, only : dp
    implicit none
    ! Params
    integer, intent(in)   :: asz
    real(dp), intent(in)  :: array(asz)
    character(256)        :: fnam
    ! Locals
    real(dp), allocatable :: array2(:)
    integer               :: i, iun
     
    allocate ( array2(asz) )

    iun = 13
    call myopen(fnam, iun, i)

    do i = 1, asz
        read(iun) array2(i)
        if ( array(i) .ne. array2(i) ) then
            write(*,*) 'Mismatch found: ', array2(i), ', ', array(i)
            call exit(2)
        end if
    end do

    close(iun)
end subroutine do_verify


subroutine seq_verify(array, asz, fnam)
    ! verify array is EXACTLY the same as that stored in fnam (sequential)
    use kinds, only : dp
    implicit none
    ! Params
    integer, intent(in)   :: asz
    real(dp), intent(in)  :: array(asz)
    character(256)        :: fnam
    ! Locals
    real(dp), allocatable :: array2(:)
    integer               :: i, iun
     
    allocate ( array2(asz) )

    iun = 13
    call myopen(fnam, iun, i)

    read(iun) array2
    do i = 1, asz
        if ( array(i) .ne. array2(i) ) then
            write(*,*) 'Mismatch found: ', array2(i), ', ', array(i)
            call exit(2)
        end if
    end do

    close(iun)
end subroutine seq_verify