C
C
      SUBROUTINE MODCOVAR ( N, IP, X, P, A, ISTAT, C, D, R )
C
C   Fast algorithm for the solution of the modified covariance least
C   squares 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) through X(N)
C
C   Output Parameters:
C
C     P    - Real linear prediction variance at order IP
C     A    - Array of complex linear prediction coefficients
C     ISTAT - Integer status indicator at time of exit
C             0 for normal exit (no numerical ill-conditioning)
C             1 if P is not a positive value
C             2 if DELTA' and GAMMA' do not lie in the range 0 to 1
C             3 if P' is not a positive value
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 array X must be dimensioned .GE. N and array A must be
C     dimensioned .GE. IP in the calling program.   Internal array  R
C     must be dimensioned .GE. IP  and arrays C,D must be dimensioned
C     .GE. IP+1 .
C
      COMPLEX X(1),A(1),C(IP+1),D(IP+1),R(IP),LAMBDA,THETA,PSI,XI
      COMPLEX EF,EB,C1,C2,C3,C4,SAVE1,SAVE2,SAVE3,SAVE4
C
C   Initialization
C
      R1=0.
      DO 10 K=2,N-1
10      R1=R1+2.*(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
      R4=1./(R1+2.*(R2+R3))
      P=R1+R2+R3
      DELTA=1.-R2*R4
      GAMMA=1.-R3*R4
      LAMBDA=CONJG(X(1)*X(N))*R4
      C(1)=X(N)*R4
      D(1)=CONJG(X(1))*R4
      ISTAT=0
      M=0
      IF (IP .NE. 0)  GO TO 1000
      P=(.5*R1+R2+R3)/FLOAT(N)
      RETURN
C
C   Main loop
C
1000  M=M+1
      SAVE1=(0.,0.)
      DO 20 K=M+1,N
20      SAVE1=SAVE1+X(K)*CONJG(X(K-M))
      SAVE1=2.*SAVE1
      R(M)=CONJG(SAVE1)
      THETA=X(N)*D(1)
      PSI=X(N)*C(1)
      XI=CONJG(X(1))*D(1)
      IF (M .EQ. 1)  GO TO 40
      DO 30 K=1,M-1
        THETA=THETA+X(N-K)*D(K+1)                  
        PSI=PSI+X(N-K)*C(K+1)                      
        XI=XI+CONJG(X(K+1))*D(K+1)                 
        R(K)=R(K)-X(N+1-M)*CONJG(X(N+1-M+K))
     *           -CONJG(X(M))*X(M-K)               
30      SAVE1=SAVE1+CONJG(R(K))*A(M-K)             
C
C   Order update of A vector
C
40    C1=-SAVE1/P
      A(M)=C1                                      
      P=P*(1.-REAL(C1)**2-AIMAG(C1)**2)            
      IF (M .EQ. 1)  GO TO 60
      DO 50 K=1,M/2
        MK=M-K
        SAVE1=A(K)
        A(K)=SAVE1+C1*CONJG(A(MK))                 
        IF (K .EQ. MK)  GO TO 50
        A(MK)=A(MK)+C1*CONJG(SAVE1)                
50      CONTINUE
60    IF (M .NE. IP)  GO TO 65
      P=.5*P/FLOAT(N-M)
      RETURN
C
C   Time update of C,D vectors and GAMMA,DELTA,LAMBDA scalars
C
65    R1=1./(DELTA*GAMMA-REAL(LAMBDA)**2-AIMAG(LAMBDA)**2)
      C1=(THETA*CONJG(LAMBDA)+PSI*DELTA)*R1
      C2=(PSI*LAMBDA+THETA*GAMMA)*R1
      C3=(XI*CONJG(LAMBDA)+THETA*DELTA)*R1
      C4=(THETA*LAMBDA+XI*GAMMA)*R1
      DO 70 K=1,(M-1)/2+1
        MK=M+1-K
        SAVE1=CONJG(C(K))
        SAVE2=CONJG(D(K))
        SAVE3=CONJG(C(MK))
        SAVE4=CONJG(D(MK))
        C(K)=C(K)+C1*SAVE3+C2*SAVE4                
        D(K)=D(K)+C3*SAVE3+C4*SAVE4                
        IF (K .EQ. MK)  GO TO 70
        C(MK)=C(MK)+C1*SAVE1+C2*SAVE2              
        D(MK)=D(MK)+C3*SAVE1+C4*SAVE2              
70      CONTINUE
      R2=REAL(PSI)**2+AIMAG(PSI)**2
      R3=REAL(THETA)**2+AIMAG(THETA)**2
      R4=REAL(XI)**2+AIMAG(XI)**2
      R5=GAMMA-(R2*DELTA+R3*GAMMA+2.*REAL(PSI*LAMBDA*CONJG(THETA)))*R1
      R2=DELTA-(R3*DELTA+R4*GAMMA+2.*REAL(THETA*LAMBDA*CONJG(XI)))*R1
      GAMMA=R5                                     
      DELTA=R2                                     
      LAMBDA=LAMBDA+C3*CONJG(PSI)+C4*CONJG(THETA)  
      IF (P .GT. 0.)  GO TO 80
      ISTAT=1
      RETURN
80    IF (DELTA .GT. 0. .AND. DELTA .LE. 1. .AND. GAMMA .GT. 0.
     *    .AND. GAMMA .LE. 1.)  GO TO 90
      ISTAT=2
      RETURN
C
C   Time update of A vector; order updates of C,D vectors aand GAMMA,
C   DELTA,LAMBDA scalars
C
90    R1=1./P
      R2=1./(DELTA*GAMMA-REAL(LAMBDA)**2-AIMAG(LAMBDA)**2) 
      EF=X(M+1)
      EB=X(N-M)
      DO 100 K=1,M
        EF=EF+A(K)*X(M+1-K)                        
100     EB=EB+CONJG(A(K))*X(N-M+K)                 
      C1=EB*R1                                     
      C2=CONJG(EF)*R1                              
      C3=(CONJG(EB)*DELTA+EF*LAMBDA)*R2
      C4=(EF*GAMMA+CONJG(EB*LAMBDA))*R2
      DO 110 K=M,1,-1
        SAVE1=A(K)
        A(K)=SAVE1+C3*C(K)+C4*D(K)                 
        C(K+1)=C(K)+C1*SAVE1                       
110     D(K+1)=D(K)+C2*SAVE1                       
      C(1)=C1
      D(1)=C2
      R3=REAL(EB)**2+AIMAG(EB)**2
      R4=REAL(EF)**2+AIMAG(EF)**2
      P=P-(R3*DELTA+R4*GAMMA+2.*REAL(EF*EB*LAMBDA))*R2  
      DELTA=DELTA-R4*R1                            
      GAMMA=GAMMA-R3*R1                            
      LAMBDA=LAMBDA+CONJG(EF*EB)*R1                
      IF (P .GT. 0.)  GO TO 120
      ISTAT=3
      RETURN
120   IF (DELTA .GT. 0. .AND. DELTA .LE. 1. .AND. GAMMA .GT. 0.
     *    .AND. GAMMA .LE. 1.)  GO TO 1000
      ISTAT=4
      RETURN
      END
