!     Last change:  DOS  26 Jul 2000    8:45 pm
!                                       *** copyright 2000 ***
! *** filename sweep.f95                *** John F. Monahan **
!                                       **********************
      program psweep
!  test of sweep on Brownlee data, Exercise 13.1, p462
      implicit none
      INTEGER, parameter :: p = 5
      INTEGER, parameter :: n = 22
      REAL, DIMENSION(p,p) :: xx
      REAL, DIMENSION(p) :: row
      integer i,j,k,kk
!
      interface
      subroutine sweep(a,m,k)
      real, dimension(:,:), intent(in out) :: a
      integer, intent(in) :: m,k
      end subroutine sweep
      end interface
!
  22  format(14x,f5.3,f5.0,f4.0,f6.2)
  23  format(5f15.8)
  24  format(' after step kk=',i2,'  sweeping k=',i2)      
!
      open( unit=5, file='brown463.dat')
      OPEN( UNIT=6, FILE='sweep.out' )
!
      xx = 0.
!                   read in the Brownlee data, Exercise 13.1
!
!                   last variable is response (Y)
      do i = 1,n
      row(1) = 1.
      read(5,22) row(2),row(3),row(4),row(5)
!                   form the X'X matrix
      do j = 1,p
      do k = 1,p
      xx(j,k) = xx(j,k) + row(j)*row(k)
      end do  ! loop on k
      end do  ! loop on j
!
      end do  ! loop i
!                   write out X'X
      write(6,24) 0,0
      do  i = 1,p
      write(6,23) (xx(i,j),j=1,p)
      end do  ! loop on i
!                   sweep in variable k and print it all out
      do kk = 1,8
      k = 1 + mod(kk-1,4) 
!                   doing it twice, second sweep deletes variable
      write(6,24) kk,k  
      call sweep(xx,p,k)
      do i = 1,p
      write(6,23) (xx(i,j),j=1,p)
      end do  ! loop on i
      end do  ! loop on kk
      stop
      end program psweep
      subroutine sweep(a,m,k)
! sweep operator on the sum of squares and cross products mx A
!  a  real matrix m by m  (dimensioned 10 by 10 for convenience)
!  m  integer size of matrix p+1 for one dependent and p regressors
!  k  integer variable to be swept
!
      implicit none
      real, dimension(:,:), intent(in out) :: a
      integer, intent(in) :: m,k
      real b,d
      integer i,j
!                   sweeping in/out variable k
      d = a(k,k)
      do j = 1,m
      a(k,j) = a(k,j)/d
      end do  ! loop on j
!                   fixed row k, now rest of matrix
      do i = 1,m
      if( i .ne. k ) then
      b = a(i,k) 
      do j = 1,m
      a(i,j) = a(i,j) - b*a(k,j)
      end do  ! loop on j
      a(i,k) = -b/d
                     END if
      end do  ! loop on i
      a(k,k) = 1.d0/d 
!                   done
      return
      end subroutine sweep
! *** end of filename sweep.f95         *********************