      subroutine lfit (nord, ndat, xx, dd, sigma, lsig, mm, ss, chi2)
c
c     subroutine to perform a general linear least squares fit
c     order of fit determined by user
c
c     nord = order (2 = line, etc)
c     ndat = number of points in inversion
c     xx = array of basis function
c     dd = array of data
c     sigma = array of errors
c     lsig = true if the errrors are independent
c     mm = estimated parameters
c     ss = errors
c     chi2 - misfit
c
      parameter (npmax = 65536)
      parameter (ncol = 5)
c
      real*4 xx(*), dd(*), sigma(*), mm(*), ss(*)
      real*4 am(npmax,ncol), u(npmax,ncol), v(ncol,ncol)
      real*4 w(ncol), cov(ncol,ncol)
      real*4 data(npmax), xm(ncol), sum, chi2
c
      logical lsig
c
c     index controls the svd operation
c
      index = 1
c
c     alpha is the ls damping parameter
c
      alpha = 0.0
c
      if (nord .gt. ncol) then
        print*,' parameter ncol exceeded in lfit'
        return
      elseif (ndat .gt. npmax) then
        print*,' parameter npmax exceeded in lfit'
        return
      endif
c
      nfil = ndat
c
c     form the A or design matrix - weighted by the input errors
c
      if (nord .gt. 0) then
        ninv = nord
        do ii = 1, nfil
          s = 1.0/sigma(ii)
          do jj = 1, nord
            am(ii,jj) = s*(xx(ii)**(jj-1))
          end do
          data(ii) = s*dd(ii)
cccc
c         print*, ii, s, am(ii,1), am(ii,2), data(ii)
cccc
        end do 
      else
        ninv = 1
        do ii = 1, nfil
          s = 1.0/sigma(ii)
          am(ii,1) = s*(xx(ii)**(abs(nord)-1))
          data(ii) = s*dd(ii)
        end do
      endif
c
c     compute singular value decomposition 
c
      call svd1(npmax,ncol,nfil,ninv,am,u,v,w,index)
c
c     preform backsubstitution
c
      call inv(npmax,ncol,nfil,ninv,u,v,w,data,xm,alpha,cov)
c
c     evaluate the chi-square
c
      chi2 = 0.0
      if (nord .gt. 0) then
        do ii = 1, nfil
          sum = 0.0
          do jj = 1, ninv
            sum = sum + xm(jj)*(xx(ii)**(jj-1))
          end do
          chi2 = chi2 + ((data(ii) - sum)/sigma(ii))**2
        end do
      else
        do ii = 1, nfil
          ax = xm(1)*(xx(ii)**(abs(nord)-1))
          chi2 = chi2 + ((data(ii) - ax)/sigma(ii))**2
        end do
      end if
c
c     test whether the errors are independent
c
      if (lsig) then 
        do ii = 1, ninv
          ss(ii) = sqrt(cov(ii,ii))
          mm(ii) = xm(ii)
        end do
      else
        scale = sqrt(chi2/(nfil-ninv))
        do ii = 1, ninv
          ss(ii) = sqrt(cov(ii,ii))*scale
          mm(ii) = xm(ii)
        end do
      endif
c
c     all done! - go home!
c
      return
      end
