c....*...1.........2.........3.........4.........5.........6.........7.*.......8
c     regr9    1/11/00
c
c     purpose
c     multiple regression
c
c     usage
c     call regr9(y,x,b,e,m,n,k,g,s,r,eps,ier)
c
c     arguments
c     y   - an n by m input vector stored columnwise (storage mode of 0)
c           elements of y are real*8
c     x   - an n by k input matrix stored columnwise (storage mode of 0)
c           elements of x are real*8
c     b   - a  k by m vector of regression coefficients stored columnwise 
c           (storage mode of 0)
c           elements of b are real*8
c     e   - an n by m vector of residuals stored columnwise 
c           (storage mode of 0)
c           elements of e are real*8
c     m   - number of columns in y and b 
c           integer
c     n   - number of rows in x, y, and e
c           integer
c     k   - number of columns in x, rows of b, number of rows and
c           columns in g.  n must be .ge. k.
c           integer
c     g   - a k by k matrix containing the inverse of x'x stored 
c           columnwise  (storage mode of 0)
c           elements of g are real*8
c     s   - estimated m by m variance matrix on n-k+ier degrees freedom 
c           stored columnwise (storage mode of 0)
c           elements of s are real*8
c     r   - vector of length m containing the uncorrected R squared
c           elements of r are real*8
c     eps - input constant used as a relative tolerince in testing for 
c           degenerate rank.  a resonable value for eps is 1.d-13.
c           real*8
c     ier - error parameter coded as follows:
c           ier=0  no error.
c           ier>0  x is not full rank (rank of x is k-ier).
c           integer
c
c     subroutines called
c     dgmapa, dgmabp, dsweep, dgmprd, dgmsub
c
c     remarks
c     on return:
c     b=inverse(x'x)x'y
c     e=y-xb
c     s=sum(e'e)/(n-k)
c     g=inverse(x'x)
c  
c     this routine just does the indicated algebra and is therefore
c     not numerically stable
c
      subroutine regr9(y,x,b,e,m,n,k,g,s,r,eps,ier)
      implicit none
      integer*4 m,n,k,ier
      real*8 y(n,m),x(n,k),b(k,m),e(n,m),s(m,m),g(k,k),r(m),eps
      real*8 df
      integer*4 i,j
      call dgmapa(y,s,n,m)
      do i=1,m
        r(i)=s(i,i)
      end do
      call dgmapa(x,g,n,k)
      call dgmapb(x,y,e,n,k,m)
      call dsweep(g,k,eps,ier)
      call dgmprd(g,e,b,k,k,m)
      call dgmprd(x,b,e,n,k,m)
      call dgmsub(y,e,e,n,m)
      call dgmapa(e,s,n,m)
      do i=1,m
        r(i)=1.d0-s(i,i)/r(i)
      end do
      df=n-k+ier
      do i=1,m
      do j=1,m
        s(i,j)=s(i,j)/df
      end do
      end do
      end