C***********************************************************************
C                                                                      *
C        CHEBYSHEV RATIONAL APPROXIMATION                              *
C                                                                      *
C***********************************************************************
C
C
C
C     To obtain the rational approximation
C
C     rT(x) = (p0*T0 + p1*T1 +...+ Pn*Tn) / (q0*T0 + q1*T1 +...+ qm*Tm)
C
C     for a given function f(x):
C
C     INPUT  nonnegative integers m and n.
C
C     OUTPUT  coefficients q0, q1, ... , qm, p0, p1, ... , pn.
C
C     The coefficients of the Chebyshev expansion a0, a1,  ... could
C     be calculated instead of input as is assumed in this program.
C
C     Note that a0 is to be doubled.
      DIMENSION A(11,12),AA(11),NROW(11),P(7),Q(7)
      INTEGER PP,BN
      CHARACTER NAME*14,NAME1*14,AAA*1
      INTEGER INP,OUP,FLAG
      LOGICAL OK
      OPEN(UNIT=5,FILE='CON',ACCESS='SEQUENTIAL')
      OPEN(UNIT=6,FILE='CON',ACCESS='SEQUENTIAL')
      ZERO = 1.0E-20
      OK = .FALSE.
      WRITE(6,*) 'This is Chebyshev Approximation.'
11    IF (OK) GOTO 12
      WRITE(6,*) 'Input m and n'
      WRITE(6,*) ' '
      READ(5,*) LM,LN
      BN = LM+LN
      IF( (LM.GE.0) .AND. (LN.GE.0)) THEN
         OK = .TRUE.
      ELSE
         WRITE(6,*) 'm and n must be nonnegative'
      ENDIF
      GOTO 11
12    IF ( (LM.EQ.0) .AND. (LN.EQ.0)) THEN
         OK = .FALSE.
         WRITE(6,*) 'Not both m and n can be zero.'
         GOTO 11
      ENDIF
      OK = .FALSE.
21    IF ( .NOT. OK) THEN
         WRITE(6,*) 'The Chebyshev coefficients a(0),...,a(N)'
         WRITE(6,*) 'are to be input.'
         WRITE(6,*) 'Choice of input method: '
         WRITE(6,*) '1. Input entry by entry from keyboard '
         WRITE(6,*) '2. Input data from a text file '
         WRITE(6,*) 'Choose 1 or 2 please '
         WRITE(6,*) ' '
         READ(5,*)  FLAG
         IF( ( FLAG .GE. 1 ) .AND. ( FLAG .LE. 2 )) OK = .TRUE.
         GOTO 21
      ENDIF
      IF (FLAG .EQ. 1) THEN
         DO 31 I = 0, BN
            WRITE(6,*) 'Input a(',I,')'
            WRITE(6,*) ' '
            READ(5,*) AA(I+1)
31       CONTINUE
      ENDIF
      IF (FLAG .EQ. 2) THEN
         WRITE(6,*) 'Has a text file been created.'
         WRITE(6,*) 'Enter Y or N - letter within quotes '
         WRITE(6,*) ' '
         READ(5,*)  AAA
         IF (( AAA .EQ. 'Y' ) .OR.( AAA .EQ. 'y' )) THEN
            WRITE(6,*) 'Input the file name in the form - '
            WRITE(6,*) 'drive:name.ext  contained in quotes'
            WRITE(6,*) 'as example:   ''A:DATA.DTA'' '
            WRITE(6,*) ' '
            READ(5,*)  NAME
            INP = 4
            OPEN(UNIT=INP,FILE=NAME,ACCESS='SEQUENTIAL')
            READ(4,*) (AA(I+1),I=0,BN)
            CLOSE(UNIT=4)
         ELSE
            WRITE(6,*) 'The program will end so the input file can '
            WRITE(6,*) 'be created. '
            OK = .FALSE.
         ENDIF
      ENDIF
      IF(.NOT.OK) GOTO 400
      WRITE(6,*) 'Select output destinations: '
      WRITE(6,*) '1. Screen '
      WRITE(6,*) '2. Text file '
      WRITE(6,*) 'Enter 1 or 2 '
      WRITE(6,*) ' '
      READ(5,*) FLAG
      IF ( FLAG .EQ. 2 ) THEN
         WRITE(6,*) 'Input the file name in the form - '
         WRITE(6,*) 'drive:name.ext'
         WRITE(6,*) 'with the name contained within quotes'
         WRITE(6,*) 'as example:   ''A:OUTPUT.DTA'' '
         WRITE(6,*) ' '
         READ(5,*) NAME1
         OUP = 3
         OPEN(UNIT=OUP,FILE=NAME1,STATUS='NEW')
      ELSE
         OUP = 6
      ENDIF
      WRITE(OUP,*) 'PADE APPROXIMATION'
      LL1 = BN+1
      LL2 = BN + LM
      DO 111 I = LL1,LL2
111      AA(I+1) = 0.0
      OK = .TRUE.
C     STEP 1
      N = BN
      M = N + 1
C     STEP 2  -  performed in input
      DO 20 I = 1, M
20       NROW(I) = I
C     initialize row pointer for linear system
      NN = N - 1
C     STEP 3
      Q(1) = 1.0
C     STEP 4
C     Set up a linear system, but use A[i,j] instead of B[i,j].
      DO 30 I = 0, N
C        STEP 5
         DO 40 J = 0, I
            IF (J .LE. LN) A(I+1,J+1) = 0.0
40       CONTINUE
C        STEP 6
         IF (I .LE. LN) A(I+1,I+1) = 1.0
C        STEP 7
         LLL = I + 1
         IF (LLL .LE. LN) THEN
         DO 50 J = LLL, LN
50          A(I+1,J+1) = 0.0
         ENDIF
C        STEP 8
         LLL = LN+1
         IF (LLL .LE. N) THEN
         DO 60 J = LLL, N
            IF (I .NE. 0) THEN
               PP = I - J + LN
               IF (PP .LT. 0) PP = -PP
               A(I+1,J+1) = -(AA(I+J-LN+1)+AA(PP+1))/2.0
            ELSE
               A(I+1,J+1) = -AA(J-LN+1)/2.0
            ENDIF
60       CONTINUE
         ENDIF
         A(I+1,N+2) = AA(I+1)
30    CONTINUE
C        STEP 9
      A(1,N+2) = A(1,N+2)/2.0
      ICHG = 0
C     Solve the linear system using partial pivoting.
      I = LN+2
      NNN = N+1
      MM = M+1
C     STEP 10
80    IF ((.NOT. OK ) .OR. ( I .GT. N )) GOTO 90
C        STEP 11
         IMAX = NROW(I)
         AMAX = ABS( A(IMAX,I) )
         IMAX = I
         JJ = I + 1
         DO 100 IP = JJ, NNN
            JP = NROW(IP)
            IF ( ABS( A(JP,I) ) .GT. AMAX ) THEN
               AMAX = ABS( A(JP,I) )
               IMAX = IP
            ENDIF
100      CONTINUE
C        STEP 12
         IF ( AMAX .LE. ZERO ) THEN
            OK = .FALSE.
         ELSE
C           STEP 13
C           simulate row interchange
            IF ( NROW(I) .NE. NROW(IMAX) ) THEN
               ICHG = ICHG + 1
               NCOPY = NROW(I)
               NROW(I) = NROW(IMAX)
               NROW(IMAX) = NCOPY
            ENDIF
            I1 = NROW(I)
C           STEP 14
C           Perform elimination.
            DO 110 J = JJ, MM
               J1 = NROW(J)
C              STEP 15
               XM = A(J1,I) / A(I1,I)
C              STEP 16
               DO 120 K = JJ, MM
120               A(J1,K) = A(J1,K) - XM * A(I1,K)
C              STEP 17
               A(J1,I) = XM
110         CONTINUE
         ENDIF
         I = I + 1
      GOTO 80
90    IF ( OK ) THEN
C        STEP 18
         N1 = NROW(NNN)
         IF ( ABS( A(N1,NNN) ) .LE. ZERO ) THEN
            OK = .FALSE.
C           system has no unique solution
         ELSE
C           STEP 19
C           Start backward substitution.
            IF (LM .GT. 0) THEN
               Q(LM+1) = A(N1,MM) / A(N1,NNN)
               A(N1,MM) = Q(LM+1)
            ENDIF
            PP = 1
C           STEP 20
            LLL = LN+2
            IF (LLL .LE. NNN) THEN
            DO 130 K = LLL, N
               I = N - K + LLL
               JJ = I + 1
               N2 = NROW(I)
               SUM = A(N2,MM)
               IF (JJ .LE. NNN) THEN
               DO 140 KK = JJ, NNN
                  LL = NROW(KK)
                  SUM = SUM - A(N2,KK) * A(LL,MM)
140            CONTINUE
               ENDIF
               A(N2,MM) = SUM / A(N2,I)
               Q(LM-PP+1) = A(N2,MM)
               PP = PP + 1
130         CONTINUE
            ENDIF
C           STEP 21
            LLL = LN+1
            DO 150 K = 1, LLL
               I = LLL - K +1
               N2 = NROW(I)
               SUM = A(N2,MM)
               LLLL = LN+2
               IF (LLLL .LE. NNN) THEN
               DO 160 KK = LLLL, NNN
                  LL = NROW(KK)
                  SUM = SUM - A(N2,KK) * A(LL,MM)
160            CONTINUE
               ENDIF
               A(N2,MM) = SUM
               P(LN-K+2) = A(N2,MM)
150         CONTINUE
C           STEP 22
C           procedure completed successfully
         ENDIF
      ENDIF
      IF ( .NOT. OK ) THEN
         WRITE(OUP,*) 'System has no unique solution '
      ELSE
         WRITE(OUP,*)  'Coefficients Q[0], ..., Q[M]'
         DO 170 I = 0, LM
170         WRITE(OUP,1) Q(I+1)
         WRITE(OUP,*)  'Coefficients P[0], ..., P[N]'
         DO 180 I = 0, LN
180         WRITE(OUP,1) P(I+1)
      ENDIF
400   CLOSE(UNIT=5)
      CLOSE(UNIT=OUP)
      IF(OUP.NE.6) CLOSE(UNIT=6)
      STOP
1     FORMAT(1X,E15.8)
      END
