![[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|
Return to dlasv2.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 1: SUBROUTINE DLASV2( F, G, H, SSMIN, SSMAX, SNR, CSR, SNL, CSL )
2: *
3: * -- LAPACK auxiliary routine (version 3.2) --
4: * -- LAPACK is a software package provided by Univ. of Tennessee, --
5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6: * November 2006
7: *
8: * .. Scalar Arguments ..
9: DOUBLE PRECISION CSL, CSR, F, G, H, SNL, SNR, SSMAX, SSMIN
10: * ..
11: *
12: * Purpose
13: * =======
14: *
15: * DLASV2 computes the singular value decomposition of a 2-by-2
16: * triangular matrix
17: * [ F G ]
18: * [ 0 H ].
19: * On return, abs(SSMAX) is the larger singular value, abs(SSMIN) is the
20: * smaller singular value, and (CSL,SNL) and (CSR,SNR) are the left and
21: * right singular vectors for abs(SSMAX), giving the decomposition
22: *
23: * [ CSL SNL ] [ F G ] [ CSR -SNR ] = [ SSMAX 0 ]
24: * [-SNL CSL ] [ 0 H ] [ SNR CSR ] [ 0 SSMIN ].
25: *
26: * Arguments
27: * =========
28: *
29: * F (input) DOUBLE PRECISION
30: * The (1,1) element of the 2-by-2 matrix.
31: *
32: * G (input) DOUBLE PRECISION
33: * The (1,2) element of the 2-by-2 matrix.
34: *
35: * H (input) DOUBLE PRECISION
36: * The (2,2) element of the 2-by-2 matrix.
37: *
38: * SSMIN (output) DOUBLE PRECISION
39: * abs(SSMIN) is the smaller singular value.
40: *
41: * SSMAX (output) DOUBLE PRECISION
42: * abs(SSMAX) is the larger singular value.
43: *
44: * SNL (output) DOUBLE PRECISION
45: * CSL (output) DOUBLE PRECISION
46: * The vector (CSL, SNL) is a unit left singular vector for the
47: * singular value abs(SSMAX).
48: *
49: * SNR (output) DOUBLE PRECISION
50: * CSR (output) DOUBLE PRECISION
51: * The vector (CSR, SNR) is a unit right singular vector for the
52: * singular value abs(SSMAX).
53: *
54: * Further Details
55: * ===============
56: *
57: * Any input parameter may be aliased with any output parameter.
58: *
59: * Barring over/underflow and assuming a guard digit in subtraction, all
60: * output quantities are correct to within a few units in the last
61: * place (ulps).
62: *
63: * In IEEE arithmetic, the code works correctly if one matrix element is
64: * infinite.
65: *
66: * Overflow will not occur unless the largest singular value itself
67: * overflows or is within a few ulps of overflow. (On machines with
68: * partial overflow, like the Cray, overflow may occur if the largest
69: * singular value is within a factor of 2 of overflow.)
70: *
71: * Underflow is harmless if underflow is gradual. Otherwise, results
72: * may correspond to a matrix modified by perturbations of size near
73: * the underflow threshold.
74: *
75: * =====================================================================
76: *
77: * .. Parameters ..
78: DOUBLE PRECISION ZERO
79: PARAMETER ( ZERO = 0.0D0 )
80: DOUBLE PRECISION HALF
81: PARAMETER ( HALF = 0.5D0 )
82: DOUBLE PRECISION ONE
83: PARAMETER ( ONE = 1.0D0 )
84: DOUBLE PRECISION TWO
85: PARAMETER ( TWO = 2.0D0 )
86: DOUBLE PRECISION FOUR
87: PARAMETER ( FOUR = 4.0D0 )
88: * ..
89: * .. Local Scalars ..
90: LOGICAL GASMAL, SWAP
91: INTEGER PMAX
92: DOUBLE PRECISION A, CLT, CRT, D, FA, FT, GA, GT, HA, HT, L, M,
93: $ MM, R, S, SLT, SRT, T, TEMP, TSIGN, TT
94: * ..
95: * .. Intrinsic Functions ..
96: INTRINSIC ABS, SIGN, SQRT
97: * ..
98: * .. External Functions ..
99: DOUBLE PRECISION DLAMCH
100: EXTERNAL DLAMCH
101: * ..
102: * .. Executable Statements ..
103: *
104: FT = F
105: FA = ABS( FT )
106: HT = H
107: HA = ABS( H )
108: *
109: * PMAX points to the maximum absolute element of matrix
110: * PMAX = 1 if F largest in absolute values
111: * PMAX = 2 if G largest in absolute values
112: * PMAX = 3 if H largest in absolute values
113: *
114: PMAX = 1
115: SWAP = ( HA.GT.FA )
116: IF( SWAP ) THEN
117: PMAX = 3
118: TEMP = FT
119: FT = HT
120: HT = TEMP
121: TEMP = FA
122: FA = HA
123: HA = TEMP
124: *
125: * Now FA .ge. HA
126: *
127: END IF
128: GT = G
129: GA = ABS( GT )
130: IF( GA.EQ.ZERO ) THEN
131: *
132: * Diagonal matrix
133: *
134: SSMIN = HA
135: SSMAX = FA
136: CLT = ONE
137: CRT = ONE
138: SLT = ZERO
139: SRT = ZERO
140: ELSE
141: GASMAL = .TRUE.
142: IF( GA.GT.FA ) THEN
143: PMAX = 2
144: IF( ( FA / GA ).LT.DLAMCH( 'EPS' ) ) THEN
145: *
146: * Case of very large GA
147: *
148: GASMAL = .FALSE.
149: SSMAX = GA
150: IF( HA.GT.ONE ) THEN
151: SSMIN = FA / ( GA / HA )
152: ELSE
153: SSMIN = ( FA / GA )*HA
154: END IF
155: CLT = ONE
156: SLT = HT / GT
157: SRT = ONE
158: CRT = FT / GT
159: END IF
160: END IF
161: IF( GASMAL ) THEN
162: *
163: * Normal case
164: *
165: D = FA - HA
166: IF( D.EQ.FA ) THEN
167: *
168: * Copes with infinite F or H
169: *
170: L = ONE
171: ELSE
172: L = D / FA
173: END IF
174: *
175: * Note that 0 .le. L .le. 1
176: *
177: M = GT / FT
178: *
179: * Note that abs(M) .le. 1/macheps
180: *
181: T = TWO - L
182: *
183: * Note that T .ge. 1
184: *
185: MM = M*M
186: TT = T*T
187: S = SQRT( TT+MM )
188: *
189: * Note that 1 .le. S .le. 1 + 1/macheps
190: *
191: IF( L.EQ.ZERO ) THEN
192: R = ABS( M )
193: ELSE
194: R = SQRT( L*L+MM )
195: END IF
196: *
197: * Note that 0 .le. R .le. 1 + 1/macheps
198: *
199: A = HALF*( S+R )
200: *
201: * Note that 1 .le. A .le. 1 + abs(M)
202: *
203: SSMIN = HA / A
204: SSMAX = FA*A
205: IF( MM.EQ.ZERO ) THEN
206: *
207: * Note that M is very tiny
208: *
209: IF( L.EQ.ZERO ) THEN
210: T = SIGN( TWO, FT )*SIGN( ONE, GT )
211: ELSE
212: T = GT / SIGN( D, FT ) + M / T
213: END IF
214: ELSE
215: T = ( M / ( S+T )+M / ( R+L ) )*( ONE+A )
216: END IF
217: L = SQRT( T*T+FOUR )
218: CRT = TWO / L
219: SRT = T / L
220: CLT = ( CRT+SRT*M ) / A
221: SLT = ( HT / FT )*SRT / A
222: END IF
223: END IF
224: IF( SWAP ) THEN
225: CSL = SRT
226: SNL = CRT
227: CSR = SLT
228: SNR = CLT
229: ELSE
230: CSL = CLT
231: SNL = SLT
232: CSR = CRT
233: SNR = SRT
234: END IF
235: *
236: * Correct signs of SSMAX and SSMIN
237: *
238: IF( PMAX.EQ.1 )
239: $ TSIGN = SIGN( ONE, CSR )*SIGN( ONE, CSL )*SIGN( ONE, F )
240: IF( PMAX.EQ.2 )
241: $ TSIGN = SIGN( ONE, SNR )*SIGN( ONE, CSL )*SIGN( ONE, G )
242: IF( PMAX.EQ.3 )
243: $ TSIGN = SIGN( ONE, SNR )*SIGN( ONE, SNL )*SIGN( ONE, H )
244: SSMAX = SIGN( SSMAX, TSIGN )
245: SSMIN = SIGN( SSMIN, TSIGN*SIGN( ONE, F )*SIGN( ONE, H ) )
246: RETURN
247: *
248: * End of DLASV2
249: *
250: END