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