1*> \brief \b DLASV2 computes the singular value decomposition of a 2-by-2 triangular matrix.
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DLASV2 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasv2.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasv2.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasv2.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE DLASV2( F, G, H, SSMIN, SSMAX, SNR, CSR, SNL, CSL )
22*
23*       .. Scalar Arguments ..
24*       DOUBLE PRECISION   CSL, CSR, F, G, H, SNL, SNR, SSMAX, SSMIN
25*       ..
26*
27*
28*> \par Purpose:
29*  =============
30*>
31*> \verbatim
32*>
33*> DLASV2 computes the singular value decomposition of a 2-by-2
34*> triangular matrix
35*>    [  F   G  ]
36*>    [  0   H  ].
37*> On return, abs(SSMAX) is the larger singular value, abs(SSMIN) is the
38*> smaller singular value, and (CSL,SNL) and (CSR,SNR) are the left and
39*> right singular vectors for abs(SSMAX), giving the decomposition
40*>
41*>    [ CSL  SNL ] [  F   G  ] [ CSR -SNR ]  =  [ SSMAX   0   ]
42*>    [-SNL  CSL ] [  0   H  ] [ SNR  CSR ]     [  0    SSMIN ].
43*> \endverbatim
44*
45*  Arguments:
46*  ==========
47*
48*> \param[in] F
49*> \verbatim
50*>          F is DOUBLE PRECISION
51*>          The (1,1) element of the 2-by-2 matrix.
52*> \endverbatim
53*>
54*> \param[in] G
55*> \verbatim
56*>          G is DOUBLE PRECISION
57*>          The (1,2) element of the 2-by-2 matrix.
58*> \endverbatim
59*>
60*> \param[in] H
61*> \verbatim
62*>          H is DOUBLE PRECISION
63*>          The (2,2) element of the 2-by-2 matrix.
64*> \endverbatim
65*>
66*> \param[out] SSMIN
67*> \verbatim
68*>          SSMIN is DOUBLE PRECISION
69*>          abs(SSMIN) is the smaller singular value.
70*> \endverbatim
71*>
72*> \param[out] SSMAX
73*> \verbatim
74*>          SSMAX is DOUBLE PRECISION
75*>          abs(SSMAX) is the larger singular value.
76*> \endverbatim
77*>
78*> \param[out] SNL
79*> \verbatim
80*>          SNL is DOUBLE PRECISION
81*> \endverbatim
82*>
83*> \param[out] CSL
84*> \verbatim
85*>          CSL is DOUBLE PRECISION
86*>          The vector (CSL, SNL) is a unit left singular vector for the
87*>          singular value abs(SSMAX).
88*> \endverbatim
89*>
90*> \param[out] SNR
91*> \verbatim
92*>          SNR is DOUBLE PRECISION
93*> \endverbatim
94*>
95*> \param[out] CSR
96*> \verbatim
97*>          CSR is DOUBLE PRECISION
98*>          The vector (CSR, SNR) is a unit right singular vector for the
99*>          singular value abs(SSMAX).
100*> \endverbatim
101*
102*  Authors:
103*  ========
104*
105*> \author Univ. of Tennessee
106*> \author Univ. of California Berkeley
107*> \author Univ. of Colorado Denver
108*> \author NAG Ltd.
109*
110*> \date September 2012
111*
112*> \ingroup auxOTHERauxiliary
113*
114*> \par Further Details:
115*  =====================
116*>
117*> \verbatim
118*>
119*>  Any input parameter may be aliased with any output parameter.
120*>
121*>  Barring over/underflow and assuming a guard digit in subtraction, all
122*>  output quantities are correct to within a few units in the last
123*>  place (ulps).
124*>
125*>  In IEEE arithmetic, the code works correctly if one matrix element is
126*>  infinite.
127*>
128*>  Overflow will not occur unless the largest singular value itself
129*>  overflows or is within a few ulps of overflow. (On machines with
130*>  partial overflow, like the Cray, overflow may occur if the largest
131*>  singular value is within a factor of 2 of overflow.)
132*>
133*>  Underflow is harmless if underflow is gradual. Otherwise, results
134*>  may correspond to a matrix modified by perturbations of size near
135*>  the underflow threshold.
136*> \endverbatim
137*>
138*  =====================================================================
139      SUBROUTINE DLASV2( F, G, H, SSMIN, SSMAX, SNR, CSR, SNL, CSL )
140*
141*  -- LAPACK auxiliary routine (version 3.4.2) --
142*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
143*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
144*     September 2012
145*
146*     .. Scalar Arguments ..
147      DOUBLE PRECISION   CSL, CSR, F, G, H, SNL, SNR, SSMAX, SSMIN
148*     ..
149*
150* =====================================================================
151*
152*     .. Parameters ..
153      DOUBLE PRECISION   ZERO
154      PARAMETER          ( ZERO = 0.0D0 )
155      DOUBLE PRECISION   HALF
156      PARAMETER          ( HALF = 0.5D0 )
157      DOUBLE PRECISION   ONE
158      PARAMETER          ( ONE = 1.0D0 )
159      DOUBLE PRECISION   TWO
160      PARAMETER          ( TWO = 2.0D0 )
161      DOUBLE PRECISION   FOUR
162      PARAMETER          ( FOUR = 4.0D0 )
163*     ..
164*     .. Local Scalars ..
165      LOGICAL            GASMAL, SWAP
166      INTEGER            PMAX
167      DOUBLE PRECISION   A, CLT, CRT, D, FA, FT, GA, GT, HA, HT, L, M,
168     $                   MM, R, S, SLT, SRT, T, TEMP, TSIGN, TT
169*     ..
170*     .. Intrinsic Functions ..
171      INTRINSIC          ABS, SIGN, SQRT
172*     ..
173*     .. External Functions ..
174      DOUBLE PRECISION   DLAMCH
175      EXTERNAL           DLAMCH
176*     ..
177*     .. Executable Statements ..
178*
179      FT = F
180      FA = ABS( FT )
181      HT = H
182      HA = ABS( H )
183*
184*     PMAX points to the maximum absolute element of matrix
185*       PMAX = 1 if F largest in absolute values
186*       PMAX = 2 if G largest in absolute values
187*       PMAX = 3 if H largest in absolute values
188*
189      PMAX = 1
190      SWAP = ( HA.GT.FA )
191      IF( SWAP ) THEN
192         PMAX = 3
193         TEMP = FT
194         FT = HT
195         HT = TEMP
196         TEMP = FA
197         FA = HA
198         HA = TEMP
199*
200*        Now FA .ge. HA
201*
202      END IF
203      GT = G
204      GA = ABS( GT )
205      IF( GA.EQ.ZERO ) THEN
206*
207*        Diagonal matrix
208*
209         SSMIN = HA
210         SSMAX = FA
211         CLT = ONE
212         CRT = ONE
213         SLT = ZERO
214         SRT = ZERO
215      ELSE
216         GASMAL = .TRUE.
217         IF( GA.GT.FA ) THEN
218            PMAX = 2
219            IF( ( FA / GA ).LT.DLAMCH( 'EPS' ) ) THEN
220*
221*              Case of very large GA
222*
223               GASMAL = .FALSE.
224               SSMAX = GA
225               IF( HA.GT.ONE ) THEN
226                  SSMIN = FA / ( GA / HA )
227               ELSE
228                  SSMIN = ( FA / GA )*HA
229               END IF
230               CLT = ONE
231               SLT = HT / GT
232               SRT = ONE
233               CRT = FT / GT
234            END IF
235         END IF
236         IF( GASMAL ) THEN
237*
238*           Normal case
239*
240            D = FA - HA
241            IF( D.EQ.FA ) THEN
242*
243*              Copes with infinite F or H
244*
245               L = ONE
246            ELSE
247               L = D / FA
248            END IF
249*
250*           Note that 0 .le. L .le. 1
251*
252            M = GT / FT
253*
254*           Note that abs(M) .le. 1/macheps
255*
256            T = TWO - L
257*
258*           Note that T .ge. 1
259*
260            MM = M*M
261            TT = T*T
262            S = SQRT( TT+MM )
263*
264*           Note that 1 .le. S .le. 1 + 1/macheps
265*
266            IF( L.EQ.ZERO ) THEN
267               R = ABS( M )
268            ELSE
269               R = SQRT( L*L+MM )
270            END IF
271*
272*           Note that 0 .le. R .le. 1 + 1/macheps
273*
274            A = HALF*( S+R )
275*
276*           Note that 1 .le. A .le. 1 + abs(M)
277*
278            SSMIN = HA / A
279            SSMAX = FA*A
280            IF( MM.EQ.ZERO ) THEN
281*
282*              Note that M is very tiny
283*
284               IF( L.EQ.ZERO ) THEN
285                  T = SIGN( TWO, FT )*SIGN( ONE, GT )
286               ELSE
287                  T = GT / SIGN( D, FT ) + M / T
288               END IF
289            ELSE
290               T = ( M / ( S+T )+M / ( R+L ) )*( ONE+A )
291            END IF
292            L = SQRT( T*T+FOUR )
293            CRT = TWO / L
294            SRT = T / L
295            CLT = ( CRT+SRT*M ) / A
296            SLT = ( HT / FT )*SRT / A
297         END IF
298      END IF
299      IF( SWAP ) THEN
300         CSL = SRT
301         SNL = CRT
302         CSR = SLT
303         SNR = CLT
304      ELSE
305         CSL = CLT
306         SNL = SLT
307         CSR = CRT
308         SNR = SRT
309      END IF
310*
311*     Correct signs of SSMAX and SSMIN
312*
313      IF( PMAX.EQ.1 )
314     $   TSIGN = SIGN( ONE, CSR )*SIGN( ONE, CSL )*SIGN( ONE, F )
315      IF( PMAX.EQ.2 )
316     $   TSIGN = SIGN( ONE, SNR )*SIGN( ONE, CSL )*SIGN( ONE, G )
317      IF( PMAX.EQ.3 )
318     $   TSIGN = SIGN( ONE, SNR )*SIGN( ONE, SNL )*SIGN( ONE, H )
319      SSMAX = SIGN( SSMAX, TSIGN )
320      SSMIN = SIGN( SSMIN, TSIGN*SIGN( ONE, F )*SIGN( ONE, H ) )
321      RETURN
322*
323*     End of DLASV2
324*
325      END
326