Statistiques
| Révision :

root / src / lapack / double / dlas2.f @ 10

Historique | Voir | Annoter | Télécharger (3,55 ko)

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