魏凤英老师的ACTA的fortran程序为什么一运行总是显示错误

青衣

魏凤英老师的ACTA的fortran程序为什么一运行总是显示错误，之后我把第18行的OPEN(4,FILE='GT.RES',STATUS='NEW')，改成了OPEN(4,FILE='e:\aa.txt',STATUS='NEW'),结果显示没有错误了，aa.TXT生成了一级差分和二级差分序列，以及CSC序列，但其他的结果就显示运行出错没有结果，不知道怎么回事，先把程序和数据粘出来希望能得到帮助。
        PROGRAM ACTA
        PARAMETER (LY=100,LZ=100,LX=50)
        DIMENSION X(100),X1(100),X2(100),X3(100),X4(100),
     &        XX1(100),XX2(100),XX3(100)
        COMMON/AA/XX(3,LZ)
        COMMON/BB/PMGF(LZ,LY)
        COMMON/CC/XY(LZ,LX)
        COMMON/DD/A(LX,LX)
        COMMON/EE/EX(LX,LX)
        COMMON/FF/COE(LX,LX)
        COMMON/GG/FX(LZ,LX)
        WRITE(*,10)
  10        FORMAT(5X,'N=?,NN=?,IN=?')
        READ(*,*)N,NN,IN
        N1=N+NN
        OPEN(11,FILE='GT.dat')
        READ(11,*)(X(I),I=1,N)
        OPEN(4,FILE='GT.RES',STATUS='NEW')
        CALL ACTA (X,X1,N,NN,N1,IN)
        W=0.0
        DO 20 I=1,N
  20        W=W+(X(I)-X1(I))*(X(I)-X1(I))
        V=SQRT(W/N)
        WRITE(4,30)V
  30        FORMAT(30X,'THE FINAL RESULT OF FITTING AND FORECAST'/
     &        (5X,'RMSE=',F15.2))
        DO 40 I=1,N1
        WRITE(4,50)I,X(I),X1(I)
  50        FORMAT(1X,I5,F10.2,F10.2)
  40        CONTINUE
        STOP
        END
C--------------------------------------END-----------------------------------
        SUBROUTINE ACTA(X,X1,N,NN,N1,IN)
C        THIS IS A SUBROUTINE FOR TIME SERIES ANALYSIS AND FORECAST
C        BASED ON MEAN GENERATING FUNCTION OF THE SERIES AND
C        DIFFERENCE AND ADDITION
        integer    ly, lz, n, nn, in, i, j, mm, m, num, ip,
     &        n1, l, k, kk, kl
        real   fx, x, x1, x2, f, xxx, csca, sy, b, pmgf,xy,
     &        x0, kf, xx,a,ex,coe
        parameter ( ly = 100, lz = 100,lx=50)

        DIMENSION X(lz),X1(lz),X2(lz),F(lz),NUM(lz),
     &            IP(lz),XXX(lz),CSCA(lz),SY(lz),B(lz)
        COMMON/AA/XX(3,lz)
        COMMON/BB/PMGF(lz, ly)
        COMMON/CC/XY(lz, lx)
        COMMON/DD/A(lx,lx)
        COMMON/EE/EX(lx,lx)
        COMMON/FF/COE(lx,lx)
        COMMON/GG/FX(lz,lx)
111        format(5x,'mark1')
        DO 30 I=1,N-1
  30           X1(I)=X(I+1)-X(I)
        X1(I)=0.0
        DO 40 I=1,N-2
  40           X2(I)=X1(I+1)-X1(I)
        X2(N-1)=0.0
        X2(N)=0.0
        WRITE(4,50)(X1(I),I=1,N)
  50        FORMAT(30X,'DIFFERENCE DATA-1:'/(5X,10F7.2))
        WRITE(4,60)(X2(I),I=1,N)
  60        FORMAT(30X,'DIFFERENCE DATA-2:'/(5X,10F7.2))
        write(*,222)
222        format(5x,'mark2')
        DO 70 I=1,N
           XX(1,I)=X(I)
           XX(2,I)=X1(I)
  70             XX(3,I)=X2(I)
        N1=N+NN
        MM=INT(N/3)
        M=MM*3
        CALL MMGF(N,3,MM,N1)
        M=M-3
        DO 71 I=MM,2*MM-2
        DO 71 J=1,N1
  71           FX(J,I-MM+1)=PMGF(J,I)
        DO 72 I=1,MM-1
           PMGF(1,I+M)=X(1)
           X0=X(1)
           DO 73 J=2,N1
              X0=X0+FX(J-1,I)
  73              PMGF(J,I+M)=X0
  72           CONTINUE
        M=M+MM-1
C        DO 74 I=1,M
C        WRITE(4,75)(PMGF(J,I),J=1,N)
C 75        FORMAT(30X,'MGF'/(5X,10F7.2))
C 74        CONTINUE
        DO 80 I=1,M
           DO 90 J=1,N
  90              F(J)=PMGF(J,I)
           IF(I.LT.MM)THEN
              NUM(I)=0
             IP(I)=I+1
C        WRITE(4,100)NUM(I),IP(I)
100        FORMAT(//30X,'MGF--NO.',I3,'-',I3)
           ELSE IF(I.GE.MM.AND.I.LT.2*MM-1)THEN
              NUM(I)=1
              IP(I)=I+1-MM+1
C        WRITE(4,100)NUM(I),IP(I)
           ELSE IF(I.GE.2*MM.AND.I.LT.3*MM-1)THEN
              NUM(I)=2
              IP(I)=I+2-2*MM+1
C        WRITE(4,100)NUM(I),IP(I)
           ELSE
              NUM(I)=3
              IP(I)=I+3-3*MM+1
C        WRITE(4,100)NUM(I),IP(I)
           END IF
           CALL RA1(N,X,F,XXX)
           CALL CSC(N,X,XXX,CSC3,1)
           CSCA(I)=CSC3
  80           CONTINUE
        write(*,333)
333        format(5x,'mark3')
        KF=13.39
        WRITE(4,110)
110        FORMAT(//10X,'NO.',10X,'PERIOD',10X,'CSC'/)
        L=0
        DO 120 I=1,M
           IF(CSCA(I).GT.KF)THEN
              L=L+1
              WRITE(4,130)L,NUM(I),IP(I),CSCA(I)
130        FORMAT(/10X,I3,5X,I4,'-',I4,10X,F10.2)
              CSCA(L)=CSCA(I)
              DO 140 J=1,N1
140                 pmgf(J,L)=PMGF(J,I)
           END IF
120           CONTINUE

        DO 141 I=1,L
141       IP(I)=I

        DO 142 I=1,L-1
           DO 143 J=I+1,L
              IF(CSCA(I).GE.CSCA(J))GOTO 143
              W=CSCA(I)
              CSCA(I)=CSCA(J)
              CSCA(J)=W
              K=IP(I)
              IP(I)=IP(J)
              IP(J)=K
143              CONTINUE
142           CONTINUE

        WRITE(4,144)(IP(I),I=1,L)
144        FORMAT(//30X,'MAX-MIN-N0:'/(5X,10I5))
        WRITE(4,145)(CSCA(I),I=1,L)
145        FORMAT(//30X,'MAX-MIN-CSC:'/(5X,10F7.2))
        IF(L.LT.IN)IN=L
        DO 146 I=1,IN
           KL=IP(I)
           DO 147 J=1,N1
147              FX(J,I)=pmgf(J,KL)
146           CONTINUE
        DO 148 I=1,IN
        DO 148 J=1,N
148         XY(J,I)=FX(J,I)
        DO 149 I=1,N
149        XY(I,IN+1)=X(I)
        IIN=IN+1
        CALL SREG(X1,N,N1,IN,IIN)
        RETURN
        END
C----------------------------------END---------------------------------------

        SUBROUTINE MMGF(N,M,MM,N1)
        integer  ly, lz, i, j, l, n, kk, m, mm, n1, k, kkk
        real sb, xx, pmgf
        parameter (ly = 100, lz = 100)        
        DIMENSION SB(lz)
        COMMON/AA/XX(3,lz)
        COMMON/BB/PMGF(lz,ly)
        DO 10 K=1,M
           DO 20 KK=2,MM
              DO 30 J=1,KK
   30                 SB(J)=0.0
              DO 40 J=1,KK
                 I=0
                 DO 50 L=J,N,KK
                    SB(J)=SB(J)+XX(K,L)
                    I=I+1
   50                    CONTINUE
                 SB(J)=SB(J)/I
   40                 CONTINUE
              DO 60 L=1,N1
                 J=L-L/KK*KK
                 IF(J.EQ.0)J=KK
                 KKK=(K-1)*(MM-1)+KK-1
                 PMGF(L,KKK)=SB(J)
   60                 CONTINUE
   20              CONTINUE
   10           CONTINUE

C        DO 70 I=1,M*MM-M
C        WRITE(4,80)I,(PMGF(J,I),J=1,N1)
C   80        FORMAT(30X,'MGF'/5X,I5/(3X,10F7.1))
C   70        CONTINUE
        RETURN
        END
C-------------------------------------END------------------------------------

C        SUBROUTINE RA1 FOR REGRESSION OF ONE VARIABLE
C        N--SAMPLE NUMBER
C        X--TIME SERIES
C        F--MEAN GENERATING FUNCTION SERIES
C        X1--FITTING SERIES OF X
        SUBROUTINE RA1(N,X,F,X1)
        integer lz, i, n
        real x, x1, c, d, fc, xc, fxc,df, dx, dfx, b, b0, r, f
        parameter ( lz = 100)
        DIMENSION X(lz),F(lz),X1(lz)
        C=0.0
        D=0.0
        FC=0.0
        XC=0.0
        FXC=0.0
        DO 10 I=1,N
           C=C+F(I)
           D=D+X(I)
           FC=FC+F(I)*F(I)
           XC=XC+X(I)*X(I)
  10           FXC=FXC+F(I)*X(I)
        DF=FC-C*C/N
        DX=XC-D*D/N
        DFX=FXC-D*C/N
C        WRITE(4,50)C,D,FC,XC,FXC,DF,DX,DFX
  50        FORMAT(2X,8F9.2)
        B=DFX/DF
        B0=(D-B*C)/N
        R=DFX/SQRT(DF*DX)
        DO 20 I=1,N
  20           X1(I)=B0+B*F(I)
C        WRITE(4,30)B0,B
  30        FORMAT(//5X,'B0=',F15.7,5X,'B1=',F15.7/)
C        WRITE(4,40)(X1(I),I=1,N)
  40        FORMAT(//30X,'FITTING OF X'/(5X,10F7.1)/)
        RETURN
        END
C--------------------------------------END-----------------------------------
C        SUBROUTINE CALCULATING CSC
C        N--SAMPLE NUMBER
C        X--TIME SERIES
C        X1--FITTING SERIES OF X
C        K--FACTOR NUMBER
        SUBROUTINE CSC(N,X,X1,CSC3,K)
        integer lz, n, k, i, j, le
        real x, x1, p, q, xp, xp1, u, v, r1, r2, r3,
     &        s1, s2, xm1, qk, qx, csc3
        parameter (lz = 100)
        DIMENSION X(lz),X1(lz),LE(3,3),P(3),Q(3)
        U=0.0
        V=0.0
        DO 10 I=1,N-1
           U=U+ABS(X1(I+1)-X1(I))/(N-1)
  10           V=V+ABS(X(I+1)-X(I))/(N-1)
        DO 20 I=1,3
        DO 20 J=1,3
  20           LE(I,J)=0.0
        DO 30 I=1,N-1
           XP=X(I+1)-X(I)
           XP1=X1(I+1)-X1(I)
        IF(XP1.GE.U.AND.XP.GE.V)LE(1,1)=LE(1,1)+1
        IF(XP1.GE.U.AND.XP.LE.0)LE(1,3)=LE(1,3)+1
        IF(XP1.GE.U.AND.XP.LT.V.AND.XP.GT.0)LE(1,2)=LE(1,2)+1
        IF(XP1.LT.U.AND.XP1.GT.0.AND.XP.GE.V)LE(2,1)=LE(2,1)+1
        IF(XP1.LT.U.AND.XP1.GT.0.AND.XP.LT.V.AND.XP.GT.0)LE(2,2)=LE(2,2)+1
        IF(XP1.LT.U.AND.XP1.GT.0.AND.XP.LE.0)LE(2,3)=LE(2,3)+1
        IF(XP1.LE.0.AND.XP.GE.V)LE(3,1)=LE(3,1)+1
        IF(XP1.LE.0.AND.XP.GT.0.AND.XP.LT.V)LE(3,2)=LE(3,2)+1
        IF(XP1.LE.0.AND.XP.LE.0)LE(3,3)=LE(3,3)+1
  30           CONTINUE
C        WRITE(4,40)((LE(I,J),J=1,3),I=1,3)
c  40        FORMAT(//30X,'CONTINGENCY TABLE'/(5X,3I4/))

        R1=0.0
        R2=0.0
        R3=0.0
        DO 50 I=1,3
           P(I)=0.0
  50           Q(I)=0.0
        DO 55 I=1,3
        DO 55 J=1,3       
           IF(LE(I,J).NE.0.0)R1=R1+LE(I,J)*LOG(FLOAT(LE(I,J)))
  55           CONTINUE
        DO 60 I=1,3
        DO 60 J=1,3
           P(I)=P(I)+LE(I,J)
  60           Q(I)=Q(I)+LE(J,I)
        DO 70 I=1,3
           IF(Q(I).NE.0.0)R2=R2+Q(I)*LOG(Q(I))
           IF(P(I).NE.0.0)R3=R3+P(I)*LOG(P(I))
  70           CONTINUE
        S1=2*(R1+(N-1)*LOG(FLOAT(N-1))-(R2+R3))
        XM1=0.0
        DO 80 I=1,N
  80           XM1=XM1+X(I)/FLOAT(N)
        QK=0.0
        QX=0.0
        DO 90 I=1,N
           QK=QK+(X(I)-X1(I))*(X(I)-X1(I))/N
  90           QX=QX+(X(I)-XM1)*(X(I)-XM1)/N
        S2=(N-K)*(1-QK/QX)
        CSC3=S1+S2
C        WRITE(4,100)S1,S2,CSC3
C 100        FORMAT(//5X,'S1=',F10.2,5X,'S2=',F10.2,5X,'CSC=',F10.2//)
        RETURN
        END
C------------------------------------END--------------------------------------
        subroutine sreg(x1,n,n1,ip,iip)
        integer n,n1,ip,iip,lx,ly,lz
        parameter(lx=50,lz=100)
          dimension xym(lx),ks(2000),ma(2000),cc(lx,lx),xx(lz,lx),xf(lx,lz)
        dimension ihh(lx),rss(lx),x1(lz),x2(lz),x(lz),csca(lx)
        common/ee/ex(lx,lx)
        common/ff/coe(lx,lx)
        common/cc/xy(lz,lx)
        common/dd/a(lx,lx)
        common/gg/fx(lz,lx)
        integer ts
        do 10 j=1,ip+1
        s=0
        do 20 i=1,n
   20   s=s+xy(i,j)
        xym(j)=s/n
   10   continue
        write(4,11)(xym(i),i=1,iip)
   11        format('mean=',10f10.2)
        do 30 i=1,ip+1
        do 30 j=1,ip+1
        s=0
        do 40 k=1,n
   40   s=s+xy(k,i)*xy(k,j)
        s=s-n*xym(i)*xym(j)
        a(i,j)=s
        a(j,i)=s
   30   continue
        sto=a(ip+1,ip+1)
        ih=2**(ip-1)-1
        ks(1)=1
        do 50 k=2,ip-1
        j=2**(k-1)
        ks(j)=k
        do 55 i=1,j-1
        m=ks(j-i)
        ks(j+i)=-m
   55   continue
   50        continue
        nb=0
        do 60 i=1,ip
        rss(i)=10.0**20
        ma(i)=0
  60    continue
        ts=0
        do 70 m=1,ih
        it=abs(ks(m))
C        write(4,71)it
C  71        format(2x,'IT=',I3)         
        iip=ip+1
        if(ks(m).gt.0) then
        call reg350(K,ts,ma,nb,rss,iip,it,ip)
        else
        call reg430(K,ts,ma,nb,rss,iip,it,ip)
        end if
  70    continue
        it=ip
        call reg350(K,ts,ma,nb,rss,iip,it,ip)
        do 80 m=1,ih
        m1=(-1)*ks(ih+1-m)
        it=abs(m1)
        if(m1.gt.0) then
         call reg350(K,ts,ma,nb,rss,iip,it,ip)
        else
        call reg430(K,ts,ma,nb,rss,iip,it,ip)
        end if
   80   continue
        do 90 k=1,ip
        write(4,91)k
   91        format(//30x,i2,'ORDER REGRESSION SET ')
        do 92 j=1,ip
        if( ex(k,j).eq.0) goto 92
        write(4,93)j
   93        format(5x,'X',i2,2x)
   92        continue
        s=rss(k)
        c=0.0
        do 100 j=1,ip
        if(ex(k,j).eq.0) goto 100
        c=c+coe(k,j)*xym(j)
  100   continue
        c=xym(ip+1)-c
        write(4,101)c
  101        format(/30x,'REGRESSION EQUATION'/2X,'Y=',f10.3)
        do 102 j=1,ip
        if(ex(k,j).eq.0)goto 102
        write(4,104)coe(k,j),j
  104        format('+',f10.3,'x',i2)
  102        continue       
        write(4,105) sto
        write(4,106) sto-s
        write(4,107) s
        write(4,108) sqrt((sto-s)/sto)
  105   format(/5x,'DEPARTURE SQUARE SUM =',f10.4)
  106   format(5X,'REGRESSION SQUARE SUM=',f10.4)
  107        format(5x,'RESIDUES SQUARE SUM=',f10.4)
  108        format(5x,'MULTIPLE CORRELATION COEFFICIENT=',F10.4)
        mm=0.0
        do 301 j=1,ip
        if(ex(k,j).eq.0.0)goto 301
        mm=mm+1
        cc(k,mm)=coe(k,j)
        do 300 i=1,n1
  300        xx(i,mm)=fx(i,j)
  301        continue
        do 302 i=1,n1
        xf(k,i)=c
        do 303 j=1,mm
  303        xf(k,i)=xf(k,i)+cc(k,j)*xx(i,j)
        x(i)=xy(i,iip)
        x1(i)=xf(k,i)
  302        continue
        do 401 j=1,n
  401        x2(j)=xf(k,j)-x(j)
        w=0.0
        do 402 j=1,n
  402        w=w+x2(j)*x2(j)
        v=sqrt(w/n)
        call csc(n,x,x1,csc3,k)
        csca(k)=csc3
        write(4,308)xym(iip),v
  308        format(5x,'MEAN--Y=',f10.2,2x,'RMSE=',f10.2)
        write(4,304)(x1(i),i=1,n)
  304        format(/30x,'FITTING'/(5x,10f8.2))
        write(4,305)(x1(i),i=n+1,n1)
  305        format(/30x,'FORECAST'/(5x,10f8.2))
   90   continue
        aa=0.0
        do 201 i=1,ip
        if(csca(i).gt.aa)aa=csca(i)
  201        continue
        do 202 i=1,ip
        if(csca(i).eq.aa)goto203
  202        continue
  203        ihh(1)=i
        WRITE(4,181)ihh(1)
  181        format(//5x,'REGRESSION SET IS--'I3/)
        KK=ihh(1)
        do 180 i=1,ip
        if(ex(KK,i).eq.0) goto 180
        write(4,182)i
  180        continue
  182        format(5x,'X',i2,2x)
        do 188 j=1,ip
        if(j.EQ.KK) THEN
        DO 185 I=1,N1
        X1(I)=XF(J,I)
  185        CONTINUE
        ELSE
        GOTO 188
        END IF
  188        continue
        return
        end
C---------------------------------END----------------------------------------
        subroutine reg350(K,ts,ma,nb,rss,iip,it,ip)
        PARAMETER(lx=50)
        integer ts
        dimension ma(1000),rss(ip)
        common /ee/ex(lx,lx)
        common/ff/coe(lx,lx)
        common/dd/a(lx,lx)
        k=it
        call reg300(k,iip)
        ts=ts+1
C       write(4,15)ts
C  15   format(10x,i5)
        ma(it)=1
        nb=nb+1
        if(a(iip,iip).gt.rss(nb)) goto 111
        rss(nb)=a(iip,iip)
        do 10 j=1,ip
        ex(nb,j)=ma(j)
        coe(nb,j)=a(j,iip)
   10   continue
  111        return
        end
C----------------------------------END---------------------------------------

        subroutine reg430(K,ts,ma,nb,rss,iip,it,ip)
        integer ts
        parameter(lx=50)
        dimension ma(1000),rss(ip)
        common /ee/ex(lx,lx)
        common/ff/coe(lx,lx)
        common/dd/a(lx,lx)
        k=it
        call reg300(k,iip)
        ts=ts+1
C        write(*,15) ts
C  15   format(10x,i5)
        ma(it)=0
        nb=nb-1
        if(a(iip,iip).gt.rss(nb)) goto 111
        rss(nb)=a(iip,iip)
        do 10 j=1,ip
        ex(nb,j)=ma(j)
        coe(nb,j)=a(j,iip)
   10   continue
  111        return
        end
C-------------------------------------END------------------------------------
        subroutine reg300(k,iip)
        parameter(lx=500)
        common/dd/a(lx,lx)
        if(abs(a(k,k)).gt.10**(-12)) goto 200
        write(*,*) 'error stop'
        goto 999
  200   do 10 i=1,iip
        do 10 j=1,iip
        if(i.eq.k.or.j.eq.k) goto 10
        a(i,j)=a(i,j)-a(i,k)*a(k,j)/a(k,k)
   10   continue
        do 20 j=1,iip
        if(j.eq.k) goto 20
        a(k,j)=a(k,j)/a(k,k)
        a(j,k)=-a(j,k)/a(k,k)
   20   continue
        a(k,k)=1/a(k,k)
        return
  999        end