c     ..................................................................
c
c        subroutine dgels
c
c        purpose
c           to solve a system of simultaneous linear equations with
c           symmetric coefficient matrix upper triangular part of which
c           is assumed to be stored columnwise.
c
c        usage
c           call dgels(r,a,m,n,eps,ier,aux)
c
c        description of parameters
c           r      - double precision m by n right hand side matrix
c                    (destroyed). on return r contains the solution of
c                    the equations.
c           a      - upper triangular part of the symmetric double
c                    precision m by m coefficient matrix.  (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 main diagonal
c                             element of matrix a.
c           aux    - double precision auxiliary storage array
c                    with dimension m-1.
c
c        remarks
c           upper triangular part of matrix a is assumed to be stored
c           columnwise in m*(m+1)/2 successive storage locations, right
c           hand side matrix r columnwise in n*m successive storage
c           locations. on return solution matrix r is stored columnwise
c           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           error parameter ier=-1 does not necessarily mean that
c           matrix a is singular, as only main diagonal elements
c           are used as pivot elements. possibly subroutine dgelg (which
c           works with total pivoting) would be able to find a solution.
c
c        subroutines and function subprograms required
c           none
c
c        method
c           solution is done by means of gauss-elimination with
c           pivoting in main diagonal, in order to preserve
c           symmetry in remaining coefficient matrices.
c
c     ..................................................................
c
      subroutine dgels (r,a,m,n,eps,ier,aux)
c
      dimension a(1), r(1), aux(1)
      double precision r, a, aux, piv, tb, tol, pivi
      if (m) 240,240,10
c
c     search for greatest main diagonal element
   10 ier=0
      piv=0.d0
      l=0
      do 30 k=1,m
      l=l+k
      tb=dabs(a(l))
      if (tb-piv) 30,30,20
   20 piv=tb
      i=l
      j=k
   30 continue
      tol=eps*piv
c     main diagonal element a(i)=a(j,j) is first pivot element.
c     piv contains the absolute value of a(i).
c
c     start elimination loop
      lst=0
      nm=n*m
      lend=m-1
      do 180 k=1,m
c
c     test on usefulness of symmetric algorithm
      if (piv) 240,240,40
   40 if (ier) 70,50,70
   50 if (piv-tol) 60,60,70
   60 ier=k-1
   70 lt=j-k
      lst=lst+k
c
c     pivot row reduction and row interchange in right hand side r
      pivi=1.d0/a(i)
      do 80 l=k,nm,m
      ll=l+lt
      tb=pivi*r(ll)
      r(ll)=r(l)
   80 r(l)=tb
c
c     is elimination terminated
      if (k-m) 90,190,190
c
c     row and column interchange and pivot row reduction in matrix a.
c     elements of pivot column are saved in auxiliary vector aux.
   90 lr=lst+(lt*(k+j-1))/2
      ll=lr
      l=lst
      do 140 ii=k,lend
      l=l+ii
      ll=ll+1
      if (l-lr) 120,100,110
  100 a(ll)=a(lst)
      tb=a(l)
      go to 130
  110 ll=l+lt
  120 tb=a(ll)
      a(ll)=a(l)
  130 aux(ii)=tb
  140 a(l)=pivi*tb
c
c     save column interchange information
      a(lst)=lt
c
c     element reduction and search for next pivot
      piv=0.d0
      llst=lst
      lt=0
      do 180 ii=k,lend
      pivi=-aux(ii)
      ll=llst
      lt=lt+1
      do 150 lld=ii,lend
      ll=ll+lld
      l=ll+lt
  150 a(l)=a(l)+pivi*a(ll)
      llst=llst+ii
      lr=llst+lt
      tb=dabs(a(lr))
      if (tb-piv) 170,170,160
  160 piv=tb
      i=lr
      j=ii+1
  170 do 180 lr=k,nm,m
      ll=lr+lt
  180 r(ll)=r(ll)+pivi*r(lr)
c     end of elimination loop
c
c     back substitution and back interchange
  190 if (lend) 240,230,200
  200 ii=m
      do 220 i=2,m
      lst=lst-ii
      ii=ii-1
      l=a(lst)+.5d0
      do 220 j=ii,nm,m
      tb=r(j)
      ll=j
      k=lst
      do 210 lt=ii,lend
      ll=ll+1
      k=k+lt
  210 tb=tb-a(k)*r(ll)
      k=j+l
      r(j)=r(k)
  220 r(k)=tb
  230 return
c
c     error return
  240 ier=-1
      return
      end