!---------------------------------------------------------------------------
!  program:  serial_lu_solve, Fortran 90 version
!  filename: ser_lu.f
!
!  description:  
!    Example use of LAPACK LU factor and solve routines dgetrf and 
!    dgetrs to solve the system of linear equations Ax = b. The 
!    factorization step is timed and the result reported upon completion.
!
!  associated files: Makefile, lu_solve.in
!
!---------------------------------------------------------------------------
!
!---------------------------------------------------------------------------
!$$$$$$$$$$$$$$$$$$$       Main Program        $$$$$$$$$$$$$$$$$$$$$$$$$$$$$
!---------------------------------------------------------------------------

program serial_lu_solve

implicit none

! the working precision of real variables, define to give double precision
integer, parameter :: wp = kind(1.d0)

! system matrix A and inhomgenous term b
real(kind=wp), allocatable, dimension(:,:) :: A
real(kind=wp), allocatable, dimension(:) :: b

integer :: sys_order

real(kind=wp) :: wall_time

!---------------------------------------------------------------------------
 
! read input and set number of processors
call get_parameters(sys_order)
 
! allocate space for matrices
allocate(A(sys_order,sys_order),b(sys_order))

! initialize and solve system, return timing results
call init_system(A,b)
call solver(A,b,wall_time)

! output timing results
call put_results(sys_order,wall_time)
    
! deallocate space
deallocate(A,b)

!--------------------       subroutines        -----------------------------
contains

!===========================================================================

subroutine get_parameters(sys_order)

! argument declarations
integer, intent(out) :: sys_order

!---------------------------------------------------------------------------

! get input parameters from input file lu_solve.in
open(10,file='lu_solve.in',status='old')

read(10,*) sys_order

close(10)

end subroutine get_parameters

!===========================================================================


subroutine init_system(A,b)

implicit none

! argument declarations
real(kind=wp), dimension(:,:), intent(out) :: A
real(kind=wp), dimension(:), intent(out) :: b

! local variables
integer :: sys_order, i, j

!---------------------------------------------------------------------------

sys_order = size(b)

! initialize system matrix A and ihomogenous term b
do j = 1, sys_order
  do i = 1, sys_order
    A(i,j) = 1.0_wp + sqrt(dble(abs(i-j))/dble(sys_order))
  end do
  b(j) = j
end do

end subroutine init_system

!===========================================================================

subroutine solver(A,b,wall_time)

! argument declarations
real(kind=wp), dimension(:,:), intent(inout) :: A
real(kind=wp), dimension(:), intent(inout) :: b
real(kind=wp), intent(out) :: wall_time

! local variables
integer :: sys_order, nrhs, info
integer, dimension(:), allocatable :: ipvt
real(kind=wp) :: rtc
real(kind=wp), dimension(2) :: wtime

!---------------------------------------------------------------------------

! set number of right hand sides and order of system
nrhs=1
sys_order = size(b)

! allocate space for pivots
allocate(ipvt(sys_order))

!  Call LU factorization routine
wtime(1) = rtc()
call dgetrf(sys_order,sys_order,A,sys_order,ipvt,info)
wtime(2) = rtc()

if (info .ne. 0) print*,'dgetrf info: ',info

!  Call LU solver routine
call dgetrs('N',sys_order,nrhs,A,sys_order,ipvt,b,sys_order,info)

if (info .ne. 0) print*,'dgetrs info: ',info

wall_time = wtime(2) - wtime(1)

deallocate(ipvt)

end subroutine solver

!===========================================================================

subroutine put_header

!---------------------------------------------------------------------------

write(*,*)
write(*,fmt='(a)') 'dgetrf timing'
write(*,fmt='(a)') '============='
write(*,*)

write(*,fmt='(a)') '  system      wall'
write(*,fmt='(a)') '   size       time'
write(*,fmt='(a)') '  ------      ----'

end subroutine put_header

!===========================================================================

subroutine put_results(sys_order,wall_time)

! argument declarations
integer, intent(in) :: sys_order
real(kind=wp), intent(in) :: wall_time

!---------------------------------------------------------------------------

call put_header

write(*,fmt='(2x,i5,5x,f6.2)') sys_order,wall_time

end subroutine put_results

!---------------------------------------------------------------------------

end program serial_lu_solve