source: trunk/SOURCES/tridiag_mod-0.3.f90 @ 23

!> \file tridiag_mod-0.3.f90
!! Module pour la resolution d'un systeme tri-diagonal
!<

!> \namespace TRIDIAGMOD
!! Module pour la resolution d'un systeme tri-diagonal
!! \author ...
!! \date ...
!<
MODULE TRIDIAGMOD

CONTAINS

  !***************************************************************
  !
  ! tridiag : routine de resolution d'un systeme tri-diagonal
  !  M * U = R   ou M est tridiag.  (ssdiag :A, diag:B, surdiag : C)
  !  U est la solution recherchee
  !
  !******************************************************************

  !> SUBROUTINE: TRIDIAG()
  !! \author ... 
  !! \date ...
  !! @note Cette routine permet de resoudre d'un systeme tri-diagonal
  !! @note  M * U = R   ou M est tridiag.  (ssdiag :A, diag:B, surdiag : C)
  !! @note  U est la solution recherchee
  !! \param A   sous-diagonale l'indice 1 ne sert pas
  !! \param B   diagonale
  !! \param C   sur-diagonale, l'indice n ne sert pas
  !! \param R   vecteur membre de droite
  !! \param U   solution 
  !! \param n  
  !! \param iFAIL     permet de detecter une erreur
  !<


  subroutine TRIDIAG (A,B,C,R,U,n,iFAIL)


    IMPLICIT NONE

    INTEGER, intent(in) :: n
    REAL,dimension(n), intent(in) :: a       !< sous-diagonale l'indice 1 ne sert pas
    REAL,dimension(n), intent(inout) :: b       !< diagonale
    REAL,dimension(n), intent(in) :: c       !< sur-diagonale, l'indice n ne sert pas
    REAL,dimension(n), intent(in) :: r       !< vecteur membre de droite
    REAL,dimension(n), intent(out) :: u       !< solution 
    REAL,dimension(n) :: gam       !< tableau de travail.
    REAL :: BET
    INTEGER,intent(out) :: iFAIL !< permet de detecter une erreur
    INTEGER :: JJ




    iFAIL=0 

    if (abs(B(1)).lt.1.e-20) then 
       B(1)=1.0 
       iFAIL=1 
    endif


    BET=1./B(1) 
    U(1)=R(1)*BET


    !BET=B(1)
    !U(1)=R(1)/BET

    !GAM(:)=EOSHIFT(C,shift=-1)

    do jj=2,n
       GAM(jj)=C(jj-1)*BET
       BET=1./(B(jj)-A(jj)*GAM(jj))

       !   GAM(JJ)=C(JJ-1)/BET
       !   BET=B(JJ)-A(JJ)*GAM(JJ) 


       if (abs(BET)>1.e20) then
          !   if (abs(BET)<1.e-20) then
          BET=1.0 
          iFAIL=1 
       endif

       U(jj) = (R(JJ)-A(JJ)*U(JJ-1))*BET
       !  U(jj) = (R(JJ)-A(JJ)*U(JJ-1))/BET  
    end do

    do JJ=N-1,1,-1
       U(JJ)=U(JJ)-GAM(JJ+1)*U(JJ+1) 
    end do

    if (iFAIL.eq.1) then 
       write(6,*) 'a B element is nul. TRIDIAG cannot work.'
    endif

  end subroutine TRIDIAG

end module tridiagmod