您所在位置：主题：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]
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]
