Code: GAUSS

!  The below subroutine GAUSS_ELIM solves the input linear equations
!  via gaussian elimination with partial pivoting.

       program gauss

       interface
         subroutine gauss_elim(c,n)
           integer n
           real, dimension(1:n, 1:n+1) :: c
         end subroutine
       end interface

       parameter(n =3)
       integer i
       real, dimension(n*(n+1))::vector
       real, dimension(n, n+1)::augmented

       vector = [ 1.0, 1.0, 1.0, 2.0, -1.0, 2.0, -3.0, 32.0, &
                 3.0, 0.0, -4.0, 17.0 ]
       augmented = reshape( vector, [3,4], order= [2, 1])

       print *, ' Original augmented matrix is : '
       do i = 1, n
         print *, augmented(i, 1:n+1)
       end do

       call gauss_elim(augmented,n)

       print *, '   '

       print *, ' Final augmented matrix (solution in right hand column: '
       do i = 1, n
         print *, augmented(i, 1:n+1)
       end do


       end

!**********************************************************************
   subroutine gauss_elim(mat,n)
   integer n
   real, dimension (1:n, 1:n+1)  :: mat
   ! local var
   integer i,j,pivot_index
   integer, dimension (1) :: pivot_index_array
   real, dimension (1:n+1) :: temprow

!  main loop for gauss elimination

   do i = 1,n
      ! find pivot index
      pivot_index_array = maxloc(abs(mat(i:n,i)))
      pivot_index = pivot_index_array(1) + i - 1
      ! if pivot row not already in place, make it so
      if( pivot_index .ne. i) then
         temprow = mat(i,1:n+1)
         mat(i,1:n+1) = mat(pivot_index,1:n+1)
         mat(pivot_index,1:n+1) = temprow
      endif
      ! check for singular matrix
      if(abs(mat(i,i)) .lt. 1.0e-6) then
         write(*,*) 'error: matrix is singular to machine precision'
         stop
      endif
      mat(i,i:n+1) = mat(i,i:n+1)/mat(i,i)
      do j = i+1,n
         mat(j,i:n+1) = mat(j,i:n+1) - mat(j,i) * mat(i,i:n+1)
      enddo
   enddo

!  back substitution

   do i = n-1,1,-1
      mat(i,n+1) = mat(i,n+1) - sum(mat(i,i+1:n) * mat(i+1:n,n+1))
      mat(i,i+1:n) = 0.0
   enddo

end