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