	subroutine reduce(h,z1,noref)
	implicit none
	real*4 z1(3,3),h(3,3),z2(3,3),t1(3),temp,sgn
	integer*4 i,j,noref,nminus
c
	do 99 i=1,3
	  do 98 j=1,3
	    z2(j,i)=0.0
98	  continue
99	continue
100	if(z1(1,1).gt.z1(2,2).or.(z1(1,1).eq.z1(2,2).and.
	1 abs(z1(3,2)).gt.abs(z1(3,1))))go to 102
	go to 110
c		case 1 - S11 > S22 or S11 = S22 and abs(S23) > abs(S31)
102	  z2(1,2)=1.0
	  z2(2,1)=1.0
	  z2(3,3)=-1.0
	  call transf(z1,z2,h)
	  go to 100
c	end if
110	if(z1(2,2).gt.z1(3,3).or.(z1(2,2).eq.z1(3,3).and.
	1 abs(z1(3,1)).gt.abs(z1(2,1))))go to 112
	go to 120
c		case 2 - S22 > S33 or S22 = S33 and abs(S31) > abs(S12)
112	  z2(1,1)=-1.0
	  z2(2,3)=1.0
	  z2(3,2)=1.0
	  call transf(z1,z2,h)
	  go to 100
c	end if
120	do i=1,3
	    z2(i,i)=1.0
	enddo
	if(z1(3,2)/sqrt(z1(3,3)*z1(2,2)).gt.1.0e-3)then
	      z2(3,3)=-z2(3,3)
	      z2(2,2)=-z2(2,2)
	endif
	if(z1(3,1)/sqrt(z1(3,3)*z1(1,1)).gt.1.0e-3)then
		z2(3,3)=-z2(3,3)
		z2(1,1)=-z2(1,1)
	endif
	if(z1(2,1)/sqrt(z1(1,1)*z1(2,2)).gt.1.0e-3)then
		z2(1,1)=-z2(1,1)
		z2(2,2)=-z2(2,2)
	endif
	if(abs(z2(1,1)+z2(2,2)+z2(3,3)).ne.3.0)then
c		case 3,4 - All Off-Diagonal Terms Must be < 0 or > 0
   	    call transf(z1,z2,h)
	end if
130	z2(1,1)=1.0
	z2(2,2)=1.0
	z2(3,3)=1.0
	if(2.0*abs(z1(3,2)).gt.z1(2,2).or.
	1  (2.0*z1(3,2).eq.z1(2,2).and.2.0*z1(3,1).lt.z1(2,1)).or.
	2  (2.0*z1(3,2).eq.-z1(2,2).and.z1(2,1).lt.0.0))go to 132
	go to 140
c		case 5
132	   z2(3,2)=-sgn(z1(3,2)/sqrt(z1(3,3)*z1(2,2)))
	   call transf(z1,z2,h)
	   go to 100
c	end if
140	if(2.0*abs(z1(3,1)).gt.z1(1,1).or.
	1  (2.0*z1(3,1).eq.z1(1,1).and.2.0*z1(3,2).lt.z1(2,1)).or.
	2  (2.0*z1(3,1).eq.-z1(1,1).and.z1(2,1).lt.0.0))go to 142
	go to 150
c		case 6
142	   z2(3,1)=-sgn(z1(3,1)/sqrt(z1(1,1)*z1(3,3)))
	   call transf(z1,z2,h)
	   go to 100
c	end if
150	if(2.0*abs(z1(2,1)).gt.z1(1,1).or.
	1  (2.0*z1(2,1).eq.z1(1,1).and.2.0*z1(3,2).lt.z1(3,1)).or.
	2  (2.0*z1(2,1).eq.-z1(1,1).and.z1(2,1).lt.0.0))go to 152
	go to 160
c		case 7
152	   z2(2,1)=-sgn(z1(2,1)/sqrt(z1(1,1)*z1(2,2)))
	   call transf(z1,z2,h)
	   go to 100
c	end if
160	temp=2.0*(z1(2,1)+z1(3,1)+z1(3,2))+z1(1,1)+z1(2,2)
	if(temp.lt.0.0.or.(temp.eq.0.0.and.2.0*(z1(1,1)+2.0*z1(3,1))
	1+2.0*z1(2,1).gt.0.0))go to 162
	go to 170
c		case 8
162	  z2(1,1)=1.0
	  z2(2,2)=1.0
	  do 164 i=1,3
	    z2(3,i)=1.0
164	  continue
	  call transf(z1,z2,h)
	  go to 100
c	end if
170	do i=1,3
	    z2(i,i)=1.0
	enddo
	nminus=0
	if(sgn(z1(3,2)/sqrt(z1(3,3)*z1(2,2))).lt.0.0)nminus=nminus+1
	if(sgn(z1(3,1)/sqrt(z1(3,3)*z1(1,1))).lt.0.0)nminus=nminus+1
	if(sgn(z1(2,1)/sqrt(z1(1,1)*z1(2,2))).lt.0.0)nminus=nminus+1
	if(nminus.eq.0.or.nminus.eq.3)return
	if(sgn(z1(3,2)/sqrt(z1(3,3)*z1(2,2))).gt.0.0)then
	      z2(1,1)=-z2(1,1)
	      z2(2,2)=-z2(2,2)
	endif
	if(sgn(z1(3,1)/sqrt(z1(3,3)*z1(1,1))).gt.0.0)then
		z2(3,3)=-z2(3,3)
		z2(2,2)=-z2(2,2)
	endif
	if(sgn(z1(2,1)/sqrt(z1(1,1)*z1(2,2))).gt.0.0)then
		z2(1,1)=-z2(1,1)
		z2(3,3)=-z2(3,3)
	endif
	if(abs(z2(1,1)+z2(2,2)+z2(3,3)).ne.3.0)then
c		case 3,4 - All Off-Diagonal Terms Must be < 0 or > 0
   	    call transf(z1,z2,h)
	end if
	return
	end
	subroutine transf(z1,z2,b)
	dimension z1(3,3),z2(3,3),b(3,3),t1(3,3)
c
	do 120 i=1,3
	  do 110 j=1,3
	    t1(i,j)=0.0
	    do 100 k=1,3
		t1(i,j)=t1(i,j)+z2(i,k)*z1(k,j)
100	    continue
110	  continue
120	continue
	do 150 i=1,3
	  do 140 j=1,3
	    z1(i,j)=0.0
	    do 130 k=1,3
		z1(i,j)=z1(i,j)+t1(i,k)*z2(j,k)
130	    continue
140	  continue
150	continue
	do 180 i=1,3
	  do 170 j=1,3
	    t1(i,j)=0.0
	    do 160 k=1,3
		t1(i,j)=t1(i,j)+z2(i,k)*b(k,j)
160	    continue
170	  continue
180	continue
	do 200 i=1,3
	  do 190 j=1,3
	    z2(i,j)=0.0
	    b(i,j)=t1(i,j)
190	  continue
200	continue
	return
	end
