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 / dlas2.f < prev    next >
Text File  |  1996-09-28  |  4KB  |  123 lines
  1.       SUBROUTINE DLAS2( F, G, H, SSMIN, SSMAX )
  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. *     September 30, 1994
  7. *
  8. *     .. Scalar Arguments ..
  9.       DOUBLE PRECISION   F, G, H, SSMAX, SSMIN
  10. *     ..
  11. *
  12. *  Purpose
  13. *  =======
  14. *
  15. *  DLAS2  computes the singular values of the 2-by-2 matrix
  16. *     [  F   G  ]
  17. *     [  0   H  ].
  18. *  On return, SSMIN is the smaller singular value and SSMAX is the
  19. *  larger singular value.
  20. *
  21. *  Arguments
  22. *  =========
  23. *
  24. *  F       (input) DOUBLE PRECISION
  25. *          The (1,1) element of the 2-by-2 matrix.
  26. *
  27. *  G       (input) DOUBLE PRECISION
  28. *          The (1,2) element of the 2-by-2 matrix.
  29. *
  30. *  H       (input) DOUBLE PRECISION
  31. *          The (2,2) element of the 2-by-2 matrix.
  32. *
  33. *  SSMIN   (output) DOUBLE PRECISION
  34. *          The smaller singular value.
  35. *
  36. *  SSMAX   (output) DOUBLE PRECISION
  37. *          The larger singular value.
  38. *
  39. *  Further Details
  40. *  ===============
  41. *
  42. *  Barring over/underflow, all output quantities are correct to within
  43. *  a few units in the last place (ulps), even in the absence of a guard
  44. *  digit in addition/subtraction.
  45. *
  46. *  In IEEE arithmetic, the code works correctly if one matrix element is
  47. *  infinite.
  48. *
  49. *  Overflow will not occur unless the largest singular value itself
  50. *  overflows, or is within a few ulps of overflow. (On machines with
  51. *  partial overflow, like the Cray, overflow may occur if the largest
  52. *  singular value is within a factor of 2 of overflow.)
  53. *
  54. *  Underflow is harmless if underflow is gradual. Otherwise, results
  55. *  may correspond to a matrix modified by perturbations of size near
  56. *  the underflow threshold.
  57. *
  58. *  ====================================================================
  59. *
  60. *     .. Parameters ..
  61.       DOUBLE PRECISION   ZERO
  62.       PARAMETER          ( ZERO = 0.0D0 )
  63.       DOUBLE PRECISION   ONE
  64.       PARAMETER          ( ONE = 1.0D0 )
  65.       DOUBLE PRECISION   TWO
  66.       PARAMETER          ( TWO = 2.0D0 )
  67. *     ..
  68. *     .. Local Scalars ..
  69.       DOUBLE PRECISION   AS, AT, AU, C, FA, FHMN, FHMX, GA, HA
  70. *     ..
  71. *     .. Intrinsic Functions ..
  72.       INTRINSIC          ABS, MAX, MIN, SQRT
  73. *     ..
  74. *     .. Executable Statements ..
  75. *
  76.       FA = ABS( F )
  77.       GA = ABS( G )
  78.       HA = ABS( H )
  79.       FHMN = MIN( FA, HA )
  80.       FHMX = MAX( FA, HA )
  81.       IF( FHMN.EQ.ZERO ) THEN
  82.          SSMIN = ZERO
  83.          IF( FHMX.EQ.ZERO ) THEN
  84.             SSMAX = GA
  85.          ELSE
  86.             SSMAX = MAX( FHMX, GA )*SQRT( ONE+
  87.      $              ( MIN( FHMX, GA ) / MAX( FHMX, GA ) )**2 )
  88.          END IF
  89.       ELSE
  90.          IF( GA.LT.FHMX ) THEN
  91.             AS = ONE + FHMN / FHMX
  92.             AT = ( FHMX-FHMN ) / FHMX
  93.             AU = ( GA / FHMX )**2
  94.             C = TWO / ( SQRT( AS*AS+AU )+SQRT( AT*AT+AU ) )
  95.             SSMIN = FHMN*C
  96.             SSMAX = FHMX / C
  97.          ELSE
  98.             AU = FHMX / GA
  99.             IF( AU.EQ.ZERO ) THEN
  100. *
  101. *              Avoid possible harmful underflow if exponent range
  102. *              asymmetric (true SSMIN may not underflow even if
  103. *              AU underflows)
  104. *
  105.                SSMIN = ( FHMN*FHMX ) / GA
  106.                SSMAX = GA
  107.             ELSE
  108.                AS = ONE + FHMN / FHMX
  109.                AT = ( FHMX-FHMN ) / FHMX
  110.                C = ONE / ( SQRT( ONE+( AS*AU )**2 )+
  111.      $             SQRT( ONE+( AT*AU )**2 ) )
  112.                SSMIN = ( FHMN*C )*AU
  113.                SSMIN = SSMIN + SSMIN
  114.                SSMAX = GA / ( C+C )
  115.             END IF
  116.          END IF
  117.       END IF
  118.       RETURN
  119. *
  120. *     End of DLAS2
  121. *
  122.       END
  123.