SUBROUTINE CHEB(M, N, MDIM, NDIM, A, B, TOL, RELERR, X, RANK, CHE 10 * 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