LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches
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
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*> \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 dlasv2( 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 DOUBLE PRECISION CSL, CSR, F, G, H, SNL, SNR, SSMAX, SSMIN
145* ..
146*
147* =====================================================================
148*
149* .. Parameters ..
150 DOUBLE PRECISION ZERO
151 parameter( zero = 0.0d0 )
152 DOUBLE PRECISION HALF
153 parameter( half = 0.5d0 )
154 DOUBLE PRECISION ONE
155 parameter( one = 1.0d0 )
156 DOUBLE PRECISION TWO
157 parameter( two = 2.0d0 )
158 DOUBLE PRECISION FOUR
159 parameter( four = 4.0d0 )
160* ..
161* .. Local Scalars ..
162 LOGICAL GASMAL, SWAP
163 INTEGER PMAX
164 DOUBLE PRECISION 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 DOUBLE PRECISION DLAMCH
172 EXTERNAL dlamch
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.dlamch( '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 DLASV2
321*
322 END
subroutine dlasv2(F, G, H, SSMIN, SSMAX, SNR, CSR, SNL, CSL)
DLASV2 computes the singular value decomposition of a 2-by-2 triangular matrix.
Definition: dlasv2.f:138