      subroutine inv(npmax,ncol,nfil,ninv,u,v,w,d,x,alpha,cov)
c
c     performs backsubstitution, given svd decomposition
c
c     single precision version
c
      parameter (ncolm=550)
c
      real*4 u(npmax,ncol), v(ncol,ncol), w(ncol)
      real*4 cov(ncol,ncol)
      real*4 d(npmax), x(ncol)
      real*4 wmax, wmin, s, alpha
      real*4 w2(ncolm), temp(ncolm)
c
c     check array bounds
c
      if (ncol .gt. ncolm) then
        print*,' array bound ncolm exceeded in inv'
      end if
c
c     zero the x and temp arrays for good measure
c
      do i = 1, ncol
        x(i) = 0.0
        temp(i) = 0.0
      end do
c
c     set threshold for smallest eigenvalues
c
      wmax = 0.0
      do j = 1, ninv
        if (w(j) .gt. wmax) then 
          wmax = w(j)
        endif
      end do
      wmin = wmax*1.0e-5
      do j = 1, ninv
        if (w(j) .lt. wmin) then 
          print*,' zeroing a column # ',j,w(j)
          w(j) = 0.0
        endif
      end do
c
c     perform back substitution
c
      do j = 1, ninv
        s = 0.0
        if (w(j) .ne. 0.0) then
          do i = 1, nfil
            s = s + u(i,j)*d(i)
          end do
          s = s*(w(j)/(w(j)**2 + alpha))
        endif
        temp(j) = s
      end do
      do j = 1, ninv
        s = 0.0
        do jj = 1, ninv
          s = s + v(j,jj)*temp(jj)
        end do
        x(j) = s 
      end do
c       
c     estimate the variance associated with this model
c
c     form the square of the singular values
c
      do j = 1, ninv
        w2(j) = 0.d0
        if ((w(j) + alpha) .ne. 0.0) then
          s = w(j)**2
          w2(j) = s / (s + alpha)**2
        endif
      end do
c
c     now form the variance matrix
c
      do ii = 1, ninv
        do jj = 1, ii
          s = 0.0
          do kk = 1, ninv
            s = s + v(ii,kk)*v(jj,kk)*w2(kk)
          end do
          cov(ii,jj) = s
          cov(jj,ii) = s
        end do
      end do
c
      return
      end
