home *** CD-ROM | disk | FTP | other *** search
/ Frostbyte's 1980s DOS Shareware Collection / floppyshareware.zip / floppyshareware / DOOG / PCSSP1.ZIP / ITRPAPSM.ZIP / DAPFS.FOR < prev    next >
Text File  |  1985-11-29  |  7KB  |  190 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE DAPFS
  5. C
  6. C        PURPOSE
  7. C           PERFORM SYMMETRIC FACTORIZATION OF THE MATRIX OF THE NORMAL
  8. C           EQUATIONS FOLLOWED BY CALCULATION OF THE LEAST SQUARES FIT
  9. C           OPTIONALLY
  10. C
  11. C        USAGE
  12. C           CALL DAPFS(WORK,IP,IRES,IOP,EPS,ETA,IER)
  13. C
  14. C        DESCRIPTION OF PARAMETERS
  15. C           WORK  - GIVEN SYMMETRIC COEFFICIENT MATRIX, STORED
  16. C                   COMPRESSED, I.E UPPER TRIANGULAR PART COLUMNWISE.
  17. C                   THE GIVEN RIGHT HAND SIDE OCCUPIES THE NEXT IP
  18. C                   LOCATIONS IN WORK. THE VERY LAST COMPONENT OF WORK
  19. C                   CONTAINS THE SQUARE SUM OF FUNCTION VALUES E0
  20. C                   THIS SCHEME OF STORAGE ALLOCATION IS PRODUCED E.G.
  21. C                   BY SUBROUTINE APLL.
  22. C                   THE GIVEN MATRIX IS FACTORED IN THE FORM
  23. C                   TRANSPOSE(T)*T AND THE GIVEN RIGHT HAND SIDE IS
  24. C                   DIVIDED BY TRANSPOSE(T).
  25. C                   THE UPPER TRIANGULAR FACTOR T IS RETURNED IN WORK IF
  26. C                   IOP EQUALS ZERO.
  27. C                   IN CASE OF NONZERO IOP THE CALCULATED SOLUTIONS ARE
  28. C                   STORED IN THE COLUMNS OF TRIANGULAR ARRAY WORK OF
  29. C                   CORRESPONDING DIMENSION AND E0  IS REPLACED BY THE
  30. C                   SQUARE SUM OF THE ERRORS FOR FIT OF DIMENSION IRES.
  31. C                   THE TOTAL DIMENSION OF WORK IS (IP+1)*(IP+2)/2
  32. C                   WORK MUST BE OF DOUBLE PRECISION
  33. C           IP    - NUMBER OF FUNDAMENTAL FUNCTIONS USED FOR LEAST
  34. C                   SQUARES FIT
  35. C           IRES  - DIMENSION OF CALCULATED LEAST SQUARES FIT.
  36. C                   LET N1, N2, DENOTE THE FOLLOWING NUMBERS
  37. C                   N1 = MAXIMAL DIMENSION FOR WHICH NO LOSS OF
  38. C                        SIGNIFICANCE WAS INDICATED DURING FACTORIZATION
  39. C                   N2 = SMALLEST DIMENSION FOR WHICH THE SQUARE SUM OF
  40. C                        THE ERRORS DOES NOT EXCEED TEST=ABS(ETA*FSQ)
  41. C                   THEN IRES=MINO(IP,N1) IF IOP IS NONNEGATIVE
  42. C                   AND  IRES=MINO(IP,N1,N2) IF IOP IS NEGATIVE
  43. C           IOP   - INPUT PARAMETER FOR SELECTION OF OPERATION
  44. C                   IOP = 0 MEANS TRIANGULAR FACTORIZATION, DIVISION OF
  45. C                           THE RIGHT HAND SIDE BY TRANSPOSE(T) AND
  46. C                           CALCULATION OF THE SQUARE SUM OF ERRORS IS
  47. C                           PERFORMED ONLY
  48. C                   IOP = +1 OR -1 MEANS THE SOLUTION OF DIMENSION IRES
  49. C                           IS CALCULATED ADDITIONALLY
  50. C                   IOP = +2 OR -2 MEANS ALL SOLUTIONS FOR DIMENSION ONE
  51. C                           UP TO IRES ARE CALCULATED ADDITIONALLY
  52. C           EPS   - RELATIVE TOLERANCE FOR TEST ON LOSS OF SIGNIFICANCE.
  53. C                   A SENSIBLE VALUE IS BETWEEN 1.E-10 AND 1.E-15
  54. C           ETA   - RELATIVE TOLERANCE FOR TOLERATED SQUARE SUM OF
  55. C                   ERRORS. A REALISTIC VALUE IS BETWEEN 1.E0 AND 1.E-15
  56. C           IER   - RESULTANT ERROR PARAMETER
  57. C                   IER =-1 MEANS NONPOSITIVE IP
  58. C                   IER = 0 MEANS NO LOSS OF SIGNIFICANCE DETECTED
  59. C                           AND SPECIFIED TOLERANCE OF ERRORS REACHED
  60. C                   IER = 1 MEANS LOSS OF SIGNIFICANCE DETECTED OR
  61. C                           SPECIFIED TOLERANCE OF ERRORS NOT REACHED
  62. C
  63. C        REMARKS
  64. C           THE ABSOLUTE TOLERANCE USED INTERNALLY FOR TEST ON LOSS OF
  65. C           SIGNIFICANCE IS TOL=ABS(EPS*SNGL(WORK(1))).
  66. C           THE ABSOLUTE TOLERANCE USED INTERNALLY FOR THE SQUARE SUM OF
  67. C           ERRORS IS ABS(ETA*SNGL(FSQ)).
  68. C           IOP GREATER THAN 2 HAS THE SAME EFFECT AS IOP = 2.
  69. C           IOP LESS THAN -2 HAS THE SAME EFFECT AS IOP =-2.
  70. C           IRES = 0 MEANS THE ABSOLUTE VALUE OF EPS IS NOT LESS THAN
  71. C           ONE AND/OR WORK(1) IS NOT POSITIVE AND/OR IP IS NOT POSITIVE
  72. C
  73. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  74. C           NONE
  75. C
  76. C        METHOD
  77. C           CALCULATION OF THE LEAST SQUARES FITS IS DONE USING
  78. C           CHOLESKYS SQUARE ROOT METHOD FOR SYMMETRIC FACTORIZATION.
  79. C           THE INCORPORATED TEST ON LOSS OF SIGNIFICANCE MEANS EACH
  80. C           RADICAND MUST BE GREATER THAN THE INTERNAL ABSOLUTE
  81. C           TOLERANCE TOL.
  82. C           IN CASE OF LOSS OF SIGNIFICANCE IN THE ABOVE SENSE ONLY A
  83. C           SUBSYSTEM OF THE NORMAL EQUATIONS IS SOLVED.
  84. C           IN CASE OF NEGATIVE IOP THE TRIANGULAR FACTORIZATION IS
  85. C           TERMINATED PREMATURELY EITHER IF THE SQUARE SUM OF THE
  86. C           ERRORS DOES NOT EXCEED ETA*FSQ OR IF THERE IS INDICATION
  87. C           FOR LOSS OF SIGNIFICANCE
  88. C
  89. C     ..................................................................
  90. C
  91.       SUBROUTINE DAPFS(WORK,IP,IRES,IOP,EPS,ETA,IER)
  92. C
  93. C
  94. C        DIMENSIONED DUMMY VARIABLES
  95.       DIMENSION WORK(1)
  96.       DOUBLE PRECISION WORK,SUM,PIV
  97.       IRES=0
  98. C
  99. C        TEST OF SPECIFIED DIMENSION
  100.       IF(IP)1,1,2
  101. C
  102. C        ERROR RETURN IN CASE OF ILLEGAL DIMENSION
  103.     1 IER=-1
  104.       RETURN
  105. C
  106. C        INITIALIZE FACTORIZATION PROCESS
  107.     2 IPIV=0
  108.       IPP1=IP+1
  109.       IER=1
  110.       ITE=IP*IPP1/2
  111.       IEND=ITE+IPP1
  112.       TOL=ABS(EPS*SNGL(WORK(1)))
  113.       TEST=ABS(ETA*SNGL(WORK(IEND)))
  114. C
  115. C        START LOOP OVER ALL ROWS OF WORK
  116.       DO 11 I=1,IP
  117.       IPIV=IPIV+I
  118.       JA=IPIV-IRES
  119.       JE=IPIV-1
  120. C
  121. C        FORM SCALAR PRODUCT NEEDED TO MODIFY CURRENT ROW ELEMENTS
  122.       JK=IPIV
  123.       DO 9 K=I,IPP1
  124.       SUM=0.D0
  125.       IF(IRES)5,5,3
  126.     3 JK=JK-IRES
  127.       DO 4 J=JA,JE
  128.       SUM=SUM+WORK(J)*WORK(JK)
  129.     4 JK=JK+1
  130.     5 IF(JK-IPIV)6,6,8
  131. C
  132. C        TEST FOR LOSS OF SIGNIFICANCE
  133.     6 SUM=WORK(IPIV)-SUM
  134.       IF(SNGL(SUM)-TOL)12,12,7
  135.     7 SUM=DSQRT(SUM)
  136.       WORK(IPIV)=SUM
  137.       PIV=1.D0/SUM
  138.       GOTO 9
  139. C
  140. C        UPDATE OFF-DIAGONAL TERMS
  141.     8 SUM=(WORK(JK)-SUM)*PIV
  142.       WORK(JK)=SUM
  143.     9 JK=JK+K
  144. C
  145. C        UPDATE SQUARE SUM OF ERRORS
  146.       WORK(IEND)=WORK(IEND)-SUM*SUM
  147. C
  148. C        RECORD ADDRESS OF LAST PIVOT ELEMENT
  149.       IRES=IRES+1
  150.       IADR=IPIV
  151. C
  152. C        TEST FOR TOLERABLE ERROR IF SPECIFIED
  153.       IF(IOP)10,11,11
  154.    10 IF(SNGL(WORK(IEND))-TEST)13,13,11
  155.    11 CONTINUE
  156.       IF(IOP)12,22,12
  157. C
  158. C        PERFORM BACK SUBSTITUTION IF SPECIFIED
  159.    12 IF(IOP)14,23,14
  160.    13 IER=0
  161.    14 IPIV=IRES
  162.    15 IF(IPIV)23,23,16
  163.    16 SUM=0.D0
  164.       JA=ITE+IPIV
  165.       JJ=IADR
  166.       JK=IADR
  167.       K=IPIV
  168.       DO 19 I=1,IPIV
  169.       WORK(JK)=(WORK(JA)-SUM)/WORK(JJ)
  170.       IF(K-1)20,20,17
  171.    17 JE=JJ-1
  172.       SUM=0.D0
  173.       DO 18 J=K,IPIV
  174.       SUM=SUM+WORK(JK)*WORK(JE)
  175.       JK=JK+1
  176.    18 JE=JE+J
  177.       JK=JE-IPIV
  178.       JA=JA-1
  179.       JJ=JJ-K
  180.    19 K=K-1
  181.    20 IF(IOP/2)21,23,21
  182.    21 IADR=IADR-IPIV
  183.       IPIV=IPIV-1
  184.       GOTO 15
  185. C
  186. C        NORMAL RETURN
  187.    22 IER=0
  188.    23 RETURN
  189.       END
  190.