C PROGRAM AKAIKE_IRC C This program was originally described by Prof. H. Akaike of C the Institute of Statistical Mathematics, Tokyo, Japan, and C was slightly modified by T. Morikawa, M. SIGEMORI and T. Wada. C ------------------------------------------------------------------ C 90.7.16 MODIFIED BY SIGEMORI, WADA C 96.4.5 MODIFIED BY MORIKAWA C ------------------------------------------------------------------ PARAMETER (NVAR=3, MAXL=5) C WE ADVISE YOU NOT TO CHANGE THE PARAMETERS WHICH FOLLOW. PARAMETER (LAGH=15,MAXV=5,ML=50,NOP=1000,IPLS=1) C NVAR: NUMBER OF VARIABLES C MAXL: UPPER LIMIT OF MODEL ORDER C LAGH: MAXIMUM LAG FOR COMPUTING AUTO-CORRELATION C MAXV: LIMITATION OF NUMBER OF VARIABLES C ML: LENGTH OF IMPULSE RESPONSE CURVES C NOP: MAXIMUM NUMBER OF POINTS ADDED WITH 1 C IPLS: SIZE OF IMPULSE (0: SIZE=2SD OF INNOVATION; 1:SIZE=1) IMPLICIT REAL*8(A-H,O-Z) CHARACTER*4 TITLE(MAXV) CHARACTER*5 UNIT REAL*4 XX1 REAL*4 X11 REAL*4 WS,X,Y DIMENSION R1(LAGH,MAXV,MAXV),OSD(MAXV,MAXV),SM(MAXV) DIMENSION AO(LAGH,MAXV,MAXV) DIMENSION X11(NOP,MAXV) DIMENSION A1(LAGH,MAXV,MAXV),B1(LAGH,MAXV,MAXV) DIMENSION A(LAGH),C(LAGH) DIMENSION AR(MAXV,MAXV,LAGH) DIMENSION C0(MAXV) DIMENSION SD(MAXV,MAXV),SE(MAXV,MAXV),SF(MAXV,MAXV) DIMENSION XSD(MAXV,MAXV),XSF(MAXV,MAXV) DIMENSION D(MAXV,MAXV),E(MAXV,MAXV),Z1(MAXV,MAXV) DIMENSION SDL(MAXV,MAXV) DIMENSION Z2(MAXV,MAXV) DIMENSION R(MAXV,MAXV),Z(MAXV,MAXV) DIMENSION SD1(MAXV,MAXV) DIMENSION IDS(MAXV) DIMENSION DE(MAXV,MAXV,100),SM1(MAXV,MAXV,100) DIMENSION WS(NOP) DIMENSION X(NOP),Y(NOP) DIMENSION PD(MAXV,MAXV,100),PI(MAXV,MAXV,100) COMMON /COM11/XX1(1000,5) L=MAXL IL=0 C *********INITIAL CONDITION INPUT AND OUTPUT********** CALL ARPREPD2(N,UNIT,NVAR,TITLE,MAXV) K=NVAR+IL DO 303 II=1,K DO 202 JJ=1,N X11(JJ,II)=XX1(JJ,II) 202 CONTINUE 303 CONTINUE IP=NVAR+IL LH=(N-NVAR-1)/(IP+IP) IF(L.GT.LH) L=LH L1=L+1 IP0=K 909 CALL NSMLCO(R1,N,LAGH,K,L,SM,TITLE,MAXV,C0,X,Y,NOP) C FPEC OR MFPE COMPUTATION CALL NFPEC(R1,OSD,AO,L,NVAR,IL,N,IFPEC,LAGH,MAXV,A1,B1, 1 SD,SE,SF,XSD,XSF,D,E,Z1,SDL,Z2,R,Z,SD1,IDS) CALL SDECNV(N,IFPEC,NVAR,IL,OSD,AO,ML,TITLE,SM,UNIT,IPLS, 1 LAGH,A,C,AR,MAXV,DE,SM1,X,Y,NOP,PD,PI) STOP END C sssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss SUBROUTINE ARPREPD2(N,UNIT,NVAR,TITL,MAXV) REAL*4 X C COMMON /COM11/X(1000,5) C CHARACTER*4 TITL(MAXV) CHARACTER*5 UNIT C ?????????????????????? 1:VARIABLE NAMES ?????????????????????????? TITL(1)='A' TITL(2)='B' TITL(3)='C' C TITL(4)='' C TITL(5)='' C ?????????????????????? 2:MODE OF OUTPUTS ????????????????????????? M=NVAR C:TIME UNITS(DAYS, WEEKS etc) UNIT='sec' C:INITIALIZATION N=0 I=1 C ???????????????????? 3:NAME OF DATA FILE ????????????????????????? OPEN(1,STATUS='OLD',FILE='AR.DAT', + ACCESS='SEQUENTIAL',FORM='FORMATTED') OPEN (6,FILE='AIRC98.TXT',STATUS='NEW', + ACCESS='SEQUENTIAL',FORM='FORMATTED') 16 READ(1,17,END=99) (X(I,J),J=1,M) N=N+1 I=I+1 GO TO 16 99 CONTINUE CLOSE(1,STATUS='KEEP') C ????????????????????? 4:FORMAT OF DATA ??????????????????????????? CC 17 FORMAT(9X,3F10.4) CCC 17 FORMAT(4X,F6.4,4X,F6.4,4X,F6.4) 17 FORMAT(3F10.5) WRITE(6,111) DO 18 J=1,M WRITE(6,4) TITL(J) WRITE(6,1) N WRITE(6,50) (X(I,J),I=1,N) 18 CONTINUE RETURN 1 FORMAT(10I5) 4 FORMAT(5X,A4) 50 FORMAT(6F13.5) 111 FORMAT(1H ,20X,'*0*0*0*0* ORIGINAL DATA *0*0*0*0*') END C sssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss SUBROUTINE NSMLCO(R1,N,LAGH,K,L,SM,TITL,MAXV,C0,X,Y,NOP) REAL*8 SM REAL*4 XX1 DOUBLE PRECISION C1,C2,CN1,CN2,C0,CX0,CY0 DOUBLE PRECISION R1 CHARACTER*4 TITL(MAXV) DIMENSION X(NOP),Y(NOP) DIMENSION C1(101),C2(101),CN1(101),CN2(101) DIMENSION SM(MAXV),C0(MAXV) DIMENSION R1(LAGH,MAXV,MAXV) C COMMON /COM11/XX1(1000,5) C L1=L+1 LAGH1=LAGH+1 C MEAN DELETION DO 300 II=1,K DO 310 I=1,N X(I)=XX1(I,II) 310 CONTINUE CALL SMEADL(X,N,XMEAN) DO 320 I=1,N XX1(I,II)=X(I) 320 CONTINUE SM(II)=XMEAN 300 CONTINUE C COVARIANCE COMPUTATION DO 10 II=1,K DO 110 I=1,N X(I)=XX1(I,II) 110 CONTINUE C AUTO COVARIANCE COMPUTATION CALL CROSCO(X,X,N,C1,LAGH1) C NORMALIZATION C0(II)=C1(1) CX0=C0(II) CALL CORNOM(C1,CN1,LAGH1,CX0,CX0) DO 20 I=1,L1 R1(I,II,II)=C1(I) 20 CONTINUE IF(II.EQ.1) GO TO 10 IIM1=II-1 DO 11 JJ=1,IIM1 DO 120 I=1,N Y(I)=XX1(I,JJ) 120 CONTINUE C CROSS COVARIANCE COMPUTATION CALL CROSCO(X,Y,N,C1,LAGH1) CALL CROSCO(Y,X,N,C2,LAGH1) C NORMALIZATION CX0=C0(II) CY0=C0(JJ) CALL CORNOM(C1,CN1,LAGH1,CX0,CY0) CALL CORNOM(C2,CN2,LAGH1,CX0,CY0) C CROSS COVARIANCE PRINT OUT DO 30 I=1,L1 R1(I,II,JJ)=C1(I) R1(I,JJ,II)=C2(I) 30 CONTINUE 11 CONTINUE 10 CONTINUE WRITE (6,601) 601 FORMAT(' ***** ORIGINAL DATA ***** ') 408 DO 303 II=1,K WRITE (6,602) II WRITE (6,603) TITL(II) WRITE (6,604) SM(II),R1(1,II,II) 602 FORMAT(1X,'J = ',I3) 603 FORMAT(1X ,A4) 604 FORMAT(1X ,'MEAN = ',E12.5,5X,'VARIANCE = ',E12.5) 303 CONTINUE 404 CONTINUE RETURN END C sssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss SUBROUTINE CORNOM(C,CN,LAGH1,CX0,CY0) C NORMALIZATION OF COVARIANCE IMPLICIT REAL*8(A-H,O-Z) DIMENSION C(LAGH1),CN(LAGH1) C CST1=1.0D-00 DS=CST1/DSQRT(CX0*CY0) DO 10 I=1,LAGH1 CN(I)=C(I)*DS 10 CONTINUE RETURN END C sssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss SUBROUTINE CROSCO(X,Y,N,C,LAGH1) C THIS SUBROUTINE COMPUTES C(L)=COVARIANCE(X(S+L),Y(S)) C (L=0,1,...,LAGH1-1). DOUBLE PRECISION C,T DOUBLE PRECISION AN,BN1,BN,CT0 DIMENSION X(N),Y(N),C(LAGH1) C AN=N BN1=1.0D-00 BN=BN1/AN CT0=0.0D-00 DO 10 II=1,LAGH1 I=II-1 T=CT0 IL=N-I DO 20 J=1,IL J1=J+I T=T+DBLE(X(J1))*DBLE(Y(J)) 20 CONTINUE C(II)=T*BN 10 CONTINUE RETURN END C sssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss SUBROUTINE SMEADL(X,N,XMEAN) C MEAN DELETION DIMENSION X(N) C AN=N XMEAN=SUMF(X,N)/AN DO 10 I=1,N X(I)=X(I)-XMEAN 10 CONTINUE RETURN END C fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff FUNCTION SUMF(X,N) DIMENSION X(N) C SUMF=0.0 DO 10 I=1,N SUMF=SUMF+X(I) 10 CONTINUE RETURN END C sssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss SUBROUTINE NFPEC(R1,OSD,AO,L,NVAR,IL,N,IFPEC,LAGH,MAXV,A1,B1, 1 SD,SE,SF,XSD,XSF,D,E,Z1,SDL,Z2,R,Z,SD1,IDS) C PROGRAM 5.3.2 FPEC(AR-MODEL FITTING FOR CONTROL) IMPLICIT REAL*8(A-H,O-Z) COMMON /COMA/AIC DIMENSION R1(LAGH,MAXV,MAXV) DIMENSION A1(LAGH,MAXV,MAXV),B1(LAGH,MAXV,MAXV),AO(LAGH,MAXV,MAXV) DIMENSION SD(MAXV,MAXV),SE(MAXV,MAXV),SF(MAXV,MAXV),OSD(MAXV,MAXV) DIMENSION XSD(MAXV,MAXV),XSF(MAXV,MAXV) DIMENSION D(MAXV,MAXV),E(MAXV,MAXV),Z1(MAXV,MAXV) DIMENSION SDL(MAXV,MAXV) C IP=NVAR+IL L1=L+1 C INITIAL CONDITION AND COVARIANCE PRINT OUT WRITE(6,39) WRITE(6,40) WRITE(6,41) N,L,NVAR,IL WRITE(6,42) CALL PRMAT3(R1,L1,IP,IP,1,LAGH,MAXV,MAXV) C INITIAL SD, SF, SE COMPUTATION DO 330 II=1,IP DO 330 JJ=1,IP SD(II,JJ)=R1(1,II,JJ) SDL(II,JJ)=SD(II,JJ) SF(II,JJ)=SD(II,JJ) SE(II,JJ)=R1(2,II,JJ) XSD(II,JJ)=SD(II,JJ) XSF(II,JJ)=SF(II,JJ) 330 CONTINUE C 0-TH STEP COMPUTATION IFPEC=0 MS=0 C OFPEC, ORFPEC COMPUTATION CALL SFPEC(SD,N,IP,NVAR,MS,OFPEC,ORFPEC,OOFPEC,MAXV,SD1) C OFPEC, ORFPEC PRINT OUT WRITE(6,600) WRITE(6,264) MS,OFPEC,ORFPEC,AIC C ITERATION M=1 TO L DO 400 M=1,L C INVERSE OF SD, SF COMPUTATION CALL INVDET(XSD,SDDET,IP,MAXV,IDS) CALL INVDET(XSF,SFDET,IP,MAXV,IDS) C D, E, SD, SF COMPUTATION CALL MULPLY(SE,XSF,D,IP,IP,IP,MAXV,MAXV,MAXV) CALL TRAMDL(SE,XSD,E,IP,IP,IP,MAXV,MAXV,MAXV) CALL TRAMDR(D,SE,Z1,IP,IP,IP,MAXV,MAXV,MAXV) CALL SUBTAL(SD,Z1,IP,IP,MAXV,MAXV) CALL MULPLY(E,SE,Z1,IP,IP,IP,MAXV,MAXV,MAXV) CALL SUBTAL(SF,Z1,IP,IP,MAXV,MAXV) MS=M DO 410 II=1,IP DO 410 JJ=1,IP XSD(II,JJ)=SD(II,JJ) SDL(II,JJ)=SD(II,JJ) XSF(II,JJ)=SF(II,JJ) 410 CONTINUE C FPEC,RFPEC COMPUTATION CALL SFPEC(SD,N,IP,NVAR,MS,FPEC,RFPEC,OOFPEC,MAXV,SD1) C FPEC,RFPEC PRINT OUT WRITE(6,264) MS,FPEC,RFPEC,AIC C FORWARD AND BACKWARD PREDICTOR COMPUTATION CALL COEFAB(A1,B1,D,E,MS,IP,LAGH,MAXV,Z1,Z2) C MIN.FPEC, MIN.RFPEC COMPUTATION IF(M.EQ.1) GO TO 444 IF(OFPEC.LE.FPEC) GO TO 440 444 OFPEC=FPEC ORFPEC=RFPEC IFPEC=M DO 560 II=1,NVAR DO 560 JJ=1,NVAR OSD(II,JJ)=SD(II,JJ) 560 CONTINUE DO 561 I=1,M DO 562 II=1,NVAR DO 562 JJ=1,IP AO(I,II,JJ)=A1(I,II,JJ) 562 CONTINUE 561 CONTINUE 440 IF(M.EQ.L) GO TO 400 C SE COMPUTATION CALL NEWSE(A1,R1,SE,MS,IP,LAGH,MAXV,R,Z) 400 CONTINUE C MIN.FPEC, MIN.RFPEC PRINT OUT WRITE(6,607) OFPEC,ORFPEC,IFPEC C OSD, AO PRINT AND PUNCH OUT WRITE(6,608) CALL SUBMPR(OSD,NVAR,NVAR,MAXV,MAXV) IF(IFPEC.LE.0) GO TO 699 WRITE(6,609) CALL PRMAT3(AO,IFPEC,NVAR,IP,0,LAGH,MAXV,MAXV) 699 RETURN 1 FORMAT(10I5) 2 FORMAT(4D20.10) 39 FORMAT(///1H0,50HPROGRAM 5.3.2 FPEC(AR-MODEL FITTING FOR CONTROL 1)) 40 FORMAT(1X,17HINITIAL CONDITION) 41 FORMAT(1X,2HN=,I5,5X,2HL=,I5,5X,5HNVAR=,I5,5X,3HIL=,I5) 42 FORMAT(//1X,17HCOVARIANCE MATRIX) 264 FORMAT(1X ,I5,2X,3D14.5) 600 FORMAT(///1X,4X,1HI,12X,4HFPEC,9X,5HRFPEC,11X,3HAIC) 607 FORMAT(1X,13HMINIMUM FPEC=,D12.5,2X,14HMINIMUM RFPEC=,D12.5,2X,14 AHATTAINED AT M=,I5) 608 FORMAT(//1X,10X,10HOSD(II,JJ)) 609 FORMAT(//1X,10X,10H(A(I)B(I))) END C sssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss SUBROUTINE COEFAB(A1,B1,D,E,MS,K,LAGH,MAXV,Z1,Z2) C THIS SUBROUTINE COMPUTES FORWARD(A) AND BACKWARD(B) PREDICTOR C COEFFICIENTS. IMPLICIT REAL*8(A-H,O-Z) DIMENSION A1(LAGH,MAXV,MAXV),B1(LAGH,MAXV,MAXV) DIMENSION D(MAXV,MAXV),E(MAXV,MAXV) DIMENSION A(5,5),B(5,5),Z1(MAXV,MAXV),Z2(MAXV,MAXV) C IF(MS.EQ.1) GO TO 40 MSM1=MS-1 DO 10 I=1,MSM1 MMI=MS-I DO 20 II=1,K DO 20 JJ=1,K A(II,JJ)=A1(I,II,JJ) B(II,JJ)=B1(MMI,II,JJ) 20 CONTINUE CALL MULPLY(D,B,Z1,K,K,K,MAXV,MAXV,MAXV) CALL MULPLY(E,A,Z2,K,K,K,MAXV,MAXV,MAXV) CALL SUBTAL(A,Z1,K,K,MAXV,MAXV) CALL SUBTAL(B,Z2,K,K,MAXV,MAXV) DO 21 II=1,K DO 21 JJ=1,K A1(I,II,JJ)=A(II,JJ) B1(MMI,II,JJ)=B(II,JJ) 21 CONTINUE 10 CONTINUE 40 DO 30 II=1,K DO 30 JJ=1,K A1(MS,II,JJ)=D(II,JJ) B1(MS,II,JJ)=E(II,JJ) 30 CONTINUE RETURN END C sssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss SUBROUTINE NEWSE(A1,R1,SE,MS,K,LAGH,MAXV,R,Z) C SE COMPUTATION IMPLICIT REAL*8(A-H,O-Z) DIMENSION A1(LAGH,MAXV,MAXV),R1(LAGH,MAXV,MAXV) DIMENSION SE(MAXV,MAXV) DIMENSION A(5,5),R(MAXV,MAXV),Z(MAXV,MAXV) C CST0=0.0D-00 DO 10 II=1,K DO 10 JJ=1,K Z(II,JJ)=CST0 10 CONTINUE MSP2=MS+2 DO 11 I=1,MS MMI=MSP2-I DO 12 II=1,K DO 12 JJ=1,K A(II,JJ)=A1(I,II,JJ) R(II,JJ)=R1(MMI,II,JJ) 12 CONTINUE CALL MULPLY(A,R,SE,K,K,K,MAXV,MAXV,MAXV) CALL MATADL(Z,SE,K,K,MAXV,MAXV) 11 CONTINUE DO 14 II=1,K DO 14 JJ=1,K R(II,JJ)=R1(MSP2,II,JJ) 14 CONTINUE CALL SUBTAC(R,Z,SE,K,K,MAXV,MAXV) RETURN END C sssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss SUBROUTINE SFPEC(SD,N,K,NVAF,MS,Z,RZ,OOZ,MAXV,SD1) C FPEC COMPUTATION IMPLICIT REAL*8(A-H,O-Z) COMMON /COMA/AIC DIMENSION SD(MAXV,MAXV) DIMENSION SD1(MAXV,MAXV) C AN=N KM=K*MS ANP=N+1+KM ANM=N-1-KM AP=ANP/ANM APR=AP**NVAF CST1=1.0D-00 DO 9 I=1,NVAF DO 9 J=1,NVAF SD1(I,J)=SD(I,J) 9 CONTINUE CALL SUBDET(SD1,SDRM,NVAF,MAXV) Z=APR*SDRM ARM2=2*MS*K*NVAF AIC=AN*DLOG(SDRM)+ARM2 IF(MS.NE.0) GO TO 10 OOZ=CST1/Z 10 RZ=Z*OOZ RETURN END C sssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss SUBROUTINE SUBDET(X,XDETMI,MM,MAXV) IMPLICIT REAL*8(X) DIMENSION X(MAXV,MAXV) C CST0=0.0D-00 CST1=1.0D-00 XDETMI=CST1 IF(MM.EQ.1) GO TO 18 MM1=MM-1 DO 10 I=1,MM1 XDETMI=XDETMI*X(I,I) XC=CST1/X(I,I) I1=I+1 DO 15 J=I1,MM XXC=X(J,I)*XC DO 16 K=I1,MM X(J,K)=X(J,K)-X(I,K)*XXC 16 CONTINUE 15 CONTINUE 10 CONTINUE 18 XDETMI=XDETMI*X(MM,MM) 17 RETURN END C sssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss SUBROUTINE SUBTAC(X,Y,Z,MM,NN,MJ1,MJ2) C MATRIX SUBTRACTION C Z=X-Y IMPLICIT REAL*8(A-H,O-Z) DIMENSION X(MJ1,MJ2),Y(MJ1,MJ2),Z(MJ1,MJ2) C DO 10 I=1,MM DO 10 J=1,NN Z(I,J)=X(I,J)-Y(I,J) 10 CONTINUE RETURN END C sssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss SUBROUTINE SUBTAL(X,Y,MM,NN,MJ1,MJ2) C MATRIX SUBTRACTION C X=X-Y IMPLICIT REAL*8(A-H,O-Z) DIMENSION X(MJ1,MJ2),Y(MJ1,MJ2) C DO 10 I=1,MM DO 10 J=1,NN X(I,J)=X(I,J)-Y(I,J) 10 CONTINUE RETURN END C sssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss SUBROUTINE TRAMDR(X,Y,Z,MM,NN,NC,MJ1,MJ2,MJ3) C TRANSPOSE MULTIPLY (RIGHT) C Z=X*Y' IMPLICIT REAL*8(A-H,O-Z) DIMENSION X(MJ1,MJ2),Y(MJ3,MJ2),Z(MJ1,MJ3) C CST0=0.0D-00 DO 10 I=1,MM DO 11 J=1,NC SUM=CST0 DO 12 K=1,NN SUM=SUM+X(I,K)*Y(J,K) 12 CONTINUE Z(I,J)=SUM 11 CONTINUE 10 CONTINUE RETURN END C sssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss SUBROUTINE INVDET(X,XDET,MM,MAXV,IDS) IMPLICIT REAL*8(X) DIMENSION X(MAXV,MAXV) DIMENSION IDS(MAXV) C CST0=0.0D-00 CST1=1.0D-00 XDET=CST1 DO 10 L=1,MM C PIVOTING AT L-TH STAGE XMAXP=0.10000D-10 MAXI=0 DO 110 I=L,MM C FOR COMPLEX VERSION NEXT STATEMENT SHOULD BE REPLACED BY C IF(CDABS(XMAXP).GE.CDABS(X(I,L))) GO TO 110 IF(DABS(XMAXP).GE.DABS(X(I,L))) GO TO 110 XMAXP=X(I,L) MAXI=I 110 CONTINUE IDS(L)=MAXI IF(MAXI.EQ.L) GO TO 120 IF(MAXI.GT.0) GO TO 121 XDET=CST0 GO TO 140 C ROW INTERCHANGE 121 DO 14 J=1,MM XC=X(MAXI,J) X(MAXI,J)=X(L,J) X(L,J)=XC 14 CONTINUE XDET=-XDET 120 XDET=XDET*XMAXP XC=CST1/XMAXP X(L,L)=CST1 DO 11 J=1,MM X(L,J)=X(L,J)*XC 11 CONTINUE DO 12 I=1,MM IF(I.EQ.L) GO TO 12 XC=X(I,L) X(I,L)=CST0 DO 13 J=1,MM X(I,J)=X(I,J)-XC*X(L,J) 13 CONTINUE 12 CONTINUE 10 CONTINUE IF(MM.GT.1) GO TO 123 GO TO 140 C COLUMN INTERCHANGE 123 MM1=MM-1 DO 130 J=1,MM1 MMJ=MM-J JJ=IDS(MMJ) IF(JJ.EQ.MMJ) GO TO 130 DO 131 I=1,MM XC=X(I,JJ) X(I,JJ)=X(I,MMJ) X(I,MMJ)=XC 131 CONTINUE 130 CONTINUE 140 RETURN END C sssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss SUBROUTINE MATADL(X,Y,MM,NN,MJ1,MJ2) C MATRIX ADDITION C X=X+Y IMPLICIT REAL*8(A-H,O-Z) DIMENSION X(MJ1,MJ2),Y(MJ1,MJ2) C DO 10 I=1,MM DO 10 J=1,NN X(I,J)=X(I,J)+Y(I,J) 10 CONTINUE RETURN END C sssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss SUBROUTINE MULPLY(X,Y,Z,MM,NN,NC,MJ1,MJ2,MJ3) C MATRIX MULTIPLICATION C Z=X*Y IMPLICIT REAL*8(A-H,O-Z) DIMENSION X(MJ1,MJ2),Y(MJ2,MJ3),Z(MJ1,MJ3) C CST0=0.0D-00 DO 10 I=1,MM DO 11 J=1,NC SUM=CST0 DO 12 K=1,NN SUM=SUM+X(I,K)*Y(K,J) 12 CONTINUE Z(I,J)=SUM 11 CONTINUE 10 CONTINUE RETURN END C sssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss SUBROUTINE TRAMDL(X,Y,Z,MM,NN,NC,MJ1,MJ2,MJ3) C TRANSPOSE MULTIPLY (LEFT) C Z=X'*Y IMPLICIT REAL*8(A-H,O-Z) DIMENSION X(MJ1,MJ2),Y(MJ1,MJ3),Z(MJ2,MJ3) C CST0=0.0D-00 DO 10 I=1,NN DO 11 J=1,NC SUM=CST0 DO 12 K=1,MM SUM=SUM+X(K,I)*Y(K,J) 12 CONTINUE Z(I,J)=SUM 11 CONTINUE 10 CONTINUE RETURN END C sssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss SUBROUTINE PRMAT3(X1,MM,NN,NC,ISHIFT,LAGH,MJ1,MJ2) C 3-DIMENSIONAL MATRIX X1(MM,NN,NC) PRINT OUT C (LAGH,MJ1,MJ2): ABSOLUTE DIMENSION OF X1 IN THE MAIN ROUTINE C NN,NC: SHOULD BE LESS THAN 11 IMPLICIT REAL*8(A-H,O-Z) DIMENSION X1(LAGH,MJ1,MJ2) DIMENSION X(5,5) C WRITE(6,60) MM DO 61 I=1,MM DO 62 II=1,NN DO 62 JJ=1,NC X(II,JJ)=X1(I,II,JJ) 62 CONTINUE I1=I-ISHIFT WRITE(6,63) I1 CALL SUBMPR(X,NN,NC,MJ1,MJ2) 61 CONTINUE RETURN 60 FORMAT(1X,3HLL=,I5) 63 FORMAT(/1X,2HI=,I5) END C sssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss SUBROUTINE SUBMPR(RM,MR,LC,MJ1,MJ2) C THIS SUBROUTINE PRINTS OUT UPPER LEFT MR X LC OF RM. C (MJ1,MJ2): ABSOLUTE DIMENSION OF X IN THE MAIN ROUTINE IMPLICIT REAL*8(A-H,O-Z) DIMENSION RM(MJ1,MJ2) C WRITE(6,39) MR,LC WRITE(6,20) DO 10 I=1,MR J1=0 J2=0 JC=0 J1=J2+1 16 CONTINUE J2=J1+9 IF(J2.LE.LC) GO TO 13 J2=LC 13 JC=JC+1 IF(JC.GT.1) GO TO 19 18 WRITE(6,21) I,(RM(I,J),J=J1,J2) GO TO 17 19 WRITE(6,22) (RM(I,J),J=J1,J2) 17 IF(J2.LT.LC) GO TO 16 10 CONTINUE RETURN 39 FORMAT(1X,6HMATRIX,5X,I5,4X,1HX,I5) 20 FORMAT(16X,1H1,11X,1H2,11X,1H3,11X,1H4,11X,1H5,11X,1H6,11X,1 AH7,11X,1H8,11X,1H9,10X,2H10) 21 FORMAT(1X,I5,4X,10D12.5) 22 FORMAT(1X ,9X,10D12.5) END C sssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss SUBROUTINE SDECNV(N,M0,NVAF,IL,OSD,A1,ML,TITL,SM,UNIT,IPLS,LAGH, 1 A,C,AR,MAXV,DE,SM1,X,Y,NOP,PD,PI) C PROGRAM 5.3.4 DECONVOLUTION IMPLICIT REAL*8(A-H,O-Z) REAL*4 X,Y DIMENSION A1(LAGH,MAXV,MAXV),A(LAGH),C(LAGH),SA(101),C1(101) DIMENSION OSD(MAXV,MAXV),SM(MAXV),P3(100) DIMENSION DE(MAXV,MAXV,100),SM1(MAXV,MAXV,100),P1(100),P2(100) DIMENSION AR(MAXV,MAXV,LAGH) DIMENSION PD(MAXV,MAXV,100),PI(MAXV,MAXV,100) CHARACTER*4 TITL(MAXV) CHARACTER*5 UNIT C C ABSOLUTE DIMENSIONS USED FOR SUBROUTINE CALL CST0=0.0D-00 IP=NVAF+IL C INITIAL CONDITION PRINT OUT WRITE(6,60) WRITE(6,61) WRITE(6,62) ML WRITE(6,63) N,M0,NVAF,IL WRITE(6,64) CALL PRMAT3(A1,M0,NVAF,IP,0,LAGH,MAXV,MAXV) C COMPUTATION START ML=ML+1 MLP1=ML+1 DO 10 II=1,IP C C ARRANGEMENT DO 11 L=1,M0 C(L)=A1(L,II,II) 11 CONTINUE DO 20 JJ=1,IP IF(II.EQ.JJ) GO TO 20 C A ARRANGEMENT DO 22 L=1,M0 A(L)=A1(L,II,JJ) 22 CONTINUE C INITIAL CONDITIONING DO 30 I=1,MLP1 C1(I)=CST0 SA(I)=CST0 30 CONTINUE C SA COMPUTATION SA(2)=A(1) DO 40 I=2,ML IM1=I-1 IF(IM1.GT.M0) GO TO 42 C1(I)=C(IM1) 42 IP1=I+1 CALL CONVOL(C1,SA,SA2,IP1) IF(I.LE.M0) SA2=A(I)+SA2 SA(IP1)=SA2 DE(II,JJ,1)=SA(1) DE(II,JJ,2)=SA(2) DE(II,JJ,IP1)=SA(IP1) 40 CONTINUE 20 CONTINUE 10 CONTINUE CALL SMCON1(M0,NVAF,IL,A1,ML,SM1,LAGH,AR,MAXV,X,Y,NOP) COS2=2.0 DO 111 II=1,IP DO 222 JJ=1,IP WN0=SQRT(OSD(JJ,JJ)) WN=COS2*SQRT(OSD(JJ,JJ)) IF(IPLS.NE.0) WN=1.0*IPLS DO 333 LL=1,ML P2(LL)=WN*SM1(II,JJ,LL) PI(II,JJ,LL)=WN*SM1(II,JJ,LL) 333 CONTINUE DO 444 LL=1,ML P1(LL)=WN*DE(II,JJ,LL) PD(II,JJ,LL)=WN*DE(II,JJ,LL) 444 CONTINUE CALL PRCOL3(P1,P2,1,ML,-4,TITL,II,JJ,WN,SM,UNIT,OSD,WN0, 1 IPLS,MAXV) 222 CONTINUE 111 CONTINUE WRITE(6,1111) DO 3001 I=1,IP WRITE(6,661) I WRITE(6,663) '1','2','3','4','5' DO 3002 L=1,ML IF(J.EQ.I) GO TO 3002 WRITE(6,662) (PD(J,I,L),J=1,IP) 3002 CONTINUE 3001 CONTINUE WRITE(6,1212) DO 3003 I=1,IP WRITE(6,661) I WRITE(6,663) '1','2','3','4','5' DO 3004 L=1,ML WRITE(6,662) (PI(J,I,L),J=1,IP) 3004 CONTINUE 3003 CONTINUE RETURN 1 FORMAT(10I5) 60 FORMAT(29HPROGRAM 5.3.4 DECONVOLUTION) 61 FORMAT(1X,17HINITIAL CONDITION) 62 FORMAT(1X,3HML=,I5) 63 FORMAT(1X ,2HN=,I5,5X,3HM0=,I5,5X,5HNVAR=,I5,5X,3HIL=,I5) 64 FORMAT(1X,7HA1(I,J)) 160 FORMAT(2X,2HI=,I5,5X,2HJ=,I5) 1111 FORMAT(5X,'*** IRC DIRECT VALUE ***') 1212 FORMAT(5X,'*** IRC INDIRECT VALUE ***') 661 FORMAT(7H INPUT=,I2) 662 FORMAT(5F10.5) 663 FORMAT(1X,7HOUTPUT=,A1,4(9X,A1)) END C sssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss SUBROUTINE CONVOL(A,B,SUM,K1) C THIS SUBROUTINE COMPUTES CONVOLUTION. C C(K)=A(0)B(K)+A(1)B(K-1)+...+A(K-1)B(1)+A(K)B(0) C K1: K PLUS 1 IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(K1),B(K1) C CST0=0.0D-00 K2=K1+1 SUM=CST0 DO 10 I=1,K1 KI=K2-I SUM=SUM+A(I)*B(KI) 10 CONTINUE RETURN END C sssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss SUBROUTINE SMCON1(K,NVAR,IL,AO,ML,SM1,LAGH,AR,MAXV,X,Y,NOP) IMPLICIT REAL*8(A-H,O-Z) REAL*4 X,Y DIMENSION AR(MAXV,MAXV,LAGH),AO(LAGH,MAXV,MAXV) DIMENSION X(NOP),Y(NOP),YY(100) DIMENSION SM1(MAXV,MAXV,100) C SCST1=1.0 CST0=0.0D-00 DO 100 II=1,MAXV DO 100 JJ=1,MAXV DO 100 KK=1,LAGH AR(II,JJ,KK)=CST0 100 CONTINUE NSH=0 N1=ML+NSH ID=NVAR+IL C ..... SIGN CHANGE ........ DO 150 II=1,K DO 160 I=1,ID DO 160 J=1,ID AR(I,J,II)=-AO(II,I,J) 160 CONTINUE 150 CONTINUE C -------- IMPLS=1 3210 CONTINUE JJ=IMPLS KK=IMPLS CALL SUBIMP(X,N1,ID,IMPLS,NOP) CALL MLARMA(AR,Y,X,K,N1,ID,MAXV,LAGH,NOP) N1ID=N1*ID NID=ML*ID DO 590 IX=1,ID I=0 DO 595 IC=IX,NID,ID I=I+1 YY(I)=Y(IC) SM1(IX,JJ,I)=YY(I) 595 CONTINUE 590 CONTINUE IF(IMPLS.GE.ID) GO TO 4500 IMPLS=IMPLS+1 GO TO 3210 4500 CONTINUE RETURN END C sssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss SUBROUTINE SUBIMP(WS,N1,ID,IMPLS,NOP) REAL*4 WS,SCST0 DIMENSION WS(NOP) C N1ID=N1*ID IF(IMPLS.LE.0) RETURN SCST0=0.0 DO 340 II=1,N1ID WS(II)=SCST0 340 CONTINUE WS(IMPLS)=1.0 RETURN END C sssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss SUBROUTINE MLARMA(AR,Y,X,K,N,ID,MJ1,LAGH,NOP) REAL*8 AR,Z,CST0 REAL*4 X,Y DIMENSION AR(MJ1,MJ1,LAGH) DIMENSION X(NOP),Y(NOP) C SCST0=0.0 CST0=0.0 NID=N*ID DO 10 I=1,NID Y(I)=SCST0 10 CONTINUE IMD=0 DO 100 I=1,N C ----- AR(1)Y(I-1)+...+AR(K)Y(I-K) COMPUTATION ----- DO 200 J=1,K IJ=I-J IMG=(IJ-1)*ID IF(IJ.LE.0) GO TO 200 DO 210 IX=1,ID Z=CST0 DO 220 I2=1,ID IG=IMG+I2 Z=Z+AR(IX,I2,J)*Y(IG) 220 CONTINUE IC=IMD+IX Y(IC)=Y(IC)-Z 210 CONTINUE 200 CONTINUE DO 420 IX=1,ID IC=IMD+IX Y(IC)=Y(IC)+X(IC) 420 CONTINUE IMD=IMD+ID 100 CONTINUE RETURN END C sssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss SUBROUTINE PRCOL3(P1,P2,INDI,INDL,ISHIFT,TITL,II,JJ,WN,SM,UNIT, 1 SD,WN0,IPLS,MAXV) C THIS SUBROUTINE PRINTS VECTOR (P1(I),I=INDI,INDL) C COLUMNWISE IN THE FORMAT (I-ISHIFT,P1(I),I=INDI,INDL). IMPLICIT REAL*8(A-H,O-Z) CHARACTER Y*1,UNIT*5,Z*1 CHARACTER*4 TITL(20) DIMENSION P1(100),P2(100),Y(41),SM(MAXV) DIMENSION Z(41),SD(MAXV,MAXV) C S5=0.5 SD1=SQRT(SD(II,II)) SDD=S5*SQRT(SD(II,II)) WN1=IPLS*(1.0/WN0) CNT40=0.4D+02 C SEARCH MAX(ABS(P1)) P1MAX=0.0D-00 DO 1000 I=INDI,INDL P1ABS=ABS(P1(I)) IF(P1MAX.LT.P1ABS) P1MAX=P1ABS 1000 CONTINUE DO 1001 I=INDI,INDL P2ABS=ABS(P2(I)) IF(P1MAX.LT.P2ABS) P1MAX=P2ABS IF(P1MAX.LT.SDD) P1MAX=SDD P2MAX=P1MAX 1001 CONTINUE IF(IPLS.NE.0) GO TO 1011 WRITE(6,5500) II,TITL(II),SM(II),JJ,TITL(JJ),WN 5500 FORMAT(/5X,7HOUTPUT=,I3,3H ; ,A3,2X,7H( MEAN=,F10.4,2H ), 1 /5X,6HINPUT=,I4,3H ; ,A3,2X,8H( IMPLS:,F9.4,7H =2SD )) GO TO 1022 1011 CONTINUE WRITE(6,5505) II,TITL(II),SM(II),JJ,TITL(JJ),WN,WN1 5505 FORMAT(/5X,7HOUTPUT=,I3,3H ; ,A3,2X,7H( MEAN=,F10.4,2H ),10X, 1 /5X,6HINPUT=,I4,3H ; ,A3,2X,8H( IMPLS:,F6.3,2H =,F8.4,4HSD )) 1022 CONTINUE C FIRST LINE WRITE P1MAXM=-P1MAX 5011 WRITE(6,5033)SDD,P1MAXM,P1MAX C WRITE(6,5010) WRITE(6,5010) 5033 FORMAT(47X,10H( # 0.5SD=,F9.4,1H)/ 1 2X,1HN,7X,4HD(N),8X,4HI(N),F9.4,15X,1H0,11X,F10.4) P0A=SDD+P1MAX P0A=P0A/P1MAX*CNT40 P0IA=AINT(P0A) IP0A=INT(P0A) IP0A=(IP0A+3)/2 P00A=-SDD+P1MAX P00A=P00A/P1MAX*CNT40 P00IA=AINT(P00A) IP00A=INT(P00A) IP00A=(IP00A+3)/2 C WRITE GRAPH DO 2000 I=INDI,INDL IM1=I-ISHIFT MODIM1=MOD(IM1,3) IF(MODIM1.EQ.0) GO TO 2010 C SET ' ' DO 2020 J=1,41 Y(J)=' ' Z(J)=' ' 2020 CONTINUE GO TO 2030 C SET '~~~~~~~~~' 2010 DO 2040 J=1,41 Y(J)='-' Z(J)='-' 2040 CONTINUE C SET '|' 2030 DO 2050 J=1,41,10 Y(J)='|' Z(J)='|' 2050 CONTINUE C SET '#' Y(IP0A)='#' Y(IP00A)='#' C SET 'D','I','+' P1A=P1(I)+P1MAX P1A=P1A/P1MAX*CNT40 P1IA=AINT(P1A) IP1A=INT(P1A) IP1A=(IP1A+3)/2 Y(IP1A)='D' P2A=P2(I)+P1MAX P2A=P2A/P1MAX*CNT40 P2IA=AINT(P2A) IP2A=INT(P2A) IP2A=(IP2A+3)/2 IF(IP2A.EQ.IP1A) THEN Y(IP2A)='+' ELSE Y(IP2A)='I' ENDIF C SET 'D','I','+','*' I1=I-1 WRITE(6,5020)I1,P1(I),P2(I),(Y(J),J=1,41) GO TO 2000 2000 CONTINUE WRITE(6,5030) UNIT RETURN 5010 FORMAT(30X,'|---------|---------|---------|---------|') 5020 FORMAT(1X ,I3,2F12.5,2X,41A1) 5060 FORMAT(1X ,I3,2D12.4,2X,41A1,3X,41A1) 5030 FORMAT(3X,4H( X ,A5,2H )) END