The Computer Language
23.03 Benchmarks Game

spectral-norm Fortran gfortran #3 program

source code

! The Computer Language Benchmarks Game
! https://salsa.debian.org/benchmarksgame-team/benchmarksgame/
!
! Original C contributed by Sebastien Loisel
! Conversion to C++ by Jon Harrop
! OpenMP parallelize by The Anh Tran
! Add SSE by The Anh Tran
! Reconversion into C by Dan Farina
! Conversion to Fortran by Brian Taylor

program main
!$ use omp_lib
implicit none

character(len=6) :: argv
integer :: n
real*8, allocatable :: u(:), v(:), tmp(:)
integer :: n2, r_begin, r_end
real*8 uv, vv
integer :: i, tid, tcount, chunk, ite

call get_command_argument(1, argv)
read (argv, *) n

n2 = n / 2

allocate(u(0:n-1), v(0:n-1), tmp(0:n-1))

uv = 0.d0
vv = 0.d0

!$omp parallel default(shared) private(i,tid,tcount,chunk,r_begin,r_end)

!$omp do schedule(static)
do i = 0, n - 1
  u(i) = 1.d0
end do

tid = omp_get_thread_num()
tcount = omp_get_num_threads()
chunk = n / tcount

r_begin = tid * chunk
if (tid < tcount - 1) then
  r_end = r_begin + chunk - 1
else
  r_end = n - 1
end if

do i = 1, 10
  call eval_AtA_times_u(r_begin, r_end, u, v)
  call eval_AtA_times_u(r_begin, r_end, v, u)
end do

!$omp do schedule(static) reduction(+:uv) reduction(+:vv)
do i = 0, n - 1
  uv = uv + u(i) * v(i)
  vv = vv + v(i) * v(i)
end do
!$omp end do nowait

!$omp end parallel

write (*, "(f0.9)") sqrt(uv / vv)

contains


! Return element (i,j) of matrix A
pure function eval_A(i, j)
real*8 :: eval_A
integer, intent(in) :: i, j
real*8 :: di, dj
integer :: d
di = real(i,8)
dj = real(j,8)
eval_A = 1.d0 / (0.5d0 * ((di + dj) * (di + dj + 1.d0)) + di + 1.d0)
end function


subroutine eval_A_times_u(r_begin, r_end, src, dest)
integer, intent(in) :: r_begin, r_end
real*8, intent(in) :: src(0:)
real*8, intent(out) :: dest(0:)
real*8 sum1
integer :: i, j
do i = r_begin, r_end
  sum1 = 0.d0
  do j = 0, n - 1
    sum1 = sum1 + src(j) * eval_A(i, j)
  end do
  dest(i) = sum1
end do
end subroutine


subroutine eval_At_times_u(r_begin, r_end, src, dest)
integer, intent(in) :: r_begin, r_end
real*8, intent(in) :: src(0:)
real*8, intent(out) :: dest(0:)
real*8 sum1
integer :: i, j
do i = r_begin, r_end
  sum1 = 0.d0
  do j = 0, n - 1
    sum1 = sum1 + src(j) * eval_A(j, i)
  end do
  dest(i) = sum1
end do
end subroutine


subroutine eval_AtA_times_u(r_begin, r_end, src, dest)
integer, intent(in) :: r_begin, r_end
real*8, intent(in) :: src(0:)
real*8, intent(out) :: dest(0:)
call eval_A_times_u(r_begin, r_end, src, tmp)
!$omp barrier
call eval_At_times_u(r_begin, r_end, tmp, dest)
!$omp barrier
end subroutine

end program
    

notes, command-line, and program output

NOTES:
64-bit Ubuntu quad core
Ubuntu 12.2.0-3ubuntu1


Wed, 25 Jan 2023 03:32:20 GMT

MAKE:
/usr/bin/gfortran -pipe -O3 -fomit-frame-pointer -march=ivybridge  -fopenmp spectralnorm.gfortran-3.f95 -o spectralnorm.gfortran-3.gfortran_run 
rm spectralnorm.gfortran-3.f95

2.55s to complete and log all make actions

COMMAND LINE:
./spectralnorm.gfortran-3.gfortran_run 5500

PROGRAM OUTPUT:
1.274224153