
      SUBROUTINE CHEB(M, N, MDIM, NDIM, A, B, TOL, RELERR, X, RANK,     CHE   1
0
     * RESMAX, ITER, OCODE)
C THIS SUBROUTINE USES A MODIFICATION OF THE SIMPLEX METHOD
C OF LINEAR PROGRAMMING TO CALCULATE A CHEBYSHEV SOLUTION TO
C AN OVER-DETERMINED SYSTEM OF LINEAR EQUATIONS.
C DESCRIPTION OF PARAMETERS.
C M      NUMBER OF EQUATIONS.
C N      NUMBER OF UNKNOWNS (N MUST NOT EXCEED M).
C MDIM   THE NUMBER OF COLUMNS OF A, AT LEAST M+1.
C NDIM   THE NUMBER OF ROWS OF A, AT LEAST N+3.
C A      TWO DIMENSIONAL REAL ARRAY OF SIZE (NDIM,MDIM).
C        ON ENTRY,THE TRANSPOSE OF THE MATRIX OF
C        COEFFICIENTS OF THE OVER-DETERMINED SYSTEM MUST
C        BE STORED IN THE FIRST M COLUMNS AND N ROWS OF A.
C        THESE VALUES ARE DESTROYED BY THE SUBROUTINE.
C B      ONE DIMENSIONAL REAL ARRAY OF SIZE MDIM. ON ENTRY,
C        B MUST CONTAIN THE RIGHT-HAND SIDES OF THE
C        EQUATIONS IN ITS FIRST M LOCATIONS. ON EXIT, B
C        CONTAINS THE RESIDUALS FOR THE EQUATIONS IN ITS
C        FIRST M LOCATIONS (SEE DESCRIPTION).
C TOL    A SMALL POSITIVE TOLERANCE. EMPIRICAL EVIDENCE
C        SUGGESTS TOL=10**(-D+1) WHERE D REPRESENTS THE
C        NUMBER OF DECIMAL DIGITS OF ACCURACY AVAILABLE
C        (SEE DESCRIPTION).
C RELERR A REAL VARIABLE WHICH ON ENTRY MUST HAVE THE VALUE
C        0.0 IF A CHEBYSHEV SOLUTION IS REQUIRED. IF RELERR
C        IS POSITIVE, THE SUBROUTINE CALCULATES AN
C        APPROXIMATE SOLUTION WITH RELERR AS AN UPPER BOUND
C        ON THE RELATIVE ERROR OF ITS LARGEST RESIDUAL (SEE
C        DESCRIPTION-INEQUALITY (2)). ON EXIT, THE VALUE OF
C        RELERR GIVES A SMALLER UPPER BOUND FOR THIS
C        RELATIVE ERROR.
C X      ONE DIMENSIONAL REAL ARRAY OF SIZE NDIM. ON EXIT,
C        THIS ARRAY CONTAINS A SOLUTION TO THE PROBLEM IN
C        ITS FIRST N LOCATIONS.
C RANK   AN INTEGER WHICH GIVES ON EXIT THE RANK OF THE
C        MATRIX OF COEFFICIENTS.
C RESMAX THE LARGEST RESIDUAL IN MAGNITUDE.
C ITER   THE NUMBER OF SIMPLEX ITERATIONS PERFORMED.
C OCODE  AN EXIT CODE WITH VALUES..
C              0 - OPTIMAL SOLUTION WHICH IS PROBABLY
C                  NON-UNIQUE (SEE DESCRIPTION).
C              1 - UNIQUE OPTIMAL SOLUTION.
C              2 - CALCULATIONS TERMINATED PREMATURELY
C                  DUE TO ROUNDING ERRORS.
C IF YOUR FORTRAN COMPILER PERMITS A SINGLE COLUMN OF A TWO
C DIMENSIONAL ARRAY TO BE PASSED TO A ONE DIMENSIONAL ARRAY
C THROUGH A SUBROUTINE CALL, CONSIDERABLE SAVINGS IN
C EXECUTION TIME MAY BE ACHIEVED THROUGH THE USE OF THE
C FOLLOWING SUBROUTINE WHICH OPERATES ON COLUMN VECTORS.
C     SUBROUTINE COL (V1,V2,MLT,NOTROW,I1,NP2)
C THIS SUBROUTINE SUBTRACTS FROM THE VECTOR V1 A MULTIPLE OF
C THE VECTOR V2 STARTING AT THE I1*TH ELEMENT UP TO THE
C NP2*TH ELEMENT, EXCEPT FOR THE NOTROW*TH ELEMENT.
C     REAL V1(NP2),V2(NP2),MLT
C     DO 1 I=I1,NP2
C       IF(I.EQ.NOTROW) GO TO 1
C       V1(I)=V1(I)-MLT*V2(I)
C   1   CONTINUE
C     RETURN
C     END
C SEE COMMENTS FOLLOWING STATEMENT NUMBER 340 FOR
C INSTRUCTIONS ON THE IMPLEMENTATION OF THIS MODIFICATION.
      DIMENSION A(NDIM,MDIM), B(MDIM), X(NDIM)
      INTEGER PROW, PCOL, RANK, RANKP1, OCODE
C BIG MUST BE SET EQUAL TO ANY VERY LARGE REAL CONSTANT.
C ITS VALUE HERE IS APPROPRIATE FOR THE IBM/370/145.
      DATA BIG /1.E+75/
C INITIALIZATION.
      MP1 = M + 1
      NP1 = N + 1
      NP2 = N + 2
      NP3 = N + 3
      NP1MR = 1
      RANK = N
      RELTMP = RELERR
      RELERR = 0.
      DO 10 J=1,M
        A(NP1,J) = 1.
        A(NP2,J) = -B(J)
        A(NP3,J) = N + J
   10 CONTINUE
      A(NP1,MP1) = 0.
      ITER = 0
      OCODE = 1
      DO 20 I=1,N
        X(I) = 0.
        A(I,MP1) = I
   20 CONTINUE
C LEVEL 1.
      LEV = 1
      K = 0
   30 K = K + 1
      KP1 = K + 1
      NP1MK = NP1 - K
      MODE = 0
      DO 40 J=K,M
        B(J) = 1.
   40 CONTINUE
C DETERMINE THE VECTOR TO ENTER THE BASIS.
   50 D = -BIG
      DO 60 J=K,M
        IF (B(J).EQ.0.) GO TO 60
        DD = ABS(A(NP2,J))
        IF (DD.LE.D) GO TO 60
        PCOL = J
        D = DD
   60 CONTINUE
      IF (K.GT.1) GO TO 70
C TEST FOR ZERO RIGHT-HAND SIDE.
      IF (D.GT.TOL) GO TO 70
      RESMAX = 0.
      MODE = 2
      GO TO 380
C DETERMINE THE VECTOR TO LEAVE THE BASIS.
   70 D = TOL
      DO 80 I=1,NP1MK
        DD = ABS(A(I,PCOL))
        IF (DD.LE.D) GO TO 80
        PROW = I
        D = DD
   80 CONTINUE
      IF (D.GT.TOL) GO TO 330
C CHECK FOR LINEAR DEPENDENCE IN LEVEL 1.
      B(PCOL) = 0.
      IF (MODE.EQ.1) GO TO 50
      DO 100 J=K,M
        IF (B(J).EQ.0.) GO TO 100
        DO 90 I=1,NP1MK
          IF (ABS(A(I,J)).LE.TOL) GO TO 90
          MODE = 1
          GO TO 50
   90   CONTINUE
  100 CONTINUE
      RANK = K - 1
      NP1MR = NP1 - RANK
      OCODE = 0
      GO TO 160
  110 IF (PCOL.EQ.K) GO TO 130
C INTERCHANGE COLUMNS IN LEVEL 1.
      DO 120 I=1,NP3
        D = A(I,PCOL)
        A(I,PCOL) = A(I,K)
        A(I,K) = D
  120 CONTINUE
  130 IF (PROW.EQ.NP1MK) GO TO 150
C INTERCHANGE ROWS IN LEVEL 1.
      DO 140 J=1,MP1
        D = A(PROW,J)
        A(PROW,J) = A(NP1MK,J)
        A(NP1MK,J) = D
  140 CONTINUE
  150 IF (K.LT.N) GO TO 30
  160 IF (RANK.EQ.M) GO TO 380
      RANKP1 = RANK + 1
C LEVEL 2.
      LEV = 2
C DETERMINE THE VECTOR TO ENTER THE BASIS
      D = TOL
      DO 170 J=RANKP1,M
        DD = ABS(A(NP2,J))
        IF (DD.LE.D) GO TO 170
        PCOL = J
        D = DD
  170 CONTINUE
C COMPARE CHEBYSHEV ERROR WITH TOL.
      IF (D.GT.TOL) GO TO 180
      RESMAX = 0.
      MODE = 3
      GO TO 380
  180 IF (A(NP2,PCOL).LT.-TOL) GO TO 200
      A(NP1,PCOL) = 2. - A(NP1,PCOL)
      DO 190 I=NP1MR,NP3
        IF (I.EQ.NP1) GO TO 190
        A(I,PCOL) = -A(I,PCOL)
  190 CONTINUE
C ARRANGE FOR ALL ENTRIES IN PIVOT COLUMN
C (EXCEPT PIVOT) TO BE NEGATIVE.
  200 DO 220 I=NP1MR,N
        IF (A(I,PCOL).LT.TOL) GO TO 220
        DO 210 J=1,M
          A(NP1,J) = A(NP1,J) + 2.*A(I,J)
          A(I,J) = -A(I,J)
  210   CONTINUE
        A(I,MP1) = -A(I,MP1)
  220 CONTINUE
      PROW = NP1
      GO TO 330
  230 IF (RANKP1.EQ.M) GO TO 380
      IF (PCOL.EQ.M) GO TO 250
C INTERCHANGE COLUMNS IN LEVEL 2.
      DO 240 I=NP1MR,NP3
        D = A(I,PCOL)
        A(I,PCOL) = A(I,M)
        A(I,M) = D
  240 CONTINUE
  250 MM1 = M - 1
C LEVEL 3.
      LEV = 3
C DETERMINE THE VECTOR TO ENTER THE BASIS.
  260 D = -TOL
      VAL = 2.*A(NP2,M)
      DO 280 J=RANKP1,MM1
        IF (A(NP2,J).GE.D) GO TO 270
        PCOL = J
        D = A(NP2,J)
        MODE = 0
        GO TO 280
  270   DD = VAL - A(NP2,J)
        IF (DD.GE.D) GO TO 280
        MODE = 1
        PCOL = J
        D = DD
  280 CONTINUE
      IF (D.GE.-TOL) GO TO 380
      DD = -D/A(NP2,M)
      IF (DD.GE.RELTMP) GO TO 290
      RELERR = DD
      MODE = 4
      GO TO 380
  290 IF (MODE.EQ.0) GO TO 310
      DO 300 I=NP1MR,NP1
        A(I,PCOL) = 2.*A(I,M) - A(I,PCOL)
  300 CONTINUE
      A(NP2,PCOL) = D
      A(NP3,PCOL) = -A(NP3,PCOL)
C DETERMINE THE VECTOR TO LEAVE THE BASIS.
  310 D = BIG
      DO 320 I=NP1MR,NP1
        IF (A(I,PCOL).LE.TOL) GO TO 320
        DD = A(I,M)/A(I,PCOL)
        IF (DD.GE.D) GO TO 320
        PROW = I
        D = DD
  320 CONTINUE
      IF (D.LT.BIG) GO TO 330
      OCODE = 2
      GO TO 380
C PIVOT ON A(PROW,PCOL).
  330 PIVOT = A(PROW,PCOL)
      DO 340 J=1,M
        A(PROW,J) = A(PROW,J)/PIVOT
  340 CONTINUE
C IF PERMITTED, USE SUBROUTINE COL IN THE DESCRIPTION
C SECTION AND REPLACE THE FOLLOWING EIGHT STATEMENTS DOWN TO
C AND INCLUDING STATEMENT NUMBER 360 BY..
C     DO 360 J=1,M
C       IF(J.EQ.PCOL) GO TO 360
C       CALL COL (A(1,J),A(1,PCOL),A(PROW,J),PROW,NP1MR,NP2)
C 360 CONTINUE
      DO 360 J=1,M
        IF (J.EQ.PCOL) GO TO 360
        D = A(PROW,J)
        DO 350 I=NP1MR,NP2
          IF (I.EQ.PROW) GO TO 350
          A(I,J) = A(I,J) - D*A(I,PCOL)
  350   CONTINUE
  360 CONTINUE
      TPIVOT = -PIVOT
      DO 370 I=NP1MR,NP2
        A(I,PCOL) = A(I,PCOL)/TPIVOT
  370 CONTINUE
      A(PROW,PCOL) = 1./PIVOT
      D = A(PROW,MP1)
      A(PROW,MP1) = A(NP3,PCOL)
      A(NP3,PCOL) = D
      ITER = ITER + 1
      GO TO (110, 230, 260), LEV
C PREPARE OUTPUT.
  380 DO 390 J=1,M
        B(J) = 0.
  390 CONTINUE
      IF (MODE.EQ.2) GO TO 450
      DO 400 J=1,RANK
        K = A(NP3,J)
        X(K) = A(NP2,J)
  400 CONTINUE
      IF (MODE.EQ.3 .OR. RANK.EQ.M) GO TO 450
      DO 410 I=NP1MR,NP1
        K = ABS(A(I,MP1)) - FLOAT(N)
        B(K) = A(NP2,M)*SIGN(1.,A(I,MP1))
  410 CONTINUE
      IF (RANKP1.EQ.M) GO TO 430
      DO 420 J=RANKP1,MM1
        K = ABS(A(NP3,J)) - FLOAT(N)
        B(K) = (A(NP2,M)-A(NP2,J))*SIGN(1.,A(NP3,J))
  420 CONTINUE
C TEST FOR NON-UNIQUE SOLUTION.
  430 DO 440 I=NP1MR,NP1
        IF (ABS(A(I,M)).GT.TOL) GO TO 440
        OCODE = 0
        GO TO 450
  440 CONTINUE
  450 IF (MODE.NE.2 .AND. MODE.NE.3) RESMAX = A(NP2,M)
      IF (RANK.EQ.M) RESMAX = 0.
      IF (MODE.EQ.4) RESMAX = RESMAX - D
      RETURN
      END

