IMPLICIT NONE
	INTEGER NMAX,N,I,J,P,Q
	PARAMETER(NMAX=50)
	REAL A(NMAX,NMAX),AMAX,TEMP,ZEMP,COO,SII,CO,SI,APP,AQQ,APQ,API,AQI
	REAL R(NMAX,NMAX),RIP,RIQ
	CHARACTER NAME*12,NAMEO*12,CHR*1
 ! 從文件中讀入實對稱矩陣A
  
	WRITE(*,*)'輸入實對稱矩陣維數n(n<51):'
	READ(*,*) N
	WRITE(*,*)'輸入矩陣文件:'
	READ(*,*) NAME
	OPEN(6,FILE=NAME)
	DO I=1,N
	 READ(6,*) (A(I,J),J=1,I)
	 DO J=1,I
	 A(J,I)=A(I,J)
	 ENDDO
	ENDDO
	CLOSE(6)
	!R矩陣中存放正交變換矩陣U，在這先初始化，即單位矩陣				  
	DO I=1,N
	DO J=1,N
	R(I,J)=0
	ENDDO
	R(I,I)=1
	ENDDO
	!矩陣A非主對角線元素中,找出按模最大的元素Apq
100	AMAX=ABS(A(2,1))
	P=2
	Q=1
	DO I=2,N
	 DO J=1,I-1
	 IF(ABS(A(I,J)).GT.AMAX) THEN
	 AMAX=ABS(A(I,J))
	 P=I
	 Q=J
	 ENDIF
	 ENDDO
	ENDDO

	do i=1,n
		WRITE(*,*)(A(I,J),J=1,N)
	enddo
!當非主對角元素化為0,及小于給定精度時, 輸出特征值與特征向量
	IF(AMAX.LE.1.0E-7) THEN
	WRITE(*,*) 'A的特征值為:'
	WRITE(*,*) (A(I,I),I=1,N)
	WRITE(*,*) 'A的特征向量為:'
	WRITE(*,*) '  X1      X2     X3     ...:'
	DO I=1,N
	WRITE(*,*)(R(I,J),J=1,N)
	ENDDO

	WRITE(*,*) '是否將結果存入文件(Y/N)?'
	READ(*,*) CHR
	IF(CHR.EQ.'Y'.OR.CHR.EQ.'y') THEN
	WRITE(*,*) '輸入文件名(小于12字符):'
	READ(*,*) NAMEO
	OPEN(8,FILE=NAMEO)
	WRITE(8,*) 'A的特征值為:'
	WRITE(8,*) (A(I,I),I=1,N)
	WRITE(8,*) 'A的特征向量為:'
	WRITE(8,*) '  X1      X2     X3     ...:'
	DO I=1,N
	WRITE(8,*)(R(I,J),J=1,N)
	ENDDO
	CLOSE(8)
	ENDIF

	STOP
	ENDIF
!開始準備計算平面旋轉矩陣 u
	TEMP=2*A(P,Q)/(A(P,P)-A(Q,Q)+1.0e-30)
	ZEMP=(A(P,P)-A(Q,Q))/(2*A(P,Q))
	
	IF(ABS(TEMP).LT.1.0) THEN
	COO=(1+TEMP**2)**(-0.5)
	SII=TEMP*(1+TEMP**2)**(-0.5)
	ELSE
	COO=ABS(ZEMP)*(1+ZEMP**2)**(-0.5)
	SII=SIGN(1.0,ZEMP)*(1+ZEMP**2)**(-0.5)
	ENDIF
	
	CO=SQRT(0.5*(1+COO))
	SI=SII/(2*CO)
! 計算平面旋轉矩陣U
	DO I=1,N
	RIP=R(I,P)*CO+R(I,Q)*SI
	RIQ=-R(I,P)*SI+R(I,Q)*CO
	R(I,P)=RIP
	R(I,Q)=RIQ
	ENDDO
!對A進行變換	
	APP=A(P,P)*CO**2+A(Q,Q)*SI**2+2*A(P,Q)*CO*SI
	AQQ=A(P,P)*SI**2+A(Q,Q)*CO**2-2*A(P,Q)*CO*SI
	APQ=0.5*(A(Q,Q)-A(P,P))*SII+A(P,Q)*COO
	A(P,P)=APP
	A(Q,Q)=AQQ
	A(P,Q)=APQ
	A(Q,P)=A(P,Q)
	
	DO I=1,N
	IF(I.EQ.P.OR.I.EQ.Q) THEN
	ELSE
	API=A(P,I)*CO+A(Q,I)*SI
	AQI=-A(P,I)*SI+A(Q,I)*CO
	A(P,I)=API
	A(Q,I)=AQI
	A(I,P)=A(P,I)
	A(I,Q)=A(Q,I)
	ENDIF
	ENDDO
	
      GOTO 100  
	END