'********************************************************************* '* This program solves a Toeplitz linear system: * '* N * '* S R x = y (i=1,...,N) * '* j=1 N+i-j j i * '* ----------------------------------------------------------------- * '* SAMPLE RUN: * '* N=10 * '* R = (-10,-12,3,4,5,6,7,8,9,5,11,12,13,14,15,16,17,18,19) * '* Y = (1,1,1,1,1,1,1,1,30,1) * '* * '* The solution vector is: * '* 1 4.145257687310531 * '* 2 -7.202165439584252 * '* 3 .1361628410567307 * '* 4 .4738847986141193 * '* 5 .811606756171504 * '* 6 1.149328713728888 * '* 7 1.48705067128628 * '* 8 1.824772628843657 * '* 9 -3.637505413598958 * '* , 10 2.500216543958427 * '* * '* ----------------------------------------------------------------- * '* Reference: "Numerical Recipes by W.H. Press, B.P. Flannery, S.A. * '* Teukolsky, W.T. Vetterling, Cambridge University * '* Press, 1987" * '* * '* Basic Version By J-P Moreau, Paris * '* (www.jpmoreau.fr) * '********************************************************************* DEFINT I, N DEFDBL A-H, O-Z OPTION BASE 1 DIM R(19), X(10), Y(10) N = 10 DATA -10,-12,3,4,5,6,7,8,9,5,11,12,13,14,15,16,17,18,19 FOR i = 1 TO (2 * N - 1) READ R(i) NEXT i DATA 1,1,1,1,1,1,1,1,30,1 FOR i = 1 TO N READ Y(i) NEXT i CALL TOEPLZ(R(), X(), Y(), N) PRINT PRINT " The solution vector is:" FOR i = 1 TO N PRINT " "; i; " "; X(i) NEXT i END 'of main program SUB TOEPLZ (R(), X(), Y(), N) '--------------------------------------------------------------------- '* This program solves a Toeplitz linear system: * '* N * '* S R x = y (i=1,...,N) * '* j=1 N+i-j j i * '* * '* Input consists of vectors R of size 2N-1 and Y of size N. * '* Output is solution vector X of size N. * '* N is the size of the linear system. * '--------------------------------------------------------------------- NMAX = 100 ZERO = 0 ' NMAX is the maximum expected value of N. DIM G(NMAX), H(NMAX) 'Dim PP, PT1, PT2, QQ, QT1, QT2, SD, SGD, SGN, SHN, SXN IF R(N) = ZERO THEN GOTO 99 X(1) = Y(1) / R(N) 'initialize for the recursion. IF N = 1 THEN GOTO 100 G(1) = R(N - 1) / R(N) H(1) = R(N + 1) / R(N) FOR M = 1 TO N 'main loop over the recursion. M1 = M + 1 SXN = -Y(M1) SD = -R(N) FOR J = 1 TO M SXN = SXN + R(N + M1 - J) * X(J) SD = SD + R(N + M1 - J) * G(M - J + 1) NEXT J IF SD = ZERO THEN GOTO 99 X(M1) = SXN / SD 'whence x. FOR J = 1 TO M X(J) = X(J) - X(M1) * G(M - J + 1) NEXT J IF M1 = N THEN GOTO 100 'normal exit END IF SGN1 = -R(N - M1) 'compute numerator and denominator for G and H. SHN = -R(N + M1) SGD = -R(N) FOR J = 1 TO M SGN1 = SGN1 + R(N + J - M1) * G(J) SHN = SHN + R(N - J + M1) * H(J) SGD = SGD + R(N + J - M1) * H(M - J + 1) NEXT J IF SD = ZERO OR SGD = ZERO THEN GOTO 99 G(M1) = SGN1 / SGD 'whence G and H H(M1) = SHN / SD K = M M2 = (M + 1) / 2 PP = G(M1) QQ = H(M1) FOR J = 1 TO M2 PT1 = G(J) PT2 = G(K) QT1 = H(J) QT2 = H(K) G(J) = PT1 - PP * QT2 G(K) = PT2 - PP * QT1 H(J) = QT1 - QQ * PT2 H(K) = QT2 - QQ * PT1 K = K - 1 NEXT J NEXT M 99 PRINT PRINT " Levinson method fails: singular principal minor." 100 END SUB 'end of file toeplitz.bas