*TITLES1EURK -- DECON MATRIX INVERSION ROUTINE 00000001 C*********************************************************************** C COPYRIGHT ATLANTIC RICHFIELD COMPANY 1991 * C*********************************************************************** *A AUTHOR HOOGSTRAAT 00000002 *A DESIGNER HOOGSTRAAT 00000003 *A LANGUAGE S/370 ASSEMBLER F 00000004 *A WRITTEN UNKNOWN 00000005 * REVISED MO-DA-YR 00000006 *A 00000007 *A 00000008 * 00000009 * 00000010 S1EURK CSECT 00000020 * 00000030 * TOPLITZ MATRIX INVERSION ROUTINE AS DESCRIBED BY ROBINSON 00000040 * IN HIS BOOK ABOUT DIGITAL FILTERING, UNDER THE NAME EUREKA. 00000050 * 00000060 * THE FOLLOWING IS THE FORTRAN EQUIVALENT. 00000070 * 00000080 *A CALL S1EURK (LR, R, G, F, A, E) 00000090 *A INPUT LR = LENGTH OF R I4 00000100 *A INPUT R = AUTOCORRELATION ARRAY R4 00000110 *A INPUT G = SAME AS R (USED FOR PREDICTIVE R4 00000120 *A DECON) 00000130 *A OUTPUT F = FILTER R4 00000140 *A INPUT A = WORK ARRAY R4 00000141 *A OUTPUT E = ERROR TRACE OPERATOR R4 00000142 *A 00000143 *A 00000144 *A TOPLITZ MATRIX INVERSION ROUTINE USED FOR DEVELOPMENT 00000145 *A OF A DECON OPERATOR. 00000146 *A 00000147 *A 00000148 * E(1) = 100. 00000150 * A(1) = 1. 00000160 * V = R(1) 00000170 * D = R(2) 00000180 * F(1) = G(1) / V 00000190 * Q = F(1) * R(2) 00000200 * IF ( LR .EQ. 1 ) RETURN 00000210 *C 00000220 * DO 4 L = 2, LR 00000230 * A(L) = - D / V 00000240 * IF ( L .EQ. 2 ) GO TO 2 00000250 *C 00000260 * L1 = ( L - 2 ) / 2 00000270 * L2 = L1 + 1 00000280 * IF ( L2 .LE. 1 ) GO TO 18 00000290 *C 00000300 * DO 1 J = 2, L2 00000310 * HOLD = A(J) 00000320 * K = L - J + 1 00000330 * A(J) = A(J) + A(L) * A(K) 00000340 *1 A(K) = A(K) + A(L) * HOLD 00000350 *C 00000360 *18 IF ( 2 * L1 .EQ. L - 2 ) GO TO 2 00000370 * A(L2+1) = A(L2+1) + A(L) * A(L2+1) 00000380 *C 00000390 *2 V = V + A(L) * D 00000400 * F(L) = ( G(L) - Q ) / V 00000410 * L3 = L - 1 00000420 *C 00000430 * DO 3 J = 1, L3 00000440 * K = L - J + 1 00000450 *3 F(J) = F(J) + F(L) * A(K) 00000460 *C 00000470 * IF ( L .EQ. LR ) GO TO 10 00000480 *C 00000490 * D = 0. 00000500 * Q = 0. 00000510 *C 00000520 * DO 4 I = 1, L 00000530 * E(L) = V * 100. / R(1) 00000540 * K = L - I + 2 00000550 * D = D + A(I) * R(K) 00000560 *4 Q = Q + F(I) * R(K) 00000570 *C 00000580 *10 E(L) = V * 100. / R(1) 00000590 * RETURN 00000600 * END 00000610 * 00000620 * ====================================================== 00000630 * * PROGRAM DEVELOPED BY 00000640 * * H * HOOGSTRAAT PROGRAMMING SERVICES LTD. 00000650 * * P S * BOX 20, SITE 7, SS #1 00000660 * * * * * * * * CALGARY, ALTA. CANADA PH 288-8088 00000670 * ====================================================== 00000680 * 00000690 ** 00000700 * SUBROUTINE S1EURK (LR,R,G,F,A) 00000710 * REAL R(LR),G(LR),F(LR),A(LR) 00000720 * 00000730 COPY S1REG 00000740 J EQU R1 00000750 I EQU R1 00000760 R EQU R2 00000770 G EQU R3 00000780 F EQU R4 00000790 A EQU R5 00000800 LR EQU R7 00000810 L2 EQU R9 00000820 L3 EQU R9 00000830 L EQU R11 00000840 K EQU R12 00000850 L1 EQU R13 00000860 BASE EQU R15 00000870 * 00000880 V EQU FR2 00000890 D EQU FR4 00000900 Q EQU FR6 00000910 HOLD EQU V 00000920 AL EQU D 00000930 AK EQU Q 00000940 * 00000950 USING *,BASE 00000960 STM R14,R12,12(R13) 00000970 ST R13,SAVE+4 00000980 * 00000990 LM R,A,4(R1) 00001000 L LR,0(R1) 00001010 L LR,0(LR) 00001020 SLL LR,2 00001030 * 00001040 LA R6,4 00001050 LR R8,R6 00001060 LR R10,R6 00001070 SR LR,R6 00001080 * 00001090 * V = R(1) 00001100 * 00001110 LE V,0(R) 00001120 * 00001130 LER FR0,V 00001140 ME FR0,MIN 00001150 STE FR0,LOW 00001160 * 00001170 * D = R(2) 00001180 * 00001190 LE D,4(R) 00001200 * 00001210 * A(1) = 1. 00001220 * 00001230 LE FR0,ONE 00001240 STE FR0,0(A) 00001250 * 00001260 * F(1) = G(1) / V 00001270 * 00001280 LE Q,0(G) 00001290 DER Q,V 00001300 STE Q,0(F) 00001310 * 00001320 * Q = F(1) * R(2) 00001330 * 00001340 ME Q,4(R) 00001350 * 00001360 * 00001370 * IF (LR.EQ.1) RETURN 00001380 * 00001390 LTR LR,LR 00001400 BZ RETURN 00001410 * 00001420 * DO 4 L = 2, LR 00001430 * 00001440 LA L,4 00001450 * 00001460 * A(L) = - D / V 00001470 * 00001480 LOOP41 LCER FR0,D 00001490 DER FR0,V 00001500 STE FR0,0(A,L) 00001510 * 00001520 * IF (L.EQ.2) GO TO 2 00001530 * 00001540 CR L,R6 00001550 BE STAT2 00001560 * 00001570 * L1 = (L - 2 ) / 2 00001580 * 00001590 LR L1,L 00001600 SR L1,R6 00001610 SRL L1,3 00001620 SLL L1,2 00001630 * 00001640 * L2 = L1 + 1 00001650 * 00001660 LR L2,L1 00001670 * 00001680 * IF (L2.LE.1) GO TO 5 00001690 * 00001700 LTR L1,L1 00001710 BZ STAT5 00001720 * 00001730 * DO 1 J = 2, L2 00001740 * 00001750 LR J,R6 00001760 * 00001770 * K = L - J + 1 00001780 * 00001790 LR K,L 00001800 SR K,R6 00001810 * 00001820 STE V,VHOLD 00001830 STE D,DHOLD 00001840 STE Q,QHOLD 00001850 * 00001860 * HOLD = A(J) 00001870 * 00001880 LOOP1 LE HOLD,0(A,J) 00001890 LE AL,0(A,L) 00001900 LE AK,0(A,K) 00001910 * 00001920 * A(J) = A(J) + A(L) * A(K) 00001930 * 00001940 LER FR0,AL 00001950 MER FR0,AK 00001960 AER FR0,HOLD 00001970 STE FR0,0(A,J) 00001980 * 00001990 *1 A(K) = A(K) + A(L) * HOLD 00002000 * 00002010 LER FR0,AL 00002020 MER FR0,HOLD 00002030 AER FR0,AK 00002040 STE FR0,0(A,K) 00002050 * 00002060 SR K,R6 00002070 * 00002080 BXLE J,L2-1,LOOP1 00002090 * 00002100 LE V,VHOLD 00002110 LE D,DHOLD 00002120 LE Q,QHOLD 00002130 * 00002140 *5 IF (2*L1.EQ.L-2) GO TO 2 00002150 * 00002160 STAT5 LR J,L1 00002170 SLL J,1 00002180 AR J,R6 00002190 CR J,L 00002200 BE STAT2 00002210 * 00002220 * A(L2+1) = A(L2+1) + A(L) * A(L2+1) 00002230 * 00002240 AR L2,R6 00002250 LE FR0,0(A,L) 00002260 ME FR0,0(A,L2) 00002270 AE FR0,0(A,L2) 00002280 STE FR0,0(A,L2) 00002290 * 00002300 * 00002310 *2 V = V + A(L) * D 00002320 * 00002330 STAT2 LE FR0,0(A,L) 00002340 MER FR0,D 00002350 AER V,FR0 00002360 * 00002370 CE V,LOW 00002380 BH NOROUND 00002390 LE V,LOW 00002400 * 00002410 * F(L) = ( G(L) - Q ) / V 00002420 * 00002430 NOROUND LE FR0,0(G,L) 00002440 SER FR0,Q 00002450 DER FR0,V 00002460 STE FR0,0(F,L) 00002470 * 00002480 * L3 = L - 1 00002490 * 00002500 LR L3,L 00002510 SR L3,R6 00002520 * 00002530 * DO 3 J = 1, L3 00002540 * 00002550 LA J,0 00002560 * 00002570 * K = L - J + 1 00002580 * 00002590 LR K,L 00002600 * 00002610 *3 F(J) = F(J) + F(L) * A(K) 00002620 * 00002630 LOOP3 LE FR0,0(A,K) 00002640 ME FR0,0(F,L) 00002650 AE FR0,0(F,J) 00002660 STE FR0,0(F,J) 00002670 * 00002680 SR K,R6 00002690 * 00002700 BXLE J,L3-1,LOOP3 00002710 * 00002720 * IF (L.EQ.LR) RETURN 00002730 * 00002740 CR L,LR 00002750 BE RETURN 00002760 * 00002770 * D = 0.0 00002780 * 00002790 LE D,ZERO 00002800 * 00002810 * Q = 0.0 00002820 * 00002830 LE Q,ZERO 00002840 * 00002850 * DO 4 I = 1, L 00002860 * 00002870 LA I,0 00002880 * 00002890 * K = L - I + 2 00002900 * 00002910 LR K,L 00002920 AR K,R6 00002930 * 00002940 * D = D + A(I) * R(K) 00002950 * 00002960 LOOP42 LE FR0,0(A,I) 00002970 ME FR0,0(R,K) 00002980 AER D,FR0 00002990 * 00003000 *4 Q = Q + F(I) * R(K) 00003010 * 00003020 LE FR0,0(F,I) 00003030 ME FR0,0(R,K) 00003040 AER Q,FR0 00003050 * 00003060 SR K,R6 00003070 * 00003080 BXLE I,L-1,LOOP42 00003090 * 00003100 BXLE L,LR-1,LOOP41 00003110 * 00003120 * END 00003130 * 00003140 RETURN L R13,SAVE+4 00003150 LM R14,R12,12(R13) 00003160 BR R14 00003170 * 00003180 ZERO DC E'0.0' 00003190 ONE DC E'1.0' 00003200 M8 DC F'8' 00003210 VHOLD DS 1D 00003220 DHOLD DS 1D 00003230 QHOLD DS 1D 00003240 SAVE DS 9D 00003250 MIN DC E'0.00001' 00003260 LOW DC E'0.00001' 00003270 LTORG 00003280 END 00003290