'*********************************************************************
'* 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