主题:float point error 新帖
iamyangel@126
[专家分:0] 发布于 2010-04-06 20:13:00
PROGRAM MAIN
REAL,DIMENSION(1:6)::XR,XI,A
REAL,DIMENSION(1:3)::CL,CM,CN,SP,CK,AA,BB,CC,CHECK
REAL::EPS
READ(*,*) SX,SY,SZ,SXY,SYZ,SXZ
C 求解空间任一点的主应力和主应力方向余弦
AI1=SX+SY+SZ
AI2=SX*SY+SY*SZ+SX*SZ-SXY**2-SYZ**2-SXZ**2
AI3=SX*SY*SZ+2*SXY*SYZ*SXZ-SX*SYZ**2-SY*SXZ**2-SZ*SXY**2
A(1)=-AI1
A(2)=AI2
A(3)=-AI3
EPS=1.0D-10
C 其中的AI1.AI2,AI3是中间量,这步应该没什么问题
CALL DQRRT(A,3,XR,XI,EPS,B,L)
IF (L.NE.0) THEN
SP(1)=XR(1)
SP(2)=XR(2)
SP(3)=XR(3)
ENDIF
CALL HTOL(SP,3)
DO 12 J=1,3
CHECK(J)=1.0
AA(J)=(SY-SP(J))*(SZ-SP(J))-SYZ**2
BB(J)=-(SXY*(SZ-SP(J))-SXZ*SYZ)
CC(J)=SXY*SYZ-(SXZ*(SY-SP(J)))
CK(J)=1.0/SQRT(AA(J)**2+BB(J)**2+CC(J)**2)
CL(J)=AA(J)*CK(J)
CM(J)=BB(J)*CK(J)
CN(J)=CC(J)*CK(J)
CL,CM,CN 是计算出的方向余弦,即是最后结果
CHECK(J)=CL(J)**2+CM(J)**2+CN(J)**2
SSMAX=(SP(1)-SP(3))/2.0
WRITE(*,315) J,SP(J),CL(J),CM(J),CN(J)
315 FORMAT(1X,I3,5X,D15.8,3(5X,D15.8))
12 CONTINUE
END
C 主应力大小求解
SUBROUTINE DQRRT(A,N,XR,XI,EPS,B,L)
DIMENSION B(N,N),XR(N),XI(N),A(N)
DOUBLE PRECISION B,XR,XI,A,EPS
DO 10 J=1,N
10 B(1,J)=-A(J)
DO 20 I=2,N
DO 20 J=1,N
20 B(I,J)=0.0
DO 30 I=1,N-1
30 B(I+1,I)=1.0
CALL CHHQR(B,N,XR,XI,EPS,L)
RETURN
END
C 解三元一次方程
C [color=FF0000]这个子程序是我从一本书上直接抄过来的,所以里面很多代码变量我也不清楚,麻烦各位C 高人前辈帮我找找里面是否有问题[/color]
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,EPS
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
C 由大到小排序
SUBROUTINE HTOL(A,N)
REAL*8 A(N)
DO 20 I=1,N-1
DO 30 J=1,N-I
IF(A(J).LT.A(J+1)) THEN
C=A(J)
A(J)=A(J+1)
A(J+1)=C
ENDIF
30 CONTINUE
20 CONTINUE
END
在输入[color=FF00FF]六个应力分量即SX=.1500E+03,SY=-.1000E+03,SZ=.8000E+02 SXY=-.5000E+02,SYZ=.4000E+02,SXZ=.3000E+02 [/color]
后显示float point error:invlalid 求高手指点,谢谢
P.S:谢谢上个帖子里雪球的指点,我是才接触这种论坛,多多指教,谢谢大家的帮忙
[em11]
最后更新于:2010-04-07 09:54:00
回复列表 (共9个回复)
沙发
asymptotic [专家分:16630] 发布于 2010-04-06 23:58:00
养成好的编程习惯很重要,为什么不用 implicit none 关掉所有隐式声明呢?
CALL DQRRT(A,3,XR,XI,EPS,B,L) 编译器告诉你,这一句有问题;因为 caller 是单精度,而 callee 是双精度。
既然是菜鸟,那就一步一步做吧,先把所有变量都申明了。
板凳
臭石头雪球 [专家分:23030] 发布于 2010-04-07 08:15:00
还是那个数组越界问题,自己好好查吧
3 楼
臭石头雪球 [专家分:23030] 发布于 2010-04-07 08:30:00
1 0.16574303D+03 0.95315790D+00 -0.13889344D+00 0.26869822D+00
2 0.84085464D+02 -0.22566131D+00 0.26499623D+00 0.93747211D+00
3 -0.11982851D+03 0.20141263D+00 0.95419371D+00 -0.22124034D+00
请按任意键继续. . .
这是我的答案
4 楼
臭石头雪球 [专家分:23030] 发布于 2010-04-07 08:41:00
正如渐近线前辈所言,老程序最大的问题就是不定义变量,大量的变量使用了,却没有定义。
这就造成实参,虚参,变量名字,变量类型,变量精度非常混淆。
强烈建议所有新手,老手,熟手,僵尸手全部使用 Implicit None!!!强制定义所有变量。
实参是 REAL ,那么对应的虚参也必须是 REAL。如果是 REAL*8 那么也要统一!!
5 楼
iamyangel@126 [专家分:0] 发布于 2010-04-07 09:52:00
嗯,雪球前辈的答案和我那本参考书上的一样[em28][em28],请问您是怎么得到的?我还是没法克服那个问题,是不是就是您说的那个“数组越界问题”啊,其实您说的这个越界问题我还是不很清楚,请您再说详细些,好么?
非常感谢各位前辈的指点,小女子受益很大[em28][em28][em28]
6 楼
iamyangel@126 [专家分:0] 发布于 2010-04-07 09:58:00
谢谢前辈的指导,呵呵,我会记住您说的
7 楼
臭石头雪球 [专家分:23030] 发布于 2010-04-07 13:09:00
如果我的结果没错的话,你可以这样改:
PROGRAM MAIN
Implicit None
REAL*8,DIMENSION(1:6)::XR,XI,A,B[color=FF0000] !//部分变量要变成 Real*8 的,以和子程序相匹配[/color]
REAL,DIMENSION(1:3)::CL,CM,CN,SP,CK,AA,BB,CC,CHECK
REAL*8::EPS [color=FF0000] !//这里也是[/color]
Real SX,SY,SZ,SXY,SYZ,SXZ,AI1,AI2,AI3,ssmax
Integer L,j
!READ(*,*) SX,SY,SZ,SXY,SYZ,SXZ
SX=.1500E+03
SY=-.1000E+03
SZ=.8000E+02
SXY=-.5000E+02
SYZ=.4000E+02
SXZ=.3000E+02
C 求解空间任一点的主应力和主应力方向余弦
AI1=SX+SY+SZ
AI2=SX*SY+SY*SZ+SX*SZ-SXY**2-SYZ**2-SXZ**2
AI3=SX*SY*SZ+2*SXY*SYZ*SXZ-SX*SYZ**2-SY*SXZ**2-SZ*SXY**2
A(1)=-AI1
A(2)=AI2
A(3)=-AI3
EPS=1.0D-10
C 其中的AI1.AI2,AI3是中间量,这步应该没什么问题
CALL DQRRT(A,3,XR,XI,EPS,B,L)
IF (L.NE.0) THEN
SP(1)=XR(1)
SP(2)=XR(2)
SP(3)=XR(3)
ENDIF
CALL HTOL(SP,3)
DO 12 J=1,3
CHECK(J)=1.0
AA(J)=(SY-SP(J))*(SZ-SP(J))-SYZ**2
BB(J)=-(SXY*(SZ-SP(J))-SXZ*SYZ)
CC(J)=SXY*SYZ-(SXZ*(SY-SP(J)))
CK(J)=1.0/SQRT(AA(J)**2+BB(J)**2+CC(J)**2)
CL(J)=AA(J)*CK(J)
CM(J)=BB(J)*CK(J)
CN(J)=CC(J)*CK(J)
CL,CM,CN 是计算出的方向余弦,即是最后结果
CHECK(J)=CL(J)**2+CM(J)**2+CN(J)**2
SSMAX=(SP(1)-SP(3))/2.0
WRITE(*,315) J,SP(J),CL(J),CM(J),CN(J)
315 FORMAT(1X,I3,5X,D15.8,3(5X,D15.8))
12 CONTINUE
END
8 楼
臭石头雪球 [专家分:23030] 发布于 2010-04-07 13:09:00
C 主应力大小求解
SUBROUTINE DQRRT(A,N,XR,XI,EPS,B,L)
Implicit None
Integer L,i,j,n [color=FF0000]!//强制定义变量,并定义[/color]
DIMENSION B(N,N),XR(N),XI(N),A(N)
DOUBLE PRECISION B,XR,XI,A,EPS
DO 10 J=1,N
10 B(1,J)=-A(J)
DO 20 I=2,N
DO 20 J=1,N
20 B(I,J)=0.0
DO 30 I=1,N-1
30 B(I+1,I)=1.0
CALL CHHQR(B,N,XR,XI,EPS,L)
RETURN
END
C 解三元一次方程
C 这个子程序是我从一本书上直接抄过来的,所以里面很多代码变量我也不清楚,麻烦各位C 高人前辈帮我找找里面是否有问题
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,EPS
JT=1
M=N
IT=0
10 IF (M.EQ.0) RETURN
L=M+1
40 L=L-1
IF (L.gt.1) then [color=FF0000]!//这里改一下,防止越界[/color]
If (ABS(A(L,L-1)).GT.EPS*(ABS(A(L-1,L-1))+ABS(A(L,L)))) then
goto 40
End If
EndIf
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
C 由大到小排序
SUBROUTINE HTOL(A,N)
Implicit None [color=FF0000]!//这里加上它,强制变量定义[/color]
REAL A(N),c
Integer i,j,n [color=FF0000] !//定义变量[/color]
DO 20 I=1,N-1
DO 30 J=1,N-I
IF(A(J).LT.A(J+1)) THEN
C=A(J)
A(J)=A(J+1)
A(J+1)=C
ENDIF
30 CONTINUE
20 CONTINUE
END
9 楼
iamyangel@126 [专家分:0] 发布于 2010-04-07 20:43:00
谢谢雪球前辈,辛苦了,小女子真是感激涕零[em28]
我来回复