3 * =========== DOCUMENTATION ===========
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
9 *> Download SBDSDC + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sbdsdc.f">
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sbdsdc.f">
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sbdsdc.f">
21 * SUBROUTINE SBDSDC( UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ,
24 * .. Scalar Arguments ..
25 * CHARACTER COMPQ, UPLO
26 * INTEGER INFO, LDU, LDVT, N
28 * .. Array Arguments ..
29 * INTEGER IQ( * ), IWORK( * )
30 * REAL D( * ), E( * ), Q( * ), U( LDU, * ),
31 * $ VT( LDVT, * ), WORK( * )
40 *> SBDSDC computes the singular value decomposition (SVD) of a real
41 *> N-by-N (upper or lower) bidiagonal matrix B: B = U * S * VT,
42 *> using a divide and conquer method, where S is a diagonal matrix
43 *> with non-negative diagonal elements (the singular values of B), and
44 *> U and VT are orthogonal matrices of left and right singular vectors,
45 *> respectively. SBDSDC can be used to compute all singular values,
46 *> and optionally, singular vectors or singular vectors in compact form.
48 *> This code makes very mild assumptions about floating point
49 *> arithmetic. It will work on machines with a guard digit in
50 *> add/subtract, or on those binary machines without guard digits
51 *> which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
52 *> It could conceivably fail on hexadecimal or decimal machines
53 *> without guard digits, but we know of none. See SLASD3 for details.
55 *> The code currently calls SLASDQ if singular values only are desired.
56 *> However, it can be slightly modified to compute singular values
57 *> using the divide and conquer method.
65 *> UPLO is CHARACTER*1
66 *> = 'U': B is upper bidiagonal.
67 *> = 'L': B is lower bidiagonal.
72 *> COMPQ is CHARACTER*1
73 *> Specifies whether singular vectors are to be computed
75 *> = 'N': Compute singular values only;
76 *> = 'P': Compute singular values and compute singular
77 *> vectors in compact form;
78 *> = 'I': Compute singular values and singular vectors.
84 *> The order of the matrix B. N >= 0.
89 *> D is REAL array, dimension (N)
90 *> On entry, the n diagonal elements of the bidiagonal matrix B.
91 *> On exit, if INFO=0, the singular values of B.
96 *> E is REAL array, dimension (N-1)
97 *> On entry, the elements of E contain the offdiagonal
98 *> elements of the bidiagonal matrix whose SVD is desired.
99 *> On exit, E has been destroyed.
104 *> U is REAL array, dimension (LDU,N)
105 *> If COMPQ = 'I', then:
106 *> On exit, if INFO = 0, U contains the left singular vectors
107 *> of the bidiagonal matrix.
108 *> For other values of COMPQ, U is not referenced.
114 *> The leading dimension of the array U. LDU >= 1.
115 *> If singular vectors are desired, then LDU >= max( 1, N ).
120 *> VT is REAL array, dimension (LDVT,N)
121 *> If COMPQ = 'I', then:
122 *> On exit, if INFO = 0, VT**T contains the right singular
123 *> vectors of the bidiagonal matrix.
124 *> For other values of COMPQ, VT is not referenced.
130 *> The leading dimension of the array VT. LDVT >= 1.
131 *> If singular vectors are desired, then LDVT >= max( 1, N ).
136 *> Q is REAL array, dimension (LDQ)
137 *> If COMPQ = 'P', then:
138 *> On exit, if INFO = 0, Q and IQ contain the left
139 *> and right singular vectors in a compact form,
140 *> requiring O(N log N) space instead of 2*N**2.
141 *> In particular, Q contains all the REAL data in
142 *> LDQ >= N*(11 + 2*SMLSIZ + 8*INT(LOG_2(N/(SMLSIZ+1))))
143 *> words of memory, where SMLSIZ is returned by ILAENV and
144 *> is equal to the maximum size of the subproblems at the
145 *> bottom of the computation tree (usually about 25).
146 *> For other values of COMPQ, Q is not referenced.
151 *> IQ is INTEGER array, dimension (LDIQ)
152 *> If COMPQ = 'P', then:
153 *> On exit, if INFO = 0, Q and IQ contain the left
154 *> and right singular vectors in a compact form,
155 *> requiring O(N log N) space instead of 2*N**2.
156 *> In particular, IQ contains all INTEGER data in
157 *> LDIQ >= N*(3 + 3*INT(LOG_2(N/(SMLSIZ+1))))
158 *> words of memory, where SMLSIZ is returned by ILAENV and
159 *> is equal to the maximum size of the subproblems at the
160 *> bottom of the computation tree (usually about 25).
161 *> For other values of COMPQ, IQ is not referenced.
166 *> WORK is REAL array, dimension (MAX(1,LWORK))
167 *> If COMPQ = 'N' then LWORK >= (4 * N).
168 *> If COMPQ = 'P' then LWORK >= (6 * N).
169 *> If COMPQ = 'I' then LWORK >= (3 * N**2 + 4 * N).
174 *> IWORK is INTEGER array, dimension (8*N)
180 *> = 0: successful exit.
181 *> < 0: if INFO = -i, the i-th argument had an illegal value.
182 *> > 0: The algorithm failed to compute a singular value.
183 *> The update process of divide and conquer failed.
189 *> \author Univ. of Tennessee
190 *> \author Univ. of California Berkeley
191 *> \author Univ. of Colorado Denver
196 *> \ingroup auxOTHERcomputational
198 *> \par Contributors:
201 *> Ming Gu and Huan Ren, Computer Science Division, University of
202 *> California at Berkeley, USA
204 * =====================================================================
205 SUBROUTINE SBDSDC( UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ,
206 $ WORK, IWORK, INFO )
208 * -- LAPACK computational routine (version 3.6.1) --
209 * -- LAPACK is a software package provided by Univ. of Tennessee, --
210 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
213 * .. Scalar Arguments ..
214 CHARACTER COMPQ, UPLO
215 INTEGER INFO, LDU, LDVT, N
217 * .. Array Arguments ..
218 INTEGER IQ( * ), IWORK( * )
219 REAL D( * ), E( * ), Q( * ), U( LDU, * ),
220 $ VT( LDVT, * ), WORK( * )
223 * =====================================================================
224 * Changed dimension statement in comment describing E from (N) to
225 * (N-1). Sven, 17 Feb 05.
226 * =====================================================================
230 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0, TWO = 2.0E+0 )
232 * .. Local Scalars ..
233 INTEGER DIFL, DIFR, GIVCOL, GIVNUM, GIVPTR, I, IC,
234 $ ICOMPQ, IERR, II, IS, IU, IUPLO, IVT, J, K, KK,
235 $ MLVL, NM1, NSIZE, PERM, POLES, QSTART, SMLSIZ,
236 $ SMLSZP, SQRE, START, WSTART, Z
237 REAL CS, EPS, ORGNRM, P, R, SN
239 * .. External Functions ..
243 EXTERNAL SLAMCH, SLANST, ILAENV, LSAME
245 * .. External Subroutines ..
246 EXTERNAL SCOPY, SLARTG, SLASCL, SLASD0, SLASDA, SLASDQ,
247 $ SLASET, SLASR, SSWAP, XERBLA
249 * .. Intrinsic Functions ..
250 INTRINSIC REAL, ABS, INT, LOG, SIGN
252 * .. Executable Statements ..
254 * Test the input parameters.
259 IF( LSAME( UPLO, 'U' ) )
261 IF( LSAME( UPLO, 'L' ) )
263 IF( LSAME( COMPQ, 'N' ) ) THEN
265 ELSE IF( LSAME( COMPQ, 'P' ) ) THEN
267 ELSE IF( LSAME( COMPQ, 'I' ) ) THEN
272 IF( IUPLO.EQ.0 ) THEN
274 ELSE IF( ICOMPQ.LT.0 ) THEN
276 ELSE IF( N.LT.0 ) THEN
278 ELSE IF( ( LDU.LT.1 ) .OR. ( ( ICOMPQ.EQ.2 ) .AND. ( LDU.LT.
281 ELSE IF( ( LDVT.LT.1 ) .OR. ( ( ICOMPQ.EQ.2 ) .AND. ( LDVT.LT.
286 CALL XERBLA( 'SBDSDC', -INFO )
290 * Quick return if possible
294 SMLSIZ = ILAENV( 9, 'SBDSDC', ' ', 0, 0, 0, 0 )
296 IF( ICOMPQ.EQ.1 ) THEN
297 Q( 1 ) = SIGN( ONE, D( 1 ) )
298 Q( 1+SMLSIZ*N ) = ONE
299 ELSE IF( ICOMPQ.EQ.2 ) THEN
300 U( 1, 1 ) = SIGN( ONE, D( 1 ) )
303 D( 1 ) = ABS( D( 1 ) )
308 * If matrix lower bidiagonal, rotate to be upper bidiagonal
309 * by applying Givens rotations on the left
313 IF( ICOMPQ.EQ.1 ) THEN
314 CALL SCOPY( N, D, 1, Q( 1 ), 1 )
315 CALL SCOPY( N-1, E, 1, Q( N+1 ), 1 )
317 IF( IUPLO.EQ.2 ) THEN
321 CALL SLARTG( D( I ), E( I ), CS, SN, R )
324 D( I+1 ) = CS*D( I+1 )
325 IF( ICOMPQ.EQ.1 ) THEN
328 ELSE IF( ICOMPQ.EQ.2 ) THEN
335 * If ICOMPQ = 0, use SLASDQ to compute the singular values.
337 IF( ICOMPQ.EQ.0 ) THEN
338 * Ignore WSTART, instead using WORK( 1 ), since the two vectors
339 * for CS and -SN above are added only if ICOMPQ == 2,
340 * and adding them exceeds documented WORK size of 4*n.
341 CALL SLASDQ( 'U', 0, N, 0, 0, 0, D, E, VT, LDVT, U, LDU, U,
342 $ LDU, WORK( 1 ), INFO )
346 * If N is smaller than the minimum divide size SMLSIZ, then solve
347 * the problem with another solver.
349 IF( N.LE.SMLSIZ ) THEN
350 IF( ICOMPQ.EQ.2 ) THEN
351 CALL SLASET( 'A', N, N, ZERO, ONE, U, LDU )
352 CALL SLASET( 'A', N, N, ZERO, ONE, VT, LDVT )
353 CALL SLASDQ( 'U', 0, N, N, N, 0, D, E, VT, LDVT, U, LDU, U,
354 $ LDU, WORK( WSTART ), INFO )
355 ELSE IF( ICOMPQ.EQ.1 ) THEN
358 CALL SLASET( 'A', N, N, ZERO, ONE, Q( IU+( QSTART-1 )*N ),
360 CALL SLASET( 'A', N, N, ZERO, ONE, Q( IVT+( QSTART-1 )*N ),
362 CALL SLASDQ( 'U', 0, N, N, N, 0, D, E,
363 $ Q( IVT+( QSTART-1 )*N ), N,
364 $ Q( IU+( QSTART-1 )*N ), N,
365 $ Q( IU+( QSTART-1 )*N ), N, WORK( WSTART ),
371 IF( ICOMPQ.EQ.2 ) THEN
372 CALL SLASET( 'A', N, N, ZERO, ONE, U, LDU )
373 CALL SLASET( 'A', N, N, ZERO, ONE, VT, LDVT )
378 ORGNRM = SLANST( 'M', N, D, E )
381 CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, IERR )
382 CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, NM1, 1, E, NM1, IERR )
384 EPS = SLAMCH( 'Epsilon' )
386 MLVL = INT( LOG( REAL( N ) / REAL( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1
389 IF( ICOMPQ.EQ.1 ) THEN
398 GIVNUM = POLES + 2*MLVL
407 IF( ABS( D( I ) ).LT.EPS ) THEN
408 D( I ) = SIGN( EPS, D( I ) )
416 IF( ( ABS( E( I ) ).LT.EPS ) .OR. ( I.EQ.NM1 ) ) THEN
418 * Subproblem found. First determine its size and then
419 * apply divide and conquer on it.
423 * A subproblem with E(I) small for I < NM1.
425 NSIZE = I - START + 1
426 ELSE IF( ABS( E( I ) ).GE.EPS ) THEN
428 * A subproblem with E(NM1) not too small but I = NM1.
430 NSIZE = N - START + 1
433 * A subproblem with E(NM1) small. This implies an
434 * 1-by-1 subproblem at D(N). Solve this 1-by-1 problem
437 NSIZE = I - START + 1
438 IF( ICOMPQ.EQ.2 ) THEN
439 U( N, N ) = SIGN( ONE, D( N ) )
441 ELSE IF( ICOMPQ.EQ.1 ) THEN
442 Q( N+( QSTART-1 )*N ) = SIGN( ONE, D( N ) )
443 Q( N+( SMLSIZ+QSTART-1 )*N ) = ONE
445 D( N ) = ABS( D( N ) )
447 IF( ICOMPQ.EQ.2 ) THEN
448 CALL SLASD0( NSIZE, SQRE, D( START ), E( START ),
449 $ U( START, START ), LDU, VT( START, START ),
450 $ LDVT, SMLSIZ, IWORK, WORK( WSTART ), INFO )
452 CALL SLASDA( ICOMPQ, SMLSIZ, NSIZE, SQRE, D( START ),
453 $ E( START ), Q( START+( IU+QSTART-2 )*N ), N,
454 $ Q( START+( IVT+QSTART-2 )*N ),
455 $ IQ( START+K*N ), Q( START+( DIFL+QSTART-2 )*
456 $ N ), Q( START+( DIFR+QSTART-2 )*N ),
457 $ Q( START+( Z+QSTART-2 )*N ),
458 $ Q( START+( POLES+QSTART-2 )*N ),
459 $ IQ( START+GIVPTR*N ), IQ( START+GIVCOL*N ),
460 $ N, IQ( START+PERM*N ),
461 $ Q( START+( GIVNUM+QSTART-2 )*N ),
462 $ Q( START+( IC+QSTART-2 )*N ),
463 $ Q( START+( IS+QSTART-2 )*N ),
464 $ WORK( WSTART ), IWORK, INFO )
475 CALL SLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, IERR )
478 * Use Selection Sort to minimize swaps of singular vectors
485 IF( D( J ).GT.P ) THEN
493 IF( ICOMPQ.EQ.1 ) THEN
495 ELSE IF( ICOMPQ.EQ.2 ) THEN
496 CALL SSWAP( N, U( 1, I ), 1, U( 1, KK ), 1 )
497 CALL SSWAP( N, VT( I, 1 ), LDVT, VT( KK, 1 ), LDVT )
499 ELSE IF( ICOMPQ.EQ.1 ) THEN
504 * If ICOMPQ = 1, use IQ(N,1) as the indicator for UPLO
506 IF( ICOMPQ.EQ.1 ) THEN
507 IF( IUPLO.EQ.1 ) THEN
514 * If B is lower bidiagonal, update U by those Givens rotations
515 * which rotated B to be upper bidiagonal
517 IF( ( IUPLO.EQ.2 ) .AND. ( ICOMPQ.EQ.2 ) )
518 $ CALL SLASR( 'L', 'V', 'B', N, N, WORK( 1 ), WORK( N ), U, LDU )