program tridiag 

!
! Program to solve
!
!        
! | 1 4 0 0 | | x1 |    | 8 |
! | 3 4 1 0 | | x2 | =  | 4 |
! | 0 2 3 4 | | x3 |    | 2 |
! | 0 0 1 3 | | x4 |    | 4 |
! 
! to compile on cobra
!
! gfortran tridiag.f90 /usr/lib64/liblapack.so.3
!

parameter (MAXDIM=4)

! Define program variables

real (kind=8), dimension(MAXDIM) :: A, ASUB, ASUP, ASUP2, B
integer, dimension(MAXDIM) ::  IPIV
integer ::  INFO, N

N = MAXDIM

! Fill the diagonal, superdiagonal and subdiagonal matrix elements  
 
A(1) = 1.0
A(2) = 4.0
A(3) = 3.0
A(4) = 3.0
ASUP(1) = 4.0
ASUP(2) = 1.0
ASUP(3) = 4.0
ASUB(1) = 3.0
ASUB(2) = 2.0  
ASUB(3) = 1.0  

B(1) = 8.0
B(2) = 4.0
B(3) = 2.0
B(4) = 4.0

! Call the appropriate LAPACK routines to solve the tridiagonal system

CALL DGTTRF( N, ASUB, A, ASUP, ASUP2, IPIV, INFO )
CALL DGTTRS( 'N', N, 1, ASUB, A, ASUP, ASUP2, IPIV, B, N, INFO)  
      
! Now compute the values for the various coefficient arrays

do  i = 1, N
  print *, B(i)
enddo

end program tridiag