program test_gauss

  implicit none
  integer,parameter :: n=3,p=1
  real*8,dimension(n,n+p) :: A=0

      A(1,1) = 4
      A(1,2) = 8
      A(1,3) = 12
      A(2,1) = 3
      A(2,2) = 8
      A(2,3) = 13
      A(3,1) = 2
      A(3,2) = 9
      A(3,3) = 18

      A(1,4)=4
      A(2,4)=5
      A(3,4)=11


  call gauss(A,n,p)


end program test_gauss


subroutine gauss(A,n,p)
  implicit none
  integer :: n,p
  real*8,dimension(n,n+p) :: A
  real*8 :: pivot,coef
  integer :: i,k,r,j

  do r=1,n-1
     pivot=A(r,r)
     if (pivot.eq.0) then
        write(*,*) 'Pivot nul'
        stop
     end if
     do i=r+1,n
        coef=A(i,r)/pivot
        do k=1,n+p
           A(i,k)=A(i,k)-coef*A(r,k)
        end do
     end do
  end do

  ! algorithme de remontée
  do k=n+1,n+p
     A(n,k)=A(n,k)/A(n,n)  ! x(n)=b(n)/a(n,n)
  end do

  do i=n-1,1,-1
     do k=n+1,n+p
        do j=i+1,n
           A(i,k)=A(i,k)-A(i,j)*A(j,k)
        end do
        A(i,k)=A(i,k)/A(i,i)
     end do
  end do
 
  write(*,*) 'solution'
  write(*,*) 'A='
  do i=1,n
     write(*,*) A(i,1:n+p)
  end do

end subroutine gauss
