szjs.f90

带状矩阵特征值jocobi算法的fortran源码

!***************************************************************************!
!                                                                           !
!  PROGRAM: SZJS                                                            !
!                                                                           !
!  PURPOSE: 计算求解数值问题                                                !
!                                                                           !
!***************************************************************************!

	PROGRAM SZJS

	INTEGER P,Q
	
	DOUBLE PRECISION C(1:5,1:501)
	
	DOUBLE PRECISION MAX
	
	CALL INIT(C)

	DO 150, WHILE(MAX(C).GT.1E-12)
		CALL MAXPQ(P,Q,C)
		CALL CAL(P,Q,C)
		PRINT *,MAX(C)	
150 CONTINUE

	OPEN(UNIT=1,FILE='RESULTS.TXT',FORM='FORMATTED')

	WRITE(1,*)C(1,1),C(1,2),C(1,3),C(1,4),C(1,499),C(1,500),C(1,501),'$$$$$'
	WRITE(1,*)C(2,1),C(2,2),C(2,3),C(2,4),C(2,499),C(2,500),C(2,501),'$$$$$'
	WRITE(1,*)C(3,1),C(3,2),C(3,3),C(3,4),C(3,499),C(3,500),C(3,501),'$$$$$'
	WRITE(1,*)C(4,1),C(4,2),C(4,3),C(4,4),C(4,499),C(4,500),C(4,501),'$$$$$'
	WRITE(1,*)C(5,1),C(5,2),C(5,3),C(5,4),C(5,499),C(5,500),C(5,501),'$$$$$'
	WRITE(1,*)MAX(C)

	CLOSE(1)

	END PROGRAM SZJS

!************************初始化矩阵*****************************************!

	SUBROUTINE INIT(A)

	DOUBLE PRECISION A(1:5,1:501)

	DO 10 I=1,499
 		A(1,I+2)=-0.064
10	CONTINUE

	DO 20 I=1,500
		A(2,I+1)=0.16
20	CONTINUE

	DO 30 I=1,501
		A(3,I)=(1.64-0.024*I)*SIN(0.2*I)-0.64*EXP(0.1/I)
30	CONTINUE

	DO 40 I=1,500
		A(4,I)=0.16
40	CONTINUE

	DO 50 I=1,499
		A(5,I)=-0.064
50	CONTINUE

	END

!***************************************************************************!


!****************寻找矩阵A非主对角线元素中按模最大的元素********************!

	SUBROUTINE MAXPQ(X,Y,A)
	INTEGER X,Y
	DOUBLE PRECISION A(1:5,1:501)

	DOUBLE PRECISION S
	S=0.0
	DO 60 I=1,499
	
		IF((ABS(A(1,I+2))).GT.S) THEN
			S=ABS(A(1,I+2))
			X=I
			Y=I+2
		END IF
60 CONTINUE

	DO 70 I=1,500
		IF((ABS(A(2,I+1))).GT.S) THEN
			S=ABS(A(2,I+1))
			X=I
			Y=I+1
		END IF
70	CONTINUE

		
	END

!***************************************************************************!

!********************计算转换矩阵A1*****************************************!

	SUBROUTINE CAL(X,Y,A)
	INTEGER X,Y
	DOUBLE PRECISION A(1:5,1:501)

	DOUBLE PRECISION ALF
	DOUBLE PRECISION B(1:5,1:501)

	ALF=ATAN(2*A(X-Y+3,Y)/(A(3,X)-A(3,Y)))/2

	DO 80 I=1,5
		DO 90 J=1,501
			B(I,J)=A(I,J)
90		CONTINUE
80	CONTINUE

	DO 100 I=1,5
		B(I,X)=A(I,X)*COS(ALF)+A(I,Y)*SIN(ALF)
		B(I,Y)=-A(I,X)*SIN(ALF)+A(I,Y)*COS(ALF)
		IF(((I+X-3).GE.1).AND.((I+X-3).LE.501)) THEN
			B((6-I),(I+X-3))=B(I,X)
		END IF
		IF(((I+Y-3).GE.1).AND.((I+Y-3).LE.501)) THEN
			B((6-I),(I+Y-3))=B(I,Y)
		END IF
100	CONTINUE

	B(3,X)=A(3,X)*COS(ALF)**2+A(3,Y)*SIN(ALF)**2+2*A(X-Y+3,Y)*COS(ALF)*SIN(ALF)
	B(3,Y)=A(3,X)*SIN(ALF)**2+A(3,Y)*COS(ALF)**2-2*A(X-Y+3,Y)*COS(ALF)*SIN(ALF)

	B(X-Y+3,Y)=0.5*(A(3,Y)-A(3,X))*SIN(2*ALF)+A(X-Y+3,Y)*COS(2*ALF)
	B(Y-X+3,X)=0.5*(A(3,Y)-A(3,X))*SIN(2*ALF)+A(X-Y+3,Y)*COS(2*ALF)

	DO 110 I=1,5
		DO 120 J=1,501
			A(I,J)=B(I,J)
120		CONTINUE
110	CONTINUE

	END
!***************************************************************************!

!****************计算经过变换后矩阵中的按模取元素的最大值*******************!

	FUNCTION MAX(A)
	DOUBLE PRECISION MAX
	DOUBLE PRECISION A(1:5,1:501)

	DO 130 I=1,499
		IF(ABS(A(1,I+2)).GT.MAX) THEN
			MAX=ABS(A(1,I+2))
		END IF
130	CONTINUE

	DO 140 I=1,500
		IF(ABS(A(2,I+1)).GT.MAX) THEN
			MAX=ABS(A(2,I+1))
		END IF
140	CONTINUE

	RETURN

	END
!***************************************************************************!