LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
dlasv2.f
Go to the documentation of this file.
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