home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / octave-1.1.1p1-src.tgz / tar.out / fsf / octave / libcruft / minpack / r1updt.f < prev   
Text File  |  1996-09-28  |  6KB  |  208 lines
  1.       SUBROUTINE R1UPDT(M,N,S,LS,U,V,W,SING)
  2.       INTEGER M,N,LS
  3.       LOGICAL SING
  4.       DOUBLE PRECISION S(LS),U(M),V(N),W(M)
  5. C     **********
  6. C
  7. C     SUBROUTINE R1UPDT
  8. C
  9. C     GIVEN AN M BY N LOWER TRAPEZOIDAL MATRIX S, AN M-VECTOR U,
  10. C     AND AN N-VECTOR V, THE PROBLEM IS TO DETERMINE AN
  11. C     ORTHOGONAL MATRIX Q SUCH THAT
  12. C
  13. C                   T
  14. C           (S + U*V )*Q
  15. C
  16. C     IS AGAIN LOWER TRAPEZOIDAL.
  17. C
  18. C     THIS SUBROUTINE DETERMINES Q AS THE PRODUCT OF 2*(N - 1)
  19. C     TRANSFORMATIONS
  20. C
  21. C           GV(N-1)*...*GV(1)*GW(1)*...*GW(N-1)
  22. C
  23. C     WHERE GV(I), GW(I) ARE GIVENS ROTATIONS IN THE (I,N) PLANE
  24. C     WHICH ELIMINATE ELEMENTS IN THE I-TH AND N-TH PLANES,
  25. C     RESPECTIVELY. Q ITSELF IS NOT ACCUMULATED, RATHER THE
  26. C     INFORMATION TO RECOVER THE GV, GW ROTATIONS IS RETURNED.
  27. C
  28. C     THE SUBROUTINE STATEMENT IS
  29. C
  30. C       SUBROUTINE R1UPDT(M,N,S,LS,U,V,W,SING)
  31. C
  32. C     WHERE
  33. C
  34. C       M IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
  35. C         OF ROWS OF S.
  36. C
  37. C       N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
  38. C         OF COLUMNS OF S. N MUST NOT EXCEED M.
  39. C
  40. C       S IS AN ARRAY OF LENGTH LS. ON INPUT S MUST CONTAIN THE LOWER
  41. C         TRAPEZOIDAL MATRIX S STORED BY COLUMNS. ON OUTPUT S CONTAINS
  42. C         THE LOWER TRAPEZOIDAL MATRIX PRODUCED AS DESCRIBED ABOVE.
  43. C
  44. C       LS IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN
  45. C         (N*(2*M-N+1))/2.
  46. C
  47. C       U IS AN INPUT ARRAY OF LENGTH M WHICH MUST CONTAIN THE
  48. C         VECTOR U.
  49. C
  50. C       V IS AN ARRAY OF LENGTH N. ON INPUT V MUST CONTAIN THE VECTOR
  51. C         V. ON OUTPUT V(I) CONTAINS THE INFORMATION NECESSARY TO
  52. C         RECOVER THE GIVENS ROTATION GV(I) DESCRIBED ABOVE.
  53. C
  54. C       W IS AN OUTPUT ARRAY OF LENGTH M. W(I) CONTAINS INFORMATION
  55. C         NECESSARY TO RECOVER THE GIVENS ROTATION GW(I) DESCRIBED
  56. C         ABOVE.
  57. C
  58. C       SING IS A LOGICAL OUTPUT VARIABLE. SING IS SET TRUE IF ANY
  59. C         OF THE DIAGONAL ELEMENTS OF THE OUTPUT S ARE ZERO. OTHERWISE
  60. C         V IS SET FALSE.
  61. C
  62. C     SUBPROGRAMS CALLED
  63. C
  64. C       MINPACK-SUPPLIED ... DPMPAR
  65. C
  66. C       FORTRAN-SUPPLIED ... DABS,DSQRT
  67. C
  68. C     MINPACK. VERSION OF DECEMBER 1978.
  69. C     BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE,
  70. C     JOHN L. NAZARETH
  71. C
  72. C     **********
  73.       INTEGER I,J,JJ,L,NMJ,NM1
  74.       DOUBLE PRECISION COS,COTAN,GIANT,ONE,P5,P25,SIN,TAN,TAU,TEMP,
  75.      *                 ZERO
  76.       DOUBLE PRECISION DPMPAR
  77.       DATA ONE,P5,P25,ZERO /1.0D0,5.0D-1,2.5D-1,0.0D0/
  78. C
  79. C     GIANT IS THE LARGEST MAGNITUDE.
  80. C
  81.       GIANT = DPMPAR(3)
  82. C
  83. C     INITIALIZE THE DIAGONAL ELEMENT POINTER.
  84. C
  85.       JJ = (N*(2*M - N + 1))/2 - (M - N)
  86. C
  87. C     MOVE THE NONTRIVIAL PART OF THE LAST COLUMN OF S INTO W.
  88. C
  89.       L = JJ
  90.       DO 10 I = N, M
  91.          W(I) = S(L)
  92.          L = L + 1
  93.    10    CONTINUE
  94. C
  95. C     ROTATE THE VECTOR V INTO A MULTIPLE OF THE N-TH UNIT VECTOR
  96. C     IN SUCH A WAY THAT A SPIKE IS INTRODUCED INTO W.
  97. C
  98.       NM1 = N - 1
  99.       IF (NM1 .LT. 1) GO TO 70
  100.       DO 60 NMJ = 1, NM1
  101.          J = N - NMJ
  102.          JJ = JJ - (M - J + 1)
  103.          W(J) = ZERO
  104.          IF (V(J) .EQ. ZERO) GO TO 50
  105. C
  106. C        DETERMINE A GIVENS ROTATION WHICH ELIMINATES THE
  107. C        J-TH ELEMENT OF V.
  108. C
  109.          IF (DABS(V(N)) .GE. DABS(V(J))) GO TO 20
  110.             COTAN = V(N)/V(J)
  111.             SIN = P5/DSQRT(P25+P25*COTAN**2)
  112.             COS = SIN*COTAN
  113.             TAU = ONE
  114.             IF (DABS(COS)*GIANT .GT. ONE) TAU = ONE/COS
  115.             GO TO 30
  116.    20    CONTINUE
  117.             TAN = V(J)/V(N)
  118.             COS = P5/DSQRT(P25+P25*TAN**2)
  119.             SIN = COS*TAN
  120.             TAU = SIN
  121.    30    CONTINUE
  122. C
  123. C        APPLY THE TRANSFORMATION TO V AND STORE THE INFORMATION
  124. C        NECESSARY TO RECOVER THE GIVENS ROTATION.
  125. C
  126.          V(N) = SIN*V(J) + COS*V(N)
  127.          V(J) = TAU
  128. C
  129. C        APPLY THE TRANSFORMATION TO S AND EXTEND THE SPIKE IN W.
  130. C
  131.          L = JJ
  132.          DO 40 I = J, M
  133.             TEMP = COS*S(L) - SIN*W(I)
  134.             W(I) = SIN*S(L) + COS*W(I)
  135.             S(L) = TEMP
  136.             L = L + 1
  137.    40       CONTINUE
  138.    50    CONTINUE
  139.    60    CONTINUE
  140.    70 CONTINUE
  141. C
  142. C     ADD THE SPIKE FROM THE RANK 1 UPDATE TO W.
  143. C
  144.       DO 80 I = 1, M
  145.          W(I) = W(I) + V(N)*U(I)
  146.    80    CONTINUE
  147. C
  148. C     ELIMINATE THE SPIKE.
  149. C
  150.       SING = .FALSE.
  151.       IF (NM1 .LT. 1) GO TO 140
  152.       DO 130 J = 1, NM1
  153.          IF (W(J) .EQ. ZERO) GO TO 120
  154. C
  155. C        DETERMINE A GIVENS ROTATION WHICH ELIMINATES THE
  156. C        J-TH ELEMENT OF THE SPIKE.
  157. C
  158.          IF (DABS(S(JJ)) .GE. DABS(W(J))) GO TO 90
  159.             COTAN = S(JJ)/W(J)
  160.             SIN = P5/DSQRT(P25+P25*COTAN**2)
  161.             COS = SIN*COTAN
  162.             TAU = ONE
  163.             IF (DABS(COS)*GIANT .GT. ONE) TAU = ONE/COS
  164.             GO TO 100
  165.    90    CONTINUE
  166.             TAN = W(J)/S(JJ)
  167.             COS = P5/DSQRT(P25+P25*TAN**2)
  168.             SIN = COS*TAN
  169.             TAU = SIN
  170.   100    CONTINUE
  171. C
  172. C        APPLY THE TRANSFORMATION TO S AND REDUCE THE SPIKE IN W.
  173. C
  174.          L = JJ
  175.          DO 110 I = J, M
  176.             TEMP = COS*S(L) + SIN*W(I)
  177.             W(I) = -SIN*S(L) + COS*W(I)
  178.             S(L) = TEMP
  179.             L = L + 1
  180.   110       CONTINUE
  181. C
  182. C        STORE THE INFORMATION NECESSARY TO RECOVER THE
  183. C        GIVENS ROTATION.
  184. C
  185.          W(J) = TAU
  186.   120    CONTINUE
  187. C
  188. C        TEST FOR ZERO DIAGONAL ELEMENTS IN THE OUTPUT S.
  189. C
  190.          IF (S(JJ) .EQ. ZERO) SING = .TRUE.
  191.          JJ = JJ + (M - J + 1)
  192.   130    CONTINUE
  193.   140 CONTINUE
  194. C
  195. C     MOVE W BACK INTO THE LAST COLUMN OF THE OUTPUT S.
  196. C
  197.       L = JJ
  198.       DO 150 I = N, M
  199.          S(L) = W(I)
  200.          L = L + 1
  201.   150    CONTINUE
  202.       IF (S(JJ) .EQ. ZERO) SING = .TRUE.
  203.       RETURN
  204. C
  205. C     LAST CARD OF SUBROUTINE R1UPDT.
  206. C
  207.       END
  208.