c                                                                      |
c - * -- 1 -- * -- 2 -- * -- 3 -- * -- 4 -- * -- 5 -- * -- 6 -- * -- 7 +
c                                                                      |
      function gsnelm(mat, vec, dim, mxdm)
c
      implicit none
      integer :: mxdm, dim, gsnelm
      double precision :: mat(mxdm,mxdm), vec(mxdm)
c
      integer :: pivot, swprws, sclrw, elmrw
      integer :: i, j, k, dmk, ip, rc
      double precision :: eps
c
      eps = 5.0d-16
      if (dim .gt. mxdm) then
        gsnelm = -1
        return
      end if
c
      do dmk = 1, dim - 1
        k = dim - dmk + 1
c
        ip = pivot(mat, k, dim, mxdm)
        if (ip .ne. k) then
          rc = swprws(mat, vec, ip, k, dim, mxdm)
        end if
c
        if (abs(mat(k,k)) .lt. eps) then
          gsnelm = k
          return
        end if
c
        rc = sclrw(mat, vec, k, dim, mxdm)
        do i = 1, k - 1
          rc = elmrw(mat, vec, k, i, dim, mxdm)
        end do
c
      end do
c
      if (abs(mat(1,1)) .lt. eps) then
        gsnelm = 1
        return
      end if
c
      vec(1) = vec(1) / mat(1,1)
      do k = 2, dim
        do j = 1, k - 1
          vec(k) = vec(k) - mat(k, j) * vec(j)
        end do
      end do
c
      gsnelm = 0
c
      end
c                                                                      |
c - * -- 1 -- * -- 2 -- * -- 3 -- * -- 4 -- * -- 5 -- * -- 6 -- * -- 7 +
c                                                                      |
      function pivot(mat, k, dim, mxdm)
c
      implicit none
      integer :: mxdm, dim, k, pivot
      double precision :: mat(mxdm,mxdm)
c
      integer :: i
c
      pivot = k
      do i = 1, k - 1
        if(abs(mat(pivot,k)) .lt. abs(mat(i,k))) then
          pivot = i
        end if
      end do
c
      end
c                                                                      |
c - * -- 1 -- * -- 2 -- * -- 3 -- * -- 4 -- * -- 5 -- * -- 6 -- * -- 7 +
c                                                                      |
      function swprws(mat, vec, ip, k, dim, mxdm)
c
      implicit none
      integer :: mxdm, dim, ip, k, swprws
      double precision :: mat(mxdm,mxdm), vec(mxdm)
c
      integer :: j
      double precision :: swp
c
      do j = 1, k
        swp = mat(ip,j)
        mat(ip,j) = mat(k,j)
        mat(k,j) = swp
      end do
c
      swp = vec(ip)
      vec(ip) = vec(k)
      vec(k) = swp
c
      swprws = 0
      end
c                                                                      |
c - * -- 1 -- * -- 2 -- * -- 3 -- * -- 4 -- * -- 5 -- * -- 6 -- * -- 7 +
c                                                                      |
      function sclrw(mat, vec, ip, dim, mxdm)
c
      implicit none
      integer :: mxdm, dim, ip, sclrw
      double precision :: mat(mxdm,mxdm), vec(mxdm)
      integer :: j
      double precision :: pv
c
      pv = mat(ip,ip)
      do j = 1, ip - 1
        mat(ip,j) = mat(ip,j) / pv
      end do
      mat(ip,ip) = 1.0d+0
      vec(ip) = vec(ip) / pv
c
      sclrw = 0
      end
c                                                                      |
c - * -- 1 -- * -- 2 -- * -- 3 -- * -- 4 -- * -- 5 -- * -- 6 -- * -- 7 +
c                                                                      |
      function elmrw(mat, vec, ip, i, dim, mxdm)
c
      implicit none
      integer :: mxdm, dim, ip, i, elmrw
      double precision :: mat(mxdm,mxdm), vec(mxdm)
c
      integer :: j
      double precision :: pv
c
      pv = mat(i,ip)
      do j = 1, ip - 1
        mat(i,j) = mat(i,j) - mat(ip,j) * pv
      end do
      mat(i,ip) = 0.0d+0
      vec(i) = vec(i) - vec(ip) * pv
c
      elmrw = 0
      end
c                                                                      |
c - * -- 1 -- * -- 2 -- * -- 3 -- * -- 4 -- * -- 5 -- * -- 6 -- * -- 7 +
c                                                                      |
