! 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 *********************