!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! RayleighBernard2D 
! solvers.f90 module
! James Sandham 
! 11 Mar 2014
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!



module solvers
use def


contains
  subroutine tridiag(d,x,a,b,c)
   complex, dimension(n), intent(in) :: d
   complex, dimension(n), intent(inout) :: x
   real, dimension(n), intent(in) :: a,b,c

   integer :: i
   complex, dimension(n) :: dp
   real, dimension(n) :: cp
   real :: m,a0,c0

   !initialize c-prime and d-prime
   cp(1) = c(1)/b(1)
   dp(1) = d(1)/b(1)
   
   !for constant lower and upper off-diagonals
   a0 = a(2)
   c0 = c(2)

   !solve for vectors c-prime and d-prime 
   do i=2,n-1
      m = b(i)-cp(i-1)*a0
      cp(i) = c0/m
      dp(i) = (d(i)-dp(i-1)*a0)/m
   enddo
   m = b(n)-cp(n-1)*a(n)
   cp(n) = c(n)/m
   dp(n) = (d(n)-dp(n-1)*a(n))/m

   !initialize x
   x(n) = dp(n)

   !solve for x from the vectors c-prime and d-prime
   do i=n-1,1,-1
      x(i) = dp(i)-cp(i)*x(i+1)
   enddo
   end subroutine tridiag



   !subroutine tridiag(d,x,a,b,c)
   !complex, dimension(n), intent(in) :: d
   !complex, dimension(n), intent(inout) :: x
   !real, dimension(n), intent(in) :: a,b,c

   !integer :: i
   !complex, dimension(n) :: dp
   !real, dimension(n) :: cp
   !real :: m

   !initialize c-prime and d-prime
   !cp(1) = c(1)/b(1)
   !dp(1) = d(1)/b(1)

   !solve for vectors c-prime and d-prime 
   !do i=2,n
   !   m = b(i)-cp(i-1)*a(i)
   !   cp(i) = c(i)/m
   !   dp(i) = (d(i)-dp(i-1)*a(i))/m
   !enddo

   !initialize x
   !x(n) = dp(n)

   !solve for x from the vectors c-prime and d-prime
   !do i=n-1,1,-1
   !   x(i) = dp(i)-cp(i)*x(i+1)
   !enddo
   !end subroutine tridiag

end module solvers
