C C .................................................................. C C SUBROUTINE APMM 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 APMM(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 - RESULT VECTOR OF DIMENSION M CONTAINING C THE VALUES OF GIVEN CONTINUOUS FUNCTIONS C FOR GIVEN ARGUMENT X C X - 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 CNP, CSP, LEP, C LAP, HEP, 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 - 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 - 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 - AUXILIARY VECTOR OF DIMENSION (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. APMM 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 APMM(FCT,N,M,TOP,IHE,PIV,T,ITER,IER) C C DIMENSION TOP(1),IHE(1),PIV(1),T(1) DOUBLE PRECISION DSUM 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. K=1 J=L+1 DO 5 I=1,L T(K)=1. 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. 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.E38 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./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+DBLE(HELP) PIV(IDO+1)=HELP 41 IDO=IDO-1 PIV(L)=-DSUM PIV(1)=1. 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. C C START MULTIPLICATION-LOOP DO 46 I=1,L DSUM=DSUM+DBLE(PIV(I)*T(IDO)) TOL=ABS(SNGL(DSUM)) IF(TOL-HELP) 46,46,45 45 HELP=TOL 46 IDO=IDO+L C END OF MULTIPLICATION-LOOP C TOL=1.E-5*HELP IF(ABS(SNGL(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*DBLE(SAVE) TOL=1.E-5*ABS(SNGL(DSUM)) TOP(J)=TOP(J)+SNGL(DSUM) IF(ABS(TOP(J))-TOL) 52,52,53 52 TOP(J)=0. 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. 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. 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. HELP=0. 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