|
2329
|
1 SUBROUTINE DLASQ1( N, D, E, WORK, INFO ) |
|
|
2 * |
|
|
3 * -- LAPACK routine (version 2.0) -- |
|
|
4 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., |
|
|
5 * Courant Institute, Argonne National Lab, and Rice University |
|
|
6 * September 30, 1994 |
|
|
7 * |
|
|
8 * .. Scalar Arguments .. |
|
|
9 INTEGER INFO, N |
|
|
10 * .. |
|
|
11 * .. Array Arguments .. |
|
|
12 DOUBLE PRECISION D( * ), E( * ), WORK( * ) |
|
|
13 * .. |
|
|
14 * |
|
|
15 * Purpose |
|
|
16 * ======= |
|
|
17 * |
|
|
18 * DLASQ1 computes the singular values of a real N-by-N bidiagonal |
|
|
19 * matrix with diagonal D and off-diagonal E. The singular values are |
|
|
20 * computed to high relative accuracy, barring over/underflow or |
|
|
21 * denormalization. The algorithm is described in |
|
|
22 * |
|
|
23 * "Accurate singular values and differential qd algorithms," by |
|
|
24 * K. V. Fernando and B. N. Parlett, |
|
|
25 * Numer. Math., Vol-67, No. 2, pp. 191-230,1994. |
|
|
26 * |
|
|
27 * See also |
|
|
28 * "Implementation of differential qd algorithms," by |
|
|
29 * K. V. Fernando and B. N. Parlett, Technical Report, |
|
|
30 * Department of Mathematics, University of California at Berkeley, |
|
|
31 * 1994 (Under preparation). |
|
|
32 * |
|
|
33 * Arguments |
|
|
34 * ========= |
|
|
35 * |
|
|
36 * N (input) INTEGER |
|
|
37 * The number of rows and columns in the matrix. N >= 0. |
|
|
38 * |
|
|
39 * D (input/output) DOUBLE PRECISION array, dimension (N) |
|
|
40 * On entry, D contains the diagonal elements of the |
|
|
41 * bidiagonal matrix whose SVD is desired. On normal exit, |
|
|
42 * D contains the singular values in decreasing order. |
|
|
43 * |
|
|
44 * E (input/output) DOUBLE PRECISION array, dimension (N) |
|
|
45 * On entry, elements E(1:N-1) contain the off-diagonal elements |
|
|
46 * of the bidiagonal matrix whose SVD is desired. |
|
|
47 * On exit, E is overwritten. |
|
|
48 * |
|
|
49 * WORK (workspace) DOUBLE PRECISION array, dimension (2*N) |
|
|
50 * |
|
|
51 * INFO (output) INTEGER |
|
|
52 * = 0: successful exit |
|
|
53 * < 0: if INFO = -i, the i-th argument had an illegal value |
|
|
54 * > 0: if INFO = i, the algorithm did not converge; i |
|
|
55 * specifies how many superdiagonals did not converge. |
|
|
56 * |
|
|
57 * ===================================================================== |
|
|
58 * |
|
|
59 * .. Parameters .. |
|
|
60 DOUBLE PRECISION MEIGTH |
|
|
61 PARAMETER ( MEIGTH = -0.125D0 ) |
|
|
62 DOUBLE PRECISION ZERO |
|
|
63 PARAMETER ( ZERO = 0.0D0 ) |
|
|
64 DOUBLE PRECISION ONE |
|
|
65 PARAMETER ( ONE = 1.0D0 ) |
|
|
66 DOUBLE PRECISION TEN |
|
|
67 PARAMETER ( TEN = 10.0D0 ) |
|
|
68 DOUBLE PRECISION HUNDRD |
|
|
69 PARAMETER ( HUNDRD = 100.0D0 ) |
|
|
70 DOUBLE PRECISION TWO56 |
|
|
71 PARAMETER ( TWO56 = 256.0D0 ) |
|
|
72 * .. |
|
|
73 * .. Local Scalars .. |
|
|
74 LOGICAL RESTRT |
|
|
75 INTEGER I, IERR, J, KE, KEND, M, NY |
|
|
76 DOUBLE PRECISION DM, DX, EPS, SCL, SFMIN, SIG1, SIG2, SIGMN, |
|
|
77 $ SIGMX, SMALL2, THRESH, TOL, TOL2, TOLMUL |
|
|
78 * .. |
|
|
79 * .. External Functions .. |
|
|
80 DOUBLE PRECISION DLAMCH |
|
|
81 EXTERNAL DLAMCH |
|
|
82 * .. |
|
|
83 * .. External Subroutines .. |
|
|
84 EXTERNAL DCOPY, DLAS2, DLASCL, DLASQ2, DLASRT, XERBLA |
|
|
85 * .. |
|
|
86 * .. Intrinsic Functions .. |
|
|
87 INTRINSIC ABS, DBLE, MAX, MIN, SQRT |
|
|
88 * .. |
|
|
89 * .. Executable Statements .. |
|
|
90 INFO = 0 |
|
|
91 IF( N.LT.0 ) THEN |
|
|
92 INFO = -2 |
|
|
93 CALL XERBLA( 'DLASQ1', -INFO ) |
|
|
94 RETURN |
|
|
95 ELSE IF( N.EQ.0 ) THEN |
|
|
96 RETURN |
|
|
97 ELSE IF( N.EQ.1 ) THEN |
|
|
98 D( 1 ) = ABS( D( 1 ) ) |
|
|
99 RETURN |
|
|
100 ELSE IF( N.EQ.2 ) THEN |
|
|
101 CALL DLAS2( D( 1 ), E( 1 ), D( 2 ), SIGMN, SIGMX ) |
|
|
102 D( 1 ) = SIGMX |
|
|
103 D( 2 ) = SIGMN |
|
|
104 RETURN |
|
|
105 END IF |
|
|
106 * |
|
|
107 * Estimate the largest singular value |
|
|
108 * |
|
|
109 SIGMX = ZERO |
|
|
110 DO 10 I = 1, N - 1 |
|
|
111 SIGMX = MAX( SIGMX, ABS( E( I ) ) ) |
|
|
112 10 CONTINUE |
|
|
113 * |
|
|
114 * Early return if sigmx is zero (matrix is already diagonal) |
|
|
115 * |
|
|
116 IF( SIGMX.EQ.ZERO ) |
|
|
117 $ GO TO 70 |
|
|
118 * |
|
|
119 DO 20 I = 1, N |
|
|
120 D( I ) = ABS( D( I ) ) |
|
|
121 SIGMX = MAX( SIGMX, D( I ) ) |
|
|
122 20 CONTINUE |
|
|
123 * |
|
|
124 * Get machine parameters |
|
|
125 * |
|
|
126 EPS = DLAMCH( 'EPSILON' ) |
|
|
127 SFMIN = DLAMCH( 'SAFE MINIMUM' ) |
|
|
128 * |
|
|
129 * Compute singular values to relative accuracy TOL |
|
|
130 * It is assumed that tol**2 does not underflow. |
|
|
131 * |
|
|
132 TOLMUL = MAX( TEN, MIN( HUNDRD, EPS**( -MEIGTH ) ) ) |
|
|
133 TOL = TOLMUL*EPS |
|
|
134 TOL2 = TOL**2 |
|
|
135 * |
|
|
136 THRESH = SIGMX*SQRT( SFMIN )*TOL |
|
|
137 * |
|
|
138 * Scale matrix so the square of the largest element is |
|
|
139 * 1 / ( 256 * SFMIN ) |
|
|
140 * |
|
|
141 SCL = SQRT( ONE / ( TWO56*SFMIN ) ) |
|
|
142 SMALL2 = ONE / ( TWO56*TOLMUL**2 ) |
|
|
143 CALL DCOPY( N, D, 1, WORK( 1 ), 1 ) |
|
|
144 CALL DCOPY( N-1, E, 1, WORK( N+1 ), 1 ) |
|
|
145 CALL DLASCL( 'G', 0, 0, SIGMX, SCL, N, 1, WORK( 1 ), N, IERR ) |
|
|
146 CALL DLASCL( 'G', 0, 0, SIGMX, SCL, N-1, 1, WORK( N+1 ), N-1, |
|
|
147 $ IERR ) |
|
|
148 * |
|
|
149 * Square D and E (the input for the qd algorithm) |
|
|
150 * |
|
|
151 DO 30 J = 1, 2*N - 1 |
|
|
152 WORK( J ) = WORK( J )**2 |
|
|
153 30 CONTINUE |
|
|
154 * |
|
|
155 * Apply qd algorithm |
|
|
156 * |
|
|
157 M = 0 |
|
|
158 E( N ) = ZERO |
|
|
159 DX = WORK( 1 ) |
|
|
160 DM = DX |
|
|
161 KE = 0 |
|
|
162 RESTRT = .FALSE. |
|
|
163 DO 60 I = 1, N |
|
|
164 IF( ABS( E( I ) ).LE.THRESH .OR. WORK( N+I ).LE.TOL2* |
|
|
165 $ ( DM / DBLE( I-M ) ) ) THEN |
|
|
166 NY = I - M |
|
|
167 IF( NY.EQ.1 ) THEN |
|
|
168 GO TO 50 |
|
|
169 ELSE IF( NY.EQ.2 ) THEN |
|
|
170 CALL DLAS2( D( M+1 ), E( M+1 ), D( M+2 ), SIG1, SIG2 ) |
|
|
171 D( M+1 ) = SIG1 |
|
|
172 D( M+2 ) = SIG2 |
|
|
173 ELSE |
|
|
174 KEND = KE + 1 - M |
|
|
175 CALL DLASQ2( NY, D( M+1 ), E( M+1 ), WORK( M+1 ), |
|
|
176 $ WORK( M+N+1 ), EPS, TOL2, SMALL2, DM, KEND, |
|
|
177 $ INFO ) |
|
|
178 * |
|
|
179 * Return, INFO = number of unconverged superdiagonals |
|
|
180 * |
|
|
181 IF( INFO.NE.0 ) THEN |
|
|
182 INFO = INFO + I |
|
|
183 RETURN |
|
|
184 END IF |
|
|
185 * |
|
|
186 * Undo scaling |
|
|
187 * |
|
|
188 DO 40 J = M + 1, M + NY |
|
|
189 D( J ) = SQRT( D( J ) ) |
|
|
190 40 CONTINUE |
|
|
191 CALL DLASCL( 'G', 0, 0, SCL, SIGMX, NY, 1, D( M+1 ), NY, |
|
|
192 $ IERR ) |
|
|
193 END IF |
|
|
194 50 CONTINUE |
|
|
195 M = I |
|
|
196 IF( I.NE.N ) THEN |
|
|
197 DX = WORK( I+1 ) |
|
|
198 DM = DX |
|
|
199 KE = I |
|
|
200 RESTRT = .TRUE. |
|
|
201 END IF |
|
|
202 END IF |
|
|
203 IF( I.NE.N .AND. .NOT.RESTRT ) THEN |
|
|
204 DX = WORK( I+1 )*( DX / ( DX+WORK( N+I ) ) ) |
|
|
205 IF( DM.GT.DX ) THEN |
|
|
206 DM = DX |
|
|
207 KE = I |
|
|
208 END IF |
|
|
209 END IF |
|
|
210 RESTRT = .FALSE. |
|
|
211 60 CONTINUE |
|
|
212 KEND = KE + 1 |
|
|
213 * |
|
|
214 * Sort the singular values into decreasing order |
|
|
215 * |
|
|
216 70 CONTINUE |
|
|
217 CALL DLASRT( 'D', N, D, INFO ) |
|
|
218 RETURN |
|
|
219 * |
|
|
220 * End of DLASQ1 |
|
|
221 * |
|
|
222 END |
