C
C
      SUBROUTINE COVAR ( N, IP, X, PF, AF, PB, AB, ISTAT, C, D, R )
C
C   Fast algorithm for the solution of the covariance least squares
C   normal equations.
C
C   Input Parameters:
C
C     N    - Number of data samples (integer)
C     IP   - Order of linear prediction model (integer)
C     X    - Array of complex data samples X(1) to X(N)
C
C   Output Parameters:
C
C     PF   - Real forward linear prediction variance at order IP
C     AF   - Array of complex forward linear prediction coefficients
C     PB   - Real backward linear prediction variance at order IP
C     AB    - Array of complex backward linear prediction coefficients
C     ISTAT - Integer indicating status of algorithm at time of exit:
C             0 for normal exit if no numerical ill-conditioning
C             1 if PF and PB are not positive values
C             2 if DELTA'' and GAMMA' do not lie in the range 0 to 1
C             3 if PF' and PB'' are not positive values
C             4 if DELTA and GAMMA do not lie in the range 0 to 1
C     C     - Complex working vector ( length .ge. IP+1 )
C     D     - Complex working vector ( length .ge. IP+1 )
C     R     - Comples working vector ( length .ge. IP )
C
C   Notes:
C
C     External arrays AF,AB must be dimensioned .GE. IP and array X must
C     be dimensioned .GE. N in the calling program.  Internal array  R
C     must be dimensioned .GE. IP and arrays  C,D  must be dimensioned
C     .GE. IP+1.  Program elements C(k+1) and D(k+1) correspond to text
C     elements c(k) and d(k), for k=0 to IP .
C
      COMPLEX X(1),AF(1),AB(1)
      COMPLEX C(IP+1),D(IP+1),R(IP),EF,EB,THETA,TEMP,SAVE,C1,C2,C3,C4
C
C   Initialization
C
      R1=0.
      DO 10 K=2,N-1
10      R1=R1+REAL(X(K))**2+AIMAG(X(K))**2
      R2=REAL(X(1))**2+AIMAG(X(1))**2
      R3=REAL(X(N))**2+AIMAG(X(N))**2
      PF=R1+R3
      PB=R1+R2
      R4=R1+R2+R3
      R5=1./R4
      DELTA=1.-R2*R5
      GAMMA=1.-R3*R5
      C(1)=CONJG(X(N))*R5
      D(1)=CONJG(X(1))*R5
      ISTAT=0
      M=0
      IF (IP .NE. 0)  GO TO 1000
      PF=R4/FLOAT(N)
      PB=PF
      RETURN
C
C   MAIN LOOP
C
1000  M=M+1
      R1=1./PF
      R2=1./PB
      R3=1./DELTA
      R4=1./GAMMA
C
C   Order update: AF and AB vectors ; time update: C and D vectors
C
      TEMP=(0.,0.)
      DO 20 K=M+1,N
20      TEMP=TEMP+X(K)*CONJG(X(K-M))
      R(M)=CONJG(TEMP)
      THETA=X(1)*C(M)
      IF (M .EQ. 1)  GO TO 40
      DO 30 K=1,M-1
        THETA=THETA+X(M+1-K)*C(K)              
        R(K)=R(K)-X(N+1-M)*CONJG(X(N+1-M+K))   
30      TEMP=TEMP+CONJG(R(K))*AF(M-K)
40    C1=-TEMP*R2
      C2=-CONJG(TEMP)*R1
      C3=THETA*R3
      C4=CONJG(THETA)*R4
      AF(M)=C1                                 
      AB(M)=C2                                 
      SAVE=C(M)
      C(M)=SAVE+C3*D(M)
      D(M)=D(M)+C4*SAVE
      IF (M .EQ. 1)  GO TO 60
      DO 50 K=1,M-1
        SAVE=AF(K)
        AF(K)=SAVE+C1*AB(M-K)                  
        AB(M-K)=AB(M-K)+C2*SAVE                
        SAVE=C(K)
        C(K)=SAVE+C3*D(K)                      
50      D(K)=D(K)+C4*SAVE                      
60    R5=REAL(TEMP)**2+AIMAG(TEMP)**2
      PF=PF-R5*R2                              
      PB=PB-R5*R1                              
      R5=REAL(THETA)**2+AIMAG(THETA)**2
      DELTA=DELTA-R5*R4                        
      GAMMA=GAMMA-R5*R3                        
      IF (M .NE. IP)  GO TO 65
      PF=PF/FLOAT(N-M)
      PB=PB/FLOAT(N-M)
      RETURN
65    IF (PF .GT. 0. .AND. PB .GT. 0.)  GO TO 70
      ISTAT=1
      RETURN
70    IF (DELTA .GT. 0. .AND. DELTA .LE. 1. .AND. GAMMA .GT. 0.
     *    .AND. GAMMA .LE. 1)  GO TO 80
      ISTAT=2
      RETURN
C
C   Time update:  AF and AB vectors; order update:  C and D vectors
C
80    R1=1./PF
      R2=1./PB
      R3=1./DELTA
      R4=1./GAMMA
      EF=X(M+1)
      EB=X(N-M)
      DO 90 K=1,M
        EF=EF+AF(K)*X(M+1-K)                   
90      EB=EB+AB(K)*X(N-M+K)                   
      C1=EF*R3
      C2=EB*R4
      C3=CONJG(EB)*R2
      C4=CONJG(EF)*R1
      DO 100 K=M,1,-1
        SAVE=AF(K)
        AF(K)=SAVE+C1*D(K)                     
        D(K+1)=D(K)+C4*SAVE                    
        SAVE=AB(K)
        AB(K)=SAVE+C2*C(M+1-K)                 
100     C(M+1-K)=C(M+1-K)+C3*SAVE              
      C(M+1)=C3
      D(1)=C4
      R5=REAL(EF)**2+AIMAG(EF)**2
      PF=PF-R5*R3                              
      DELTA=DELTA-R5*R1                        
      R5=REAL(EB)**2+AIMAG(EB)**2
      PB=PB-R5*R4                              
      GAMMA=GAMMA-R5*R2                        
      IF (PF .GT. 0. .AND. PB .GT. 0.)  GO TO 110
      ISTAT=3
      RETURN
110   IF (DELTA .GT. 0. .AND. DELTA .LE. 1. .AND. GAMMA .GT. 0.
     *    .AND. GAMMA .LE. 1.)  GO TO 1000
      ISTAT=4
      RETURN
      END
