home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / INFO / FORTRAN / LLSQ.ZIP / HFTI.FOR < prev    next >
Encoding:
Text File  |  1984-02-23  |  10.9 KB  |  137 lines
  1.       SUBROUTINE HFTI (A,MDA,M,N,B,MDB,NB,TAU,KRANK,RNORM,H,G,IP)       HFT00100
  2. C     C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY, 1973 JUN 12 HFT00200
  3. C     TO APPEAR IN 'SOLVING LEAST SQUARES PROBLEMS', PRENTICE-HALL, 1974HFT00300
  4. C          SOLVE LEAST SQUARES PROBLEM USING ALGORITHM, HFTI.           HFT00400
  5. C                                                                       HFT00500
  6.       DIMENSION A(MDA,N),B(MDB,NB),H(N),G(N),RNORM(NB)                  HFT00600
  7.       INTEGER IP(N)                                                     HFT00700
  8.       DOUBLE PRECISION SM,DZERO                                         HFT00800
  9.       SZERO=0.                                                          HFT00900
  10.       DZERO=0.D0                                                        HFT01000
  11.       FACTOR=0.001                                                      HFT01100
  12. C                                                                       HFT01200
  13.       K=0                                                               HFT01300
  14.       LDIAG=MIN0(M,N)                                                   HFT01400
  15.       IF (LDIAG.LE.0) GO TO 270                                         HFT01500
  16.           DO 80 J=1,LDIAG                                               HFT01600
  17.           IF (J.EQ.1) GO TO 20                                          HFT01700
  18. C                                                                       HFT01800
  19. C     UPDATE SQUARED COLUMN LENGTHS AND FIND LMAX                       HFT01900
  20. C    ..                                                                 HFT02000
  21.           LMAX=J                                                        HFT02100
  22.               DO 10 L=J,N                                               HFT02200
  23.               H(L)=H(L)-A(J-1,L)**2                                     HFT02300
  24.               IF (H(L).GT.H(LMAX)) LMAX=L                               HFT02400
  25.    10         CONTINUE                                                  HFT02500
  26.           IF(DIFF(HMAX+FACTOR*H(LMAX),HMAX)) 20,20,50                   HFT02600
  27. C                                                                       HFT02700
  28. C     COMPUTE SQUARED COLUMN LENGTHS AND FIND LMAX                      HFT02800
  29. C    ..                                                                 HFT02900
  30.    20     LMAX=J                                                        HFT03000
  31.               DO 40 L=J,N                                               HFT03100
  32.               H(L)=0.                                                   HFT03200
  33.                   DO 30 I=J,M                                           HFT03300
  34.    30             H(L)=H(L)+A(I,L)**2                                   HFT03400
  35.               IF (H(L).GT.H(LMAX)) LMAX=L                               HFT03500
  36.    40         CONTINUE                                                  HFT03600
  37.           HMAX=H(LMAX)                                                  HFT03700
  38. C    ..                                                                 HFT03800
  39. C     LMAX HAS BEEN DETERMINED                                          HFT03900
  40. C                                                                       HFT04000
  41. C     DO COLUMN INTERCHANGES IF NEEDED.                                 HFT04100
  42. C    ..                                                                 HFT04200
  43.    50     CONTINUE                                                      HFT04300
  44.           IP(J)=LMAX                                                    HFT04400
  45.           IF (IP(J).EQ.J) GO TO 70                                      HFT04500
  46.               DO 60 I=1,M                                               HFT04600
  47.               TMP=A(I,J)                                                HFT04700
  48.               A(I,J)=A(I,LMAX)                                          HFT04800
  49.    60         A(I,LMAX)=TMP                                             HFT04900
  50.           H(LMAX)=H(J)                                                  HFT05000
  51. C                                                                       HFT05100
  52. C     COMPUTE THE J-TH TRANSFORMATION AND APPLY IT TO A AND B.          HFT05200
  53. C    ..                                                                 HFT05300
  54.    70     CALL H12 (1,J,J+1,M,A(1,J),1,H(J),A(1,J+1),1,MDA,N-J)         HFT05400
  55.    80     CALL H12 (2,J,J+1,M,A(1,J),1,H(J),B,1,MDB,NB)                 HFT05500
  56. C                                                                       HFT05600
  57. C     DETERMINE THE PSEUDORANK, K, USING THE TOLERANCE, TAU.            HFT05700
  58. C    ..                                                                 HFT05800
  59.           DO 90 J=1,LDIAG                                               HFT05900
  60.           IF (ABS(A(J,J)).LE.TAU) GO TO 100                             HFT06000
  61.    90     CONTINUE                                                      HFT06100
  62.       K=LDIAG                                                           HFT06200
  63.       GO TO 110                                                         HFT06300
  64.   100 K=J-1                                                             HFT06400
  65.   110 KP1=K+1                                                           HFT06500
  66. C                                                                       HFT06600
  67. C     COMPUTE THE NORMS OF THE RESIDUAL VECTORS.                        HFT06700
  68. C                                                                       HFT06800
  69.       IF (NB.LE.0) GO TO 140                                            HFT06900
  70.           DO 130 JB=1,NB                                                HFT07000
  71.           TMP=SZERO                                                     HFT07100
  72.           IF (KP1.GT.M) GO TO 130                                       HFT07200
  73.               DO 120 I=KP1,M                                            HFT07300
  74.   120         TMP=TMP+B(I,JB)**2                                        HFT07400
  75.   130     RNORM(JB)=SQRT(TMP)                                           HFT07500
  76.   140 CONTINUE                                                          HFT07600
  77. C                                           SPECIAL FOR PSEUDORANK = 0  HFT07700
  78.       IF (K.GT.0) GO TO 160                                             HFT07800
  79.       IF (NB.LE.0) GO TO 270                                            HFT07900
  80.           DO 150 JB=1,NB                                                HFT08000
  81.               DO 150 I=1,N                                              HFT08100
  82.   150         B(I,JB)=SZERO                                             HFT08200
  83.       GO TO 270                                                         HFT08300
  84. C                                                                       HFT08400
  85. C     IF THE PSEUDORANK IS LESS THAN N COMPUTE HOUSEHOLDER              HFT08500
  86. C     DECOMPOSITION OF FIRST K ROWS.                                    HFT08600
  87. C    ..                                                                 HFT08700
  88.   160 IF (K.EQ.N) GO TO 180                                             HFT08800
  89.           DO 170 II=1,K                                                 HFT08900
  90.           I=KP1-II                                                      HFT09000
  91.   170     CALL H12 (1,I,KP1,N,A(I,1),MDA,G(I),A,MDA,1,I-1)              HFT09100
  92.   180 CONTINUE                                                          HFT09200
  93. C                                                                       HFT09300
  94. C                                                                       HFT09400
  95.       IF (NB.LE.0) GO TO 270                                            HFT09500
  96.           DO 260 JB=1,NB                                                HFT09600
  97. C                                                                       HFT09700
  98. C     SOLVE THE K BY K TRIANGULAR SYSTEM.                               HFT09800
  99. C    ..                                                                 HFT09900
  100.               DO 210 L=1,K                                              HFT10000
  101.               SM=DZERO                                                  HFT10100
  102.               I=KP1-L                                                   HFT10200
  103.               IF (I.EQ.K) GO TO 200                                     HFT10300
  104.               IP1=I+1                                                   HFT10400
  105.                   DO 190 J=IP1,K                                        HFT10500
  106.   190             SM=SM+A(I,J)*DBLE(B(J,JB))                            HFT10600
  107.   200         SM1=SM                                                    HFT10700
  108.   210         B(I,JB)=(B(I,JB)-SM1)/A(I,I)                              HFT10800
  109. C                                                                       HFT10900
  110. C     COMPLETE COMPUTATION OF SOLUTION VECTOR.                          HFT11000
  111. C    ..                                                                 HFT11100
  112.           IF (K.EQ.N) GO TO 240                                         HFT11200
  113.               DO 220 J=KP1,N                                            HFT11300
  114.   220         B(J,JB)=SZERO                                             HFT11400
  115.               DO 230 I=1,K                                              HFT11500
  116.   230         CALL H12 (2,I,KP1,N,A(I,1),MDA,G(I),B(1,JB),1,MDB,1)      HFT11600
  117. C                                                                       HFT11700
  118. C      RE-ORDER THE SOLUTION VECTOR TO COMPENSATE FOR THE               HFT11800
  119. C      COLUMN INTERCHANGES.                                             HFT11900
  120. C    ..                                                                 HFT12000
  121.   240         DO 250 JJ=1,LDIAG                                         HFT12100
  122.               J=LDIAG+1-JJ                                              HFT12200
  123.               IF (IP(J).EQ.J) GO TO 250                                 HFT12300
  124.               L=IP(J)                                                   HFT12400
  125.               TMP=B(L,JB)                                               HFT12500
  126.               B(L,JB)=B(J,JB)                                           HFT12600
  127.               B(J,JB)=TMP                                               HFT12700
  128.   250         CONTINUE                                                  HFT12800
  129.   260     CONTINUE                                                      HFT12900
  130. C    ..                                                                 HFT13000
  131. C     THE SOLUTION VECTORS, X, ARE NOW                                  HFT13100
  132. C     IN THE FIRST  N  ROWS OF THE ARRAY B(,).                          HFT13200
  133. C                                                                       HFT13300
  134.   270 KRANK=K                                                           HFT13400
  135.       RETURN                                                            HFT13500
  136.       END                                                               HFT13600
  137.