|
7072
|
1 SUBROUTINE DLASDQ( UPLO, SQRE, N, NCVT, NRU, NCC, D, E, VT, LDVT, |
|
|
2 $ U, LDU, C, LDC, WORK, INFO ) |
|
|
3 * |
|
|
4 * -- LAPACK auxiliary routine (version 3.1) -- |
|
|
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. |
|
|
6 * November 2006 |
|
|
7 * |
|
|
8 * .. Scalar Arguments .. |
|
|
9 CHARACTER UPLO |
|
|
10 INTEGER INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU, SQRE |
|
|
11 * .. |
|
|
12 * .. Array Arguments .. |
|
|
13 DOUBLE PRECISION C( LDC, * ), D( * ), E( * ), U( LDU, * ), |
|
|
14 $ VT( LDVT, * ), WORK( * ) |
|
|
15 * .. |
|
|
16 * |
|
|
17 * Purpose |
|
|
18 * ======= |
|
|
19 * |
|
|
20 * DLASDQ computes the singular value decomposition (SVD) of a real |
|
|
21 * (upper or lower) bidiagonal matrix with diagonal D and offdiagonal |
|
|
22 * E, accumulating the transformations if desired. Letting B denote |
|
|
23 * the input bidiagonal matrix, the algorithm computes orthogonal |
|
|
24 * matrices Q and P such that B = Q * S * P' (P' denotes the transpose |
|
|
25 * of P). The singular values S are overwritten on D. |
|
|
26 * |
|
|
27 * The input matrix U is changed to U * Q if desired. |
|
|
28 * The input matrix VT is changed to P' * VT if desired. |
|
|
29 * The input matrix C is changed to Q' * C if desired. |
|
|
30 * |
|
|
31 * See "Computing Small Singular Values of Bidiagonal Matrices With |
|
|
32 * Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan, |
|
|
33 * LAPACK Working Note #3, for a detailed description of the algorithm. |
|
|
34 * |
|
|
35 * Arguments |
|
|
36 * ========= |
|
|
37 * |
|
|
38 * UPLO (input) CHARACTER*1 |
|
|
39 * On entry, UPLO specifies whether the input bidiagonal matrix |
|
|
40 * is upper or lower bidiagonal, and wether it is square are |
|
|
41 * not. |
|
|
42 * UPLO = 'U' or 'u' B is upper bidiagonal. |
|
|
43 * UPLO = 'L' or 'l' B is lower bidiagonal. |
|
|
44 * |
|
|
45 * SQRE (input) INTEGER |
|
|
46 * = 0: then the input matrix is N-by-N. |
|
|
47 * = 1: then the input matrix is N-by-(N+1) if UPLU = 'U' and |
|
|
48 * (N+1)-by-N if UPLU = 'L'. |
|
|
49 * |
|
|
50 * The bidiagonal matrix has |
|
|
51 * N = NL + NR + 1 rows and |
|
|
52 * M = N + SQRE >= N columns. |
|
|
53 * |
|
|
54 * N (input) INTEGER |
|
|
55 * On entry, N specifies the number of rows and columns |
|
|
56 * in the matrix. N must be at least 0. |
|
|
57 * |
|
|
58 * NCVT (input) INTEGER |
|
|
59 * On entry, NCVT specifies the number of columns of |
|
|
60 * the matrix VT. NCVT must be at least 0. |
|
|
61 * |
|
|
62 * NRU (input) INTEGER |
|
|
63 * On entry, NRU specifies the number of rows of |
|
|
64 * the matrix U. NRU must be at least 0. |
|
|
65 * |
|
|
66 * NCC (input) INTEGER |
|
|
67 * On entry, NCC specifies the number of columns of |
|
|
68 * the matrix C. NCC must be at least 0. |
|
|
69 * |
|
|
70 * D (input/output) DOUBLE PRECISION array, dimension (N) |
|
|
71 * On entry, D contains the diagonal entries of the |
|
|
72 * bidiagonal matrix whose SVD is desired. On normal exit, |
|
|
73 * D contains the singular values in ascending order. |
|
|
74 * |
|
|
75 * E (input/output) DOUBLE PRECISION array. |
|
|
76 * dimension is (N-1) if SQRE = 0 and N if SQRE = 1. |
|
|
77 * On entry, the entries of E contain the offdiagonal entries |
|
|
78 * of the bidiagonal matrix whose SVD is desired. On normal |
|
|
79 * exit, E will contain 0. If the algorithm does not converge, |
|
|
80 * D and E will contain the diagonal and superdiagonal entries |
|
|
81 * of a bidiagonal matrix orthogonally equivalent to the one |
|
|
82 * given as input. |
|
|
83 * |
|
|
84 * VT (input/output) DOUBLE PRECISION array, dimension (LDVT, NCVT) |
|
|
85 * On entry, contains a matrix which on exit has been |
|
|
86 * premultiplied by P', dimension N-by-NCVT if SQRE = 0 |
|
|
87 * and (N+1)-by-NCVT if SQRE = 1 (not referenced if NCVT=0). |
|
|
88 * |
|
|
89 * LDVT (input) INTEGER |
|
|
90 * On entry, LDVT specifies the leading dimension of VT as |
|
|
91 * declared in the calling (sub) program. LDVT must be at |
|
|
92 * least 1. If NCVT is nonzero LDVT must also be at least N. |
|
|
93 * |
|
|
94 * U (input/output) DOUBLE PRECISION array, dimension (LDU, N) |
|
|
95 * On entry, contains a matrix which on exit has been |
|
|
96 * postmultiplied by Q, dimension NRU-by-N if SQRE = 0 |
|
|
97 * and NRU-by-(N+1) if SQRE = 1 (not referenced if NRU=0). |
|
|
98 * |
|
|
99 * LDU (input) INTEGER |
|
|
100 * On entry, LDU specifies the leading dimension of U as |
|
|
101 * declared in the calling (sub) program. LDU must be at |
|
|
102 * least max( 1, NRU ) . |
|
|
103 * |
|
|
104 * C (input/output) DOUBLE PRECISION array, dimension (LDC, NCC) |
|
|
105 * On entry, contains an N-by-NCC matrix which on exit |
|
|
106 * has been premultiplied by Q' dimension N-by-NCC if SQRE = 0 |
|
|
107 * and (N+1)-by-NCC if SQRE = 1 (not referenced if NCC=0). |
|
|
108 * |
|
|
109 * LDC (input) INTEGER |
|
|
110 * On entry, LDC specifies the leading dimension of C as |
|
|
111 * declared in the calling (sub) program. LDC must be at |
|
|
112 * least 1. If NCC is nonzero, LDC must also be at least N. |
|
|
113 * |
|
|
114 * WORK (workspace) DOUBLE PRECISION array, dimension (4*N) |
|
|
115 * Workspace. Only referenced if one of NCVT, NRU, or NCC is |
|
|
116 * nonzero, and if N is at least 2. |
|
|
117 * |
|
|
118 * INFO (output) INTEGER |
|
|
119 * On exit, a value of 0 indicates a successful exit. |
|
|
120 * If INFO < 0, argument number -INFO is illegal. |
|
|
121 * If INFO > 0, the algorithm did not converge, and INFO |
|
|
122 * specifies how many superdiagonals did not converge. |
|
|
123 * |
|
|
124 * Further Details |
|
|
125 * =============== |
|
|
126 * |
|
|
127 * Based on contributions by |
|
|
128 * Ming Gu and Huan Ren, Computer Science Division, University of |
|
|
129 * California at Berkeley, USA |
|
|
130 * |
|
|
131 * ===================================================================== |
|
|
132 * |
|
|
133 * .. Parameters .. |
|
|
134 DOUBLE PRECISION ZERO |
|
|
135 PARAMETER ( ZERO = 0.0D+0 ) |
|
|
136 * .. |
|
|
137 * .. Local Scalars .. |
|
|
138 LOGICAL ROTATE |
|
|
139 INTEGER I, ISUB, IUPLO, J, NP1, SQRE1 |
|
|
140 DOUBLE PRECISION CS, R, SMIN, SN |
|
|
141 * .. |
|
|
142 * .. External Subroutines .. |
|
|
143 EXTERNAL DBDSQR, DLARTG, DLASR, DSWAP, XERBLA |
|
|
144 * .. |
|
|
145 * .. External Functions .. |
|
|
146 LOGICAL LSAME |
|
|
147 EXTERNAL LSAME |
|
|
148 * .. |
|
|
149 * .. Intrinsic Functions .. |
|
|
150 INTRINSIC MAX |
|
|
151 * .. |
|
|
152 * .. Executable Statements .. |
|
|
153 * |
|
|
154 * Test the input parameters. |
|
|
155 * |
|
|
156 INFO = 0 |
|
|
157 IUPLO = 0 |
|
|
158 IF( LSAME( UPLO, 'U' ) ) |
|
|
159 $ IUPLO = 1 |
|
|
160 IF( LSAME( UPLO, 'L' ) ) |
|
|
161 $ IUPLO = 2 |
|
|
162 IF( IUPLO.EQ.0 ) THEN |
|
|
163 INFO = -1 |
|
|
164 ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN |
|
|
165 INFO = -2 |
|
|
166 ELSE IF( N.LT.0 ) THEN |
|
|
167 INFO = -3 |
|
|
168 ELSE IF( NCVT.LT.0 ) THEN |
|
|
169 INFO = -4 |
|
|
170 ELSE IF( NRU.LT.0 ) THEN |
|
|
171 INFO = -5 |
|
|
172 ELSE IF( NCC.LT.0 ) THEN |
|
|
173 INFO = -6 |
|
|
174 ELSE IF( ( NCVT.EQ.0 .AND. LDVT.LT.1 ) .OR. |
|
|
175 $ ( NCVT.GT.0 .AND. LDVT.LT.MAX( 1, N ) ) ) THEN |
|
|
176 INFO = -10 |
|
|
177 ELSE IF( LDU.LT.MAX( 1, NRU ) ) THEN |
|
|
178 INFO = -12 |
|
|
179 ELSE IF( ( NCC.EQ.0 .AND. LDC.LT.1 ) .OR. |
|
|
180 $ ( NCC.GT.0 .AND. LDC.LT.MAX( 1, N ) ) ) THEN |
|
|
181 INFO = -14 |
|
|
182 END IF |
|
|
183 IF( INFO.NE.0 ) THEN |
|
|
184 CALL XERBLA( 'DLASDQ', -INFO ) |
|
|
185 RETURN |
|
|
186 END IF |
|
|
187 IF( N.EQ.0 ) |
|
|
188 $ RETURN |
|
|
189 * |
|
|
190 * ROTATE is true if any singular vectors desired, false otherwise |
|
|
191 * |
|
|
192 ROTATE = ( NCVT.GT.0 ) .OR. ( NRU.GT.0 ) .OR. ( NCC.GT.0 ) |
|
|
193 NP1 = N + 1 |
|
|
194 SQRE1 = SQRE |
|
|
195 * |
|
|
196 * If matrix non-square upper bidiagonal, rotate to be lower |
|
|
197 * bidiagonal. The rotations are on the right. |
|
|
198 * |
|
|
199 IF( ( IUPLO.EQ.1 ) .AND. ( SQRE1.EQ.1 ) ) THEN |
|
|
200 DO 10 I = 1, N - 1 |
|
|
201 CALL DLARTG( D( I ), E( I ), CS, SN, R ) |
|
|
202 D( I ) = R |
|
|
203 E( I ) = SN*D( I+1 ) |
|
|
204 D( I+1 ) = CS*D( I+1 ) |
|
|
205 IF( ROTATE ) THEN |
|
|
206 WORK( I ) = CS |
|
|
207 WORK( N+I ) = SN |
|
|
208 END IF |
|
|
209 10 CONTINUE |
|
|
210 CALL DLARTG( D( N ), E( N ), CS, SN, R ) |
|
|
211 D( N ) = R |
|
|
212 E( N ) = ZERO |
|
|
213 IF( ROTATE ) THEN |
|
|
214 WORK( N ) = CS |
|
|
215 WORK( N+N ) = SN |
|
|
216 END IF |
|
|
217 IUPLO = 2 |
|
|
218 SQRE1 = 0 |
|
|
219 * |
|
|
220 * Update singular vectors if desired. |
|
|
221 * |
|
|
222 IF( NCVT.GT.0 ) |
|
|
223 $ CALL DLASR( 'L', 'V', 'F', NP1, NCVT, WORK( 1 ), |
|
|
224 $ WORK( NP1 ), VT, LDVT ) |
|
|
225 END IF |
|
|
226 * |
|
|
227 * If matrix lower bidiagonal, rotate to be upper bidiagonal |
|
|
228 * by applying Givens rotations on the left. |
|
|
229 * |
|
|
230 IF( IUPLO.EQ.2 ) THEN |
|
|
231 DO 20 I = 1, N - 1 |
|
|
232 CALL DLARTG( D( I ), E( I ), CS, SN, R ) |
|
|
233 D( I ) = R |
|
|
234 E( I ) = SN*D( I+1 ) |
|
|
235 D( I+1 ) = CS*D( I+1 ) |
|
|
236 IF( ROTATE ) THEN |
|
|
237 WORK( I ) = CS |
|
|
238 WORK( N+I ) = SN |
|
|
239 END IF |
|
|
240 20 CONTINUE |
|
|
241 * |
|
|
242 * If matrix (N+1)-by-N lower bidiagonal, one additional |
|
|
243 * rotation is needed. |
|
|
244 * |
|
|
245 IF( SQRE1.EQ.1 ) THEN |
|
|
246 CALL DLARTG( D( N ), E( N ), CS, SN, R ) |
|
|
247 D( N ) = R |
|
|
248 IF( ROTATE ) THEN |
|
|
249 WORK( N ) = CS |
|
|
250 WORK( N+N ) = SN |
|
|
251 END IF |
|
|
252 END IF |
|
|
253 * |
|
|
254 * Update singular vectors if desired. |
|
|
255 * |
|
|
256 IF( NRU.GT.0 ) THEN |
|
|
257 IF( SQRE1.EQ.0 ) THEN |
|
|
258 CALL DLASR( 'R', 'V', 'F', NRU, N, WORK( 1 ), |
|
|
259 $ WORK( NP1 ), U, LDU ) |
|
|
260 ELSE |
|
|
261 CALL DLASR( 'R', 'V', 'F', NRU, NP1, WORK( 1 ), |
|
|
262 $ WORK( NP1 ), U, LDU ) |
|
|
263 END IF |
|
|
264 END IF |
|
|
265 IF( NCC.GT.0 ) THEN |
|
|
266 IF( SQRE1.EQ.0 ) THEN |
|
|
267 CALL DLASR( 'L', 'V', 'F', N, NCC, WORK( 1 ), |
|
|
268 $ WORK( NP1 ), C, LDC ) |
|
|
269 ELSE |
|
|
270 CALL DLASR( 'L', 'V', 'F', NP1, NCC, WORK( 1 ), |
|
|
271 $ WORK( NP1 ), C, LDC ) |
|
|
272 END IF |
|
|
273 END IF |
|
|
274 END IF |
|
|
275 * |
|
|
276 * Call DBDSQR to compute the SVD of the reduced real |
|
|
277 * N-by-N upper bidiagonal matrix. |
|
|
278 * |
|
|
279 CALL DBDSQR( 'U', N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, |
|
|
280 $ LDC, WORK, INFO ) |
|
|
281 * |
|
|
282 * Sort the singular values into ascending order (insertion sort on |
|
|
283 * singular values, but only one transposition per singular vector) |
|
|
284 * |
|
|
285 DO 40 I = 1, N |
|
|
286 * |
|
|
287 * Scan for smallest D(I). |
|
|
288 * |
|
|
289 ISUB = I |
|
|
290 SMIN = D( I ) |
|
|
291 DO 30 J = I + 1, N |
|
|
292 IF( D( J ).LT.SMIN ) THEN |
|
|
293 ISUB = J |
|
|
294 SMIN = D( J ) |
|
|
295 END IF |
|
|
296 30 CONTINUE |
|
|
297 IF( ISUB.NE.I ) THEN |
|
|
298 * |
|
|
299 * Swap singular values and vectors. |
|
|
300 * |
|
|
301 D( ISUB ) = D( I ) |
|
|
302 D( I ) = SMIN |
|
|
303 IF( NCVT.GT.0 ) |
|
|
304 $ CALL DSWAP( NCVT, VT( ISUB, 1 ), LDVT, VT( I, 1 ), LDVT ) |
|
|
305 IF( NRU.GT.0 ) |
|
|
306 $ CALL DSWAP( NRU, U( 1, ISUB ), 1, U( 1, I ), 1 ) |
|
|
307 IF( NCC.GT.0 ) |
|
|
308 $ CALL DSWAP( NCC, C( ISUB, 1 ), LDC, C( I, 1 ), LDC ) |
|
|
309 END IF |
|
|
310 40 CONTINUE |
|
|
311 * |
|
|
312 RETURN |
|
|
313 * |
|
|
314 * End of DLASDQ |
|
|
315 * |
|
|
316 END |
