program test
  implicit none
  integer,parameter :: n=80
  ! matrice
  real*8 :: L(2:n),u(1:n-1),d(1:n)
  real*8 :: x(n),f(n),erreur,t
  real*8 :: longueur=1,pression=1,eps=1,dx,pi=4*atan(1.d0)
  integer :: i

  dx=longueur/(n+1)

  do i=1,n
     t=i*dx
     f(i)=dx**2*(-0.8D1 * cos(0.2D1 * 0.3141592654D1 * t) ** 2 * 0.3141592654D1 ** 2 + &
          0.8D1 * sin(0.2D1 * 0.3141592654D1 * t) ** 2 * 0.3141592654D1 ** 2)
  end do

  d=2
  L=-1
  u=-1

  call factolu(n,L,d,u,f,x)


  open(101,file='num')
  write(101,*) 0,0
  do i=1,n
     write(101,*) i*dx,x(i)
  end do
  write(101,*) longueur,0
  close(101)

  open(102,file='exa')
  do i=0,n+1
     write(102,*) i*dx,sin(2*pi*i*dx)**2
  end do
  close(102)

  erreur=0
  do i=1,n
     erreur=erreur + &
          (x(i)-sin(2*pi*i*dx)**2)**2
  end do
  write(*,*) n,erreur
  


end program test


subroutine factolu(n,L,d,u,f,x)
  implicit none
  real*8 :: L(2:n),d(1:n),u(1:n-1),f(1:n),x(1:n)
  integer :: i,n
  logical,save :: facto = .true.


  if (facto) then
     facto = .false.
     write(*,*) 'factorisation'
     ! factorisation
     do i=2,n
        L(i)=L(i)/d(i-1)
        d(i)=d(i)-L(i)*u(i-1)
     end do
  end if

  ! descente
  x(1)=f(1)
  do i=2,n
     x(i)=f(i)-L(i)*x(i-1)
  end do

  ! remontée
  x(n)=x(n)/d(n)
  do i=n-1,1,-1
     x(i)=(x(i)-u(i)*x(i+1))/d(i)
  end do

end subroutine factolu




!!$program test
!!$  implicit none
!!$  integer,parameter :: n=10
!!$  ! matrice et sa copie
!!$  real*8 :: L(2:n),u(1:n-1),d(1:n)
!!$
!!$  real*8 :: x(n),f(n),xp(n)
!!$
!!$  integer :: i
!!$
!!$
!!$  d(1)=6+rand()
!!$  u(1)=-1+rand()
!!$  xp(1)=5+rand()
!!$  do i=2,n-1
!!$     d(i)=6+rand()
!!$     L(i)=-2+rand()
!!$     u(i)=-1+rand()
!!$     xp(i)=5+rand()
!!$  end do
!!$  xp(n)=5+rand()
!!$  d(n)=6+rand()
!!$  L(n)=-2+rand()
!!$
!!$  f(1)=d(1)*xp(1)+u(1)*xp(2)
!!$  do i=2,n-1
!!$     f(i)=L(i)*xp(i-1)+d(i)*xp(i)+u(i)*xp(i+1)
!!$  end do
!!$  f(n)=L(n)*xp(n-1)+d(n)*xp(n)
!!$
!!$
!!$
!!$  write(*,*) 'L=',L
!!$  write(*,*) 'd=',d
!!$  write(*,*) 'u=',u
!!$  write(*,*) 'f=',f
!!$
!!$
!!$  call factolu(n,L,d,u,f,x)
!!$
!!$
!!$  write(*,*) 'erreur=',xp-x
!!$
!!$
!!$  f=5
!!$
!!$  call factolu(n,L,d,u,f,x)
!!$
!!$
!!$end program test
