module solvers_1D

! objetivo: resolver sistema linear de equações algébricas
!           originado de problemas unidimensionais

use variaveis

implicit none

contains

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

  ! Método direto eliminação de Gauss

  subroutine EG (N,ap,b,c,d,T)

    implicit none

    integer :: i      ! número do nó
	integer :: ii, j
	real*8  :: mji, S

    integer,intent(in) :: N ! número de nós

    real*8,dimension(:,:),allocatable :: A ! matriz de coeficientes

    real*8,intent(in), dimension(N) :: ap ! coeficiente aP
    real*8,intent(in), dimension(N) :: b  ! coeficiente aW
    real*8,intent(in), dimension(N) :: c  ! coeficiente aE
    real*8,intent(in), dimension(N) :: d  ! termo fonte bP

    real*8,intent(out),dimension(N) :: T  ! incógnita

    allocate ( A(N,N+1) )
	
    A = 0.0d0

    ! gera a matriz de coeficientes A com o termo independente

    do i = 1,N

       if (i > 1) A(i,i-1) = -b(i)

       A(i,i) = ap(i)

       if (i < N) A(i,i+1) = -c(i)

       A(i,N+1) = d(i)

    end do

	! Escalonamento

	do i = 1,N-1
	
		do ii = i+1,N 
		
			mji= A(ii,i) / A(i,i)

			do j = i,N+1

				A(ii,j) = A(ii,j) - mji*A(i,j)
			
			end do

		end do

	end do 

	!Substituicao retroativa

	T(N) = A(N,N+1) / A(N,N)

	do i = N-1,1,-1

		S = 0
		
		do j = i+1,N

			S = S + A(i,j)*T(j)	

		end do

		T(i) = (A(i,N+1) - S) / A(i,i)

	end do 

    deallocate ( A )

end subroutine EG

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

  ! método iterativo de Gauss-Seidel

  subroutine GS (N,ite,tol,a,b,c,d,T)

    implicit none

    integer :: i   ! número do nó
    integer :: it  ! número da iteração
    integer :: ite ! número de iterações

    integer,intent(in) :: N   ! número de nós

    real*8, intent(in) :: Tol ! tolerância sobre R

    real*8,intent(in), dimension(N) :: a ! coeficiente aP
    real*8,intent(in), dimension(N) :: b ! coeficiente aW
    real*8,intent(in), dimension(N) :: c ! coeficiente aE
    real*8,intent(in), dimension(N) :: d ! termo fonte bP

    real*8,intent(out),dimension(N) :: T ! incógnita

    real*8 :: R ! norma l1 média dos resíduos

    T = 0.0d0

    do it = 1, ite

       T(1) = ( c(1)*T(2) + d(1) ) / a(1)

       do i = 2, N-1

          T(i) = ( b(i)*T(i-1) + c(i)*T(i+1) + d(i) ) / a(i)

       end do

       T(N) = ( b(N)*T(N-1) + d(N) ) / a(N)

       call norma (N,a,b,c,d,T,R)

       if ( R <= Tol ) then
          write(10,1) it, Tol, R
          1 format(//, 'O processo iterativo convergiu em', I8, ' iterações', &
                    /, 'para a tolerância de', 1pe10.2, &
                    /, 'A norma l1 média dos resíduos é', 1pe10.2, /)
          exit
       end if

    end do

    if ( R > Tol ) then
       write(10,2) ite, Tol, R
       2 format(//, 'O processo iterativo NÃO convergiu em', I8, ' iterações', &
                 /, 'para a tolerância de', 1pe10.2, &
                 /, 'A norma l1 média dos resíduos é', 1pe10.2, /)
    end if

  end subroutine GS

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

  ! calcula a norma l1 média do resíduo das equações

  subroutine norma (N,a,b,c,d,T,R)

    implicit none

    integer :: i ! número do nó

    integer,intent(in) :: N ! número de nós

    real*8,intent(in), dimension(N) :: a ! coeficiente aP
    real*8,intent(in), dimension(N) :: b ! coeficiente aW
    real*8,intent(in), dimension(N) :: c ! coeficiente aE
    real*8,intent(in), dimension(N) :: d ! termo fonte bP

    real*8,intent(out),dimension(N) :: T ! incógnita

    real*8,intent(inout) :: R ! norma l1 média dos resíduos

    R = 0.0d0

    R = R + dabs ( c(1)*T(2) + d(1) - T(1)*a(1) )

    do i = 2, N-1

       R = R + dabs ( b(i)*T(i-1) + c(i)*T(i+1) + d(i) - T(i)*a(i) )

    end do

    R = R + dabs ( b(N)*T(N-1) + d(N) - T(N)*a(N) )

    R = R / N

  end subroutine norma

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

  ! método direto Tri-Diagonal Matrix Algorithm (TDMA)

  subroutine TDMA (N,a,b,c,d,T)

    implicit none

    integer :: i   ! número do nó
    real*8  :: div ! variável auxiliar

    integer,intent(in) :: N ! número de nós

    real*8,dimension(:),allocatable :: P ! coeficiente do tdma
    real*8,dimension(:),allocatable :: Q ! coeficiente do tdma

    real*8,intent(in), dimension(N) :: a ! coeficiente aP
    real*8,intent(in), dimension(N) :: b ! coeficiente aW
    real*8,intent(in), dimension(N) :: c ! coeficiente aE
    real*8,intent(in), dimension(N) :: d ! termo fonte bP

    real*8,intent(out),dimension(N) :: T ! incógnita

    allocate(P(N),Q(N))

    P(1) = c(1) / a(1)
    Q(1) = d(1) / a(1)

    do i = 2, N
      div  = a(i) - b(i)*P(i-1)
      P(i) = c(i) / div
      Q(i) = (d(i) + b(i)*Q(i-1))/div
    end do

    T(N) = Q(N)

    do i = N-1, 1, -1
      T(i) = P(i)*T(i+1) + Q(i)
    end do

    deallocate(P,Q)

  end subroutine tdma

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

end module solvers_1D
