CTITLEMTOEPL -- SOLVE LINEAR SYSTEM OF COMPLEX BLOCK TOEPLITZ EQUATIONS
C***********************************************************************
C                 COPYRIGHT ATLANTIC RICHFIELD COMPANY 1991            *
C***********************************************************************
C***********************************************************************
CABS  MTOEPL - SOLVE LINEAR SYSTEM OF COMPLEX BLOCK TOEPLITZ EQUATIONS
C
C     COPYRIGHT ATLANTIC RICHFIELD COMPANY 1989.
C
C     ALL RIGHTS RESERVED.  NO PART OF THIS PROGRAM MAY BE PHOTOCOPIED,
C     REPRODUCED, OR  TRANSLATED TO ANOTHER PROGRAM LANGUAGE WITHOUT THE
C     PRIOR CONSENT OF ATLANTIC RICHFIELD COMPANY.
C
C
CA    DESIGNER             D CORRIGAN
CA    AUTHOR               D CORRIGAN
CA    LANGUAGE             VS FORTRAN
CA    SYSTEM               IBM/CRAY
CA    WRITTEN              12-13-88
CA
CA    REVISED      10-05-90 D.CORRIGAN
CA                               TO AVOID ASSUMPTIONS ABOUT
CA                               THE OVERALL MAGNITUDE OF THE
CA                               CROSS CORRELATIONS BY SETTING
CA                               ZTOL = ETOL * MAX(ABS(R))
CA    REVISED      10-19-90 CLJ  RENAME FROM BCLEV, ADD TO CRAY
CA
CA    PURPOSE OF PROGRAM:
CA
CA    SOLVE A LINEAR SYSTEM WHOSE MATRIX IS
CA    COMPLEX AND BLOCK TOEPLITZ
CA
CA    AFTER E.A. ROBINSON
CA          MULTICHANNEL TIME SERIES ANALYSIS
CA          SEC 6.3 / SUBROUTINE NORMEQ
CA          (THIS REFERENCE NEEDED MODIFICATIONS:
CA                1. THE MATRICES ARE COMPLEX (HERMITEAN)
CA                2. I SOLVE FOR COLUMN VECTOR SOLUTIONS
CA
CA           SEE CLAERBOUT, FUNDAMENTALS OF GEOPHYSICAL
CA                          DATA PROCESSING, PP 140)
CA
CA    CALLING PROCEDURE:
CA       SUBROUTINE MTOEPL( N,L,R,RX,G,F,A,AP,B,BP,
CA    *                    VA,VB,DA,DB,CA,CB,CF,CI,GA,TM,* )
CA
CA    CALLING ARGUMENTS:
CA
CA    ARGUMENTS ( INPUT):
CA
CA    N   - NUMBER OF CHANNELS
CA    L   - LENGTH OF EACH N CHANNEL FILTER
CA    R   - L ORDER N TOEPLITZ MATRICES (N,N,L)
CA    RX  - MAGNITUDE OF MAXIMUM DIAGONAL ELEMENT
CA    G   - L N-COMPONENT VECTORS - THE RHS
CA
CA    ARGUMENTS (OUTPUT):
CA
CA    F   - L N-COMPONENT VECTORS - THE FILTERS
CA
CA    ARGUMENTS (WORKSPACE)
CA
CA    A   - L N X N MATRICES
CA    B   - L N X N MATRICES
CA    AP  - L N X N MATRICES
CA    BP  - L N X N MATRICES
CA    VA  - N X N MATRIX
CA    VB  - N X N MATRIX
CA    DA  - N X N MATRIX
CA    DB  - N X N MATRIX
CA    CA  - N X N MATRIX
CA    CB  - N X N MATRIX
CA    CF  - N-COMPONENT VECTOR
CA    CI  - N-COMPONENT VECTOR
CA    GA  - N-COMPONENT VECTOR
CA    TM  - N X N MATRIX
CA
C     SUBROUTINES CALLED:
C
C                              ARMVE
C                              ARSET
C                              CMFUFS (MATH ADV.)
C
CEND
C***********************************************************************
C
      SUBROUTINE MTOEPL( N,L,R,RX,G,F,A,AP,B,BP,
     *                  VA,VB,DA,DB,CA,CB,CF,CI,GA,TM,* )
C
      IMPLICIT   INTEGER(A-Z)
C
C  ---------------------------------------------------------------------
C
      COMPLEX R(N,N,L),G(N,L)
      COMPLEX F(N,L)
      COMPLEX A(N,N,L),AP(N,N,L),B(N,N,L),BP(N,N,L)
      COMPLEX VA(N,N),VB(N,N),DA(N,N),DB(N,N)
      COMPLEX CA(N,N),CB(N,N)
      COMPLEX CI(N),CF(N),GA(N)
      COMPLEX TM(N,N)
C
      REAL     ETOL,ZTOL
      REAL     RX
C
      DATA ETOL/ 1.E-5 /
C
C  ---------------------------------------------------------------------
C
C            INITIALIZE
C
      ZTOL = RX*ETOL
      NN  = N*N
      CALL ARSET( A,2*NN*L,0. )
      CALL ARSET( B,2*NN*L,0. )
      CALL ARSET( F, 2*N*L,0. )
C
      DO 200 I = 1,N
        DO 100 J = 1,N
          VA(I,J)  = R(I,J,1)
100       VB(I,J)  = R(I,J,1)
        A(I,I,1) = CMPLX(1.,0.)
200   B(I,I,1) = CMPLX(1.,0.)
C
C          FILTERS FOR STAGE 1
C
      CALL ARMVE( R(1,1,1),TM,2*NN )
      CALL CMFUFS( TM,N,N,G,N,1,ZTOL,CI,F,N,IERR )
      IF(IERR.NE.0) RETURN 1
      IF( L.EQ.1 ) RETURN
C
C  ---------------------------------------------------------------------
C
C          ITERATE TO BUILD UP LENGTH LF FILTERS
C
      DO 1000 M = 2,L
C
        CALL ARSET( DA,2*NN,0. )
        CALL ARSET( DB,2*NN,0. )
        CALL ARMVE( G(1,M),GA,2*N )
C
        DO 500 LI = 1,M
          LD  = M - LI + 1
C
          DO 400 I  = 1,N
C
            DO 350 J  = 1,N
              DO 300 K  = 1,N
                DA(I,J)   = DA(I,J) + CONJG(R(K,I,LD))*A(K,J,LI)
                DB(I,J)   = DB(I,J) + R(I,K,LD)*B(K,J,LI)
300           CONTINUE
350         CONTINUE
C
            DO 380 K  = 1,N
              GA(I)    = GA(I) - CONJG(R(K,I,LD))*F(K,LI)
380         CONTINUE
C
400       CONTINUE
C
500     CONTINUE
C
C  ---------------------------------------------------------------------
C
        CALL ARMVE( VB,TM,2*NN )
        CALL CMFUFS( TM,N,N,DA,N,N,ZTOL,CI,CA,N,IERR )
        IF(IERR.NE.0) RETURN 1
C
        CALL ARMVE( VA,TM,2*NN )
        CALL CMFUFS( TM,N,N,DB,N,N,ZTOL,CI,CB,N,IERR )
        IF(IERR.NE.0) RETURN 1
C
        CALL ARMVE( A,AP,2*NN*M )
        CALL ARMVE( B,BP,2*NN*M )
C
        DO 800 LI = 1,M
          LD  = M - LI + 1
          DO 700 I  = 1,N
            DO 650 J  = 1,N
              DO 600 K  = 1,N
                A(I,J,LI) = A(I,J,LI) - BP(I,K,LD)*CA(K,J)
                B(I,J,LI) = B(I,J,LI) - AP(I,K,LD)*CB(K,J)
600           CONTINUE
650         CONTINUE
700       CONTINUE
800     CONTINUE
C
        DO 860 I  = 1,N
          DO 840 J  = 1,N
            DO 820 K  = 1,N
              VA(I,J)   = VA(I,J)   - DB(I,K)*CA(K,J)
              VB(I,J)   = VB(I,J)   - DA(I,K)*CB(K,J)
820         CONTINUE
840       CONTINUE
860     CONTINUE
C
        CALL ARMVE( VB,TM,2*NN )
        CALL CMFUFS( TM,N,N,GA,N,1,ZTOL,CI,CF,N,IERR )
        IF(IERR.NE.0) RETURN 1
C
        DO 900 LI = 1,M
          LD  = M - LI + 1
          DO 890 J = 1,N
            DO 880 K = 1,N
              F(J,LI)  = F(J,LI) + B(J,K,LD)*CF(K)
880         CONTINUE
890       CONTINUE
900     CONTINUE
C
1000  CONTINUE
C
C  ---------------------------------------------------------------------
C
      RETURN
      END