Commit a23c5b28 authored by Rafal Gandecki's avatar Rafal Gandecki
Browse files

Added OpenMP + Fortran 90 code sample for LU decomposition (Doolittle algorithm)

parent 7e79c47f
cmake_minimum_required(VERSION 2.8.7 FATAL_ERROR)
include(${CMAKE_CURRENT_SOURCE_DIR}/../../../cmake/common.cmake)
# ==================================================================================================
if ("${DWARF_PREFIX}" STREQUAL "")
set(DWARF_PREFIX dense_linear_algebra)
endif()
find_package(Common)
find_package(OpenMP)
enable_language (Fortran)
set(NAME ${DWARF_PREFIX}_lud_openmp_fortran)
if (OPENMP_FOUND)
set (CMAKE_Fortran_FLAGS "${CMAKE_CXX_FLAGS} -fopenmp")
add_executable(${NAME} src/lud.f90)
install(TARGETS ${NAME} DESTINATION bin)
else()
message("## Skipping '${NAME}_omp': no OpenMP support found")
dummy_install(${NAME} "OpenMP")
endif()
unset(NAME)
=======
README
=======
# 1. Code sample name
lud_openmp
# 2. Description of the code sample package
This example demonstrates the use of Fortran 90 with OpenMP for LU decomposition (Doolittle algorithm).
# 3. Release date
25 August 2016
# 4. Version history
1.0
# 6. Copyright / License of the code sample
Apache Version 2.0
# 5. Contributor (s) / Maintainer(s)
Rafal Gandecki <rafal.gandecki@pwr.edu.pl>
# 7. Language(s)
Fortran 90
# 8. Parallelisation Implementation(s)
OpenMP
# 9. Level of the code sample complexity
basic
# 10. Instructions on how to compile the code
Uses the CodeVault CMake infrastructure, see main README.md
# 11. Instructions on how to run the code
Just run compiled executable
# 12. Sample input(s)
andomly generated on sample code execution
# 13. Sample output(s)
execution time of algorthims with and without OpenMP
! =================================================================================================
! This file is part of the CodeVault project. The project is licensed under Apache Version 2.0.
! CodeVault is part of the EU-project PRACE-4IP (WP7.3.C).
!
! Author(s):
! Rafal Gandecki <rafal.gandecki@pwr.edu.pl>
!
! This example demonstrates the use of Fortran 90 with OpenMP for LU decomposition (Doolittle algorithm).
!
! See [http://www.openmp.org/] for the full OpenMP documentation.
!
! =================================================================================================
program lud_openmp
integer, parameter :: n = 3000
real, allocatable :: A(:)
real, allocatable :: L(:)
real, allocatable :: U(:)
double precision :: t_config
integer :: t1, t2, clock_rate, clock_max
allocate(A(n*n))
allocate(L(n*n))
allocate(U(n*n))
call fill_random(A, n, n)
call system_clock (t1, clock_rate, clock_max )
call lud_algorithm(A, L, U, n)
call system_clock (t2, clock_rate, clock_max )
t_config = real ( t2 - t1 ) / real ( clock_rate )
print '("Time without OpenMp: ",f6.3," seconds.")', t_config
call system_clock (t1, clock_rate, clock_max )
call lud_openmp_algorithm(A, L, U, n)
call system_clock (t2, clock_rate, clock_max )
t_config = real ( t2 - t1 ) / real ( clock_rate )
print '("Time with OpenMp: ",f6.3," seconds.")', t_config
deallocate(A)
deallocate(L)
deallocate(U)
end program lud_openmp
subroutine fill_random(a, n, m)
implicit none
integer, intent(in) :: n, m
real, dimension(n*m), intent(inout) :: a
integer :: i, j
do i=1, n
do j=1, m
a((i-1)*m+j) = rand()
end do
end do
end subroutine fill_random
subroutine lud_algorithm(A, L, U, n)
integer, intent(in) :: n
real, dimension(n*n), intent(in) :: A
real, dimension(n*n), intent(inout) :: L
real, dimension(n*n), intent(inout) :: U
integer :: i, j, k
do i=1, n
do j=1, n
if (j>i) then
U((j-1)*n+i) = 0
end if
U((i-1)*n+j) = A((i-1)*n+j)
do k=1, i-1
U((i-1)*n+j) = U((i-1)*n+j) - (U((k-1)*n+j) * L((i-1)*n+k))
end do
end do
do j=1, n
if(i>j) then
L((j-1)*n+i) = 0
else if (j==i) then
L((j-1)*n+i) = 1
else
L((j-1)*n+i) = A((j-1)*n+i) / U((i-1)*n+i)
do k=1, i-1
L((j-1)*n+i) = L((j-1)*n+i) - ((U((k-1)*n+i) * L((j-1)*n+k)) / U((i-1)*n+i))
end do
end if
end do
end do
end subroutine
subroutine lud_openmp_algorithm(A, L, U, n)
integer, intent(in) :: n
real, dimension(n*n), intent(in) :: A
real, dimension(n*n), intent(inout) :: L
real, dimension(n*n), intent(inout) :: U
integer :: i, j, k
!$OMP PARALLEL DO DEFAULT (SHARED) PRIVATE(i,j,k)
do i=1, n
do j=1, n
if (j>i) then
U((j-1)*n+i) = 0
end if
U((i-1)*n+j) = A((i-1)*n+j)
do k=1, i-1
U((i-1)*n+j) = U((i-1)*n+j) - (U((k-1)*n+j) * L((i-1)*n+k))
end do
end do
do j=1, n
if(i>j) then
L((j-1)*n+i) = 0
else if (j==i) then
L((j-1)*n+i) = 1
else
L((j-1)*n+i) = A((j-1)*n+i) / U((i-1)*n+i)
do k=1, i-1
L((j-1)*n+i) = L((j-1)*n+i) - ((U((k-1)*n+i) * L((j-1)*n+k)) / U((i-1)*n+i))
end do
end if
end do
end do
!$OMP END PARALLEL DO
end subroutine
subroutine print_array(a, n, m)
implicit none
integer, intent(in) :: n, m
real, dimension(n*m), intent(in) :: a
integer :: i, j
do i=1, n
do j=1, m
print '(F8.4,$)', a((i-1)*m+j)
end do
print *, ''
end do
end subroutine print_array
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment