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