program gaussmp
  use omp_lib
  implicit none
  integer,parameter :: n=1024
  integer :: i,j,k,k0,nbt
  real*8,allocatable :: a(:,:),ac(:,:),xc(:),x(:),b(:)
  real*8 :: akk,aik,temp,tend,tstart



  allocate(a(n,n+1))
  allocate(ac(n,n+1))
  allocate(xc(n))
  allocate(x(n))
  allocate(b(n))

  !$omp parallel
  nbt=omp_get_num_threads()
  !$omp end parallel
  write(*,*) 'Il y a ',nbt,' threads'

! construction d'une matrice arbitraire
! le second membre sera dans a(1:n,n+1)
  do i=1,n
     do j=1,n+1
        !a(i,j)=1.d0/(abs(i-j)+1)
        a(i,j)=1.d0/(i+j-1)
        if (i.eq.j) a(i,j)=a(i,j)+10
     end do
  end do

! solution exacte pour les tests
  do i=1,n
     xc(i)=a(i,n+1)
  end do

! second membre correspondant à la solution
! exacte
  do i=1,n
     a(i,n+1)=0
     do j=1,n
        a(i,n+1)=a(i,n+1)+a(i,j)*xc(j)
     end do
  end do


  write(*,*) 'étape ',0
  !call display(a,n)

! on lance le chronomètre
  tstart=omp_get_wtime()

  ! étape de la méthode du pivot
  do k=1,n-1
     ! rajouter ici la méthode du pivot max
     akk=a(k,k)
     do i=k+1,n
        ! opération Li <- Li - aik/akk * Lk
        aik=a(i,k)
        a(i,1:n+1)=a(i,1:n+1)-aik/akk*a(k,1:n+1)
     end do
     if (mod(k,n/10).eq.0) write(*,*) 'étape ',k,' sur ',n
     !call display(a,n)
  end do

! temps intermédiaire
  tend=omp_get_wtime()
  write(*,*) 'triangulation:',tend-tstart,' s'


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

! arrêt chrono
  tend=omp_get_wtime()


! vérification
  write(*,*)
  write(*,*) maxval(abs(x-xc))
  write(*,*) 'temps=',tend-tstart

end program gaussmp

subroutine display(a,n)
  implicit none
  integer :: n,i
  real*8 :: a(n,n+1)
  do i=1,n
     write(*,*) a(i,1:n+1)
  end do
end subroutine display
