! LU decomposition of tridiagonal matrix A
! d(i)=A(i,i)  i=1..n
! L(i)=A(i,i-1) i=2..n
! s(i)=A(i,i+1) i=1..(n-1)
! if ityp=1, returns the LU decomposition of A
! in d(1:n), s(1:n-1) and L(2:n)
! if ityp=2 returns the solution of Ax=f in x
! (a previous call with ityp=2 has to be performed)
subroutine factolu(n,d,L,s,x,f,ityp)
  implicit none
  real*8,dimension(n) :: d,x,f
  real*8,dimension(:),allocatable :: y
  real*8,dimension(1:n-1) :: s
  real*8,dimension(2:n) :: L
  integer :: i,n,ityp

  if (ityp.eq.1) then ! factorization
     do i=2,n
        L(i)=L(i)/d(i-1)
        d(i)=d(i)-L(i)*s(i-1)
     end do
  end if

  if (ityp.eq.2) then
     allocate(y(n))
     ! top-bottom algorithm
     y(1)=f(1)
     do i=2,n
        y(i)=f(i)-L(i)*y(i-1)
     end do
     ! bottom-top algorithm
     x(n)=y(n)/d(n)
     do i=n-1,1,-1
        x(i)=(y(i)-s(i)*x(i+1))/d(i)
     end do
     deallocate(y)
  end if

end subroutine factolu

! simple test of the factolu subroutine
program test
  implicit none
  integer,parameter  :: n=40
  real*8,dimension(1:n) :: d,f,x
  real*8,dimension(2:n) :: L
  real*8,dimension(1:n-1) :: s
  integer    :: i
  real*8    :: dx,xexact,xi

  dx=1.d0/(n+1)

  do i=1,n
     d(i)=2.d0
  end do
  do i=2,n
     L(i)=-1.d0
  end do
  do i=1,n-1
     s(i)=-1.d0
  end do

  do i=1,n
     xi=i*dx
     f(i)=1.d0*dx**2
  end do

  call factolu(n,d,l,s,x(1),f,1)
  call factolu(n,d,l,s,x(1),f,2)

  open(1,file='resu')
  do i=1,n
     write(1,*) i*dx,x(i)
  end do
  write(1,*)
  do i=1,n
     write(1,*) i*dx,xexact(i*dx)
  end do

  close(1)

end program test

function xexact(x)
  implicit none
  integer :: p=3
  real*8 :: x,xexact,pi
  pi=4.d0*atan(1.d0)
  xexact=log(1.d0+x)*sin(2*p*pi*x)
  xexact=x*(1.d0-x)/2
end function xexact
