c
c
c***********************************************************************
c
      subroutine ginv ( a, mr, nr, nc, tau, u, atemp )
c
c-----------------------------------------------------------------------
c
c     upon entry   a = a matrix with nr rows and nc columns 
c     after the return the original a is destroyed and the array a
c     contains the transpose of  t   the transformation matrix
c        (x = t*b) which solves the pair of equations
c
c                  a*x = b    and    tau i*x = 0
c
c          in the least square sense.
c
c     mr = 1st dimension no. of array a in the calling programs
c     tau is a nonnegative constant which controls the error
c     propigation of the transformation matrix.   if tau = 0.   then the
c     resulting transformation matrix is the ordinary least squares
c     transformation matrix.   (the inverse if nr = nc)
c     atemp and u are used for working space by the algorithm and do not
c     necessarily contain any relavant numbers at the conclusion.
c     atemp must be dimensioned at least nc by the calling programs
c     u must be dimensioned at least nc*nc by the calling programs
c     no. of multiplications = nc**2 (5/2 nr + 2/3 nc)
c
c-----------------------------------------------------------------------
c
      dimension a(mr,nc),u(nc,nc),atemp(nc)
c
c------------------------------------------------------------------------
c
      tausq = tau**2
c
c     place unit matrix in u
c
      do  5  i = 1,nc
          do  4  j = 1,nc
    4         u(i,j) = 0.0
    5     u(i,i) = 1.0
c
c     orthogonalize combined matrix (a above u) by gramm-schmidt-hilbert
c        method with first nr rows weighted with 1  and the other
c        nc rows weighted with 1/tau.  then reorthogonalize to
c        lessen roundoff error.
c
      do  20  i = 1,nc
          ii = i - 1
          if  (ii)  2,11,2
    2     do  10  ll = 1,2
              do  10  j = 1,ii
                  dot = 0.0
                  dot2 = 0.0
                  do  3  k = 1,nr
    3                 dot = a(k,i)*a(k,j) + dot
                  do  6  k = 1,j
    6                 dot2 = u(k,i)*u(k,j) + dot2 
                  dot = dot + dot2*tausq
                  do  8  k = 1,j
    8                 u(k,i) = u(k,i) - dot*u(k,j)
                  do  10  k = 1,nr
   10                 a(k,i) = a(k,i) - dot*a(k,j)
c
c     normalize the column  i of the combined matrix
c
   11     dot = 0.0 
          dot2 = 0.0
          do  12  k = 1,nr
   12         dot = a(k,i)*a(k,i) + dot 
          do  14  k = 1,i
   14         dot2 = u(k,i)*u(k,i) + dot2
          dot = dot + dot2*tausq
          dot = sqrt(dot)
          do 17  k = 1,i
   17         u(k,i) = u(k,i)/dot
          do  19  k = 1,nr
   19         a(k,i) = a(k,i)/dot
   20 continue
c
c     calculation of the transpose of the transformation
c        matrix   t(trans.) = a * u(tr.)
c
      do  50  i = 1,nr
          do  45  j = 1,nc
              atemp(j) = 0.0
              do  45  k = j,nc
                  atemp(j) = a(i,k)*u(j,k) + atemp(j)
   45         continue
          do  50  j = 1,nc
              a(i,j) = atemp(j)
   50         continue
c
      return
c
      end 
