LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
slasv2.f
Go to the documentation of this file.
1 *> \brief \b SLASV2 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 SLASV2 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slasv2.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slasv2.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slasv2.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SLASV2( F, G, H, SSMIN, SSMAX, SNR, CSR, SNL, CSL )
22 *
23 * .. Scalar Arguments ..
24 * REAL CSL, CSR, F, G, H, SNL, SNR, SSMAX, SSMIN
25 * ..
26 *
27 *
28 *> \par Purpose:
29 * =============
30 *>
31 *> \verbatim
32 *>
33 *> SLASV2 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 REAL
51 *> The (1,1) element of the 2-by-2 matrix.
52 *> \endverbatim
53 *>
54 *> \param[in] G
55 *> \verbatim
56 *> G is REAL
57 *> The (1,2) element of the 2-by-2 matrix.
58 *> \endverbatim
59 *>
60 *> \param[in] H
61 *> \verbatim
62 *> H is REAL
63 *> The (2,2) element of the 2-by-2 matrix.
64 *> \endverbatim
65 *>
66 *> \param[out] SSMIN
67 *> \verbatim
68 *> SSMIN is REAL
69 *> abs(SSMIN) is the smaller singular value.
70 *> \endverbatim
71 *>
72 *> \param[out] SSMAX
73 *> \verbatim
74 *> SSMAX is REAL
75 *> abs(SSMAX) is the larger singular value.
76 *> \endverbatim
77 *>
78 *> \param[out] SNL
79 *> \verbatim
80 *> SNL is REAL
81 *> \endverbatim
82 *>
83 *> \param[out] CSL
84 *> \verbatim
85 *> CSL is REAL
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 REAL
93 *> \endverbatim
94 *>
95 *> \param[out] CSR
96 *> \verbatim
97 *> CSR is REAL
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 *> \ingroup OTHERauxiliary
111 *
112 *> \par Further Details:
113 * =====================
114 *>
115 *> \verbatim
116 *>
117 *> Any input parameter may be aliased with any output parameter.
118 *>
119 *> Barring over/underflow and assuming a guard digit in subtraction, all
120 *> output quantities are correct to within a few units in the last
121 *> place (ulps).
122 *>
123 *> In IEEE arithmetic, the code works correctly if one matrix element is
124 *> infinite.
125 *>
126 *> Overflow will not occur unless the largest singular value itself
127 *> overflows or is within a few ulps of overflow. (On machines with
128 *> partial overflow, like the Cray, overflow may occur if the largest
129 *> singular value is within a factor of 2 of overflow.)
130 *>
131 *> Underflow is harmless if underflow is gradual. Otherwise, results
132 *> may correspond to a matrix modified by perturbations of size near
133 *> the underflow threshold.
134 *> \endverbatim
135 *>
136 * =====================================================================
137  SUBROUTINE slasv2( F, G, H, SSMIN, SSMAX, SNR, CSR, SNL, CSL )
138 *
139 * -- LAPACK auxiliary routine --
140 * -- LAPACK is a software package provided by Univ. of Tennessee, --
141 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
142 *
143 * .. Scalar Arguments ..
144  REAL CSL, CSR, F, G, H, SNL, SNR, SSMAX, SSMIN
145 * ..
146 *
147 * =====================================================================
148 *
149 * .. Parameters ..
150  REAL ZERO
151  parameter( zero = 0.0e0 )
152  REAL HALF
153  parameter( half = 0.5e0 )
154  REAL ONE
155  parameter( one = 1.0e0 )
156  REAL TWO
157  parameter( two = 2.0e0 )
158  REAL FOUR
159  parameter( four = 4.0e0 )
160 * ..
161 * .. Local Scalars ..
162  LOGICAL GASMAL, SWAP
163  INTEGER PMAX
164  REAL A, CLT, CRT, D, FA, FT, GA, GT, HA, HT, L, M,
165  $ MM, R, S, SLT, SRT, T, TEMP, TSIGN, TT
166 * ..
167 * .. Intrinsic Functions ..
168  INTRINSIC abs, sign, sqrt
169 * ..
170 * .. External Functions ..
171  REAL SLAMCH
172  EXTERNAL slamch
173 * ..
174 * .. Executable Statements ..
175 *
176  ft = f
177  fa = abs( ft )
178  ht = h
179  ha = abs( h )
180 *
181 * PMAX points to the maximum absolute element of matrix
182 * PMAX = 1 if F largest in absolute values
183 * PMAX = 2 if G largest in absolute values
184 * PMAX = 3 if H largest in absolute values
185 *
186  pmax = 1
187  swap = ( ha.GT.fa )
188  IF( swap ) THEN
189  pmax = 3
190  temp = ft
191  ft = ht
192  ht = temp
193  temp = fa
194  fa = ha
195  ha = temp
196 *
197 * Now FA .ge. HA
198 *
199  END IF
200  gt = g
201  ga = abs( gt )
202  IF( ga.EQ.zero ) THEN
203 *
204 * Diagonal matrix
205 *
206  ssmin = ha
207  ssmax = fa
208  clt = one
209  crt = one
210  slt = zero
211  srt = zero
212  ELSE
213  gasmal = .true.
214  IF( ga.GT.fa ) THEN
215  pmax = 2
216  IF( ( fa / ga ).LT.slamch( 'EPS' ) ) THEN
217 *
218 * Case of very large GA
219 *
220  gasmal = .false.
221  ssmax = ga
222  IF( ha.GT.one ) THEN
223  ssmin = fa / ( ga / ha )
224  ELSE
225  ssmin = ( fa / ga )*ha
226  END IF
227  clt = one
228  slt = ht / gt
229  srt = one
230  crt = ft / gt
231  END IF
232  END IF
233  IF( gasmal ) THEN
234 *
235 * Normal case
236 *
237  d = fa - ha
238  IF( d.EQ.fa ) THEN
239 *
240 * Copes with infinite F or H
241 *
242  l = one
243  ELSE
244  l = d / fa
245  END IF
246 *
247 * Note that 0 .le. L .le. 1
248 *
249  m = gt / ft
250 *
251 * Note that abs(M) .le. 1/macheps
252 *
253  t = two - l
254 *
255 * Note that T .ge. 1
256 *
257  mm = m*m
258  tt = t*t
259  s = sqrt( tt+mm )
260 *
261 * Note that 1 .le. S .le. 1 + 1/macheps
262 *
263  IF( l.EQ.zero ) THEN
264  r = abs( m )
265  ELSE
266  r = sqrt( l*l+mm )
267  END IF
268 *
269 * Note that 0 .le. R .le. 1 + 1/macheps
270 *
271  a = half*( s+r )
272 *
273 * Note that 1 .le. A .le. 1 + abs(M)
274 *
275  ssmin = ha / a
276  ssmax = fa*a
277  IF( mm.EQ.zero ) THEN
278 *
279 * Note that M is very tiny
280 *
281  IF( l.EQ.zero ) THEN
282  t = sign( two, ft )*sign( one, gt )
283  ELSE
284  t = gt / sign( d, ft ) + m / t
285  END IF
286  ELSE
287  t = ( m / ( s+t )+m / ( r+l ) )*( one+a )
288  END IF
289  l = sqrt( t*t+four )
290  crt = two / l
291  srt = t / l
292  clt = ( crt+srt*m ) / a
293  slt = ( ht / ft )*srt / a
294  END IF
295  END IF
296  IF( swap ) THEN
297  csl = srt
298  snl = crt
299  csr = slt
300  snr = clt
301  ELSE
302  csl = clt
303  snl = slt
304  csr = crt
305  snr = srt
306  END IF
307 *
308 * Correct signs of SSMAX and SSMIN
309 *
310  IF( pmax.EQ.1 )
311  $ tsign = sign( one, csr )*sign( one, csl )*sign( one, f )
312  IF( pmax.EQ.2 )
313  $ tsign = sign( one, snr )*sign( one, csl )*sign( one, g )
314  IF( pmax.EQ.3 )
315  $ tsign = sign( one, snr )*sign( one, snl )*sign( one, h )
316  ssmax = sign( ssmax, tsign )
317  ssmin = sign( ssmin, tsign*sign( one, f )*sign( one, h ) )
318  RETURN
319 *
320 * End of SLASV2
321 *
322  END
subroutine slasv2(F, G, H, SSMIN, SSMAX, SNR, CSR, SNL, CSL)
SLASV2 computes the singular value decomposition of a 2-by-2 triangular matrix.
Definition: slasv2.f:138