Mercurial > octave

comparison libcruft/lapack/slasv2.f @ 7789:82be108cc558

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