LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
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 wether 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 *> \date September 2012
201 *
202 *> \ingroup auxOTHERauxiliary
203 *
204 *> \par Contributors:
205 * ==================
206 *>
207 *> Ming Gu and Huan Ren, Computer Science Division, University of
208 *> California at Berkeley, USA
209 *>
210 * =====================================================================
211  SUBROUTINE dlasdq( UPLO, SQRE, N, NCVT, NRU, NCC, D, E, VT, LDVT,
212  $ u, ldu, c, ldc, work, info )
213 *
214 * -- LAPACK auxiliary routine (version 3.4.2) --
215 * -- LAPACK is a software package provided by Univ. of Tennessee, --
216 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
217 * September 2012
218 *
219 * .. Scalar Arguments ..
220  CHARACTER uplo
221  INTEGER info, ldc, ldu, ldvt, n, ncc, ncvt, nru, sqre
222 * ..
223 * .. Array Arguments ..
224  DOUBLE PRECISION c( ldc, * ), d( * ), e( * ), u( ldu, * ),
225  $ vt( ldvt, * ), work( * )
226 * ..
227 *
228 * =====================================================================
229 *
230 * .. Parameters ..
231  DOUBLE PRECISION zero
232  parameter( zero = 0.0d+0 )
233 * ..
234 * .. Local Scalars ..
235  LOGICAL rotate
236  INTEGER i, isub, iuplo, j, np1, sqre1
237  DOUBLE PRECISION cs, r, smin, sn
238 * ..
239 * .. External Subroutines ..
240  EXTERNAL dbdsqr, dlartg, dlasr, dswap, xerbla
241 * ..
242 * .. External Functions ..
243  LOGICAL lsame
244  EXTERNAL lsame
245 * ..
246 * .. Intrinsic Functions ..
247  INTRINSIC max
248 * ..
249 * .. Executable Statements ..
250 *
251 * Test the input parameters.
252 *
253  info = 0
254  iuplo = 0
255  IF( lsame( uplo, 'U' ) )
256  $ iuplo = 1
257  IF( lsame( uplo, 'L' ) )
258  $ iuplo = 2
259  IF( iuplo.EQ.0 ) THEN
260  info = -1
261  ELSE IF( ( sqre.LT.0 ) .OR. ( sqre.GT.1 ) ) THEN
262  info = -2
263  ELSE IF( n.LT.0 ) THEN
264  info = -3
265  ELSE IF( ncvt.LT.0 ) THEN
266  info = -4
267  ELSE IF( nru.LT.0 ) THEN
268  info = -5
269  ELSE IF( ncc.LT.0 ) THEN
270  info = -6
271  ELSE IF( ( ncvt.EQ.0 .AND. ldvt.LT.1 ) .OR.
272  $ ( ncvt.GT.0 .AND. ldvt.LT.max( 1, n ) ) ) THEN
273  info = -10
274  ELSE IF( ldu.LT.max( 1, nru ) ) THEN
275  info = -12
276  ELSE IF( ( ncc.EQ.0 .AND. ldc.LT.1 ) .OR.
277  $ ( ncc.GT.0 .AND. ldc.LT.max( 1, n ) ) ) THEN
278  info = -14
279  END IF
280  IF( info.NE.0 ) THEN
281  CALL xerbla( 'DLASDQ', -info )
282  RETURN
283  END IF
284  IF( n.EQ.0 )
285  $ RETURN
286 *
287 * ROTATE is true if any singular vectors desired, false otherwise
288 *
289  rotate = ( ncvt.GT.0 ) .OR. ( nru.GT.0 ) .OR. ( ncc.GT.0 )
290  np1 = n + 1
291  sqre1 = sqre
292 *
293 * If matrix non-square upper bidiagonal, rotate to be lower
294 * bidiagonal. The rotations are on the right.
295 *
296  IF( ( iuplo.EQ.1 ) .AND. ( sqre1.EQ.1 ) ) THEN
297  DO 10 i = 1, n - 1
298  CALL dlartg( d( i ), e( i ), cs, sn, r )
299  d( i ) = r
300  e( i ) = sn*d( i+1 )
301  d( i+1 ) = cs*d( i+1 )
302  IF( rotate ) THEN
303  work( i ) = cs
304  work( n+i ) = sn
305  END IF
306  10 CONTINUE
307  CALL dlartg( d( n ), e( n ), cs, sn, r )
308  d( n ) = r
309  e( n ) = zero
310  IF( rotate ) THEN
311  work( n ) = cs
312  work( n+n ) = sn
313  END IF
314  iuplo = 2
315  sqre1 = 0
316 *
317 * Update singular vectors if desired.
318 *
319  IF( ncvt.GT.0 )
320  $ CALL dlasr( 'L', 'V', 'F', np1, ncvt, work( 1 ),
321  $ work( np1 ), vt, ldvt )
322  END IF
323 *
324 * If matrix lower bidiagonal, rotate to be upper bidiagonal
325 * by applying Givens rotations on the left.
326 *
327  IF( iuplo.EQ.2 ) THEN
328  DO 20 i = 1, n - 1
329  CALL dlartg( d( i ), e( i ), cs, sn, r )
330  d( i ) = r
331  e( i ) = sn*d( i+1 )
332  d( i+1 ) = cs*d( i+1 )
333  IF( rotate ) THEN
334  work( i ) = cs
335  work( n+i ) = sn
336  END IF
337  20 CONTINUE
338 *
339 * If matrix (N+1)-by-N lower bidiagonal, one additional
340 * rotation is needed.
341 *
342  IF( sqre1.EQ.1 ) THEN
343  CALL dlartg( d( n ), e( n ), cs, sn, r )
344  d( n ) = r
345  IF( rotate ) THEN
346  work( n ) = cs
347  work( n+n ) = sn
348  END IF
349  END IF
350 *
351 * Update singular vectors if desired.
352 *
353  IF( nru.GT.0 ) THEN
354  IF( sqre1.EQ.0 ) THEN
355  CALL dlasr( 'R', 'V', 'F', nru, n, work( 1 ),
356  $ work( np1 ), u, ldu )
357  ELSE
358  CALL dlasr( 'R', 'V', 'F', nru, np1, work( 1 ),
359  $ work( np1 ), u, ldu )
360  END IF
361  END IF
362  IF( ncc.GT.0 ) THEN
363  IF( sqre1.EQ.0 ) THEN
364  CALL dlasr( 'L', 'V', 'F', n, ncc, work( 1 ),
365  $ work( np1 ), c, ldc )
366  ELSE
367  CALL dlasr( 'L', 'V', 'F', np1, ncc, work( 1 ),
368  $ work( np1 ), c, ldc )
369  END IF
370  END IF
371  END IF
372 *
373 * Call DBDSQR to compute the SVD of the reduced real
374 * N-by-N upper bidiagonal matrix.
375 *
376  CALL dbdsqr( 'U', n, ncvt, nru, ncc, d, e, vt, ldvt, u, ldu, c,
377  $ ldc, work, info )
378 *
379 * Sort the singular values into ascending order (insertion sort on
380 * singular values, but only one transposition per singular vector)
381 *
382  DO 40 i = 1, n
383 *
384 * Scan for smallest D(I).
385 *
386  isub = i
387  smin = d( i )
388  DO 30 j = i + 1, n
389  IF( d( j ).LT.smin ) THEN
390  isub = j
391  smin = d( j )
392  END IF
393  30 CONTINUE
394  IF( isub.NE.i ) THEN
395 *
396 * Swap singular values and vectors.
397 *
398  d( isub ) = d( i )
399  d( i ) = smin
400  IF( ncvt.GT.0 )
401  $ CALL dswap( ncvt, vt( isub, 1 ), ldvt, vt( i, 1 ), ldvt )
402  IF( nru.GT.0 )
403  $ CALL dswap( nru, u( 1, isub ), 1, u( 1, i ), 1 )
404  IF( ncc.GT.0 )
405  $ CALL dswap( ncc, c( isub, 1 ), ldc, c( i, 1 ), ldc )
406  END IF
407  40 CONTINUE
408 *
409  RETURN
410 *
411 * End of DLASDQ
412 *
413  END