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 December 2016
111: *
112: *> \ingroup OTHERauxiliary
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.7.0) --
142: * -- LAPACK is a software package provided by Univ. of Tennessee, --
143: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
144: * December 2016
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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlasv2.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack