! Copyright (C) HPC WORLD (Prometech Software, Inc. and GDEP Solutions, Inc.) All rights reserved.
!
! Permission is hereby granted, free of charge, to any person obtaining a
! copy of this software and associated documentation files (the "Software"),
! to deal in the Software without restriction, including without limitation
! the rights to use, copy, modify, merge, publish, distribute, sublicense,
! and/or sell copies of the Software, and to permit persons to whom the
! Software is furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in
! all copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
! THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
! FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
! DEALINGS IN THE SOFTWARE.

program main
   implicit none
   integer, parameter :: m = 5      ! row number of matrix A
   integer, parameter :: n = m      ! column number of matrix A
   integer, parameter :: nrhs = 2   ! number of vector b
   integer, parameter :: lda = m
   integer, parameter :: ldb = n
   real(8), parameter :: pi = datan(1.0d0)*4.0d0
   integer :: i, j, k, min_mn
   real(8) :: a(lda, n)       ! m x m matrix A
   real(8) :: b(ldb, nrhs)    ! n x nrhs RHS (set of vectors b)
   real(8) :: a1(lda, n)      ! initial A
   real(8) :: b1(ldb, nrhs)   ! initial b
   real(8) :: b2(ldb, nrhs)   ! check
   real(8) :: theta(m)        ! variable theta
   real(8) :: init_L(m, m)    ! initial matrix L
   real(8) :: init_U(m, m)    ! initial matrix U
   real(8) :: temp
   real(8) :: x1, x2
   real(8) :: sum0, relative_error
   real(8) :: phase, detA

   ! LAPACK
   integer, allocatable :: ipiv(:)
   integer :: info_f, info_s

   ! Elapsed time function initialized
   integer ic, icr, icm
   real*8 cpu0, cpu1, elapsed, t_ac
   real*8 second
   external second

   write (*, *) "Dgetrs (LAPACK)"

   ! 時間の計測の精度
   call system_clock(ic, icr, icm)
   t_ac = 1.0d0/real(icr, 8)
   print '(1x,a,e10.5)', 'Time measurement accuracy : ', t_ac

   write (*, *) 'lda,ldb', lda, ldb

   ! elements of matrix A
   ! ここでは、LとUを作り、それらをかけて行列Aを作る。
   ! calculation of theta
   temp = pi/(dble(m)*2.0d0)
   do i = 1, m
      theta(i) = dble(i)*temp   ! 適当な数
   end do

   ! initial matrix L
   ! (ここでは、init_Lを転置で定義しているので注意。後のloopでの演算で利点がある。)
   ! initial matrix U
   ! case: i = j
!   write (*, *) "case: i=j"
   do i = 1, m
      init_U(i, i) = 1.0d0   ! det(A) = det(U) = 1
      init_L(i, i) = 1.0d0
   end do

   ! case: i > j (lower)
!   write (*, *) "i>j"
   do j = 1, m - 1
      do i = j + 1, m
         init_U(i, j) = 0.0d0
         init_L(i, j) = 0.0d0
      end do
   end do

   ! case: i < j (upper)
!   write (*, *) "i<j"
   do i = 1, m - 1
      x1 = dsin(theta(i))   ! 適当な数
      do j = i + 1, m
         x2 = dcos(theta(j))   ! 適当な数
         init_U(i, j) = x1 + x2
         init_L(i, j) = init_U(i, j)
      end do
   end do

   ! A = L * U ()
   write (*, *) "A=LU"
   do i = 1, m
      do j = 1, m
         a(i, j) = 0.0d0   ! initialization
         do k = 1, m
            a(i, j) = a(i, j) + init_L(k, i)*init_U(k, j)
         end do
      end do
   end do

   ! generate vector b
   write (*, *) "generate b"
   do k = 1, nrhs
      x1 = k*0.1d0   ! 適当な数
      do i = 1, m
         b(i, k) = dsin(theta(i)) + x1   ! 適当な数
      end do
   end do

   ! display
   print *, "Matrix A"
   do i = 1, m
      write (*, "(5(f8.2,2x))") a(i, :)
   end do
   print *, "RHS"
   do i = 1, m
      write (*, "(5(f8.2,2x))") b(i, :)
   end do

   ! stock A and b
   a1 = a
   b1 = b

   min_mn = min(m, n)        ! mとnの小さい方をmin_mnに代入
   allocate (ipiv(min_mn))   ! allocateにより、配列ipivの次元をここで定義

   cpu0 = second()   ! 時間の計測を開始

   ! LU factorization
   call Dgetrf(m, n, a, lda, ipiv, info_f)
   ! solve simultaneous equation
   ! 'N': non-transpose of the matrix A
   call Dgetrs('N', n, nrhs, a, lda, ipiv, b, ldb, info_s)

   cpu1 = second() ! 時間の計測を終了

   ! Printout Elapsed time
   elapsed = (cpu1 - cpu0)*t_ac
   print *, "time: ", elapsed, "sec"

   ! check
   if (info_f /= 0) write (*, *) 'Dgetrf failed'
   if (info_s /= 0) write (*, *) 'Dgetrs failed'

   ! output
   print *, "LU factorization"
   do i = 1, lda
      write (*, "(5(f8.2,2x))") a(i, :)
   end do
   print *, "Solution"
   do i = 1, ldb
      write (*, '(5(e10.3,2x))') b(i, :)
   end do

   ! considering pivot (or permutation) matrix
   k = 0 ! initialization
   do i = 1, lda
      if (ipiv(i) .gt. i) k = k + 1
   end do
   phase = (-1.0d0)**k

   ! calculation of det(A)
   detA = phase ! initialization
   do i = 1, lda
      detA = detA*a(i, i)
   end do
   write (*, *) 'det(A):', detA

end program main

function second() result(rtime)   !-- SofTek

   implicit none
   integer :: ic, ir, im
   real(8) :: rtime

   call system_clock(ic, ir, im)
   rtime = real(ic, 8)

   return
end function