|
 |
Software Release
Disclaimer:
// EXEC FORT1CLG,REGION=512K 00000010
//* 00000020
/*JOBPARM COPIES=4 00000030
C 00000040
C PROGRAM GFLS: 00000050
C EQUATION NUMBERS IN COMMENT STATEMENTS REFER TO 00000060
C ''FLEXIBLE LEAST SQUARES FOR APPROXIMATELY LINEAR SYSTEMS'' 00000070
C BY R. KALABA AND L. TESFATSION, IEEE TRANSACTIONS ON SYSTEMS, 00000080
C MAN, AND CYBERNETICS 20 (1990), PAGES 978-989. 00000090
C LAST UPDATED: 14 JUNE 1994 00000100
C 00000110
IMPLICIT REAL*8(A-H,O-Z) 00000120
INTEGER T,TCAP,TCAP1 00000130
REAL*8 M 00000140
C 00000150
C THIS PROGRAM IS CURRENTLY DIMENSIONED FOR A MAXIMUM OF TCAP=110 00000160
C OBSERVATION VECTORS Y OF MAXIMUM DIMENSION MOBS = 15 WITH STATE 00000170
C VECTORS X OF MAXIMUM DIMENSION N = 15. 00000180
C 00000190
DIMENSION QO(15,15),PO(15,15),RO(15,15),QZERO(15,15) 00000200
DIMENSION PZERO(15,15),RZERO(15,15) 00000210
DIMENSION F(15,15),A(15,15),H(15,15),B(15,15),D(15,15) 00000220
DIMENSION M(15,15),Y(15,15),TRUEX(15,110),YY(15,110) 00000230
DIMENSION HT(15,15),U(15,15),C(15,15),W(15,15),V(15,15) 00000240
DIMENSION E(15,15),Z(15,15),G(15,15),QNEW(15,15),PNEW(15,15) 00000250
DIMENSION S(15,15),RNEW(15,15),XTCAP(15,15),X(15,110) 00000260
DIMENSION AA(15,15),BB(15,15),CC(15,15),DD(15,15),EE(15,15) 00000270
DIMENSION FF(15,15),HH(15,15),OO(15,15),PP(15,15),QQ(15,15) 00000280
DIMENSION RR(15,15),TT(15,15) 00000290
C 00000300
C ADDITIONAL ARRAYS IF SMOOTHED ESTIMATES ARE TO BE CALCULATED 00000310
C FOR INTERMEDIATE X VALUES (I.E., IF IFLAGS IS SET AT 1) 00000320
C 00000330
DIMENSION GG(15,15,110),SS(15,110) 00000340
C 00000350
C THE FOLLOWING SUBROUTINE INITIALIZES THE PENALTY WEIGHT AMU, 00000360
C THE NUMBER OF OBSERVATIONS TCAP, THE STATE VECTOR DIMENSION 00000370
C N, THE OBSERVATION VECTOR DIMENSION MOBS, AND THE INITIAL COST 00000380
C FUNCTION CHARACTERISTICS QZERO, PZERO, AND RZERO. IT ALSO SETS 00000390
C THE VALUE FOR A FLAG "IFLAGR" TO DETERMINE IF RNEW IS TO BE 00000400
C CALCULATED (1) OR NOT (0) AND A FLAG "IFLAGS" TO DETERMINE IF 00000410
C SMOOTHED ESTIMATES FOR INTERMEDIATE X VALUES ARE TO BE CALCULATED 00000420
C (1) OR NOT (0). 00000430
C 00000440
CALL INPUT(AMU,TCAP,N,MOBS,QZERO,PZERO,RZERO,IFLAGR,IFLAGS) 00000450
CALL SHIFT(N,N,QZERO,QO) 00000460
CALL SHIFT(N,1,PZERO,PO) 00000470
CALL SHIFT(1,1,RZERO,RO) 00000480
C 00000490
C ENTERING THE MAIN DO LOOP FOR GENERATING Q,P,AND R FOR 00000500
C SUCCESSIVE TIMES T = 1,TCAP USING EQS.(24), (26), AND (27). 00000510
C 00000520
DO 50 T=1,TCAP 00000530
CALL MODEL(T,F,A,H,B,D,M,Y,TRUEX) 00000540
DO 5 I=1,MOBS 00000550
YY(I,T) = Y(I,1) 00000560
5 CONTINUE 00000570
C 00000580
C GETTING U=HT*M*H + QO IN EQ.(28) 00000590
C 00000600
CALL MUL(MOBS,MOBS,N,M,H,AA) 00000610
CALL TRANS(MOBS,N,H,HT) 00000620
CALL MUL(N,MOBS,N,HT,AA,BB) 00000630
CALL ADD(N,N,BB,QO,U) 00000640
C 00000650
C GETTING C=FT*D 00000660
C 00000670
CALL TRANS(N,N,F,AA) 00000680
CALL MUL(N,N,N,AA,D,C) 00000690
C 00000700
C GETTING W=AMU*C*F+U 00000710
C 00000720
CALL MUL(N,N,N,C,F,AA) 00000730
CALL MULCON(N,N,AMU,AA,BB) 00000740
CALL ADD(N,N,BB,U,W) 00000750
C 00000760
C GETTING V=WINV IN EQ.(17) 00000770
C 00000780
CALL INV(N,W,V) 00000790
C 00000800
C GETTING E = (Y-B) 00000810
C 00000820
CALL SUB(MOBS,1,Y,B,E) 00000830
C 00000840
C GETTING Z = HT*M*E + PO IN EQ.(29) 00000850
C 00000860
CALL MUL(MOBS,MOBS,1,M,E,AA) 00000870
CALL MUL(N,MOBS,1,HT,AA,BB) 00000880
CALL ADD(N,1,BB,PO,Z) 00000890
C 00000900
C GETTING G = AMU*V*C IN EQ.(20) 00000910
C 00000920
CALL MUL(N,N,N,V,C,AA) 00000930
CALL MULCON(N,N,AMU,AA,G) 00000940
IF(IFLAGS.EQ.0) GO TO 110 00000950
C 00000960
C STORE G FOR CALCULATION OF SMOOTHED ESTIMATES 00000970
C 00000980
DO 10 I=1,N 00000990
DO 20 J=1,N 00001000
GG(I,J,T)=G(I,J) 00001010
20 CONTINUE 00001020
10 CONTINUE 00001030
110 CONTINUE 00001040
C 00001050
C GETTING QNEW = AMU*D*(I-F*G) IN EQ.(24) 00001060
C 00001070
CALL MUL(N,N,N,F,G,AA) 00001080
CALL IDEN(N,BB) 00001090
CALL SUB(N,N,BB,AA,CC) 00001100
CALL MUL(N,N,N,D,CC,DD) 00001110
CALL MULCON(N,N,AMU,DD,QNEW) 00001120
C 00001130
C GETTING PNEW = GT*Z+QNEWT*A IN EQ.(26) 00001140
C 00001150
CALL TRANS(N,N,G,AA) 00001160
CALL MUL(N,N,1,AA,Z,BB) 00001170
CALL TRANS(N,N,QNEW,CC) 00001180
CALL MUL(N,N,1,CC,A,DD) 00001190
CALL ADD(N,1,BB,DD,PNEW) 00001200
C 00001210
C GETTING S = V*(Z - AMU*C*A) IN EQ.(19) 00001220
C 00001230
CALL MUL(N,N,1,C,A,BB) 00001240
CALL MULCON(N,1,AMU,BB,CC) 00001250
CALL SUB(N,1,Z,CC,DD) 00001260
CALL MUL(N,N,1,V,DD,S) 00001270
IF(IFLAGS.EQ.0) GO TO 210 00001280
C 00001290
C STORE S FOR CALCULATION OF SMOOTHED ESTIMATES 00001300
C 00001310
DO 30 I=1,N 00001320
SS(I,T)=S(I,1) 00001330
30 CONTINUE 00001340
210 CONTINUE 00001350
IF(IFLAGR.EQ.0) GO TO 310 00001360
C 00001370
C GETTING RNEW = RO + ET*M*E + AMU*AT*D*A - ST*W*S IN EQ.(27) 00001380
C 00001390
CALL MUL(MOBS,MOBS,1,M,E,AA) 00001400
CALL TRANS(MOBS,1,E,BB) 00001410
CALL MUL(1,MOBS,1,BB,AA,CC) 00001420
CALL ADD(1,1,RO,CC,DD) 00001430
CALL MUL(N,N,1,D,A,EE) 00001440
CALL TRANS(N,1,A,FF) 00001450
CALL MUL(1,N,1,FF,EE,HH) 00001460
CALL MULCON(1,1,AMU,HH,OO) 00001470
CALL ADD(1,1,DD,OO,PP) 00001480
CALL MUL(N,N,1,W,S,QQ) 00001490
CALL TRANS(N,1,S,RR) 00001500
CALL MUL(1,N,1,RR,QQ,TT) 00001510
CALL SUB(1,1,PP,TT,RNEW) 00001520
310 CONTINUE 00001530
IF(T.EQ.TCAP) GO TO 50 00001540
C 00001550
C UPDATING QO,PO, AND RO 00001560
C 00001570
CALL SHIFT(N,N,QNEW,QO) 00001580
CALL SHIFT(N,1,PNEW,PO) 00001590
IF(IFLAGR.EQ.0) GO TO 50 00001600
CALL SHIFT(1,1,RNEW,RO) 00001610
50 CONTINUE 00001620
C 00001630
C GETTING THE FLS FILTER ESTIMATE FOR XTCAP = UINV*Z IN EQ.(30) 00001640
C 00001650
CALL INV(N,U,AA) 00001660
CALL MUL(N,N,1,AA,Z,XTCAP) 00001670
DO 65 I=1,N 00001680
X(I,TCAP)=XTCAP(I,1) 00001690
65 CONTINUE 00001700
IF (IFLAGS.EQ.1) GOTO 410 00001710
C 00001720
C PRINTING OUT THE FLS FILTER ESTIMATE FOR XTCAP 00001730
C 00001740
CALL OUTPUT(TCAP,N,X,TRUEX) 00001750
IF(IFLAGS.EQ.0) GOTO 510 00001760
410 CONTINUE 00001770
C 00001780
C GETTING SMOOTHED ESTIMATES FOR X1,..., XTCAP-1 IN EQS.(33A) 00001790
C 00001800
TCAP1=TCAP-1 00001810
DO 70 T=1,TCAP1 00001820
L=TCAP-T 00001830
DO 80 I=1,N 00001840
X(I,L)=SS(I,L) 00001850
DO 90 J=1,N 00001860
X(I,L)=X(I,L)+GG(I,J,L)*X(J,L+1) 00001870
90 CONTINUE 00001880
80 CONTINUE 00001890
70 CONTINUE 00001900
C 00001910
C PRINTING OUT THE FLS ESTIMATES FOR X1,...,XTCAP 00001920
C 00001930
DO 150 T=1,TCAP 00001940
CALL OUTPUT(T,N,X,TRUEX) 00001950
150 CONTINUE 00001960
C VALIDATION TEST: HOW WELL DO THE FLS ESTIMATES SATISFY THE 00001970
C FIRST-ORDER CONDITIONS FOR THE COST MINIMIZATION PROBLEM (5) 00001980
CALL FOCTST(X,YY) 00001990
510 CONTINUE 00002000
STOP 00002010
END 00002020
C 00002030
C MATRIX SUBROUTINES FOR ADDITION, MULTIPLICATION, TRANSPOSITION, 00002040
C SUBTRACTION, INVERSION, MULTIPLICATION BY A SCALAR, SHIFT, AND 00002050
C FORMATION OF AN IDENTITY MATRIX 00002060
C 00002070
C OBTAINING THE SUM C=A+B OF TWO NROW X MCOL MATRICES A AND B 00002080
C 00002090
SUBROUTINE ADD(NROW,MCOL,A,B,C) 00002100
IMPLICIT REAL*8(A-H,O-Z) 00002110
DIMENSION A(15,15),B(15,15),C(15,15) 00002120
DO 10 I=1,NROW 00002130
DO 20 J=1,MCOL 00002140
C(I,J)=A(I,J)+B(I,J) 00002150
20 CONTINUE 00002160
10 CONTINUE 00002170
RETURN 00002180
END 00002190
C 00002200
C OBTAINING THE PRODUCT C=A*B OF AN NROW X L MATRIX A AND AN 00002210
C L X MCOL MATRIX B 00002220
C 00002230
SUBROUTINE MUL(NROW,L,MCOL,A,B,C) 00002240
IMPLICIT REAL*8(A-H,O-Z) 00002250
DIMENSION A(15,15),B(15,15),C(15,15) 00002260
DO 10 I=1,NROW 00002270
DO 20 J=1,MCOL 00002280
SUM=0.0D+00 00002290
DO 30 K=1,L 00002300
SUM=SUM+A(I,K)*B(K,J) 00002310
30 CONTINUE 00002320
C(I,J)=SUM 00002330
20 CONTINUE 00002340
10 CONTINUE 00002350
RETURN 00002360
END 00002370
C 00002380
C OBTAINING THE TRANSPOSE B OF AN NROW X MCOL MATRIX A 00002390
C 00002400
SUBROUTINE TRANS(NROW,MCOL,A,B) 00002410
IMPLICIT REAL*8(A-H,O-Z) 00002420
DIMENSION A(15,15),B(15,15) 00002430
DO 10 I=1,NROW 00002440
DO 20 J=1,MCOL 00002450
B(J,I)=A(I,J) 00002460
20 CONTINUE 00002470
10 CONTINUE 00002480
RETURN 00002490
END 00002500
C 00002510
C OBTAINING THE DIFFERENCE C=A-B BETWEEN NROW X MCOL MATRICES 00002520
C A AND B 00002530
C 00002540
SUBROUTINE SUB(NROW,MCOL,A,B,C) 00002550
IMPLICIT REAL*8(A-H,O-Z) 00002560
DIMENSION A(15,15),B(15,15),C(15,15) 00002570
DO 10 I=1,NROW 00002580
DO 20 J=1,MCOL 00002590
C(I,J)=A(I,J)-B(I,J) 00002600
20 CONTINUE 00002610
10 CONTINUE 00002620
RETURN 00002630
END 00002640
C 00002650
C OBTAINING THE INVERSE C OF A K X K MATRIX A 00002660
C 00002670
SUBROUTINE INV(K,A,C) 00002680
IMPLICIT REAL*8(A-H,O-Z) 00002690
DIMENSION A(15,15),B(15,30),C(15,15) 00002700
DO 5 J=1,K 00002710
DO 6 I=1,K 00002720
B(I,J)=A(I,J) 00002730
6 CONTINUE 00002740
5 CONTINUE 00002750
K2=K*2 00002760
DO 7 J=1,K 00002770
DO 8 I=1,K 00002780
B(I,K+J)=0.0D+00 00002790
IF(I.EQ.J) B(I,K+J)=1.0D+00 00002800
8 CONTINUE 00002810
7 CONTINUE 00002820
C 00002830
C THE PIVOT OPERATION STARTS HERE 00002840
C 00002850
DO 9 L=1,K 00002860
PIVOT = B(L,L) 00002870
DO 13 J=L,K2 00002880
B(L,J)=B(L,J)/PIVOT 00002890
13 CONTINUE 00002900
C 00002910
C TO IMPROVE THE ROWS 00002920
C 00002930
DO 14 I=1,K 00002940
IF(I.EQ.L) GO TO 14 00002950
AIL=B(I,L) 00002960
DO 15 J=L,K2 00002970
B(I,J)=B(I,J)-AIL*B(L,J) 00002980
15 CONTINUE 00002990
14 CONTINUE 00003000
9 CONTINUE 00003010
DO 45 I=1,K 00003020
DO 46 J=1,K 00003030
C(I,J)=B(I,K+J) 00003040
46 CONTINUE 00003050
45 CONTINUE 00003060
RETURN 00003070
END 00003080
C 00003090
C OBTAINING THE PRODUCT C*A OF A SCALAR C AND AN NROW X MCOL 00003100
C MATRIX A 00003110
C 00003120
SUBROUTINE MULCON(NROW,MCOL,C,A,CA) 00003130
IMPLICIT REAL*8(A-H,O-Z) 00003140
DIMENSION A(15,15),CA(15,15) 00003150
DO 10 I=1,NROW 00003160
DO 20 J=1,MCOL 00003170
CA(I,J)=C*A(I,J) 00003180
20 CONTINUE 00003190
10 CONTINUE 00003200
RETURN 00003210
END 00003220
C 00003230
C PUTTING AN NROW X MCOL MATRIX A INTO AN NROW X MCOL MATRIX B 00003240
C 00003250
SUBROUTINE SHIFT(NROW,MCOL,A,B) 00003260
IMPLICIT REAL*8(A-H,O-Z) 00003270
DIMENSION A(15,15),B(15,15) 00003280
DO 10 I=1,NROW 00003290
DO 20 J=1,MCOL 00003300
B(I,J)=A(I,J) 00003310
20 CONTINUE 00003320
10 CONTINUE 00003330
RETURN 00003340
END 00003350
C 00003360
C FORMING THE N X N IDENTITY MATRIX E 00003370
C 00003380
SUBROUTINE IDEN(N,E) 00003390
IMPLICIT REAL*8(A-H,O-Z) 00003400
DIMENSION E(15,15) 00003410
ZERO=0.0D+00 00003420
ONE=1.0D+00 00003430
DO 10 I=1,N 00003440
DO 20 J=1,N 00003450
E(I,J)=ZERO 00003460
20 CONTINUE 00003470
10 CONTINUE 00003480
DO 30 L=1,N 00003490
E(L,L)=ONE 00003500
30 CONTINUE 00003510
RETURN 00003520
END 00003530
C 00003540
SUBROUTINE INPUT(AMU,TCAP,N,MOBS,QZERO,PZERO,RZERO,IFLAGR,IFLAGS) 00003550
IMPLICIT REAL*8(A-H,O-Z) 00003560
INTEGER TCAP 00003570
DIMENSION QZERO(15,15),PZERO(15,15),RZERO(15,15) 00003580
AMU = 1.0D+00 00003590
TCAP = 30 00003600
N = 2 00003610
MOBS = 1 00003620
DO 10 J = 1,N 00003630
DO 20 I = 1,N 00003640
QZERO(I,J) = 0.0D+00 00003650
PZERO(I,J) = 0.0D+00 00003660
RZERO(I,J) = 0.0D+00 00003670
20 CONTINUE 00003680
10 CONTINUE 00003690
IFLAGR=1 00003700
IFLAGS=1 00003710
RETURN 00003720
END 00003730
C 00003740
C 00003750
SUBROUTINE MODEL(T,F,A,H,B,D,M,Y,TRUEX) 00003760
IMPLICIT REAL*8(A-H,O-Z) 00003770
REAL*8 M 00003780
INTEGER T,TCAP 00003790
DIMENSION F(15,15),A(15,15),H(15,15),B(15,15),D(15,15),M(15,15) 00003800
DIMENSION Y(15,15),TRUEX(15,110),ZERO(15,15) 00003810
DIMENSION QZERO(15,15),PZERO(15,15),RZERO(15,15) 00003820
C 00003830
C REGIME SHIFT EXPERIMENT WITH POSSIBLY NOISY OBSERVATIONS 00003840
C 00003850
CALL INPUT(AMU,TCAP,N,MOBS,QZERO,PZERO,RZERO,IFLAGR,IFLAGS) 00003860
SIGMA = 0.00D+00 00003870
DO 10 I=1,15 00003880
DO 20 J=1,15 00003890
ZERO(I,J) = 0.0D+00 00003900
20 CONTINUE 00003910
10 CONTINUE 00003920
CALL IDEN(N,F) 00003930
CALL SHIFT(N,1,ZERO,A) 00003940
H(1,1)=1.0D+00 00003950
H(1,2)=1.0D+00 00003960
AT=DFLOAT(T) 00003970
IF(T.EQ.1) GO TO 200 00003980
H(1,1)=DSIN(10.0D+00+(AT))+0.01D+00 00003990
H(1,2)=DCOS(10.0D+00+(AT)) 00004000
200 CONTINUE 00004010
CALL SHIFT(MOBS,1,ZERO,B) 00004020
CALL IDEN(N,D) 00004030
CALL IDEN(MOBS,M) 00004040
IF (T.GT.15) GOTO 150 00004050
TRUEX(1,T) = 2.0D+00 00004060
TRUEX(2,T) = 3.0D+00 00004070
GOTO 175 00004080
150 TRUEX(1,T) = 4.0D+00 00004090
TRUEX(2,T) = 5.0D+00 00004100
175 CONTINUE 00004110
C THE DISCREPANCY TERM UU IN THE NEXT LINE IS SET TO ZERO;
C IT CAN BE REPLACED BY A CALL TO A PSEUDO-RANDOM NUMBER GENERATOR
C TO TEST GFLS IN THE PRESENCE OF NOISE IN THE OBSERVATIONS
UU = 0.0D+00 00004120
Y(1,1)=H(1,1)*TRUEX(1,T) + H(1,2)*TRUEX(2,T) + UU 00004130
RETURN 00004140
END 00004150
C 00004160
SUBROUTINE OUTPUT(T,N,X,TRUEX) 00004170
IMPLICIT REAL*8(A-H,O-Z) 00004180
INTEGER T 00004190
DIMENSION X(15,110),TRUEX(15,110) 00004200
L = T 00004210
WRITE(6,100) L,(X(I,L),I=1,N) 00004220
100 FORMAT(1H0,'TIME EQUALS',I3/1X,'FLS ESTIMATES',7X,2D25.10) 00004230
WRITE(6,200) (TRUEX(I,L),I=1,N) 00004240
200 FORMAT(1X,'TRUE X VALUES',7X,2D25.10) 00004250
RETURN 00004260
END 00004270
C 00004280
C VALIDATION TEST: HOW WELL DO THE FLS ESTIMATES SATISFY THE 00004290
C FIRST-ORDER CONDITIONS FOR THE COST MINIMIZATION PROBLEM (5) 00004300
C 00004310
SUBROUTINE FOCTST(X,YY) 00004320
IMPLICIT REAL*8(A-H,O-Z) 00004330
INTEGER T,TP1,TCAP,TCAP1 00004340
REAL*8 M,MH 00004350
DIMENSION QZERO(15,15),PZERO(15,15),RZERO(15,15) 00004360
DIMENSION XT(15,15),X(15,110),XTT(15,15),E(15,15) 00004370
DIMENSION PZEROT(15,15),EE(15,15),CO(15,15),YT(15,15),YY(15,110) 00004380
DIMENSION F(15,15),A(15,15),H(15,15),B(15,15),D(15,15),M(15,15) 00004390
DIMENSION Y(15,15),TRUEX(15,110) 00004400
DIMENSION MH(15,15),EM(15,15),EMT(15,15),W(15,15),XTP1(15,15) 00004410
DIMENSION ED(15,15),EDT(15,15),U(15,15),V(15,15),FOCD(15,15) 00004420
C = -1.0D+00 00004430
C FORM THE STATE VECTOR FOR TIME T = 1 00004440
CALL INPUT(AMU,TCAP,N,MOBS,QZERO,PZERO,RZERO,IFLAGR,IFLAGS) 00004450
DO 100 I=1,N 00004460
XT(I,1) = X(I,1) 00004470
100 CONTINUE 00004480
C FORM THE INITIAL INCREMENTAL COST CO = -(X1'QO - PO') 00004490
CALL TRANS(N,1,XT,XTT) 00004500
CALL MUL(1,N,N,XTT,QZERO,E) 00004510
CALL TRANS(N,1,PZERO,PZEROT) 00004520
CALL SUB(1,N,E,PZEROT,EE) 00004530
CALL MULCON(1,N,C,EE,CO) 00004540
C DO LOOP FOR THE SEQUENTIAL CHECK OF THE FOC FOR T=1,TCAP 00004550
DO 200 T=1,TCAP 00004560
C FORM THE TIME-T STATE VECTOR XT 00004570
DO 300 I=1,N 00004580
XT(I,1) = X(I,T) 00004590
300 CONTINUE 00004600
C FORM THE TIME-T OBSERVATION VECTOR YT 00004610
DO 400 J=1,MOBS 00004620
YT(J,1) = YY(J,T) 00004630
400 CONTINUE 00004640
CALL MODEL(T,F,A,H,B,D,M,Y,TRUEX) 00004650
C FORM W = (YT - H(T)XT - B(T))'M(T)H(T) 00004660
CALL MUL(MOBS,MOBS,N,M,H,MH) 00004670
CALL RME(N,MOBS,YT,XT,H,B,EM) 00004680
CALL TRANS(MOBS,1,EM,EMT) 00004690
CALL MUL(1,MOBS,N,EMT,MH,W) 00004700
IF(T.EQ.TCAP) GOTO 600 00004710
C FORM THE TIME-T+1 STATE VECTOR XTP1 00004720
TP1 = T + 1 00004730
DO 500 I=1,N 00004740
XTP1(I,1) = X(I,TP1) 00004750
500 CONTINUE 00004760
C FORM U = AMU*(XTP1 - F(T)XT - A(T))'*D(T) 00004770
CALL RDE(N,XTP1,XT,F,A,ED) 00004780
CALL TRANS(N,1,ED,EDT) 00004790
CALL MUL(1,N,N,EDT,D,E) 00004800
CALL MULCON(1,N,AMU,E,U) 00004810
C FORM V = U*F 00004820
CALL MUL(1,N,N,U,F,V) 00004830
GOTO 800 00004840
600 CONTINUE 00004850
DO 700 I=1,N 00004860
V(1,I) = 0.0D+00 00004870
700 CONTINUE 00004880
800 CONTINUE 00004890
C DETERMINE THE FOC DISCREPANCIES FOR TIME T 00004900
C GIVEN BY FOCD = CO + V + W 00004910
CALL ADD(1,N,CO,V,E) 00004920
CALL ADD(1,N,E,W,FOCD) 00004930
C PRINT OUT THE FOC DISCREPANCIES FOCD FOR TIME T 00004940
WRITE (6,36) T 00004950
36 FORMAT(1H0,'FOC DISCREPANCIES FOR TIME',I3) 00004960
WRITE (6,37) (FOCD(1,I),I=1,N) 00004970
37 FORMAT(1X,13D10.2) 00004980
C UPDATE THE INITIAL INCREMENTAL COST CO 00004990
CALL MULCON(1,N,C,U,CO) 00005000
200 CONTINUE 00005010
RETURN 00005020
END 00005030
C 00005040
C SUBROUTINE FOR EVALUATING THE MEASUREMENT SPECIFICATION ERROR 00005050
C EM = (YT - H(T)XT - B(T)) FOR TIME T 00005060
C 00005070
SUBROUTINE RME(N,MOBS,YT,XT,H,B,EM) 00005080
IMPLICIT REAL*8(A-H,O-Z) 00005090
DIMENSION YT(15,15),XT(15,15),H(15,15),B(15,15),EM(15,15) 00005100
DIMENSION HX(15,15),HXPB(15,15) 00005110
CALL MUL(MOBS,N,1,H,XT,HX) 00005120
CALL ADD(MOBS,1,HX,B,HXPB) 00005130
CALL SUB(MOBS,1,YT,HXPB,EM) 00005140
RETURN 00005150
END 00005160
C 00005170
C SUBROUTINE FOR EVALUATING THE DYNAMIC SPECIFICATION ERROR 00005180
C ED = (XTP1 - F(T)XT - A(T)) FOR TIME T 00005190
C 00005200
SUBROUTINE RDE(N,XTP1,XT,F,A,ED) 00005210
IMPLICIT REAL*8(A-H,O-Z) 00005220
DIMENSION XTP1(15,15),XT(15,15),F(15,15),A(15,15),ED(15,15) 00005230
DIMENSION FXT(15,15),FXTPA(15,15) 00005240
CALL MUL(N,N,1,F,XT,FXT) 00005250
CALL ADD(N,1,FXT,A,FXTPA) 00005260
CALL SUB(N,1,XTP1,FXTPA,ED) 00005270
RETURN 00005280
END 00005290