徐士良老师书中程序，遇到问题
SUBROUTINE CHHBG(A,N)DIMENSION A(N,N)
DOUBLE PRECISION A,D,T
DO 100 K=2,N-1
D=0.0
DO 10 J=K,N
IF (ABS(A(J,K-1)).GT.ABS(D)) THEN
D=A(J,K-1)
I=J
END IF
10 CONTINUE
IF (ABS(D)+1.0.NE.1.0) THEN
IF (I.NE.K) THEN
DO 20 J=K-1,N
T=A(I,J)
A(I,J)=A(K,J)
20 A(K,J)=T
DO 30 J=1,N
T=A(J,I)
A(J,I)=A(J,K)
30 A(J,K)=T
END IF
DO 90 I=K+1,N
T=A(I,K-1)/D
A(I,K-1)=0.0
DO 40 J=K,N
40 A(I,J)=A(I,J)-T*A(K,J)
DO 50 J=1,N
50 A(J,K)=A(J,K)+T*A(J,I)
90 CONTINUE
END IF
100 CONTINUE
RETURN
END
SUBROUTINE CHHQR(A,N,U,V,EPS,JT)
DIMENSION A(N,N),U(N),V(N)
DOUBLE PRECISION A,U,V,W,B,C,X,Y,XY,P,Q,R,E,F,Z,G
JT=1
M=N
IT=0
10 IF (M.EQ.0) RETURN
L=M+1
40 L=L-1
IF (ABS(A(L,L-1)).GT.EPS*(ABS(A(L-1,L-1))+ABS(A(L,L))).AND.
* L.GT.1) GOTO 40
IF (L.EQ.M) THEN
U(M)=A(M,M)
V(M)=0.0
M=M-1
IT=0
GOTO 10
END IF
IF (L.EQ.M-1) THEN
B=-(A(M,M)+A(M-1,M-1))
C=A(M,M)*A(M-1,M-1)-A(M,M-1)*A(M-1,M)
W=B*B-4*C
Y=SQRT(ABS(W))
IF (W.GT.0.0) THEN
XY=1.0
IF (B.LT.0.0) XY=-1.0
U(M)=(-B-XY*Y)/2.0
U(M-1)=C/U(M)
V(M)=0.0
V(M-1)=0.0
ELSE
U(M)=-B/2.0
U(M-1)=U(M)
V(M)=Y/2.0
V(M-1)=-V(M)
END IF
M=M-2
IT=0
GOTO 10
END IF
IF (IT.GE.60) THEN
WRITE(*,50)
50 FORMAT(1X,'FAIL')
JT=0
RETURN
END IF
IT=IT+1
DO 60 J=L+2,M
60 A(J,J-2)=0.0
DO 70 J=L+3,M
70 A(J,J-3)=0.0
DO 150 K=L,M-1
IF (K.NE.L) THEN
P=A(K,K-1)
Q=A(K+1,K-1)
R=0.0
IF (K.NE.M-1) R=A(K+2,K-1)
ELSE
X=A(M,M)+A(M-1,M-1)
Y=A(M-1,M-1)*A(M,M)-A(M-1,M)*A(M,M-1)
P=A(L,L)*(A(L,L)-X)+A(L,L+1)*A(L+1,L)+Y
Q=A(L+1,L)*(A(L,L)+A(L+1,L+1)-X)
R=A(L+1,L)*A(L+2,L+1)
END IF
IF (ABS(P)+ABS(Q)+ABS(R).NE.0.0) THEN
XY=1.0
IF (P.LT.0.0) XY=-1.0
S=XY*SQRT(P*P+Q*Q+R*R)
IF (K.NE.L) A(K,K-1)=-S
E=-Q/S
F=-R/S
X=-P/S
Y=-X-F*R/(P+S)
G=E*R/(P+S)
Z=-X-E*Q/(P+S)
DO 110 J=K,M
P=X*A(K,J)+E*A(K+1,J)
Q=E*A(K,J)+Y*A(K+1,J)
R=F*A(K,J)+G*A(K+1,J)
IF (K.NE.M-1) THEN
P=P+F*A(K+2,J)
Q=Q+G*A(K+2,J)
R=R+Z*A(K+2,J)
A(K+2,J)=R
END IF
A(K+1,J)=Q
A(K,J)=P
110 CONTINUE
J=K+3
IF (J.GE.M) J=M
DO 120 I=L,J
P=X*A(I,K)+E*A(I,K+1)
Q=E*A(I,K)+Y*A(I,K+1)
R=F*A(I,K)+G*A(I,K+1)
IF (K.NE.M-1) THEN
P=P+F*A(I,K+2)
Q=Q+G*A(I,K+2)
R=R+Z*A(I,K+2)
A(I,K+2)=R
END IF
A(I,K+1)=Q
A(I,K)=P
120 CONTINUE
END IF
150 CONTINUE
GOTO 10
END
DIMENSION A(5,5),U(5),V(5)
DOUBLE PRECISION A,U,V
DATA A/1.0,8.0,-2.0,-13.0,17.0,6.0,-15.0,11.0,2.0,
* 22.0,-3.0,18.0,9.0,21.0,-5.0,-1.0,5.0,15.0,
* 30.0,3.0,7.0,4.0,20.0,-6.0,6.0/
CALL CHHBG(A,5)
WRITE(*,10) ((A(I,J),J=1,5),I=1,5)
10 FORMAT(1X,5D15.6)
EPS=1.0E-06
CALL CHHQR(A,5,U,V,EPS,JT)
IF (JT.NE.0) THEN
DO 20 I=1,5
WRITE(*,30) U(I),V(I)
20 CONTINUE
END IF
30 FORMAT(1X,D15.6,2X,'+J',1X,D15.6)
END
什么问题？ necrohan 发表于 2020-12-23 16:09
什么问题？
您好，程序会报错，是数组越界的问题
第45，46 行
IF (ABS(A(L,L-1)).GT.EPS*(ABS(A(L-1,L-1))+ABS(A(L,L))).AND.
* L.GT.1) GOTO 40
替换为
if( l > 1 ) then
IF (ABS(A(L,L-1))>EPS*(ABS(A(L-1,L-1))+ABS(A(L,L)))) GOTO 40
end if 布衣龙共 发表于 2020-12-23 22:05
第45，46 行
IF (ABS(A(L,L-1)).GT.EPS*(ABS(A(L-1,L-1))+ABS(A(L,L))).AND.
* L.GT.1)...
谢谢你
页:
[1]