!=======================================================================
!
!      Course       : CPS615
!
!      Program Name : Linear.f
!
!      Author       : Kim, Dongmin (096-08-3342)
!
!      This is a program to solve the linear equation A x = b
!
!      where A is N x N matrix and x and b are N x 1 vectors.
!
!      A and b have the initial value, and the goal of this program
!
!      is to find out the vector x using these value.
!
!      In this program, two kind of methods are used :
!
!      one is guassian elimination and the other is Crout algorithm.
!
!      In this program, 1 dimension is used as a column and
!
!      2 dimension is used as a row e.g.
!
!      A(1,1)    A(1,2)   . . .  A(1,Max)
!      A(2,1)    A(2,2)   . . .  A(2,Max)
!        .         .               .
!        .         .               .
!        .         .               .
!
!      A(Max,1)  A(Max,2) . . .  A(Max,Max)
!
!=======================================================================

program Solve_Linear_Equation

      implicit none
      integer, parameter :: Max = 3             ! The maximun size of the problem
      real, dimension( 1 : Max, 1 : Max ) :: A1 ! A matrix
      real, dimension( 1 : Max ) :: x, b1       ! x and b vector
      real, dimension( 1 : Max, 1 : Max ) :: A2 ! A matrix
      real, dimension( 1 : Max ) ::    b2       ! x and b vector
      logical :: error                          ! error check flag
      integer :: i, j
      integer :: operation 

!HPF$ DISTRIBUTE A1(CYCLIC, *)
!HPF$ DISTRIBUTE A1(CYCLIC, *)

      ! Part 1 :
      !   x1 +  x2 +  x3 =  2
      !  -x1 + 2x2 - 3x3 = 32
      !  3x1       - 4x3 = 17

      ! initialize the matrix A1 and the vector b1
      data (( A1( i, j ), j = 1, Max ), i = 1, Max ) &
      / 1.0 , 1.0, 1.0, -1.0, 2.0, -3.0, 3.0, 0.0, -4.0 /
      data ( b1( i ), i = 1, Max ) &
      / 2.0, 32.0, 17.0 /

      print *
      print *
      print *, "A1 matrix"
      do i = 1, Max
         print *, A1( i, :)
      end do

      print *
      print *, "b1 vector"
      do i = 1, Max
         print *, b1(i)
      end do

      print *
      print *, "Guassian elimination"
      print *, "===================="

      ! Guassian elimination
      call Gauss( A1, x, b1, error, Max, operation )

      if ( error ) then
         print *, "Non-linear Equation - Cannot solve"
      else
!         print *, "Float Point Operation : ", operation
         print *
         print *, "Result"
         print *, "======"
         print *, x
      end if

      print *
      print *, "Crout Algorithm"
      print *, "==============="

      ! Crout Algorithm
      call Crout( A1, x, b1, error, Max, operation )

      if ( error ) then
         print *, "Non-linear Equation - Cannot solve"
      else
!         print *, "Float Point Operation : ", operation
         print *
         print *, "Result"
         print *, "======"
         print *, x
      end if

      ! Part 2 :
      !  3x1 + 5x2 + 2x3 =  8
      !        8x2 + 2x3 = -7
      !  6x1 + 2x2 + 8x3 = 26

      ! initialize the matrix A2 and the vector b2
      data (( A2( i, j ), j = 1, Max ), i = 1, Max ) &
      / 3.0, 5.0, 2.0, 0.0, 8.0, 2.0, 6.0, 2.0, 8.0 /
      data ( b2( i ), i = 1, Max ) &
      / 8.0, -7.0, 26.0 /

      print *
      print *
      print *, "A2 matrix"
      do i = 1, Max
         print *, A2( i, :)
      end do

      print *
      print *, "b2 vector"
      do i = 1, Max
         print *, b2(i)
      end do

      print *
      print *, "Guassian elimination"
      print *, "===================="

      ! Guassian elimination
      call Gauss( A2, x, b2, error, Max, operation )

      if ( error ) then
         print *, "Non-linear Equation - Cannot solve"
      else
!         print *, "Float Point Operation : ", operation
         print *
         print *, "Result"
         print *, "======"
         print *, x
      end if

      print *
      print *, "Crout Algorithm"
      print *, "==============="

      ! Crout Algorithm
      call Crout( A2, x, b2, error, Max, operation )

      if ( error ) then
         print *, "Non-linear Equation - Cannot solve"
      else
!         print *, "Float Point Operation : ", operation
         print *
         print *, "Result"
         print *, "======"
         print *, x
      end if

end program Solve_Linear_Equation


subroutine Gauss( A, x, b, error, Max, operation )

      implicit none
      integer :: Max
      real, dimension( 1 : Max, 1 : Max ) :: A    ! A matrix
      real, dimension( 1 : Max) :: x, b           ! x and b vector
      logical :: error                            ! error check flag
      integer :: operation

      real, dimension(size(b), size(b) + 1) :: M  ! array that handles A and b
      real, dimension(size(b) + 1) :: temp_row    ! temporary row
      integer, dimension(1) :: max_location       ! pivot location
      integer :: n, k

      n = size(b)

      operation = 0

      M( 1 : n, 1 : n ) = A  ! copy A to M
      M( 1 : n, n + 1 ) = b  ! copy b to M

      error = .FALSE.

      t_loop: do k = 1, n - 1

         ! find the pivot location in k-th column.
         ! the actual row # is k - 1 + max_location,
         ! since maxloc returns the number that is counting from k.
         max_location = maxloc( abs( M( k : n, k ) ) )

         ! swap the pivot row with k-th row
         temp_row( k : n + 1 ) = M( k, k : n + 1 )
         M( k, k : n + 1 ) = M( k - 1 + max_location(1), k : n + 1 )
         M( k - 1 + max_location(1), k : n + 1 ) = temp_row( k : n + 1 )

         ! check that the chosen pivot is 0
         if ( M( k, k) .EQ. 0 ) then
            error = .TRUE.
            exit t_loop
         else
            ! calculating
            M( k, k + 1 : n + 1 ) = M( k, k + 1 : n + 1 ) / M( k, k)
            operation = operation + 1
            M( k, k ) = 1.
            M( k + 1 : n , k + 1 : n + 1 ) = M( k + 1: n , k + 1 : n + 1 ) - &
            spread( M( k, k + 1 : n + 1 ), 1, n - k) * &
            spread( M( k + 1 : n, k ), 2, n - k + 1 )
            M( k + 1 : n, k ) = 0
         end if
      end do t_loop

      if ( error ) then
         x = 0.0
      else
         ! back substitution
         do k = n, 1, -1
            x( k ) = ( M( k, n + 1 ) - sum( M( k, k + 1 : n ) * & 
            x ( k + 1 : n ) ) ) / M( k, k )
            operation = operation + 1
         end do
      end if

!      return operation

end subroutine Gauss


subroutine Crout( A, x, b, error, Max, operation )

      implicit none
      integer :: Max
      real, dimension( 1 : Max, 1 : Max ) :: A    ! A matrix
      real, dimension( 1 : Max) :: x, b           ! x and b vector
      logical :: error                            ! error check flag
      integer :: operation

      real, dimension(size(b), size(b) + 1) :: M  ! array that has L and U
      real, dimension(size(b) + 1) :: temp_row    ! temporary row
      real, dimension(size(b)) :: y               ! temporary vector
      real, dimension(size(b), size(b)) :: L, U   ! L and U matrix
      integer, dimension(1) :: max_location       ! pivot location
      integer :: n, i, j, k

      n = size(b)

      operation = 0

      M( 1 : n, 1 : n ) = A  ! copy A to M
      M( 1 : n, n + 1 ) = b  ! copy b to M

      error = .FALSE.

      ! LU decomposition using Crout algorithm

      do k = 1, n - 1

         ! find the pivot location in k-th column.
         ! the actual row # is k - 1 + max_location,
         ! since maxloc returns the number that is counting from k.
         max_location = maxloc( abs( M( k : n, k ) ) )

         ! swap the pivot row with k-th row
         temp_row( k : n + 1 ) = M( k, k : n + 1 )
         M( k, k : n + 1 ) = M( k - 1 + max_location(1), k : n + 1 )
         M( k - 1 + max_location(1), k : n + 1 ) = temp_row( k : n + 1 )

         ! check that the chosen pivot is 0
         if ( M( k, k) .EQ. 0 ) then
            error = .TRUE.
            exit
         else
            ! calculating
! Doolittle Algorithm
!           do i = k + 1, n
!               M( i, k ) = M( i, k) / M( k, k)
!               M( i, k + 1 : Max ) = M( i, k + 1 : Max ) - &
!               M( i, k ) * M( k, k + 1 : Max )
!               
!               operation = operation + 2
!            end do

! Crout Algorithm
            forall ( j = k + 1, n )
               M( k, j ) = M( k, j ) / M( k, k)
               M( k + 1 : Max, j ) = M( k + 1 : Max, j ) - &
               M( k, j ) * M( k + 1 : Max, k )
               
               operation = operation + 2
            end forall


         end if
      end do

      ! print out L and U of A
      print *

      do i = 1, Max
         L( i, i + 1 : Max ) = 0.0
         U( i, i ) = 1.0
         U( i + 1 : Max, i ) = 0.0
      end do

      print *, "L part of A"
      print *, "==========="
      do i = 1, Max
         print *, L( i, :)
      end do

      print *
      print *, "U part of A"
      print *, "==========="
      do i = 1, Max
         print *, U( i, :)
      end do

      print *

      if ( error ) then
         x = 0.0
      else
         ! forward substitution
         L = M( :, 1 : n )
         y = M( :, n + 1 )
         do k = 1, n
            y( k ) = ( y( k ) - sum( L( k, 1 : k - 1) * y( 1 : k - 1 ) ) ) &
            / L( k, k )
            operation = operation + 1
         end do

         ! back substitution
         U = M( :, 1 : n )
         do i = 1, n
            U( i, i ) = 1.0
         end do
         x = y
         do k = n - 1, 1, -1
            x( k ) = ( y( k ) - sum( U( k, k + 1 : n ) * x( k + 1 : n ) ) )
            operation = operation + 1
         end do
      end if

!      return operation

end subroutine Crout
