← Back to branch summary

~ubuntu-branches/ubuntu/karmic/scilab/karmic

« back to all changes in this revision

Viewing changes to routines/lapack/dlasv2.f

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2002-03-21 16:57:43 UTC
  • Revision ID: james.westby@ubuntu.com-20020321165743-e9mv12c1tb1plztg
Tags: upstream-2.6
Import upstream version 2.6

Show diffs side-by-side

added added

removed removed

Lines of Context:
routines/lapack/dlasv2.f
 
1
      SUBROUTINE DLASV2( F, G, H, SSMIN, SSMAX, SNR, CSR, SNL, CSL )
 
2
*
 
3
*  -- LAPACK auxiliary routine (version 2.0) --
 
4
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
 
5
*     Courant Institute, Argonne National Lab, and Rice University
 
6
*     October 31, 1992
 
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