1 SUBROUTINE DBDSDC( UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ,
2 $ WORK, IWORK, INFO )
3 *
4 * -- LAPACK routine (version 3.3.1) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * -- April 2011 --
8 *
9 * .. Scalar Arguments ..
10 CHARACTER COMPQ, UPLO
11 INTEGER INFO, LDU, LDVT, N
12 * ..
13 * .. Array Arguments ..
14 INTEGER IQ( * ), IWORK( * )
15 DOUBLE PRECISION D( * ), E( * ), Q( * ), U( LDU, * ),
16 $ VT( LDVT, * ), WORK( * )
17 * ..
18 *
19 * Purpose
20 * =======
21 *
22 * DBDSDC computes the singular value decomposition (SVD) of a real
23 * N-by-N (upper or lower) bidiagonal matrix B: B = U * S * VT,
24 * using a divide and conquer method, where S is a diagonal matrix
25 * with non-negative diagonal elements (the singular values of B), and
26 * U and VT are orthogonal matrices of left and right singular vectors,
27 * respectively. DBDSDC can be used to compute all singular values,
28 * and optionally, singular vectors or singular vectors in compact form.
29 *
30 * This code makes very mild assumptions about floating point
31 * arithmetic. It will work on machines with a guard digit in
32 * add/subtract, or on those binary machines without guard digits
33 * which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
34 * It could conceivably fail on hexadecimal or decimal machines
35 * without guard digits, but we know of none. See DLASD3 for details.
36 *
37 * The code currently calls DLASDQ if singular values only are desired.
38 * However, it can be slightly modified to compute singular values
39 * using the divide and conquer method.
40 *
41 * Arguments
42 * =========
43 *
44 * UPLO (input) CHARACTER*1
45 * = 'U': B is upper bidiagonal.
46 * = 'L': B is lower bidiagonal.
47 *
48 * COMPQ (input) CHARACTER*1
49 * Specifies whether singular vectors are to be computed
50 * as follows:
51 * = 'N': Compute singular values only;
52 * = 'P': Compute singular values and compute singular
53 * vectors in compact form;
54 * = 'I': Compute singular values and singular vectors.
55 *
56 * N (input) INTEGER
57 * The order of the matrix B. N >= 0.
58 *
59 * D (input/output) DOUBLE PRECISION array, dimension (N)
60 * On entry, the n diagonal elements of the bidiagonal matrix B.
61 * On exit, if INFO=0, the singular values of B.
62 *
63 * E (input/output) DOUBLE PRECISION array, dimension (N-1)
64 * On entry, the elements of E contain the offdiagonal
65 * elements of the bidiagonal matrix whose SVD is desired.
66 * On exit, E has been destroyed.
67 *
68 * U (output) DOUBLE PRECISION array, dimension (LDU,N)
69 * If COMPQ = 'I', then:
70 * On exit, if INFO = 0, U contains the left singular vectors
71 * of the bidiagonal matrix.
72 * For other values of COMPQ, U is not referenced.
73 *
74 * LDU (input) INTEGER
75 * The leading dimension of the array U. LDU >= 1.
76 * If singular vectors are desired, then LDU >= max( 1, N ).
77 *
78 * VT (output) DOUBLE PRECISION array, dimension (LDVT,N)
79 * If COMPQ = 'I', then:
80 * On exit, if INFO = 0, VT**T contains the right singular
81 * vectors of the bidiagonal matrix.
82 * For other values of COMPQ, VT is not referenced.
83 *
84 * LDVT (input) INTEGER
85 * The leading dimension of the array VT. LDVT >= 1.
86 * If singular vectors are desired, then LDVT >= max( 1, N ).
87 *
88 * Q (output) DOUBLE PRECISION array, dimension (LDQ)
89 * If COMPQ = 'P', then:
90 * On exit, if INFO = 0, Q and IQ contain the left
91 * and right singular vectors in a compact form,
92 * requiring O(N log N) space instead of 2*N**2.
93 * In particular, Q contains all the DOUBLE PRECISION data in
94 * LDQ >= N*(11 + 2*SMLSIZ + 8*INT(LOG_2(N/(SMLSIZ+1))))
95 * words of memory, where SMLSIZ is returned by ILAENV and
96 * is equal to the maximum size of the subproblems at the
97 * bottom of the computation tree (usually about 25).
98 * For other values of COMPQ, Q is not referenced.
99 *
100 * IQ (output) INTEGER array, dimension (LDIQ)
101 * If COMPQ = 'P', then:
102 * On exit, if INFO = 0, Q and IQ contain the left
103 * and right singular vectors in a compact form,
104 * requiring O(N log N) space instead of 2*N**2.
105 * In particular, IQ contains all INTEGER data in
106 * LDIQ >= N*(3 + 3*INT(LOG_2(N/(SMLSIZ+1))))
107 * words of memory, where SMLSIZ is returned by ILAENV and
108 * is equal to the maximum size of the subproblems at the
109 * bottom of the computation tree (usually about 25).
110 * For other values of COMPQ, IQ is not referenced.
111 *
112 * WORK (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
113 * If COMPQ = 'N' then LWORK >= (4 * N).
114 * If COMPQ = 'P' then LWORK >= (6 * N).
115 * If COMPQ = 'I' then LWORK >= (3 * N**2 + 4 * N).
116 *
117 * IWORK (workspace) INTEGER array, dimension (8*N)
118 *
119 * INFO (output) INTEGER
120 * = 0: successful exit.
121 * < 0: if INFO = -i, the i-th argument had an illegal value.
122 * > 0: The algorithm failed to compute a singular value.
123 * The update process of divide and conquer failed.
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 * Changed dimension statement in comment describing E from (N) to
134 * (N-1). Sven, 17 Feb 05.
135 * =====================================================================
136 *
137 * .. Parameters ..
138 DOUBLE PRECISION ZERO, ONE, TWO
139 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 )
140 * ..
141 * .. Local Scalars ..
142 INTEGER DIFL, DIFR, GIVCOL, GIVNUM, GIVPTR, I, IC,
143 $ ICOMPQ, IERR, II, IS, IU, IUPLO, IVT, J, K, KK,
144 $ MLVL, NM1, NSIZE, PERM, POLES, QSTART, SMLSIZ,
145 $ SMLSZP, SQRE, START, WSTART, Z
146 DOUBLE PRECISION CS, EPS, ORGNRM, P, R, SN
147 * ..
148 * .. External Functions ..
149 LOGICAL LSAME
150 INTEGER ILAENV
151 DOUBLE PRECISION DLAMCH, DLANST
152 EXTERNAL LSAME, ILAENV, DLAMCH, DLANST
153 * ..
154 * .. External Subroutines ..
155 EXTERNAL DCOPY, DLARTG, DLASCL, DLASD0, DLASDA, DLASDQ,
156 $ DLASET, DLASR, DSWAP, XERBLA
157 * ..
158 * .. Intrinsic Functions ..
159 INTRINSIC ABS, DBLE, INT, LOG, SIGN
160 * ..
161 * .. Executable Statements ..
162 *
163 * Test the input parameters.
164 *
165 INFO = 0
166 *
167 IUPLO = 0
168 IF( LSAME( UPLO, 'U' ) )
169 $ IUPLO = 1
170 IF( LSAME( UPLO, 'L' ) )
171 $ IUPLO = 2
172 IF( LSAME( COMPQ, 'N' ) ) THEN
173 ICOMPQ = 0
174 ELSE IF( LSAME( COMPQ, 'P' ) ) THEN
175 ICOMPQ = 1
176 ELSE IF( LSAME( COMPQ, 'I' ) ) THEN
177 ICOMPQ = 2
178 ELSE
179 ICOMPQ = -1
180 END IF
181 IF( IUPLO.EQ.0 ) THEN
182 INFO = -1
183 ELSE IF( ICOMPQ.LT.0 ) THEN
184 INFO = -2
185 ELSE IF( N.LT.0 ) THEN
186 INFO = -3
187 ELSE IF( ( LDU.LT.1 ) .OR. ( ( ICOMPQ.EQ.2 ) .AND. ( LDU.LT.
188 $ N ) ) ) THEN
189 INFO = -7
190 ELSE IF( ( LDVT.LT.1 ) .OR. ( ( ICOMPQ.EQ.2 ) .AND. ( LDVT.LT.
191 $ N ) ) ) THEN
192 INFO = -9
193 END IF
194 IF( INFO.NE.0 ) THEN
195 CALL XERBLA( 'DBDSDC', -INFO )
196 RETURN
197 END IF
198 *
199 * Quick return if possible
200 *
201 IF( N.EQ.0 )
202 $ RETURN
203 SMLSIZ = ILAENV( 9, 'DBDSDC', ' ', 0, 0, 0, 0 )
204 IF( N.EQ.1 ) THEN
205 IF( ICOMPQ.EQ.1 ) THEN
206 Q( 1 ) = SIGN( ONE, D( 1 ) )
207 Q( 1+SMLSIZ*N ) = ONE
208 ELSE IF( ICOMPQ.EQ.2 ) THEN
209 U( 1, 1 ) = SIGN( ONE, D( 1 ) )
210 VT( 1, 1 ) = ONE
211 END IF
212 D( 1 ) = ABS( D( 1 ) )
213 RETURN
214 END IF
215 NM1 = N - 1
216 *
217 * If matrix lower bidiagonal, rotate to be upper bidiagonal
218 * by applying Givens rotations on the left
219 *
220 WSTART = 1
221 QSTART = 3
222 IF( ICOMPQ.EQ.1 ) THEN
223 CALL DCOPY( N, D, 1, Q( 1 ), 1 )
224 CALL DCOPY( N-1, E, 1, Q( N+1 ), 1 )
225 END IF
226 IF( IUPLO.EQ.2 ) THEN
227 QSTART = 5
228 WSTART = 2*N - 1
229 DO 10 I = 1, N - 1
230 CALL DLARTG( D( I ), E( I ), CS, SN, R )
231 D( I ) = R
232 E( I ) = SN*D( I+1 )
233 D( I+1 ) = CS*D( I+1 )
234 IF( ICOMPQ.EQ.1 ) THEN
235 Q( I+2*N ) = CS
236 Q( I+3*N ) = SN
237 ELSE IF( ICOMPQ.EQ.2 ) THEN
238 WORK( I ) = CS
239 WORK( NM1+I ) = -SN
240 END IF
241 10 CONTINUE
242 END IF
243 *
244 * If ICOMPQ = 0, use DLASDQ to compute the singular values.
245 *
246 IF( ICOMPQ.EQ.0 ) THEN
247 CALL DLASDQ( 'U', 0, N, 0, 0, 0, D, E, VT, LDVT, U, LDU, U,
248 $ LDU, WORK( WSTART ), INFO )
249 GO TO 40
250 END IF
251 *
252 * If N is smaller than the minimum divide size SMLSIZ, then solve
253 * the problem with another solver.
254 *
255 IF( N.LE.SMLSIZ ) THEN
256 IF( ICOMPQ.EQ.2 ) THEN
257 CALL DLASET( 'A', N, N, ZERO, ONE, U, LDU )
258 CALL DLASET( 'A', N, N, ZERO, ONE, VT, LDVT )
259 CALL DLASDQ( 'U', 0, N, N, N, 0, D, E, VT, LDVT, U, LDU, U,
260 $ LDU, WORK( WSTART ), INFO )
261 ELSE IF( ICOMPQ.EQ.1 ) THEN
262 IU = 1
263 IVT = IU + N
264 CALL DLASET( 'A', N, N, ZERO, ONE, Q( IU+( QSTART-1 )*N ),
265 $ N )
266 CALL DLASET( 'A', N, N, ZERO, ONE, Q( IVT+( QSTART-1 )*N ),
267 $ N )
268 CALL DLASDQ( 'U', 0, N, N, N, 0, D, E,
269 $ Q( IVT+( QSTART-1 )*N ), N,
270 $ Q( IU+( QSTART-1 )*N ), N,
271 $ Q( IU+( QSTART-1 )*N ), N, WORK( WSTART ),
272 $ INFO )
273 END IF
274 GO TO 40
275 END IF
276 *
277 IF( ICOMPQ.EQ.2 ) THEN
278 CALL DLASET( 'A', N, N, ZERO, ONE, U, LDU )
279 CALL DLASET( 'A', N, N, ZERO, ONE, VT, LDVT )
280 END IF
281 *
282 * Scale.
283 *
284 ORGNRM = DLANST( 'M', N, D, E )
285 IF( ORGNRM.EQ.ZERO )
286 $ RETURN
287 CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, IERR )
288 CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, NM1, 1, E, NM1, IERR )
289 *
290 EPS = (0.9D+0)*DLAMCH( 'Epsilon' )
291 *
292 MLVL = INT( LOG( DBLE( N ) / DBLE( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1
293 SMLSZP = SMLSIZ + 1
294 *
295 IF( ICOMPQ.EQ.1 ) THEN
296 IU = 1
297 IVT = 1 + SMLSIZ
298 DIFL = IVT + SMLSZP
299 DIFR = DIFL + MLVL
300 Z = DIFR + MLVL*2
301 IC = Z + MLVL
302 IS = IC + 1
303 POLES = IS + 1
304 GIVNUM = POLES + 2*MLVL
305 *
306 K = 1
307 GIVPTR = 2
308 PERM = 3
309 GIVCOL = PERM + MLVL
310 END IF
311 *
312 DO 20 I = 1, N
313 IF( ABS( D( I ) ).LT.EPS ) THEN
314 D( I ) = SIGN( EPS, D( I ) )
315 END IF
316 20 CONTINUE
317 *
318 START = 1
319 SQRE = 0
320 *
321 DO 30 I = 1, NM1
322 IF( ( ABS( E( I ) ).LT.EPS ) .OR. ( I.EQ.NM1 ) ) THEN
323 *
324 * Subproblem found. First determine its size and then
325 * apply divide and conquer on it.
326 *
327 IF( I.LT.NM1 ) THEN
328 *
329 * A subproblem with E(I) small for I < NM1.
330 *
331 NSIZE = I - START + 1
332 ELSE IF( ABS( E( I ) ).GE.EPS ) THEN
333 *
334 * A subproblem with E(NM1) not too small but I = NM1.
335 *
336 NSIZE = N - START + 1
337 ELSE
338 *
339 * A subproblem with E(NM1) small. This implies an
340 * 1-by-1 subproblem at D(N). Solve this 1-by-1 problem
341 * first.
342 *
343 NSIZE = I - START + 1
344 IF( ICOMPQ.EQ.2 ) THEN
345 U( N, N ) = SIGN( ONE, D( N ) )
346 VT( N, N ) = ONE
347 ELSE IF( ICOMPQ.EQ.1 ) THEN
348 Q( N+( QSTART-1 )*N ) = SIGN( ONE, D( N ) )
349 Q( N+( SMLSIZ+QSTART-1 )*N ) = ONE
350 END IF
351 D( N ) = ABS( D( N ) )
352 END IF
353 IF( ICOMPQ.EQ.2 ) THEN
354 CALL DLASD0( NSIZE, SQRE, D( START ), E( START ),
355 $ U( START, START ), LDU, VT( START, START ),
356 $ LDVT, SMLSIZ, IWORK, WORK( WSTART ), INFO )
357 ELSE
358 CALL DLASDA( ICOMPQ, SMLSIZ, NSIZE, SQRE, D( START ),
359 $ E( START ), Q( START+( IU+QSTART-2 )*N ), N,
360 $ Q( START+( IVT+QSTART-2 )*N ),
361 $ IQ( START+K*N ), Q( START+( DIFL+QSTART-2 )*
362 $ N ), Q( START+( DIFR+QSTART-2 )*N ),
363 $ Q( START+( Z+QSTART-2 )*N ),
364 $ Q( START+( POLES+QSTART-2 )*N ),
365 $ IQ( START+GIVPTR*N ), IQ( START+GIVCOL*N ),
366 $ N, IQ( START+PERM*N ),
367 $ Q( START+( GIVNUM+QSTART-2 )*N ),
368 $ Q( START+( IC+QSTART-2 )*N ),
369 $ Q( START+( IS+QSTART-2 )*N ),
370 $ WORK( WSTART ), IWORK, INFO )
371 END IF
372 IF( INFO.NE.0 ) THEN
373 RETURN
374 END IF
375 START = I + 1
376 END IF
377 30 CONTINUE
378 *
379 * Unscale
380 *
381 CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, IERR )
382 40 CONTINUE
383 *
384 * Use Selection Sort to minimize swaps of singular vectors
385 *
386 DO 60 II = 2, N
387 I = II - 1
388 KK = I
389 P = D( I )
390 DO 50 J = II, N
391 IF( D( J ).GT.P ) THEN
392 KK = J
393 P = D( J )
394 END IF
395 50 CONTINUE
396 IF( KK.NE.I ) THEN
397 D( KK ) = D( I )
398 D( I ) = P
399 IF( ICOMPQ.EQ.1 ) THEN
400 IQ( I ) = KK
401 ELSE IF( ICOMPQ.EQ.2 ) THEN
402 CALL DSWAP( N, U( 1, I ), 1, U( 1, KK ), 1 )
403 CALL DSWAP( N, VT( I, 1 ), LDVT, VT( KK, 1 ), LDVT )
404 END IF
405 ELSE IF( ICOMPQ.EQ.1 ) THEN
406 IQ( I ) = I
407 END IF
408 60 CONTINUE
409 *
410 * If ICOMPQ = 1, use IQ(N,1) as the indicator for UPLO
411 *
412 IF( ICOMPQ.EQ.1 ) THEN
413 IF( IUPLO.EQ.1 ) THEN
414 IQ( N ) = 1
415 ELSE
416 IQ( N ) = 0
417 END IF
418 END IF
419 *
420 * If B is lower bidiagonal, update U by those Givens rotations
421 * which rotated B to be upper bidiagonal
422 *
423 IF( ( IUPLO.EQ.2 ) .AND. ( ICOMPQ.EQ.2 ) )
424 $ CALL DLASR( 'L', 'V', 'B', N, N, WORK( 1 ), WORK( N ), U, LDU )
425 *
426 RETURN
427 *
428 * End of DBDSDC
429 *
430 END
2 $ WORK, IWORK, INFO )
3 *
4 * -- LAPACK routine (version 3.3.1) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * -- April 2011 --
8 *
9 * .. Scalar Arguments ..
10 CHARACTER COMPQ, UPLO
11 INTEGER INFO, LDU, LDVT, N
12 * ..
13 * .. Array Arguments ..
14 INTEGER IQ( * ), IWORK( * )
15 DOUBLE PRECISION D( * ), E( * ), Q( * ), U( LDU, * ),
16 $ VT( LDVT, * ), WORK( * )
17 * ..
18 *
19 * Purpose
20 * =======
21 *
22 * DBDSDC computes the singular value decomposition (SVD) of a real
23 * N-by-N (upper or lower) bidiagonal matrix B: B = U * S * VT,
24 * using a divide and conquer method, where S is a diagonal matrix
25 * with non-negative diagonal elements (the singular values of B), and
26 * U and VT are orthogonal matrices of left and right singular vectors,
27 * respectively. DBDSDC can be used to compute all singular values,
28 * and optionally, singular vectors or singular vectors in compact form.
29 *
30 * This code makes very mild assumptions about floating point
31 * arithmetic. It will work on machines with a guard digit in
32 * add/subtract, or on those binary machines without guard digits
33 * which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
34 * It could conceivably fail on hexadecimal or decimal machines
35 * without guard digits, but we know of none. See DLASD3 for details.
36 *
37 * The code currently calls DLASDQ if singular values only are desired.
38 * However, it can be slightly modified to compute singular values
39 * using the divide and conquer method.
40 *
41 * Arguments
42 * =========
43 *
44 * UPLO (input) CHARACTER*1
45 * = 'U': B is upper bidiagonal.
46 * = 'L': B is lower bidiagonal.
47 *
48 * COMPQ (input) CHARACTER*1
49 * Specifies whether singular vectors are to be computed
50 * as follows:
51 * = 'N': Compute singular values only;
52 * = 'P': Compute singular values and compute singular
53 * vectors in compact form;
54 * = 'I': Compute singular values and singular vectors.
55 *
56 * N (input) INTEGER
57 * The order of the matrix B. N >= 0.
58 *
59 * D (input/output) DOUBLE PRECISION array, dimension (N)
60 * On entry, the n diagonal elements of the bidiagonal matrix B.
61 * On exit, if INFO=0, the singular values of B.
62 *
63 * E (input/output) DOUBLE PRECISION array, dimension (N-1)
64 * On entry, the elements of E contain the offdiagonal
65 * elements of the bidiagonal matrix whose SVD is desired.
66 * On exit, E has been destroyed.
67 *
68 * U (output) DOUBLE PRECISION array, dimension (LDU,N)
69 * If COMPQ = 'I', then:
70 * On exit, if INFO = 0, U contains the left singular vectors
71 * of the bidiagonal matrix.
72 * For other values of COMPQ, U is not referenced.
73 *
74 * LDU (input) INTEGER
75 * The leading dimension of the array U. LDU >= 1.
76 * If singular vectors are desired, then LDU >= max( 1, N ).
77 *
78 * VT (output) DOUBLE PRECISION array, dimension (LDVT,N)
79 * If COMPQ = 'I', then:
80 * On exit, if INFO = 0, VT**T contains the right singular
81 * vectors of the bidiagonal matrix.
82 * For other values of COMPQ, VT is not referenced.
83 *
84 * LDVT (input) INTEGER
85 * The leading dimension of the array VT. LDVT >= 1.
86 * If singular vectors are desired, then LDVT >= max( 1, N ).
87 *
88 * Q (output) DOUBLE PRECISION array, dimension (LDQ)
89 * If COMPQ = 'P', then:
90 * On exit, if INFO = 0, Q and IQ contain the left
91 * and right singular vectors in a compact form,
92 * requiring O(N log N) space instead of 2*N**2.
93 * In particular, Q contains all the DOUBLE PRECISION data in
94 * LDQ >= N*(11 + 2*SMLSIZ + 8*INT(LOG_2(N/(SMLSIZ+1))))
95 * words of memory, where SMLSIZ is returned by ILAENV and
96 * is equal to the maximum size of the subproblems at the
97 * bottom of the computation tree (usually about 25).
98 * For other values of COMPQ, Q is not referenced.
99 *
100 * IQ (output) INTEGER array, dimension (LDIQ)
101 * If COMPQ = 'P', then:
102 * On exit, if INFO = 0, Q and IQ contain the left
103 * and right singular vectors in a compact form,
104 * requiring O(N log N) space instead of 2*N**2.
105 * In particular, IQ contains all INTEGER data in
106 * LDIQ >= N*(3 + 3*INT(LOG_2(N/(SMLSIZ+1))))
107 * words of memory, where SMLSIZ is returned by ILAENV and
108 * is equal to the maximum size of the subproblems at the
109 * bottom of the computation tree (usually about 25).
110 * For other values of COMPQ, IQ is not referenced.
111 *
112 * WORK (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
113 * If COMPQ = 'N' then LWORK >= (4 * N).
114 * If COMPQ = 'P' then LWORK >= (6 * N).
115 * If COMPQ = 'I' then LWORK >= (3 * N**2 + 4 * N).
116 *
117 * IWORK (workspace) INTEGER array, dimension (8*N)
118 *
119 * INFO (output) INTEGER
120 * = 0: successful exit.
121 * < 0: if INFO = -i, the i-th argument had an illegal value.
122 * > 0: The algorithm failed to compute a singular value.
123 * The update process of divide and conquer failed.
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 * Changed dimension statement in comment describing E from (N) to
134 * (N-1). Sven, 17 Feb 05.
135 * =====================================================================
136 *
137 * .. Parameters ..
138 DOUBLE PRECISION ZERO, ONE, TWO
139 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 )
140 * ..
141 * .. Local Scalars ..
142 INTEGER DIFL, DIFR, GIVCOL, GIVNUM, GIVPTR, I, IC,
143 $ ICOMPQ, IERR, II, IS, IU, IUPLO, IVT, J, K, KK,
144 $ MLVL, NM1, NSIZE, PERM, POLES, QSTART, SMLSIZ,
145 $ SMLSZP, SQRE, START, WSTART, Z
146 DOUBLE PRECISION CS, EPS, ORGNRM, P, R, SN
147 * ..
148 * .. External Functions ..
149 LOGICAL LSAME
150 INTEGER ILAENV
151 DOUBLE PRECISION DLAMCH, DLANST
152 EXTERNAL LSAME, ILAENV, DLAMCH, DLANST
153 * ..
154 * .. External Subroutines ..
155 EXTERNAL DCOPY, DLARTG, DLASCL, DLASD0, DLASDA, DLASDQ,
156 $ DLASET, DLASR, DSWAP, XERBLA
157 * ..
158 * .. Intrinsic Functions ..
159 INTRINSIC ABS, DBLE, INT, LOG, SIGN
160 * ..
161 * .. Executable Statements ..
162 *
163 * Test the input parameters.
164 *
165 INFO = 0
166 *
167 IUPLO = 0
168 IF( LSAME( UPLO, 'U' ) )
169 $ IUPLO = 1
170 IF( LSAME( UPLO, 'L' ) )
171 $ IUPLO = 2
172 IF( LSAME( COMPQ, 'N' ) ) THEN
173 ICOMPQ = 0
174 ELSE IF( LSAME( COMPQ, 'P' ) ) THEN
175 ICOMPQ = 1
176 ELSE IF( LSAME( COMPQ, 'I' ) ) THEN
177 ICOMPQ = 2
178 ELSE
179 ICOMPQ = -1
180 END IF
181 IF( IUPLO.EQ.0 ) THEN
182 INFO = -1
183 ELSE IF( ICOMPQ.LT.0 ) THEN
184 INFO = -2
185 ELSE IF( N.LT.0 ) THEN
186 INFO = -3
187 ELSE IF( ( LDU.LT.1 ) .OR. ( ( ICOMPQ.EQ.2 ) .AND. ( LDU.LT.
188 $ N ) ) ) THEN
189 INFO = -7
190 ELSE IF( ( LDVT.LT.1 ) .OR. ( ( ICOMPQ.EQ.2 ) .AND. ( LDVT.LT.
191 $ N ) ) ) THEN
192 INFO = -9
193 END IF
194 IF( INFO.NE.0 ) THEN
195 CALL XERBLA( 'DBDSDC', -INFO )
196 RETURN
197 END IF
198 *
199 * Quick return if possible
200 *
201 IF( N.EQ.0 )
202 $ RETURN
203 SMLSIZ = ILAENV( 9, 'DBDSDC', ' ', 0, 0, 0, 0 )
204 IF( N.EQ.1 ) THEN
205 IF( ICOMPQ.EQ.1 ) THEN
206 Q( 1 ) = SIGN( ONE, D( 1 ) )
207 Q( 1+SMLSIZ*N ) = ONE
208 ELSE IF( ICOMPQ.EQ.2 ) THEN
209 U( 1, 1 ) = SIGN( ONE, D( 1 ) )
210 VT( 1, 1 ) = ONE
211 END IF
212 D( 1 ) = ABS( D( 1 ) )
213 RETURN
214 END IF
215 NM1 = N - 1
216 *
217 * If matrix lower bidiagonal, rotate to be upper bidiagonal
218 * by applying Givens rotations on the left
219 *
220 WSTART = 1
221 QSTART = 3
222 IF( ICOMPQ.EQ.1 ) THEN
223 CALL DCOPY( N, D, 1, Q( 1 ), 1 )
224 CALL DCOPY( N-1, E, 1, Q( N+1 ), 1 )
225 END IF
226 IF( IUPLO.EQ.2 ) THEN
227 QSTART = 5
228 WSTART = 2*N - 1
229 DO 10 I = 1, N - 1
230 CALL DLARTG( D( I ), E( I ), CS, SN, R )
231 D( I ) = R
232 E( I ) = SN*D( I+1 )
233 D( I+1 ) = CS*D( I+1 )
234 IF( ICOMPQ.EQ.1 ) THEN
235 Q( I+2*N ) = CS
236 Q( I+3*N ) = SN
237 ELSE IF( ICOMPQ.EQ.2 ) THEN
238 WORK( I ) = CS
239 WORK( NM1+I ) = -SN
240 END IF
241 10 CONTINUE
242 END IF
243 *
244 * If ICOMPQ = 0, use DLASDQ to compute the singular values.
245 *
246 IF( ICOMPQ.EQ.0 ) THEN
247 CALL DLASDQ( 'U', 0, N, 0, 0, 0, D, E, VT, LDVT, U, LDU, U,
248 $ LDU, WORK( WSTART ), INFO )
249 GO TO 40
250 END IF
251 *
252 * If N is smaller than the minimum divide size SMLSIZ, then solve
253 * the problem with another solver.
254 *
255 IF( N.LE.SMLSIZ ) THEN
256 IF( ICOMPQ.EQ.2 ) THEN
257 CALL DLASET( 'A', N, N, ZERO, ONE, U, LDU )
258 CALL DLASET( 'A', N, N, ZERO, ONE, VT, LDVT )
259 CALL DLASDQ( 'U', 0, N, N, N, 0, D, E, VT, LDVT, U, LDU, U,
260 $ LDU, WORK( WSTART ), INFO )
261 ELSE IF( ICOMPQ.EQ.1 ) THEN
262 IU = 1
263 IVT = IU + N
264 CALL DLASET( 'A', N, N, ZERO, ONE, Q( IU+( QSTART-1 )*N ),
265 $ N )
266 CALL DLASET( 'A', N, N, ZERO, ONE, Q( IVT+( QSTART-1 )*N ),
267 $ N )
268 CALL DLASDQ( 'U', 0, N, N, N, 0, D, E,
269 $ Q( IVT+( QSTART-1 )*N ), N,
270 $ Q( IU+( QSTART-1 )*N ), N,
271 $ Q( IU+( QSTART-1 )*N ), N, WORK( WSTART ),
272 $ INFO )
273 END IF
274 GO TO 40
275 END IF
276 *
277 IF( ICOMPQ.EQ.2 ) THEN
278 CALL DLASET( 'A', N, N, ZERO, ONE, U, LDU )
279 CALL DLASET( 'A', N, N, ZERO, ONE, VT, LDVT )
280 END IF
281 *
282 * Scale.
283 *
284 ORGNRM = DLANST( 'M', N, D, E )
285 IF( ORGNRM.EQ.ZERO )
286 $ RETURN
287 CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, IERR )
288 CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, NM1, 1, E, NM1, IERR )
289 *
290 EPS = (0.9D+0)*DLAMCH( 'Epsilon' )
291 *
292 MLVL = INT( LOG( DBLE( N ) / DBLE( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1
293 SMLSZP = SMLSIZ + 1
294 *
295 IF( ICOMPQ.EQ.1 ) THEN
296 IU = 1
297 IVT = 1 + SMLSIZ
298 DIFL = IVT + SMLSZP
299 DIFR = DIFL + MLVL
300 Z = DIFR + MLVL*2
301 IC = Z + MLVL
302 IS = IC + 1
303 POLES = IS + 1
304 GIVNUM = POLES + 2*MLVL
305 *
306 K = 1
307 GIVPTR = 2
308 PERM = 3
309 GIVCOL = PERM + MLVL
310 END IF
311 *
312 DO 20 I = 1, N
313 IF( ABS( D( I ) ).LT.EPS ) THEN
314 D( I ) = SIGN( EPS, D( I ) )
315 END IF
316 20 CONTINUE
317 *
318 START = 1
319 SQRE = 0
320 *
321 DO 30 I = 1, NM1
322 IF( ( ABS( E( I ) ).LT.EPS ) .OR. ( I.EQ.NM1 ) ) THEN
323 *
324 * Subproblem found. First determine its size and then
325 * apply divide and conquer on it.
326 *
327 IF( I.LT.NM1 ) THEN
328 *
329 * A subproblem with E(I) small for I < NM1.
330 *
331 NSIZE = I - START + 1
332 ELSE IF( ABS( E( I ) ).GE.EPS ) THEN
333 *
334 * A subproblem with E(NM1) not too small but I = NM1.
335 *
336 NSIZE = N - START + 1
337 ELSE
338 *
339 * A subproblem with E(NM1) small. This implies an
340 * 1-by-1 subproblem at D(N). Solve this 1-by-1 problem
341 * first.
342 *
343 NSIZE = I - START + 1
344 IF( ICOMPQ.EQ.2 ) THEN
345 U( N, N ) = SIGN( ONE, D( N ) )
346 VT( N, N ) = ONE
347 ELSE IF( ICOMPQ.EQ.1 ) THEN
348 Q( N+( QSTART-1 )*N ) = SIGN( ONE, D( N ) )
349 Q( N+( SMLSIZ+QSTART-1 )*N ) = ONE
350 END IF
351 D( N ) = ABS( D( N ) )
352 END IF
353 IF( ICOMPQ.EQ.2 ) THEN
354 CALL DLASD0( NSIZE, SQRE, D( START ), E( START ),
355 $ U( START, START ), LDU, VT( START, START ),
356 $ LDVT, SMLSIZ, IWORK, WORK( WSTART ), INFO )
357 ELSE
358 CALL DLASDA( ICOMPQ, SMLSIZ, NSIZE, SQRE, D( START ),
359 $ E( START ), Q( START+( IU+QSTART-2 )*N ), N,
360 $ Q( START+( IVT+QSTART-2 )*N ),
361 $ IQ( START+K*N ), Q( START+( DIFL+QSTART-2 )*
362 $ N ), Q( START+( DIFR+QSTART-2 )*N ),
363 $ Q( START+( Z+QSTART-2 )*N ),
364 $ Q( START+( POLES+QSTART-2 )*N ),
365 $ IQ( START+GIVPTR*N ), IQ( START+GIVCOL*N ),
366 $ N, IQ( START+PERM*N ),
367 $ Q( START+( GIVNUM+QSTART-2 )*N ),
368 $ Q( START+( IC+QSTART-2 )*N ),
369 $ Q( START+( IS+QSTART-2 )*N ),
370 $ WORK( WSTART ), IWORK, INFO )
371 END IF
372 IF( INFO.NE.0 ) THEN
373 RETURN
374 END IF
375 START = I + 1
376 END IF
377 30 CONTINUE
378 *
379 * Unscale
380 *
381 CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, IERR )
382 40 CONTINUE
383 *
384 * Use Selection Sort to minimize swaps of singular vectors
385 *
386 DO 60 II = 2, N
387 I = II - 1
388 KK = I
389 P = D( I )
390 DO 50 J = II, N
391 IF( D( J ).GT.P ) THEN
392 KK = J
393 P = D( J )
394 END IF
395 50 CONTINUE
396 IF( KK.NE.I ) THEN
397 D( KK ) = D( I )
398 D( I ) = P
399 IF( ICOMPQ.EQ.1 ) THEN
400 IQ( I ) = KK
401 ELSE IF( ICOMPQ.EQ.2 ) THEN
402 CALL DSWAP( N, U( 1, I ), 1, U( 1, KK ), 1 )
403 CALL DSWAP( N, VT( I, 1 ), LDVT, VT( KK, 1 ), LDVT )
404 END IF
405 ELSE IF( ICOMPQ.EQ.1 ) THEN
406 IQ( I ) = I
407 END IF
408 60 CONTINUE
409 *
410 * If ICOMPQ = 1, use IQ(N,1) as the indicator for UPLO
411 *
412 IF( ICOMPQ.EQ.1 ) THEN
413 IF( IUPLO.EQ.1 ) THEN
414 IQ( N ) = 1
415 ELSE
416 IQ( N ) = 0
417 END IF
418 END IF
419 *
420 * If B is lower bidiagonal, update U by those Givens rotations
421 * which rotated B to be upper bidiagonal
422 *
423 IF( ( IUPLO.EQ.2 ) .AND. ( ICOMPQ.EQ.2 ) )
424 $ CALL DLASR( 'L', 'V', 'B', N, N, WORK( 1 ), WORK( N ), U, LDU )
425 *
426 RETURN
427 *
428 * End of DBDSDC
429 *
430 END