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 / lapack / dlasv2.f < prev    next >
Text File  |  1996-09-28  |  7KB  |  251 lines
  1.       SUBROUTINE DLASV2( F, G, H, SSMIN, SSMAX, SNR, CSR, SNL, CSL )
  2. *
  3. *  -- LAPACK auxiliary routine (version 2.0) --
  4. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  5. *     Courant Institute, Argonne National Lab, and Rice University
  6. *     October 31, 1992
  7. *
  8. *     .. Scalar Arguments ..
  9.       DOUBLE PRECISION   CSL, CSR, F, G, H, SNL, SNR, SSMAX, SSMIN
  10. *     ..
  11. *
  12. *  Purpose
  13. *  =======
  14. *
  15. *  DLASV2 computes the singular value decomposition of a 2-by-2
  16. *  triangular matrix
  17. *     [  F   G  ]
  18. *     [  0   H  ].
  19. *  On return, abs(SSMAX) is the larger singular value, abs(SSMIN) is the
  20. *  smaller singular value, and (CSL,SNL) and (CSR,SNR) are the left and
  21. *  right singular vectors for abs(SSMAX), giving the decomposition
  22. *
  23. *     [ CSL  SNL ] [  F   G  ] [ CSR -SNR ]  =  [ SSMAX   0   ]
  24. *     [-SNL  CSL ] [  0   H  ] [ SNR  CSR ]     [  0    SSMIN ].
  25. *
  26. *  Arguments
  27. *  =========
  28. *
  29. *  F       (input) DOUBLE PRECISION
  30. *          The (1,1) element of the 2-by-2 matrix.
  31. *
  32. *  G       (input) DOUBLE PRECISION
  33. *          The (1,2) element of the 2-by-2 matrix.
  34. *
  35. *  H       (input) DOUBLE PRECISION
  36. *          The (2,2) element of the 2-by-2 matrix.
  37. *
  38. *  SSMIN   (output) DOUBLE PRECISION
  39. *          abs(SSMIN) is the smaller singular value.
  40. *
  41. *  SSMAX   (output) DOUBLE PRECISION
  42. *          abs(SSMAX) is the larger singular value.
  43. *
  44. *  SNL     (output) DOUBLE PRECISION
  45. *  CSL     (output) DOUBLE PRECISION
  46. *          The vector (CSL, SNL) is a unit left singular vector for the
  47. *          singular value abs(SSMAX).
  48. *
  49. *  SNR     (output) DOUBLE PRECISION
  50. *  CSR     (output) DOUBLE PRECISION
  51. *          The vector (CSR, SNR) is a unit right singular vector for the
  52. *          singular value abs(SSMAX).
  53. *
  54. *  Further Details
  55. *  ===============
  56. *
  57. *  Any input parameter may be aliased with any output parameter.
  58. *
  59. *  Barring over/underflow and assuming a guard digit in subtraction, all
  60. *  output quantities are correct to within a few units in the last
  61. *  place (ulps).
  62. *
  63. *  In IEEE arithmetic, the code works correctly if one matrix element is
  64. *  infinite.
  65. *
  66. *  Overflow will not occur unless the largest singular value itself
  67. *  overflows or is within a few ulps of overflow. (On machines with
  68. *  partial overflow, like the Cray, overflow may occur if the largest
  69. *  singular value is within a factor of 2 of overflow.)
  70. *
  71. *  Underflow is harmless if underflow is gradual. Otherwise, results
  72. *  may correspond to a matrix modified by perturbations of size near
  73. *  the underflow threshold.
  74. *
  75. * =====================================================================
  76. *
  77. *     .. Parameters ..
  78.       DOUBLE PRECISION   ZERO
  79.       PARAMETER          ( ZERO = 0.0D0 )
  80.       DOUBLE PRECISION   HALF
  81.       PARAMETER          ( HALF = 0.5D0 )
  82.       DOUBLE PRECISION   ONE
  83.       PARAMETER          ( ONE = 1.0D0 )
  84.       DOUBLE PRECISION   TWO
  85.       PARAMETER          ( TWO = 2.0D0 )
  86.       DOUBLE PRECISION   FOUR
  87.       PARAMETER          ( FOUR = 4.0D0 )
  88. *     ..
  89. *     .. Local Scalars ..
  90.       LOGICAL            GASMAL, SWAP
  91.       INTEGER            PMAX
  92.       DOUBLE PRECISION   A, CLT, CRT, D, FA, FT, GA, GT, HA, HT, L, M,
  93.      $                   MM, R, S, SLT, SRT, T, TEMP, TSIGN, TT
  94. *     ..
  95. *     .. Intrinsic Functions ..
  96.       INTRINSIC          ABS, SIGN, SQRT
  97. *     ..
  98. *     .. External Functions ..
  99.       DOUBLE PRECISION   DLAMCH
  100.       EXTERNAL           DLAMCH
  101. *     ..
  102. *     .. Executable Statements ..
  103. *
  104.       FT = F
  105.       FA = ABS( FT )
  106.       HT = H
  107.       HA = ABS( H )
  108. *
  109. *     PMAX points to the maximum absolute element of matrix
  110. *       PMAX = 1 if F largest in absolute values
  111. *       PMAX = 2 if G largest in absolute values
  112. *       PMAX = 3 if H largest in absolute values
  113. *
  114.       PMAX = 1
  115.       SWAP = ( HA.GT.FA )
  116.       IF( SWAP ) THEN
  117.          PMAX = 3
  118.          TEMP = FT
  119.          FT = HT
  120.          HT = TEMP
  121.          TEMP = FA
  122.          FA = HA
  123.          HA = TEMP
  124. *
  125. *        Now FA .ge. HA
  126. *
  127.       END IF
  128.       GT = G
  129.       GA = ABS( GT )
  130.       IF( GA.EQ.ZERO ) THEN
  131. *
  132. *        Diagonal matrix
  133. *
  134.          SSMIN = HA
  135.          SSMAX = FA
  136.          CLT = ONE
  137.          CRT = ONE
  138.          SLT = ZERO
  139.          SRT = ZERO
  140.       ELSE
  141.          GASMAL = .TRUE.
  142.          IF( GA.GT.FA ) THEN
  143.             PMAX = 2
  144.             IF( ( FA / GA ).LT.DLAMCH( 'EPS' ) ) THEN
  145. *
  146. *              Case of very large GA
  147. *
  148.                GASMAL = .FALSE.
  149.                SSMAX = GA
  150.                IF( HA.GT.ONE ) THEN
  151.                   SSMIN = FA / ( GA / HA )
  152.                ELSE
  153.                   SSMIN = ( FA / GA )*HA
  154.                END IF
  155.                CLT = ONE
  156.                SLT = HT / GT
  157.                SRT = ONE
  158.                CRT = FT / GT
  159.             END IF
  160.          END IF
  161.          IF( GASMAL ) THEN
  162. *
  163. *           Normal case
  164. *
  165.             D = FA - HA
  166.             IF( D.EQ.FA ) THEN
  167. *
  168. *              Copes with infinite F or H
  169. *
  170.                L = ONE
  171.             ELSE
  172.                L = D / FA
  173.             END IF
  174. *
  175. *           Note that 0 .le. L .le. 1
  176. *
  177.             M = GT / FT
  178. *
  179. *           Note that abs(M) .le. 1/macheps
  180. *
  181.             T = TWO - L
  182. *
  183. *           Note that T .ge. 1
  184. *
  185.             MM = M*M
  186.             TT = T*T
  187.             S = SQRT( TT+MM )
  188. *
  189. *           Note that 1 .le. S .le. 1 + 1/macheps
  190. *
  191.             IF( L.EQ.ZERO ) THEN
  192.                R = ABS( M )
  193.             ELSE
  194.                R = SQRT( L*L+MM )
  195.             END IF
  196. *
  197. *           Note that 0 .le. R .le. 1 + 1/macheps
  198. *
  199.             A = HALF*( S+R )
  200. *
  201. *           Note that 1 .le. A .le. 1 + abs(M)
  202. *
  203.             SSMIN = HA / A
  204.             SSMAX = FA*A
  205.             IF( MM.EQ.ZERO ) THEN
  206. *
  207. *              Note that M is very tiny
  208. *
  209.                IF( L.EQ.ZERO ) THEN
  210.                   T = SIGN( TWO, FT )*SIGN( ONE, GT )
  211.                ELSE
  212.                   T = GT / SIGN( D, FT ) + M / T
  213.                END IF
  214.             ELSE
  215.                T = ( M / ( S+T )+M / ( R+L ) )*( ONE+A )
  216.             END IF
  217.             L = SQRT( T*T+FOUR )
  218.             CRT = TWO / L
  219.             SRT = T / L
  220.             CLT = ( CRT+SRT*M ) / A
  221.             SLT = ( HT / FT )*SRT / A
  222.          END IF
  223.       END IF
  224.       IF( SWAP ) THEN
  225.          CSL = SRT
  226.          SNL = CRT
  227.          CSR = SLT
  228.          SNR = CLT
  229.       ELSE
  230.          CSL = CLT
  231.          SNL = SLT
  232.          CSR = CRT
  233.          SNR = SRT
  234.       END IF
  235. *
  236. *     Correct signs of SSMAX and SSMIN
  237. *
  238.       IF( PMAX.EQ.1 )
  239.      $   TSIGN = SIGN( ONE, CSR )*SIGN( ONE, CSL )*SIGN( ONE, F )
  240.       IF( PMAX.EQ.2 )
  241.      $   TSIGN = SIGN( ONE, SNR )*SIGN( ONE, CSL )*SIGN( ONE, G )
  242.       IF( PMAX.EQ.3 )
  243.      $   TSIGN = SIGN( ONE, SNR )*SIGN( ONE, SNL )*SIGN( ONE, H )
  244.       SSMAX = SIGN( SSMAX, TSIGN )
  245.       SSMIN = SIGN( SSMIN, TSIGN*SIGN( ONE, F )*SIGN( ONE, H ) )
  246.       RETURN
  247. *
  248. *     End of DLASV2
  249. *
  250.       END
  251.