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