C
C     ..................................................................
C
C        SUBROUTINE DAPMM
C
C        PURPOSE
C           APPROXIMATE A FUNCTION TABULATED IN N POINTS BY ANY LINEAR
C           COMBINATION OF M GIVEN CONTINUOUS FUNCTIONS IN THE SENSE
C           OF CHEBYSHEV.
C
C        USAGE
C           CALL DAPMM(FCT,N,M,TOP,IHE,PIV,T,ITER,IER)
C           PARAMETER FCT REQUIRES AN EXTERNAL STATEMENT IN THE
C           CALLING PROGRAM.
C
C        DESCRIPTION OF PARAMETERS
C           FCT    - NAME OF SUBROUTINE TO BE SUPPLIED BY THE USER.
C                    IT COMPUTES VALUES OF M GIVEN FUNCTIONS FOR
C                    ARGUMENT VALUE X.
C                    USAGE
C                       CALL FCT(Y,X,K)
C                    DESCRIPTION OF PARAMETERS
C                       Y   - DOUBLE PRECISION RESULT VECTOR OF DIMEN-
C                             SION M CONTAINING THE VALUES OF GIVEN
C                             CONTINUOUS FUNCTIONS FOR GIVEN ARGUMENT X
C                       X   - DOUBLE PRECISON ARGUMENT VALUE
C                       K   - AN INTEGER VALUE WHICH IS EQUAL TO M-1
C                    REMARKS
C                       IF APPROXIMATION BY NORMAL CHEBYSHEV, SHIFTED
C                       CHEBYSHEV, LEGENDRE, LAGUERRE, HERMITE POLYNO-
C                       MIALS IS DESIRED SUBROUTINES DCNP,DCSP,DLEP,
C                       DLAP,DHEP, RESPECTIVELY FROM SSP COULD BE USED.
C           N      - NUMBER OF DATA POINTS DEFINING THE FUNCTION WHICH
C                    IS TO BE APPROXIMATED
C           M      - NUMBER OF GIVEN CONTINUOUS FUNCTIONS FROM WHICH
C                    THE APPROXIMATING FUNCTION IS CONSTRUCTED.
C           TOP    - DOUBLE PRECISION VECTOR OF DIMENSION 3*N.
C                    ON ENTRY IT MUST CONTAIN FROM TOP(1) UP TO TOP(N)
C                    THE GIVEN N FUNCTION VALUES AND FROM TOP(N+1) UP
C                    TO TOP(2*N) THE CORRESPONDING NODES
C                    ON RETURN TOP CONTAINS FROM TOP(1) UP TO TOP(N)
C                    THE ERRORS AT THOSE N NODES.
C                    OTHER VALUES OF TOP ARE SCRATCH.
C           IHE    - INTEGER VECTOR OF DIMENSION 3*M+4*N+6
C           PIV    - DOUBLE PRECISION VECTOR OF DIMENSION 3*M+6.
C                    ON RETURN PIV CONTAINS AT PIV(1) UP TO PIV(M) THE
C                    RESULTING COEFFICIENTS OF LINEAR APPROXIMATION.
C           T      - DOUBLE PRECISION AUXILIARY VECTOR OF DIMENSION
C                    (M+2)*(M+2)
C           ITER   - RESULTANT INTEGER WHICH SPECIFIES THE NUMBER OF
C                    ITERATIONS NEEDED
C           IER    - RESULTANT ERROR PARAMETER CODED IN THE FOLLOWING
C                    FORM
C                     IER=0  - NO ERROR
C                     IER=1  - THE NUMBER OF ITERATIONS HAS REACHED
C                              THE INTERNAL MAXIMUM N+M
C                     IER=-1 - NO RESULT BECAUSE OF WRONG INPUT PARA-
C                              METER M OR N OR SINCE AT SOME ITERATION
C                              NO SUITABLE PIVOT COULD BE FOUND
C
C        REMARKS
C           NO ACTION BESIDES ERROR MESSAGE IN CASE M LESS THAN 1 OR
C           N LESS THAN 2.
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           THE EXTERNAL SUBROUTINE FCT MUST BE FURNISHED BY THE USER.
C
C        METHOD
C           THE PROBLEM OF APPROXIMATION A TABULATED FUNCTION BY ANY
C           LINEAR COMBINATION OF GIVEN FUNCTIONS IN THE SENSE OF
C           CHEBYSHEV (I.E. TO MINIMIZE THE MAXIMUM ERROR) IS TRANS-
C           FORMED INTO A LINEAR PROGRAMMING PROBLEM. DAPMM USES A
C           REVISED SIMPLEX METHOD TO SOLVE A CORRESPONDING DUAL
C           PROBLEM. FOR REFERENCE, SEE
C           I.BARRODALE/A.YOUNG, ALGORITHMS FOR BEST L-SUB-ONE AND
C           L-SUB-INFINITY, LINEAR APPROXIMATIONS ON A DISCRETE SET,
C           NUMERISCHE MATHEMATIK, VOL.8, ISS.3 (1966), PP.295-306.
C
C     ..................................................................
C
      SUBROUTINE DAPMM(FCT,N,M,TOP,IHE,PIV,T,ITER,IER)
C
C
      DIMENSION TOP(1),IHE(1),PIV(1),T(1)
      DOUBLE PRECISION DSUM,TOP,PIV,T,SAVE,HELP,REPI,TOL
C
C        TEST ON WRONG INPUT PARAMETERS N AND M
      IER=-1
      IF (N-1) 81,81,1
    1 IF(M) 81,81,2
C
C        INITIALIZE CHARACTERISTIC VECTORS FOR THE TABLEAU
    2 IER=0
C
C        PREPARE TOP-ROW TOP
      DO 3 I=1,N
      K=I+N
      J=K+N
      TOP(J)=TOP(K)
    3 TOP(K)=-TOP(I)
C
C        PREPARE INVERSE TRANSFORMATION MATRIX T
      L=M+2
      LL=L*L
      DO 4 I=1,LL
    4 T(I)=0.D0
      K=1
      J=L+1
      DO 5 I=1,L
      T(K)=1.D0
    5 K=K+J
C
C        PREPARE INDEX-VECTOR IHE
      DO 6 I=1,L
      K=I+L
      J=K+L
      IHE(I)=0
      IHE(K)=I
    6 IHE(J)=1-I
      NAN=N+N
      K=L+L+L
      J=K+NAN
      DO 7 I=1,NAN
      K=K+1
      IHE(K)=I
      J=J+1
    7 IHE(J)=I
C
C        SET COUNTER ITER FOR ITERATION-STEPS
      ITER=-1
    8 ITER=ITER+1
C
C        TEST FOR MAXIMUM ITERATION-STEPS
      IF(N+M-ITER) 9,9,10
    9 IER=1
      GO TO 69
C
C        DETERMINE THE COLUMN WITH THE MOST POSITIVE ELEMENT IN TOP
   10 ISE=0
      IPIV=0
      K=L+L+L
      SAVE=0.D0
C
C        START TOP-LOOP
      DO 14 I=1,NAN
      IDO=K+I
      HELP=TOP(I)
      IF(HELP-SAVE) 12,12,11
   11 SAVE=HELP
      IPIV=I
   12 IF(IHE(IDO)) 14,13,14
   13 ISE=I
   14 CONTINUE
C        END OF TOP-LOOP
C
C        IS OPTIMAL TABLEAU REACHED
      IF(IPIV) 69,69,15
C
C        DETERMINE THE PIVOT-ELEMENT FOR THE COLUMN CHOSEN UPOVE
   15 ILAB=1
      IND=0
      J=ISE
      IF(J) 21,21,34
C
C        TRANSFER K-TH COLUMN FROM T TO PIV
   16 K=(K-1)*L
      DO 17 I=1,L
      J=L+I
      K=K+1
   17 PIV(J)=T(K)
C
C        IS ANOTHER COLUMN NEEDED FOR SEARCH FOR PIVOT-ELEMENT
   18 IF(ISE) 22,22,19
   19 ISE=-ISE
C
C        TRANSFER COLUMNS IN PIV
      J=L+1
      IDO=L+L
      DO 20 I=J,IDO
      K=I+L
   20 PIV(K)=PIV(I)
   21 J=IPIV
      GO TO 34
C
C        SEARCH PIVOT-ELEMENT PIV(IND)
   22 SAVE=1.D38
      IDO=0
      K=L+1
      LL=L+L
      IND=0
C
C        START PIVOT-LOOP
      DO 29 I=K,LL
      J=I+L
      HELP=PIV(I)
      IF(HELP) 29,29,23
   23 HELP=-HELP
      IF(ISE) 26,24,26
   24 IF(IHE(J)) 27,25,27
   25 IDO=I
      GO TO 29
   26 HELP=-PIV(J)/HELP
   27 IF(HELP-SAVE) 28,29,29
   28 SAVE=HELP
      IND=I
   29 CONTINUE
C        END OF PIVOT-LOOP
C
C        TEST FOR SUITABLE PIVOT-ELEMENT
      IF(IND) 30,30,32
   30 IF(IDO) 68,68,31
   31 IND=IDO
C        PIVOT-ELEMENT IS STORED IN PIV(IND)
C
C        COMPUTE THE RECIPROCAL OF THE PIVOT-ELEMENT REPI
   32 REPI=1.D0/PIV(IND)
      IND=IND-L
C
C        UPDATE THE TOP-ROW TOP OF THE TABLEAU
      ILAB=0
      SAVE=-TOP(IPIV)*REPI
      TOP(IPIV)=SAVE
C
C        INITIALIZE J AS COUNTER FOR TOP-LOOP
      J=NAN
   33 IF(J-IPIV) 34,53,34
   34 K=0
C
C        SEARCH COLUMN IN TRANSFORMATION-MATRIX T
      DO 36 I=1,L
      IF(IHE(I)-J) 36,35,36
   35 K=I
      IF(ILAB) 50,50,16
   36 CONTINUE
C
C        GENERATE COLUMN USING SUBROUTINE FCT AND TRANSFORMATION-MATRIX
      I=L+L+L+NAN+J
      I=IHE(I)-N
      IF(I) 37,37,38
   37 I=I+N
      K=1
   38 I=I+NAN
C
C        CALL SUBROUTINE FCT
      CALL FCT(PIV,TOP(I),M-1)
C
C        PREPARE THE CALLED VECTOR PIV
      DSUM=0.D0
      IDO=M
      DO 41 I=1,M
      HELP=PIV(IDO)
      IF(K) 39,39,40
   39 HELP=-HELP
   40 DSUM=DSUM+HELP
      PIV(IDO+1)=HELP
   41 IDO=IDO-1
      PIV(L)=-DSUM
      PIV(1)=1.D0
C
C        TRANSFORM VECTOR PIV WITH ROWS OF MATRIX T
      IDO=IND
      IF(ILAB) 44,44,42
   42 K=1
   43 IDO=K
   44 DSUM=0.D0
      HELP=0.D0
C
C        START MULTIPLICATION-LOOP
      DO 46 I=1,L
      DSUM=DSUM+PIV(I)*T(IDO)
      TOL=DABS(DSUM)
      IF(TOL-HELP) 46,46,45
   45 HELP=TOL
   46 IDO=IDO+L
C        END OF MULTIPLICATION-LOOP
C
      TOL=1.D-14*HELP
      IF(DABS(DSUM)-TOL) 47,47,48
   47 DSUM=0.D0
   48 IF(ILAB) 51,51,49
   49 I=K+L
      PIV(I)=DSUM
C
C        TEST FOR LAST COLUMN-TERM
      K=K+1
      IF(K-L) 43,43,18
   50 I=(K-1)*L+IND
      DSUM=T(I)
C
C        COMPUTE NEW TOP-ELEMENT
   51 DSUM=DSUM*SAVE
      TOL=1.D-14*DABS(DSUM)
      TOP(J)=TOP(J)+DSUM
      IF(DABS(TOP(J))-TOL) 52,52,53
   52 TOP(J)=0.D0
C
C        TEST FOR LAST TOP-TERM
   53 J=J-1
      IF(J) 54,54,33
C        END OF TOP-LOOP
C
C        TRANSFORM PIVOT-COLUMN
   54 I=IND+L
      PIV(I)=-1.D0
      DO 55 I=1,L
      J=I+L
   55 PIV(I)=-PIV(J)*REPI
C
C        UPDATE TRANSFORMATION-MATRIX T
      J=0
      DO 57 I=1,L
      IDO=J+IND
      SAVE=T(IDO)
      T(IDO)=0.D0
      DO 56 K=1,L
      ISE=K+J
   56 T(ISE)=T(ISE)+SAVE*PIV(K)
   57 J=J+L
C
C        UPDATE INDEX-VECTOR IHE
C        INITIALIZE CHARACTERISTICS
      J=0
      K=0
      ISE=0
      IDO=0
C
C        START QUESTION-LOOP
      DO 61 I=1,L
      LL=I+L
      ILAB=IHE(LL)
      IF(IHE(I)-IPIV) 59,58,59
   58 ISE=I
      J=ILAB
   59 IF(ILAB-IND) 61,60,61
   60 IDO=I
      K=IHE(I)
   61 CONTINUE
C        END OF QUESTION-LOOP
C
C        START MODIFICATION
      IF(K) 62,62,63
   62 IHE(IDO)=IPIV
      IF(ISE) 67,67,65
   63 IF(IND-J) 64,66,64
   64 LL=L+L+L+NAN
      K=K+LL
      I=IPIV+LL
      ILAB=IHE(K)
      IHE(K)=IHE(I)
      IHE(I)=ILAB
      IF(ISE) 67,67,65
   65 IDO=IDO+L
      I=ISE+L
      IHE(IDO)=J
      IHE(I)=IND
   66 IHE(ISE)=0
   67 LL=L+L
      J=LL+IND
      I=LL+L+IPIV
      ILAB=IHE(I)
      IHE(I)=IHE(J)
      IHE(J)=ILAB
C        END OF MODIFICATION
C
      GO TO 8
C
C        SET ERROR PARAMETER IER=-1 SINCE NO SUITABLE PIVOT IS FOUND
   68 IER=-1
C
C        EVALUATE FINAL TABLEAU
C        COMPUTE SAVE AS MAXIMUM ERROR OF APPROXIMATION AND
C        HELP AS ADDITIVE CONSTANCE FOR RESULTING COEFFICIENTS
   69 SAVE=0.D0
      HELP=0.D0
      K=L+L+L
      DO 73 I=1,NAN
      IDO=K+I
      J=IHE(IDO)
      IF(J) 71,70,73
   70 SAVE=-TOP(I)
   71 IF(M+J+1) 73,72,73
   72 HELP=TOP(I)
   73 CONTINUE
C
C        PREPARE T,TOP,PIV
      T(1)=SAVE
      IDO=NAN+1
      J=NAN+N
      DO 74 I=IDO,J
   74 TOP(I)=SAVE
      DO 75 I=1,M
   75 PIV(I)=HELP
C
C        COMPUTE COEFFICIENTS OF RESULTING POLYNOMIAL IN PIV(1) UP TO PI
C        AND CALCULATE ERRORS AT GIVEN NODES IN TOP(1) UP TO TOP(N)
      DO 79 I=1,NAN
      IDO=K+I
      J=IHE(IDO)
      IF(J) 76,79,77
   76 J=-J
      PIV(J)=HELP-TOP(I)
      GO TO 79
   77 IF(J-N) 78,78,79
   78 J=J+NAN
      TOP(J)=SAVE+TOP(I)
   79 CONTINUE
      DO 80 I=1,N
      IDO=NAN+I
   80 TOP(I)=TOP(IDO)
   81 RETURN
      END
