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: *> \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
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