!     Last change:  DOS  26 Jul 2000    8:38 pm
!                                       *** copyright 2000 ***
! *** filename chex57.f95               *** John F. Monahan **
!                                       **********************
      program chex57
!  check Example 5.7 using ROT534 to do Givens' transformations
      implicit none
      integer, parameter :: n = 4
      integer, parameter :: p = 2
      real, dimension(n,p) :: x
      real, dimension(n,n) :: q         ! orthogonal matrix
      real a,b,c,s,rnu,y,z
      integer i,j,ii,jj,m
!
  24  format(/'cosine and sine parts of rotation',2f8.4)
  25  format(/'results after rotating in rows',i2,' and',i2     &
     &       /10x,'X',38x,'Q')
  26  format(2f10.6,4x,4f12.6)
!
!  use Example 5.7, simple linear regression (p=2) with n=4 obs
      data x/ 5*1., 2., 3., 4. /
!  start with identity matrix for Q
      data q/ 1., 4*0., 1., 4*0., 1., 4*0., 1./
!
      open( unit=6, file='chex57.out' )
!
!  use the second method of looping -- by observations
!
      do j = 2,n              ! zero in this row
      m = min( j-1, p )
      do i = 1,m            ! put ss in x(i,i)
!                            select elements
      a = x(i,i)
      b = x(j,i)
!                            find Givens' tfm
      call ROT734(a,b,c,s,rnu)
      write(6,24) c,s
!                   apply transformations to X and Q
!                   first to X
      do jj = 1,p
      y = c*x(i,jj) + s*x(j,jj)
      z = c*x(j,jj) - s*x(i,jj)
      x(i,jj) = y
      x(j,jj) = z
      end do  ! loop on jj
!                   now to Q
      do jj = 1,n
      y = c*q(i,jj) + s*q(j,jj)
      z = c*q(j,jj) - s*q(i,jj)
      q(i,jj) = y
      q(j,jj) = z
      end do  ! loop on jj
!  write out current x and q
      write(6,25) i,j
      do ii = 1,n
      write(6,26) (x(ii,jj),jj=1,p),(q(ii,jj),jj=1,n)
      end do  ! loop on jj
!  all done with changing rows i and j
      end do  ! loop on i
      end do  ! loop on j
      stop
      end program chex57
      SUBROUTINE ROT734(ALPHA,BETA,GAMA,SIGMA,RNU)
!  GIVEN ALPHA AND BETA, PRODUCES GIVENS TRANSFORMATION VALUES
!
!      (  GAMA       SIGMA  )  (  ALPHA  )       (  RNU  )
!      (                    )  (         )   =   (       )
!      ( -SIGMA       GAMA  )  (  BETA   )       ( ZERO  )
!
!  ALGORITHM 7.3.4 FROM STEWART
!
      IMPLICIT NONE
      REAL, INTENT(IN) :: ALPHA,BETA
      REAL, INTENT(OUT) :: GAMA,SIGMA,RNU
      REAL ETA,DELTA
!
      ETA = MAX( ABS(ALPHA), ABS(BETA) )
!                     IF BOTH ALPHA AND BETA ARE ZERO THEN
!                     USE IDENTITY TRANSFORMATION
      IF( ETA .EQ. 0. ) then
!                     SET UP IDENTITY TRANSFORMATION IF NO WORK DONE
      RNU = 0.
      GAMA = 1.                                                         
      SIGMA = 0.                                                        
      RETURN
                        else
      DELTA = (ALPHA/ETA)*(ALPHA/ETA) + (BETA/ETA)*(BETA/ETA)           
      DELTA = SQRT(DELTA)                                               
      RNU = ETA * DELTA                                                 
      GAMA = ALPHA / RNU                                                
      SIGMA = BETA / RNU
                        END if
      RETURN                                                            
      END SUBROUTINE ROT734
! *** end of filename chex57.f95        *********************