home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / fortran / library / ssp / itrpapsm / apmm.for < prev    next >
Encoding:
Text File  |  1985-11-29  |  10.4 KB  |  411 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE APMM
  5. C
  6. C        PURPOSE
  7. C           APPROXIMATE A FUNCTION TABULATED IN N POINTS BY ANY LINEAR
  8. C           COMBINATION OF M GIVEN CONTINUOUS FUNCTIONS IN THE SENSE
  9. C           OF CHEBYSHEV.
  10. C
  11. C        USAGE
  12. C           CALL APMM(FCT,N,M,TOP,IHE,PIV,T,ITER,IER)
  13. C           PARAMETER FCT REQUIRES AN EXTERNAL STATEMENT IN THE
  14. C           CALLING PROGRAM.
  15. C
  16. C        DESCRIPTION OF PARAMETERS
  17. C           FCT    - NAME OF SUBROUTINE TO BE SUPPLIED BY THE USER.
  18. C                    IT COMPUTES VALUES OF M GIVEN FUNCTIONS FOR
  19. C                    ARGUMENT VALUE X.
  20. C                    USAGE
  21. C                       CALL FCT(Y,X,K)
  22. C                    DESCRIPTION OF PARAMETERS
  23. C                       Y   - RESULT VECTOR OF DIMENSION M CONTAINING
  24. C                             THE VALUES OF GIVEN CONTINUOUS FUNCTIONS
  25. C                             FOR GIVEN ARGUMENT X
  26. C                       X   - ARGUMENT VALUE
  27. C                       K   - AN INTEGER VALUE WHICH IS EQUAL TO M-1
  28. C                    REMARKS
  29. C                       IF APPROXIMATION BY NORMAL CHEBYSHEV, SHIFTED
  30. C                       CHEBYSHEV, LEGENDRE, LAGUERRE, HERMITE POLYNO-
  31. C                       MIALS IS DESIRED SUBROUTINES CNP, CSP, LEP,
  32. C                       LAP, HEP, RESPECTIVELY FROM SSP COULD BE USED.
  33. C           N      - NUMBER OF DATA POINTS DEFINING THE FUNCTION WHICH
  34. C                    IS TO BE APPROXIMATED
  35. C           M      - NUMBER OF GIVEN CONTINUOUS FUNCTIONS FROM WHICH
  36. C                    THE APPROXIMATING FUNCTION IS CONSTRUCTED.
  37. C           TOP    - VECTOR OF DIMENSION 3*N.
  38. C                    ON ENTRY IT MUST CONTAIN FROM TOP(1) UP TO TOP(N)
  39. C                    THE GIVEN N FUNCTION VALUES AND FROM TOP(N+1) UP
  40. C                    TO TOP(2*N) THE CORRESPONDING NODES
  41. C                    ON RETURN TOP CONTAINS FROM TOP(1) UP TO TOP(N)
  42. C                    THE ERRORS AT THOSE N NODES.
  43. C                    OTHER VALUES OF TOP ARE SCRATCH.
  44. C           IHE    - INTEGER VECTOR OF DIMENSION 3*M+4*N+6
  45. C           PIV    - VECTOR OF DIMENSION 3*M+6.
  46. C                    ON RETURN PIV CONTAINS AT PIV(1) UP TO PIV(M) THE
  47. C                    RESULTING COEFFICIENTS OF LINEAR APPROXIMATION.
  48. C           T      - AUXILIARY VECTOR OF DIMENSION (M+2)*(M+2)
  49. C           ITER   - RESULTANT INTEGER WHICH SPECIFIES THE NUMBER OF
  50. C                    ITERATIONS NEEDED
  51. C           IER    - RESULTANT ERROR PARAMETER CODED IN THE FOLLOWING
  52. C                    FORM
  53. C                     IER=0  - NO ERROR
  54. C                     IER=1  - THE NUMBER OF ITERATIONS HAS REACHED
  55. C                              THE INTERNAL MAXIMUM N+M
  56. C                     IER=-1 - NO RESULT BECAUSE OF WRONG INPUT PARA-
  57. C                              METER M OR N OR SINCE AT SOME ITERATION
  58. C                              NO SUITABLE PIVOT COULD BE FOUND
  59. C
  60. C        REMARKS
  61. C           NO ACTION BESIDES ERROR MESSAGE IN CASE M LESS THAN 1 OR
  62. C           N LESS THAN 2.
  63. C
  64. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  65. C           THE EXTERNAL SUBROUTINE FCT MUST BE FURNISHED BY THE USER.
  66. C
  67. C        METHOD
  68. C           THE PROBLEM OF APPROXIMATION A TABULATED FUNCTION BY ANY
  69. C           LINEAR COMBINATION OF GIVEN FUNCTIONS IN THE SENSE OF
  70. C           CHEBYSHEV (I.E. TO MINIMIZE THE MAXIMUM ERROR) IS TRANS-
  71. C           FORMED INTO A LINEAR PROGRAMMING PROBLEM. APMM USES A
  72. C           REVISED SIMPLEX METHOD TO SOLVE A CORRESPONDING DUAL
  73. C           PROBLEM. FOR REFERENCE, SEE
  74. C           I.BARRODALE/A.YOUNG, ALGORITHMS FOR BEST L-SUB-ONE AND
  75. C           L-SUB-INFINITY, LINEAR APPROXIMATIONS ON A DISCRETE SET,
  76. C           NUMERISCHE MATHEMATIK, VOL.8, ISS.3 (1966), PP.295-306.
  77. C
  78. C     ..................................................................
  79. C
  80.       SUBROUTINE APMM(FCT,N,M,TOP,IHE,PIV,T,ITER,IER)
  81. C
  82. C
  83.       DIMENSION TOP(1),IHE(1),PIV(1),T(1)
  84.       DOUBLE PRECISION DSUM
  85. C
  86. C        TEST ON WRONG INPUT PARAMETERS N AND M
  87.       IER=-1
  88.       IF (N-1) 81,81,1
  89.     1 IF(M) 81,81,2
  90. C
  91. C        INITIALIZE CHARACTERISTIC VECTORS FOR THE TABLEAU
  92.     2 IER=0
  93. C
  94. C        PREPARE TOP-ROW TOP
  95.       DO 3 I=1,N
  96.       K=I+N
  97.       J=K+N
  98.       TOP(J)=TOP(K)
  99.     3 TOP(K)=-TOP(I)
  100. C
  101. C        PREPARE INVERSE TRANSFORMATION MATRIX T
  102.       L=M+2
  103.       LL=L*L
  104.       DO 4 I=1,LL
  105.     4 T(I)=0.
  106.       K=1
  107.       J=L+1
  108.       DO 5 I=1,L
  109.       T(K)=1.
  110.     5 K=K+J
  111. C
  112. C        PREPARE INDEX-VECTOR IHE
  113.       DO 6 I=1,L
  114.       K=I+L
  115.       J=K+L
  116.       IHE(I)=0
  117.       IHE(K)=I
  118.     6 IHE(J)=1-I
  119.       NAN=N+N
  120.       K=L+L+L
  121.       J=K+NAN
  122.       DO 7 I=1,NAN
  123.       K=K+1
  124.       IHE(K)=I
  125.       J=J+1
  126.     7 IHE(J)=I
  127. C
  128. C        SET COUNTER ITER FOR ITERATION-STEPS
  129.       ITER=-1
  130.     8 ITER=ITER+1
  131. C
  132. C        TEST FOR MAXIMUM ITERATION-STEPS
  133.       IF(N+M-ITER) 9,9,10
  134.     9 IER=1
  135.       GO TO 69
  136. C
  137. C        DETERMINE THE COLUMN WITH THE MOST POSITIVE ELEMENT IN TOP
  138.    10 ISE=0
  139.       IPIV=0
  140.       K=L+L+L
  141.       SAVE=0.
  142. C
  143. C        START TOP-LOOP
  144.       DO 14 I=1,NAN
  145.       IDO=K+I
  146.       HELP=TOP(I)
  147.       IF(HELP-SAVE) 12,12,11
  148.    11 SAVE=HELP
  149.       IPIV=I
  150.    12 IF(IHE(IDO)) 14,13,14
  151.    13 ISE=I
  152.    14 CONTINUE
  153. C        END OF TOP-LOOP
  154. C
  155. C        IS OPTIMAL TABLEAU REACHED
  156.       IF(IPIV) 69,69,15
  157. C
  158. C        DETERMINE THE PIVOT-ELEMENT FOR THE COLUMN CHOSEN UPOVE
  159.    15 ILAB=1
  160.       IND=0
  161.       J=ISE
  162.       IF(J) 21,21,34
  163. C
  164. C        TRANSFER K-TH COLUMN FROM T TO PIV
  165.    16 K=(K-1)*L
  166.       DO 17 I=1,L
  167.       J=L+I
  168.       K=K+1
  169.    17 PIV(J)=T(K)
  170. C
  171. C        IS ANOTHER COLUMN NEEDED FOR SEARCH FOR PIVOT-ELEMENT
  172.    18 IF(ISE) 22,22,19
  173.    19 ISE=-ISE
  174. C
  175. C        TRANSFER COLUMNS IN PIV
  176.       J=L+1
  177.       IDO=L+L
  178.       DO 20 I=J,IDO
  179.       K=I+L
  180.    20 PIV(K)=PIV(I)
  181.    21 J=IPIV
  182.       GO TO 34
  183. C
  184. C        SEARCH PIVOT-ELEMENT PIV(IND)
  185.    22 SAVE=1.E38
  186.       IDO=0
  187.       K=L+1
  188.       LL=L+L
  189.       IND=0
  190. C
  191. C        START PIVOT-LOOP
  192.       DO 29 I=K,LL
  193.       J=I+L
  194.       HELP=PIV(I)
  195.       IF(HELP) 29,29,23
  196.    23 HELP=-HELP
  197.       IF(ISE) 26,24,26
  198.    24 IF(IHE(J)) 27,25,27
  199.    25 IDO=I
  200.       GO TO 29
  201.    26 HELP=-PIV(J)/HELP
  202.    27 IF(HELP-SAVE) 28,29,29
  203.    28 SAVE=HELP
  204.       IND=I
  205.    29 CONTINUE
  206. C        END OF PIVOT-LOOP
  207. C
  208. C        TEST FOR SUITABLE PIVOT-ELEMENT
  209.       IF(IND) 30,30,32
  210.    30 IF(IDO) 68,68,31
  211.    31 IND=IDO
  212. C        PIVOT-ELEMENT IS STORED IN PIV(IND)
  213. C
  214. C        COMPUTE THE RECIPROCAL OF THE PIVOT-ELEMENT REPI
  215.    32 REPI=1./PIV(IND)
  216.       IND=IND-L
  217. C
  218. C        UPDATE THE TOP-ROW TOP OF THE TABLEAU
  219.       ILAB=0
  220.       SAVE=-TOP(IPIV)*REPI
  221.       TOP(IPIV)=SAVE
  222. C
  223. C        INITIALIZE J AS COUNTER FOR TOP-LOOP
  224.       J=NAN
  225.    33 IF(J-IPIV) 34,53,34
  226.    34 K=0
  227. C
  228. C        SEARCH COLUMN IN TRANSFORMATION-MATRIX T
  229.       DO 36 I=1,L
  230.       IF(IHE(I)-J) 36,35,36
  231.    35 K=I
  232.       IF(ILAB) 50,50,16
  233.    36 CONTINUE
  234. C
  235. C        GENERATE COLUMN USING SUBROUTINE FCT AND TRANSFORMATION-MATRIX
  236.       I=L+L+L+NAN+J
  237.       I=IHE(I)-N
  238.       IF(I) 37,37,38
  239.    37 I=I+N
  240.       K=1
  241.    38 I=I+NAN
  242. C
  243. C        CALL SUBROUTINE FCT
  244.       CALL FCT(PIV,TOP(I),M-1)
  245. C
  246. C        PREPARE THE CALLED VECTOR PIV
  247.       DSUM=0.D0
  248.       IDO=M
  249.       DO 41 I=1,M
  250.       HELP=PIV(IDO)
  251.       IF(K) 39,39,40
  252.    39 HELP=-HELP
  253.    40 DSUM=DSUM+DBLE(HELP)
  254.       PIV(IDO+1)=HELP
  255.    41 IDO=IDO-1
  256.       PIV(L)=-DSUM
  257.       PIV(1)=1.
  258. C
  259. C        TRANSFORM VECTOR PIV WITH ROWS OF MATRIX T
  260.       IDO=IND
  261.       IF(ILAB) 44,44,42
  262.    42 K=1
  263.    43 IDO=K
  264.    44 DSUM=0.D0
  265.       HELP=0.
  266. C
  267. C        START MULTIPLICATION-LOOP
  268.       DO 46 I=1,L
  269.       DSUM=DSUM+DBLE(PIV(I)*T(IDO))
  270.       TOL=ABS(SNGL(DSUM))
  271.       IF(TOL-HELP) 46,46,45
  272.    45 HELP=TOL
  273.    46 IDO=IDO+L
  274. C        END OF MULTIPLICATION-LOOP
  275. C
  276.       TOL=1.E-5*HELP
  277.       IF(ABS(SNGL(DSUM))-TOL) 47,47,48
  278.    47 DSUM=0.D0
  279.    48 IF(ILAB) 51,51,49
  280.    49 I=K+L
  281.       PIV(I)=DSUM
  282. C
  283. C        TEST FOR LAST COLUMN-TERM
  284.       K=K+1
  285.       IF(K-L) 43,43,18
  286.    50 I=(K-1)*L+IND
  287.       DSUM=T(I)
  288. C
  289. C        COMPUTE NEW TOP-ELEMENT
  290.    51 DSUM=DSUM*DBLE(SAVE)
  291.       TOL=1.E-5*ABS(SNGL(DSUM))
  292.       TOP(J)=TOP(J)+SNGL(DSUM)
  293.       IF(ABS(TOP(J))-TOL) 52,52,53
  294.    52 TOP(J)=0.
  295. C
  296. C        TEST FOR LAST TOP-TERM
  297.    53 J=J-1
  298.       IF(J) 54,54,33
  299. C        END OF TOP-LOOP
  300. C
  301. C        TRANSFORM PIVOT-COLUMN
  302.    54 I=IND+L
  303.       PIV(I)=-1.
  304.       DO 55 I=1,L
  305.       J=I+L
  306.    55 PIV(I)=-PIV(J)*REPI
  307. C
  308. C        UPDATE TRANSFORMATION-MATRIX T
  309.       J=0
  310.       DO 57 I=1,L
  311.       IDO=J+IND
  312.       SAVE=T(IDO)
  313.       T(IDO)=0.
  314.       DO 56 K=1,L
  315.       ISE=K+J
  316.    56 T(ISE)=T(ISE)+SAVE*PIV(K)
  317.    57 J=J+L
  318. C
  319. C        UPDATE INDEX-VECTOR IHE
  320. C        INITIALIZE CHARACTERISTICS
  321.       J=0
  322.       K=0
  323.       ISE=0
  324.       IDO=0
  325. C
  326. C        START QUESTION-LOOP
  327.       DO 61 I=1,L
  328.       LL=I+L
  329.       ILAB=IHE(LL)
  330.       IF(IHE(I)-IPIV) 59,58,59
  331.    58 ISE=I
  332.       J=ILAB
  333.    59 IF(ILAB-IND) 61,60,61
  334.    60 IDO=I
  335.       K=IHE(I)
  336.    61 CONTINUE
  337. C        END OF QUESTION-LOOP
  338. C
  339. C        START MODIFICATION
  340.       IF(K) 62,62,63
  341.    62 IHE(IDO)=IPIV
  342.       IF(ISE) 67,67,65
  343.    63 IF(IND-J) 64,66,64
  344.    64 LL=L+L+L+NAN
  345.       K=K+LL
  346.       I=IPIV+LL
  347.       ILAB=IHE(K)
  348.       IHE(K)=IHE(I)
  349.       IHE(I)=ILAB
  350.       IF(ISE) 67,67,65
  351.    65 IDO=IDO+L
  352.       I=ISE+L
  353.       IHE(IDO)=J
  354.       IHE(I)=IND
  355.    66 IHE(ISE)=0
  356.    67 LL=L+L
  357.       J=LL+IND
  358.       I=LL+L+IPIV
  359.       ILAB=IHE(I)
  360.       IHE(I)=IHE(J)
  361.       IHE(J)=ILAB
  362. C        END OF MODIFICATION
  363. C
  364.       GO TO 8
  365. C
  366. C        SET ERROR PARAMETER IER=-1 SINCE NO SUITABLE PIVOT IS FOUND
  367.    68 IER=-1
  368. C
  369. C        EVALUATE FINAL TABLEAU
  370. C        COMPUTE SAVE AS MAXIMUM ERROR OF APPROXIMATION AND
  371. C        HELP AS ADDITIVE CONSTANCE FOR RESULTING COEFFICIENTS
  372.    69 SAVE=0.
  373.       HELP=0.
  374.       K=L+L+L
  375.       DO 73 I=1,NAN
  376.       IDO=K+I
  377.       J=IHE(IDO)
  378.       IF(J) 71,70,73
  379.    70 SAVE=-TOP(I)
  380.    71 IF(M+J+1) 73,72,73
  381.    72 HELP=TOP(I)
  382.    73 CONTINUE
  383. C
  384. C        PREPARE T,TOP,PIV
  385.       T(1)=SAVE
  386.       IDO=NAN+1
  387.       J=NAN+N
  388.       DO 74 I=IDO,J
  389.    74 TOP(I)=SAVE
  390.       DO 75 I=1,M
  391.    75 PIV(I)=HELP
  392. C
  393. C        COMPUTE COEFFICIENTS OF RESULTING POLYNOMIAL IN PIV(1) UP TO PI
  394. C        AND CALCULATE ERRORS AT GIVEN NODES IN TOP(1) UP TO TOP(N)
  395.       DO 79 I=1,NAN
  396.       IDO=K+I
  397.       J=IHE(IDO)
  398.       IF(J) 76,79,77
  399.    76 J=-J
  400.       PIV(J)=HELP-TOP(I)
  401.       GO TO 79
  402.    77 IF(J-N) 78,78,79
  403.    78 J=J+NAN
  404.       TOP(J)=SAVE+TOP(I)
  405.    79 CONTINUE
  406.       DO 80 I=1,N
  407.       IDO=NAN+I
  408.    80 TOP(I)=TOP(IDO)
  409.    81 RETURN
  410.       END
  411.