subroutine qr(m, n, a, b, x, resq)
      implicit double precision (a-h, o-z)                              *doub
c$$$$ calls no other routines
c  Relies on FORTRAN77 do-loop conventions!
c  Solves over-determined least-squares problem  ax ~ b
c  where  a  is an  m by n  matrix,  b  is an m-vector .
c  resq  is the sum of squared residuals of optimal solution.  Also used
c  to signal error conditions - if -2 , system is underdetermined,  if
c  -1,  system is singular.
c  Method - successive Householder rotations.  See Lawson & Hanson - 
c  Solving Least Squares Problems (1974).
c  Routine will also work when m=n.
c*****   CAUTION -  a and b  are overwritten by this routine.
      dimension a(m,n),b(m),x(n)
      double precision sum, dot
c
      resq=-2.0
      if (m .lt. n) return
      resq=-1.0
c   Loop ending on 1800 rotates  a  into upper triangular form.
      do 1800 j=1, n
c   Find constants for rotation and diagonal entry.
        sq=0.0
        do 1100 i=j, m
          sq=a(i,j)**2 + sq
 1100   continue
        if (sq .eq. 0.0) return
        qv1=-sign(sqrt(sq), a(j,j))
        u1=a(j,j) - qv1
        a(j,j)=qv1
        j1=j + 1
c  Rotate remaining columns of sub-matrix.
        do 1400 jj=j1, n
          dot=u1*a(j,jj)
          do 1200 i=j1, m
            dot=a(i,jj)*a(i,j) + dot
 1200     continue
          const=dot/abs(qv1*u1)
          do 1300 i=j1, m
            a(i,jj)=a(i,jj) - const*a(i,j)
 1300     continue
          a(j,jj)=a(j,jj) - const*u1
 1400   continue
c  Rotate  b  vector.
        dot=u1*b(j)
        do 1600 i=j1, m
          dot=b(i)*a(i,j) + dot
 1600   continue
        const=dot/abs(qv1*u1)
        b(j)=b(j) - const*u1
        do 1700 i=j1, m
          b(i)=b(i) - const*a(i,j)
 1700   continue
 1800 continue
c  Solve triangular system by back-substitution.
      do 2200 ii=1, n
        i=n-ii+1
        sum=b(i)
        do 2100 j=i+1, n
          sum=sum - a(i,j)*x(j)
 2100   continue
        if (a(i,i).eq. 0.0) return
         x(i)=sum/a(i,i)
 2200 continue
c  Find residual in overdetermined case.
      resq=0.0
      do 2300 i=n+1, m
        resq=b(i)**2 + resq
 2300 continue
      return
      end                                                               qr
c______________________________________________________________________