home *** CD-ROM | disk | FTP | other *** search
/ Frostbyte's 1980s DOS Shareware Collection / floppyshareware.zip / floppyshareware / DOOG / PCSSP1.ZIP / ITRPAPSM.ZIP / DAPLL.FOR < prev    next >
Text File  |  1985-11-29  |  5KB  |  125 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE DAPLL
  5. C        PURPOSE
  6. C           SET UP NORMAL EQUATIONS FOR A LINEAR LEAST SQUARES FIT
  7. C           TO A GIVEN DISCRETE FUNCTION
  8. C
  9. C        USAGE
  10. C           CALL DAPLL(FFCT,N,IP,P,WORK,DATI,IER)
  11. C           SUBROUTINE FFCT REQUIRES AN EXTERNAL STATEMENT
  12. C
  13. C        DESCRIPTION OF PARAMETERS
  14. C           FFCT  - USER CODED SUBROUTINE WHICH MUST BE DECLARED
  15. C                   EXTERNAL IN THE MAIN PROGRAM. IT IS CALLED
  16. C                   CALL FFCT(I,N,IP,P,DATI,WGT,IER) AND RETURNS
  17. C                   THE VALUES OF THE FUNDAMENTAL FUNCTIONS FOR
  18. C                   THE I-TH ARGUMENT IN P(1) UP TO P(IP)
  19. C                   FOLLOWED BY THE I-TH FUNCTION VALUE IN P(IP+1)
  20. C                   N IS THE NUMBER OF ALL POINTS
  21. C                   P,DATI,WGT MUST BE OF DOUBLE PRECISION.
  22. C                   DATI IS A DUMMY PARAMETER WHICH IS USED AS ARRAY
  23. C                   NAME. THE GIVEN DATA SET MAY BE ALLOCATED IN DATI
  24. C                   WGT IS THE WEIGHT FACTOR FOR THE I-TH POINT
  25. C                   IER IS USED AS RESULTANT ERROR PARAMETER IN FFCT
  26. C           N     - NUMBER OF GIVEN POINTS
  27. C           IP    - NUMBER OF FUNDAMENTAL FUNCTIONS USED FOR LEAST
  28. C                   SQUARES FIT
  29. C                   IP SHOULD NOT EXCEED N
  30. C           P     - WORKING STORAGE OF DIMENSION IP+1, WHICH
  31. C                   IS USED AS INTERFACE BETWEEN APLL AND THE USER
  32. C                   CODED SUBROUTINE FFCT
  33. C                   P MUST BE OF DOUBLE PRECISION.
  34. C           WORK  - WORKING STORAGE OF DIMENSION (IP+1)*(IP+2)/2.
  35. C                   ON RETURN WORK CONTAINS THE SYMMETRIC COEFFICIENT
  36. C                   MATRIX OF THE NORMAL EQUATIONS IN COMPRESSED FORM,
  37. C                   I.E. UPPER TRINGULAR PART ONLY STORED COLUMNWISE.
  38. C                   THE FOLLOWING IP POSITIONS CONTAIN THE RIGHT
  39. C                   HAND SIDE AND WORK((IP+1)*(IP+2)/2) CONTAINS
  40. C                   THE WEIGHTED SQUARE SUM OF THE FUNCTION VALUES
  41. C                   WORK MUST BE OF DOUBLE PRECISION.
  42. C           DATI  - DUMMY ENTRY TO COMMUNICATE AN ARRAY NAME BETWEEN
  43. C                   MAIN LINE AND SUBROUTINE FFCT.
  44. C                   DATI MUST BE OF DOUBLE PRECISION.
  45. C           IER   - RESULTING ERROR PARAMETER
  46. C                   IER =-1 MEANS FORMAL ERRORS IN SPECIFIED DIMENSIONS
  47. C                   IER = 0 MEANS NO ERRORS
  48. C                   IER = 1 MEANS ERROR IN EXTERNAL SUBROUTINE FFCT
  49. C
  50. C        REMARKS
  51. C           TO ALLOW FOR EASY COMMUNICATION OF INTEGER VALUES
  52. C           BETWEEN MAINLINE AND EXTERNAL SUBROUTINE FFCT, THE ERROR
  53. C           PARAMETER IER IS TREATED AS A VECTOR OF DIMENSION 1 WITHIN
  54. C           SUBROUTINE DAPLL. ADDITIONAL COMPONENTS OF IER MAY BE
  55. C           INTRODUCED BY THE USER FOR COMMUNICATION BACK AND FORTH.
  56. C           IN THIS CASE, HOWEVER, THE USER MUST SPECIFY IER AS A
  57. C           VECTOR IN HIS MAINLINE.
  58. C           EXECUTION OF SUBROUTINE DAPLL IS A PREPARATORY STEP FOR
  59. C           CALCULATION OF THE LINEAR LEAST SQUARES FIT.
  60. C           NORMALLY IT IS FOLLOWED BY EXECUTION OF SUBROUTINE DAPFS
  61. C
  62. C       SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  63. C           THE EXTERNAL SUBROUTINE FFCT MUST BE FURNISHED BY THE USER
  64. C
  65. C        METHOD
  66. C           HANDLING OF THE GIVEN DATA SET (ARGUMENTS,FUNCTION VALUES
  67. C           AND WEIGHTS) IS COMPLETELY LEFT TO THE USER
  68. C           ESSENTIALLY HE HAS THREE CHOICES
  69. C           (1) THE I-TH VALUES OF ARGUMENT, FUNCTION VALUE AND WEIGHT
  70. C               ARE CALCULATED WITHIN SUBROUTINE FFCT FOR GIVEN I.
  71. C           (2) THE I-TH VALUES OF ARGUMENT, FUNCTION VALUE AND WEIGHT
  72. C               ARE DETERMINED BY TABLE LOOK UP. THE STORAGE LOCATIONS
  73. C               REQUIRED ARE ALLOCATED WITHIN THE DUMMY ARRAY DATI
  74. C               (POSSIBLY IN P TOO, IN EXCESS OF THE SPECIFIED IP + 1
  75. C               LOCATIONS).
  76. C               ANOTHER POSSIBILITY WOULD BE TO USE COMMON AS INTERFACE
  77. C               BETWEEN MAIN LINE AND SUBROUTINE FFCT AND TO ALLOCATE
  78. C               STORAGE FOR THE DATA SET IN COMMON.
  79. C           (3) THE I-TH VALUES OF ARGUMENT, FUNCTION VALUE AND WEIGHT
  80. C               ARE READ IN FROM AN EXTERNAL DEVICE. THIS MAY BE EASILY
  81. C               ACCOMPLISHED SINCE I IS USED STRICTLY INCREASING FROM
  82. C               ONE UP TO N WITHIN APLL
  83. C
  84. C     ..................................................................
  85. C
  86.       SUBROUTINE DAPLL(FFCT,N,IP,P,WORK,DATI,IER)
  87. C
  88. C
  89. C        DIMENSIONED DUMMY VARIABLES
  90.       DIMENSION P(1),WORK(1),DATI(1),IER(1)
  91.       DOUBLE PRECISION P,WORK,DATI,WGT,AUX
  92. C
  93. C        CHECK FOR FORMAL ERRORS IN SPECIFIED DIMENSIONS
  94.       IF(N)10,10,1
  95.     1 IF(IP)10,10,2
  96.     2 IF(N-IP)10,3,3
  97. C
  98. C        SET WORKING STORAGE AND RIGHT HAND SIDE TO ZERO
  99.     3 IPP1=IP+1
  100.       M=IPP1*(IP+2)/2
  101.       IER(1)=0
  102.       DO 4 I=1,M
  103.     4 WORK(I)=0.D0
  104. C
  105. C        START GREAT LOOP OVER ALL GIVEN POINTS
  106.       DO 8 I=1,N
  107.       CALL FFCT(I,N,IP,P,DATI,WGT,IER)
  108.       IF(IER(1))9,5,9
  109.     5 J=0
  110.       DO 7 K=1,IPP1
  111.       AUX=P(K)*WGT
  112.       DO 6 L=1,K
  113.       J=J+1
  114.     6 WORK(J)=WORK(J)+P(L)*AUX
  115.     7 CONTINUE
  116.     8 CONTINUE
  117. C
  118. C        NORMAL RETURN
  119.     9 RETURN
  120. C
  121. C        ERROR RETURN IN CASE OF FORMAL ERRORS
  122.    10 IER(1)=-1
  123.       RETURN
  124.       END
  125.