program tridiagonale
  ! risolve un sistema tridiagonale Ax=f
  ! mediante fattorizzazione LU
  ! i vettori a,b,c contengono gli elementi duagonali,sotto-diagonali e
  ! sopra-diagonali
  implicit none
  real(kind(0.d0)),dimension(:),allocatable :: a,b,c,x,f
  real(kind(0.d0))::m
  integer(4) :: n,i
  
  write(*,*)"n="
  read(*,*) n
  allocate(a(n),b(n-1),c(n-1),x(n),f(n))
  
  a=2
  b=-1
  c=-1
  f=0
  f(1)=1
  f(n)=1
  
  ! fattorizzazione LU
  do i=1,n-1
     m=b(i)/a(i)
     a(i+1)=a(i+1)-m*c(i)
     f(i+1)=f(i+1)-m*f(i)
  end do
  
  !sostituzione
  x(n)=f(n)/a(n)
  do i=n-1,1,-1
     x(i)=(f(i)-c(i)*x(i+1))/a(i)
  end do
  
  write(*,*)"errore=",maxval(abs(x-1))

end program tridiagonale
