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
    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)
    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     解三元一次方程
    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

在输入相应数据后显示float point error:invlalid     求高手指点，谢谢