LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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 lasdq
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 xerbla(srname, info)
Definition cblat2.f:3285
subroutine dbdsqr(uplo, n, ncvt, nru, ncc, d, e, vt, ldvt, u, ldu, c, ldc, work, info)
DBDSQR
Definition dbdsqr.f:241
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
Definition dlartg.f90:111
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 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 dswap(n, dx, incx, dy, incy)
DSWAP
Definition dswap.f:82