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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlasdq.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack