text1.f90

本程序采用fortran编写

	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