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: *> \date June 2016
201: *
202: *> \ingroup OTHERauxiliary
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.7.0) --
215: * -- LAPACK is a software package provided by Univ. of Tennessee, --
216: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
217: * June 2016
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
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