Directory of image this file is from
This file as a plain text file
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