c     ..................................................................
c
c        subroutine dgelg
c
c        purpose
c           to solve a general system of simultaneous linear equations.
c
c        usage
c           call dgelg(r,a,m,n,eps,ier)
c
c        description of parameters
c           r      - double precision m by n right hand side matrix
c                    (destroyed). on return r contains the solutions
c                    of the equations.
c           a      - double precision m by m coefficient matrix
c                    (destroyed).
c           m      - the number of equations in the system.
c           n      - the number of right hand side vectors.
c           eps    - single precision input constant which is used as
c                    relative tolerance for test on loss of
c                    significance.
c           ier    - resulting error parameter coded as follows
c                    ier=0  - no error,
c                    ier=-1 - no result because of m less than 1 or
c                             pivot element at any elimination step
c                             equal to 0,
c                    ier=k  - warning due to possible loss of signifi-
c                             cance indicated at elimination step k+1,
c                             where pivot element was less than or
c                             equal to the internal tolerance eps times
c                             absolutely greatest element of matrix a.
c
c        remarks
c           input matrices r and a are assumed to be stored columnwise
c           in m*n resp. m*m successive storage locations. on return
c           solution matrix r is stored columnwise too.
c           the procedure gives results if the number of equations m is
c           greater than 0 and pivot elements at all elimination steps
c           are different from 0. however warning ier=k - if given -
c           indicates possible loss of significance. in case of a well
c           scaled matrix a and appropriate tolerance eps, ier=k may be
c           interpreted that matrix a has the rank k. no warning is
c           given in case m=1.
c
c        subroutines and function subprograms required
c           none
c
c        method
c           solution is done by means of gauss-elimination with
c           complete pivoting.
c
c     ..................................................................
c
      subroutine dgelg (r,a,m,n,eps,ier)
c
      dimension a(1), r(1)
      double precision r, a, piv, tb, tol, pivi
      if (m) 230,230,10
c
c     search for greatest element in matrix a
   10 ier=0
      piv=0.d0
      mm=m*m
      nm=n*m
      do 30 l=1,mm
      tb=dabs(a(l))
      if (tb-piv) 30,30,20
   20 piv=tb
      i=l
   30 continue
      tol=eps*piv
c     a(i) is pivot element. piv contains the absolute value of a(i).
c
c     start elimination loop
      lst=1
      do 170 k=1,m
c
c     test on singularity
      if (piv) 230,230,40
   40 if (ier) 70,50,70
   50 if (piv-tol) 60,60,70
   60 ier=k-1
   70 pivi=1.d0/a(i)
      j=(i-1)/m
      i=i-j*m-k
      j=j+1-k
c     i+k is row-index, j+k column-index of pivot element
c
c     pivot row reduction and row interchange in right hand side r
      do 80 l=k,nm,m
      ll=l+i
      tb=pivi*r(ll)
      r(ll)=r(l)
   80 r(l)=tb
c
c     is elimination terminated
      if (k-m) 90,180,180
c
c     column interchange in matrix a
   90 lend=lst+m-k
      if (j) 120,120,100
  100 ii=j*m
      do 110 l=lst,lend
      tb=a(l)
      ll=l+ii
      a(l)=a(ll)
  110 a(ll)=tb
c
c     row interchange and pivot row reduction in matrix a
  120 do 130 l=lst,mm,m
      ll=l+i
      tb=pivi*a(ll)
      a(ll)=a(l)
  130 a(l)=tb
c
c     save column interchange information
      a(lst)=j
c
c     element reduction and next pivot search
      piv=0.d0
      lst=lst+1
      j=0
      do 160 ii=lst,lend
      pivi=-a(ii)
      ist=ii+m
      j=j+1
      do 150 l=ist,mm,m
      ll=l-j
      a(l)=a(l)+pivi*a(ll)
      tb=dabs(a(l))
      if (tb-piv) 150,150,140
  140 piv=tb
      i=l
  150 continue
      do 160 l=k,nm,m
      ll=l+j
  160 r(ll)=r(ll)+pivi*r(l)
  170 lst=lst+m
c     end of elimination loop
c
c     back substitution and back interchange
  180 if (m-1) 230,220,190
  190 ist=mm+m
      lst=m+1
      do 210 i=2,m
      ii=lst-i
      ist=ist-lst
      l=ist-m
      l=a(l)+.5d0
      do 210 j=ii,nm,m
      tb=r(j)
      ll=j
      do 200 k=ist,mm,m
      ll=ll+1
  200 tb=tb-a(k)*r(ll)
      k=j+l
      r(j)=r(k)
  210 r(k)=tb
  220 return
c
c     error return
  230 ier=-1
      return
      end