LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
dlasdq.f
Go to the documentation of this file.
1 *> \brief \b DLASDQ computes the SVD of a real bidiagonal matrix with diagonal d and off-diagonal e. Used by sbdsdc.
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download DLASDQ + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasdq.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasdq.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasdq.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DLASDQ( UPLO, SQRE, N, NCVT, NRU, NCC, D, E, VT, LDVT,
22 * U, LDU, C, LDC, WORK, INFO )
23 *
24 * .. Scalar Arguments ..
25 * CHARACTER UPLO
26 * INTEGER INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU, SQRE
27 * ..
28 * .. Array Arguments ..
29 * DOUBLE PRECISION C( LDC, * ), D( * ), E( * ), U( LDU, * ),
30 * $ VT( LDVT, * ), WORK( * )
31 * ..
32 *
33 *
34 *> \par Purpose:
35 * =============
36 *>
37 *> \verbatim
38 *>
39 *> DLASDQ computes the singular value decomposition (SVD) of a real
40 *> (upper or lower) bidiagonal matrix with diagonal D and offdiagonal
41 *> E, accumulating the transformations if desired. Letting B denote
42 *> the input bidiagonal matrix, the algorithm computes orthogonal
43 *> matrices Q and P such that B = Q * S * P**T (P**T denotes the transpose
44 *> of P). The singular values S are overwritten on D.
45 *>
46 *> The input matrix U is changed to U * Q if desired.
47 *> The input matrix VT is changed to P**T * VT if desired.
48 *> The input matrix C is changed to Q**T * C if desired.
49 *>
50 *> See "Computing Small Singular Values of Bidiagonal Matrices With
51 *> Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,
52 *> LAPACK Working Note #3, for a detailed description of the algorithm.
53 *> \endverbatim
54 *
55 * Arguments:
56 * ==========
57 *
58 *> \param[in] UPLO
59 *> \verbatim
60 *> UPLO is CHARACTER*1
61 *> On entry, UPLO specifies whether the input bidiagonal matrix
62 *> is upper or lower bidiagonal, and whether it is square are
63 *> not.
64 *> UPLO = 'U' or 'u' B is upper bidiagonal.
65 *> UPLO = 'L' or 'l' B is lower bidiagonal.
66 *> \endverbatim
67 *>
68 *> \param[in] SQRE
69 *> \verbatim
70 *> SQRE is INTEGER
71 *> = 0: then the input matrix is N-by-N.
72 *> = 1: then the input matrix is N-by-(N+1) if UPLU = 'U' and
73 *> (N+1)-by-N if UPLU = 'L'.
74 *>
75 *> The bidiagonal matrix has
76 *> N = NL + NR + 1 rows and
77 *> M = N + SQRE >= N columns.
78 *> \endverbatim
79 *>
80 *> \param[in] N
81 *> \verbatim
82 *> N is INTEGER
83 *> On entry, N specifies the number of rows and columns
84 *> in the matrix. N must be at least 0.
85 *> \endverbatim
86 *>
87 *> \param[in] NCVT
88 *> \verbatim
89 *> NCVT is INTEGER
90 *> On entry, NCVT specifies the number of columns of
91 *> the matrix VT. NCVT must be at least 0.
92 *> \endverbatim
93 *>
94 *> \param[in] NRU
95 *> \verbatim
96 *> NRU is INTEGER
97 *> On entry, NRU specifies the number of rows of
98 *> the matrix U. NRU must be at least 0.
99 *> \endverbatim
100 *>
101 *> \param[in] NCC
102 *> \verbatim
103 *> NCC is INTEGER
104 *> On entry, NCC specifies the number of columns of
105 *> the matrix C. NCC must be at least 0.
106 *> \endverbatim
107 *>
108 *> \param[in,out] D
109 *> \verbatim
110 *> D is DOUBLE PRECISION array, dimension (N)
111 *> On entry, D contains the diagonal entries of the
112 *> bidiagonal matrix whose SVD is desired. On normal exit,
113 *> D contains the singular values in ascending order.
114 *> \endverbatim
115 *>
116 *> \param[in,out] E
117 *> \verbatim
118 *> E is DOUBLE PRECISION array.
119 *> dimension is (N-1) if SQRE = 0 and N if SQRE = 1.
120 *> On entry, the entries of E contain the offdiagonal entries
121 *> of the bidiagonal matrix whose SVD is desired. On normal
122 *> exit, E will contain 0. If the algorithm does not converge,
123 *> D and E will contain the diagonal and superdiagonal entries
124 *> of a bidiagonal matrix orthogonally equivalent to the one
125 *> given as input.
126 *> \endverbatim
127 *>
128 *> \param[in,out] VT
129 *> \verbatim
130 *> VT is DOUBLE PRECISION array, dimension (LDVT, NCVT)
131 *> On entry, contains a matrix which on exit has been
132 *> premultiplied by P**T, dimension N-by-NCVT if SQRE = 0
133 *> and (N+1)-by-NCVT if SQRE = 1 (not referenced if NCVT=0).
134 *> \endverbatim
135 *>
136 *> \param[in] LDVT
137 *> \verbatim
138 *> LDVT is INTEGER
139 *> On entry, LDVT specifies the leading dimension of VT as
140 *> declared in the calling (sub) program. LDVT must be at
141 *> least 1. If NCVT is nonzero LDVT must also be at least N.
142 *> \endverbatim
143 *>
144 *> \param[in,out] U
145 *> \verbatim
146 *> U is DOUBLE PRECISION array, dimension (LDU, N)
147 *> On entry, contains a matrix which on exit has been
148 *> postmultiplied by Q, dimension NRU-by-N if SQRE = 0
149 *> and NRU-by-(N+1) if SQRE = 1 (not referenced if NRU=0).
150 *> \endverbatim
151 *>
152 *> \param[in] LDU
153 *> \verbatim
154 *> LDU is INTEGER
155 *> On entry, LDU specifies the leading dimension of U as
156 *> declared in the calling (sub) program. LDU must be at
157 *> least max( 1, NRU ) .
158 *> \endverbatim
159 *>
160 *> \param[in,out] C
161 *> \verbatim
162 *> C is DOUBLE PRECISION array, dimension (LDC, NCC)
163 *> On entry, contains an N-by-NCC matrix which on exit
164 *> has been premultiplied by Q**T dimension N-by-NCC if SQRE = 0
165 *> and (N+1)-by-NCC if SQRE = 1 (not referenced if NCC=0).
166 *> \endverbatim
167 *>
168 *> \param[in] LDC
169 *> \verbatim
170 *> LDC is INTEGER
171 *> On entry, LDC specifies the leading dimension of C as
172 *> declared in the calling (sub) program. LDC must be at
173 *> least 1. If NCC is nonzero, LDC must also be at least N.
174 *> \endverbatim
175 *>
176 *> \param[out] WORK
177 *> \verbatim
178 *> WORK is DOUBLE PRECISION array, dimension (4*N)
179 *> Workspace. Only referenced if one of NCVT, NRU, or NCC is
180 *> nonzero, and if N is at least 2.
181 *> \endverbatim
182 *>
183 *> \param[out] INFO
184 *> \verbatim
185 *> INFO is INTEGER
186 *> On exit, a value of 0 indicates a successful exit.
187 *> If INFO < 0, argument number -INFO is illegal.
188 *> If INFO > 0, the algorithm did not converge, and INFO
189 *> specifies how many superdiagonals did not converge.
190 *> \endverbatim
191 *
192 * Authors:
193 * ========
194 *
195 *> \author Univ. of Tennessee
196 *> \author Univ. of California Berkeley
197 *> \author Univ. of Colorado Denver
198 *> \author NAG Ltd.
199 *
200 *> \ingroup OTHERauxiliary
201 *
202 *> \par Contributors:
203 * ==================
204 *>
205 *> Ming Gu and Huan Ren, Computer Science Division, University of
206 *> California at Berkeley, USA
207 *>
208 * =====================================================================
209  SUBROUTINE dlasdq( UPLO, SQRE, N, NCVT, NRU, NCC, D, E, VT, LDVT,
210  $ U, LDU, C, LDC, WORK, INFO )
211 *
212 * -- LAPACK auxiliary routine --
213 * -- LAPACK is a software package provided by Univ. of Tennessee, --
214 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
215 *
216 * .. Scalar Arguments ..
217  CHARACTER UPLO
218  INTEGER INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU, SQRE
219 * ..
220 * .. Array Arguments ..
221  DOUBLE PRECISION C( LDC, * ), D( * ), E( * ), U( LDU, * ),
222  $ vt( ldvt, * ), work( * )
223 * ..
224 *
225 * =====================================================================
226 *
227 * .. Parameters ..
228  DOUBLE PRECISION ZERO
229  parameter( zero = 0.0d+0 )
230 * ..
231 * .. Local Scalars ..
232  LOGICAL ROTATE
233  INTEGER I, ISUB, IUPLO, J, NP1, SQRE1
234  DOUBLE PRECISION CS, R, SMIN, SN
235 * ..
236 * .. External Subroutines ..
237  EXTERNAL dbdsqr, dlartg, dlasr, dswap, xerbla
238 * ..
239 * .. External Functions ..
240  LOGICAL LSAME
241  EXTERNAL lsame
242 * ..
243 * .. Intrinsic Functions ..
244  INTRINSIC max
245 * ..
246 * .. Executable Statements ..
247 *
248 * Test the input parameters.
249 *
250  info = 0
251  iuplo = 0
252  IF( lsame( uplo, 'U' ) )
253  $ iuplo = 1
254  IF( lsame( uplo, 'L' ) )
255  $ iuplo = 2
256  IF( iuplo.EQ.0 ) THEN
257  info = -1
258  ELSE IF( ( sqre.LT.0 ) .OR. ( sqre.GT.1 ) ) THEN
259  info = -2
260  ELSE IF( n.LT.0 ) THEN
261  info = -3
262  ELSE IF( ncvt.LT.0 ) THEN
263  info = -4
264  ELSE IF( nru.LT.0 ) THEN
265  info = -5
266  ELSE IF( ncc.LT.0 ) THEN
267  info = -6
268  ELSE IF( ( ncvt.EQ.0 .AND. ldvt.LT.1 ) .OR.
269  $ ( ncvt.GT.0 .AND. ldvt.LT.max( 1, n ) ) ) THEN
270  info = -10
271  ELSE IF( ldu.LT.max( 1, nru ) ) THEN
272  info = -12
273  ELSE IF( ( ncc.EQ.0 .AND. ldc.LT.1 ) .OR.
274  $ ( ncc.GT.0 .AND. ldc.LT.max( 1, n ) ) ) THEN
275  info = -14
276  END IF
277  IF( info.NE.0 ) THEN
278  CALL xerbla( 'DLASDQ', -info )
279  RETURN
280  END IF
281  IF( n.EQ.0 )
282  $ RETURN
283 *
284 * ROTATE is true if any singular vectors desired, false otherwise
285 *
286  rotate = ( ncvt.GT.0 ) .OR. ( nru.GT.0 ) .OR. ( ncc.GT.0 )
287  np1 = n + 1
288  sqre1 = sqre
289 *
290 * If matrix non-square upper bidiagonal, rotate to be lower
291 * bidiagonal. The rotations are on the right.
292 *
293  IF( ( iuplo.EQ.1 ) .AND. ( sqre1.EQ.1 ) ) THEN
294  DO 10 i = 1, n - 1
295  CALL dlartg( d( i ), e( i ), cs, sn, r )
296  d( i ) = r
297  e( i ) = sn*d( i+1 )
298  d( i+1 ) = cs*d( i+1 )
299  IF( rotate ) THEN
300  work( i ) = cs
301  work( n+i ) = sn
302  END IF
303  10 CONTINUE
304  CALL dlartg( d( n ), e( n ), cs, sn, r )
305  d( n ) = r
306  e( n ) = zero
307  IF( rotate ) THEN
308  work( n ) = cs
309  work( n+n ) = sn
310  END IF
311  iuplo = 2
312  sqre1 = 0
313 *
314 * Update singular vectors if desired.
315 *
316  IF( ncvt.GT.0 )
317  $ CALL dlasr( 'L', 'V', 'F', np1, ncvt, work( 1 ),
318  $ work( np1 ), vt, ldvt )
319  END IF
320 *
321 * If matrix lower bidiagonal, rotate to be upper bidiagonal
322 * by applying Givens rotations on the left.
323 *
324  IF( iuplo.EQ.2 ) THEN
325  DO 20 i = 1, n - 1
326  CALL dlartg( d( i ), e( i ), cs, sn, r )
327  d( i ) = r
328  e( i ) = sn*d( i+1 )
329  d( i+1 ) = cs*d( i+1 )
330  IF( rotate ) THEN
331  work( i ) = cs
332  work( n+i ) = sn
333  END IF
334  20 CONTINUE
335 *
336 * If matrix (N+1)-by-N lower bidiagonal, one additional
337 * rotation is needed.
338 *
339  IF( sqre1.EQ.1 ) THEN
340  CALL dlartg( d( n ), e( n ), cs, sn, r )
341  d( n ) = r
342  IF( rotate ) THEN
343  work( n ) = cs
344  work( n+n ) = sn
345  END IF
346  END IF
347 *
348 * Update singular vectors if desired.
349 *
350  IF( nru.GT.0 ) THEN
351  IF( sqre1.EQ.0 ) THEN
352  CALL dlasr( 'R', 'V', 'F', nru, n, work( 1 ),
353  $ work( np1 ), u, ldu )
354  ELSE
355  CALL dlasr( 'R', 'V', 'F', nru, np1, work( 1 ),
356  $ work( np1 ), u, ldu )
357  END IF
358  END IF
359  IF( ncc.GT.0 ) THEN
360  IF( sqre1.EQ.0 ) THEN
361  CALL dlasr( 'L', 'V', 'F', n, ncc, work( 1 ),
362  $ work( np1 ), c, ldc )
363  ELSE
364  CALL dlasr( 'L', 'V', 'F', np1, ncc, work( 1 ),
365  $ work( np1 ), c, ldc )
366  END IF
367  END IF
368  END IF
369 *
370 * Call DBDSQR to compute the SVD of the reduced real
371 * N-by-N upper bidiagonal matrix.
372 *
373  CALL dbdsqr( 'U', n, ncvt, nru, ncc, d, e, vt, ldvt, u, ldu, c,
374  $ ldc, work, info )
375 *
376 * Sort the singular values into ascending order (insertion sort on
377 * singular values, but only one transposition per singular vector)
378 *
379  DO 40 i = 1, n
380 *
381 * Scan for smallest D(I).
382 *
383  isub = i
384  smin = d( i )
385  DO 30 j = i + 1, n
386  IF( d( j ).LT.smin ) THEN
387  isub = j
388  smin = d( j )
389  END IF
390  30 CONTINUE
391  IF( isub.NE.i ) THEN
392 *
393 * Swap singular values and vectors.
394 *
395  d( isub ) = d( i )
396  d( i ) = smin
397  IF( ncvt.GT.0 )
398  $ CALL dswap( ncvt, vt( isub, 1 ), ldvt, vt( i, 1 ), ldvt )
399  IF( nru.GT.0 )
400  $ CALL dswap( nru, u( 1, isub ), 1, u( 1, i ), 1 )
401  IF( ncc.GT.0 )
402  $ CALL dswap( ncc, c( isub, 1 ), ldc, c( i, 1 ), ldc )
403  END IF
404  40 CONTINUE
405 *
406  RETURN
407 *
408 * End of DLASDQ
409 *
410  END
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
Definition: dlartg.f90:113
subroutine dlasr(SIDE, PIVOT, DIRECT, M, N, C, S, A, LDA)
DLASR applies a sequence of plane rotations to a general rectangular matrix.
Definition: dlasr.f:199
subroutine dlasdq(UPLO, SQRE, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)
DLASDQ computes the SVD of a real bidiagonal matrix with diagonal d and off-diagonal e....
Definition: dlasdq.f:211
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine dbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)
DBDSQR
Definition: dbdsqr.f:241
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
Definition: dswap.f:82