      SUBROUTINE LEVINSON (M,T0,T,P,A,ISTAT)
C
C   Solves the set of complex linear simultaneous equations
C                           TA = P
C   by the Levinson algorithm.  T is M+1 by M+1 Hermitian Toeplitz
C   matrix,  A is column vector with elements 1,a(1),...,a(M), and
C   P  is a column vector with  p as the top element and zeros for
C   the remaining elements.
C
C   Input Parameters:
C
C      M  - Number of unknown A vector components;  dimension of A
C           is M. (This is the autoregressive model order; see
C           Chap. 7.) (integer)
C      T0 - Real scalar corresponding to matrix element t(0)
C           (Hermitian symmetry forces this element to be real.)
C      T  - Array of M complex matrix elements t(1),...,t(M) from
C           the left column of the Toeplitz matrix (These are the
C           autocorrelation lags of the autoregressive model; see
C           Chap. 7.)
C
C   Output Parameters:
C
C      P  - Real scalar of top element of right-hand-side vector
C           (This is the autoregressive model driving noise
C           variance; see Chap. 7.)
C      A  - Array of M complex elements of solution vector. (These
C           are the autoregressive coefficients; see Chap. 7.)
C      ISTAT - Integer status indicator at time of exit
C              0 for normal exit
C              1 if P=0 (singular matrix)
C
C   Notes:
C
C     External arrays A and T must be dimensioned .GE. M by calling
C     program.
C
      COMPLEX T(1),A(1),TEMP,SAVE
      REAL T0,P
      P=T0
      ISTAT=0
C   Handle M=0 as special case
      IF (M .EQ. 0)  RETURN
C
C   Main Levinson recursion
C
      K=0
100   K=K+1
      SAVE=T(K)
      IF (K .EQ. 1)  GO TO 20
      DO 10 J=1,K-1
10      SAVE=SAVE+A(J)*T(K-J)                 ! Eq. (3.136)
20    TEMP=-SAVE/P
      P=P*(1.-REAL(TEMP)**2-AIMAG(TEMP)**2)   ! Eq. (3.158)
      IF (P .GT. 0.)  GO TO 30
      ISTAT=1
      RETURN
30    A(K)=TEMP                               ! Eq. (3.139)
      IF (K .EQ. 1)  GO TO 50
      KHALF=K/2
      DO 40 J=1,KHALF
        KJ=K-J
        SAVE=A(J)
        A(J)=SAVE+TEMP*CONJG(A(KJ))           ! Eqs. (3.147),(3.156)
        IF (J .EQ. KJ)  GO TO 40
        A(KJ)=A(KJ)+TEMP*CONJG(SAVE)          ! Eqs. (3.147),(3.156)
40      CONTINUE
50    IF (K .LT. M)  GO TO 100
      RETURN
      END
