История изменений
Исправление grem, (текущая версия) :
Вот тебе для коллекции fortran.f03:
!! gfortran -std=f2003 -Ofast -march=native
program matrix_mult
implicit none
integer, parameter :: dp = kind(1.d0)
character(100) :: argv
integer :: n, IOSTATUS
real(dp), dimension(:, :), allocatable :: AA, BB, CC
real(dp) :: start, finish
if(command_argument_count() .eq. 1) then
call get_command_argument(1, argv)
read(argv, *, IOSTAT=IOSTATUS) n
if (IOSTATUS .ne. 0) then
write(*, *) 'Error: Comman-line argument must be integer. Stop.'
stop
endif
else
write(*, *) 'Error: One command-line argument is required. Stop.'
stop
end if
allocate(AA(n, n))
allocate(BB(n, n))
allocate(CC(n, n))
call matrix_fill(AA, n)
call matrix_fill(BB, n)
call matrix_print(AA, n)
call cpu_time(start)
call matrix_prod(AA, BB, CC, n)
call cpu_time(finish)
print '("time: ", f10.3, " seconds")', (finish - start)
write(*, '(f12.6)') CC(n/2 + 1, n/2 + 1)
deallocate(AA)
deallocate(BB)
deallocate(CC)
contains
subroutine matrix_fill(A, n)
implicit none
real(dp), dimension(:, :), intent(out) :: A
integer, intent(in) :: n
real(dp) :: tmp
integer :: i, j
tmp = 1.0_dp / n / n
do j = 0, n-1
do i = 0, n-1
A(i+1, j+1) = tmp * (j - i) * (j + i)
end do
end do
end subroutine matrix_fill
subroutine matrix_print(A, n)
implicit none
real(dp), dimension(:, :), intent(in) :: A
integer, intent(in) :: n
integer :: j
if ( n .le. 10 ) then
do j = 1, n
write(*, '(*(f12.6))') A( : , j )
end do
else
print '("Matrix is too large to print, n = ", (i4))', n
end if
end subroutine matrix_print
subroutine matrix_prod(A, B, C, n)
implicit none
real(dp), dimension(:, :), intent(in) :: A, B
real(dp), dimension(:, :), intent(out) :: C
integer, intent(in) :: n
real(dp), dimension(:, :), allocatable :: T
integer :: i, j
allocate(T(n, n))
T = TRANSPOSE(B)
do j = 1, n
do i = 1, n
C(i, j) = sum( A(:, i) * T(:, j) )
end do
end do
deallocate(T)
end subroutine matrix_prod
end program matrix_mult
Для gfortran из mingw64 при n = 3000 работает либо так же, либо на 0.5 секунды медленнее варианта на C, при временах ~13 секунд. Если использовать для перемножения встроенную функцию matmult, то ответ выдаёт через ~3 секунды.
Если хочется распаралеллить, то внутри реализации процедуры subroutine matrix_prod(A, B, C, n) вместо do j = 1,n пишешь do concurrent (j = 1:n), а при сборке вместо флага -std=2003 указываешь флаг -std=2008 -ftree-parallelize-loops=4, где 4 - колличество потоков. Будет выполняться секунд за 9.
Исходная версия grem, :
Вот тебе для коллекции fortran.f03:
!! gfortran -std=f2003 -Ofast -march=native
program matrix_mult
implicit none
integer, parameter :: dp = kind(1.d0)
character(100) :: argv
integer :: n, IOSTATUS
real(dp), dimension(:, :), allocatable :: AA, BB, CC
real(dp) :: start, finish
if(command_argument_count() .eq. 1) then
call get_command_argument(1, argv)
read(argv, *, IOSTAT=IOSTATUS) n
if (IOSTATUS .ne. 0) then
write(*, *) 'Error: Comman-line argument must be integer. Stop.'
stop
endif
else
write(*, *) 'Error: One command-line argument is required. Stop.'
stop
end if
allocate(AA(n, n))
allocate(BB(n, n))
allocate(CC(n, n))
call matrix_fill(AA, n)
call matrix_fill(BB, n)
call matrix_print(AA, n)
call cpu_time(start)
call matrix_prod(AA, BB, CC, n)
call cpu_time(finish)
print '("time: ", f10.3, " seconds")', (finish - start)
write(*, '(f12.6)') CC(n/2 + 1, n/2 + 1)
deallocate(AA)
deallocate(BB)
deallocate(CC)
contains
subroutine matrix_fill(A, n)
implicit none
real(dp), dimension(:, :), intent(inout) :: A
integer, intent(in) :: n
real(dp) :: tmp
integer :: i, j
tmp = 1.0_dp / n / n
do j = 0, n-1
do i = 0, n-1
A(i+1, j+1) = tmp * (j - i) * (j + i)
end do
end do
end subroutine matrix_fill
subroutine matrix_print(A, n)
implicit none
real(dp), dimension(:, :), intent(in) :: A
integer, intent(in) :: n
integer :: j
if ( n .le. 10 ) then
do j = 1, n
write(*, '(*(f12.6))') A( : , j )
end do
else
print '("Matrix is too large to print, n = ", (i4))', n
end if
end subroutine matrix_print
subroutine matrix_prod(A, B, C, n)
implicit none
real(dp), dimension(:, :), intent(in) :: A, B
real(dp), dimension(:, :), intent(out) :: C
integer, intent(in) :: n
real(dp), dimension(:, :), allocatable :: T
integer :: i, j
allocate(T(n, n))
T = TRANSPOSE(B)
do j = 1, n
do i = 1, n
C(i, j) = sum( A(:, i) * T(:, j) )
end do
end do
deallocate(T)
end subroutine matrix_prod
end program matrix_mult
Для gfortran из mingw64 при n = 3000 работает либо так же, либо на 0.5 секунды медленнее варианта на C, при временах ~13 секунд. Если использовать для перемножения встроенную функцию matmult, то ответ выдаёт через ~3 секунды.
Если хочется распаралеллить, то внутри реализации процедуры subroutine matrix_prod(A, B, C, n) вместо do j = 1,n пишешь do concurrent (j = 1:n), а при сборке вместо флага -std=2003 указываешь флаг -std=2008 -ftree-parallelize-loops=4, где 4 - колличество потоков. Будет выполняться секунд за 9.