[Fortran] 纯文本查看 复制代码
SUBROUTINE DOOLITTLE(A,L,U,N)
!Doolittle分解
IMPLICIT REAL*8(A-Z)
INTEGER::N,I,K,R
REAL*8::A(N,N),L(N,N),U(N,N)
N = 501
!U的第一行
U(1,:)=A(1,:)
!L的第一列
L(:,1)=A(:,1)/U(1,1)
DO K=2,N
L(K,K)=1
DO J=K,N
S=0
DO M=1,K-1
S=S+L(K,M)*U(M,J)
END DO
U(K,J)=A(K,J)-S
END DO
DO I=K+1,N
S=0
DO M=1,K-1
S=S+L(I,M)*U(M,K)
END DO
L(I,K)=(A(I,K)-S)/U(K,K)
END DO
END DO
END SUBROUTINE DOOLITTLE
SUBROUTINE UPTRI(A,B,X,N)
!解上三角方程组
IMPLICIT REAL*8(A-Z)
INTEGER::I,J,K,N
REAL*8::A(N,N),B(N),X(N)
N = 501
X(N)=B(N)/A(N,N)
!回代部分
DO I=N-1,1,-1
X(I)=B(I)
DO J=I+1,N
X(I)=X(I)-A(I,J)*X(J)
END DO
X(I)=X(I)/A(I,I)
END DO
END SUBROUTINE UPTRI
SUBROUTINE DOWNTRI(A,B,X,N)
!解下三角方程组
IMPLICIT REAL*8(A-Z)
INTEGER::I,J,N
REAL*8::A(N,N),B(N),X(N)
N = 501
X(1)=B(1)/A(1,1)
DO K=2,N
X(K)=B(K)
DO I=1,K-1
X(K)=X(K)-A(K,I)*X(I)
END DO
X(K)=X(K)/A(K,K)
END DO
END SUBROUTINE DOWNTRI
SUBROUTINE UNMIFA(A,NAMDA4)
IMPLICIT NONE
INTEGER::N,I,K,L,S,J,M,K1
DIMENSION A(501,501),U(501),U0(501),V(501),U1(501,501),L1(501,501),Y(501)
DOUBLE PRECISION MA,A,NAMDA4,U,U0,V,M0,M1,TOL,U1,L1,Y
N = 501
DO I=1,N
U0(N)=10D0
END DO
!给出u0
U=U0
M0=0
!------------------------------------------------
CALL DOOLITTLE(A,L1,U1,N)
!------------------------------------------------
DO K1=1,500
CALL DOWNTRI(L1,U,Y,N)
CALL UPTRI(U1,Y,V,N)
!------------------------------------------------
!用解方程组的方法求矩阵v
!求最大值
M1=DABS(V(1))
DO I=2,N
IF (DABS(V(I))>M1) THEN
M1=DABS(V(I))
L = I !用k记录指标，但是ma不取r(i)
END IF
END DO
M1 = V(L)
!找出v中的最大值，返回到m1 SUBROUTINE MAX_ROU(R,N,MA)
U = V/M1
!IF (DABS(M1-M0)<1E-20) EXIT
IF (DABS(DABS(M1)-DABS(M0))/DABS(M1)<TOL) EXIT
M0=M1
END DO
NAMDA4=M1
NAMDA4=1/NAMDA4
END
PROGRAM MAIN
IMPLICIT NONE
INTEGER I,J,N
DIMENSION A(501,501)
DOUBLE PRECISION A,NAMDA4
!给a赋值!!!!!!!!!!!!!!!!!!!!!!!!!
DO I = 1,501
DO J = 1,501
IF(I == J) THEN
A(I,J) =(1.64-0.024*I)*SIN(0.2*I)-0.64*EXP(0.1/I)
ELSE IF(I == J+1) THEN
A(I,J) = 0.16
ELSE IF(I+1 == J) THEN
A(I,J) = 0.16
ELSE IF(I == J+2) THEN
A(I,J) = -0.064
ELSE IF(I+2 == J) THEN
A(I,J) = -0.064
!ELSE
! A(I,J) = 0
END IF
END DO
END DO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
CALL UNMIFA(A,NAMDA4)
WRITE(*,*) NAMDA4
END